diff --git a/src/ASM/ASMs2D.C b/src/ASM/ASMs2D.C index e9cff1aa..6ecd0034 100644 --- a/src/ASM/ASMs2D.C +++ b/src/ASM/ASMs2D.C @@ -1663,6 +1663,9 @@ bool ASMs2D::integrate (Integrand& integrand, fe.detJxW = utl::Jacobian(Jac,fe.dNdX,Xnod,dNdu); if (fe.detJxW == 0.0) continue; // skip singular points + // Store tangent vectors in fe.G for shells + if (nsd > 2) fe.G = Jac; + // Cartesian coordinates of current integration point X.assign(Xnod * fe.N); X.t = time.t; @@ -1710,6 +1713,8 @@ bool ASMs2D::integrate (Integrand& integrand, // Compute G-matrix if (integrand.getIntegrandType() & Integrand::G_MATRIX) utl::getGmat(Jac,dXidu,fe.G); + else if (nsd > 2) + fe.G = Jac; // Store tangent vectors in fe.G for shells #if SP_DEBUG > 4 if (iel == dbgElm || iel == -dbgElm || dbgElm == 0) @@ -1718,6 +1723,8 @@ bool ASMs2D::integrate (Integrand& integrand, <<"\nN ="<< fe.N <<"dNdX ="<< fe.dNdX; if (!fe.d2NdX2.empty()) std::cout <<"d2NdX2 ="<< fe.d2NdX2; + if (!fe.G.empty()) + std::cout <<"G ="<< fe.G; } #endif @@ -1745,7 +1752,7 @@ bool ASMs2D::integrate (Integrand& integrand, A->destruct(); #ifdef SP_DEBUG - if (iel == -dbgElm) break; // Skipping all elements, except for -dbgElm + if (iel == -dbgElm) break; // Skipping all elements, except for -dbgElm #endif } } @@ -1900,6 +1907,8 @@ bool ASMs2D::integrate (Integrand& integrand, // Compute G-matrix if (integrand.getIntegrandType() & Integrand::G_MATRIX) utl::getGmat(Jac,dXidu,fe.G); + else if (nsd > 2) + fe.G = Jac; // Store tangent vectors in fe.G for shells #if SP_DEBUG > 4 if (iel == dbgElm || iel == -dbgElm || dbgElm == 0) @@ -1932,7 +1941,7 @@ bool ASMs2D::integrate (Integrand& integrand, A->destruct(); #ifdef SP_DEBUG - if (iel == -dbgElm) break; // Skipping all elements, except for -dbgElm + if (iel == -dbgElm) break; // Skipping all elements, except for -dbgElm #endif } } @@ -2056,6 +2065,9 @@ bool ASMs2D::integrate (Integrand& integrand, if (edgeDir < 0) normal *= -1.0; + // Store tangent vectors in fe.G for shells + if (nsd > 2) fe.G = Jac; + // Cartesian coordinates of current integration point X.assign(Xnod * fe.N); X.t = time.t; @@ -2253,6 +2265,8 @@ bool ASMs2D::integrate (Integrand& integrand, int lIndex, // Compute G-matrix if (integrand.getIntegrandType() & Integrand::G_MATRIX) utl::getGmat(Jac,dXidu,fe.G); + else if (nsd > 2) + fe.G = Jac; // Store tangent vectors in fe.G for shells // Cartesian coordinates of current integration point X.assign(Xnod * fe.N); @@ -2643,6 +2657,9 @@ bool ASMs2D::evalSolution (Matrix& sField, const IntegrandBase& integrand, if (!utl::Hessian(Hess,fe.d2NdX2,Jac,Xtmp,d2Ndu2,fe.dNdX)) continue; + // Store tangent vectors in fe.G for shells + if (nsd > 2) fe.G = Jac; + // Now evaluate the solution field if (!integrand.evalSol(solPt,fe,Xtmp*fe.N,ip)) return false; diff --git a/src/ASM/LR/ASMu2D.C b/src/ASM/LR/ASMu2D.C index 7115323e..454e2298 100644 --- a/src/ASM/LR/ASMu2D.C +++ b/src/ASM/LR/ASMu2D.C @@ -337,7 +337,7 @@ bool ASMu2D::checkElementSize (int elmId, bool globalNum) const elmId = it - MLGE.begin(); } - if (elmId < lrspline->nElements()) + if (elmId >= 0 && elmId < lrspline->nElements()) return lrspline->getElement(elmId)->area() > aMin+1.0e-12; else return false; @@ -406,7 +406,7 @@ bool ASMu2D::raiseOrder (int ru, int rv) bool ASMu2D::evaluateBasis (FiniteElement& fe, int derivs) const { - LR::Element* el = lrspline->getElement(fe.iel-1); + const LR::Element* el = lrspline->getElement(fe.iel-1); if (!el) return false; fe.xi = 2.0*(fe.u - el->umin()) / (el->umax() - el->umin()) - 1.0; @@ -870,7 +870,7 @@ double ASMu2D::getParametricLength (int iel, int dir) const } #endif - LR::Element* el = lrspline->getElement(iel-1); + const LR::Element* el = lrspline->getElement(iel-1); switch (dir) { case 1: return el->vmax() - el->vmin(); @@ -894,7 +894,7 @@ bool ASMu2D::getElementCoordinates (Matrix& X, int iel) const } #endif - LR::Element* el = lrspline->getElement(iel-1); + const LR::Element* el = lrspline->getElement(iel-1); X.resize(nsd,el->nBasisFunctions()); int n = 1; @@ -979,7 +979,7 @@ double ASMu2D::getElementCorners (int iel, Vec3Vec& XC) const } #endif - LR::Element* el = lrspline->getElement(iel-1); + const LR::Element* el = lrspline->getElement(iel-1); double u[4] = { el->umin(), el->umax(), el->umin(), el->umax() }; double v[4] = { el->vmin(), el->vmin(), el->vmax(), el->vmax() }; @@ -1036,12 +1036,12 @@ bool ASMu2D::integrate (Integrand& integrand, continue; int iel = threadGroups[0][t][e] + 1; - #ifdef SP_DEBUG if (dbgElm < 0 && iel != -dbgElm) continue; // Skipping all elements, except for -dbgElm #endif - FiniteElement fe(MNPC[iel-1].size()); + + FiniteElement fe; fe.iel = MLGE[iel-1]; Matrix dNdu, Xnod, Jac; Matrix3D d2Ndu2, Hess; @@ -1083,21 +1083,22 @@ bool ASMu2D::integrate (Integrand& integrand, double u0 = 0.5*(gpar[0].front() + gpar[0].back()); double v0 = 0.5*(gpar[1].front() + gpar[1].back()); lrspline->point(X0,u0,v0); - for (unsigned char i = 0; i < nsd; i++) - X[i] = X0[i]; + X = SplineUtils::toVec3(X0,nsd); } if (integrand.getIntegrandType() & Integrand::G_MATRIX) { // Element size in parametric space - dXidu[0] = lrspline->getElement(iel-1)->umax()-lrspline->getElement(iel-1)->umin(); - dXidu[1] = lrspline->getElement(iel-1)->vmax()-lrspline->getElement(iel-1)->vmin(); + const LR::Element* el = lrspline->getElement(iel-1); + dXidu[0] = el->umax() - el->umin(); + dXidu[1] = el->vmax() - el->vmin(); } // Initialize element quantities - LocalIntegral* A = integrand.getLocalIntegral(fe.N.size(),fe.iel); + LocalIntegral* A = integrand.getLocalIntegral(MNPC[iel-1].size(),fe.iel); if (!integrand.initElement(MNPC[iel-1],fe,X,nRed*nRed,*A)) { + A->destruct(); ok = false; continue; } @@ -1124,6 +1125,10 @@ bool ASMu2D::integrate (Integrand& integrand, // Compute Jacobian inverse and derivatives fe.detJxW = utl::Jacobian(Jac,fe.dNdX,Xnod,dNdu); + if (fe.detJxW == 0.0) continue; // skip singular points + + // Store tangent vectors in fe.G for shells + if (nsd > 2) fe.G = Jac; // Cartesian coordinates of current integration point X.assign(Xnod * fe.N); @@ -1132,10 +1137,7 @@ bool ASMu2D::integrate (Integrand& integrand, // Compute the reduced integration terms of the integrand fe.detJxW *= 0.25*dA*wr[i]*wr[j]; if (!integrand.reducedInt(*A,fe,X)) - { ok = false; - continue; - } } } @@ -1197,10 +1199,12 @@ bool ASMu2D::integrate (Integrand& integrand, // Compute G-matrix if (integrand.getIntegrandType() & Integrand::G_MATRIX) utl::getGmat(Jac,dXidu,fe.G); + else if (nsd > 2) + fe.G = Jac; // Store tangent vectors in fe.G for shells #if SP_DEBUG > 4 if (iel == dbgElm || iel == -dbgElm || dbgElm == 0) - std::cout <<"\nN ="<< fe.N <<"dNdX ="<< fe.dNdX; + std::cout <<"\n"<< fe; #endif // Cartesian coordinates of current integration point @@ -1209,33 +1213,26 @@ bool ASMu2D::integrate (Integrand& integrand, // Evaluate the integrand and accumulate element contributions fe.detJxW *= 0.25*dA*wg[i]*wg[j]; +#ifndef USE_OPENMP PROFILE3("Integrand::evalInt"); +#endif if (!integrand.evalInt(*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(); #ifdef SP_DEBUG if (iel == -dbgElm) - continue; // Skipping all elements, except for -dbgElm + break; // Skipping all elements, except for -dbgElm #endif } } @@ -1281,7 +1278,8 @@ bool ASMu2D::integrate (Integrand& integrand, if (dbgElm < 0 && iel != -dbgElm) continue; // Skipping all elements, except for -dbgElm #endif - FiniteElement fe(MNPC[iel-1].size()); + + FiniteElement fe; fe.iel = MLGE[iel-1]; Matrix dNdu, Xnod, Jac; Matrix3D d2Ndu2, Hess; @@ -1313,7 +1311,7 @@ bool ASMu2D::integrate (Integrand& integrand, } // Initialize element quantities - LocalIntegral* A = integrand.getLocalIntegral(fe.N.size(),fe.iel); + LocalIntegral* A = integrand.getLocalIntegral(MNPC[iel-1].size(),fe.iel); if (!integrand.initElement(MNPC[iel-1],fe,X,0,*A)) { ok = false; @@ -1332,7 +1330,7 @@ bool ASMu2D::integrate (Integrand& integrand, fe.u = elmPts[ip][0]; fe.v = elmPts[ip][1]; - // Compute basis function derivatives at current integration point + // Compute basis function derivatives at current integration point if (integrand.getIntegrandType() & Integrand::SECOND_DERIVATIVES) { Go::BasisDerivsSf2 spline; lrspline->computeBasis(fe.u,fe.v,spline,iel-1); @@ -1356,10 +1354,12 @@ bool ASMu2D::integrate (Integrand& integrand, continue; } + // Store tangent vectors in fe.G for shells + if (nsd > 2) fe.G = Jac; + #if SP_DEBUG > 4 if (iel == dbgElm || iel == -dbgElm || dbgElm == 0) - std::cout <<"\niel, ip = "<< iel <<" "<< ip - <<"\nN ="<< fe.N <<"dNdX ="<< fe.dNdX; + std::cout <<"\n"<< fe; #endif // Cartesian coordinates of current integration point @@ -1369,29 +1369,27 @@ bool ASMu2D::integrate (Integrand& integrand, // Evaluate the integrand and accumulate element contributions fe.detJxW *= 0.25*dA*elmPts[ip][2]; +#ifndef USE_OPENMP PROFILE3("Integrand::evalInt"); +#endif if (!integrand.evalInt(*A,fe,time,X)) - { ok = false; - continue; - } } // Finalize the element quantities - if (!integrand.finalizeElement(*A,time,firstIp+MPitg[iel])) - { + if (ok && !integrand.finalizeElement(*A,time,firstIp+MPitg[iel])) 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(); + +#ifdef SP_DEBUG + if (iel == -dbgElm) + break; // Skipping all elements, except for -dbgElm +#endif } } @@ -1438,6 +1436,7 @@ bool ASMu2D::integrate (Integrand& integrand, int lIndex, std::map::const_iterator iit = firstBp.find(lIndex%10); size_t firstp = iit == firstBp.end() ? 0 : iit->second; + FiniteElement fe; Matrix dNdu, Xnod, Jac; double param[3] = { 0.0, 0.0, 0.0 }; Vec4 X(param); @@ -1473,10 +1472,10 @@ bool ASMu2D::integrate (Integrand& integrand, int lIndex, if (!this->getElementCoordinates(Xnod,iel)) return false; // Initialize element quantities - FiniteElement fe((**el).nBasisFunctions()); fe.iel = MLGE[iel-1]; - fe.xi = fe.eta = edgeDir < 0 ? -1.0 : 1.0; - LocalIntegral* A = integrand.getLocalIntegral(fe.N.size(),fe.iel,true); + fe.xi = fe.eta = edgeDir < 0 ? -1.0 : 1.0; + LocalIntegral* A = integrand.getLocalIntegral(MNPC[iel-1].size(), + fe.iel,true); if (!integrand.initElementBou(MNPC[iel-1],*A)) return false; // Get integration gauss points over this element @@ -1512,6 +1511,9 @@ bool ASMu2D::integrate (Integrand& integrand, int lIndex, if (edgeDir < 0) normal *= -1.0; + // Store tangent vectors in fe.G for shells + if (nsd > 2) fe.G = Jac; + // Cartesian coordinates of current integration point X.assign(Xnod * fe.N); X.t = time.t; @@ -1666,6 +1668,7 @@ bool ASMu2D::tesselate (ElementBlock& grid, const int* npe) const grid.unStructResize(nElements * nSubElPerElement, nElements * nNodesPerElement); + Go::Point pt; std::vector::iterator el; int inod = 0; int iel = 0; @@ -1676,15 +1679,12 @@ bool ASMu2D::tesselate (ElementBlock& grid, const int* npe) const double vmin = (**el).vmin(); double vmax = (**el).vmax(); for(int iv=0; ivpoint(pt, u,v, iel, iu!=npe[0]-1, iv!=npe[1]-1); + grid.setCoor(inod, SplineUtils::toVec3(pt,nsd)); grid.setParams(inod, u, v); - for(int dim=0; dimgetElementContaining(gpar[0][i],gpar[1][i]); - FiniteElement fe(lrspline->getElement(iel)->nBasisFunctions()); fe.iel = iel + 1; fe.u = gpar[0][i]; fe.v = gpar[1][i]; @@ -1881,19 +1881,19 @@ bool ASMu2D::evalSolution (Matrix& sField, const IntegrandBase& integrand, if (nPoints != gpar[1].size()) return false; + FiniteElement fe(0,firstIp); Vector solPt; Matrix dNdu, Jac, Xnod; Matrix3D d2Ndu2, Hess; // Evaluate the secondary solution field at each point - for (size_t i = 0; i < nPoints; i++) + for (size_t i = 0; i < nPoints; i++, fe.iGP++) { // Fetch element containing evaluation point // sadly, points are not always ordered in the same way as the elements int iel = lrspline->getElementContaining(gpar[0][i],gpar[1][i]); // Evaluate the basis functions at current parametric point - FiniteElement fe(lrspline->getElement(iel)->nBasisFunctions(),firstIp+i); if (use2ndDer) { Go::BasisDerivsSf2 spline; @@ -1918,6 +1918,9 @@ bool ASMu2D::evalSolution (Matrix& sField, const IntegrandBase& integrand, if (!utl::Hessian(Hess,fe.d2NdX2,Jac,Xnod,d2Ndu2,dNdu)) continue; + // Store tangent vectors in fe.G for shells + if (nsd > 2) fe.G = Jac; + // Now evaluate the solution field if (!integrand.evalSol(solPt,fe,Xnod*fe.N,MNPC[iel])) return false; diff --git a/src/Utility/CoordinateMapping.C b/src/Utility/CoordinateMapping.C index 4959bb1c..a117f6dc 100644 --- a/src/Utility/CoordinateMapping.C +++ b/src/Utility/CoordinateMapping.C @@ -36,6 +36,20 @@ Real utl::Jacobian (matrix& J, matrix& dNdX, // Compute the length of the tangent vector {X},u detJ = J.getColumn(1).norm2(); // |J| = sqrt(X,u*X,u + Y,u*Y,u + Z,u*Z,u) } + else if (J.cols() == 2 && J.rows() > 2) + { + // Special for two-parametric elements in 3-dimensional space (shell) + dNdX = dNdu; + // Compute the length of the normal vector {X},u x {X},v + Vec3 a1 = J.getColumn(1); + Vec3 a2 = J.getColumn(2); + detJ = Vec3(a1,a2).length(); + // Store the two tangent vectors in J + a1.normalize(); + a2.normalize(); + J.fillColumn(1,a1.ptr()); + J.fillColumn(2,a2.ptr()); + } else { // Compute the Jacobian determinant and inverse @@ -97,6 +111,15 @@ Real utl::Jacobian (matrix& J, Vec3& n, matrix& dNdX, dS = v2.normalize(); // dA = |v2| v3.normalize(); n.cross(v2,v3); // n = v2 x v3 / (|v2|*|v3|) + if (J.rows() > 2) + { + // Special treatment for two-parametric elements in 3D space (shells) + dNdX = dNdu; + v1.normalize(); + J.fillColumn(t1,v1.ptr()); + J.fillColumn(t2,v2.ptr()); + return dS; + } } else { @@ -135,7 +158,7 @@ bool utl::Hessian (matrix3d& H, matrix3d& d2NdX2, d2NdX2.resize(0,0,0,true); return true; } - else if (Ji.cols() == 1 && Ji.rows() > 1) + else if (Ji.cols() <= 2 && Ji.rows() > Ji.cols()) { // Special treatment for one-parametric elements in multi-dimension space d2NdX2 = d2Ndu2; @@ -181,13 +204,12 @@ void utl::getGmat (const matrix& Ji, const Real* du, matrix& G) size_t nsd = Ji.cols(); G.resize(nsd,nsd,true); + Real domain = pow(2.0,nsd); for (size_t k = 1; k <= nsd; k++) for (size_t l = 1; l <= nsd; l++) { - Real scale = Real(1)/(du[k-1]*du[l-1]); + Real scale = domain/(du[k-1]*du[l-1]); for (size_t m = 1; m <= nsd; m++) G(k,l) += Ji(m,k)*Ji(m,l)*scale; } - - G *= pow(2, nsd); } diff --git a/src/Utility/Test/TestCoordinateMapping.C b/src/Utility/Test/TestCoordinateMapping.C index 7e41d13e..3222e8d0 100644 --- a/src/Utility/Test/TestCoordinateMapping.C +++ b/src/Utility/Test/TestCoordinateMapping.C @@ -208,6 +208,51 @@ TEST(TestCoordinateMapping, Hessian2D) } +class SIMShell : public SIM2D +{ +public: + SIMShell() { nsd = 3; } + virtual ~SIMShell() {} +}; + + +TEST(TestCoordinateMapping, JacobianShell) +{ + SIMShell sim; + ASSERT_TRUE(sim.createDefaultModel()); + ASMs2D& p = static_cast(*sim.getPatch(1)); + ASSERT_TRUE(p.uniformRefine(0,3)); + ASSERT_TRUE(p.uniformRefine(1,3)); + ASSERT_TRUE(sim.createFEMmodel()); + + Vector N; + Matrix X, dNdU; + p.extractBasis(0.25, 0.25, N, dNdU); + p.getElementCoordinates(X, 2); + + Matrix J, dNdX; + Real det = utl::Jacobian(J, dNdX, X, dNdU, true); + + const double dndx[] = { + -4.0, 4.0, 0.0, 0.0, + -4.0, 0.0, 4.0, 0.0 + }; + + EXPECT_FLOAT_EQ(det, 1.0); + EXPECT_EQ(J.rows(), 3U); + EXPECT_EQ(J.cols(), 2U); + for (size_t i = 1; i <= J.rows(); i++) + for (size_t j = 1; j <= J.cols(); j++) + EXPECT_FLOAT_EQ(J(i,j), i == j ? 1.0 : 0.0); + + size_t i = 0; + EXPECT_EQ(dNdX.rows(), 4U); + EXPECT_EQ(dNdX.cols(), 2U); + for (double v : dNdX) + EXPECT_FLOAT_EQ(dndx[i++], v); +} + + TEST(TestCoordinateMapping, Hessian2D_mixed) { SIM2D sim({1,1});