ASMsxDmx::evalSolution: support separate geometry

src/ASM/ASMbase.h

@@ -229,8 +229,9 @@ public:
   virtual Vec3 getCoord(size_t inod) const = 0;
   //! \brief Returns a matrix with all nodal coordinates within the patch.
   //! \param[out] X nsd\f$\times\f$n-matrix, where \a n is the number of nodes
+  //! \param[in] geo If true returns coordinates for geometry basis
   //! in the patch
-  virtual void getNodalCoordinates(Matrix& X) const = 0;
+  virtual void getNodalCoordinates(Matrix& X, bool geo = false) const = 0;
   //! \brief Returns a matrix with nodal coordinates for an element.
   //! \param[in] iel 1-based element index local to current patch
   //! \param[in] forceItg If true force returning integration basis coordinates

src/ASM/ASMs1D.C

@@ -715,7 +715,7 @@ bool ASMs1D::getElementNodalRotations (TensorVec& T, size_t iel) const
 }
 
 
-void ASMs1D::getNodalCoordinates (Matrix& X) const
+void ASMs1D::getNodalCoordinates (Matrix& X, bool) const
 {
   const int n1 = curv->numCoefs();
 

src/ASM/ASMs1D.h

@@ -101,7 +101,7 @@ public:
   //! \brief Returns a matrix with all nodal coordinates within the patch.
   //! \param[out] X 3\f$\times\f$n-matrix, where \a n is the number of nodes
   //! in the patch
-  virtual void getNodalCoordinates(Matrix& X) const;
+  virtual void getNodalCoordinates(Matrix& X, bool = false) const;
 
   //! \brief Updates the nodal coordinates for this patch.
   //! \param[in] displ Incremental displacements to update the coordinates with

src/ASM/ASMs1DLag.C

@@ -149,7 +149,7 @@ bool ASMs1DLag::getElementCoordinates (Matrix& X, int iel, bool) const
 }
 
 
-void ASMs1DLag::getNodalCoordinates (Matrix& X) const
+void ASMs1DLag::getNodalCoordinates (Matrix& X, bool) const
 {
   X.resize(nsd,coord.size());
 

src/ASM/ASMs1DLag.h

@@ -60,7 +60,7 @@ public:
   //! \brief Returns a matrix with all nodal coordinates within the patch.
   //! \param[out] X 3\f$\times\f$n-matrix, where \a n is the number of nodes
   //! in the patch
-  virtual void getNodalCoordinates(Matrix& X) const;
+  virtual void getNodalCoordinates(Matrix& X, bool = false) const;
 
   //! \brief Returns the global coordinates for the given node.
   //! \param[in] inod 1-based node index local to current patch

src/ASM/ASMs2D.C

@@ -1430,20 +1430,21 @@ bool ASMs2D::getElementCoordinatesPrm (Matrix& X, double u, double v) const
 }
 
 
-void ASMs2D::getNodalCoordinates (Matrix& X) const
+void ASMs2D::getNodalCoordinates (Matrix& X, bool geo) const
 {
-  const int n1 = surf->numCoefs_u();
-  const int n2 = surf->numCoefs_v();
+  const Go::SplineSurface* spline = geo ? this->getBasis(ASM::GEOMETRY_BASIS) : surf;
+  const int n1 = spline->numCoefs_u();
+  const int n2 = spline->numCoefs_v();
   X.resize(nsd,n1*n2);
 
-  RealArray::const_iterator cit = surf->coefs_begin();
+  RealArray::const_iterator cit = spline->coefs_begin();
   size_t inod = 1;
   for (int i2 = 0; i2 < n2; i2++)
     for (int i1 = 0; i1 < n1; i1++, inod++)
     {
-      int ip = (i2*n1 + i1)*surf->dimension();
+      int ip = (i2*n1 + i1)*spline->dimension();
       for (size_t i = 0; i < nsd; i++)
-	X(i+1,inod) = *(cit+(ip+i));
+        X(i+1,inod) = *(cit+(ip+i));
     }
 }
 

