added: ability to read mixed bases from streams

src/ASM/ASMbase.h

@@ -116,7 +116,7 @@ public:
   virtual bool empty() const = 0;
 
   //! \brief Creates an instance by reading the given input stream.
-  virtual bool read(std::istream& is) = 0;
+  virtual bool read(std::istream& is, int basis = 0) = 0;
   //! \brief Writes the geometry/basis of the patch to the given stream.
   virtual bool write(std::ostream& os, int basis = 0) const = 0;
 

src/ASM/ASMs1D.C

@@ -58,7 +58,7 @@ ASMs1D::ASMs1D (const ASMs1D& patch, unsigned char n_f)
 }
 
 
-bool ASMs1D::read (std::istream& is)
+bool ASMs1D::read (std::istream& is, int)
 {
   if (shareFE) return true;
   if (curv) delete curv;

src/ASM/ASMs1D.h

@@ -49,7 +49,7 @@ public:
   // ============================
 
   //! \brief Creates an instance by reading the given input stream.
-  virtual bool read(std::istream&);
+  virtual bool read(std::istream&, int = 0);
   //! \brief Writes the geometry of the SplineCurve object to given stream.
   virtual bool write(std::ostream&, int = 0) const;
 

src/ASM/ASMs2D.C

@@ -111,7 +111,7 @@ void ASMs2D::copyParameterDomain (const ASMbase* other)
 }
 
 
-bool ASMs2D::read (std::istream& is)
+bool ASMs2D::read (std::istream& is, int)
 {
   if (shareFE) return true;
   if (surf) delete surf;

src/ASM/ASMs2D.h

@@ -135,9 +135,9 @@ public:
   // ============================
 
   //! \brief Creates an instance by reading the given input stream.
-  virtual bool read(std::istream&);
+  virtual bool read(std::istream&, int = 0);
   //! \brief Writes the geometry of the SplineSurface object to given stream.
-  virtual bool write(std::ostream&, int = 0) const;
+  virtual bool write(std::ostream&, int) const;
 
   //! \brief Generates the finite element topology data for the patch.
   //! \details The data generated are the element-to-node connectivity array,

src/ASM/ASMs2Dmx.C

@@ -11,6 +11,7 @@
 //!
 //==============================================================================
 
+#include "GoTools/geometry/ObjectHeader.h"
 #include "GoTools/geometry/SplineSurface.h"
 
 #include "ASMs2Dmx.h"
@@ -72,6 +73,28 @@ Go::SplineCurve* ASMs2Dmx::getBoundary (int dir, int basis)
 }
 
 
+bool ASMs2Dmx::read (std::istream& is, int basis)
+{
+  if (basis == 0)
+    return this->ASMs2D::read(is);
+
+  if (basis < 0 || basis > static_cast<int>(nfx.size()))
+    return false;
+
+  if (m_basis.empty()) {
+    m_basis.resize(nfx.size());
+    nb.resize(nfx.size(), 0);
+  }
+
+  Go::ObjectHeader head;
+  m_basis[basis-1] = std::make_shared<Go::SplineSurface>();
+  is >> head >> *m_basis[basis-1];
+  nb[basis-1] = m_basis[basis-1]->numCoefs_u()*m_basis[basis-1]->numCoefs_v();
+
+  return true;
+}
+
+
 bool ASMs2Dmx::write (std::ostream& os, int basis) const
 {
   if (basis == -1)

src/ASM/ASMs2Dmx.h

@@ -50,6 +50,8 @@ public:
   // Methods for model generation
   // ============================
 
+  //! \brief Creates an instance by reading the given input stream.
+  virtual bool read(std::istream&, int basis = 0);
   //! \brief Writes the geometry/basis of the patch to given stream.
   virtual bool write(std::ostream& os, int basis) const;
 

src/ASM/ASMs3D.C

@@ -89,7 +89,7 @@ void ASMs3D::copyParameterDomain (const ASMbase* other)
 }
 
 
