Added: Support for integration of dirac-delta functions.

This can be used as interface for element point loads.

src/ASM/ASMbase.h

@@ -440,6 +440,15 @@ public:
                          const TimeDomain& time,
                          const ASM::InterfaceChecker& iChk) { return true; }
 
+  //! \brief Integrates a spatial dirac-delta function over a patch.
+  //! \param integrand Object with problem-specific data and methods
+  //! \param glbInt The integrated quantity
+  //! \param[in] u Parameters of the non-zero point of the dirac-delta function
+  //! \param[in] p Function value at the specified point
+  virtual bool diracPoint(Integrand& integrand, GlobalIntegral& glbInt,
+                          const double* u, const Vec3& p) = 0;
+
+
   // Post-processing methods
   // =======================
 
@@ -449,6 +458,10 @@ public:
   //! \param[out] X The Cartesian coordinates of the point
   //! \return Local node number within the patch that matches the point
   virtual int evalPoint(const double* xi, double* param, Vec3& X) const = 0;
+  //! \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 = 0;
 
   //! \brief Creates a standard FE model of this patch for visualization.
   //! \param[out] grid The generated finite element grid

src/ASM/Integrand.h

@@ -239,6 +239,16 @@ public:
     return this->evalIntMx(elmInt,fe,X,normal);
   }
 
+  //! \brief Evaluates the dirac-delta integrand at a specified point.
+  //! \param elmInt The local integral object to receive the contributions
+  //! \param[in] fe Finite element data of current integration point
+  //! \param[in] pval Function value at the specified point
+  //!
+  //! \details This interface can be used for implementation calculation
+  //! of consistent load vectors due to a point load on the element.
+  virtual bool evalPoint(LocalIntegral& elmInt, const FiniteElement& fe,
+                         const Vec3& pval) { return false; }
+
   //! \brief Finalizes the element quantities after the numerical integration.
   //! \param elmInt The local integral object to receive the contributions
   //! \param[in] fe Nodal and integration point data for current element

src/ASM/LR/ASMu2D.C

@@ -1576,6 +1576,28 @@ bool ASMu2D::integrate (Integrand& integrand, int lIndex,
 }
 
 
+bool ASMu2D::diracPoint (Integrand& integrand, GlobalIntegral& glInt,
+                         const double* u, const Vec3& pval)
+{
+  if (!lrspline) return false;
+
+  FiniteElement fe;
+  fe.iel = 1 + lrspline->getElementContaining(u[0],u[1]);
+  fe.u   = u[0];
+  fe.v   = u[1];
+  if (!this->evaluateBasis(fe)) return false;
+
+  LocalIntegral* A = integrand.getLocalIntegral(MNPC[fe.iel-1].size(),
+                                                fe.iel,true);
+
+  bool ok = integrand.evalPoint(*A,fe,pval) && glInt.assemble(A,fe.iel);
+
+  A->destruct();
+
+  return ok;
+}
+
+
 int ASMu2D::evalPoint (const double* xi, double* param, Vec3& X) const
 {
   if (!lrspline) return -2;
@@ -1601,6 +1623,12 @@ int ASMu2D::evalPoint (const double* xi, double* param, Vec3& X) const
 }
 
 
+int ASMu2D::findElementContaining (const double* param) const
+{
+  return lrspline ? 1 + lrspline->getElementContaining(param[0],param[1]) : -2;
+}
+
+
 bool ASMu2D::getGridParameters (RealArray& prm, int dir, int nSegPerSpan) const
 {
 #ifdef SP_DEBUG

src/ASM/LR/ASMu2D.h

@@ -271,6 +271,14 @@ public:
   virtual bool integrate(Integrand& integrand, int lIndex,
                          GlobalIntegral& glbInt, const TimeDomain& time);
 
+  //! \brief Integrates a spatial dirac-delta function over a patch.
+  //! \param integrand Object with problem-specific data and methods
+  //! \param glbInt 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& integrand, GlobalIntegral& glbInt,
+                          const double* u, const Vec3& pval);
+
   //! \brief Updates the time-dependent in-homogeneous Dirichlet coefficients.
   //! \param[in] func Scalar property fields
   //! \param[in] vfunc Vector property fields
