Point searching for 2D lagrange meshes

CMakeLists.txt

@@ -100,8 +100,17 @@ target_sources(
     src/ASM/BDFMats.C
     src/ASM/BlockElmMats.C
     src/ASM/DomainDecomposition.C
-    src/ASM/dslbln.f
-    src/ASM/dsolv3.f
+    src/ASM/Ftn/dcros3.f
+    src/ASM/Ftn/dmuax3.f
+    src/ASM/Ftn/dnorm3.f
+    src/ASM/Ftn/dprint.f
+    src/ASM/Ftn/dslbln.f
+    src/ASM/Ftn/dsolv3.f
+    src/ASM/Ftn/errmsg.f
+    src/ASM/Ftn/q4_inv.f
+    src/ASM/Ftn/shape2.f
+    src/ASM/Ftn/shltri.f
+    src/ASM/Ftn/t3_inv.f
     src/ASM/DualField.C
     src/ASM/ElmMats.C
     src/ASM/Field.C

src/ASM/ASMLagBase.C

@@ -12,6 +12,7 @@
 //==============================================================================
 
 #include "ASMLagBase.h"
+#include "Vec3Oper.h"
 
 
 ASMLagBase::ASMLagBase (const ASMLagBase& patch, bool cpCoord) : coord(myCoord)
@@ -20,15 +21,16 @@ ASMLagBase::ASMLagBase (const ASMLagBase& patch, bool cpCoord) : coord(myCoord)
 }
 
 
-bool ASMLagBase::nodalField (Matrix& field, const Vector& sol, size_t nno) const
+bool ASMLagBase::nodalField (Matrix& field,
+                             const RealArray& sol, size_t nno) const
 {
   size_t nPnts = coord.size();
   size_t nComp = sol.size() / nno;
   if (nno < nPnts || nComp*nno != sol.size())
     return false;
 
   field.resize(nComp,nPnts);
-  const double* u = sol.ptr();
+  const double* u = sol.data();
   for (size_t iPt = 1; iPt <= nPnts; iPt++, u += nComp)
     field.fillColumn(iPt,u);
 
@@ -48,7 +50,22 @@ Vec3 ASMLagBase::getGeometricCenter (const std::vector<int>& MNPC) const
 }
 
 