-bool ASMs3D::read (std::istream& is)
+bool ASMs3D::read (std::istream& is, int)
 {
   if (shareFE) return true;
   if (svol) delete svol;

src/ASM/ASMs3D.h

@@ -155,9 +155,9 @@ public:
   // ============================
 
   //! \brief Creates an instance by reading the given input stream.
-  virtual bool read(std::istream&);
+  virtual bool read(std::istream&, int = 0);
   //! \brief Writes the geometry of the SplineVolume object to given stream.
-  virtual bool write(std::ostream&, int = 0) const;
+  virtual bool write(std::ostream&, int) const;
 
   //! \brief Generates the finite element topology data for the patch.
   //! \details The data generated are the element-to-node connectivity array,

src/ASM/ASMs3DLagrecovery.C

@@ -0,0 +1,143 @@
+// $Id$
+//==============================================================================
+//!
+//! \file ASMs3DLagrecovery.C
+//!
+//! \date Jul 7 2020
+//!
+//! \author Arne Morten Kvarving / SINTEF
+//!
+//! \brief Recovery of secondary solutions for structured 3D FE models.
+//!
+//==============================================================================
+
+#include "ASMs3DLag.h"
+#include "ItgPoint.h"
+#include "Field.h"
+#include "CoordinateMapping.h"
+#include "GaussQuadrature.h"
+#include "Lagrange.h"
+#include "SparseMatrix.h"
+#include "SplineUtils.h"
+#include "Utilities.h"
+#include "Profiler.h"
+#include <array>
+
+
+/*
+bool ASMs3DLag::assembleL2matrices (SparseMatrix& A, StdVector& B,
+                                    const IntegrandBase& integrand,
+                                    bool continuous) const
+{
+  const size_t nnod = this->getNoNodes(1);
+
+  const int nel1 = (nx-1)/(p1-1);
+  const int nel2 = (ny-1)/(p2-1);
+  const int nel3 = (nz-1)/(p3-1);
+
+  int pmax = p1 > p2 ? p1 : p2;
+  if (pmax < p3) pmax = p3;
+
+  // Get Gaussian quadrature point coordinates (and weights if continuous)
+  const int ng1 = continuous ? this->getNoGaussPt(pmax,true) : p1 - 1;
+  const int ng2 = continuous ? ng1 : p2 - 1;
+  const int ng3 = continuous ? ng2 : p3 - 1;
+  const double* xg = GaussQuadrature::getCoord(ng1);
+  const double* yg = GaussQuadrature::getCoord(ng2);
+  const double* zg = GaussQuadrature::getCoord(ng3);
+  const double* wg = continuous ? GaussQuadrature::getWeight(ng1) : 0;
+  if (!xg || !yg || !zg) return false;
+  if (continuous && !wg) return false;
+
+  // Compute parameter values of the Gauss points over the whole patch
+  Matrix gp;
+  std::array<RealArray,3> gpar;
+  gpar[0] = this->getGaussPointParameters(gp,0,ng1,xg);
+  gpar[1] = this->getGaussPointParameters(gp,1,ng2,yg);
+  gpar[2] = this->getGaussPointParameters(gp,2,ng3,zg);
+
+  // Evaluate the secondary solution at all integration points
+  Matrix sField;
+  if (!this->evalSolution(sField,integrand,gpar.data()))
+  {
+    std::cerr <<" *** ASMs3DLag::assembleL2matrices: Failed for patch "<< idx+1
+              <<" nPoints="<< gpar[0].size()*gpar[1].size()*gpar[2].size()
+              << std::endl;
+    return false;
+  }
+
+  double dV = 1.0;
+  Vector phi(p1*p2*p3);
+  Matrix dNdu, Xnod, J;
+
+
+  // === Assembly loop over all elements in the patch ==========================
+
+  int iel = 0;
+  for (int i3 = 0; i3 < nel3; i3++)
+    for (int i2 = 0; i2 < nel2; i2++)
+      for (int i1 = 0; i1 < nel1; i1++, iel++)
+      {
+        if (MLGE[iel] < 1) continue; // zero-volume element
+
+        if (continuous)
+        {
+          // Set up control point (nodal) coordinates for current element
+          if (!this->getElementCoordinates(Xnod,1+iel))
+            return false;
+          else if ((dV = 0.125*this->getParametricVolume(1+iel)) < 0.0)
+            return false; // topology error (probably logic error)
+        }
+
+        int ip = ((i3*ng2*nel2 + i2)*ng1*nel1 + i1)*ng3;
+        const IntVec& mnpc = MNPC[iel];
+
+        // --- Integration loop over all Gauss points in each direction --------
+
+        Matrix eA(p1*p2*p3, p1*p2*p3);
+        Vectors eB(sField.rows(), Vector(p1*p2*p3));
+        for (int k = 0; k < ng3; k++)
+          for (int j = 0; j < ng2; j++)
+            for (int i = 0; i < ng1; i++, ip++)
+            {
+              // Compute basis function derivatives at current point
+              // using tensor product of one-dimensional Lagrange polynomials
+              if (continuous) {
+                if (!Lagrange::computeBasis(phi, dNdu,
+                                            p1,xg[i],p2,xg[j],p3,xg[k]))
+                  return false;
+              } else {
+                if (!Lagrange::computeBasis(phi,
+                                            p1,xg[i],p2,xg[j],p3,xg[k]))
+                  return false;
+              }
+
+              // Compute the Jacobian inverse and derivatives
+              double dJw = dV;
+              if (continuous)
+              {
+                dJw *= wg[i]*wg[j]*wg[k]*utl::Jacobian(J,dNdu,Xnod,dNdu,false);
+                if (dJw == 0.0) continue; // skip singular points
+              }
+
+              // Integrate the mass matrix
+              eA.outer_product(phi, phi, true, dJw);
+
+              // Integrate the rhs vector B
+              for (size_t r = 1; r <= sField.rows(); r++)
+                eB[r-1].add(phi,sField(r,ip+1)*dJw);
+            }
+
+        for (int i = 0; i < p1*p2*p3; ++i) {
+          for (int j = 0; j < p1*p2*p3; ++j)
+            A(mnpc[i]+1, mnpc[j]+1) += eA(i+1, j+1);
+
+          int jp = mnpc[i]+1;
+          for (size_t r = 0; r < sField.rows(); r++, jp += nnod)
+            B(jp) += eB[r](1+i);
+        }
+      }
+
+  return true;
+}
+*/

src/ASM/ASMs3Dmx.C

@@ -11,6 +11,7 @@
 //!
 //==============================================================================
 
+#include "GoTools/geometry/ObjectHeader.h"
 #include "GoTools/trivariate/SplineVolume.h"
 #include "GoTools/trivariate/VolumeInterpolator.h"
 
@@ -71,6 +72,30 @@ Go::SplineSurface* ASMs3Dmx::getBoundary (int dir, int basis)
 }
 
 
