Added: Dirac-delta integration for structured patches

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;

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

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
 {

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

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
 {

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

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<std::set<int>>& nodes,
                                    int group, bool ignoreGlobalLM)
 {
+  if (group < 0 || group >= (int)nodes.size())
+    return false;
+
 #if SP_DEBUG > 1
-  auto nit = nodes[group].begin();
+  int n = 0;
-  std::cout <<"\n\t   nodes: "<< *(nit++);
+  std::cout <<"\n\t   nodes:";
-  for (int k = 1; nit != nodes[group].end(); ++nit, k++)
+  for (int node : nodes[group])
-    std::cout << (k%10 > 0 ? " " : "\n\t          ") << *nit;
+    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;
 }

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.
+  //! \return \e true if the groups pass checks, otherwise \e false
-  //! \details This checks that no nodes exist on several threads
+  //!
+  //! \details This checks that no nodes exist on several threads.
   bool checkThreadGroups(const std::vector<std::set<int>>& 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