diff --git a/src/ASM/ASMbase.h b/src/ASM/ASMbase.h index 382bfd34..9aa806cf 100644 --- a/src/ASM/ASMbase.h +++ b/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; diff --git a/src/ASM/ASMs1D.C b/src/ASM/ASMs1D.C index 74d5448a..12e494e7 100644 --- a/src/ASM/ASMs1D.C +++ b/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; diff --git a/src/ASM/ASMs1D.h b/src/ASM/ASMs1D.h index 30114f66..101ec65a 100644 --- a/src/ASM/ASMs1D.h +++ b/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; diff --git a/src/ASM/ASMs2D.C b/src/ASM/ASMs2D.C index b2b14eab..011d9211 100644 --- a/src/ASM/ASMs2D.C +++ b/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; diff --git a/src/ASM/ASMs2D.h b/src/ASM/ASMs2D.h index d7536582..5ec7014f 100644 --- a/src/ASM/ASMs2D.h +++ b/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, diff --git a/src/ASM/ASMs2Dmx.C b/src/ASM/ASMs2Dmx.C index 428d94ba..51617d23 100644 --- a/src/ASM/ASMs2Dmx.C +++ b/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(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(); + 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) diff --git a/src/ASM/ASMs2Dmx.h b/src/ASM/ASMs2Dmx.h index cb52b0f7..ca93504b 100644 --- a/src/ASM/ASMs2Dmx.h +++ b/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; diff --git a/src/ASM/ASMs3D.C b/src/ASM/ASMs3D.C index b1ee3dfa..a3d6a089 100644 --- a/src/ASM/ASMs3D.C +++ b/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; diff --git a/src/ASM/ASMs3D.h b/src/ASM/ASMs3D.h index 0e178646..339a9b6e 100644 --- a/src/ASM/ASMs3D.h +++ b/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, diff --git a/src/ASM/ASMs3DLagrecovery.C b/src/ASM/ASMs3DLagrecovery.C new file mode 100644 index 00000000..e2e96f64 --- /dev/null +++ b/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 + + +/* +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 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; +} +*/ diff --git a/src/ASM/ASMs3Dmx.C b/src/ASM/ASMs3Dmx.C index e5bfd133..d754fd7c 100644 --- a/src/ASM/ASMs3Dmx.C +++ b/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(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(); + 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) diff --git a/src/ASM/ASMs3Dmx.h b/src/ASM/ASMs3Dmx.h index 2aa5d40c..9abb5ee2 100644 --- a/src/ASM/ASMs3Dmx.h +++ b/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; diff --git a/src/ASM/ASMsupel.C b/src/ASM/ASMsupel.C index f45ab111..615afe9a 100644 --- a/src/ASM/ASMsupel.C +++ b/src/ASM/ASMsupel.C @@ -18,7 +18,7 @@ #include -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) diff --git a/src/ASM/ASMsupel.h b/src/ASM/ASMsupel.h index 08386859..7cf40496 100644 --- a/src/ASM/ASMsupel.h +++ b/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. diff --git a/src/ASM/ASMu1DLag.C b/src/ASM/ASMu1DLag.C index f999ee12..6e1cb8f9 100644 --- a/src/ASM/ASMu1DLag.C +++ b/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': diff --git a/src/ASM/ASMu1DLag.h b/src/ASM/ASMu1DLag.h index 38f0677f..19ef629b 100644 --- a/src/ASM/ASMu1DLag.h +++ b/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); diff --git a/src/ASM/ASMu2DLag.C b/src/ASM/ASMu2DLag.C index 29af5f51..6a290540 100644 --- a/src/ASM/ASMu2DLag.C +++ b/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': diff --git a/src/ASM/ASMu2DLag.h b/src/ASM/ASMu2DLag.h index ccf58ed1..5d626e18 100644 --- a/src/ASM/ASMu2DLag.h +++ b/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. diff --git a/src/ASM/LR/ASMu2D.C b/src/ASM/LR/ASMu2D.C index 569bc772..10bca1f0 100644 --- a/src/ASM/LR/ASMu2D.C +++ b/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; diff --git a/src/ASM/LR/ASMu2D.h b/src/ASM/LR/ASMu2D.h index bc57f0dd..008f72ed 100644 --- a/src/ASM/LR/ASMu2D.h +++ b/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, diff --git a/src/ASM/LR/ASMu2Dmx.C b/src/ASM/LR/ASMu2Dmx.C index 4aa8c8f9..4b6ce18d 100644 --- a/src/ASM/LR/ASMu2Dmx.C +++ b/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(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(); + 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) diff --git a/src/ASM/LR/ASMu2Dmx.h b/src/ASM/LR/ASMu2Dmx.h index 7baeaadf..098eb3d5 100644 --- a/src/ASM/LR/ASMu2Dmx.h +++ b/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; diff --git a/src/ASM/LR/ASMu2Dnurbs.C b/src/ASM/LR/ASMu2Dnurbs.C index b87f84ff..eb51dff8 100644 --- a/src/ASM/LR/ASMu2Dnurbs.C +++ b/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())) diff --git a/src/ASM/LR/ASMu2Dnurbs.h b/src/ASM/LR/ASMu2Dnurbs.h index beb9355b..fa85ea87 100644 --- a/src/ASM/LR/ASMu2Dnurbs.h +++ b/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. diff --git a/src/ASM/LR/ASMu3D.C b/src/ASM/LR/ASMu3D.C index 8ca7539c..9a51175d 100644 --- a/src/ASM/LR/ASMu3D.C +++ b/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; diff --git a/src/ASM/LR/ASMu3D.h b/src/ASM/LR/ASMu3D.h index 081127e3..d267da85 100644 --- a/src/ASM/LR/ASMu3D.h +++ b/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, diff --git a/src/ASM/LR/ASMu3Dmx.C b/src/ASM/LR/ASMu3Dmx.C index 82f577da..c6c603b8 100644 --- a/src/ASM/LR/ASMu3Dmx.C +++ b/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(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(); + 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) diff --git a/src/ASM/LR/ASMu3Dmx.h b/src/ASM/LR/ASMu3Dmx.h index 75dcc6c1..699075b7 100644 --- a/src/ASM/LR/ASMu3Dmx.h +++ b/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;