+bool ASMs3Dmx::read (std::istream& is, int basis)
+{
+  if (basis == 0)
+    return this->ASMs3D::read(is);
+
+  if (basis < 0 || basis > static_cast<int>(nfx.size()))
+    return false;
+
+  if (m_basis.empty()) {
+    m_basis.resize(nfx.size());
+    nb.resize(nfx.size(), 0);
+  }
+
+  Go::ObjectHeader head;
+  m_basis[basis-1] = std::make_shared<Go::SplineVolume>();
+  is >> head >> *m_basis[basis-1];
+  nb[basis-1] = m_basis[basis-1]->numCoefs(0)*
+                m_basis[basis-1]->numCoefs(1)*
+                m_basis[basis-1]->numCoefs(2);
+
+  return true;
+}
+
+
 bool ASMs3Dmx::write (std::ostream& os, int basis) const
 {
   if (basis == -1)

src/ASM/ASMs3Dmx.h

@@ -70,6 +70,8 @@ public:
   //! \param[in] inod 1-based node index local to current patch
   virtual Vec3 getCoord(size_t inod) const;
 
+  //! \brief Creates an instance by reading the given input stream.
+  virtual bool read(std::istream&, int basis = 0);
   //! \brief Writes the geometry/basis of the patch to given stream.
   virtual bool write(std::ostream& os, int basis) const;
 

src/ASM/ASMsupel.C

@@ -18,7 +18,7 @@
 #include <numeric>
 
 
-bool ASMsupel::read (std::istream& is)
+bool ASMsupel::read (std::istream& is, int)
 {
   // Lambda function for reading the supernode coordinates.
   auto&& readCoord = [&is](Vec3Vec& Xsup)

src/ASM/ASMsupel.h

@@ -39,7 +39,7 @@ public:
   virtual ~ASMsupel() {}
 
   //! \brief Creates an instance by reading the given input stream.
-  virtual bool read(std::istream& is);
+  virtual bool read(std::istream& is, int = 0);
   //! \brief Writes the geometry of the patch to the given stream.
   virtual bool write(std::ostream&, int) const;
   //! \brief Generates the finite element topology data for this patch.

src/ASM/ASMu1DLag.C

@@ -37,7 +37,7 @@ ASMu1DLag::ASMu1DLag (const ASMu1DLag& p) :
 }
 
 
-bool ASMu1DLag::read (std::istream& is)
+bool ASMu1DLag::read (std::istream& is, int)
 {
   switch (fileType) {
   case 'm':

src/ASM/ASMu1DLag.h

@@ -41,7 +41,7 @@ public:
   // ============================
 
   //! \brief Creates an instance by reading the given input stream.
-  virtual bool read(std::istream& is);
+  virtual bool read(std::istream& is, int = 0);
   //! \brief Generates a beam finite element model for the patch.
   //! \param[in] Zaxis Vector defining a point in the local XZ-plane
   virtual bool generateOrientedFEModel(const Vec3& Zaxis);

src/ASM/ASMu2DLag.C

@@ -37,7 +37,7 @@ ASMu2DLag::ASMu2DLag (const ASMu2DLag& p) :
 }
 
 
-bool ASMu2DLag::read (std::istream& is)
+bool ASMu2DLag::read (std::istream& is, int)
 {
   switch (fileType) {
   case 'm':

src/ASM/ASMu2DLag.h

@@ -41,7 +41,7 @@ public:
   // ============================
 
   //! \brief Creates an instance by reading the given input stream.
-  virtual bool read(std::istream& is);
+  virtual bool read(std::istream& is, int = 0);
   //! \brief Generates the finite element topology data for the patch.
   virtual bool generateFEMTopology();
   //! \brief Checks if this patch is empty.

src/ASM/LR/ASMu2D.C

@@ -64,7 +64,7 @@ ASMu2D::ASMu2D (const ASMu2D& patch, unsigned char n_f)
 }
 
 
-bool ASMu2D::read (std::istream& is)
+bool ASMu2D::read (std::istream& is, int)
 {
   if (shareFE) return true;
 

src/ASM/LR/ASMu2D.h

@@ -103,9 +103,9 @@ public:
   // ===========================================
 
   //! \brief Creates an instance by reading the given input stream.
-  virtual bool read(std::istream&);
+  virtual bool read(std::istream&, int = 0);
   //! \brief Writes the geometry of the SplineSurface object to given stream.
-  virtual bool write(std::ostream&, int = 0) const;
+  virtual bool write(std::ostream&, int) const;
 
   //! \brief Generates the finite element topology data for the patch.
   //! \details The data generated are the element-to-node connectivity array,

src/ASM/LR/ASMu2Dmx.C

@@ -72,6 +72,28 @@ LR::LRSplineSurface* ASMu2Dmx::getBasis (int basis)
 }
 
 
+bool ASMu2Dmx::read (std::istream& is, int basis)
+{
+  if (basis == 0)
+    return this->ASMu2D::read(is);
+
+  if (basis < 0 || basis > static_cast<int>(nfx.size()))
+    return false;
+
+  if (m_basis.empty()) {
+    m_basis.resize(nfx.size());
+    nb.resize(nfx.size(), 0);
+  }
+
+  m_basis[basis-1] = std::make_shared<LR::LRSplineSurface>();
+  is >> *m_basis[basis-1];
+  nb[basis-1] = m_basis[basis-1]->nBasisFunctions();
+  m_basis[basis-1]->generateIDs();
+
+  return true;
+}
+
+
 bool ASMu2Dmx::write (std::ostream& os, int basis) const
 {
   if (basis == -1)

src/ASM/LR/ASMu2Dmx.h

@@ -47,6 +47,8 @@ public:
   // Methods for model generation
   // ============================
 
+  //! \brief Creates an instance by reading the given input stream.
+  virtual bool read(std::istream&, int basis);
   //! \brief Writes the geometry/basis of the patch to given stream.
   virtual bool write(std::ostream& os, int basis) const;
 

src/ASM/LR/ASMu2Dnurbs.C

@@ -37,7 +37,7 @@ ASMu2Dnurbs::ASMu2Dnurbs (const ASMu2Dnurbs& patch, unsigned char n_f)
 }
 
 
-bool ASMu2Dnurbs::read (std::istream& is)
+bool ASMu2Dnurbs::read (std::istream& is, int)
 {
   bool ok = this->ASMu2DC1::read(is);
   if (ok && !(tensorspline && tensorspline->rational()))

src/ASM/LR/ASMu2Dnurbs.h

@@ -33,7 +33,7 @@ public:
   virtual ~ASMu2Dnurbs() {}
 
   //! \brief Creates an instance by reading the given input stream.
-  virtual bool read(std::istream& is);
+  virtual bool read(std::istream& is, int = 0);
 
 protected:
   //! \brief Evaluates the basis functions and derivatives of an element.

src/ASM/LR/ASMu3D.C

@@ -63,7 +63,7 @@ ASMu3D::ASMu3D (const ASMu3D& patch, unsigned char n_f)
 }
 
 
-bool ASMu3D::read (std::istream& is)
+bool ASMu3D::read (std::istream& is, int)
 {
   if (shareFE) return true;
 

src/ASM/LR/ASMu3D.h

@@ -65,9 +65,9 @@ public:
   // ===========================================
 
   //! \brief Creates an instance by reading the given input stream.
-  virtual bool read(std::istream&);
+  virtual bool read(std::istream&, int = 0);
   //! \brief Writes the geometry of the SplineVolume object to given stream.
-  virtual bool write(std::ostream&, int = 0) const;
+  virtual bool write(std::ostream&, int) const;
 
   //! \brief Generates the finite element topology data for the patch.
   //! \details The data generated are the element-to-node connectivity array,

src/ASM/LR/ASMu3Dmx.C

@@ -73,6 +73,28 @@ LR::LRSplineVolume* ASMu3Dmx::getBasis (int basis)
 }
 
 
+bool ASMu3Dmx::read (std::istream& is, int basis)
+{
+  if (basis == 0)
+    return this->ASMu3D::read(is,0);
+
+  if (basis < 0 || basis > static_cast<int>(nfx.size()))
+    return false;
+
+  if (m_basis.empty()) {
+    m_basis.resize(nfx.size());
+    nb.resize(nfx.size(), 0);
+  }
+
+  m_basis[basis-1] = std::make_shared<LR::LRSplineVolume>();
+  is >> *m_basis[basis-1];
+  nb[basis-1] = m_basis[basis-1]->nBasisFunctions();
+  m_basis[basis-1]->generateIDs();
+
+  return true;
+}
+
+
 bool ASMu3Dmx::write (std::ostream& os, int basis) const
 {
   if (basis == -1)

src/ASM/LR/ASMu3Dmx.h

@@ -47,6 +47,8 @@ public:
   // Methods for model generation
   // ============================
 
+  //! \brief Creates an instance by reading the given input stream.
+  virtual bool read(std::istream&, int basis);
   //! \brief Writes the geometry/basis of the patch to given stream.
   virtual bool write(std::ostream& os, int basis) const;