From 5eb981147ec8ba1f075dc3828fea947028a7a558 Mon Sep 17 00:00:00 2001
From: Arne Morten Kvarving <arne.morten.kvarving@sintef.no>
Date: Thu, 31 Aug 2023 09:39:25 +0200
Subject: [PATCH] ASMs3D::getElementCoordinates: support a separate geometry
 basis

---
 src/ASM/ASMs2D.C |  2 +-
 src/ASM/ASMs3D.C | 59 ++++++++++++++++++++++++++++++++++++++++++++----
 src/ASM/ASMs3D.h | 10 +++++++-
 3 files changed, 64 insertions(+), 7 deletions(-)

diff --git a/src/ASM/ASMs2D.C b/src/ASM/ASMs2D.C
index 4bc7352b..e1734be7 100644
--- a/src/ASM/ASMs2D.C
+++ b/src/ASM/ASMs2D.C
@@ -1396,7 +1396,6 @@ bool ASMs2D::getElementCoordinates (Matrix& X, int iel, bool forceItg) const
 bool ASMs2D::getElementCoordinatesPrm (Matrix& X, double u, double v) const
 {
   const Go::SplineSurface* geo = this->getBasis(ASM::GEOMETRY_BASIS);
-  X.resize(nsd,geo->order_u()*geo->order_v());
 
   if (u < geo->startparam_u() ||
       u > geo->endparam_u() ||
@@ -1411,6 +1410,7 @@ bool ASMs2D::getElementCoordinatesPrm (Matrix& X, double u, double v) const
     nj = geo->basis_v().knotInterval(v) - geo->order_v() + 1;
   }
 
+  X.resize(nsd,geo->order_u()*geo->order_v());
   RealArray::const_iterator cit = geo->coefs_begin();
   for (size_t n = 0; n < X.cols(); n++)
   {
diff --git a/src/ASM/ASMs3D.C b/src/ASM/ASMs3D.C
index c35fddd7..5023458a 100644
--- a/src/ASM/ASMs3D.C
+++ b/src/ASM/ASMs3D.C
@@ -1622,21 +1622,31 @@ Vec3 ASMs3D::getCoord (size_t inod) const
 }
 
 
-bool ASMs3D::getElementCoordinates (Matrix& X, int iel, bool) const
+bool ASMs3D::getElementCoordinates (Matrix& X, int iel, bool forceItg) const
 {
 #ifdef INDEX_CHECK
   if (iel < 1 || (size_t)iel > MNPC.size())
   {
     std::cerr <<" *** ASMs3D::getElementCoordinates: Element index "<< iel
-	      <<" out of range [1,"<< MNPC.size() <<"]."<< std::endl;
+              <<" out of range [1,"<< MNPC.size() <<"]."<< std::endl;
     return false;
   }
 #endif
 
-  X.resize(3,svol->order(0)*svol->order(1)*svol->order(2));
-  int lnod0 = this->getFirstItgElmNode();
+  const Go::SplineVolume* geo = forceItg ? svol : this->getBasis(ASM::GEOMETRY_BASIS);
 
-  RealArray::const_iterator cit = svol->coefs_begin();
+  int lnod0 = this->getFirstItgElmNode();
+  if (geo != svol) {
+    const IJK& nIdx = nodeInd[MNPC[iel-1][lnod0]];
+    double u = *(svol->basis(0).begin() + nIdx.I + svol->order(0) - 1);
+    double v = *(svol->basis(1).begin() + nIdx.J + svol->order(1) - 1);
+    double w = *(svol->basis(2).begin() + nIdx.K + svol->order(2) - 1);
+    return this->getElementCoordinatesPrm(X,u,v,w);
+  }
+
+  X.resize(3,geo->order(0)*geo->order(1)*geo->order(2));
+
+  RealArray::const_iterator cit = geo->coefs_begin();
   for (size_t n = 0; n < X.cols(); n++)
   {
     int ip = this->coeffInd(MNPC[iel-1][n + lnod0])*svol->dimension();
@@ -1653,6 +1663,45 @@ bool ASMs3D::getElementCoordinates (Matrix& X, int iel, bool) const
 }
 
 
+bool ASMs3D::getElementCoordinatesPrm (Matrix& X, double u,
+                                       double v, double w) const
+{
+  const Go::SplineVolume* geo = this->getBasis(ASM::GEOMETRY_BASIS);
+
+  if (u < geo->startparam(0) || u > geo->endparam(0) ||
+      v < geo->startparam(1) || v > geo->endparam(1) ||
+      w < geo->startparam(2) || w > geo->endparam(2))
+    return false;
+
+  int ni, nj, nk;
+#pragma omp critical
+  {
+    ni = geo->basis(0).knotInterval(u) - geo->order(0) + 1;
+    nj = geo->basis(1).knotInterval(v) - geo->order(1) + 1;
+    nk = geo->basis(2).knotInterval(w) - geo->order(2) + 1;
+  }
+
+  X.resize(3,geo->order(0)*geo->order(1)*geo->order(2));
+  RealArray::const_iterator cit = geo->coefs_begin();
+  for (size_t n = 0; n < X.cols(); n++)
+  {
+    const int iu = n % geo->order(0);
+    const int iv = (n / geo->order(0)) % geo->order(1);
+    const int iw = n / (geo->order(0) * geo->order(1));
+    const int ip = (ni + iu + ((nk + iw)*geo->numCoefs(1) + (nj + iv))*geo->numCoefs(0))*geo->dimension();
+
+    for (size_t i = 0; i < 3; i++)
+      X(i+1,n+1) = *(cit+(ip+i));
+  }
+
+#if SP_DEBUG > 2
+  std::cout <<"\nCoordinates for element containing parameters ("
+            << u <<"," << v << "," << w <<"):" << X << std::endl;
+#endif
+  return true;
+}
+
+
 void ASMs3D::getNodalCoordinates (Matrix& X) const
 {
   const int n1 = svol->numCoefs(0);
diff --git a/src/ASM/ASMs3D.h b/src/ASM/ASMs3D.h
index b70f6fd8..ca4e5b2f 100644
--- a/src/ASM/ASMs3D.h
+++ b/src/ASM/ASMs3D.h
@@ -256,7 +256,8 @@ public:
   //! \param[in] iel Element index
   //! \param[out] X 3\f$\times\f$n-matrix, where \a n is the number of nodes
   //! in one element
-  virtual bool getElementCoordinates(Matrix& X, int iel, bool = false) const;
+  //! \param[in] forceItg If true return integration basis element coordinates
+  virtual bool getElementCoordinates(Matrix& X, int iel, bool forceItg = false) const;
 
   //! \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
@@ -701,6 +702,13 @@ protected:
   //! \param[in] dir Parameter direction (0,1,2)
   bool getQuasiInterplParameters(RealArray& prm, int dir) const;
 
+  //! \brief Returns a matrix with nodal coordinates for element spanning given parameters.
+  //! \param[out] X 3\f$\times\f$n-matrix, where \a n is the number of nodes
+  //! \param[in] u First parameter of point
+  //! \param[in] v Second parameter of point
+  //! \param[in] w Third parameter of point
+  bool getElementCoordinatesPrm(Matrix& X, double u, double v, double w) const;
+
   //! \brief Returns the volume in the parameter space for an element.
   //! \param[in] iel Element index
   double getParametricVolume(int iel) const;