src/ASM/ASMs2D.h

@@ -246,8 +246,9 @@ public:
 
   //! \brief Returns a matrix with all nodal coordinates within the patch.
   //! \param[out] X 3\f$\times\f$n-matrix, where \a n is the number of nodes
+  //! \param[in] geo If true return coordinates of geometry basis
   //! in the patch
-  virtual void getNodalCoordinates(Matrix& X) const;
+  virtual void getNodalCoordinates(Matrix& X, bool geo = false) const;
 
   //! \brief Returns the global coordinates for the given node.
   //! \param[in] inod 1-based node index local to current patch

src/ASM/ASMs2DLag.C

@@ -242,7 +242,7 @@ bool ASMs2DLag::getElementCoordinates (Matrix& X, int iel, bool) const
 }
 
 
-void ASMs2DLag::getNodalCoordinates (Matrix& X) const
+void ASMs2DLag::getNodalCoordinates (Matrix& X, bool) const
 {
   X.resize(nsd,coord.size());
 

src/ASM/ASMs2DLag.h

@@ -116,7 +116,7 @@ public:
   //! \brief Returns a matrix with all nodal coordinates within the patch.
   //! \param[out] X 3\f$\times\f$n-matrix, where \a n is the number of nodes
   //! in the patch
-  virtual void getNodalCoordinates(Matrix& X) const;
+  virtual void getNodalCoordinates(Matrix& X, bool = false) const;
 
   //! \brief Returns the global coordinates for the given node.
   //! \param[in] inod 1-based node index local to current patch

src/ASM/ASMs2Dmx.C

@@ -1032,29 +1032,35 @@ bool ASMs2Dmx::evalSolution (Matrix& sField, const IntegrandBase& integrand,
 {
   sField.resize(0,0);
 
+  const Go::SplineSurface* geo = this->getBasis(ASM::GEOMETRY_BASIS);
+  const bool separateGeometry = geo != surf;
+
   // Evaluate the basis functions and their derivatives at all points
-  std::vector<std::vector<Go::BasisDerivsSf>> splinex(m_basis.size());
+  std::vector<std::vector<Go::BasisDerivsSf>> splinex(m_basis.size() + separateGeometry);
   if (regular)
   {
     for (size_t b = 0; b < m_basis.size(); ++b)
       m_basis[b]->computeBasisGrid(gpar[0],gpar[1],splinex[b]);
+    if (separateGeometry)
+      geo->computeBasisGrid(gpar[0],gpar[1],splinex.back());
   }
   else if (gpar[0].size() == gpar[1].size())
   {
-    for (size_t b = 0; b < m_basis.size(); ++b) {
+    for (size_t b = 0; b < splinex.size(); ++b) {
       splinex[b].resize(gpar[0].size());
+      const Go::SplineSurface* basis = b < m_basis.size() ? m_basis[b].get() : geo;
       for (size_t i = 0; i < splinex[b].size(); i++)
-        m_basis[b]->computeBasis(gpar[0][i],gpar[1][i],splinex[b][i]);
+        basis->computeBasis(gpar[0][i],gpar[1][i],splinex[b][i]);
     }
   }
 
   // Fetch nodal (control point) coordinates
   Matrix Xnod, Xtmp;
-  this->getNodalCoordinates(Xnod);
+  this->getNodalCoordinates(Xnod,true);
 
   MxFiniteElement fe(elem_size,firstIp);
   Vector          solPt;
-  BasisValues     bfs(m_basis.size());
+  BasisValues     bfs(splinex.size());
   Matrix          Jac;
 
   // Evaluate the secondary solution field at each point
@@ -1071,20 +1077,30 @@ bool ASMs2Dmx::evalSolution (Matrix& sField, const IntegrandBase& integrand,
                  splinex[b][i].left_idx,ip[b]);
 
       // Fetch associated control point coordinates
-      if (b == static_cast<size_t>(itgBasis)-1)
+      if (!separateGeometry && b == static_cast<size_t>(itgBasis)-1)
         utl::gather(ip[itgBasis-1], nsd, Xnod, Xtmp);
 
       for (int& c : ip[b]) c += ofs;
       ipa.insert(ipa.end(), ip[b].begin(), ip[b].end());
       ofs += nb[b];
     }
+    if (separateGeometry) {
+      ip.front().clear();
+      scatterInd(geo->numCoefs_u(),geo->numCoefs_v(),
+                 geo->order_u(),geo->order_v(),
+                 splinex.back()[i].left_idx,ip.front());
+
+      utl::gather(ip.front(), nsd, Xnod, Xtmp);
+    }
 
     fe.u = splinex[0][i].param[0];
     fe.v = splinex[0][i].param[1];
 
     // Fetch basis function derivatives at current integration point
-    for (size_t b = 0; b < m_basis.size(); ++b)
-      SplineUtils::extractBasis(splinex[b][i],fe.basis(b+1),bfs[b].dNdu);
+    for (size_t b = 0; b < splinex.size(); ++b)
+      SplineUtils::extractBasis(splinex[b][i],
+                                b < m_basis.size() ? fe.basis(b+1) : bfs[b].N,
+                                bfs[b].dNdu);
 
     // Compute Jacobian inverse of the coordinate mapping and
     // basis function derivatives w.r.t. Cartesian coordinates
@@ -1092,7 +1108,8 @@ bool ASMs2Dmx::evalSolution (Matrix& sField, const IntegrandBase& integrand,
       continue; // skip singular points
 
     // Cartesian coordinates of current integration point
-    utl::Point X4(Xtmp*fe.basis(itgBasis),{fe.u,fe.v});
+    utl::Point X4(Xtmp * (separateGeometry ? bfs.back().N : fe.basis(itgBasis)),
+                  {fe.u,fe.v});
 
     // Now evaluate the solution field
     if (!integrand.evalSol(solPt,fe,X4,ipa,elem_size,nb))

src/ASM/ASMs3D.C

@@ -1702,22 +1702,23 @@ bool ASMs3D::getElementCoordinatesPrm (Matrix& X, double u,
 }
 
 
-void ASMs3D::getNodalCoordinates (Matrix& X) const
+void ASMs3D::getNodalCoordinates (Matrix& X, bool geo) const
 {
-  const int n1 = svol->numCoefs(0);
-  const int n2 = svol->numCoefs(1);
-  const int n3 = svol->numCoefs(2);
+  const Go::SplineVolume* spline = geo ? this->getBasis(ASM::GEOMETRY_BASIS) : svol;
+  const int n1 = spline->numCoefs(0);
+  const int n2 = spline->numCoefs(1);
+  const int n3 = spline->numCoefs(2);
   X.resize(3,n1*n2*n3);
 
-  RealArray::const_iterator cit = svol->coefs_begin();
+  RealArray::const_iterator cit = spline->coefs_begin();
   size_t inod = 1;
   for (int i3 = 0; i3 < n3; i3++)
     for (int i2 = 0; i2 < n2; i2++)
       for (int i1 = 0; i1 < n1; i1++, inod++)
       {
-	int ip = ((i3*n2 + i2)*n1 + i1)*svol->dimension();
-	for (size_t i = 0; i < 3; i++)
-	  X(i+1,inod) = *(cit+(ip+i));
+        int ip = ((i3*n2 + i2)*n1 + i1)*spline->dimension();
+        for (size_t i = 0; i < 3; i++)
+          X(i+1,inod) = *(cit+(ip+i));
       }
 }
 

src/ASM/ASMs3D.h

@@ -261,8 +261,9 @@ public:
 
   //! \brief Returns a matrix with all nodal coordinates within the patch.
   //! \param[out] X 3\f$\times\f$n-matrix, where \a n is the number of nodes
+  //! \param[in] geo If true returns coordinates for geometry basis
   //! in the patch
-  virtual void getNodalCoordinates(Matrix& X) const;
+  virtual void getNodalCoordinates(Matrix& X, bool geo = false) const;
 
   //! \brief Returns the global coordinates for the given node.
   //! \param[in] inod 1-based node index local to current patch

src/ASM/ASMs3DLag.C

@@ -267,7 +267,7 @@ bool ASMs3DLag::getElementCoordinates (Matrix& X, int iel, bool) const
 }
 
 
-void ASMs3DLag::getNodalCoordinates (Matrix& X) const
+void ASMs3DLag::getNodalCoordinates (Matrix& X, bool) const
 {
   X.resize(3,coord.size());
 

src/ASM/ASMs3DLag.h

@@ -116,7 +116,7 @@ public:
   //! \brief Returns a matrix with all nodal coordinates within the patch.
   //! \param[out] X 3\f$\times\f$n-matrix, where \a n is the number of nodes
   //! in the patch
-  virtual void getNodalCoordinates(Matrix& X) const;
+  virtual void getNodalCoordinates(Matrix& X, bool = false) const;
 
   //! \brief Returns the global coordinates for the given node.
   //! \param[in] inod 1-based node index local to current patch

src/ASM/ASMs3Dmx.C

@@ -1149,29 +1149,35 @@ bool ASMs3Dmx::evalSolution (Matrix& sField, const IntegrandBase& integrand,
 {
   sField.resize(0,0);
 
+  const Go::SplineVolume* geo = this->getBasis(ASM::GEOMETRY_BASIS);
+  const bool separateGeometry = geo != svol;
+
   // Evaluate the basis functions and their derivatives at all points
-  std::vector<std::vector<Go::BasisDerivs>> splinex(m_basis.size());
+  std::vector<std::vector<Go::BasisDerivs>> splinex(m_basis.size() + separateGeometry);
   if (regular)
   {
     for (size_t b = 0; b < m_basis.size(); ++b)
       m_basis[b]->computeBasisGrid(gpar[0],gpar[1],gpar[2],splinex[b]);
+    if (separateGeometry)
+      geo->computeBasisGrid(gpar[0],gpar[1],gpar[2],splinex.back());
   }
   else if (gpar[0].size() == gpar[1].size() && gpar[0].size() == gpar[2].size())
   {
-    for (size_t b = 0; b < m_basis.size(); ++b) {
+    for (size_t b = 0; b < splinex.size(); ++b) {
       splinex[b].resize(gpar[0].size());
+      const Go::SplineVolume* basis = b < m_basis.size() ? m_basis[b].get() : geo;
       for (size_t i = 0; i < splinex[b].size(); i++)
-        m_basis[b]->computeBasis(gpar[0][i],gpar[1][i],gpar[2][i],splinex[b][i]);
+        basis->computeBasis(gpar[0][i],gpar[1][i],gpar[2][i],splinex[b][i]);
     }
   }
 
   // Fetch nodal (control point) coordinates
   Matrix Xnod, Xtmp;
-  this->getNodalCoordinates(Xnod);
+  this->getNodalCoordinates(Xnod,true);
 
   MxFiniteElement fe(elem_size,firstIp);
   Vector          solPt;
-  BasisValues     bfs(m_basis.size());
+  BasisValues     bfs(splinex.size());
   Matrix          Jac;
 
   // Evaluate the secondary solution field at each point
@@ -1188,21 +1194,31 @@ bool ASMs3Dmx::evalSolution (Matrix& sField, const IntegrandBase& integrand,
                  splinex[b][i].left_idx,ip[b]);
 
       // Fetch associated control point coordinates
-      if (b == (size_t)itgBasis-1)
+      if (!separateGeometry && b == static_cast<size_t>(itgBasis)-1)
         utl::gather(ip[itgBasis-1], 3, Xnod, Xtmp);
 
       for (int& c : ip[b]) c += ofs;
       ipa.insert(ipa.end(), ip[b].begin(), ip[b].end());
       ofs += nb[b];
     }
+    if (separateGeometry) {
+      ip.front().clear();
+      scatterInd(geo->numCoefs(0),geo->numCoefs(1),geo->numCoefs(2),
+                 geo->order(0),geo->order(1),geo->order(2),
+                 splinex.back()[i].left_idx,ip.front());
+
+      utl::gather(ip.front(), 3, Xnod, Xtmp);
+    }
 
     fe.u = splinex[0][i].param[0];
     fe.v = splinex[0][i].param[1];
     fe.w = splinex[0][i].param[2];
 
     // Fetch basis function derivatives at current integration point
-    for (size_t b = 0; b < m_basis.size(); ++b)
-      SplineUtils::extractBasis(splinex[b][i],fe.basis(b+1),bfs[b].dNdu);
+    for (size_t b = 0; b < splinex.size(); ++b)
+      SplineUtils::extractBasis(splinex[b][i],
+                                b < m_basis.size() ? fe.basis(b+1) : bfs[b].N,
+                                bfs[b].dNdu);
 
     // Compute Jacobian inverse of the coordinate mapping and
     // basis function derivatives w.r.t. Cartesian coordinate
@@ -1210,7 +1226,8 @@ bool ASMs3Dmx::evalSolution (Matrix& sField, const IntegrandBase& integrand,
       continue; // skip singular points
 
     // Cartesian coordinates of current integration point
-    utl::Point X4(Xtmp * fe.basis(itgBasis),{fe.u,fe.v,fe.w});
+    utl::Point X4(Xtmp * (separateGeometry ? bfs.back().N : fe.basis(itgBasis)),
+                  {fe.u,fe.v,fe.w});
 
     // Now evaluate the solution field
     if (!integrand.evalSol(solPt,fe,X4,ipa,elem_size,nb))

src/ASM/ASMsupel.C

@@ -185,7 +185,7 @@ Vec3 ASMsupel::getCoord (size_t inod) const
 }
 
 
-void ASMsupel::getNodalCoordinates (Matrix& X) const
+void ASMsupel::getNodalCoordinates (Matrix& X, bool) const
 {
   X.resize(3,nnod);
   for (size_t inod = 1; inod <= nnod; inod++)

src/ASM/ASMsupel.h

@@ -61,7 +61,7 @@ public:
   //! \brief Returns a matrix with all nodal coordinates within the patch.
   //! \param[out] X nsd\f$\times\f$n-matrix, where \a n is the number of nodes
   //! in the patch
-  virtual void getNodalCoordinates(Matrix& X) const;
+  virtual void getNodalCoordinates(Matrix& X, bool = false) const;
   //! \brief Returns a matrix with nodal coordinates for an element.
   //! \param[in] iel 1-based element index local to current patch
   //! \param[out] X 3\f$\times\f$n-matrix, where \a n is the number of nodes

src/ASM/LR/ASMu2D.h

@@ -170,7 +170,7 @@ public:
   //! \brief Returns a matrix with all nodal coordinates within the patch.
   //! \param[out] X 3\f$\times\f$n-matrix, where \a n is the number of nodes
   //! in the patch
-  virtual void getNodalCoordinates(Matrix& X) const
+  virtual void getNodalCoordinates(Matrix& X, bool = false) const
   { this->getCoordinates(X, nsd, *lrspline); }
 
   //! \brief Returns the global coordinates for the given node.

src/ASM/LR/ASMu3D.C

@@ -713,7 +713,7 @@ bool ASMu3D::getElementCoordinates (Matrix& X, int iel, bool forceItg) const
 }
 
 
-void ASMu3D::getNodalCoordinates (Matrix& X) const
+void ASMu3D::getNodalCoordinates (Matrix& X, bool) const
 {
   X.resize(3,lrspline->nBasisFunctions());
 

src/ASM/LR/ASMu3D.h

@@ -156,7 +156,7 @@ public:
   //! \brief Returns a matrix with all nodal coordinates within the patch.
   //! \param[out] X 3\f$\times\f$n-matrix, where \a n is the number of nodes
   //! in the patch
-  virtual void getNodalCoordinates(Matrix& X) const;
+  virtual void getNodalCoordinates(Matrix& X, bool = false) const;
 
   //! \brief Returns the global coordinates for the given node.
   //! \param[in] inod 1-based node index local to current patch