-bool ASMLagBase::updateCoords (const Vector& displ, unsigned char nsd)
+double ASMLagBase::getBoundingBox (const std::vector<int>& MNPC,
+                                   Vec3Pair& bbox) const
+{
+  bbox.first = bbox.second = coord[MNPC.front()];
+  for (size_t n = 1; n < MNPC.size(); n++)
+    for (int i = 0; i < 3; i++)
+      if (double x = coord[MNPC[n]][i]; x < bbox.first[i])
+        bbox.first[i] = x;
+      else if (x > bbox.second[i])
+        bbox.second[i] = x;
+
+  return (bbox.second - bbox.first).length();
+}
+
+
+bool ASMLagBase::updateCoords (const RealArray& displ, unsigned char nsd)
 {
   if (displ.size() != nsd*coord.size())
   {
@@ -58,7 +75,7 @@ bool ASMLagBase::updateCoords (const Vector& displ, unsigned char nsd)
     return false;
   }
 
-  const double* dpt = displ.ptr();
+  const double* dpt = displ.data();
   for (Vec3& XYZ : myCoord)
   {
     XYZ += RealArray(dpt,dpt+nsd);

src/ASM/ASMLagBase.h

@@ -34,15 +34,17 @@ class ASMLagBase
   virtual ~ASMLagBase() {}
 
   //! \brief Direct nodal evaluation of a solution field.
-  bool nodalField(Matrix& field, const Vector& sol, size_t nno) const;
+  bool nodalField(Matrix& field, const RealArray& sol, size_t nno) const;
 
   //! \brief Returns the geometric center of an element.
   Vec3 getGeometricCenter(const std::vector<int>& MNPC) const;
+  //! \brief Returns the bounding box and diameter of an element.
+  double getBoundingBox(const std::vector<int>& MNPC, Vec3Pair& bbox) const;
 
   //! \brief Updates the nodal coordinates for this patch.
   //! \param[in] displ Incremental displacements to update the coordinates with
   //! \param[in] nsd Number of space dimensions
-  bool updateCoords(const Vector& displ, unsigned char nsd);
+  bool updateCoords(const RealArray& displ, unsigned char nsd);
 
 public:
   //! \brief Updates patch origin by adding a constant to all nodes.

src/ASM/ASMbase.C

@@ -1622,6 +1622,19 @@ bool ASMbase::injectNodeVec (const RealArray& nodeVec, RealArray& globRes,
 }
 
 
+bool ASMbase::findPoints (std::vector<Vec3>& points,
+                          std::vector<PointParams>& locs) const
+{
+  size_t nFound = 0;
+  locs.resize(points.size());
+  for (size_t i = 0; i < points.size(); i++)
+    if ((locs[i].dist = this->findPoint(points[i],locs[i].u)) >= 0.0)
+      ++nFound;
+
+  return nFound == points.size();
+}
+
+
 bool ASMbase::getSolution (Matrix& sField, const Vector& locSol,
                            const IntVec& nodes) const
 {
@@ -1772,7 +1785,7 @@ bool ASMbase::evalSolution (Matrix&, const IntegrandBase&,
 
 
 bool ASMbase::evalSolution (Matrix& sField, const IntegrandBase& integrand,
-                            const IntVec& elements) const
+                            const IntVec& elements, const RealArray*) const
 {
   if (elements.empty())
     return this->evalSolution(sField,integrand,nullptr,true);

src/ASM/ASMbase.h

@@ -317,6 +317,8 @@ class ASMbase
   void printNodes(std::ostream& os) const;
   //! \brief Prints out element connections of this patch to the given stream.
   void printElements(std::ostream& os) const;
+  //! \brief Prints out additional app-dependent element information.
+  virtual void printElmInfo(int, const IntegrandBase*) const {}
 
   //! \brief Increase all global node numbers by \a nshift.
   virtual void shiftGlobalNodeNums(int nshift);
@@ -582,6 +584,23 @@ class ASMbase
   //! \return Distance from the point \a X to the found point
   virtual double findPoint(Vec3& X, double* param) const = 0;
 
+  //! \brief A struct with element and associated local coordinates of a point.
+  struct PointParams
+  {
+    size_t iel;  //!< 1-based element index
+    double u[3]; //!< Parametric coordinates w.r.t. element or patch
+    double dist; //!< Distance to spatial point
+
+    //! \brief Default constructor.
+    explicit PointParams(size_t e = 0) : iel(e), u{0.0,0.0,0.0}, dist(-1.0) {}
+  };
+
+  //! \brief Finds elements and local coordinates matching spatial points.
+  //! \param points List of spatial point coordinates
+  //! \param[out] locs List of matching elements and/or parametric coordinates
+  virtual bool findPoints(std::vector<Vec3>& points,
+                          std::vector<PointParams>& locs) const;
+
   //! \brief Creates a standard FE model of this patch for visualization.
   //! \param[out] grid The generated finite element grid
   //! \param[in] npe Number of visualization nodes over each knot span
@@ -714,13 +733,17 @@ class ASMbase
   //! \param[out] sField Solution field
   //! \param[in] integrand Object with problem-specific data and methods
   //! \param[in] elements List of elements to evaluate for (all if empty)
+  //! \param[in] lpar Local parameter values of the result sampling points
   //!
   //! \details The secondary solution is derived from the primary solution,
   //! which is assumed to be stored within the \a integrand for current patch.
   //! This method evaluates the secondary solution at the center of the
   //! provided list of \a elements, or for all if the list is empty.
+  //! If \a lpar is specified, it is assumed to contain the local parameters
+  //! (w.r.t. the elements) of the evaluation points
   virtual bool evalSolution(Matrix& sField, const IntegrandBase& integrand,
-                            const IntVec& elements) const;
+                            const IntVec& elements,
+                            const RealArray* lpar = nullptr) const;
 
   //! \brief Projects the secondary solution using a (discrete) global L2-fit.
   //! \param[out] sField Secondary solution field control point values

src/ASM/ASMs2DLag.C

@@ -718,7 +718,8 @@ bool ASMs2DLag::tesselate (ElementBlock& grid, const int*) const
 }
 
 
-void ASMs2DLag::constrainEdge (int dir, bool open, int dof, int code, char basis)
+void ASMs2DLag::constrainEdge (int dir, bool open, int dof, int code,
+                               char basis)
 {
   this->ASMs2D::constrainEdge(dir, open, dof, code > 0 ? -code : code, basis);
 }
@@ -741,39 +742,79 @@ bool ASMs2DLag::evalSolution (Matrix& sField, const Vector& locSol,
   size_t nCmp = locSol.size() / this->getNoProjectionNodes();
   size_t ni   = gpar ? gpar[0].size() : nel;
   size_t nj   = gpar ? gpar[1].size() : 1;
-  size_t nen  = p1*p2;
 
   sField.resize(nCmp,ni*nj);
-  Matrix elmSol(nCmp,nen);
-  RealArray N(nen), val;
+  RealArray N, val;
 
   double xi = 0.0, eta = 0.0;
   if (!gpar && !Lagrange::computeBasis(N,p1,xi,p2,eta))
     return false;
 
-  size_t ip = 1;
+  size_t ip = 0;
   int iel = 0;
   for (size_t j = 0; j < nj; j++)
-    for (size_t i = 0; i < ni; i++, ip++)
+    for (size_t i = 0; i < ni; i++)
     {
       if (gpar)
-      {
         iel = this->findElement(gpar[0][i], gpar[1][j], &xi, &eta);
-        if (iel < 1 || iel > static_cast<int>(nel))
-          return false;
-        if (!Lagrange::computeBasis(N,p1,xi,p2,eta))
-          return false;
-      }
-      else
-        iel++;
+      else // evaluate at element centers
+        --iel;
 
-      for (size_t a = 1; a <= nen; a++)
-        elmSol.fillColumn(a, locSol.ptr() + nCmp*MNPC[iel-1][a-1]);
-
-      elmSol.multiply(N,val);
-      sField.fillColumn(ip,val);
+      if (!this->evalSolPt(iel,xi,eta,nCmp,locSol,val,N))
+        return false;
+      else if (!val.empty()) // skip elements with no results
+        sField.fillColumn(++ip,val);
     }
 
+  sField.resize(nCmp,ip);
+  return true;
+}
+
+
+bool ASMs2DLag::evalSolution (Matrix& sField, const Vector& locSol,
+                              const IntVec& elms, const RealArray* lpar) const
+{
+  const size_t nCmp = locSol.size() / this->getNoProjectionNodes();
+  const size_t nPts = lpar ? lpar[0].size() : elms.size();
+  const size_t iPts = sField.cols();
+  sField.resize(nCmp,iPts+nPts);
+
+  RealArray N, ptSol;
+  for (size_t i = 0; i < nPts; i++)
+  {
+    double xi  = lpar ? lpar[0][i] : 0.0;
+    double eta = lpar ? lpar[1][i] : 0.0;
+    if (!this->evalSolPt(elms[i],xi,eta,nCmp,locSol,ptSol,N))
+      return false;
+
+    sField.fillColumn(iPts+i+1,ptSol);
+  }
+
+  return true;
+}
+
+
+bool ASMs2DLag::evalSolPt (int iel, double xi, double eta, size_t nCmp,
+                           const Vector& pchSol, RealArray& ptSol,
+                           RealArray& N) const
+{
+  if (iel < 0 && -iel <= static_cast<int>(nel)) // assume N is already evaluated
+    iel = -iel-1;
+  else if (iel < 1 || iel > static_cast<int>(nel))
+    return false;
+  else if (!Lagrange::computeBasis(N,p1,xi,p2,eta))
+    return false;
+  else
+    --iel;
+
+  ptSol.assign(nCmp,0.0);
+  for (size_t a = 0; a < N.size(); a++)
+  {
+    size_t ip = nCmp*MNPC[iel][a];
+    for (size_t i = 0; i < nCmp; i++)
+      ptSol[i] += N[a]*pchSol[ip++];
+  }
+
   return true;
 }
 
@@ -997,7 +1038,7 @@ int ASMs2DLag::findElement (double u, double v, double* xi, double* eta) const
   if (!surf)
   {
     std::cerr <<" *** ASMs2DLag::findElement: No spline geometry"<< std::endl;
-    return -1;
+    return 0;
   }
 
   int ku, kv;

src/ASM/ASMs2DLag.h

@@ -233,6 +233,20 @@ class ASMs2DLag : public ASMs2D, protected ASMLagBase
                             const RealArray* gpar, bool regular,
                             int, int) const;
 
+  //! \brief Evaluates the primary solution field at the given points.
+  //! \param sField Solution field
+  //! \param[in] locSol Solution vector local to current patch
+  //! \param[in] elms List of elements to evaluate at
+  //! \param[in] lpar Local parameter values of the result sampling points
+  //!
+  //! \details If \a lpar is null, the solution is evaluated at the center of
+  //! each element. Otherwise, it is assumed to contain the local ξ and η
+  //! parameters of the evaluation points, w.r.t. the provided elements.
+  //! \note The resuls are appended as extra columns to output matrix \a sField.
+  //! Therefore, any existing results is preserved.
+  bool evalSolution(Matrix& sField, const Vector& locSol,
+                    const IntVec& elms, const RealArray* lpar = nullptr) const;
+
   //! \brief Evaluates the secondary solution field at all visualization points.
   //! \param[out] sField Solution field
   //! \param[in] integrand Object with problem-specific data and methods
@@ -290,6 +304,11 @@ class ASMs2DLag : public ASMs2D, protected ASMLagBase
   static void createMNPC(size_t nx, size_t ny, int p1, int p2, IntMat& MNPC);
 
 protected:
+  //! \brief Evaluates a nodal solution field at specified point in an element.
+  virtual bool evalSolPt(int iel, double xi, double eta, size_t nCmp,
+                         const Vector& pchSol, RealArray& ptSol,
+                         RealArray& N) const;
+
   size_t nx; //!< Number of nodes in first parameter direction
   size_t ny; //!< Number of nodes in second parameter direction
   int    p1; //!< Polynomial order in first parameter direction