diff --git a/src/ASM/ASMs2D.C b/src/ASM/ASMs2D.C index 7cc8cfb5..6a7556ee 100644 --- a/src/ASM/ASMs2D.C +++ b/src/ASM/ASMs2D.C @@ -1544,6 +1544,8 @@ bool ASMs2D::integrate (Integrand& integrand, for (size_t t = 0; t < threadGroups[g].size(); t++) { FiniteElement fe(p1*p2); + fe.p = p1 - 1; + fe.q = p2 - 1; Matrix dNdu, Xnod, Jac; Matrix3D d2Ndu2, Hess; double dXidu[2]; @@ -1707,8 +1709,12 @@ bool ASMs2D::integrate (Integrand& integrand, // Compute Hessian of coordinate mapping and 2nd order derivatives if (use2ndDer) + { if (!utl::Hessian(Hess,fe.d2NdX2,Jac,Xnod,d2Ndu2,fe.dNdX)) ok = false; + else if (nsd > 2) + utl::Hessian(Hess,fe.H); + } // Compute G-matrix if (integrand.getIntegrandType() & Integrand::G_MATRIX) @@ -1725,6 +1731,8 @@ bool ASMs2D::integrate (Integrand& integrand, std::cout <<"d2NdX2 ="<< fe.d2NdX2; if (!fe.G.empty()) std::cout <<"G ="<< fe.G; + if (!fe.H.empty()) + std::cout <<"H ="<< fe.H; } #endif @@ -1816,6 +1824,8 @@ bool ASMs2D::integrate (Integrand& integrand, for (size_t t = 0; t < threadGroups[g].size(); t++) { FiniteElement fe(p1*p2); + fe.p = p1 - 1; + fe.q = p2 - 1; Matrix dNdu, Xnod, Jac; Matrix3D d2Ndu2, Hess; double dXidu[2]; @@ -1901,8 +1911,12 @@ bool ASMs2D::integrate (Integrand& integrand, // Compute Hessian of coordinate mapping and 2nd order derivatives if (use2ndDer) + { if (!utl::Hessian(Hess,fe.d2NdX2,Jac,Xnod,d2Ndu2,fe.dNdX)) ok = false; + else if (nsd > 2) + utl::Hessian(Hess,fe.H); + } // Compute G-matrix if (integrand.getIntegrandType() & Integrand::G_MATRIX) @@ -2175,6 +2189,8 @@ bool ASMs2D::integrate (Integrand& integrand, int lIndex, size_t firstp = iit == firstBp.end() ? 0 : iit->second; FiniteElement fe(p1*p2); + fe.p = p1 - 1; + fe.q = p2 - 1; fe.xi = fe.eta = edgeDir < 0 ? -1.0 : 1.0; fe.u = gpar[0](1,1); fe.v = gpar[1](1,1); @@ -2617,9 +2633,11 @@ bool ASMs2D::evalSolution (Matrix& sField, const IntegrandBase& integrand, this->getNodalCoordinates(Xnod); FiniteElement fe(p1*p2,firstIp); - Vector solPt; - Matrix dNdu, Jac; - Matrix3D d2Ndu2, Hess; + fe.p = p1 - 1; + fe.q = p2 - 1; + Vector solPt; + Matrix dNdu, Jac; + Matrix3D d2Ndu2, Hess; // Evaluate the secondary solution field at each point for (size_t i = 0; i < nPoints; i++, fe.iGP++) @@ -2654,8 +2672,12 @@ bool ASMs2D::evalSolution (Matrix& sField, const IntegrandBase& integrand, // Compute Hessian of coordinate mapping and 2nd order derivatives if (use2ndDer) + { if (!utl::Hessian(Hess,fe.d2NdX2,Jac,Xtmp,d2Ndu2,fe.dNdX)) continue; + else if (nsd > 2) + utl::Hessian(Hess,fe.H); + } // Store tangent vectors in fe.G for shells if (nsd > 2) fe.G = Jac; diff --git a/src/ASM/FiniteElement.C b/src/ASM/FiniteElement.C index d12f44d4..2644a838 100644 --- a/src/ASM/FiniteElement.C +++ b/src/ASM/FiniteElement.C @@ -16,7 +16,8 @@ std::ostream& FiniteElement::write (std::ostream& os) const { - os <<"FiniteElement: iel="<< iel <<" iGP="<< iGP <<" p="<< p + os <<"FiniteElement: iel="<< iel <<" iGP="<< iGP + <<"\n p, q, r: "<< p <<" "<< q <<" "<< r <<"\n u, v, w: "<< u <<" "<< v <<" "<< w <<"\n xi, eta, zeta: "<< xi <<" "<< eta <<" "<< zeta <<"\n h, detJxW: "<< h <<" "<< detJxW << std::endl; @@ -24,6 +25,7 @@ std::ostream& FiniteElement::write (std::ostream& os) const if (!dNdX.empty()) os <<"dNdX:"<< dNdX; if (!d2NdX2.empty()) os <<"d2NdX2:"<< d2NdX2; if (!G.empty()) os <<"G:"<< G; + if (!H.empty()) os <<"H:"<< H; if (!Navg.empty()) os <<"Navg:"<< Navg; if (!Xn.empty()) os <<"Xn:"<< Xn; if (!Te.isZero(0.0)) os <<"Te:\n"<< Te; diff --git a/src/ASM/FiniteElement.h b/src/ASM/FiniteElement.h index 94088974..fda71bb7 100644 --- a/src/ASM/FiniteElement.h +++ b/src/ASM/FiniteElement.h @@ -28,19 +28,19 @@ class FiniteElement public: //! \brief Default constructor. explicit FiniteElement(size_t n = 0, size_t i = 0) : iGP(i), N(n), Te(3) - { iel = p = 0; u = v = w = xi = eta = zeta = h = 0.0; detJxW = 1.0; } + { iel = p = q = r = 0; u = v = w = xi = eta = zeta = h = 0.0; detJxW = 1.0; } //! \brief Empty destructor. virtual ~FiniteElement() {} + //! \brief Returns the number of bases. + virtual size_t getNoBasis() const { return 1; } + //! \brief Returns a const reference to the basis function values. virtual const Vector& basis(char) const { return N; } //! \brief Returns a reference to the basis function values. virtual Vector& basis(char) { return N; } - //! \brief Returns number of bases - virtual size_t getNoBasis() const { return 1; } - //! \brief Returns a const reference to the basis function derivatives. virtual const Matrix& grad(char) const { return dNdX; } //! \brief Returns a reference to the basis function derivatives. @@ -74,11 +74,14 @@ public: Vector N; //!< Basis function values Matrix dNdX; //!< First derivatives (gradient) of the basis functions Matrix3D d2NdX2; //!< Second derivatives of the basis functions - Matrix G; //!< Matrix used for stabilized methods + Matrix G; //!< Covariant basis / Matrix used for stabilized methods + Matrix H; //!< Hessian // Element quantities int iel; //!< Element identifier - short int p; //!< Polynomial order of the basis functions + short int p; //!< Polynomial order of the basis in u-direction + short int q; //!< Polynomial order of the basis in v-direction + short int r; //!< Polynomial order of the basis in r-direction double h; //!< Characteristic element size/diameter Vec3Vec XC; //!< Array with element corner coordinate vectors Vector Navg; //!< Volume-averaged basis function values @@ -101,14 +104,14 @@ public: //! \brief Empty destructor. virtual ~MxFiniteElement() {} + //! \brief Returns the number of bases. + virtual size_t getNoBasis() const { return 1+Nx.size(); } + //! \brief Returns a const reference to the basis function values. virtual const Vector& basis(char b) const { return b == 1 ? N : Nx[b-2]; } //! \brief Returns a reference to the basis function values. virtual Vector& basis(char b) { return b == 1 ? N : Nx[b-2]; } - //! \brief Returns number of bases - virtual size_t getNoBasis() const { return 1+Nx.size(); } - //! \brief Returns a const reference to the basis function derivatives. virtual const Matrix& grad(char b) const { return b == 1 ? dNdX : dNxdX[b-2]; } //! \brief Returns a reference to the basis function derivatives. diff --git a/src/ASM/LR/ASMu2D.C b/src/ASM/LR/ASMu2D.C index 454e2298..875ae2d1 100644 --- a/src/ASM/LR/ASMu2D.C +++ b/src/ASM/LR/ASMu2D.C @@ -1043,6 +1043,8 @@ bool ASMu2D::integrate (Integrand& integrand, FiniteElement fe; fe.iel = MLGE[iel-1]; + fe.p = lrspline->order(0) - 1; + fe.q = lrspline->order(1) - 1; Matrix dNdu, Xnod, Jac; Matrix3D d2Ndu2, Hess; double dXidu[2]; @@ -1190,11 +1192,15 @@ bool ASMu2D::integrate (Integrand& integrand, // Compute Hessian of coordinate mapping and 2nd order derivatives if (integrand.getIntegrandType() & Integrand::SECOND_DERIVATIVES) + { if (!utl::Hessian(Hess,fe.d2NdX2,Jac,Xnod,d2Ndu2,dNdu)) { ok = false; continue; } + else if (nsd > 2) + utl::Hessian(Hess,fe.H); + } // Compute G-matrix if (integrand.getIntegrandType() & Integrand::G_MATRIX) @@ -1281,6 +1287,8 @@ bool ASMu2D::integrate (Integrand& integrand, FiniteElement fe; fe.iel = MLGE[iel-1]; + fe.p = lrspline->order(0) - 1; + fe.q = lrspline->order(1) - 1; Matrix dNdu, Xnod, Jac; Matrix3D d2Ndu2, Hess; Vec4 X; @@ -1348,11 +1356,15 @@ bool ASMu2D::integrate (Integrand& integrand, // Compute Hessian of coordinate mapping and 2nd order derivatives if (integrand.getIntegrandType() & Integrand::SECOND_DERIVATIVES) + { if (!utl::Hessian(Hess,fe.d2NdX2,Jac,Xnod,d2Ndu2,dNdu)) { ok = false; continue; } + else if (nsd > 2) + utl::Hessian(Hess,fe.H); + } // Store tangent vectors in fe.G for shells if (nsd > 2) fe.G = Jac; @@ -1473,6 +1485,8 @@ bool ASMu2D::integrate (Integrand& integrand, int lIndex, // Initialize element quantities fe.iel = MLGE[iel-1]; + fe.p = lrspline->order(0) - 1; + fe.q = lrspline->order(1) - 1; fe.xi = fe.eta = edgeDir < 0 ? -1.0 : 1.0; LocalIntegral* A = integrand.getLocalIntegral(MNPC[iel-1].size(), fe.iel,true); @@ -1742,6 +1756,8 @@ bool ASMu2D::evalSolution (Matrix& sField, const Vector& locSol, return false; FiniteElement fe; + fe.p = lrspline->order(0) - 1; + fe.q = lrspline->order(1) - 1; Vector ptSol; Matrix dNdu, dNdX, Jac, Xnod, eSol, ptDer; Matrix3D d2Ndu2, d2NdX2, Hess, ptDer2; @@ -1882,6 +1898,8 @@ bool ASMu2D::evalSolution (Matrix& sField, const IntegrandBase& integrand, return false; FiniteElement fe(0,firstIp); + fe.p = lrspline->order(0) - 1; + fe.q = lrspline->order(1) - 1; Vector solPt; Matrix dNdu, Jac, Xnod; Matrix3D d2Ndu2, Hess; @@ -1915,8 +1933,12 @@ bool ASMu2D::evalSolution (Matrix& sField, const IntegrandBase& integrand, // Compute Hessian of coordinate mapping and 2nd order derivatives if (use2ndDer) + { if (!utl::Hessian(Hess,fe.d2NdX2,Jac,Xnod,d2Ndu2,dNdu)) continue; + else if (nsd > 2) + utl::Hessian(Hess,fe.H); + } // Store tangent vectors in fe.G for shells if (nsd > 2) fe.G = Jac; diff --git a/src/Utility/CoordinateMapping.C b/src/Utility/CoordinateMapping.C index a117f6dc..dbc155f7 100644 --- a/src/Utility/CoordinateMapping.C +++ b/src/Utility/CoordinateMapping.C @@ -41,14 +41,7 @@ Real utl::Jacobian (matrix& J, matrix& dNdX, // 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()); + detJ = Vec3(J.getColumn(1),J.getColumn(2)).length(); } else { @@ -115,9 +108,6 @@ Real utl::Jacobian (matrix& J, Vec3& n, matrix& dNdX, { // 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; } } @@ -161,6 +151,7 @@ bool utl::Hessian (matrix3d& H, matrix3d& d2NdX2, else if (Ji.cols() <= 2 && Ji.rows() > Ji.cols()) { // Special treatment for one-parametric elements in multi-dimension space + // as well as two-parametric elements in 3D space (shells) d2NdX2 = d2Ndu2; return true; } @@ -199,6 +190,24 @@ bool utl::Hessian (matrix3d& H, matrix3d& d2NdX2, } +void utl::Hessian (const matrix3d& Hess, matrix& H) +{ + size_t i, n = Hess.dim(2); + if (n == Hess.dim(3)) + { + H.resize(Hess.dim(1),n*(n+1)/2); + for (i = 1; i <= n; i++) + H.fillColumn(i,Hess.getColumn(i,i)); + for (i = 2; i <= n; i++) + H.fillColumn(n+i-1,Hess.getColumn(1,i)); + if (n > 2) + H.fillColumn(n+n,Hess.getColumn(2,3)); + } + else + H.clear(); +} + + void utl::getGmat (const matrix& Ji, const Real* du, matrix& G) { size_t nsd = Ji.cols(); diff --git a/src/Utility/CoordinateMapping.h b/src/Utility/CoordinateMapping.h index 7108d315..df86b1bc 100644 --- a/src/Utility/CoordinateMapping.h +++ b/src/Utility/CoordinateMapping.h @@ -64,13 +64,17 @@ namespace utl //! \param[in] X Matrix of element nodal coordinates //! \param[in] d2Ndu2 Second order derivatives of basis functions //! \param[in] dNdX First order derivatives of basis functions - //! \param[in] geoMapping If true, calculate geometry mapping + //! \param[in] geoMapping If \e true, calculate geometry mapping //! \return \e false if matrix dimensions are incompatible, otherwise \e true - //! \details If geoMapping is true, H is output, else H is input and assumed to be already calculated. + //! + //! \details If geoMapping is \e true, H is output, + //! else H is input and assumed to be already calculated in a previous call. bool Hessian(matrix3d& H, matrix3d& d2NdX2, const matrix& Ji, const matrix& X, const matrix3d& d2Ndu2, const matrix& dNdX, - bool geoMapping=true); + bool geoMapping = true); + //! \brief Convert a Hessian from a matrix3d to a matrix assuming symmetry. + void Hessian(const matrix3d& Hess, matrix& H); //! \brief Compute the stabilization matrix \b G from the Jacobian inverse. //! \param[in] Ji The inverse of the Jacobian matrix diff --git a/src/Utility/Test/TestCoordinateMapping.C b/src/Utility/Test/TestCoordinateMapping.C index 3222e8d0..f33d7e08 100644 --- a/src/Utility/Test/TestCoordinateMapping.C +++ b/src/Utility/Test/TestCoordinateMapping.C @@ -221,35 +221,58 @@ TEST(TestCoordinateMapping, JacobianShell) SIMShell sim; ASSERT_TRUE(sim.createDefaultModel()); ASMs2D& p = static_cast(*sim.getPatch(1)); + ASSERT_TRUE(p.raiseOrder(1,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); + Matrix3D d2Ndu2; + p.extractBasis(0.25, 0.25, N, dNdU, d2Ndu2); p.getElementCoordinates(X, 2); - Matrix J, dNdX; + Matrix J, dNdX, Hs; + Matrix3D H, d2NdX2; Real det = utl::Jacobian(J, dNdX, X, dNdU, true); + utl::Hessian(H, d2NdX2, J, X, d2Ndu2, dNdX); + utl::Hessian(H, Hs); - const double dndx[] = { - -4.0, 4.0, 0.0, 0.0, - -4.0, 0.0, 4.0, 0.0 + const double j[] = { + 1.0, 0.0, 0.0, + 0.0, 0.5, 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); + const double h[] = { + 0.0, 0.0, 0.0, + 0.0, 2.0, 0.0, + 0.0, 0.0, 0.0 + }; + + const double dndx[] = { + -2.0, 2.0, 0.0,-2.0, 2.0, 0.0, 0.0, 0.0, 0.0, + -2.0,-2.0, 0.0, 2.0, 2.0, 0.0, 0.0, 0.0, 0.0 + }; + + EXPECT_FLOAT_EQ(det, 0.5); size_t i = 0; - EXPECT_EQ(dNdX.rows(), 4U); + EXPECT_EQ(J.rows(), 3U); + EXPECT_EQ(J.cols(), 2U); + for (double v : J) + EXPECT_FLOAT_EQ(j[i++], v); + + i = 0; + EXPECT_EQ(dNdX.rows(), 9U); EXPECT_EQ(dNdX.cols(), 2U); for (double v : dNdX) EXPECT_FLOAT_EQ(dndx[i++], v); + + i = 0; + EXPECT_EQ(Hs.rows(), 3U); + EXPECT_EQ(Hs.cols(), 3U); + for (double v : Hs) + EXPECT_FLOAT_EQ(h[i++], v); }