From 6b3954c502303bb9a7ccc426c00d76c45a79be96 Mon Sep 17 00:00:00 2001 From: Timo van Opstal Date: Thu, 17 Dec 2015 16:30:51 +0100 Subject: [PATCH] added: COMPATIBLE integrand --- Apps/Common/WeakOperators.C | 18 ++++++++++-------- Apps/Common/WeakOperators.h | 11 ++++++++--- src/ASM/ASMbase.h | 11 ++++++----- src/ASM/ASMs2Dmx.C | 14 ++++++++++---- src/ASM/ASMs2DmxLag.C | 2 +- src/ASM/ASMs3Dmx.C | 13 +++++++------ src/ASM/ASMs3DmxLag.C | 2 +- src/ASM/GlbForceVec.C | 6 +++--- src/ASM/GlbL2projector.C | 8 +++++--- src/ASM/IntegrandBase.C | 3 ++- src/ASM/IntegrandBase.h | 6 +++++- src/SIM/ForceIntegrator.C | 11 +++++++++++ src/SIM/SIMbase.C | 9 +++++---- src/Utility/HDF5Writer.C | 1 + 14 files changed, 75 insertions(+), 40 deletions(-) diff --git a/Apps/Common/WeakOperators.C b/Apps/Common/WeakOperators.C index 2182e900..e127eeae 100644 --- a/Apps/Common/WeakOperators.C +++ b/Apps/Common/WeakOperators.C @@ -140,10 +140,10 @@ namespace WeakOperators { void Laplacian(Matrix& EM, const FiniteElement& fe, double scale, size_t cmp, size_t nf, - bool stress, size_t scmp) + bool stress, size_t scmp, unsigned char basis) { Matrix A; - A.multiply(fe.dNdX,fe.dNdX,false,true); + A.multiply(fe.grad(basis),fe.grad(basis),false,true); A *= scale*fe.detJxW; addComponents(EM, A, cmp, nf, scmp); if (stress) @@ -151,7 +151,7 @@ namespace WeakOperators { for (size_t j = 1; j <= fe.N.size(); j++) for (size_t k = 1; k <= cmp; k++) for (size_t l = 1; l <= cmp; l++) - EM(nf*(j-1)+k+scmp,nf*(i-1)+l+scmp) += scale*fe.dNdX(i,k)*fe.dNdX(j,l)*fe.detJxW; + EM(nf*(j-1)+k+scmp,nf*(i-1)+l+scmp) += scale*fe.grad(basis)(i,k)*fe.grad(basis)(j,l)*fe.detJxW; } @@ -166,24 +166,26 @@ namespace WeakOperators { void Mass(Matrix& EM, const FiniteElement& fe, - double scale, size_t cmp, size_t nf, size_t scmp) + double scale, size_t cmp, size_t nf, size_t scmp, + unsigned char basis) { Matrix A; - A.outer_product(fe.N,fe.N,false); + A.outer_product(fe.basis(basis),fe.basis(basis),false); A *= scale*fe.detJxW; addComponents(EM, A, cmp, nf, scmp); } void Source(Vector& EV, const FiniteElement& fe, - double scale, size_t cmp, size_t nf, size_t scmp) + double scale, size_t cmp, size_t nf, size_t scmp, + unsigned char basis) { if (cmp == 1 && nf == 1) - EV.add(fe.N, scale*fe.detJxW); + EV.add(fe.basis(basis), scale*fe.detJxW); else { for (size_t i = 1; i <= fe.N.size(); ++i) for (size_t k = 1; k <= cmp; ++k) - EV(nf*(i-1)+k+scmp) += scale*fe.N(i)*fe.detJxW; + EV(nf*(i-1)+k+scmp) += scale*fe.basis(basis)(i)*fe.detJxW; } } diff --git a/Apps/Common/WeakOperators.h b/Apps/Common/WeakOperators.h index 2a844d5e..fc16f309 100644 --- a/Apps/Common/WeakOperators.h +++ b/Apps/Common/WeakOperators.h @@ -95,9 +95,10 @@ namespace WeakOperators //! \param[in] nf Number of fields in basis. //! \param[in] stress Whether to add extra stress formulation terms. //! \param[in] scmp Starting component. + //! \param[in] basis Basis to use. void Laplacian(Matrix& EM, const FiniteElement& fe, double scale=1.0, size_t cmp=1, size_t nf=1, - bool stress=false, size_t scmp=0); + bool stress=false, size_t scmp=0, unsigned char basis=1); //! \brief Compute a heteregenous coefficient laplacian. //! \param[out] EM The element matrix to add contribution to. @@ -114,8 +115,10 @@ namespace WeakOperators //! \param[in] cmp Number of components to add. //! \param[in] nf Number of fields in basis. //! \param[in] scmp Starting component. + //! \param[in] basis Basis to use. void Mass(Matrix& EM, const FiniteElement& fe, - double scale=1.0, size_t cmp=1, size_t nf=1, size_t scmp=0); + double scale=1.0, size_t cmp=1, size_t nf=1, size_t scmp=0, + unsigned char basis=1); //! \brief Compute a source term. //! \param[out] EV The element vector to add contribution to. @@ -124,8 +127,10 @@ namespace WeakOperators //! \param[in] cmp Number of components to add. //! \param[in] nf Number of fields in basis. //! \param[in] scmp Starting component. + //! \param[in] basis Basis to use. void Source(Vector& EV, const FiniteElement& fe, - double scale=1.0, size_t cmp=1, size_t nf=1, size_t scmp=0); + double scale=1.0, size_t cmp=1, size_t nf=1, size_t scmp=0, + unsigned char basis=1); //! \brief Compute a vector-source term. //! \param[out] EV The element vector to add contribution to. diff --git a/src/ASM/ASMbase.h b/src/ASM/ASMbase.h index f1cac1c6..67641b74 100644 --- a/src/ASM/ASMbase.h +++ b/src/ASM/ASMbase.h @@ -556,16 +556,17 @@ public: //! \param[in] code Identifier for inhomogeneous Dirichlet condition field bool add2PC(int slave, int dir, int master, int code = 0); -protected: - - // Internal methods for preprocessing of boundary conditions - // ========================================================= - //! \brief Adds a general multi-point-constraint (MPC) equation to this patch. //! \param mpc Pointer to an MPC-object //! \param[in] code Identifier for inhomogeneous Dirichlet condition field //! \param[in] silence If \e true, suppress debug print bool addMPC(MPC*& mpc, int code = 0, bool silence = false); + +protected: + + // Internal methods for preprocessing of boundary conditions + // ========================================================= + //! \brief Creates and adds a three-point constraint to this patch. //! \param[in] slave Global node number of the node to constrain //! \param[in] dir Which local DOF to constrain (1, 2, 3) diff --git a/src/ASM/ASMs2Dmx.C b/src/ASM/ASMs2Dmx.C index 91a86346..19cc5198 100644 --- a/src/ASM/ASMs2Dmx.C +++ b/src/ASM/ASMs2Dmx.C @@ -929,18 +929,24 @@ bool ASMs2Dmx::evalSolution (Matrix& sField, const IntegrandBase& integrand, // Fetch indices of the non-zero basis functions at this point std::vector ip(m_basis.size()); IntVec ipa; + size_t ofs = 0; for (size_t b = 0; b < m_basis.size(); ++b) { scatterInd(m_basis[b]->numCoefs_u(),m_basis[b]->numCoefs_v(), m_basis[b]->order_u(),m_basis[b]->order_v(),splinex[b][i].left_idx,ip[b]); + + // Fetch associated control point coordinates + if (b == (size_t)geoBasis-1) + utl::gather(ip[geoBasis-1], nsd, Xnod, Xtmp); + + for (auto& it : ip[b]) + it += ofs; ipa.insert(ipa.end(), ip[b].begin(), ip[b].end()); + ofs += nb[b]; } fe.u = splinex[0][i].param[0]; fe.v = splinex[0][i].param[1]; - // Fetch associated control point coordinates - utl::gather(ip[geoBasis-1], nsd, Xnod, Xtmp); - // 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),dNxdu[b]); @@ -962,7 +968,7 @@ bool ASMs2Dmx::evalSolution (Matrix& sField, const IntegrandBase& integrand, X = Xtmp * fe.basis(geoBasis); // Now evaluate the solution field - if (!integrand.evalSol(solPt,fe,X,ipa,elem_sizes)) + if (!integrand.evalSol(solPt,fe,X,ipa,elem_sizes,nb)) return false; else if (sField.empty()) sField.resize(solPt.size(),nPoints,true); diff --git a/src/ASM/ASMs2DmxLag.C b/src/ASM/ASMs2DmxLag.C index 3d528083..6107b370 100644 --- a/src/ASM/ASMs2DmxLag.C +++ b/src/ASM/ASMs2DmxLag.C @@ -552,7 +552,7 @@ bool ASMs2DmxLag::evalSolution (Matrix& sField, const IntegrandBase& integrand, } // Now evaluate the solution field - if (!integrand.evalSol(solPt,fe,Xnod*fe.basis(geoBasis),MNPC[iel-1],elem_size)) + if (!integrand.evalSol(solPt,fe,Xnod*fe.basis(geoBasis),MNPC[iel-1],elem_size,nb)) return false; else if (sField.empty()) sField.resize(solPt.size(),nPoints,true); diff --git a/src/ASM/ASMs3Dmx.C b/src/ASM/ASMs3Dmx.C index ef983cb4..859907c6 100644 --- a/src/ASM/ASMs3Dmx.C +++ b/src/ASM/ASMs3Dmx.C @@ -286,7 +286,8 @@ bool ASMs3Dmx::generateFEMTopology () iel = 0; if ((int)b != geoBasis-1) addBasis(m_basis[b],false); - inod += nb[b]; + else + inod += nb[b]; lnod2 += m_basis[b]->order(0)*m_basis[b]->order(1)*m_basis[b]->order(2); } @@ -391,10 +392,10 @@ Vec3 ASMs3Dmx::getCoord (size_t inod) const const int J = nodeInd[inod-1].J; const int K = nodeInd[inod-1].K; - int b = 0; - size_t nbb = 0; - while (inod < nbb) - nbb += nb[b++]; + size_t b = 0; + size_t nbb = nb[0]; + while (nbb < inod) + nbb += nb[++b]; cit = m_basis[b]->coefs_begin() + ((K*m_basis[b]->numCoefs(1)+J)*m_basis[b]->numCoefs(0)+I) * m_basis[b]->dimension(); @@ -930,7 +931,7 @@ bool ASMs3Dmx::evalSolution (Matrix& sField, const IntegrandBase& integrand, X = Xtmp * fe.basis(geoBasis); // Now evaluate the solution field - if (!integrand.evalSol(solPt,fe,X,ipa,elem_sizes)) + if (!integrand.evalSol(solPt,fe,X,ipa,elem_sizes,nb)) return false; else if (sField.empty()) sField.resize(solPt.size(),nPoints,true); diff --git a/src/ASM/ASMs3DmxLag.C b/src/ASM/ASMs3DmxLag.C index f5cd905e..19875892 100644 --- a/src/ASM/ASMs3DmxLag.C +++ b/src/ASM/ASMs3DmxLag.C @@ -622,7 +622,7 @@ bool ASMs3DmxLag::evalSolution (Matrix& sField, const IntegrandBase& integrand, } // Now evaluate the solution field - if (!integrand.evalSol(solPt,fe,Xnod*fe.basis(geoBasis),MNPC[iel-1],elem_size)) + if (!integrand.evalSol(solPt,fe,Xnod*fe.basis(geoBasis),MNPC[iel-1],elem_size,nb)) return false; else if (sField.empty()) sField.resize(solPt.size(),nPoints,true); diff --git a/src/ASM/GlbForceVec.C b/src/ASM/GlbForceVec.C index 6c573f63..0c0dc85e 100644 --- a/src/ASM/GlbForceVec.C +++ b/src/ASM/GlbForceVec.C @@ -53,19 +53,19 @@ bool GlbForceVec::assemble (const LocalIntegral* elmObj, int elmId) const ElmMats* elm = dynamic_cast(elmObj); if (!elm || elm->b.size() < 1) return false; - const Vector& ES = elm->b.front(); + const Vector& ES = elm->getRHSVector(); const size_t nfc = F.rows(); std::vector mnpc; if (!sam.getElmNodes(mnpc,elmId)) return false; - else if (ES.size() < nfc*mnpc.size()) +/* else if (ES.size() < nfc*mnpc.size()) { std::cerr <<" *** GlbForceVec::assemble: Invalid element force vector," <<" size="<< ES.size() <<" should be (at least) " << nfc*mnpc.size() << std::endl; return false; } - +*/ // Assemble the nodal forces into the Matrix F size_t i, j, k, ninod = 0; std::map::const_iterator nit; diff --git a/src/ASM/GlbL2projector.C b/src/ASM/GlbL2projector.C index 7e420268..4806c955 100644 --- a/src/ASM/GlbL2projector.C +++ b/src/ASM/GlbL2projector.C @@ -40,7 +40,8 @@ public: GlbL2& gl2Int; //!< The global L2 projection integrand LocalIntegral* elmData; //!< Element data associated with problem integrand IntVec mnpc; //!< Matrix of element nodal correspondance for bases - std::vector sizes; //!< Size of each basis + std::vector elem_sizes; //!< Size of each basis on the element + std::vector basis_sizes; //!< Size of each basis on the patch }; @@ -85,7 +86,8 @@ bool GlbL2::initElement (const IntVec& MNPC1, L2Mats& gl2 = static_cast(elmInt); gl2.mnpc = MNPC1; - gl2.sizes = elem_sizes; + gl2.elem_sizes = elem_sizes; + gl2.basis_sizes = basis_sizes; return problem.initElement(MNPC1,elem_sizes,basis_sizes,*gl2.elmData); } @@ -127,7 +129,7 @@ bool GlbL2::evalIntMx (LocalIntegral& elmInt, L2Mats& gl2 = static_cast(elmInt); Vector solPt; - if (!problem.evalSol(solPt,fe,X,gl2.mnpc,gl2.sizes)) + if (!problem.evalSol(solPt,fe,X,gl2.mnpc,gl2.elem_sizes,gl2.basis_sizes)) if (!problem.diverged(fe.iGP+1)) return false; diff --git a/src/ASM/IntegrandBase.C b/src/ASM/IntegrandBase.C index f7723293..e544e92e 100644 --- a/src/ASM/IntegrandBase.C +++ b/src/ASM/IntegrandBase.C @@ -140,7 +140,8 @@ bool IntegrandBase::evalSol (Vector& s, const FiniteElement& fe, bool IntegrandBase::evalSol (Vector& s, const MxFiniteElement& fe, const Vec3& X, const std::vector& MNPC, - const std::vector& elem_sizes) const + const std::vector& elem_sizes, + const std::vector& basis_sizes) const { std::vector MNPC1(MNPC.begin(), MNPC.begin()+elem_sizes.front()); return this->evalSol(s,fe,X,MNPC1); diff --git a/src/ASM/IntegrandBase.h b/src/ASM/IntegrandBase.h index f10bce34..e856da90 100644 --- a/src/ASM/IntegrandBase.h +++ b/src/ASM/IntegrandBase.h @@ -59,6 +59,8 @@ public: //! \brief Defines the solution mode before the element assembly is started. virtual void setMode(SIM::SolutionMode mode) { m_mode = mode; } + //! \brief Obtain current solution mode + SIM::SolutionMode getMode() const { return m_mode; } //! \brief Initializes an integration parameter for the integrand. virtual void setIntegrationPrm(unsigned short int, double) {} //! \brief Returns an integration parameter for the integrand. @@ -158,9 +160,11 @@ public: //! \param[in] X Cartesian coordinates of current point //! \param[in] MNPC Nodal point correspondance for the bases //! \param[in] elem_sizes Size of each basis on the element + //! \param[in] basis_sizes Size of each basis on the patch virtual bool evalSol(Vector& s, const MxFiniteElement& fe, const Vec3& X, const std::vector& MNPC, - const std::vector& elem_sizes) const; + const std::vector& elem_sizes, + const std::vector& basis_sizes) const; //! \brief Evaluates the analytical secondary solution at a result point. //! \param[out] s The solution field values at current point diff --git a/src/SIM/ForceIntegrator.C b/src/SIM/ForceIntegrator.C index bdc3eee6..a1c68afb 100644 --- a/src/SIM/ForceIntegrator.C +++ b/src/SIM/ForceIntegrator.C @@ -103,6 +103,17 @@ bool SIM::integrate(const Vectors& solution, SIMbase* model, int code, size_t prevPatch = 0; GlobalIntegral dummy; GlobalIntegral& frc = force ? *force : dummy; + + if (code == 0) { // special case - volume integral over the entire model + for (int p = 1; p <= model->getNoPatches() && ok; ++p) { + ASMbase* patch = model->getPatch(p); + ok = model->extractPatchSolution(solution,p-1); + model->setPatchMaterial(p); + ok &= patch->integrate(*forceInt,frc,time); + } + return ok; + } + PropertyVec::const_iterator p; for (p = model->begin_prop(); p != model->end_prop() && ok; p++) if (abs(p->pindx) == code) diff --git a/src/SIM/SIMbase.C b/src/SIM/SIMbase.C index 8b1074dc..3b970e02 100644 --- a/src/SIM/SIMbase.C +++ b/src/SIM/SIMbase.C @@ -309,6 +309,8 @@ bool SIMbase::parseBCTag (const TiXmlElement* elem) std::string set, type; utl::getAttribute(elem,"set",set); utl::getAttribute(elem,"type",type,true); + int ndir = 0; + utl::getAttribute(elem,"direction",ndir); int code = this->getUniquePropertyCode(set); if (code == 0) utl::getAttribute(elem,"code",code); if (type == "anasol") { @@ -319,14 +321,12 @@ bool SIMbase::parseBCTag (const TiXmlElement* elem) IFEM::cout <<"\tNeumann code "<< code <<" (generic)"; if (elem->FirstChild()) { this->setPropertyType(code,Property::ROBIN); - this->setNeumann(elem->FirstChild()->Value(),"expression",0,code); + this->setNeumann(elem->FirstChild()->Value(),"expression",ndir,code); } else this->setPropertyType(code,Property::NEUMANN_GENERIC); } else { - int ndir = 0; - utl::getAttribute(elem,"direction",ndir); IFEM::cout <<"\tNeumann code "<< code <<" direction "<< ndir; if (!type.empty()) IFEM::cout <<" ("<< type <<")"; if (elem->FirstChild()) @@ -1044,7 +1044,8 @@ bool SIMbase::preprocess (const IntVec& ignored, bool fixDup) (*mit)->generateThreadGroups(*myProblem,silence); for (q = myProps.begin(); q != myProps.end(); q++) - if (q->pcode == Property::NEUMANN || q->pcode == Property::NEUMANN_GENERIC) + if (q->pcode == Property::NEUMANN || q->pcode == Property::NEUMANN_GENERIC || + q->pcode == Property::ROBIN) this->generateThreadGroups(*q,silence); // Preprocess the result points diff --git a/src/Utility/HDF5Writer.C b/src/Utility/HDF5Writer.C index b54bf7e2..6a6341dd 100644 --- a/src/Utility/HDF5Writer.C +++ b/src/Utility/HDF5Writer.C @@ -465,6 +465,7 @@ void HDF5Writer::writeSIM (int level, const DataEntry& entry, if (abs(results) & DataExporter::SECONDARY) { Matrix field; if (prefix.empty()) { + sim->setMode(SIM::RECOVERY); sim->extractPatchSolution(Vectors(1,*sol), loc-1); sim->evalSecondarySolution(field,loc-1); }