added: flag in SIM for handling projection through fields

in this case, we cannot assume a global numbering for the projected
bases, so do not apply nodal averaging in SIMbase::project.
rather, store each node successively

src/ASM/ASMbase.h

@@ -473,6 +473,15 @@ public:
   virtual bool evalSolution(Matrix& sField, const Vector& locSol,
                             const int* npe, int nf = 0) const;
 
+  //! \brief Evaluates the projected solution field at all visualization points.
+  //! \param[out] sField Solution field
+  //! \param[in] locSol Solution vector local to current patch
+  //! \param[in] npe Number of visualization nodes over each knot span
+  //! \param[in] nf If nonzero, mixed evaluates nf fields on first basis
+  virtual bool evalProjSolution(Matrix& sField, const Vector& locSol,
+                                const int* npe, int nf = 0) const
+  { return this->evalSolution(sField, locSol, npe, nf); }
+
   //! \brief Evaluates the primary solution field at the given points.
   //! \param[out] sField Solution field
   //! \param[in] locSol Solution vector local to current patch

src/ASM/ASMs2Dmx.C

@@ -1231,3 +1231,44 @@ size_t ASMs2Dmx::getNoProjectionNodes() const
 {
   return projBasis->numCoefs_u() * projBasis->numCoefs_v();
 }
+
+
+bool ASMs2Dmx::evalProjSolution (Matrix& sField, const Vector& locSol,
+                                 const int* npe, int nf) const
+{
+#ifdef SP_DEBUG
+  std::cout <<"ASMu2Dmx::evalProjSolution(Matrix&,const Vector&,const int*,int)\n";
+#endif
+  if (projBasis == m_basis[0])
+    return this->evalSolution(sField, locSol, npe, nf);
+
+  // Compute parameter values of the nodal points
+  std::array<RealArray,2> gpar;
+  for (int dir = 0; dir < 2; dir++)
+    if (!this->getGridParameters(gpar[dir],dir,npe[dir]-1))
+      return false;
+
+  size_t nComp = locSol.size() / this->getNoProjectionNodes();
+  if (nComp*this->getNoProjectionNodes() != locSol.size())
+    return false;
+
+  Fields* f = this->getProjectedFields(locSol, nComp);
+
+  // Evaluate the primary solution field at each point
+  sField.resize(nComp,gpar[0].size()*gpar[1].size());
+  size_t k = 1;
+  for (size_t j = 0; j < gpar[1].size(); j++)
+    for (size_t i = 0; i < gpar[0].size(); i++)
+    {
+      Vector vals;
+      FiniteElement fe;
+      fe.u = gpar[0][i];
+      fe.v = gpar[1][j];
+      f->valueFE(fe, vals);
+      sField.fillColumn(k++, vals);
+    }
+
+  delete f;
+
+  return true;
+}

src/ASM/ASMs2Dmx.h

@@ -192,6 +192,14 @@ public:
   virtual bool evalSolution(Matrix& sField, const IntegrandBase& integrand,
 			    const RealArray* gpar, bool regular = true) const;
 
+  //! \brief Evaluates the projected solution field at all visualization points.
+  //! \param[out] sField Solution field
+  //! \param[in] locSol Solution vector local to current patch
+  //! \param[in] npe Number of visualization nodes over each knot span
+  //! \param[in] nf If nonzero, mixed evaluates nf fields on first basis
+  virtual bool evalProjSolution(Matrix& sField, const Vector& locSol,
+                                const int* npe, int nf = 0) const;
+
   //! \brief Extracts nodal results for this patch from the global vector.
   //! \param[in] globVec Global solution vector in DOF-order
   //! \param[out] nodeVec Nodal result vector for this patch

src/ASM/ASMs3Dmx.C

@@ -1375,3 +1375,46 @@ size_t ASMs3Dmx::getNoProjectionNodes() const
          projBasis->numCoefs(1) *
          projBasis->numCoefs(2);
 }