@@ -310,6 +318,10 @@ public:
   //! \return Local node number within the patch that matches the point, if any
   //! \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

src/ASM/LR/ASMu3D.C

@@ -1648,6 +1648,29 @@ bool ASMu3D::integrateEdge (Integrand& integrand, int lEdge,
 }
 
 
+bool ASMu3D::diracPoint (Integrand& integrand, GlobalIntegral& glInt,
+                         const double* u, const Vec3& pval)
+{
+  if (!lrspline) return false;
+
+  FiniteElement fe;
+  fe.iel = 1 + lrspline->getElementContaining(u[0],u[1],u[2]);
+  fe.u   = u[0];
+  fe.v   = u[1];
+  fe.w   = u[2];
+  this->evaluateBasis(fe,0,1);
+
+  LocalIntegral* A = integrand.getLocalIntegral(MNPC[fe.iel-1].size(),
+                                                fe.iel,true);
+
+  bool ok = integrand.evalPoint(*A,fe,pval) && glInt.assemble(A,fe.iel);
+
+  A->destruct();
+
+  return ok;
+}
+
+
 int ASMu3D::evalPoint (const double* xi, double* param, Vec3& X) const
 {
   std::cerr << "ASMu3D::evalPoint(...) is not properly implemented yet :(" << std::endl;
@@ -1687,6 +1710,12 @@ int ASMu3D::evalPoint (const double* xi, double* param, Vec3& X) const
 }
 
 
+int ASMu3D::findElementContaining (const double* param) const
+{
+  return lrspline ? 1 + lrspline->getElementContaining(param[0],param[1],param[2]) : -2;
+}
+
+
 bool ASMu3D::getGridParameters (RealArray& prm, int dir, int nSegPerSpan) const
 {
   if (!lrspline) return false;

src/ASM/LR/ASMu3D.h

@@ -274,6 +274,14 @@ public:
   virtual bool integrateEdge(Integrand& integrand, int lEdge,
                              GlobalIntegral& glbInt, const TimeDomain& time);
 
+  //! \brief Integrates a spatial dirac-delta function over a patch.
+  //! \param integrand Object with problem-specific data and methods
+  //! \param glbInt 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& integrand, GlobalIntegral& glbInt,
+                          const double* u, const Vec3& pval);
+
   //! \brief Updates the time-dependent in-homogeneous Dirichlet coefficients.
   //! \param[in] func Scalar property fields
   //! \param[in] vfunc Vector property fields
@@ -293,6 +301,10 @@ public:
   //! \return Local node number within the patch that matches the point, if any
   //! \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

src/SIM/SIMgeneric.C

@@ -65,3 +65,14 @@ int SIMgeneric::evalPoint (const double* xi, Vec3& X, double* param,
   int inod = pch->evalPoint(xi, param ? param : dummy, X);
   return inod > 0 && global ? pch->getNodeID(inod) : inod;
 }
+
+
+int SIMgeneric::findElementContaining (const double* param,
+                                       int patch, bool global) const
+{
+  ASMbase* pch = this->getPatch(patch,true);
+  if (!pch) return -1;
+
+  int iel = pch->findElementContaining(param);
+  return iel > 0 && global ? pch->getElmID(iel) : iel;
+}

src/SIM/SIMgeneric.h

@@ -55,6 +55,14 @@ public:
   //! \return Patch-local or global node number of node that matches the point
   int evalPoint(const double* xi, Vec3& X, double* param = nullptr,
                 int patch = 1, bool global = false) const;
+
+  //! \brief Returns the element that contains a specified spatial point.
+  //! \param[in] param The parameters of the point in the knot-span domain
+  //! \param[in] patch 1-based patch index containing the point
+  //! \param[in] global If \e true, return global number, otherwise patch-local
+  //! \return Patch-local or global number of the element containing the point
+  int findElementContaining(const double* param,
+                            int patch = 1, bool global = false) const;
 };
 
 #endif