From 2ca3f3e48273f43f6ae644d662f1527201191b3d Mon Sep 17 00:00:00 2001 From: Knut Morten Okstad Date: Fri, 7 Sep 2018 07:58:15 +0200 Subject: [PATCH 1/2] Fixed: Copy constructor that is used when shared grid. Changed: Made the 2D and 3D mixed implementations equivalent. There were some discrepancies here and there which I think were not intended, but a result of doing non-tested things at different times. --- src/ASM/ASMs2Dmx.C | 170 ++++++++++++++------------ src/ASM/ASMs2Dmx.h | 39 +++--- src/ASM/ASMs3Dmx.C | 299 +++++++++++++++++++++++---------------------- src/ASM/ASMs3Dmx.h | 20 +-- 4 files changed, 274 insertions(+), 254 deletions(-) diff --git a/src/ASM/ASMs2Dmx.C b/src/ASM/ASMs2Dmx.C index 61ad7f360..721862d7f 100644 --- a/src/ASM/ASMs2Dmx.C +++ b/src/ASM/ASMs2Dmx.C @@ -27,7 +27,6 @@ #include "Utilities.h" #include "Profiler.h" #include "Vec3Oper.h" -#include "Vec3.h" #include #include #ifdef USE_OPENMP @@ -36,13 +35,13 @@ ASMs2Dmx::ASMs2Dmx (unsigned char n_s, const CharVec& n_f) - : ASMs2D(n_s, *std::max_element(n_f.begin(),n_f.end())), ASMmxBase(n_f) + : ASMs2D(n_s,std::accumulate(n_f.begin(),n_f.end(),0)), ASMmxBase(n_f) { } ASMs2Dmx::ASMs2Dmx (const ASMs2Dmx& patch, const CharVec& n_f) - : ASMs2D(patch), ASMmxBase(n_f) + : ASMs2D(patch), ASMmxBase(n_f), m_basis(patch.m_basis) { nb = patch.nb; } @@ -59,10 +58,12 @@ Go::SplineSurface* ASMs2Dmx::getBasis (int basis) const Go::SplineCurve* ASMs2Dmx::getBoundary (int dir, int basis) { - if (dir < -2 || dir == 0 || dir > 2 || basis < 1 || basis > (int)m_basis.size()) + if (basis < 1 || basis > (int)m_basis.size()) return nullptr; int iedge = dir > 0 ? dir : 3*dir+6; + if (iedge < 0 || iedge > 3) + return nullptr; return m_basis[basis-1]->edgeCurve(iedge); } @@ -71,11 +72,11 @@ Go::SplineCurve* ASMs2Dmx::getBoundary (int dir, int basis) bool ASMs2Dmx::write (std::ostream& os, int basis) const { if (basis == -1) - os <<"200 1 0 0\n" << *proj; + os <<"200 1 0 0\n"<< *proj; else if (basis < 1 || basis > (int)m_basis.size()) - os <<"200 1 0 0\n" << *surf; + os <<"200 1 0 0\n"<< *surf; else if (m_basis[basis-1]) - os <<"200 1 0 0\n" << *m_basis[basis-1]; + os <<"200 1 0 0\n"<< *m_basis[basis-1]; else return false; @@ -105,10 +106,7 @@ size_t ASMs2Dmx::getNoNodes (int basis) const unsigned char ASMs2Dmx::getNoFields (int basis) const { - if (basis < 1 || basis > (int)nfx.size()) - return std::accumulate(nfx.begin(), nfx.end(), 0); - - return nfx[basis-1]; + return basis > 0 && basis <= (int)nfx.size() ? nfx[basis-1] : nf; } @@ -117,12 +115,12 @@ unsigned char ASMs2Dmx::getNodalDOFs (size_t inod) const if (this->isLMn(inod)) return nLag; - size_t nbc=0; - for (size_t i=0;iisLMn(inod)) return 'L'; - size_t nbc=nb.front(); + size_t nbc = nb.front(); if (inod <= nbc) return 'D'; - else for (size_t i = 1; i < nb.size(); i++) + + for (size_t i = 1; i < nb.size(); i++) if (inod <= (nbc += nb[i])) return 'O'+i; @@ -156,7 +155,7 @@ void ASMs2Dmx::extractNodeVec (const Vector& globRes, Vector& nodeVec, bool ASMs2Dmx::injectNodeVec (const Vector& nodeRes, Vector& globRes, - unsigned char, int basis) const + unsigned char, int basis) const { this->injectNodeVecMx(globRes,nodeRes,basis); return true; @@ -196,24 +195,25 @@ bool ASMs2Dmx::generateFEMTopology () elem_size.push_back(it->order_u()*it->order_v()); } + nnod = std::accumulate(nb.begin(), nb.end(), 0u); if (!nodeInd.empty() && !shareFE) { - if (nodeInd.size() == std::accumulate(nb.begin(), nb.end(), 0u)) return true; + if (nodeInd.size() == nnod) + return true; + std::cerr <<" *** ASMs2Dmx::generateFEMTopology: Inconsistency between the" <<" number of FE nodes "<< nodeInd.size() - <<"\n and the number of spline coefficients " - << std::accumulate(nb.begin(), nb.end(), 0) + <<"\n and the number of spline coefficients "<< nnod <<" in the patch."<< std::endl; return false; } - - if (shareFE == 'F') return true; + else if (shareFE == 'F') + return true; int i1, i2, j1, j2; #ifdef SP_DEBUG size_t nbasis=0; for (auto& it : m_basis) { - int i1, i2; std::cout << "Basis " << ++nbasis << ":\n"; std::cout <<"numCoefs: "<< it->numCoefs_u() <<" "<< it->numCoefs_v(); std::cout <<"\norder: "<< it->order_u() <<" "<< it->order_v(); @@ -229,15 +229,14 @@ bool ASMs2Dmx::generateFEMTopology () nel = (m_basis[geoBasis-1]->numCoefs_u()-m_basis[geoBasis-1]->order_u()+1)* (m_basis[geoBasis-1]->numCoefs_v()-m_basis[geoBasis-1]->order_v()+1); - nnod = std::accumulate(nb.begin(), nb.end(), 0); myMLGE.resize(nel,0); myMLGN.resize(nnod); myMNPC.resize(nel); myNodeInd.resize(nnod); - size_t iel, inod = 0; - for (auto& it : m_basis) { + size_t iel = 0, inod = 0; + for (auto& it : m_basis) for (i2 = 0; i2 < it->numCoefs_v(); i2++) for (i1 = 0; i1 < it->numCoefs_u(); i1++) { @@ -245,9 +244,6 @@ bool ASMs2Dmx::generateFEMTopology () myNodeInd[inod].J = i2; myMLGN[inod++] = ++gNod; } - } - - iel = 0, inod = std::accumulate(nb.begin(),nb.begin()+geoBasis-1,0u); int lnod2 = 0; int lnod3 = 0; @@ -257,6 +253,7 @@ bool ASMs2Dmx::generateFEMTopology () lnod3 += m_basis[i2]->order_u()*m_basis[i2]->order_v(); // Create nodal connectivities for bases + inod = std::accumulate(nb.begin(),nb.begin()+geoBasis-1,0u); Go::SplineSurface* b = m_basis[geoBasis-1].get(); auto knotv = b->basis_v().begin(); for (i2 = 1; i2 <= b->numCoefs_v(); i2++, ++knotv) { @@ -279,7 +276,7 @@ bool ASMs2Dmx::generateFEMTopology () // find knotspan spanning element for other bases lnod = 0; size_t lnod4 = 0; - for (size_t bas = 0; bas < m_basis.size(); ++bas) { + for (size_t bas = 0; bas < m_basis.size(); ++bas) { if (bas != (size_t)geoBasis-1) { double ku = *knotu; double kv = *knotv; @@ -313,9 +310,10 @@ bool ASMs2Dmx::connectPatch (int edge, ASM2D& neighbor, int nedge, bool revers, if (!neighMx) return false; size_t nb1 = 0, nb2 = 0; - for (size_t i = 1; i <= m_basis.size(); ++i) { + for (size_t i = 1; i <= nb.size(); i++) { if (basis == 0 || i == (size_t)basis) - if (!this->connectBasis(edge,*neighMx,nedge,revers,i,nb1,nb2,coordCheck,thick)) + if (!this->connectBasis(edge,*neighMx,nedge,revers,i,nb1,nb2, + coordCheck,thick)) return false; nb1 += nb[i-1]; @@ -330,7 +328,7 @@ bool ASMs2Dmx::connectPatch (int edge, ASM2D& neighbor, int nedge, bool revers, void ASMs2Dmx::closeBoundaries (int dir, int, int) { size_t nbi = 1; - for (size_t i = 0; i < m_basis.size(); nbi += nb[i++]) + for (size_t i = 0; i < nb.size(); nbi += nb[i++]) this->ASMs2D::closeBoundaries(dir,1+i,nbi); } @@ -413,7 +411,7 @@ Vec3 ASMs2Dmx::getCoord (size_t inod) const bool ASMs2Dmx::getSize (int& n1, int& n2, int basis) const { if (basis == 0) - return ASMs2D::getSize(n1,n2); + return this->ASMs2D::getSize(n1,n2); if (basis < 1 || basis > (int)m_basis.size()) return false; @@ -501,6 +499,7 @@ bool ASMs2Dmx::integrate (Integrand& integrand, PROFILE2("ASMs2Dmx::integrate(I)"); + bool use2ndDer = integrand.getIntegrandType() & Integrand::SECOND_DERIVATIVES; bool useElmVtx = integrand.getIntegrandType() & Integrand::ELEMENT_CORNERS; // Get Gaussian quadrature points and weights @@ -514,14 +513,22 @@ bool ASMs2Dmx::integrate (Integrand& integrand, this->getGaussPointParameters(gpar[d],d,nGauss,xg); // Evaluate basis function derivatives at all integration points - std::vector> splinex(m_basis.size()); - std::vector> splinex2(m_basis.size()); + std::vector> splinex; + std::vector> splinex2; + if (use2ndDer) + { + splinex2.resize(m_basis.size()); #pragma omp parallel for schedule(static) - for (size_t i=0;icomputeBasisGrid(gpar[0],gpar[1],splinex2[i]); - else + } + else + { + splinex.resize(m_basis.size()); +#pragma omp parallel for schedule(static) + for (size_t i = 0; i < m_basis.size(); i++) m_basis[i]->computeBasisGrid(gpar[0],gpar[1],splinex[i]); + } const int p1 = surf->order_u(); const int p2 = surf->order_v(); @@ -535,19 +542,20 @@ bool ASMs2Dmx::integrate (Integrand& integrand, // === Assembly loop over all elements in the patch ========================== - bool ok=true; - for (size_t g = 0; g < groups.size() && ok; g++) { + bool ok = true; + for (size_t g = 0; g < groups.size() && ok; g++) #pragma omp parallel for schedule(static) - for (size_t t = 0; t < groups[g].size(); t++) { + for (size_t t = 0; t < groups[g].size(); t++) + { MxFiniteElement fe(elem_size); - std::vector dNxdu(m_basis.size()); + std::vector dNxdu(m_basis.size()); std::vector d2Nxdu2(m_basis.size()); Matrix3D Hess; - double dXidu[2]; - Matrix Xnod, Jac; + double dXidu[2]; + Matrix Xnod, Jac; double param[3] = { 0.0, 0.0, 0.0 }; - Vec4 X(param); - for (size_t i = 0; i < groups[g][t].size() && ok; ++i) + Vec4 X(param); + for (size_t i = 0; i < groups[g][t].size() && ok; i++) { int iel = groups[g][t][i]; fe.iel = MLGE[iel]; @@ -557,7 +565,7 @@ bool ASMs2Dmx::integrate (Integrand& integrand, int i2 = p2 + iel / nel1; // Get element area in the parameter space - double dA = this->getParametricArea(++iel); + double dA = 0.25*this->getParametricArea(++iel); if (dA < 0.0) // topology error (probably logic error) { ok = false; @@ -609,10 +617,11 @@ bool ASMs2Dmx::integrate (Integrand& integrand, fe.v = param[1] = gpar[1](j+1,i2-p2+1); // Fetch basis function derivatives at current integration point - for (size_t b = 0; b < m_basis.size(); ++b) - if (integrand.getIntegrandType() & Integrand::SECOND_DERIVATIVES) + if (use2ndDer) + for (size_t b = 0; b < m_basis.size(); ++b) SplineUtils::extractBasis(splinex2[b][ip],fe.basis(b+1),dNxdu[b],d2Nxdu2[b]); - else + else + for (size_t b = 0; b < m_basis.size(); ++b) SplineUtils::extractBasis(splinex[b][ip],fe.basis(b+1),dNxdu[b]); // Compute Jacobian inverse of the coordinate mapping and @@ -620,12 +629,13 @@ bool ASMs2Dmx::integrate (Integrand& integrand, fe.detJxW = utl::Jacobian(Jac,fe.grad(geoBasis),Xnod, dNxdu[geoBasis-1]); if (fe.detJxW == 0.0) continue; // skip singular points + for (size_t b = 0; b < m_basis.size(); ++b) if (b != (size_t)geoBasis-1) fe.grad(b+1).multiply(dNxdu[b],Jac); // Compute Hessian of coordinate mapping and 2nd order derivatives - if (integrand.getIntegrandType() & Integrand::SECOND_DERIVATIVES) { + if (use2ndDer) { if (!utl::Hessian(Hess,fe.hess(geoBasis),Jac,Xnod, d2Nxdu2[geoBasis-1],fe.grad(geoBasis),true)) ok = false; @@ -645,7 +655,7 @@ bool ASMs2Dmx::integrate (Integrand& integrand, X.t = time.t; // Evaluate the integrand and accumulate element contributions - fe.detJxW *= 0.25*dA*wg[i]*wg[j]; + fe.detJxW *= dA*wg[i]*wg[j]; if (!integrand.evalIntMx(*A,fe,time,X)) ok = false; } @@ -661,7 +671,6 @@ bool ASMs2Dmx::integrate (Integrand& integrand, A->destruct(); } } - } return ok; } @@ -683,14 +692,14 @@ bool ASMs2Dmx::integrate (Integrand& integrand, int lIndex, if (!xg || !wg) return false; // Find the parametric direction of the edge normal {-2,-1, 1, 2} - const int edgeDir = (lIndex%10+1)/((lIndex%2) ? -2 : 2); + const int edgeDir = (lIndex%10+1) / ((lIndex%2) ? -2 : 2); const int t1 = abs(edgeDir); // Tangent direction normal to the patch edge const int t2 = 3-abs(edgeDir); // Tangent direction along the patch edge // Compute parameter values of the Gauss points along the whole patch edge std::array gpar; - for (short int d = 0; d < 2; d++) + for (int d = 0; d < 2; d++) if (-1-d == edgeDir) { gpar[d].resize(1,1); @@ -725,10 +734,11 @@ bool ASMs2Dmx::integrate (Integrand& integrand, int lIndex, fe.xi = fe.eta = edgeDir < 0 ? -1.0 : 1.0; fe.u = gpar[0](1,1); fe.v = gpar[1](1,1); + double param[3] = { fe.u, fe.v, 0.0 }; Matrices dNxdu(m_basis.size()); Matrix Xnod, Jac; - Vec4 X; + Vec4 X(param); Vec3 normal; @@ -757,7 +767,7 @@ bool ASMs2Dmx::integrate (Integrand& integrand, int lIndex, if (skipMe) continue; // Get element edge length in the parameter space - double dS = this->getParametricLength(iel,t2); + double dS = 0.5*this->getParametricLength(iel,t2); if (dS < 0.0) return false; // topology error (probably logic error) // Set up control point coordinates for current element @@ -778,16 +788,17 @@ bool ASMs2Dmx::integrate (Integrand& integrand, int lIndex, for (int i = 0; i < nGauss && ok; i++, ip++, fe.iGP++) { - // Parameter values of current integration point + // Local element coordinates and parameter values + // of current integration point if (gpar[0].size() > 1) { - fe.xi = xg[i]; - fe.u = gpar[0](i+1,i1-p1+1); + fe.xi = xg[i]; + fe.u = param[0] = gpar[0](i+1,i1-p1+1); } if (gpar[1].size() > 1) { - fe.eta = xg[i]; - fe.v = gpar[1](i+1,i2-p2+1); + fe.eta = xg[i]; + fe.v = param[1] = gpar[1](i+1,i2-p2+1); } // Fetch basis function derivatives at current integration point @@ -796,8 +807,10 @@ bool ASMs2Dmx::integrate (Integrand& integrand, int lIndex, // Compute Jacobian inverse of the coordinate mapping and // basis function derivatives w.r.t. Cartesian coordinates - fe.detJxW = utl::Jacobian(Jac,normal,fe.grad(geoBasis),Xnod,dNxdu[geoBasis-1],t1,t2); + fe.detJxW = utl::Jacobian(Jac,normal,fe.grad(geoBasis),Xnod, + dNxdu[geoBasis-1],t1,t2); if (fe.detJxW == 0.0) continue; // skip singular points + for (size_t b = 0; b < m_basis.size(); ++b) if (b != (size_t)geoBasis-1) fe.grad(b+1).multiply(dNxdu[b],Jac); @@ -805,12 +818,11 @@ bool ASMs2Dmx::integrate (Integrand& integrand, int lIndex, if (edgeDir < 0) normal *= -1.0; // Cartesian coordinates of current integration point - X .assign(Xnod * fe.basis(geoBasis)); + X.assign(Xnod * fe.basis(geoBasis)); X.t = time.t; // Evaluate the integrand and accumulate element contributions - fe.detJxW *= 0.5*dS*wg[i]; - + fe.detJxW *= dS*wg[i]; ok = integrand.evalBouMx(*A,fe,time,X,normal); } @@ -820,7 +832,7 @@ bool ASMs2Dmx::integrate (Integrand& integrand, int lIndex, // Assembly of global system integral if (ok && !glInt.assemble(A->ref(),fe.iel)) - return false; + ok = false; A->destruct(); @@ -926,7 +938,7 @@ bool ASMs2Dmx::integrate (Integrand& integrand, A_neigh->destruct(); // Get element edge length in the parameter space - double dS = this->getParametricLength(1+iel,t2); + double dS = 0.5*this->getParametricLength(1+iel,t2); if (dS < 0.0) // topology error (probably logic error) ok = false; @@ -989,7 +1001,7 @@ bool ASMs2Dmx::integrate (Integrand& integrand, } // Evaluate the integrand and accumulate element contributions - fe.detJxW *= 0.5*dS*wg[i]; + fe.detJxW *= dS*wg[i]; ok = integrand.evalIntMx(*A,fe,time,X,normal); } } @@ -1056,7 +1068,6 @@ bool ASMs2Dmx::evalSolution (Matrix& sField, const Vector& locSol, else std::copy(nfx.begin(), nfx.end(), nc.begin()); - if (std::inner_product(nb.begin(), nb.end(), nc.begin(), 0u) != locSol.size()) return false; @@ -1064,7 +1075,7 @@ bool ASMs2Dmx::evalSolution (Matrix& sField, const Vector& locSol, Vector Ytmp, Ztmp; // Evaluate the primary solution field at each point - size_t nPoints = splinex[0].size(); + size_t nPoints = splinex.front().size(); sField.resize(std::accumulate(nc.begin(), nc.end(), 0),nPoints); for (size_t i = 0; i < nPoints; i++) { @@ -1074,7 +1085,8 @@ bool ASMs2Dmx::evalSolution (Matrix& sField, const Vector& locSol, continue; IntVec ip; 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); + m_basis[b]->order_u(),m_basis[b]->order_v(), + splinex[b][i].left_idx,ip); utl::gather(ip,nc[b],locSol,Xtmp,comp); if (b == 0) @@ -1133,14 +1145,14 @@ bool ASMs2Dmx::evalSolution (Matrix& sField, const IntegrandBase& integrand, 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]); + 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; + for (int& i : ip[b]) i += ofs; ipa.insert(ipa.end(), ip[b].begin(), ip[b].end()); ofs += nb[b]; } @@ -1192,7 +1204,7 @@ void ASMs2Dmx::generateThreadGroups (const Integrand& integrand, bool silence, p2 = it->order_v(); } - ASMs2D::generateThreadGroups(p1-1, p2-1, silence, ignoreGlobalLM); + this->ASMs2D::generateThreadGroups(p1-1, p2-1, silence, ignoreGlobalLM); } @@ -1202,6 +1214,6 @@ void ASMs2Dmx::getBoundaryNodes (int lIndex, IntVec& nodes, int basis, if (basis > 0) this->ASMs2D::getBoundaryNodes(lIndex, nodes, basis, thick, 0, local); else - for (size_t b = 1; b <= this->getNoBasis(); ++b) + for (size_t b = 1; b <= m_basis.size(); b++) this->ASMs2D::getBoundaryNodes(lIndex, nodes, b, thick, 0, local); } diff --git a/src/ASM/ASMs2Dmx.h b/src/ASM/ASMs2Dmx.h index f3cfdb857..21d4ff277 100644 --- a/src/ASM/ASMs2Dmx.h +++ b/src/ASM/ASMs2Dmx.h @@ -35,7 +35,7 @@ class ASMs2Dmx : public ASMs2D, private ASMmxBase //! \brief The constructor initializes the dimension of each basis. ASMs2Dmx(unsigned char n_s, const CharVec& n_f); //! \brief Copy constructor. - ASMs2Dmx(const ASMs2Dmx& patch, const CharVec& n_f = CharVec(2,0)); + ASMs2Dmx(const ASMs2Dmx& patch, const CharVec& n_f); //! \brief Empty destructor. virtual ~ASMs2Dmx() {} @@ -50,9 +50,6 @@ class ASMs2Dmx : public ASMs2D, private ASMmxBase // Methods for model generation // ============================ - //! \brief Writes the geometry/basis of the patch to given stream. - virtual bool write(std::ostream& os, int basis = 0) const; - //! \brief Generates the finite element topology data for the patch. //! \details The data generated are the element-to-node connectivity array, //! the node-to-IJ-index array, as well as global node and element numbers. @@ -73,6 +70,9 @@ class ASMs2Dmx : public ASMs2D, private ASMmxBase //! \param[in] inod 1-based node index local to current patch virtual Vec3 getCoord(size_t inod) const; + //! \brief Writes the geometry/basis of the patch to given stream. + virtual bool write(std::ostream& os, int basis = 0) const; + //! \brief Returns the number of bases. virtual size_t getNoBasis() const { return m_basis.size(); } //! \brief Returns the total number of nodes in this patch. @@ -85,13 +85,6 @@ class ASMs2Dmx : public ASMs2D, private ASMmxBase //! \brief Returns the classification of a node. //! \param[in] inod 1-based node index local to current patch virtual char getNodeType(size_t inod) const; - //! \brief Returns the area in the parameter space for an element. - //! \param[in] iel Element index - virtual double getParametricArea(int iel) const; - //! \brief Returns boundary edge length in the parameter space for an element. - //! \param[in] iel Element index - //! \param[in] dir Local index of the boundary edge - double getParametricLength(int iel, int dir) const; //! \brief Initializes the patch level MADOF array for mixed problems. virtual void initMADOF(const int* sysMadof); @@ -105,7 +98,7 @@ class ASMs2Dmx : public ASMs2D, private ASMmxBase //! \param[in] coordCheck False to disable coordinate checks (periodic connections) //! \param[in] thick Thickness of connection virtual bool connectPatch(int edge, ASM2D& neighbor, int nedge, bool revers, - int basis = 0, bool coordCheck = true, int thick = 1); + int basis, bool coordCheck, int thick); //! \brief Makes two opposite boundary edges periodic. //! \param[in] dir Parameter direction defining the periodic edges @@ -137,7 +130,8 @@ class ASMs2Dmx : public ASMs2D, private ASMmxBase //! \param[in] time Parameters for nonlinear/time-dependent simulations //! \param[in] iChk Object checking if an element interface has contributions virtual bool integrate(Integrand& integrand, GlobalIntegral& glbInt, - const TimeDomain& time, const ASM::InterfaceChecker& iChk); + const TimeDomain& time, + const ASM::InterfaceChecker& iChk); // Post-processing methods @@ -199,18 +193,18 @@ class ASMs2Dmx : public ASMs2D, private ASMmxBase virtual void extractNodeVec(const Vector& globVec, Vector& nodeVec, unsigned char = 0, int basis = 0) const; - //! \brief Inject nodal results for this patch into a global vector. + //! \brief Injects nodal results for this patch into a global vector. //! \param[in] nodeVec Nodal result vector for this patch //! \param[out] globVec Global solution vector in DOF-order //! \param[in] basis Which basis (or 0 for both) to extract nodal values for virtual bool injectNodeVec(const Vector& nodeVec, Vector& globVec, - unsigned char = 0, int basis = 0) const; + unsigned char = 0, int basis = 0) const; using ASMs2D::generateThreadGroups; //! \brief Generates element groups for multi-threading of interior integrals. //! \param[in] integrand Object with problem-specific data and methods //! \param[in] silence If \e true, suppress threading group outprint - //! \param[in] ignoreGlobalLM If \e true, ignore global multipliers in sanity check + //! \param[in] ignoreGlobalLM Sanity check option virtual void generateThreadGroups(const Integrand& integrand, bool silence, bool ignoreGlobalLM); @@ -221,6 +215,15 @@ class ASMs2Dmx : public ASMs2D, private ASMmxBase //! \param[in] basis Which basis to return size parameters for virtual bool getSize(int& n1, int& n2, int basis = 0) const; +protected: + //! \brief Returns the area in the parameter space for an element. + //! \param[in] iel Element index + double getParametricArea(int iel) const; + //! \brief Returns boundary edge length in the parameter space for an element. + //! \param[in] iel Element index + //! \param[in] dir Local index of the boundary edge + double getParametricLength(int iel, int dir) const; + //! \brief Finds the global (or patch-local) node numbers on a patch boundary. //! \param[in] lIndex Local index of the boundary edge //! \param nodes Array of node numbers @@ -228,9 +231,9 @@ class ASMs2Dmx : public ASMs2D, private ASMmxBase //! \param[in] thick Thickness of connection //! \param[in] local If \e true return patch-local node numbers virtual void getBoundaryNodes(int lIndex, IntVec& nodes, int basis, - int thick = 1, int = 0, bool local = false) const; + int thick, int, bool local) const; -protected: +private: std::vector> m_basis; //!< Vector of bases }; diff --git a/src/ASM/ASMs3Dmx.C b/src/ASM/ASMs3Dmx.C index 5f5b2587d..c7f9ed30e 100644 --- a/src/ASM/ASMs3Dmx.C +++ b/src/ASM/ASMs3Dmx.C @@ -27,7 +27,6 @@ #include "Utilities.h" #include "Profiler.h" #include "Vec3Oper.h" -#include "Vec3.h" #include #include #ifdef USE_OPENMP @@ -36,14 +35,16 @@ ASMs3Dmx::ASMs3Dmx (const CharVec& n_f) - : ASMs3D(std::accumulate(n_f.begin(), n_f.end(), 0)), ASMmxBase(n_f) + : ASMs3D(std::accumulate(n_f.begin(),n_f.end(),0)), ASMmxBase(n_f) { } ASMs3Dmx::ASMs3Dmx (const ASMs3Dmx& patch, const CharVec& n_f) - : ASMs3D(patch), ASMmxBase(n_f), m_basis(patch.m_basis) + : ASMs3D(patch), ASMmxBase(n_f), + m_basis(patch.m_basis), projBasis(patch.projBasis) { + nb = patch.nb; } @@ -58,11 +59,13 @@ Go::SplineVolume* ASMs3Dmx::getBasis (int basis) const Go::SplineSurface* ASMs3Dmx::getBoundary (int dir, int basis) { - if (dir < -3 || dir == 0 || dir > 3) + if (basis < 1 || basis > (int)m_basis.size()) return nullptr; - // The boundary surfaces are stored internally in the SplineVolume object int iface = dir > 0 ? 2*dir-1 : -2*dir-2; + if (iface < 0 || iface > 5) + return nullptr; + return m_basis[basis-1]->getBoundarySurface(iface).get(); } @@ -70,11 +73,11 @@ Go::SplineSurface* ASMs3Dmx::getBoundary (int dir, int basis) bool ASMs3Dmx::write (std::ostream& os, int basis) const { if (basis == -1) - os <<"700 1 0 0\n" << *projBasis; + os <<"700 1 0 0\n"<< *projBasis; else if (basis < 1 || basis > (int)m_basis.size()) - os <<"700 1 0 0\n" << *svol; + os <<"700 1 0 0\n"<< *svol; else if (m_basis[basis-1]) - os <<"700 1 0 0\n" << *m_basis[basis-1]; + os <<"700 1 0 0\n"<< *m_basis[basis-1]; else return false; @@ -104,10 +107,7 @@ size_t ASMs3Dmx::getNoNodes (int basis) const unsigned char ASMs3Dmx::getNoFields (int basis) const { - if (basis < 1 || basis > (int)nfx.size()) - return std::accumulate(nfx.begin(), nfx.end(), 0); - - return nfx[basis-1]; + return basis > 0 && basis <= (int)nfx.size() ? nfx[basis-1] : nf; } @@ -116,12 +116,12 @@ unsigned char ASMs3Dmx::getNodalDOFs (size_t inod) const if (this->isLMn(inod)) return nLag; - size_t nbc=0; - for (size_t i=0;iisLMn(inod)) return 'L'; - size_t nbc=nb.front(); + size_t nbc = nb.front(); if (inod <= nbc) return 'D'; - else for (size_t i = 1; i < nb.size(); i++) + + for (size_t i = 1; i < nb.size(); i++) if (inod <= (nbc += nb[i])) return 'O'+i; @@ -184,31 +185,32 @@ bool ASMs3Dmx::generateFEMTopology () else projBasis = m_basis.front(); } - delete svol; geomB = svol = m_basis[geoBasis-1]->clone(); nb.clear(); - elem_size.clear(); nb.reserve(m_basis.size()); + elem_size.clear(); elem_size.reserve(m_basis.size()); for (auto& it : m_basis) { nb.push_back(it->numCoefs(0)*it->numCoefs(1)*it->numCoefs(2)); elem_size.push_back(it->order(0)*it->order(1)*it->order(2)); } + nnod = std::accumulate(nb.begin(), nb.end(), 0u); if (!nodeInd.empty() && !shareFE) { - if (nodeInd.size() == std::accumulate(nb.begin(), nb.end(), 0u)) return true; + if (nodeInd.size() == nnod) + return true; + std::cerr <<" *** ASMs3Dmx::generateFEMTopology: Inconsistency between the" <<" number of FE nodes "<< nodeInd.size() - <<"\n and the number of spline coefficients " - << std::accumulate(nb.begin(), nb.end(), 0) + <<"\n and the number of spline coefficients "<< nnod <<" in the patch."<< std::endl; return false; } - - if (shareFE == 'F') return true; + else if (shareFE == 'F') + return true; #ifdef SP_DEBUG for (auto it : m_basis) { @@ -230,18 +232,16 @@ bool ASMs3Dmx::generateFEMTopology () (m_basis[geoBasis-1]->numCoefs(1)-m_basis[geoBasis-1]->order(1)+1)* (m_basis[geoBasis-1]->numCoefs(2)-m_basis[geoBasis-1]->order(2)+1); - nnod = std::accumulate(nb.begin(), nb.end(), 0); - myMLGE.resize(nel,0); myMLGN.resize(nnod); myMNPC.resize(nel); myNodeInd.resize(nnod); int i1, i2, i3, j1, j2, j3; - size_t iel, inod = 0; - for (auto& it : m_basis) { + size_t iel = 0, inod = 0; + for (auto& it : m_basis) for (i3 = 0; i3 < it->numCoefs(2); i3++) - for (i2 = 0; i2 < it->numCoefs(1); i2++) + for (i2 = 0; i2 < it->numCoefs(1); i2++) for (i1 = 0; i1 < it->numCoefs(0); i1++) { myNodeInd[inod].I = i1; @@ -249,19 +249,16 @@ bool ASMs3Dmx::generateFEMTopology () myNodeInd[inod].K = i3; myMLGN[inod++] = ++gNod; } - } int lnod2 = 0; int lnod3 = 0; - - iel = 0, inod = std::accumulate(nb.begin(),nb.begin()+geoBasis-1,0u); - for (i2 = 0; i2 < geoBasis-1; ++i2) lnod2 += m_basis[i2]->order(0)*m_basis[i2]->order(1)*m_basis[i2]->order(2); for (i2 = 0; i2 < (int)m_basis.size(); ++i2) lnod3 += m_basis[i2]->order(0)*m_basis[i2]->order(1)*m_basis[i2]->order(2); // Create nodal connectivities for bases + inod = std::accumulate(nb.begin(),nb.begin()+geoBasis-1,0u); Go::SplineVolume* b = m_basis[geoBasis-1].get(); auto knotw = b->basis(2).begin(); for (i3 = 1; i3 <= b->numCoefs(2); i3++, ++knotw) { @@ -286,7 +283,7 @@ bool ASMs3Dmx::generateFEMTopology () // find knotspan spanning element for other bases lnod = 0; size_t lnod4 = 0; - for (size_t bas = 0; bas < m_basis.size(); ++bas) { + for (size_t bas = 0; bas < m_basis.size(); ++bas) { if (bas != (size_t)geoBasis-1) { double ku = *knotu; double kv = *knotv; @@ -330,9 +327,10 @@ bool ASMs3Dmx::connectPatch (int face, ASM3D& neighbor, int nface, int norient, nface = 11-nface; size_t nb1 = 0, nb2 = 0; - for (size_t i = 1; i <= m_basis.size(); ++i) { + for (size_t i = 1; i <= nb.size(); i++) { if (basis == 0 || i == (size_t)basis) - if (!this->connectBasis(face,*neighMx,nface,norient,i,nb1,nb2,coordCheck,thick)) + if (!this->connectBasis(face,*neighMx,nface,norient,i,nb1,nb2, + coordCheck,thick)) return false; nb1 += nb[i-1]; @@ -347,7 +345,7 @@ bool ASMs3Dmx::connectPatch (int face, ASM3D& neighbor, int nface, int norient, void ASMs3Dmx::closeBoundaries (int dir, int, int) { size_t nbi = 1; - for (size_t i = 0; i < m_basis.size(); nbi += nb[i++]) + for (size_t i = 0; i < nb.size(); nbi += nb[i++]) this->ASMs3D::closeBoundaries(dir,1+i,nbi); } @@ -440,12 +438,84 @@ bool ASMs3Dmx::getSize (int& n1, int& n2, int& n3, int basis) const } +#define DERR -999.99 + +double ASMs3Dmx::getParametricVolume (int iel) const +{ +#ifdef INDEX_CHECK + if (iel < 1 || (size_t)iel > MNPC.size()) + { + std::cerr <<" *** ASMs3Dmx::getParametricVolume: Element index "<< iel + <<" out of range [1,"<< MNPC.size() <<"]."<< std::endl; + return DERR; + } +#endif + if (MNPC[iel-1].empty()) + return 0.0; + + int inod1 = MNPC[iel-1][std::accumulate(elem_size.begin(), + elem_size.begin()+geoBasis, -1)]; +#ifdef INDEX_CHECK + if (inod1 < 0 || (size_t)inod1 >= nnod) + { + std::cerr <<" *** ASMs3Dmx::getParametricVolume: Node index "<< inod1 + <<" out of range [0,"<< nnod <<">."<< std::endl; + return DERR; + } +#endif + + double du = svol->knotSpan(0,nodeInd[inod1].I); + double dv = svol->knotSpan(1,nodeInd[inod1].J); + double dw = svol->knotSpan(2,nodeInd[inod1].K); + return du*dv*dw; +} + + +double ASMs3Dmx::getParametricArea (int iel, int dir) const +{ +#ifdef INDEX_CHECK + if (iel < 1 || (size_t)iel > MNPC.size()) + { + std::cerr <<" *** ASMs3Dmx::getParametricArea: Element index "<< iel + <<" out of range [1,"<< MNPC.size() <<"]."<< std::endl; + return DERR; + } +#endif + if (MNPC[iel-1].empty()) + return 0.0; + + int inod1 = MNPC[iel-1][std::accumulate(elem_size.begin(), + elem_size.begin()+geoBasis, -1)]; +#ifdef INDEX_CHECK + if (inod1 < 0 || (size_t)inod1 >= nnod) + { + std::cerr <<" *** ASMs3Dmx::getParametricArea: Node index "<< inod1 + <<" out of range [0,"<< nnod <<">."<< std::endl; + return DERR; + } +#endif + + const int ni = nodeInd[inod1].I; + const int nj = nodeInd[inod1].J; + const int nk = nodeInd[inod1].K; + switch (dir) + { + case 1: return svol->knotSpan(1,nj)*svol->knotSpan(2,nk); + case 2: return svol->knotSpan(0,ni)*svol->knotSpan(2,nk); + case 3: return svol->knotSpan(0,ni)*svol->knotSpan(1,nj); + } + + std::cerr <<" *** ASMs3Dmx::getParametricArea: Invalid face direction " + << dir << std::endl; + return DERR; +} + + bool ASMs3Dmx::integrate (Integrand& integrand, GlobalIntegral& glInt, const TimeDomain& time) { if (!svol) return true; // silently ignore empty patches - if (m_basis.empty()) return false; PROFILE2("ASMs3Dmx::integrate(I)"); @@ -463,16 +533,22 @@ bool ASMs3Dmx::integrate (Integrand& integrand, this->getGaussPointParameters(gpar[d],d,nGauss,xg); // Evaluate basis function derivatives at all integration points - std::vector> splinex(m_basis.size()); - std::vector> splinex2(m_basis.size()); + std::vector> splinex; + std::vector> splinex2; if (use2ndDer) + { + splinex2.resize(m_basis.size()); #pragma omp parallel for schedule(static) - for (size_t i=0;icomputeBasisGrid(gpar[0],gpar[1],gpar[2],splinex2[i]); + } else + { + splinex.resize(m_basis.size()); #pragma omp parallel for schedule(static) - for (size_t i=0;icomputeBasisGrid(gpar[0],gpar[1],gpar[2],splinex[i]); + } const int p1 = svol->order(0); const int p2 = svol->order(1); @@ -491,18 +567,19 @@ bool ASMs3Dmx::integrate (Integrand& integrand, // === Assembly loop over all elements in the patch ========================== bool ok = true; - for (size_t g = 0; g < groups.size() && ok; ++g) { + for (size_t g = 0; g < groups.size() && ok; g++) #pragma omp parallel for schedule(static) - for (size_t t = 0; t < groups[g].size(); ++t) { + for (size_t t = 0; t < groups[g].size(); t++) + { MxFiniteElement fe(elem_size); - std::vector dNxdu(m_basis.size()); + std::vector dNxdu(m_basis.size()); std::vector d2Nxdu2(m_basis.size()); Matrix3D Hess; - double dXidu[3]; - Matrix Xnod, Jac; + double dXidu[3]; + Matrix Xnod, Jac; double param[3]; - Vec4 X(param); - for (size_t l = 0; l < groups[g][t].size() && ok; ++l) + Vec4 X(param); + for (size_t l = 0; l < groups[g][t].size() && ok; l++) { int iel = groups[g][t][l]; fe.iel = MLGE[iel]; @@ -513,9 +590,9 @@ bool ASMs3Dmx::integrate (Integrand& integrand, int i3 = p3 + iel / (nel1*nel2); // Get element volume in the parameter space - double dV = this->getParametricVolume(++iel); + double dV = 0.125*this->getParametricVolume(++iel); if (dV < 0.0) // topology error (probably logic error) - { + { ok = false; break; } @@ -581,6 +658,7 @@ bool ASMs3Dmx::integrate (Integrand& integrand, fe.detJxW = utl::Jacobian(Jac,fe.grad(geoBasis),Xnod, dNxdu[geoBasis-1]); if (fe.detJxW == 0.0) continue; // skip singular points + for (size_t b = 0; b < m_basis.size(); ++b) if (b != (size_t)geoBasis-1) fe.grad(b+1).multiply(dNxdu[b],Jac); @@ -606,7 +684,7 @@ bool ASMs3Dmx::integrate (Integrand& integrand, X.t = time.t; // Evaluate the integrand and accumulate element contributions - fe.detJxW *= 0.125*dV*wg[i]*wg[j]*wg[k]; + fe.detJxW *= dV*wg[i]*wg[j]*wg[k]; if (!integrand.evalIntMx(*A,fe,time,X)) ok = false; } @@ -622,7 +700,6 @@ bool ASMs3Dmx::integrate (Integrand& integrand, A->destruct(); } } - } return ok; } @@ -633,7 +710,6 @@ bool ASMs3Dmx::integrate (Integrand& integrand, int lIndex, const TimeDomain& time) { if (!svol) return true; // silently ignore empty patches - if (m_basis.empty()) return false; PROFILE2("ASMs3Dmx::integrate(B)"); @@ -654,7 +730,7 @@ bool ASMs3Dmx::integrate (Integrand& integrand, int lIndex, if (!xg || !wg) return false; // Find the parametric direction of the face normal {-3,-2,-1, 1, 2, 3} - const int faceDir = (lIndex%10+1)/((lIndex%2) ? -2 : 2); + const int faceDir = (lIndex%10+1) / ((lIndex%2) ? -2 : 2); const int t1 = 1 + abs(faceDir)%3; // first tangent direction const int t2 = 1 + t1%3; // second tangent direction @@ -709,11 +785,11 @@ bool ASMs3Dmx::integrate (Integrand& integrand, int lIndex, fe.u = gpar[0](1,1); fe.v = gpar[1](1,1); fe.w = gpar[2](1,1); + double param[3] = { fe.u, fe.v, fe.w }; Matrices dNxdu(m_basis.size()); Matrix Xnod, Jac; - double param[3] = { fe.u, fe.v, fe.w }; - Vec4 X; + Vec4 X(param); Vec3 normal; for (size_t l = 0; l < threadGrp[g][t].size() && ok; ++l) { @@ -726,7 +802,7 @@ bool ASMs3Dmx::integrate (Integrand& integrand, int lIndex, int i3 = p3 + iel / (nel1*nel2); // Get element face area in the parameter space - double dA = this->getParametricArea(++iel,abs(faceDir)); + double dA = 0.25*this->getParametricArea(++iel,abs(faceDir)); if (dA < 0.0) // topology error (probably logic error) { ok = false; @@ -804,8 +880,10 @@ bool ASMs3Dmx::integrate (Integrand& integrand, int lIndex, // Compute Jacobian inverse of the coordinate mapping and // basis function derivatives w.r.t. Cartesian coordinates - fe.detJxW = utl::Jacobian(Jac,normal,fe.grad(geoBasis),Xnod,dNxdu[geoBasis-1],t1,t2); + fe.detJxW = utl::Jacobian(Jac,normal,fe.grad(geoBasis),Xnod, + dNxdu[geoBasis-1],t1,t2); if (fe.detJxW == 0.0) continue; // skip singular points + for (size_t b = 0; b < m_basis.size(); ++b) if (b != (size_t)geoBasis-1) fe.grad(b+1).multiply(dNxdu[b],Jac); @@ -813,11 +891,11 @@ bool ASMs3Dmx::integrate (Integrand& integrand, int lIndex, if (faceDir < 0) normal *= -1.0; // Cartesian coordinates of current integration point - X = Xnod * fe.basis(geoBasis); + X.assign(Xnod * fe.basis(geoBasis)); X.t = time.t; // Evaluate the integrand and accumulate element contributions - fe.detJxW *= 0.25*dA*wg[i]*wg[j]; + fe.detJxW *= dA*wg[i]*wg[j]; if (!integrand.evalBouMx(*A,fe,time,X,normal)) ok = false; } @@ -1081,8 +1159,6 @@ bool ASMs3Dmx::evalSolution (Matrix& sField, const Vector& locSol, const RealArray* gpar, bool regular, int, int nf) const { - if (m_basis.empty()) return false; - // Evaluate the basis functions at all points std::vector> splinex(m_basis.size()); if (regular) @@ -1114,7 +1190,7 @@ bool ASMs3Dmx::evalSolution (Matrix& sField, const Vector& locSol, Vector Ytmp, Ztmp; // Evaluate the primary solution field at each point - size_t nPoints = splinex[0].size(); + size_t nPoints = splinex.front().size(); sField.resize(std::accumulate(nc.begin(), nc.end(), 0),nPoints); for (size_t i = 0; i < nPoints; i++) { @@ -1148,8 +1224,6 @@ bool ASMs3Dmx::evalSolution (Matrix& sField, const IntegrandBase& integrand, { sField.resize(0,0); - if (m_basis.empty()) return false; - // Evaluate the basis functions and their derivatives at all points std::vector> splinex(m_basis.size()); if (regular) @@ -1193,8 +1267,7 @@ bool ASMs3Dmx::evalSolution (Matrix& sField, const IntegrandBase& integrand, if (b == (size_t)geoBasis-1) utl::gather(ip[geoBasis-1], 3, Xnod, Xtmp); - for (auto& it : ip[b]) - it += ofs; + for (int& i : ip[b]) i += ofs; ipa.insert(ipa.end(), ip[b].begin(), ip[b].end()); ofs += nb[b]; } @@ -1245,7 +1318,8 @@ void ASMs3Dmx::generateThreadGroups (const Integrand& integrand, bool silence, if (it->order(d) > p[d]) p[d] = it->order(d); - ASMs3D::generateThreadGroups(p[0]-1,p[1]-1,p[2]-1,silence,ignoreGlobalLM); + this->ASMs3D::generateThreadGroups(p[0]-1, p[1]-1, p[2]-1, + silence, ignoreGlobalLM); } @@ -1254,80 +1328,7 @@ void ASMs3Dmx::generateThreadGroups (char lIndex, bool silence, bool) #ifdef USE_OPENMP omp_set_num_threads(1); #endif - ASMs3D::generateThreadGroups(lIndex,silence,false); -} - - -#define DERR -999.99 - -double ASMs3Dmx::getParametricVolume (int iel) const -{ -#ifdef INDEX_CHECK - if (iel < 1 || (size_t)iel > MNPC.size()) - { - std::cerr <<" *** ASMs3Dmx::getParametricVolume: Element index "<< iel - <<" out of range [1,"<< MNPC.size() <<"]."<< std::endl; - return DERR; - } -#endif - if (MNPC[iel-1].empty()) - return 0.0; - - int inod1 = MNPC[iel-1][std::accumulate(elem_size.begin(), - elem_size.begin()+geoBasis, -1)]; -#ifdef INDEX_CHECK - if (inod1 < 0 || (size_t)inod1 >= nnod) - { - std::cerr <<" *** ASMs3Dmx::getParametricVolume: Node index "<< inod1 - <<" out of range [0,"<< nnod <<">."<< std::endl; - return DERR; - } -#endif - - double du = svol->knotSpan(0,nodeInd[inod1].I); - double dv = svol->knotSpan(1,nodeInd[inod1].J); - double dw = svol->knotSpan(2,nodeInd[inod1].K); - return du*dv*dw; -} - - -double ASMs3Dmx::getParametricArea (int iel, int dir) const -{ -#ifdef INDEX_CHECK - if (iel < 1 || (size_t)iel > MNPC.size()) - { - std::cerr <<" *** ASMs3Dmx::getParametricArea: Element index "<< iel - <<" out of range [1,"<< MNPC.size() <<"]."<< std::endl; - return DERR; - } -#endif - if (MNPC[iel-1].empty()) - return 0.0; - - int inod1 = MNPC[iel-1][std::accumulate(elem_size.begin(), - elem_size.begin()+geoBasis, -1)]; -#ifdef INDEX_CHECK - if (inod1 < 0 || (size_t)inod1 >= nnod) - { - std::cerr <<" *** ASMs3Dmx::getParametricArea: Node index "<< inod1 - <<" out of range [0,"<< nnod <<">."<< std::endl; - return DERR; - } -#endif - - const int ni = nodeInd[inod1].I; - const int nj = nodeInd[inod1].J; - const int nk = nodeInd[inod1].K; - switch (dir) - { - case 1: return svol->knotSpan(1,nj)*svol->knotSpan(2,nk); - case 2: return svol->knotSpan(0,ni)*svol->knotSpan(2,nk); - case 3: return svol->knotSpan(0,ni)*svol->knotSpan(1,nj); - } - - std::cerr <<" *** ASMs3Dmx::getParametricArea: Invalid face direction " - << dir << std::endl; - return DERR; + this->ASMs3D::generateThreadGroups(lIndex, silence, false); } @@ -1337,21 +1338,21 @@ void ASMs3Dmx::getBoundaryNodes (int lIndex, IntVec& nodes, int basis, if (basis > 0) this->ASMs3D::getBoundaryNodes(lIndex, nodes, basis, thick, 0, local); else - for (size_t b = 1; b <= this->getNoBasis(); ++b) + for (size_t b = 1; b <= m_basis.size(); b++) this->ASMs3D::getBoundaryNodes(lIndex, nodes, b, thick, 0, local); } -Fields* ASMs3Dmx::getProjectedFields(const Vector& coefs, size_t nf) const +Fields* ASMs3Dmx::getProjectedFields (const Vector& coefs, size_t nf) const { - if (projBasis != m_basis[0]) + if (projBasis != m_basis.front()) return new SplineFields3D(projBasis.get(), coefs, nf); return nullptr; } -size_t ASMs3Dmx::getNoProjectionNodes() const +size_t ASMs3Dmx::getNoProjectionNodes () const { return projBasis->numCoefs(0) * projBasis->numCoefs(1) * @@ -1365,7 +1366,7 @@ bool ASMs3Dmx::evalProjSolution (Matrix& sField, const Vector& locSol, #ifdef SP_DEBUG std::cout <<"ASMu3Dmx::evalProjSolution(Matrix&,const Vector&,const int*,int)\n"; #endif - if (projBasis == m_basis[0]) + if (projBasis == m_basis.front()) return this->evalSolution(sField, locSol, npe, nf); // Compute parameter values of the nodal points diff --git a/src/ASM/ASMs3Dmx.h b/src/ASM/ASMs3Dmx.h index 1734ab2dc..4e7b8823b 100644 --- a/src/ASM/ASMs3Dmx.h +++ b/src/ASM/ASMs3Dmx.h @@ -35,7 +35,7 @@ class ASMs3Dmx : public ASMs3D, private ASMmxBase //! \brief The constructor initializes the dimension of each basis. explicit ASMs3Dmx(const CharVec& n_f); //! \brief Copy constructor. - ASMs3Dmx(const ASMs3Dmx& patch, const CharVec& n_f = CharVec(2,0)); + ASMs3Dmx(const ASMs3Dmx& patch, const CharVec& n_f); //! \brief Empty destructor. virtual ~ASMs3Dmx() {} @@ -98,7 +98,7 @@ class ASMs3Dmx : public ASMs3D, private ASMmxBase //! \param[in] coordCheck False to disable coordinate checks (periodic connections) //! \param[in] thick Thickness of connection virtual bool connectPatch(int face, ASM3D& neighbor, int nface, int norient, - int basis = 0, bool coordCheck = true, int thick = 1); + int basis, bool coordCheck, int thick); //! \brief Makes two opposite boundary faces periodic. //! \param[in] dir Parameter direction defining the periodic faces @@ -130,7 +130,9 @@ class ASMs3Dmx : public ASMs3D, private ASMmxBase //! \param[in] time Parameters for nonlinear/time-dependent simulations //! \param[in] iChk Object checking if an element interface has contributions virtual bool integrate(Integrand& integrand, GlobalIntegral& glbInt, - const TimeDomain& time, const ASM::InterfaceChecker& iChk); + const TimeDomain& time, + const ASM::InterfaceChecker& iChk); + // Post-processing methods // ======================= @@ -157,7 +159,7 @@ class ASMs3Dmx : public ASMs3D, private ASMmxBase //! \param[in] gpar Parameter values of the result sampling points //! \param[in] regular Flag indicating how the sampling points are defined //! \param[in] deriv Derivative order to return - //! \param[in] nf If non-zero evaluates nf fields on first basis + //! \param[in] nf If nonzero, evaluate nf fields on first basis //! //! \details When \a regular is \e true, it is assumed that the parameter //! value array \a gpar forms a regular tensor-product point grid of dimension @@ -217,7 +219,7 @@ class ASMs3Dmx : public ASMs3D, private ASMmxBase //! \brief Generates element groups for multi-threading of interior integrals. //! \param[in] integrand Object with problem-specific data and methods //! \param[in] silence If \e true, suppress threading group outprint - //! \param[in] ignoreGlobalLM If \e true ignore global multipliers in sanity check + //! \param[in] ignoreGlobalLM Sanity check option virtual void generateThreadGroups(const Integrand& integrand, bool silence, bool ignoreGlobalLM); //! \brief Generates element groups for multi-threading of boundary integrals. @@ -231,13 +233,14 @@ class ASMs3Dmx : public ASMs3D, private ASMmxBase //! \param[out] n3 Number of nodes in third (w) direction //! \param[in] basis Which basis to return size parameters for virtual bool getSize(int& n1, int& n2, int& n3, int basis = 0) const; + protected: //! \brief Returns the volume in the parameter space for an element. //! \param[in] iel Element index double getParametricVolume(int iel) const; //! \brief Returns boundary face area in the parameter space for an element. //! \param[in] iel Element index - //! \param[in] dir Local face index of the boundary face + //! \param[in] dir Local index of the boundary face double getParametricArea(int iel, int dir) const; //! \brief Finds the global (or patch-local) node numbers on a patch boundary. @@ -246,8 +249,8 @@ class ASMs3Dmx : public ASMs3D, private ASMmxBase //! \param[in] basis Which basis to grab nodes for (0 for all) //! \param[in] thick Thickness of connection //! \param[in] local If \e true return patch-local node numbers - virtual void getBoundaryNodes(int lIndex, IntVec& nodes, int basis = 0, - int thick = 1, int = 0, bool local = false) const; + virtual void getBoundaryNodes(int lIndex, IntVec& nodes, int basis, + int thick, int, bool local) const; //! \brief Assembles L2-projection matrices for the secondary solution. //! \param[out] A Left-hand-side matrix @@ -258,6 +261,7 @@ class ASMs3Dmx : public ASMs3D, private ASMmxBase const IntegrandBase& integrand, bool continuous) const; +private: std::vector> m_basis; //!< Vector of bases std::shared_ptr projBasis; //!< Basis to project onto }; From d9cc3ebfa1ad95f3ac24433835b3171c3ac41def Mon Sep 17 00:00:00 2001 From: Knut Morten Okstad Date: Thu, 20 Sep 2018 17:43:26 +0200 Subject: [PATCH 2/2] Added: Overload the constrain methods for mixed to also allow specifying DOFs to constrain for all bases (when invoked with basis=0) --- src/ASM/ASMmxBase.C | 27 +++++++++-- src/ASM/ASMmxBase.h | 12 ++++- src/ASM/ASMs2Dmx.C | 26 +++++++++++ src/ASM/ASMs2Dmx.h | 16 ++++++- src/ASM/ASMs3Dmx.C | 41 +++++++++++++++++ src/ASM/ASMs3Dmx.h | 29 +++++++++++- src/ASM/LR/ASMu2Dmx.C | 62 +++++++++++++++++-------- src/ASM/LR/ASMu2Dmx.h | 20 +++++++- src/ASM/LR/ASMu3Dmx.C | 103 +++++++++++++++++++++++------------------- src/ASM/LR/ASMu3Dmx.h | 40 +++++++++++++--- src/SIM/SIMinput.C | 2 +- 11 files changed, 293 insertions(+), 85 deletions(-) diff --git a/src/ASM/ASMmxBase.C b/src/ASM/ASMmxBase.C index 4cf46ecc2..e3dafbcc9 100644 --- a/src/ASM/ASMmxBase.C +++ b/src/ASM/ASMmxBase.C @@ -12,6 +12,7 @@ //============================================================================== #include "ASMmxBase.h" +#include "Utilities.h" #include "GoTools/geometry/SplineSurface.h" #include "GoTools/geometry/SurfaceInterpolator.h" #include "GoTools/trivariate/SplineVolume.h" @@ -123,8 +124,8 @@ bool ASMmxBase::getSolutionMx (Matrix& sField, const Vector& locSol, } -ASMmxBase::SurfaceVec ASMmxBase::establishBases(Go::SplineSurface* surf, - MixedType type) +ASMmxBase::SurfaceVec ASMmxBase::establishBases (Go::SplineSurface* surf, + MixedType type) { SurfaceVec result(2); // With mixed methods we need two separate spline spaces @@ -248,8 +249,8 @@ ASMmxBase::SurfaceVec ASMmxBase::establishBases(Go::SplineSurface* surf, } -ASMmxBase::VolumeVec ASMmxBase::establishBases(Go::SplineVolume* svol, - MixedType type) +ASMmxBase::VolumeVec ASMmxBase::establishBases (Go::SplineVolume* svol, + MixedType type) { VolumeVec result(2); // With mixed methods we need two separate spline spaces @@ -498,3 +499,21 @@ Go::SplineVolume* ASMmxBase::raiseBasis (Go::SplineVolume* svol) // Project the coordinates onto the new basis (the 2nd XYZ is dummy here) return Go::VolumeInterpolator::regularInterpolation(basis[0],basis[1],basis[2],ug[0],ug[1],ug[2],XYZ,ndim,false,XYZ); } + + +int ASMmxBase::maskDOFs (int dofs, char basis) const +{ + unsigned char ofs = std::accumulate(nfx.begin(),nfx.begin()+basis-1,0u); + std::set allDofs = utl::getDigits(dofs); + dofs = 0; + + // Convert the DOF digits to local values of this basis, + // and mask off those not residing on this basis + for (int dof : allDofs) + { + dofs *= 10; + if (dof > ofs && dof <= ofs+nfx[basis-1]) + dofs += dof - ofs; + } + return dofs; +} diff --git a/src/ASM/ASMmxBase.h b/src/ASM/ASMmxBase.h index 253163024..7b8de00d0 100644 --- a/src/ASM/ASMmxBase.h +++ b/src/ASM/ASMmxBase.h @@ -74,8 +74,10 @@ class ASMmxBase static char geoBasis; //!< 1-based index of basis representing the geometry protected: - typedef std::vector> SurfaceVec; //!< Convenience type - typedef std::vector> VolumeVec; //!< Convenience type + typedef std::shared_ptr SurfacePtr; //!< Pointer to spline + typedef std::shared_ptr VolumePtr; //!< Pointer to spline + typedef std::vector SurfaceVec; //!< Convenience type + typedef std::vector VolumeVec; //!< Convenience type //! \brief Establish mixed bases in 2D. //! \param[in] surf The base basis to use @@ -94,6 +96,12 @@ class ASMmxBase //! \brief Returns a C^p-1 basis of one degree higher than \a *svol. static Go::SplineVolume* raiseBasis(Go::SplineVolume* svol); + //! \brief Mask off DOFs not residing on the specified basis + //! \param[in] dofs Encoded DOF indices (like 123456 meaning dofs 1 to 6) + //! \param[in] basis Which basis to return DOF indices for + //! \return Encoded DOF indices related to the specified basis only + int maskDOFs(int dofs, char basis) const; + private: std::vector MADOF; //!< Matrix of accumulated DOFs for this patch diff --git a/src/ASM/ASMs2Dmx.C b/src/ASM/ASMs2Dmx.C index 721862d7f..0b2cee0fa 100644 --- a/src/ASM/ASMs2Dmx.C +++ b/src/ASM/ASMs2Dmx.C @@ -303,6 +303,32 @@ bool ASMs2Dmx::generateFEMTopology () } +void ASMs2Dmx::constrainEdge (int dir, bool open, int dof, int code, char basis) +{ + if (basis > 0) + this->ASMs2D::constrainEdge(dir,open,dof,code,basis); + else for (basis = 1; basis <= (char)nfx.size(); basis++) + { + int basisDofs = this->maskDOFs(dof,basis); + if (basisDofs > 0) + this->ASMs2D::constrainEdge(dir,open,basisDofs,code,basis); + } +} + + +void ASMs2Dmx::constrainCorner (int I, int J, int dof, int code, char basis) +{ + if (basis > 0) + this->ASMs2D::constrainCorner(I,J,dof,code,basis); + else for (basis = 1; basis <= (char)nfx.size(); basis++) + { + int basisDofs = this->maskDOFs(dof,basis); + if (basisDofs > 0) + this->ASMs2D::constrainCorner(I,J,basisDofs,code,basis); + } +} + + bool ASMs2Dmx::connectPatch (int edge, ASM2D& neighbor, int nedge, bool revers, int basis, bool coordCheck, int thick) { diff --git a/src/ASM/ASMs2Dmx.h b/src/ASM/ASMs2Dmx.h index 21d4ff277..cea05d1ef 100644 --- a/src/ASM/ASMs2Dmx.h +++ b/src/ASM/ASMs2Dmx.h @@ -16,7 +16,6 @@ #include "ASMs2D.h" #include "ASMmxBase.h" -#include /*! @@ -89,6 +88,21 @@ class ASMs2Dmx : public ASMs2D, private ASMmxBase //! \brief Initializes the patch level MADOF array for mixed problems. virtual void initMADOF(const int* sysMadof); + //! \brief Constrains all DOFs on a given boundary edge. + //! \param[in] dir Parameter direction defining the edge to constrain + //! \param[in] open If \e true, exclude the end points of the edge + //! \param[in] dof Which DOFs to constrain at each node along the edge + //! \param[in] code Inhomogeneous dirichlet condition code + //! \param[in] basis Which basis to constrain edge for (0 means check all) + virtual void constrainEdge(int dir, bool open, int dof, int code, char basis); + //! \brief Constrains a corner node identified by the two parameter indices. + //! \param[in] I Parameter index in u-direction + //! \param[in] J Parameter index in v-direction + //! \param[in] dof Which DOFs to constrain at the node + //! \param[in] code Inhomogeneous dirichlet condition code + //! \param[in] basis Which basis to constrain node for (0 means check all) + virtual void constrainCorner(int I, int J, int dof, int code, char basis); + //! \brief Connects all matching nodes on two adjacent boundary edges. //! \param[in] edge Local edge index of this patch, in range [1,4] //! \param neighbor The neighbor patch diff --git a/src/ASM/ASMs3Dmx.C b/src/ASM/ASMs3Dmx.C index c7f9ed30e..9fd4e7fea 100644 --- a/src/ASM/ASMs3Dmx.C +++ b/src/ASM/ASMs3Dmx.C @@ -314,6 +314,47 @@ bool ASMs3Dmx::generateFEMTopology () } +void ASMs3Dmx::constrainFace (int dir, bool open, int dof, int code, char basis) +{ + if (basis > 0) + this->ASMs3D::constrainFace(dir,open,dof,code,basis); + else for (basis = 1; basis <= (char)nfx.size(); basis++) + { + int basisDofs = this->maskDOFs(dof,basis); + if (basisDofs > 0) + this->ASMs3D::constrainFace(dir,open,basisDofs,code,basis); + } +} + + +void ASMs3Dmx::constrainEdge (int lEdge, bool open, int dof, + int code, char basis) +{ + if (basis > 0) + this->ASMs3D::constrainEdge(lEdge,open,dof,code,basis); + else for (basis = 1; basis <= (char)nfx.size(); basis++) + { + int basisDofs = this->maskDOFs(dof,basis); + if (basisDofs > 0) + this->ASMs3D::constrainEdge(lEdge,open,basisDofs,code,basis); + } +} + + +void ASMs3Dmx::constrainCorner (int I, int J, int K, int dof, + int code, char basis) +{ + if (basis > 0) + this->ASMs3D::constrainCorner(I,J,K,dof,code,basis); + else for (basis = 1; basis <= (char)nfx.size(); basis++) + { + int basisDofs = this->maskDOFs(dof,basis); + if (basisDofs > 0) + this->ASMs3D::constrainCorner(I,J,K,basisDofs,code,basis); + } +} + + bool ASMs3Dmx::connectPatch (int face, ASM3D& neighbor, int nface, int norient, int basis, bool coordCheck, int thick) { diff --git a/src/ASM/ASMs3Dmx.h b/src/ASM/ASMs3Dmx.h index 4e7b8823b..c2493b3d3 100644 --- a/src/ASM/ASMs3Dmx.h +++ b/src/ASM/ASMs3Dmx.h @@ -16,7 +16,6 @@ #include "ASMs3D.h" #include "ASMmxBase.h" -#include /*! @@ -52,7 +51,7 @@ class ASMs3Dmx : public ASMs3D, private ASMmxBase //! \brief Generates the finite element topology data for the patch. //! \details The data generated are the element-to-node connectivity array, - //! the node-to-IJ-index array, as well as global node and element numbers. + //! the node-to-IJK-index array, as well as global node and element numbers. virtual bool generateFEMTopology(); //! \brief Clears the contents of the patch, making it empty. @@ -89,6 +88,32 @@ class ASMs3Dmx : public ASMs3D, private ASMmxBase //! \brief Initializes the patch level MADOF array for mixed problems. virtual void initMADOF(const int* sysMadof); + //! \brief Constrains all DOFs on a given boundary face. + //! \param[in] dir Parameter direction defining the face to constrain + //! \param[in] open If \e true, exclude all points along the face boundary + //! \param[in] dof Which DOFs to constrain at each node on the face + //! \param[in] code Inhomogeneous dirichlet condition code + //! \param[in] basis Which basis to constrain face for (0 means check all) + virtual void constrainFace(int dir, bool open, int dof, + int code, char basis); + //! \brief Constrains all DOFs on a given boundary edge. + //! \param[in] lEdge Local index [1,12] of the edge to constrain + //! \param[in] open If \e true, exclude the end points of the edge + //! \param[in] dof Which DOFs to constrain at each node along the edge + //! \param[in] code Inhomogeneous dirichlet condition code + //! \param[in] basis Which basis to constrain edge for (0 means check all) + virtual void constrainEdge(int lEdge, bool open, int dof, + int code, char basis); + //! \brief Constrains a corner node identified by the three parameter indices. + //! \param[in] I Parameter index in u-direction + //! \param[in] J Parameter index in v-direction + //! \param[in] K Parameter index in w-direction + //! \param[in] dof Which DOFs to constrain at the node + //! \param[in] code Inhomogeneous dirichlet condition code + //! \param[in] basis Which basis to constrain node for (0 means check all) + virtual void constrainCorner(int I, int J, int K, int dof, + int code, char basis); + //! \brief Connects all matching nodes on two adjacent boundary faces. //! \param[in] face Local face index of this patch, in range [1,6] //! \param neighbor The neighbor patch diff --git a/src/ASM/LR/ASMu2Dmx.C b/src/ASM/LR/ASMu2Dmx.C index c88000971..a4af0413f 100644 --- a/src/ASM/LR/ASMu2Dmx.C +++ b/src/ASM/LR/ASMu2Dmx.C @@ -286,6 +286,32 @@ bool ASMu2Dmx::generateFEMTopology () } +void ASMu2Dmx::constrainEdge (int dir, bool open, int dof, int code, char basis) +{ + if (basis > 0) + this->ASMu2D::constrainEdge(dir,open,dof,code,basis); + else for (basis = 1; basis <= (char)nfx.size(); basis++) + { + int basisDofs = this->maskDOFs(dof,basis); + if (basisDofs > 0) + this->ASMu2D::constrainEdge(dir,open,basisDofs,code,basis); + } +} + + +void ASMu2Dmx::constrainCorner (int I, int J, int dof, int code, char basis) +{ + if (basis > 0) + this->ASMu2D::constrainCorner(I,J,dof,code,basis); + else for (basis = 1; basis <= (char)nfx.size(); basis++) + { + int basisDofs = this->maskDOFs(dof,basis); + if (basisDofs > 0) + this->ASMu2D::constrainCorner(I,J,basisDofs,code,basis); + } +} + + bool ASMu2Dmx::integrate (Integrand& integrand, GlobalIntegral& glInt, const TimeDomain& time) @@ -315,6 +341,7 @@ bool ASMu2Dmx::integrate (Integrand& integrand, { if (!ok) continue; + int iel = groups[t][e] + 1; auto el1 = threadBasis->elementBegin()+iel-1; double uh = ((*el1)->umin()+(*el1)->umax())/2.0; @@ -339,7 +366,7 @@ bool ASMu2Dmx::integrate (Integrand& integrand, double dXidu[2]; // Get element area in the parameter space - double dA = this->getParametricArea(geoEl); + double dA = 0.25*this->getParametricArea(geoEl); if (dA < 0.0) // topology error (probably logic error) { ok = false; @@ -377,7 +404,8 @@ bool ASMu2Dmx::integrate (Integrand& integrand, continue; } - // --- Integration loop over all Gauss points in each direction ------------ + + // --- Integration loop over all Gauss points in each direction ---------- int jp = (iel-1)*nGauss*nGauss; fe.iGP = firstIp + jp; // Global integration point counter @@ -409,8 +437,10 @@ bool ASMu2Dmx::integrate (Integrand& integrand, // Compute Jacobian inverse of coordinate mapping and derivatives // basis function derivatives w.r.t. Cartesian coordinates - fe.detJxW = utl::Jacobian(Jac,fe.grad(geoBasis),Xnod,dNxdu[geoBasis-1]); + fe.detJxW = utl::Jacobian(Jac,fe.grad(geoBasis),Xnod, + dNxdu[geoBasis-1]); if (fe.detJxW == 0.0) continue; // skip singular points + for (size_t b = 0; b < m_basis.size(); ++b) if (b != (size_t)geoBasis-1) fe.grad(b+1).multiply(dNxdu[b],Jac); @@ -436,27 +466,18 @@ bool ASMu2Dmx::integrate (Integrand& integrand, X.t = time.t; // Evaluate the integrand and accumulate element contributions - fe.detJxW *= 0.25*dA*wg[i]*wg[j]; + fe.detJxW *= dA*wg[i]*wg[j]; if (!integrand.evalIntMx(*A,fe,time,X)) - { ok = false; - continue; - } } // Finalize the element quantities - if (!integrand.finalizeElement(*A,time,firstIp+jp)) - { + if (ok && !integrand.finalizeElement(*A,time,firstIp+jp)) ok = false; - continue; - } // Assembly of global system integral - if (!glInt.assemble(A->ref(),fe.iel)) - { + if (ok && !glInt.assemble(A->ref(),fe.iel)) ok = false; - continue; - } A->destruct(); } @@ -787,12 +808,13 @@ bool ASMu2Dmx::integrate (Integrand& integrand, // Compute Jacobian inverse of the coordinate mapping and // basis function derivatives w.r.t. Cartesian coordinates - fe.detJxW = utl::Jacobian(Jac2, normal, - fe.grad(geoBasis+m_basis.size()), - Xnod2,dNxdu[geoBasis-1+m_basis.size()],t1,t2); - fe.detJxW = utl::Jacobian(Jac, normal, - fe.grad(geoBasis),Xnod,dNxdu[geoBasis-1],t1,t2); + fe.detJxW = utl::Jacobian(Jac2,normal, + fe.grad(geoBasis+m_basis.size()),Xnod2, + dNxdu[geoBasis-1+m_basis.size()],t1,t2); + fe.detJxW = utl::Jacobian(Jac,normal,fe.grad(geoBasis),Xnod, + dNxdu[geoBasis-1],t1,t2); if (fe.detJxW == 0.0) continue; // skip singular points + for (size_t b = 0; b < m_basis.size(); ++b) if (b != (size_t)geoBasis-1) { fe.grad(b+1).multiply(dNxdu[b],Jac); diff --git a/src/ASM/LR/ASMu2Dmx.h b/src/ASM/LR/ASMu2Dmx.h index a72165f05..4786c5855 100644 --- a/src/ASM/LR/ASMu2Dmx.h +++ b/src/ASM/LR/ASMu2Dmx.h @@ -79,6 +79,21 @@ class ASMu2Dmx : public ASMu2D, private ASMmxBase //! \brief Initializes the patch level MADOF array for mixed problems. virtual void initMADOF(const int* sysMadof); + //! \brief Constrains all DOFs on a given boundary edge. + //! \param[in] dir Parameter direction defining the edge to constrain + //! \param[in] open If \e true, exclude the end points of the edge + //! \param[in] dof Which DOFs to constrain at each node along the edge + //! \param[in] code Inhomogeneous dirichlet condition code + //! \param[in] basis Which basis to constrain edge for (0 means check all) + virtual void constrainEdge(int dir, bool open, int dof, int code, char basis); + //! \brief Constrains a corner node identified by the two parameter indices. + //! \param[in] I Parameter index in u-direction + //! \param[in] J Parameter index in v-direction + //! \param[in] dof Which DOFs to constrain at the node + //! \param[in] code Inhomogeneous dirichlet condition code + //! \param[in] basis Which basis to constrain node for (0 means check all) + virtual void constrainCorner(int I, int J, int dof, int code, char basis); + // Methods for integration of finite element quantities. // These are the main computational methods of the ASM class hierarchy. // ==================================================================== @@ -104,7 +119,8 @@ class ASMu2Dmx : public ASMu2D, private ASMmxBase //! \param[in] time Parameters for nonlinear/time-dependent simulations //! \param[in] iChk Object checking if an element interface has contributions virtual bool integrate(Integrand& integrand, GlobalIntegral& glbInt, - const TimeDomain& time, const ASM::InterfaceChecker& iChk); + const TimeDomain& time, + const ASM::InterfaceChecker& iChk); // Post-processing methods @@ -156,7 +172,7 @@ class ASMu2Dmx : public ASMu2D, private ASMmxBase virtual void extractNodeVec(const Vector& globVec, Vector& nodeVec, unsigned char = 0, int basis = 0) const; - //! \brief Inject nodal results for this patch into a global vector. + //! \brief Injects nodal results for this patch into a global vector. //! \param[in] nodeVec Nodal result vector for this patch //! \param[out] globVec Global solution vector in DOF-order //! \param[in] basis Which basis (or 0 for both) to extract nodal values for diff --git a/src/ASM/LR/ASMu3Dmx.C b/src/ASM/LR/ASMu3Dmx.C index 52cda4243..14729e1ed 100644 --- a/src/ASM/LR/ASMu3Dmx.C +++ b/src/ASM/LR/ASMu3Dmx.C @@ -290,6 +290,47 @@ bool ASMu3Dmx::generateFEMTopology () } +void ASMu3Dmx::constrainFace (int dir, bool open, int dof, int code, char basis) +{ + if (basis > 0) + this->ASMu3D::constrainFace(dir,open,dof,code,basis); + else for (basis = 1; basis <= (char)nfx.size(); basis++) + { + int basisDofs = this->maskDOFs(dof,basis); + if (basisDofs > 0) + this->ASMu3D::constrainFace(dir,open,basisDofs,code,basis); + } +} + + +void ASMu3Dmx::constrainEdge (int lEdge, bool open, int dof, + int code, char basis) +{ + if (basis > 0) + this->ASMu3D::constrainEdge(lEdge,open,dof,code,basis); + else for (basis = 1; basis <= (char)nfx.size(); basis++) + { + int basisDofs = this->maskDOFs(dof,basis); + if (basisDofs > 0) + this->ASMu3D::constrainEdge(lEdge,open,basisDofs,code,basis); + } +} + + +void ASMu3Dmx::constrainCorner (int I, int J, int K, int dof, + int code, char basis) +{ + if (basis > 0) + this->ASMu3D::constrainCorner(I,J,K,dof,code,basis); + else for (basis = 1; basis <= (char)nfx.size(); basis++) + { + int basisDofs = this->maskDOFs(dof,basis); + if (basisDofs > 0) + this->ASMu3D::constrainCorner(I,J,K,basisDofs,code,basis); + } +} + + bool ASMu3Dmx::integrate (Integrand& integrand, GlobalIntegral& glInt, const TimeDomain& time) @@ -311,54 +352,32 @@ bool ASMu3Dmx::integrate (Integrand& integrand, int p1 = lrspline->order(0); int p2 = lrspline->order(1); int p3 = lrspline->order(2); + double u[2*p1], v[2*p2], w[2*p3]; Go::BsplineBasis basis1 = getBezierBasis(p1); Go::BsplineBasis basis2 = getBezierBasis(p2); Go::BsplineBasis basis3 = getBezierBasis(p3); - - BN[b-1].resize(p1*p2*p3, nGauss*nGauss*nGauss); + BN [b-1].resize(p1*p2*p3, nGauss*nGauss*nGauss); BdNdu[b-1].resize(p1*p2*p3, nGauss*nGauss*nGauss); BdNdv[b-1].resize(p1*p2*p3, nGauss*nGauss*nGauss); BdNdw[b-1].resize(p1*p2*p3, nGauss*nGauss*nGauss); - int ig=1; // gauss point iterator - for(int zeta=0; zeta 0) - { - xr = GaussQuadrature::getCoord(nRed); - wr = GaussQuadrature::getWeight(nRed); - if (!xr || !wr) return false; } - else if (nRed < 0) - nRed = nGauss; // The integrand needs to know nGauss ThreadGroups oneGroup; if (glInt.threadSafe()) oneGroup.oneGroup(nel); @@ -413,14 +432,9 @@ bool ASMu3Dmx::integrate (Integrand& integrand, } // Compute parameter values of the Gauss points over the whole element - std::array gpar, redpar; + std::array gpar; for (int d = 0; d < 3; d++) - { this->getGaussPointParameters(gpar[d],d,nGauss,iEl+1,xg); - if (xr) - this->getGaussPointParameters(redpar[d],d,nRed,iEl+1,xr); - } - if (integrand.getIntegrandType() & Integrand::ELEMENT_CORNERS) fe.h = this->getElementCorners(iEl+1, fe.XC); @@ -481,11 +495,6 @@ bool ASMu3Dmx::integrate (Integrand& integrand, break; } - if (xr) - { - std::cerr << "Haven't really figured out what this part does yet\n"; - exit(42142); - } // --- Integration loop over all Gauss points in each direction -------- @@ -554,8 +563,7 @@ bool ASMu3Dmx::integrate (Integrand& integrand, fe.detJxW *= 0.125*vol*wg[i]*wg[j]*wg[k]; if (!integrand.evalIntMx(*A,fe,time,X)) ok = false; - - } // end gauss integrand + } // Finalize the element quantities if (ok && !integrand.finalizeElement(*A,time,0)) @@ -655,7 +663,7 @@ bool ASMu3Dmx::integrate (Integrand& integrand, int lIndex, double dXidu[3]; // Get element face area in the parameter space - double dA = this->getParametricArea(iEl+1,abs(faceDir)); + double dA = 0.25*this->getParametricArea(iEl+1,abs(faceDir)); if (dA < 0.0) // topology error (probably logic error) { ok = false; @@ -744,7 +752,7 @@ bool ASMu3Dmx::integrate (Integrand& integrand, int lIndex, X.t = time.t; // Evaluate the integrand and accumulate element contributions - fe.detJxW *= 0.25*dA*wg[i]*wg[j]; + fe.detJxW *= dA*wg[i]*wg[j]; if (!integrand.evalBouMx(*A,fe,time,X,normal)) ok = false; } @@ -756,6 +764,7 @@ bool ASMu3Dmx::integrate (Integrand& integrand, int lIndex, // Assembly of global system integral if (ok && !glInt.assemble(A->ref(),fe.iel)) ok = false; + A->destruct(); firstp += nGP*nGP*nGP; diff --git a/src/ASM/LR/ASMu3Dmx.h b/src/ASM/LR/ASMu3Dmx.h index 73f1cb9b1..94ce7b501 100644 --- a/src/ASM/LR/ASMu3Dmx.h +++ b/src/ASM/LR/ASMu3Dmx.h @@ -79,6 +79,33 @@ class ASMu3Dmx : public ASMu3D, private ASMmxBase //! \brief Initializes the patch level MADOF array for mixed problems. virtual void initMADOF(const int* sysMadof); + //! \brief Constrains all DOFs on a given boundary face. + //! \param[in] dir Parameter direction defining the face to constrain + //! \param[in] open If \e true, exclude all points along the face boundary + //! \param[in] dof Which DOFs to constrain at each node on the face + //! \param[in] code Inhomogeneous dirichlet condition code + //! \param[in] basis Which basis to constrain face for (0 means check all) + virtual void constrainFace(int dir, bool open, int dof, + int code, char basis); + //! \brief Constrains all DOFs on a given boundary edge. + //! \param[in] lEdge Local index [1,12] of the edge to constrain + //! \param[in] open If \e true, exclude the end points of the edge + //! \param[in] dof Which DOFs to constrain at each node along the edge + //! \param[in] code Inhomogeneous dirichlet condition code + //! \param[in] basis Which basis to constrain edge for (0 means check all) + virtual void constrainEdge(int lEdge, bool open, int dof, + int code, char basis); + //! \brief Constrains a corner node identified by the three parameter indices. + //! \param[in] I Parameter index in u-direction + //! \param[in] J Parameter index in v-direction + //! \param[in] K Parameter index in w-direction + //! \param[in] dof Which DOFs to constrain at the node + //! \param[in] code Inhomogeneous dirichlet condition code + //! \param[in] basis Which basis to constrain node for (0 means check all) + virtual void constrainCorner(int I, int J, int K, int dof, + int code, char basis); + + // Methods for integration of finite element quantities. // These are the main computational methods of the ASM class hierarchy. // ==================================================================== @@ -156,7 +183,7 @@ class ASMu3Dmx : public ASMu3D, private ASMmxBase virtual void extractNodeVec(const Vector& globVec, Vector& nodeVec, unsigned char = 0, int basis = 0) const; - //! \brief Inject nodal results for this patch into a global vector. + //! \brief Injects nodal results for this patch into a global vector. //! \param[in] nodeVec Nodal result vector for this patch //! \param[out] globVec Global solution vector in DOF-order //! \param[in] basis Which basis (or 0 for both) to extract nodal values for @@ -165,10 +192,8 @@ class ASMu3Dmx : public ASMu3D, private ASMmxBase //! \brief Returns the number of projection nodes for this patch. virtual size_t getNoProjectionNodes() const; - //! \brief Returns the number of refinement nodes for this patch. virtual size_t getNoRefineNodes() const; - //! \brief Returns the number of refinement elements for this patch. virtual size_t getNoRefineElms() const; @@ -201,9 +226,12 @@ class ASMu3Dmx : public ASMu3D, private ASMmxBase bool continuous) const; private: - std::vector> m_basis; //!< Spline bases - std::shared_ptr projBasis; //!< Basis to project onto - std::shared_ptr refBasis; //!< Basis to refine based on + typedef std::shared_ptr SplinePtr; //!< Pointer to spline + + std::vector m_basis; //!< All bases + SplinePtr projBasis; //!< Basis to project onto + SplinePtr refBasis; //!< Basis to refine based on + const std::vector& bezierExtractmx; //!< Bezier extraction matrices std::vector myBezierExtractmx; //!< Bezier extraction matrices }; diff --git a/src/SIM/SIMinput.C b/src/SIM/SIMinput.C index 74b3d2c96..76c24f368 100644 --- a/src/SIM/SIMinput.C +++ b/src/SIM/SIMinput.C @@ -449,7 +449,7 @@ bool SIMinput::parseBCTag (const TiXmlElement* elem) else if (!strcasecmp(elem->Value(),"dirichlet") && !ignoreDirichlet) { const TiXmlNode* dval = nullptr; - int comp = 0, symm = 0, basis = 1; + int comp = 0, symm = 0, basis = 0; std::string set, type, axes; utl::getAttribute(elem,"set",set); utl::getAttribute(elem,"type",type,true);