+
+
+bool ASMs3Dmx::evalProjSolution (Matrix& sField, const Vector& locSol,
+                                 const int* npe, int nf) const
+{
+#ifdef SP_DEBUG
+  std::cout <<"ASMu3Dmx::evalProjSolution(Matrix&,const Vector&,const int*,int)\n";
+#endif
+  if (projBasis == m_basis[0])
+    return this->evalSolution(sField, locSol, npe, nf);
+
+  // Compute parameter values of the nodal points
+  std::array<RealArray,3> gpar;
+  for (int dir = 0; dir < 3; dir++)
+    if (!this->getGridParameters(gpar[dir],dir,npe[dir]-1))
+      return false;
+
+  size_t nComp = locSol.size() / this->getNoProjectionNodes();
+  if (nComp*this->getNoProjectionNodes() != locSol.size())
+    return false;
+
+  Fields* f = this->getProjectedFields(locSol, nComp);
+
+  // Evaluate the primary solution field at each point
+  sField.resize(nComp,gpar[0].size()*gpar[1].size()*gpar[2].size());
+  size_t l = 1;
+  for (size_t k = 0; k < gpar[2].size(); k++)
+    for (size_t j = 0; j < gpar[1].size(); j++)
+      for (size_t i = 0; i < gpar[0].size(); i++)
+      {
+        Vector vals;
+        FiniteElement fe;
+        fe.u = gpar[0][i];
+        fe.v = gpar[1][j];
+        fe.w = gpar[2][k];
+        f->valueFE(fe, vals);
+        sField.fillColumn(l++, vals);
+      }
+
+  delete f;
+
+  return true;
+}

src/ASM/ASMs3Dmx.h

@@ -184,6 +184,14 @@ public:
   virtual bool evalSolution(Matrix& sField, const IntegrandBase& integrand,
 			    const RealArray* gpar, bool regular = true) const;
 
+  //! \brief Evaluates the projected solution field at all visualization points.
+  //! \param[out] sField Solution field
+  //! \param[in] locSol Solution vector local to current patch
+  //! \param[in] npe Number of visualization nodes over each knot span
+  //! \param[in] nf If nonzero, mixed evaluates nf fields on first basis
+  virtual bool evalProjSolution(Matrix& sField, const Vector& locSol,
+                                const int* npe, int nf = 0) const;
+
   //! \brief Extracts nodal results for this patch from the global vector.
   //! \param[in] globVec Global solution vector in DOF-order
   //! \param[out] nodeVec Nodal result vector for this patch

src/ASM/IntegrandBase.h

@@ -372,6 +372,11 @@ public:
   //! \brief Returns whether projections are fed through external means.
   virtual bool hasExternalProjections() const { return false; }
 
+  //! \brief Set projected quantities as fields.
+  //! \param[in] f The field with the info
+  //! \param[in] idx Projection index
+  virtual void setProjectedFields(Fields* f, size_t idx) {}
+
 protected:
   //! \brief Initializes the projected fields for current element.
   bool initProjection(const std::vector<int>& MNPC, LocalIntegral& elmInt);

src/ASM/LR/ASMu2Dmx.C

@@ -1166,3 +1166,46 @@ Fields* ASMu2Dmx::getProjectedFields(const Vector& coefs, size_t nf) const
 
   return nullptr;
 }
+
+
+bool ASMu2Dmx::evalProjSolution (Matrix& sField, const Vector& locSol,
+                                 const int* npe, int nf) const
+{
+#ifdef SP_DEBUG
+  std::cout <<"ASMu2D::evalProjSolution(Matrix&,const Vector&,const int*,int)\n";
+#endif
+  if (projBasis == m_basis[0])
+    return this->evalSolution(sField, locSol, npe, nf);
+
+  // Compute parameter values of the result sampling points
+  std::array<RealArray,2> gpar;
+  for (int dir = 0; dir < 2; dir++)
+    if (!this->getGridParameters(gpar[dir],dir,npe[dir]-1))
+      return false;
+
+  size_t nComp = locSol.size() / this->getNoProjectionNodes();
+  if (nComp*this->getNoProjectionNodes() != locSol.size())
+    return false;
+
+  size_t nPoints = gpar[0].size();
+  if (nPoints != gpar[1].size())
+    return false;
+
+  Fields* f = this->getProjectedFields(locSol, nComp);
+
+  // Evaluate the primary solution field at each point
+  sField.resize(nComp,nPoints);
+  for (size_t i = 0; i < nPoints; i++)
+  {
+    Vector vals;
+    FiniteElement fe;
+    fe.u = gpar[0][i];
+    fe.v = gpar[1][i];
+    f->valueFE(fe, vals);
+    sField.fillColumn(1+i, vals);
+  }
+
+  delete f;
+
+  return true;
+}

