From d5d3adcdb551b10f16709a61283b6cd232cf9d58 Mon Sep 17 00:00:00 2001 From: Knut Morten Okstad Date: Sat, 28 Jul 2018 14:31:05 +0200 Subject: [PATCH] Added: Dirac-delta integration for structured patches --- src/ASM/ASMs1D.C | 23 ++++++++++++++ src/ASM/ASMs1D.h | 14 +++++++++ src/ASM/ASMs2D.C | 32 +++++++++++++++++++ src/ASM/ASMs2D.h | 19 ++++++++++-- src/ASM/ASMs3D.C | 43 ++++++++++++++++++++++++-- src/ASM/ASMs3D.h | 21 +++++++++++-- src/ASM/ASMstruct.C | 75 ++++++++++++++++++++++++++++++++++++++++++--- src/ASM/ASMstruct.h | 27 ++++++++++++++-- 8 files changed, 240 insertions(+), 14 deletions(-) diff --git a/src/ASM/ASMs1D.C b/src/ASM/ASMs1D.C index beda6a1c..88e875d3 100644 --- a/src/ASM/ASMs1D.C +++ b/src/ASM/ASMs1D.C @@ -668,6 +668,17 @@ const Vector& ASMs1D::getGaussPointParameters (Matrix& uGP, int nGauss, } +void ASMs1D::getElementBorders (int iel, double* ub) const +{ + RealArray::const_iterator uit = curv->basis().begin(); + + // Fetch parameter values of the element ends (knots) + int i = iel-1 + curv->order(); + ub[0] = uit[i-1]; + ub[1] = uit[i]; +} + + double ASMs1D::getElementEnds (int i, Vec3Vec& XC) const { RealArray::const_iterator uit = curv->basis().begin(); @@ -699,6 +710,12 @@ double ASMs1D::getElementEnds (int i, Vec3Vec& XC) const } +void ASMs1D::evaluateBasis (double u, double, double, Vector& N) const +{ + this->extractBasis(u,N); +} + + void ASMs1D::extractBasis (double u, Vector& N) const { N.resize(curv->order()); @@ -1020,6 +1037,12 @@ int ASMs1D::evalPoint (const double* xi, double* param, Vec3& X) const } +int ASMs1D::findElementContaining (const double* param) const +{ + return curv ? 2 + curv->basis().knotInterval(param[0]) - curv->order() : -1; +} + + bool ASMs1D::getGridParameters (RealArray& prm, int nSegPerSpan) const { if (!curv) return false; diff --git a/src/ASM/ASMs1D.h b/src/ASM/ASMs1D.h index d2f96775..24bea1b3 100644 --- a/src/ASM/ASMs1D.h +++ b/src/ASM/ASMs1D.h @@ -191,6 +191,11 @@ public: //! \return 0 if no node (control point) matches this point virtual int evalPoint(const double* xi, double* param, Vec3& X) const; + //! \brief Returns the element that contains a specified spatial point. + //! \param[in] param The parameter of the point in the knot-span domain + //! \return Local element number within the patch that contains the point + virtual int findElementContaining(const double* param) const; + //! \brief Creates a line element model of this patch for visualization. //! \param[out] grid The generated line grid //! \param[in] npe Number of visualization nodes over each knot span @@ -267,6 +272,10 @@ public: //! \param[in] nSegSpan Number of visualization segments over each knot-span virtual bool getGridParameters(RealArray& prm, int nSegSpan) const; + //! \brief Evaluates the basis functions at the specified point. + //! \param[in] u Parameter value of evaluation point + //! \param[out] N Basis function values + virtual void evaluateBasis(double u, double, double, Vector& N) const; //! \brief Establishes the vector with basis function values. //! \param[in] u Parameter value of current integration point //! \param[out] N Basis function values @@ -332,6 +341,11 @@ protected: //! \brief Returns the parametric length on the \a i'th knot-span. double getKnotSpan(int i) const; + //! \brief Computes the element border parameters. + //! \param[in] iel 1-based element index + //! \param[out] u Parameter values of the element borders + virtual void getElementBorders(int iel, double* u) const; + //! \brief Computes the element end coordinates. //! \param[in] i Parameter index for the knot-span element //! \param[out] XC Coordinates of the element corners diff --git a/src/ASM/ASMs2D.C b/src/ASM/ASMs2D.C index 49097c7b..87a183c6 100644 --- a/src/ASM/ASMs2D.C +++ b/src/ASM/ASMs2D.C @@ -1430,6 +1430,17 @@ const Vector& ASMs2D::getGaussPointParameters (Matrix& uGP, int dir, int nGauss, } +void ASMs2D::getElementBorders (int iel, double* ub) const +{ + int p1 = surf->order_u(); + int p2 = surf->order_v(); + int n1 = surf->numCoefs_u() - p1 + 1; + int i1 = p1 + (iel-1) % n1; + int i2 = p2 + (iel-1) / n1; + this->getElementBorders(i1-1,i2-1,ub,ub+2); +} + + void ASMs2D::getElementBorders (int i1, int i2, double* u, double* v) const { RealArray::const_iterator uit = surf->basis(0).begin(); @@ -2348,6 +2359,19 @@ int ASMs2D::evalPoint (const double* xi, double* param, Vec3& X) const } +int ASMs2D::findElementContaining (const double* param) const +{ + if (!surf) return -2; + + int p1 = surf->order_u() - 1; + int p2 = surf->order_u() - 1; + int nel1 = surf->numCoefs_u() - p1; + int uEl = surf->basis_u().knotInterval(param[0]) - p1; + int vEl = surf->basis_v().knotInterval(param[1]) - p2; + return 1 + uEl + vEl*nel1; +} + + bool ASMs2D::getGridParameters (RealArray& prm, int dir, int nSegPerSpan) const { if (!surf) return false; @@ -2792,6 +2816,14 @@ bool ASMs2D::getNoStructElms (int& n1, int& n2, int& n3) const } +void ASMs2D::evaluateBasis (double u, double v, double, Vector& N) const +{ + Go::BasisPtsSf spline; + surf->computeBasis(u,v,spline); + N = spline.basisValues; +} + + void ASMs2D::extractBasis (double u, double v, Vector& N, Matrix& dNdu, bool fromRight) const { diff --git a/src/ASM/ASMs2D.h b/src/ASM/ASMs2D.h index 1c0bf5da..dd7b0052 100644 --- a/src/ASM/ASMs2D.h +++ b/src/ASM/ASMs2D.h @@ -358,6 +358,11 @@ public: //! \return 0 if no node (control point) matches this point virtual int evalPoint(const double* xi, double* param, Vec3& X) const; + //! \brief Returns the element that contains a specified spatial point. + //! \param[in] param The parameters of the point in the knot-span domain + //! \return Local element number within the patch that contains the point + virtual int findElementContaining(const double* param) const; + //! \brief Calculates parameter values for visualization nodal points. //! \param[out] prm Parameter values in given direction for all points //! \param[in] dir Parameter direction (0,1) @@ -534,6 +539,10 @@ protected: //! \param[out] u Parameter values of the west-east borders //! \param[out] v Parameter values of the south-north borders void getElementBorders(int i1, int i2, double* u, double* v) const; + //! \brief Computes the element border parameters. + //! \param[in] iel 1-based element index + //! \param[out] u Parameter values of the element borders + virtual void getElementBorders(int iel, double* u) const; //! \brief Computes the element corner coordinates. //! \param[in] i1 Parameter index in u-direction @@ -546,7 +555,7 @@ protected: //! \brief Generates element groups for multi-threading of interior integrals. //! \param[in] integrand Object with problem-specific data and methods //! \param[in] silence If \e true, suppress threading group outprint - //! \param[in] ignoreGlobalLM If \e true ignore global multipliers in sanity check + //! \param[in] ignoreGlobalLM Sanity check option virtual void generateThreadGroups(const Integrand& integrand, bool silence, bool ignoreGlobalLM); @@ -554,7 +563,7 @@ protected: //! \param[in] strip1 Strip width in first direction //! \param[in] strip2 Strip width in second direction //! \param[in] silence If \e true, suppress threading group outprint - //! \param[in] ignoreGlobalLM If \e true ignore global multipliers in sanity check + //! \param[in] ignoreGlobalLM Sanity check option void generateThreadGroups(size_t strip1, size_t strip2, bool silence, bool ignoreGlobalLM); @@ -595,6 +604,12 @@ public: //! \brief Returns the number of elements on a boundary. virtual size_t getNoBoundaryElms(char lIndex, char ldim) const; + //! \brief Evaluates the basis functions at the specified point. + //! \param[in] u First parameter value of evaluation point + //! \param[in] v Second parameter value of evaluation point + //! \param[out] N Basis function values + virtual void evaluateBasis(double u, double v, double, Vector& N) const; + //! \brief Establishes matrices with basis functions and 1st derivatives. //! \param[in] u First parameter value of current integration point //! \param[in] v Second parameter value of current integration point diff --git a/src/ASM/ASMs3D.C b/src/ASM/ASMs3D.C index 5ffe0ff9..173c7ba8 100644 --- a/src/ASM/ASMs3D.C +++ b/src/ASM/ASMs3D.C @@ -1627,6 +1627,20 @@ const Vector& ASMs3D::getGaussPointParameters (Matrix& uGP, int dir, int nGauss, } +void ASMs3D::getElementBorders (int iel, double* ub) const +{ + int p1 = svol->order(0); + int p2 = svol->order(1); + int p3 = svol->order(2); + int n1 = svol->numCoefs(0) - p1 + 1; + int n2 = svol->numCoefs(1) - p2 + 1; + int i1 = p1 + (iel-1) % n1; + int i2 = p2 + ((iel-1) / n1) % n2; + int i3 = p3 + (iel-1) / (n1*n2); + this->getElementBorders(i1-1,i2-1,i3-1,ub,ub+2,ub+4); +} + + void ASMs3D::getElementBorders (int i1, int i2, int i3, double* u, double* v, double* w) const { @@ -2261,7 +2275,7 @@ bool ASMs3D::integrate (Integrand& integrand, int lIndex, fe.w = gpar[2](1,1); Matrix dNdu, Xnod, Jac; - double param[3] = { fe.u, fe.v, fe.w }; + double param[3] = { fe.u, fe.v, fe.w }; Vec4 X(param); Vec3 normal; double dXidu[3]; @@ -2481,7 +2495,7 @@ bool ASMs3D::integrateEdge (Integrand& integrand, int lEdge, if (gpar[2].size() == 1) fe.zeta = fe.w == svol->startparam(2) ? -1.0 : 1.0; Matrix dNdu, Xnod, Jac; - double param[3] = { 0.0, 0.0, 0.0 }; + double param[3] = { 0.0, 0.0, 0.0 }; Vec4 X(param); Vec3 tang; @@ -2605,6 +2619,23 @@ int ASMs3D::evalPoint (const double* xi, double* param, Vec3& X) const } +int ASMs3D::findElementContaining (const double* param) const +{ + if (!svol) + return -2; + + int p1 = svol->order(0) - 1; + int p2 = svol->order(1) - 1; + int p3 = svol->order(2) - 1; + int nel1 = svol->numCoefs(0) - p1; + int nel2 = svol->numCoefs(1) - p2; + int uEl = svol->basis(0).knotInterval(param[0]) - p1; + int vEl = svol->basis(1).knotInterval(param[1]) - p2; + int wEl = svol->basis(2).knotInterval(param[2]) - p3; + return 1 + uEl + (vEl + wEl*nel2)*nel1; +} + + bool ASMs3D::getGridParameters (RealArray& prm, int dir, int nSegPerSpan) const { if (!svol) return false; @@ -3131,6 +3162,14 @@ bool ASMs3D::getNoStructElms (int& n1, int& n2, int& n3) const } +void ASMs3D::evaluateBasis (double u, double v, double w, Vector& N) const +{ + Go::BasisPts spline; + svol->computeBasis(u,v,w,spline); + N = spline.basisValues; +} + + void ASMs3D::extractBasis (double u, double v, double w, Vector& N, Matrix& dNdu, bool fromRight) const { diff --git a/src/ASM/ASMs3D.h b/src/ASM/ASMs3D.h index a32dbe89..f8464a9d 100644 --- a/src/ASM/ASMs3D.h +++ b/src/ASM/ASMs3D.h @@ -423,6 +423,11 @@ public: //! \return 0 if no node (control point) matches this point virtual int evalPoint(const double* xi, double* param, Vec3& X) const; + //! \brief Returns the element that contains a specified spatial point. + //! \param[in] param The parameters of the point in the knot-span domain + //! \return Local element number within the patch that contains the point + virtual int findElementContaining(const double* param) const; + //! \brief Calculates parameter values for visualization nodal points. //! \param[out] prm Parameter values in given direction for all points //! \param[in] dir Parameter direction (0,1,2) @@ -590,6 +595,11 @@ protected: //! \param[in] dir Local face index of the boundary face double getParametricArea(int iel, int dir) const; + //! \brief Computes the element border parameters. + //! \param[in] iel 1-based element index + //! \param[out] u Parameter values of the element borders + virtual void getElementBorders(int iel, double* u) const; + //! \brief Computes the element border parameters. //! \param[in] i1 Parameter index in u-direction //! \param[in] i2 Parameter index in v-direction @@ -611,7 +621,7 @@ protected: //! \brief Generates element groups for multi-threading of interior integrals. //! \param[in] integrand Object with problem-specific data and methods //! \param[in] silence If \e true, suppress threading group outprint - //! \param[in] ignoreGlobalLM If \e true, ignore global multipliers in sanity check + //! \param[in] ignoreGlobalLM Sanity check option virtual void generateThreadGroups(const Integrand& integrand, bool silence, bool ignoreGlobalLM); //! \brief Generates element groups for multi-threading of interior integrals. @@ -619,7 +629,7 @@ protected: //! \param[in] strip2 Strip width in second direction //! \param[in] strip3 Strip width in third direction //! \param[in] silence If \e true, suppress threading group outprint - //! \param[in] ignoreGlobalLM If \e true, ignore global multipliers in sanity check + //! \param[in] ignoreGlobalLM Sanity check option void generateThreadGroups(size_t strip1, size_t strip2, size_t strip3, bool silence, bool ignoreGlobalLM); //! \brief Generates element groups for multi-threading of boundary integrals. @@ -654,6 +664,13 @@ public: //! \brief Returns the number of elements on a boundary. virtual size_t getNoBoundaryElms(char lIndex, char ldim) const; + //! \brief Evaluates the basis functions at the specified point. + //! \param[in] u First parameter value of evaluation point + //! \param[in] v Second parameter value of evaluation point + //! \param[in] w Third parameter value of evaluation point + //! \param[out] N Basis function values + virtual void evaluateBasis(double u, double v, double w, Vector& N) const; + //! \brief Establishes matrices with basis functions and 1st derivatives. //! \param[in] u First parameter value of current integration point //! \param[in] v Second parameter value of current integration point diff --git a/src/ASM/ASMstruct.C b/src/ASM/ASMstruct.C index a7347874..eee87239 100644 --- a/src/ASM/ASMstruct.C +++ b/src/ASM/ASMstruct.C @@ -12,6 +12,11 @@ //============================================================================== #include "ASMstruct.h" +#include "FiniteElement.h" +#include "GlobalIntegral.h" +#include "LocalIntegral.h" +#include "Integrand.h" + #include "GoTools/geometry/GeomObject.h" @@ -89,22 +94,82 @@ bool ASMstruct::addXNodes (unsigned short int dim, size_t nXn, IntVec& nodes) bool ASMstruct::checkThreadGroups (const std::vector>& nodes, int group, bool ignoreGlobalLM) { + if (group < 0 || group >= (int)nodes.size()) + return false; + #if SP_DEBUG > 1 - auto nit = nodes[group].begin(); - std::cout <<"\n\t nodes: "<< *(nit++); - for (int k = 1; nit != nodes[group].end(); ++nit, k++) - std::cout << (k%10 > 0 ? " " : "\n\t ") << *nit; + int n = 0; + std::cout <<"\n\t nodes:"; + for (int node : nodes[group]) + std::cout << ((++n)%10 ? " " : "\n\t ") << node; #endif bool ok = true; for (int k = 0; k < group; k++) - for (const int& node : nodes[group]) + for (int node : nodes[group]) if ((this->getLMType(node+1) != 'G' || !ignoreGlobalLM) && nodes[k].find(node) != nodes[k].end()) { std::cout <<"\n ** Warning: Node "<< node <<" is present on both" <<" thread "<< k+1 <<" and thread "<< group+1; ok = false; } + if (!ok) std::cout << std::endl; + + return ok; +} + + +bool ASMstruct::diracPoint (Integrand& integr, GlobalIntegral& glInt, + const double* u, const Vec3& pval) +{ + FiniteElement fe; + fe.iel = this->findElementContaining(u); + if (fe.iel < 1 || fe.iel > (int)nel) + { + std::cerr <<" *** ASMstruct::diracPoint: The point"; + for (unsigned char i = 0; i < ndim; i++) std::cerr <<" "<< u[i]; + std::cerr <<" is outside the patch domain."<< std::endl; + return false; + } + +#ifdef INDEX_CHECK + double uElm[6]; + this->getElementBorders(fe.iel,uElm); + for (unsigned char i = 0; i < ndim; i++) + if (u[i] < uElm[2*i] || u[i] > uElm[2*i+1]) + { + unsigned char d; + std::cerr <<" *** ASMstruct::diracPoint: The point"; + for (d = 0; d < ndim; d++) std::cerr <<" "<< u[d]; + std::cerr <<" is not within the domain of the found element "<< fe.iel + <<" which is"; + for (d = 0; d < ndim; d++) + std::cerr <<(d ? "x[":" [") << uElm[2*d] <<","<< uElm[2*d+1] <<"]"; + std::cerr << std::endl; + return false; + } +#if SP_DEBUG > 1 + else + { + unsigned char d; + std::cerr <<" * The point"; + for (d = 0; d < ndim; d++) std::cerr <<" "<< u[d]; + std::cerr <<" is within the domain of element "<< fe.iel <<" which is"; + for (d = 0; d < ndim; d++) + std::cerr << (d ? "x[":" [") << uElm[2*d] <<","<< uElm[2*d+1] <<"]"; + std::cerr << std::endl; + } +#endif +#endif + + if (ndim > 0) fe.u = u[0]; + if (ndim > 1) fe.v = u[1]; + if (ndim > 2) fe.w = u[2]; + this->evaluateBasis(fe.u,fe.v,fe.w,fe.N); + + LocalIntegral* A = integr.getLocalIntegral(MNPC[fe.iel-1].size(),fe.iel,true); + bool ok = integr.evalPoint(*A,fe,pval) && glInt.assemble(A,fe.iel); + A->destruct(); return ok; } diff --git a/src/ASM/ASMstruct.h b/src/ASM/ASMstruct.h index efc265e6..303e3112 100644 --- a/src/ASM/ASMstruct.h +++ b/src/ASM/ASMstruct.h @@ -61,11 +61,26 @@ public: //! \param[out] n3 Number of elements in third (w) direction virtual bool getNoStructElms(int& n1, int& n2, int& n3) const = 0; + //! \brief Evaluates the basis functions at the specified point. + //! \param[in] u First parameter value of evaluation point + //! \param[in] v Second parameter value of evaluation point + //! \param[in] w Third parameter value of evaluation point + //! \param[out] N Basis function values + virtual void evaluateBasis(double u, double v, double w, Vector& N) const = 0; + using ASMbase::evalSolution; //! \brief Projects the secondary solution field onto the primary basis. //! \param[in] integr Object with problem-specific data and methods virtual Go::GeomObject* evalSolution(const IntegrandBase& integr) const = 0; + //! \brief Integrates a spatial dirac-delta function over a patch. + //! \param integr Object with problem-specific data and methods + //! \param glInt The integrated quantity + //! \param[in] u Parameters of the non-zero point of the dirac-delta function + //! \param[in] pval Function value at the specified point + virtual bool diracPoint(Integrand& integr, GlobalIntegral& glInt, + const double* u, const Vec3& pval); + protected: //! \brief Adds extraordinary nodes associated with a patch boundary. //! \param[in] dim Dimension of the boundary @@ -76,15 +91,21 @@ protected: //! of the dimension-specific sub-classes. bool addXNodes(unsigned short int dim, size_t nXn, IntVec& nodes); - //! \brief Perform a sanity check on the thread groups. + //! \brief Performs a sanity check on the thread groups. //! \param[in] nodes The nodes to santiy check //! \param[in] group The group to check for //! \param[in] ignoreGlobalLM If \e true, ignore global lagrange multipliers - //! \return If \e true groups pass checks. - //! \details This checks that no nodes exist on several threads + //! \return \e true if the groups pass checks, otherwise \e false + //! + //! \details This checks that no nodes exist on several threads. bool checkThreadGroups(const std::vector>& nodes, int group, bool ignoreGlobalLM); + //! \brief Computes the element border parameters. + //! \param[in] iel 1-based element index + //! \param[out] u Parameter values of the element borders + virtual void getElementBorders(int iel, double* u) const = 0; + protected: Go::GeomObject* geo; //!< Pointer to the actual spline geometry object