Grid output of selected solution components only

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);

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,62 @@ 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;
+}
+
+
+/*!
+  If \a iel is negative, it is assumed that the basis function values aready
+  have been evaluated in \a N, and the &xi; and &eta; parameters are not used.
+  If \a iel is positive, the basis functions are evaluated at point &xi;,&eta;.
+*/
+
+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;
 }
 
@@ -995,7 +1019,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

@@ -290,6 +290,18 @@ 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.
+  //! \param[in] iel 1-based element index local to current patch
+  //! \param[in] xi First parametric coordinate, ξ ∈ [-1,1]
+  //! \param[in] eta Second parametric coordinate, η ∈ [-1,1]
+  //! \param[in] nCmp Number of solution components
+  //! \param[in] pchSol Nodal values of the field to evaluate
+  //! \param[out] ptSol Solution field values at the specified point
+  //! \param[out] N Basis function values at the sepcified point
+  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

src/ASM/ASMs3DLag.C

@@ -952,7 +952,7 @@ int ASMs3DLag::findElement (double u, double v, double w,
   if (!svol)
   {
     std::cerr <<" *** ASMs3DLag::findElement: No spline geometry"<< std::endl;
-    return -1;
+    return 0;
   }
 
   std::array<std::pair<double,int>,3> knot;
@@ -1001,41 +1001,65 @@ bool ASMs3DLag::evalSolution (Matrix& sField, const Vector& locSol,
   size_t ni   = gpar ? gpar[0].size() : nel;
   size_t nj   = gpar ? gpar[1].size() : 1;
   size_t nk   = gpar ? gpar[2].size() : 1;
-  size_t nen  = p1*p2*p3;
 
   sField.resize(nCmp,ni*nj*nk);
-  Matrix elmSol(nCmp,nen);
-  RealArray N(nen), val;
+  RealArray N, val;
 
   double xi = 0.0, eta = 0.0, zeta = 0.0;
   if (!gpar && !Lagrange::computeBasis(N,p1,xi,p2,eta,p3,zeta))
     return false;
 
-  size_t ip = 1;
+  size_t ip = 0;
   int iel = 0;
   for (size_t k = 0; k < nk; k++)
     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], gpar[2][k],
                                   &xi, &eta, &zeta);
-          if (iel < 1 || iel > static_cast<int>(nel))
-            return false;
-          if (!Lagrange::computeBasis(N,p1,xi,p2,eta,p3,zeta))
-            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,zeta,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;
+}
+
+
+/*!
+  If \a iel is negative, it is assumed that the basis function values aready
+  have been evaluated in \a N, and the &xi;, &eta; and &zeta; parameters are
+  not used. If \a iel is positive, the basis functions are evaluated at the
+  point &xi;,&eta;,&zeta;.
+*/
+
+bool ASMs3DLag::evalSolPt (int iel, double xi, double eta, double zeta,
+                           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,p3,zeta))
+    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;
 }
 

src/ASM/ASMs3DLag.h

@@ -313,6 +313,19 @@ class ASMs3DLag : public ASMs3D, protected ASMLagBase
                          int p1, int p2, int p3, IntMat& MNPC);
 
 protected:
+  //! \brief Evaluates a nodal solution field at specified point in an element.
+  //! \param[in] iel 1-based element index local to current patch
+  //! \param[in] xi First parametric coordinate in range [-1,1]
+  //! \param[in] eta Second parametric coordinate in range [-1,1]
+  //! \param[in] zeta Second parametric coordinate in range [-1,1]
+  //! \param[in] nCmp Number of solution components
+  //! \param[in] pchSol Nodal values of the field to evaluate
+  //! \param[out] ptSol Solution field values at the specified point
+  //! \param[out] N Basis function values at the sepcified point
+  virtual bool evalSolPt(int iel, double xi, double eta, double zeta,
+                         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
   size_t nz; //!< Number of nodes in third parameter direction

src/ASM/ASMu2DLag.C

@@ -21,6 +21,7 @@
 #include "IFEM.h"
 #include <numeric>
 #include <sstream>
+#include <fstream>
 
 #ifdef USE_OPENMP
 #include <omp.h>
@@ -587,3 +588,64 @@ bool ASMu2DLag::integrate (Integrand& integrand,
 
   return ok;
 }
+
+
+/*!
+  If \a iel is negative, it is assumed that the basis function values aready
+  have been evaluated in \a N, and the &xi; and &eta; parameters are not used.
+  If \a iel is positive, the basis functions are evaluated at point &xi;,&eta;.
+  The 4-noded elements  are handled by forwarding to the parent class method.
+  The 3-noded elements (linear triangles) are handled separately by setting up
+  a temporary basis function array \a N3 and not touching the incoming \a N.
+  This method will ignore elements with less than 3 nodes (\a ptSol is empty).
+*/
+
+bool ASMu2DLag::evalSolPt (int iel, double xi, double eta, size_t nCmp,
+                           const Vector& pchSol, RealArray& ptSol,
+                           RealArray& N) const
+{
+  ptSol.clear();
+  size_t jel = iel > 0 ? iel : -iel;
+  if (jel < 1 || jel > MNPC.size())
+    return false;
+  else if (MNPC[--jel].size() > 3) // 4-noded 2D element
+    return this->ASMs2DLag::evalSolPt(iel,xi,eta,nCmp,pchSol,ptSol,N);
+  else if (MNPC[jel].size() < 3)
+    return true; // no results for 0D and 1D elements
+
+  // This element is a linear triangle
+  RealArray N3(3,1.0/3.0);
+  if (iel > 0) N3 = { xi, eta, 1.0-xi-eta };
+  return this->ASMs2DLag::evalSolPt(-1-jel,xi,eta,nCmp,pchSol,ptSol,N3);
+}
+
+
+bool ASMu2DLag::writeXML (const char* fname) const
+{
+  std::ofstream os(fname);
+  if (!os) return false;
+
+  os <<"<patch>\n  <nodes>";
+  for (const Vec3& X : coord)
+    os <<"\n    "<< X;
+  os <<"\n  </nodes>";
+  for (size_t nen = 1; nen <= 4; nen++)
+  {
+    bool haveElms = false;
+    for (const IntVec& mnpc : MNPC)
+      if (mnpc.size() == nen)
+      {
+        if (!haveElms)
+          os <<"\n  <elements nenod=\""<< nen <<"\">";
+        os <<"\n   ";
+        for (int inod : mnpc)
+          os <<" "<< inod;
+        haveElms = true;
+      }
+    if (haveElms)
+      os <<"\n  </elements>";
+  }
+  os <<"\n</patch>\n";
+
+  return true;
+}

src/ASM/ASMu2DLag.h

@@ -136,7 +136,22 @@ class ASMu2DLag : public ASMs2DLag
   //! \note The number of element nodes must be set in \a grid on input.
   virtual bool tesselate(ElementBlock& grid, const int*) const;
 
+  //! \brief Dumps the mesh to the specified XML-file.
+  bool writeXML(const char* fname) const;
+
 protected:
+  //! \brief Evaluates a nodal solution field at specified point in an element.
+  //! \param[in] iel 1-based element index local to current patch
+  //! \param[in] xi First parametric coordinate in range [-1,1]
+  //! \param[in] eta Second parametric coordinate in range [-1,1]
+  //! \param[in] nCmp Number of solution components
+  //! \param[in] pchSol Nodal values of the field to evaluate
+  //! \param[out] ptSol Solution field values at the specified point
+  //! \param[out] N Basis function values at the sepcified point
+  virtual bool evalSolPt(int iel, double xi, double eta,
+                         size_t nCmp, const Vector& pchSol,
+                         RealArray& ptSol, RealArray& N) const;
+
   //! \brief Generate thread groups using multi-coloring.
   void generateThreadGroupsMultiColored(bool silence, bool separateGroup1Noded);
 

src/SIM/SIM1D.C

@@ -148,7 +148,7 @@ bool SIM1D::parseGeometryTag (const tinyxml2::XMLElement* elem)
     utl::getAttribute(elem,"offset",offset);
 
     const tinyxml2::XMLElement* child = elem->FirstChildElement("connection");
-    for (; child; child = child->NextSiblingElement())
+    for (; child; child = child->NextSiblingElement("connection"))
     {
       ASM::Interface ifc;
       if (utl::getAttribute(child,"master",ifc.master))

src/SIM/SIM2D.C

@@ -187,7 +187,7 @@ bool SIM2D::parseGeometryTag (const tinyxml2::XMLElement* elem)
 
     std::vector<Interface> top;
     const tinyxml2::XMLElement* child = elem->FirstChildElement("connection");
-    for (; child; child = child->NextSiblingElement())
+    for (; child; child = child->NextSiblingElement("connection"))
     {
       ASM::Interface ifc;
       bool rever = false, periodic = false;

src/SIM/SIM3D.C

@@ -178,7 +178,7 @@ bool SIM3D::parseGeometryTag (const tinyxml2::XMLElement* elem)
     utl::getAttribute(elem,"offset",offset);
 
     const tinyxml2::XMLElement* child = elem->FirstChildElement("connection");
-    for (; child; child = child->NextSiblingElement())
+    for (; child; child = child->NextSiblingElement("connection"))
     {
       ASM::Interface ifc;
       bool periodic = false;

src/SIM/SIMinput.C

@@ -451,7 +451,7 @@ bool SIMinput::parseBCTag (const tinyxml2::XMLElement* elem)
   else if (!strcasecmp(elem->Value(),"propertycodes"))
   {
     const tinyxml2::XMLElement* code = elem->FirstChildElement("code");
-    for (; code; code = code->NextSiblingElement())
+    for (; code; code = code->NextSiblingElement("code"))
     {
       int icode = 0;
       utl::getAttribute(code,"value",icode);