src/ASM/LR/ASMu2Dmx.h

@@ -134,6 +134,14 @@ public:
                             const RealArray* gpar, bool = false,
                             int deriv = 0, int nf = 0) const;
 
+  //! \brief Evaluates the projected solution field at all visualization points.
+  //! \param[out] sField Solution field
+  //! \param[in] locSol Solution vector local to current patch
+  //! \param[in] npe Number of visualization nodes over each knot span
+  //! \param[in] nf If nonzero, mixed evaluates nf fields on first basis
+  virtual bool evalProjSolution(Matrix& sField, const Vector& locSol,
+                                const int* npe, int nf = 0) const;
+
   //! \brief Evaluates the secondary solution field at the given points.
   //! \param[out] sField Solution field
   //! \param[in] integrand Object with problem-specific data and methods

src/ASM/LR/ASMu3Dmx.C

@@ -943,3 +943,47 @@ Fields* ASMu3Dmx::getProjectedFields(const Vector& coefs, size_t nf) const
 
   return nullptr;
 }
+
+
+bool ASMu3Dmx::evalProjSolution (Matrix& sField, const Vector& locSol,
+                                 const int* npe, int nf) const
+{
+#ifdef SP_DEBUG
+  std::cout <<"ASMu3Dmx::evalProjSolution(Matrix&,const Vector&,const int*,int)\n";
+#endif
+  if (projBasis == m_basis[0])
+    return this->evalSolution(sField, locSol, npe, nf);
+
+  // Compute parameter values of the result sampling points
+  std::array<RealArray,3> gpar;
+  for (int dir = 0; dir < 3; dir++)
+    if (!this->getGridParameters(gpar[dir],dir,npe[dir]-1))
+      return false;
+
+  size_t nComp = locSol.size() / this->getNoProjectionNodes();
+  if (nComp*this->getNoProjectionNodes() != locSol.size())
+    return false;
+
+  size_t nPoints = gpar[0].size();
+  if (nPoints != gpar[1].size())
+    return false;
+
+  Fields* f = this->getProjectedFields(locSol, nComp);
+
+  // Evaluate the primary solution field at each point
+  sField.resize(nComp,nPoints);
+  for (size_t i = 0; i < nPoints; i++)
+  {
+    Vector vals;
+    FiniteElement fe;
+    fe.u = gpar[0][i];
+    fe.v = gpar[1][i];
+    fe.w = gpar[2][i];
+    f->valueFE(fe, vals);
+    sField.fillColumn(1+i, vals);
+  }
+
+  delete f;
+
+  return true;
+}

src/ASM/LR/ASMu3Dmx.h

@@ -141,6 +141,14 @@ public:
   virtual bool evalSolution(Matrix& sField, const IntegrandBase& integrand,
                             const RealArray* gpar, bool = false) const;
 
+  //! \brief Evaluates the projected solution field at all visualization points.
+  //! \param[out] sField Solution field
+  //! \param[in] locSol Solution vector local to current patch
+  //! \param[in] npe Number of visualization nodes over each knot span
+  //! \param[in] nf If nonzero, mixed evaluates nf fields on first basis
+  virtual bool evalProjSolution(Matrix& sField, const Vector& locSol,
+                                const int* npe, int nf = 0) const;
+
   //! \brief Extracts nodal results for this patch from the global vector.
   //! \param[in] globVec Global solution vector in DOF-order
   //! \param[out] nodeVec Nodal result vector for this patch

src/SIM/SIMbase.C

@@ -1332,6 +1332,7 @@ bool SIMbase::solutionNorms (const TimeDomain& time,
   size_t k, lp = 0;
   ASMbase* pch = nullptr;
   PropertyVec::const_iterator p;
+  size_t projOfs = 0;
   for (p = myProps.begin(); p != myProps.end() && ok; ++p)
     if (p->pcode == Property::MATERIAL)
       if (!(pch = this->getPatch(p->patch)))
@@ -1343,10 +1344,20 @@ bool SIMbase::solutionNorms (const TimeDomain& time,
 	for (k = 0; k < ssol.size(); k++)
           if (ssol[k].empty())
             norm->getProjection(k).clear();
-          else
+          else if (this->fieldProjections()) {
+            size_t ndof = pch->getNoProjectionNodes()*myProblem->getNoFields(2);
+            Vector c(ndof);
+            std::copy(ssol[k].begin()+projOfs,
+                ssol[k].begin()+projOfs+ndof, c.begin());
+            Fields* f = pch->getProjectedFields(c, myProblem->getNoFields(2));
+            norm->setProjectedFields(f, k);
+            projOfs += ndof;
+          } else
             this->extractPatchSolution(ssol[k],norm->getProjection(k),lp-1,nCmp,1);
+
         if (mySol)
           mySol->initPatch(pch->idx);
+
 	ok &= pch->integrate(*norm,globalNorm,time);
         if (norm->getIntegrandType() & IntegrandBase::INTERFACE_TERMS) {
           ASM::InterfaceChecker* iChk = this->getInterfaceChecker(pch->idx);
@@ -1368,10 +1379,20 @@ bool SIMbase::solutionNorms (const TimeDomain& time,
       for (k = 0; k < ssol.size(); k++)
         if (ssol[k].empty())
           norm->getProjection(k).clear();
-        else
+        else if (this->fieldProjections()) {
+          size_t ndof = myModel[i]->getNoProjectionNodes()*myProblem->getNoFields(2);
+          Vector c(ndof);
+          std::copy(ssol[k].begin()+projOfs,
+                    ssol[k].begin()+projOfs+ndof, c.begin());
+          Fields* f = myModel[i]->getProjectedFields(c, myProblem->getNoFields(2));
+          norm->setProjectedFields(f, k);
+          projOfs += ndof;
+        } else
           this->extractPatchSolution(ssol[k],norm->getProjection(k),i,nCmp,1);
+
       if (mySol)
         mySol->initPatch(myModel[i]->idx);
+
       ok &= myModel[i]->integrate(*norm,globalNorm,time);
       if (norm->getIntegrandType() & IntegrandBase::INTERFACE_TERMS) {
         ASM::InterfaceChecker* iChk = this->getInterfaceChecker(myModel[i]->idx);
@@ -1583,6 +1604,17 @@ bool SIMbase::project (Matrix& ssol, const Vector& psol,
     myProblem->initIntegration(time,psol);
   }
 
+  // no nodal averaging - full (potentially discontinuous) representation
+  if (this->fieldProjections()) {
+    ngNodes = 0;
+    for (i = 0; i < myModel.size(); i++)
+    {
+      if (myModel[i]->empty()) continue; // skip empty patches
+      ngNodes += myModel[i]->getNoProjectionNodes();
+    }
+  }
+
+  size_t ofs = 0;
   for (i = 0; i < myModel.size(); i++)
   {
     if (myModel[i]->empty()) continue; // skip empty patches
@@ -1658,30 +1690,41 @@ bool SIMbase::project (Matrix& ssol, const Vector& psol,
       return false;
     }
 
-    size_t nComps = values.rows();
-    size_t nNodes = values.cols();
-    if (ssol.empty())
-      ssol.resize(nComps,ngNodes);
+    // If true, we cannot assume we have a multi-patch numbering for
+    // patch projections, so simply append each vector successively.
+    if (this->fieldProjections()) {
+      if (ssol.empty())
+        ssol.resize(myProblem->getNoFields(2),ngNodes);
 
-    // Nodal averaging for nodes referred to by two or more patches
-    // (these are typically the interface nodes)
-    for (n = 1; n <= nNodes; n++)
-      if (count.empty())
-	ssol.fillColumn(myModel[i]->getNodeID(n),values.getColumn(n));
-      else
-      {
-	int inod = myModel[i]->getNodeID(n);
-	for (j = 1; j <= nComps; j++)
-	  ssol(j,inod) += values(j,n);
-	count(inod) ++;
-      }
+      ssol.fillBlock(values, 1, ofs+1);
+      ofs += myModel[i]->getNoProjectionNodes()*myProblem->getNoFields(2);
+    } else {
+      size_t nComps = values.rows();
+      size_t nNodes = values.cols();
+      if (ssol.empty())
+        ssol.resize(nComps,ngNodes);
+
+      // Nodal averaging for nodes referred to by two or more patches
+      // (these are typically the interface nodes)
+      for (n = 1; n <= nNodes; n++)
+        if (count.empty())
+          ssol.fillColumn(myModel[i]->getNodeID(n),values.getColumn(n));
+        else
+        {
+          int inod = myModel[i]->getNodeID(n);
+          for (j = 1; j <= nComps; j++)
+            ssol(j,inod) += values(j,n);
+          count(inod) ++;
+        }
+    }
   }
 
   // Divide through by count(n) to get the nodal average at the interface nodes
-  for (n = 1; n <= count.size(); n++)
-    if (count(n) > 1.0)
-      for (j = 1; j <= ssol.rows(); j++)
-	ssol(j,n) /= count(n);
+  if (!this->fieldProjections())
+    for (n = 1; n <= count.size(); n++)
+      if (count(n) > 1.0)
+        for (j = 1; j <= ssol.rows(); j++)
+          ssol(j,n) /= count(n);
 
   return true;
 }

src/SIM/SIMbase.h

@@ -481,6 +481,9 @@ public:
   //! the spline basis to obtain the control point values.
   bool evalSecondarySolution(Matrix& field, int pindx) const;
 
+  //! \brief Returns whether or not projections are handled through fields.
+  virtual bool fieldProjections() const { return false; }
+
   //! \brief Returns whether an analytical solution is available or not.
   virtual bool haveAnaSol() const { return mySol ? true : false; }
 

src/SIM/SIMoutput.C

@@ -857,6 +857,7 @@ bool SIMoutput::writeGlvP (const Vector& ssol, int iStep, int& nBlock,
 
   size_t i, j;
   int geomID = myGeomID;
+  size_t projOfs = 0;
   for (i = 0; i < myModel.size(); i++)
   {
     if (myModel[i]->empty()) continue; // skip empty patches
@@ -864,9 +865,16 @@ bool SIMoutput::writeGlvP (const Vector& ssol, int iStep, int& nBlock,
     if (msgLevel > 1)
       IFEM::cout <<"Writing projected solution for patch "<< i+1 << std::endl;
 
+    if (this->fieldProjections()) {
+      size_t ndof = myModel[i]->getNoProjectionNodes() * myProblem->getNoFields(2);
+      lovec.resize(ndof);
+      std::copy(ssol.begin()+projOfs, ssol.begin()+projOfs+ndof, lovec.begin());
+      projOfs += ndof;
+    } else
+      this->extractPatchSolution(ssol,lovec,i,sID.size(),1);
+
     // Evaluate the solution variables at the visualization points
-    this->extractPatchSolution(ssol,lovec,i,sID.size(),1);
-    if (!myModel[i]->evalSolution(field,lovec,opt.nViz))
+    if (!myModel[i]->evalProjSolution(field,lovec,opt.nViz,myProblem->getNoFields(2)))
       return false;
 
     // Write out to VTF-file as scalar fields