diff --git a/Apps/FiniteDefElasticity/NonlinearElasticity.C b/Apps/FiniteDefElasticity/NonlinearElasticity.C index 7397c5b2..23a06752 100644 --- a/Apps/FiniteDefElasticity/NonlinearElasticity.C +++ b/Apps/FiniteDefElasticity/NonlinearElasticity.C @@ -13,6 +13,7 @@ #include "NonlinearElasticity.h" #include "MaterialBase.h" +#include "FiniteElement.h" #include "Utilities.h" #include "Profiler.h" @@ -75,15 +76,15 @@ void NonlinearElasticity::setMode (SIM::SolutionMode mode) } -bool NonlinearElasticity::evalInt (LocalIntegral*& elmInt, double detJW, - const Vector& N, const Matrix& dNdX, +bool NonlinearElasticity::evalInt (LocalIntegral*& elmInt, + const FiniteElement& fe, const Vec3& X) const { PROFILE3("NonlinearEl::evalInt"); // Evaluate the kinematic quantities, F and E, at this point Tensor F(nsd); - if (!this->kinematics(dNdX,F,E)) + if (!this->kinematics(fe.dNdX,F,E)) return false; // Evaluate current tangent at this point, that is @@ -106,30 +107,30 @@ bool NonlinearElasticity::evalInt (LocalIntegral*& elmInt, double detJW, PROFILE4("stiff_TL_"); if (nsd == 3) if (fullCmat) - stiff_tl3d_(N.size(),detJW,dNdX.ptr(),F.ptr(), + stiff_tl3d_(fe.N.size(),fe.detJxW,fe.dNdX.ptr(),F.ptr(), Cmat.ptr(),eKm->ptr()); else - stiff_tl3d_isoel_(N.size(),detJW,dNdX.ptr(),F.ptr(), + stiff_tl3d_isoel_(fe.N.size(),fe.detJxW,fe.dNdX.ptr(),F.ptr(), Cmat(1,1),Cmat(1,2),Cmat(4,4),eKm->ptr()); else if (nsd == 2) if (fullCmat) - stiff_tl2d_(N.size(),detJW,dNdX.ptr(),F.ptr(), + stiff_tl2d_(fe.N.size(),fe.detJxW,fe.dNdX.ptr(),F.ptr(), Cmat.ptr(),eKm->ptr()); else - stiff_tl2d_isoel_(N.size(),detJW,dNdX.ptr(),F.ptr(), + stiff_tl2d_isoel_(fe.N.size(),fe.detJxW,fe.dNdX.ptr(),F.ptr(), Cmat(1,1),Cmat(1,2),Cmat(3,3),eKm->ptr()); else if (nsd == 1) - for (a = 1; a <= N.size(); a++) - for (b = 1; b <= N.size(); b++) - (*eKm)(a,b) += dNdX(a,1)*F(1,1)*Cmat(1,1)*F(1,1)*dNdX(b,1); + for (a = 1; a <= fe.N.size(); a++) + for (b = 1; b <= fe.N.size(); b++) + (*eKm)(a,b) += fe.dNdX(a,1)*F(1,1)*Cmat(1,1)*F(1,1)*fe.dNdX(b,1); #else // This is too costly, but is basically what is done in the fortran routines PROFILE4("dNdX^t*F^t*C*F*dNdX"); unsigned short int l, m, n; SymmTensor4 D(Cmat,nsd); // fourth-order material tensor Matrix& EM = *eKm; - for (a = 1; a <= N.size(); a++) - for (b = 1; b <= N.size(); b++) + for (a = 1; a <= fe.N.size(); a++) + for (b = 1; b <= fe.N.size(); b++) for (m = 1; m <= nsd; m++) for (n = 1; n <= nsd; n++) { @@ -138,39 +139,39 @@ bool NonlinearElasticity::evalInt (LocalIntegral*& elmInt, double detJW, for (j = 1; j <= nsd; j++) for (k = 1; k <= nsd; k++) for (l = 1; l <= nsd; l++) - km += dNdX(a,i)*F(m,j)*D(i,j,k,l)*F(n,k)*dNdX(b,l); + km += fe.dNdX(a,i)*F(m,j)*D(i,j,k,l)*F(n,k)*fe.dNdX(b,l); - EM(nsd*(a-1)+m,nsd*(b-1)+n) += km*detJW; + EM(nsd*(a-1)+m,nsd*(b-1)+n) += km*fe.detJxW; } #endif } if (eKg && haveStrains) // Integrate the geometric stiffness matrix - this->formKG(*eKg,dNdX,S,detJW); + this->formKG(*eKg,fe.dNdX,S,fe.detJxW); if (eM) // Integrate the mass matrix - this->formMassMatrix(*eM,N,X,detJW); + this->formMassMatrix(*eM,fe.N,X,fe.detJxW); if (iS && haveStrains) { // Integrate the internal forces Vector& ES = *iS; - for (a = 1; a <= N.size(); a++) + for (a = 1; a <= fe.N.size(); a++) for (k = 1; k <= nsd; k++) { double f = 0.0; for (i = 1; i <= nsd; i++) for (j = 1; j <= nsd; j++) - f -= dNdX(a,i)*F(k,j)*S(i,j); - ES(nsd*(a-1)+k) += f*detJW; + f -= fe.dNdX(a,i)*F(k,j)*S(i,j); + ES(nsd*(a-1)+k) += f*fe.detJxW; } } if (eS) // Integrate the load vector due to gravitation and other body forces - this->formBodyForce(*eS,N,X,detJW); + this->formBodyForce(*eS,fe.N,X,fe.detJxW); return this->getIntegralResult(elmInt); } diff --git a/Apps/FiniteDefElasticity/NonlinearElasticity.h b/Apps/FiniteDefElasticity/NonlinearElasticity.h index 5cce868a..84d5dfe0 100644 --- a/Apps/FiniteDefElasticity/NonlinearElasticity.h +++ b/Apps/FiniteDefElasticity/NonlinearElasticity.h @@ -49,12 +49,9 @@ public: //! \brief Evaluates the integrand at an interior point. //! \param elmInt The local integral object to receive the contributions - //! \param[in] detJW Jacobian determinant times integration point weight - //! \param[in] N Basis function values - //! \param[in] dNdX Basis function gradients + //! \param[in] fe Finite element data of current integration point //! \param[in] X Cartesian coordinates of current integration point - virtual bool evalInt(LocalIntegral*& elmInt, double detJW, - const Vector& N, const Matrix& dNdX, + virtual bool evalInt(LocalIntegral*& elmInt, const FiniteElement& fe, const Vec3& X) const; //! \brief Evaluates the secondary solution at a result point. diff --git a/Apps/FiniteDefElasticity/NonlinearElasticityTL.C b/Apps/FiniteDefElasticity/NonlinearElasticityTL.C index c4b39ccd..12eb130d 100644 --- a/Apps/FiniteDefElasticity/NonlinearElasticityTL.C +++ b/Apps/FiniteDefElasticity/NonlinearElasticityTL.C @@ -13,6 +13,7 @@ #include "NonlinearElasticity.h" #include "MaterialBase.h" +#include "FiniteElement.h" #include "ElmMats.h" #include "Tensor.h" #include "Vec3Oper.h" @@ -82,15 +83,15 @@ void NonlinearElasticityTL::setMode (SIM::SolutionMode mode) } -bool NonlinearElasticityTL::evalInt (LocalIntegral*& elmInt, double detJW, - const Vector& N, const Matrix& dNdX, +bool NonlinearElasticityTL::evalInt (LocalIntegral*& elmInt, + const FiniteElement& fe, const Vec3& X) const { // Evaluate the deformation gradient, F, and the Green-Lagrange strains, E, // and compute the nonlinear strain-displacement matrix, B, from dNdX and F Tensor F(nsd); SymmTensor E(nsd), S(nsd); - if (!this->kinematics(dNdX,F,E)) + if (!this->kinematics(fe.dNdX,F,E)) return false; // Evaluate the constitutive relation @@ -105,36 +106,36 @@ bool NonlinearElasticityTL::evalInt (LocalIntegral*& elmInt, double detJW, if (eKm) { // Integrate the material stiffness matrix - CB.multiply(Cmat,Bmat).multiply(detJW); // CB = C*B*|J|*w - eKm->multiply(Bmat,CB,true,false,true); // EK += B^T * CB + CB.multiply(Cmat,Bmat).multiply(fe.detJxW); // CB = C*B*|J|*w + eKm->multiply(Bmat,CB,true,false,true); // EK += B^T * CB } if (eKg && lHaveStrains) // Integrate the geometric stiffness matrix - this->formKG(*eKg,dNdX,S,detJW); + this->formKG(*eKg,fe.dNdX,S,fe.detJxW); if (eM) // Integrate the mass matrix - this->formMassMatrix(*eM,N,X,detJW); + this->formMassMatrix(*eM,fe.N,X,fe.detJxW); if (iS && lHaveStrains) { // Integrate the internal forces - S *= -detJW; + S *= -fe.detJxW; if (!Bmat.multiply(S,*iS,true,true)) // ES -= B^T*S return false; } if (eS) // Integrate the load vector due to gravitation and other body forces - this->formBodyForce(*eS,N,X,detJW); + this->formBodyForce(*eS,fe.N,X,fe.detJxW); return this->getIntegralResult(elmInt); } -bool NonlinearElasticityTL::evalBou (LocalIntegral*& elmInt, double detJW, - const Vector& N, const Matrix& dNdX, +bool NonlinearElasticityTL::evalBou (LocalIntegral*& elmInt, + const FiniteElement& fe, const Vec3& X, const Vec3& normal) const { if (!tracFld) @@ -163,7 +164,7 @@ bool NonlinearElasticityTL::evalBou (LocalIntegral*& elmInt, double detJW, // Compute the deformation gradient, F Tensor F(nsd); SymmTensor dummy(0); - if (!this->kinematics(dNdX,F,dummy)) return false; + if (!this->kinematics(fe.dNdX,F,dummy)) return false; // Compute its inverse and determinant, J double J = F.inverse(); @@ -180,9 +181,9 @@ bool NonlinearElasticityTL::evalBou (LocalIntegral*& elmInt, double detJW, // Integrate the force vector Vector& ES = *eS; - for (size_t a = 1; a <= N.size(); a++) + for (size_t a = 1; a <= fe.N.size(); a++) for (i = 1; i <= nsd; i++) - ES(nsd*(a-1)+i) += T[i-1]*N(a)*detJW; + ES(nsd*(a-1)+i) += T[i-1]*fe.N(a)*fe.detJxW; return this->getIntegralResult(elmInt); } diff --git a/Apps/FiniteDefElasticity/NonlinearElasticityTL.h b/Apps/FiniteDefElasticity/NonlinearElasticityTL.h index 81cd0a42..070b26a3 100644 --- a/Apps/FiniteDefElasticity/NonlinearElasticityTL.h +++ b/Apps/FiniteDefElasticity/NonlinearElasticityTL.h @@ -44,26 +44,20 @@ public: //! \brief Evaluates the integrand at an interior point. //! \param elmInt The local integral object to receive the contributions - //! \param[in] detJW Jacobian determinant times integration point weight - //! \param[in] N Basis function values - //! \param[in] dNdX Basis function gradients + //! \param[in] fe Finite element data of current integration point //! \param[in] X Cartesian coordinates of current integration point - virtual bool evalInt(LocalIntegral*& elmInt, double detJW, - const Vector& N, const Matrix& dNdX, + virtual bool evalInt(LocalIntegral*& elmInt, const FiniteElement& fe, const Vec3& X) const; //! \brief Evaluates the integrand at a boundary point. //! \param elmInt The local integral object to receive the contributions - //! \param[in] detJW Jacobian determinant times integration point weight - //! \param[in] N Basis function values - //! \param[in] dNdX Basis function gradients + //! \param[in] fe Finite element data of current integration point //! \param[in] X Cartesian coordinates of current integration point //! \param[in] normal Boundary normal vector at current integration point //! //! \details This method is reimplemented in this class to account for //! possibly with-rotated traction fields in the Total-Lagrangian setting. - virtual bool evalBou(LocalIntegral*& elmInt, double detJW, - const Vector& N, const Matrix& dNdX, + virtual bool evalBou(LocalIntegral*& elmInt, const FiniteElement& fe, const Vec3& X, const Vec3& normal) const; protected: diff --git a/Apps/FiniteDefElasticity/NonlinearElasticityUL.C b/Apps/FiniteDefElasticity/NonlinearElasticityUL.C index 7e6b10a4..05948a50 100644 --- a/Apps/FiniteDefElasticity/NonlinearElasticityUL.C +++ b/Apps/FiniteDefElasticity/NonlinearElasticityUL.C @@ -13,6 +13,7 @@ #include "NonlinearElasticityUL.h" #include "MaterialBase.h" +#include "FiniteElement.h" #include "ElmMats.h" #include "ElmNorm.h" #include "TimeDomain.h" @@ -101,14 +102,14 @@ void NonlinearElasticityUL::initResultPoints () bool NonlinearElasticityUL::evalInt (LocalIntegral*& elmInt, - const TimeDomain& prm, double detJW, - const Vector& N, const Matrix& dNdX, + const FiniteElement& fe, + const TimeDomain& prm, const Vec3& X) const { // Evaluate the deformation gradient, F, and the Green-Lagrange strains, E Tensor F(nsd); SymmTensor E(nsd); - if (!this->kinematics(dNdX,F,E)) + if (!this->kinematics(fe.dNdX,F,E)) return false; bool lHaveStrains = !E.isZero(); @@ -121,12 +122,12 @@ bool NonlinearElasticityUL::evalInt (LocalIntegral*& elmInt, if (J == 0.0) return false; // Scale with J=|F| since we are integrating over current configuration - detJW *= J; + const_cast(fe.detJxW) *= J; if (eKm || iS) { // Push-forward the basis function gradients to current configuration - dNdx.multiply(dNdX,Fi); // dNdx = dNdX * F^-1 + dNdx.multiply(fe.dNdX,Fi); // dNdx = dNdX * F^-1 // Compute the small-deformation strain-displacement matrix B from dNdx if (!this->formBmatrix(dNdx)) return false; #ifdef INT_DEBUG @@ -137,7 +138,7 @@ bool NonlinearElasticityUL::evalInt (LocalIntegral*& elmInt, } else if (eKm || iS) // Initial state, no deformation yet - if (!this->formBmatrix(dNdX)) return false; + if (!this->formBmatrix(fe.dNdX)) return false; // Evaluate the constitutive relation SymmTensor sigma(nsd); @@ -151,36 +152,36 @@ bool NonlinearElasticityUL::evalInt (LocalIntegral*& elmInt, if (eKm) { // Integrate the material stiffness matrix - CB.multiply(Cmat,Bmat).multiply(detJW); // CB = C*B*|J|*w - eKm->multiply(Bmat,CB,true,false,true); // EK += B^T * CB + CB.multiply(Cmat,Bmat).multiply(fe.detJxW); // CB = C*B*|J|*w + eKm->multiply(Bmat,CB,true,false,true); // EK += B^T * CB } if (eKg && lHaveStrains) // Integrate the geometric stiffness matrix - this->formKG(*eKg,dNdx,sigma,detJW); + this->formKG(*eKg,dNdx,sigma,fe.detJxW); if (eM) // Integrate the mass matrix - this->formMassMatrix(*eM,N,X,detJW); + this->formMassMatrix(*eM,fe.N,X,fe.detJxW); if (iS && lHaveStrains) { // Integrate the internal forces - sigma *= -detJW; + sigma *= -fe.detJxW; if (!Bmat.multiply(sigma,*iS,true,true)) // ES -= B^T*sigma return false; } if (eS) // Integrate the load vector due to gravitation and other body forces - this->formBodyForce(*eS,N,X,detJW); + this->formBodyForce(*eS,fe.N,X,fe.detJxW); return this->getIntegralResult(elmInt); } -bool NonlinearElasticityUL::evalBou (LocalIntegral*& elmInt, double detJW, - const Vector& N, const Matrix& dNdX, +bool NonlinearElasticityUL::evalBou (LocalIntegral*& elmInt, + const FiniteElement& fe, const Vec3& X, const Vec3& normal) const { if (!tracFld) @@ -207,7 +208,7 @@ bool NonlinearElasticityUL::evalBou (LocalIntegral*& elmInt, double detJW, { // Compute the deformation gradient, F Tensor F(nsd); - if (!this->formDefGradient(dNdX,F)) return false; + if (!this->formDefGradient(fe.dNdX,F)) return false; // Check for with-rotated pressure load if (tracFld->isNormalPressure()) @@ -226,14 +227,14 @@ bool NonlinearElasticityUL::evalBou (LocalIntegral*& elmInt, double detJW, } else // Scale with J=|F| since we are integrating over current configuration - detJW *= F.det(); + const_cast(fe.detJxW) *= F.det(); } // Integrate the force vector Vector& ES = *eS; - for (size_t a = 1; a <= N.size(); a++) + for (size_t a = 1; a <= fe.N.size(); a++) for (i = 1; i <= nsd; i++) - ES(nsd*(a-1)+i) += T[i-1]*N(a)*detJW; + ES(nsd*(a-1)+i) += T[i-1]*fe.N(a)*fe.detJxW; return this->getIntegralResult(elmInt); } @@ -323,8 +324,8 @@ void ElasticityNormUL::initIntegration (const TimeDomain& prm) bool ElasticityNormUL::evalInt (LocalIntegral*& elmInt, - const TimeDomain& prm, double detJW, - const Vector&, const Matrix& dNdX, + const FiniteElement& fe, + const TimeDomain& prm, const Vec3& X) const { ElmNorm& pnorm = ElasticityNorm::getElmNormBuffer(elmInt); @@ -332,9 +333,9 @@ bool ElasticityNormUL::evalInt (LocalIntegral*& elmInt, NonlinearElasticityUL* ulp = static_cast(&problem); // Evaluate the deformation gradient, F, and the Green-Lagrange strains, E - Tensor F(dNdX.cols()); - SymmTensor E(dNdX.cols()); - if (!ulp->kinematics(dNdX,F,E)) + Tensor F(fe.dNdX.cols()); + SymmTensor E(fe.dNdX.cols()); + if (!ulp->kinematics(fe.dNdX,F,E)) return false; // Evaluate the 2nd Piola Kirchhoff stresses, S @@ -346,7 +347,24 @@ bool ElasticityNormUL::evalInt (LocalIntegral*& elmInt, // Integrate the energy norm a(u^h,u^h) = Int_Omega0 (S:E) dV0 if (U == 0.0) U = S.innerProd(E); - pnorm[0] += U*detJW; + pnorm[0] += U*fe.detJxW; return true; } + + +bool ElasticityNormUL::evalIntMx (LocalIntegral*& elmInt, + const MxFiniteElement& fe, + const TimeDomain& prm, + const Vec3& X) const +{ + return this->evalInt(elmInt,fe,prm,X); +} + + +bool ElasticityNormUL::evalBouMx (LocalIntegral*& elmInt, + const MxFiniteElement& fe, + const Vec3& X, const Vec3& normal) const +{ + return this->evalBou(elmInt,fe,X,normal); +} diff --git a/Apps/FiniteDefElasticity/NonlinearElasticityUL.h b/Apps/FiniteDefElasticity/NonlinearElasticityUL.h index 6b959d07..e0acb95b 100644 --- a/Apps/FiniteDefElasticity/NonlinearElasticityUL.h +++ b/Apps/FiniteDefElasticity/NonlinearElasticityUL.h @@ -51,21 +51,15 @@ public: //! \brief Evaluates the integrand at an interior point. //! \param elmInt The local integral object to receive the contributions + //! \param[in] fe Finite element data of current integration point //! \param[in] prm Nonlinear solution algorithm parameters - //! \param[in] detJW Jacobian determinant times integration point weight - //! \param[in] N Basis function values - //! \param[in] dNdX Basis function gradients //! \param[in] X Cartesian coordinates of current integration point - virtual bool evalInt(LocalIntegral*& elmInt, - const TimeDomain& prm, double detJW, - const Vector& N, const Matrix& dNdX, - const Vec3& X) const; + virtual bool evalInt(LocalIntegral*& elmInt, const FiniteElement& fe, + const TimeDomain& prm, const Vec3& X) const; //! \brief Evaluates the integrand at a boundary point. //! \param elmInt The local integral object to receive the contributions - //! \param[in] detJW Jacobian determinant times integration point weight - //! \param[in] N Basis function values - //! \param[in] dNdX Basis function gradients + //! \param[in] fe Finite element data of current integration point //! \param[in] X Cartesian coordinates of current integration point //! \param[in] normal Boundary normal vector at current integration point //! @@ -73,8 +67,7 @@ public: //! possibly with-rotated traction fields (non-conservative loads). //! For uni-directional (conservative) loads, it is similar to the //! \a LinearElasticity::evalBou method. - virtual bool evalBou(LocalIntegral*& elmInt, double detJW, - const Vector& N, const Matrix& dNdX, + virtual bool evalBou(LocalIntegral*& elmInt, const FiniteElement& fe, const Vec3& X, const Vec3& normal) const; //! \brief Returns a pointer to an Integrand for solution norm evaluation. @@ -144,45 +137,26 @@ public: //! \brief Evaluates the integrand at an interior point. //! \param elmInt The local integral object to receive the contributions + //! \param[in] fe Finite element data of current integration point //! \param[in] prm Nonlinear solution algorithm parameters - //! \param[in] detJW Jacobian determinant times integration point weight - //! \param[in] N Basis function values - //! \param[in] dNdX Basis function gradients //! \param[in] X Cartesian coordinates of current integration point - virtual bool evalInt(LocalIntegral*& elmInt, - const TimeDomain& prm, double detJW, - const Vector& N, const Matrix& dNdX, - const Vec3& X) const; + virtual bool evalInt(LocalIntegral*& elmInt, const FiniteElement& fe, + const TimeDomain& prm, const Vec3& X) const; //! \brief Evaluates the integrand at an interior point (mixed). //! \param elmInt The local integral object to receive the contributions + //! \param[in] fe Mixed finite element data of current integration point //! \param[in] prm Nonlinear solution algorithm parameters - //! \param[in] detJW Jacobian determinant times integration point weight - //! \param[in] N1 Basis function values for the displacement field - //! \param[in] dN1dX Basis function gradients for the displacement field //! \param[in] X Cartesian coordinates of current integration point - virtual bool evalInt(LocalIntegral*& elmInt, - const TimeDomain& prm, double detJW, - const Vector& N1, const Vector&, - const Matrix& dN1dX, const Matrix&, - const Vec3& X) const - { - return this->evalInt(elmInt,prm,detJW,N1,dN1dX,X); - } + virtual bool evalIntMx(LocalIntegral*& elmInt, const MxFiniteElement& fe, + const TimeDomain& prm, const Vec3& X) const; //! \brief Evaluates the integrand at a boundary point (mixed). //! \param elmInt The local integral object to receive the contributions - //! \param[in] detJW Jacobian determinant times integration point weight - //! \param[in] N1 Basis function values for the displacement field - //! \param[in] dN1dX Basis function gradients for the displacement field + //! \param[in] fe Mixed finite element data of current integration point //! \param[in] X Cartesian coordinates of current integration point //! \param[in] normal Boundary normal vector at current integration point - virtual bool evalBou(LocalIntegral*& elmInt, double detJW, - const Vector& N1, const Vector&, - const Matrix& dN1dX, const Matrix&, - const Vec3& X, const Vec3& normal) const - { - return this->ElasticityNorm::evalBou(elmInt,detJW,N1,dN1dX,X,normal); - } + virtual bool evalBouMx(LocalIntegral*& elmInt, const MxFiniteElement& fe, + const Vec3& X, const Vec3& normal) const; }; #endif diff --git a/Apps/FiniteDefElasticity/NonlinearElasticityULMX.C b/Apps/FiniteDefElasticity/NonlinearElasticityULMX.C index b323f0e9..2efe7904 100644 --- a/Apps/FiniteDefElasticity/NonlinearElasticityULMX.C +++ b/Apps/FiniteDefElasticity/NonlinearElasticityULMX.C @@ -13,6 +13,7 @@ #include "NonlinearElasticityULMX.h" #include "MaterialBase.h" +#include "FiniteElement.h" #include "ElmMats.h" #include "ElmNorm.h" #include "TimeDomain.h" @@ -140,26 +141,25 @@ bool NonlinearElasticityULMX::initElement (const std::vector& MNPC, bool NonlinearElasticityULMX::evalInt (LocalIntegral*& elmInt, - const TimeDomain&, double detJW, - const Vector& N, const Matrix& dNdX, - const Vec3& X) const + const FiniteElement& fe, + const TimeDomain&, const Vec3& X) const { if (myMats->b.size() < 3) return false; ItgPtData& ptData = myData[iP++]; #if INT_DEBUG > 0 - std::cout <<"NonlinearElasticityUL::dNdX ="<< dNdX; + std::cout <<"NonlinearElasticityUL::dNdX ="<< fe.dNdX; #endif // Evaluate the deformation gradient, Fp, at previous configuration const_cast(this)->eV = &myMats->b[2]; - if (!this->formDefGradient(dNdX,ptData.Fp)) + if (!this->formDefGradient(fe.dNdX,ptData.Fp)) return false; // Evaluate the deformation gradient, F, at current configuration const_cast(this)->eV = &myMats->b[1]; - if (!this->formDefGradient(dNdX,ptData.F)) + if (!this->formDefGradient(fe.dNdX,ptData.F)) return false; if (nsd == 2) // In 2D we always assume plane strain so set F(3,3)=1 @@ -174,10 +174,10 @@ bool NonlinearElasticityULMX::evalInt (LocalIntegral*& elmInt, if (J == 0.0) return false; // Push-forward the basis function gradients to current configuration - ptData.dNdx.multiply(dNdX,Fi); // dNdx = dNdX * F^-1 + ptData.dNdx.multiply(fe.dNdX,Fi); // dNdx = dNdX * F^-1 ptData.X.assign(X); - ptData.detJW = detJW; + ptData.detJW = fe.detJxW; // Evaluate the pressure modes (generalized coordinates) Vec3 Xg = X-X0; @@ -190,15 +190,15 @@ bool NonlinearElasticityULMX::evalInt (LocalIntegral*& elmInt, if (Hh) // Integrate the Hh-matrix - Hh->outer_product(ptData.Phi,ptData.Phi*detJW,true); + Hh->outer_product(ptData.Phi,ptData.Phi*fe.detJxW,true); if (eM) // Integrate the mass matrix - this->formMassMatrix(*eM,N,X,J*detJW); + this->formMassMatrix(*eM,fe.N,X,J*fe.detJxW); if (eS) // Integrate the load vector due to gravitation and other body forces - this->formBodyForce(*eS,N,X,J*detJW); + this->formBodyForce(*eS,fe.N,X,J*fe.detJxW); return this->getIntegralResult(elmInt); } @@ -470,12 +470,10 @@ bool ElasticityNormULMX::initElement (const std::vector& MNPC, bool ElasticityNormULMX::evalInt (LocalIntegral*& elmInt, - const TimeDomain& prm, double detJW, - const Vector& N, const Matrix& dNdX, - const Vec3& X) const + const FiniteElement& fe, + const TimeDomain& td, const Vec3& X) const { - NonlinearElasticityULMX& elp = static_cast(problem); - return elp.evalInt(elmInt,prm,detJW,N,dNdX,X); + return static_cast(problem).evalInt(elmInt,fe,td,X); } diff --git a/Apps/FiniteDefElasticity/NonlinearElasticityULMX.h b/Apps/FiniteDefElasticity/NonlinearElasticityULMX.h index 852729bd..0d8f0352 100644 --- a/Apps/FiniteDefElasticity/NonlinearElasticityULMX.h +++ b/Apps/FiniteDefElasticity/NonlinearElasticityULMX.h @@ -63,20 +63,16 @@ public: //! \brief Evaluates the integrand at an interior point. //! \param elmInt The local integral object to receive the contributions + //! \param[in] fe Finite element data of current integration point //! \param[in] prm Nonlinear solution algorithm parameters - //! \param[in] detJW Jacobian determinant times integration point weight - //! \param[in] N Basis function values - //! \param[in] dNdX Basis function gradients //! \param[in] X Cartesian coordinates of current integration point //! //! \details This method mainly stores the integration point quantities //! depending on the basis function values in the internal member \a myData. //! The actual numerical integration of the tangent stiffness and internal //! forces is then performed by the \a finalizeElement method. - virtual bool evalInt(LocalIntegral*& elmInt, - const TimeDomain& prm, double detJW, - const Vector& N, const Matrix& dNdX, - const Vec3& X) const; + virtual bool evalInt(LocalIntegral*& elmInt, const FiniteElement& fe, + const TimeDomain& prm, const Vec3& X) const; //! \brief Finalizes the element matrices after the numerical integration. //! \details This method is used to implement the multiple integration point @@ -131,20 +127,16 @@ public: //! \brief Evaluates the integrand at an interior point. //! \param elmInt The local integral object to receive the contributions + //! \param[in] fe Finite element data of current integration point //! \param[in] prm Nonlinear solution algorithm parameters - //! \param[in] detJW Jacobian determinant times integration point weight - //! \param[in] N Basis function values - //! \param[in] dNdX Basis function gradients //! \param[in] X Cartesian coordinates of current integration point //! //! \details This method just forwards to the \a evalInt method of the //! corresponding NonlinearElasticityULMX object, for which this class is //! declared as friend such that it can access the data members. The actual //! norm integration us then performed by the \a finalizeElement method. - virtual bool evalInt(LocalIntegral*& elmInt, - const TimeDomain& prm, double detJW, - const Vector& N, const Matrix& dNdX, - const Vec3& X) const; + virtual bool evalInt(LocalIntegral*& elmInt, const FiniteElement& fe, + const TimeDomain& prm, const Vec3& X) const; //! \brief Finalizes the element norms after the numerical integration. //! \details This method is used to implement the multiple integration point diff --git a/Apps/FiniteDefElasticity/NonlinearElasticityULMixed.C b/Apps/FiniteDefElasticity/NonlinearElasticityULMixed.C index e6c8d9a4..96eb15f1 100644 --- a/Apps/FiniteDefElasticity/NonlinearElasticityULMixed.C +++ b/Apps/FiniteDefElasticity/NonlinearElasticityULMixed.C @@ -13,6 +13,7 @@ #include "NonlinearElasticityULMixed.h" #include "MaterialBase.h" +#include "FiniteElement.h" #include "ElmMats.h" #include "Utilities.h" @@ -290,25 +291,22 @@ bool NonlinearElasticityULMixed::initElementBou (const std::vector& MNPC1, } -bool NonlinearElasticityULMixed::evalInt (LocalIntegral*& elmInt, - const TimeDomain& prm, double detJW, - const Vector& N1, - const Vector& N2, - const Matrix& dNdX1, - const Matrix& dNdX2, - const Vec3& X) const +bool NonlinearElasticityULMixed::evalIntMx (LocalIntegral*& elmInt, + const MxFiniteElement& fe, + const TimeDomain& prm, + const Vec3& X) const { // Evaluate the deformation gradient, F, and the Green-Lagrange strains, E Tensor F(nsd); SymmTensor E(nsd); - if (!this->kinematics(dNdX1,F,E)) + if (!this->kinematics(fe.dN1dX,F,E)) return false; bool lHaveStrains = !E.isZero(); // Evaluate the volumetric change and pressure fields - double Theta = N2.dot(myMats->b[T]) + 1.0; - double Press = N2.dot(myMats->b[P]); + double Theta = fe.N2.dot(myMats->b[T]) + 1.0; + double Press = fe.N2.dot(myMats->b[P]); #if INT_DEBUG > 0 std::cout <<"NonlinearElasticityULMixed::b_theta ="<< myMats->b[T]; std::cout <<"NonlinearElasticityULMixed::b_p ="<< myMats->b[P]; @@ -316,7 +314,7 @@ bool NonlinearElasticityULMixed::evalInt (LocalIntegral*& elmInt, #endif // Compute the mixed integration point volume - double dVol = Theta*detJW; + double dVol = Theta*fe.detJxW; // Compute the mixed model deformation gradient, F_bar Fbar = F; // notice that F_bar always has dimension 3 @@ -335,20 +333,20 @@ bool NonlinearElasticityULMixed::evalInt (LocalIntegral*& elmInt, if (J == 0.0) return false; // Push-forward the basis function gradients to current configuration - dNdx.multiply(dNdX1,Fi); // dNdx = dNdX * F^-1 + dNdx.multiply(fe.dN1dX,Fi); // dNdx = dNdX * F^-1 #if INT_DEBUG > 0 - std::cout <<"NonlinearElasticityULMixed::dNdX ="<< dNdX1; + std::cout <<"NonlinearElasticityULMixed::dNdX ="<< fe.dN1dX; std::cout <<"NonlinearElasticityULMixed::dNdx ="<< dNdx; std::cout <<"NonlinearElasticityULMixed::Fbar ="<< Fbar; #endif if (eM) // Integrate the mass matrix - this->formMassMatrix(*eM,N1,X,J*detJW); + this->formMassMatrix(*eM,fe.N1,X,J*fe.detJxW); if (eS) // Integrate the load vector due to gravitation and other body forces - this->formBodyForce(*eS,N1,X,J*detJW); + this->formBodyForce(*eS,fe.N1,X,J*fe.detJxW); // Evaluate the constitutive relation SymmTensor Eps(3), Sig(3), Sigma(nsd); @@ -384,34 +382,34 @@ bool NonlinearElasticityULMixed::evalInt (LocalIntegral*& elmInt, // Integrate the material stiffness matrix Dmat *= dVol; if (nsd == 2) - acckm2d_(N1.size(),dNdx.ptr(),Dmat.ptr(),eKm->ptr()); + acckm2d_(fe.N1.size(),dNdx.ptr(),Dmat.ptr(),eKm->ptr()); else - acckm3d_(N1.size(),dNdx.ptr(),Dmat.ptr(),eKm->ptr()); + acckm3d_(fe.N1.size(),dNdx.ptr(),Dmat.ptr(),eKm->ptr()); #endif // Integrate the volumetric change and pressure tangent terms size_t a, b; unsigned short int i, j, k; - for (a = 1; a <= N1.size(); a++) - for (b = 1; b <= N2.size(); b++) + for (a = 1; a <= fe.N1.size(); a++) + for (b = 1; b <= fe.N2.size(); b++) for (i = 1; i <= nsd; i++) { - myMats->A[Kut](nsd*(a-1)+i,b) += dNdx(a,i)*N2(b) * Dmat(i,7); + myMats->A[Kut](nsd*(a-1)+i,b) += dNdx(a,i)*fe.N2(b) * Dmat(i,7); if (nsd == 2) - myMats->A[Kut](nsd*(a-1)+i,b) += dNdx(a,3-i)*N2(b) * Dmat(4,7); + myMats->A[Kut](nsd*(a-1)+i,b) += dNdx(a,3-i)*fe.N2(b) * Dmat(4,7); else if (i < 3) { j = i + 1; k = j%3 + 1; - myMats->A[Kut](nsd*(a-1)+i,b) += dNdx(a,3-i)*N2(b) * Dmat(4,7); - myMats->A[Kut](nsd*(a-1)+j,b) += dNdx(a,5-j)*N2(b) * Dmat(5,7); - myMats->A[Kut](nsd*(a-1)+k,b) += dNdx(a,4-k)*N2(b) * Dmat(6,7); + myMats->A[Kut](nsd*(a-1)+i,b) += dNdx(a,3-i)*fe.N2(b) * Dmat(4,7); + myMats->A[Kut](nsd*(a-1)+j,b) += dNdx(a,5-j)*fe.N2(b) * Dmat(5,7); + myMats->A[Kut](nsd*(a-1)+k,b) += dNdx(a,4-k)*fe.N2(b) * Dmat(6,7); } - myMats->A[Kup](nsd*(a-1)+i,b) += dNdx(a,i)*N2(b) * dVol; + myMats->A[Kup](nsd*(a-1)+i,b) += dNdx(a,i)*fe.N2(b) * dVol; } - myMats->A[Ktt].outer_product(N2,N2*Dmat(7,7),true); // Ktt += N2*N2^T * D77 - myMats->A[Ktp].outer_product(N2,N2*(-detJW),true); // Ktp -= N2*N2^T * detJW + myMats->A[Ktt].outer_product(fe.N2,fe.N2*Dmat(7,7),true); // += N2*N2^T*D77 + myMats->A[Ktp].outer_product(fe.N2,fe.N2*(-fe.detJxW),true); // -= N2*N2^T*|J| if (lHaveStrains) { @@ -420,13 +418,28 @@ bool NonlinearElasticityULMixed::evalInt (LocalIntegral*& elmInt, // Integrate the internal forces Sigma *= -dVol; - if (!Bmat.multiply(Sigma,myMats->b[Ru],true,true)) // Ru -= B^T * sigma + if (!Bmat.multiply(Sigma,myMats->b[Ru],true,true)) // -= B^T*sigma return false; // Integrate the volumetric change and pressure forces - myMats->b[Rt].add(N2,(Press - Bpres)*detJW); // Rt += N2 * (p-pBar) - myMats->b[Rp].add(N2,(Theta - J)*detJW); // Rp += N2 * (Theta-J) + myMats->b[Rt].add(fe.N2,(Press - Bpres)*fe.detJxW); // += N2*(p-pBar)*|J| + myMats->b[Rp].add(fe.N2,(Theta - J)*fe.detJxW); // += N2*(Theta-J)*|J| } return this->getIntegralResult(elmInt); } + + +/*! + The boundary integral is the same as that of the parent class. It does not + depend on the pressure and volumetric change fields. Thus, this call is + forwarded to the corresponding single-field method of the parent class. +*/ + +bool NonlinearElasticityULMixed::evalBouMx (LocalIntegral*& elmInt, + const MxFiniteElement& fe, + const Vec3& X, + const Vec3& normal) const +{ + return this->evalBou(elmInt,fe,X,normal); +} diff --git a/Apps/FiniteDefElasticity/NonlinearElasticityULMixed.h b/Apps/FiniteDefElasticity/NonlinearElasticityULMixed.h index 4ea67f39..bd67a224 100644 --- a/Apps/FiniteDefElasticity/NonlinearElasticityULMixed.h +++ b/Apps/FiniteDefElasticity/NonlinearElasticityULMixed.h @@ -70,37 +70,19 @@ public: //! \brief Evaluates the mixed field problem integrand at an interior point. //! \param elmInt The local integral object to receive the contributions + //! \param[in] fe Mixed finite element data of current integration point //! \param[in] prm Nonlinear solution algorithm parameters - //! \param[in] detJW Jacobian determinant times integration point weight - //! \param[in] N1 Basis function values, field 1 (displacements) - //! \param[in] N2 Basis function values, field 2 (pressure and volume change) - //! \param[in] dN1dX Basis function gradients, field 1 - //! \param[in] dN2dX Basis function gradients, field 2 //! \param[in] X Cartesian coordinates of current integration point - virtual bool evalInt(LocalIntegral*& elmInt, - const TimeDomain& prm, double detJW, - const Vector& N1, const Vector& N2, - const Matrix& dN1dX, const Matrix& dN2dX, - const Vec3& X) const; + virtual bool evalIntMx(LocalIntegral*& elmInt, const MxFiniteElement& fe, + const TimeDomain& prm, const Vec3& X) const; //! \brief Evaluates the integrand at a boundary point. //! \param elmInt The local integral object to receive the contributions - //! \param[in] detJW Jacobian determinant times integration point weight - //! \param[in] N1 Basis function values, field 1 - //! \param[in] dN1dX Basis function gradients, field 1 + //! \param[in] fe Mixed finite element data of current integration point //! \param[in] X Cartesian coordinates of current integration point //! \param[in] normal Boundary normal vector at current integration point - //! - //! \details The boundary integral is the same as that of the parent class. - //! It does not depend on the pressure and volumetric change fields. - //! Thus, this call is forwarded to the single-field parent class method. - virtual bool evalBou(LocalIntegral*& elmInt, double detJW, - const Vector& N1, const Vector&, - const Matrix& dN1dX, const Matrix&, - const Vec3& X, const Vec3& normal) const - { - return this->NonlinearElasticityUL::evalBou(elmInt,detJW,N1,dN1dX,X,normal); - } + virtual bool evalBouMx(LocalIntegral*& elmInt, const MxFiniteElement& fe, + const Vec3& X, const Vec3& normal) const; protected: mutable Tensor Fbar; //!< Mixed model deformation gradient diff --git a/src/ASM/ASMs1D.C b/src/ASM/ASMs1D.C index a10a42fa..1c5e1e20 100644 --- a/src/ASM/ASMs1D.C +++ b/src/ASM/ASMs1D.C @@ -16,6 +16,7 @@ #include "ASMs1D.h" #include "TimeDomain.h" +#include "FiniteElement.h" #include "GlobalIntegral.h" #include "IntegrandBase.h" #include "CoordinateMapping.h" @@ -467,22 +468,21 @@ bool ASMs1D::integrate (Integrand& integrand, const int p1 = curv->order(); const int n1 = curv->numCoefs(); - Vector gpar; int pm1 = p1 - 1; + int nCol = n1 - pm1; + Matrix gpar(nGauss,nCol); DoubleIter uit = curv->basis().begin() + pm1; double ucurr, uprev = *(uit++); - int nCol = n1 - pm1; - gpar.reserve(nGauss*nCol); for (int j = 1; j <= nCol; uit++, j++) { ucurr = *uit; - for (int i = 0; i < nGauss; i++) - gpar.push_back(0.5*((ucurr-uprev)*xg[i] + ucurr+uprev)); + for (int i = 1; i <= nGauss; i++) + gpar(i,j) = 0.5*((ucurr-uprev)*xg[i-1] + ucurr+uprev); uprev = ucurr; } - Vector N; - Matrix dNdu, dNdX, Xnod, Jac; + FiniteElement fe(p1); + Matrix dNdu, Xnod, Jac; Vec4 X; @@ -514,21 +514,22 @@ bool ASMs1D::integrate (Integrand& integrand, for (int i = 0; i < nGauss; i++) { - // Weight of current integration point - double weight = 0.5*dL*wg[i]; + // Parameter value of current integration point + fe.u = gpar(i+1,iel); // Compute basis functions and derivatives - this->extractBasis(gpar((iel-1)*nGauss+i+1),N,dNdu); + this->extractBasis(fe.u,fe.N,dNdu); // Compute derivatives in terms of physical co-ordinates - double detJ = utl::Jacobian(Jac,dNdX,Xnod,dNdu); + fe.detJxW = utl::Jacobian(Jac,fe.dNdX,Xnod,dNdu); // Cartesian coordinates of current integration point - X = Xnod * N; + X = Xnod * fe.N; X.t = time.t; // Evaluate the integrand and accumulate element contributions - if (!integrand.evalInt(elmInt,time,detJ*weight,N,dNdX,X)) + fe.detJxW *= 0.5*dL*wg[i]; + if (!integrand.evalInt(elmInt,fe,time,X)) return false; } @@ -584,16 +585,16 @@ bool ASMs1D::integrate (Integrand& integrand, int lIndex, LocalIntegral* elmInt = locInt.empty() ? 0 : locInt[MLGE[iel-1]-1]; // Evaluate basis functions and corresponding derivatives - Vector N; + FiniteElement fe; Matrix dNdu; - this->extractBasis(param,N,dNdu); + this->extractBasis(param,fe.N,dNdu); // Cartesian coordinates of current integration point - Vec4 X(Xnod*N,time.t); + Vec4 X(Xnod*fe.N,time.t); // Compute basis function derivatives Matrix Jac, dNdX; - utl::Jacobian(Jac,dNdX,Xnod,dNdu); + utl::Jacobian(Jac,fe.dNdX,Xnod,dNdu); // Set up the normal vector Vec3 normal; @@ -603,7 +604,7 @@ bool ASMs1D::integrate (Integrand& integrand, int lIndex, normal.x = copysign(1.0,Jac(1,1)); // Evaluate the integrand and accumulate element contributions - if (!integrand.evalBou(elmInt,time,1.0,N,dNdX,X,normal)) + if (!integrand.evalBou(elmInt,fe,time,X,normal)) return false; // Assembly of global system integral @@ -623,12 +624,14 @@ bool ASMs1D::getGridParameters (RealArray& prm, int nSegPerSpan) const } DoubleIter uit = curv->basis().begin(); - double ucurr, uprev = *(uit++); + double ucurr = 0.0, uprev = *(uit++); while (uit != curv->basis().end()) { ucurr = *(uit++); if (ucurr > uprev) - for (int i = 0; i < nSegPerSpan; i++) + if (nSegPerSpan == 1) + prm.push_back(uprev); + else for (int i = 0; i < nSegPerSpan; i++) { double xg = (double)(2*i-nSegPerSpan)/(double)nSegPerSpan; prm.push_back(0.5*(ucurr*(1.0+xg) + uprev*(1.0-xg))); @@ -636,7 +639,8 @@ bool ASMs1D::getGridParameters (RealArray& prm, int nSegPerSpan) const uprev = ucurr; } - prm.push_back(curv->basis().endparam()); + if (ucurr > prm.back()) + prm.push_back(ucurr); return true; } @@ -732,7 +736,7 @@ bool ASMs1D::evalSolution (Matrix& sField, const Integrand& integrand, if (npe) { // Compute parameter values of the result sampling points - DoubleVec gpar; + RealArray gpar; if (this->getGridParameters(gpar,npe[0]-1)) // Evaluate the secondary solution at all sampling points return this->evalSolution(sField,integrand,gpar); diff --git a/src/ASM/ASMs1DLag.C b/src/ASM/ASMs1DLag.C index 77fa63b0..84f10792 100644 --- a/src/ASM/ASMs1DLag.C +++ b/src/ASM/ASMs1DLag.C @@ -16,6 +16,7 @@ #include "ASMs1DLag.h" #include "Lagrange.h" #include "TimeDomain.h" +#include "FiniteElement.h" #include "GlobalIntegral.h" #include "IntegrandBase.h" #include "CoordinateMapping.h" @@ -150,11 +151,15 @@ bool ASMs1DLag::integrate (Integrand& integrand, const double* wg = GaussQuadrature::getWeight(nGauss); if (!xg || !wg) return false; + // Get parametric coordinates of the elements + RealArray gpar; + this->getGridParameters(gpar,1); + // Order of basis (order = degree + 1) const int p1 = curv->order(); - Vector N(p1); - Matrix dNdu, dNdX, Xnod, Jac; + FiniteElement fe(p1); + Matrix dNdu, Xnod, Jac; Vec4 X; @@ -162,45 +167,48 @@ bool ASMs1DLag::integrate (Integrand& integrand, const int nel = this->getNoElms(); for (int iel = 1; iel <= nel; iel++) + { + // Set up control point coordinates for current element + if (!this->getElementCoordinates(Xnod,iel)) return false; + + // Initialize element quantities + if (!integrand.initElement(MNPC[iel-1])) return false; + + // Caution: Unless locInt is empty, we assume it points to an array of + // LocalIntegral pointers, of length at least the number of elements in + // the model (as defined by the highest number in the MLGE array). + // If the array is shorter than this, expect a segmentation fault. + LocalIntegral* elmInt = locInt.empty() ? 0 : locInt[MLGE[iel-1]-1]; + + + // --- Integration loop over all Gauss points in each direction ------------ + + for (int i = 0; i < nGauss; i++) { - // Set up control point coordinates for current element - if (!this->getElementCoordinates(Xnod,iel)) return false; + // Parameter value of current integration point + fe.u = 0.5*(gpar[iel-1]*(1.0-xg[i]) + gpar[iel]*(1.0+xg[i])); - // Initialize element quantities - if (!integrand.initElement(MNPC[iel-1])) return false; + // Compute basis function derivatives at current integration point + if (!Lagrange::computeBasis(fe.N,dNdu,p1,xg[i])) + return false; - // Caution: Unless locInt is empty, we assume it points to an array of - // LocalIntegral pointers, of length at least the number of elements in - // the model (as defined by the highest number in the MLGE array). - // If the array is shorter than this, expect a segmentation fault. - LocalIntegral* elmInt = locInt.empty() ? 0 : locInt[MLGE[iel-1]-1]; + // Compute Jacobian inverse of coordinate mapping and derivatives + fe.detJxW = utl::Jacobian(Jac,fe.dNdX,Xnod,dNdu)*wg[i]; + // Cartesian coordinates of current integration point + X = Xnod * fe.N; + X.t = time.t; - // --- Integration loop over all Gauss points in each direction ---------- - - for (int i = 0; i < nGauss; i++) - { - // Compute basis function derivatives at current integration point - if (!Lagrange::computeBasis(N,dNdu,p1,xg[i])) - return false; - - // Compute Jacobian inverse of coordinate mapping and derivatives - double detJ = utl::Jacobian(Jac,dNdX,Xnod,dNdu); - - // Cartesian coordinates of current integration point - X = Xnod * N; - X.t = time.t; - - // Evaluate the integrand and accumulate element contributions - if (!integrand.evalInt(elmInt,time,detJ*wg[i],N,dNdX,X)) - return false; - } - - // Assembly of global system integral - if (!glInt.assemble(elmInt,MLGE[iel-1])) + // Evaluate the integrand and accumulate element contributions + if (!integrand.evalInt(elmInt,fe,time,X)) return false; } + // Assembly of global system integral + if (!glInt.assemble(elmInt,MLGE[iel-1])) + return false; + } + return true; } @@ -250,16 +258,15 @@ bool ASMs1DLag::integrate (Integrand& integrand, int lIndex, LocalIntegral* elmInt = locInt.empty() ? 0 : locInt[MLGE[iel-1]-1]; // Evaluate basis functions and corresponding derivatives - Vector N; - Matrix dNdu; - if (!Lagrange::computeBasis(N,dNdu,p1,x)) return false; + FiniteElement fe; + Matrix dNdu, Jac; + if (!Lagrange::computeBasis(fe.N,dNdu,p1,x)) return false; // Cartesian coordinates of current integration point - Vec4 X(Xnod*N,time.t); + Vec4 X(Xnod*fe.N,time.t); // Compute basis function derivatives - Matrix Jac, dNdX; - utl::Jacobian(Jac,dNdX,Xnod,dNdu); + utl::Jacobian(Jac,fe.dNdX,Xnod,dNdu); // Set up the normal vector Vec3 normal; @@ -269,7 +276,7 @@ bool ASMs1DLag::integrate (Integrand& integrand, int lIndex, normal.x = copysign(1.0,Jac(1,1)); // Evaluate the integrand and accumulate element contributions - if (!integrand.evalBou(elmInt,time,1.0,N,dNdX,X,normal)) + if (!integrand.evalBou(elmInt,fe,time,X,normal)) return false; // Assembly of global system integral diff --git a/src/ASM/ASMs1DSpec.C b/src/ASM/ASMs1DSpec.C index 21f77433..3f7a61b1 100644 --- a/src/ASM/ASMs1DSpec.C +++ b/src/ASM/ASMs1DSpec.C @@ -15,6 +15,7 @@ #include "ASMs1DSpec.h" #include "TimeDomain.h" +#include "FiniteElement.h" #include "GlobalIntegral.h" #include "IntegrandBase.h" #include "CoordinateMapping.h" @@ -83,8 +84,8 @@ bool ASMs1DSpec::integrate (Integrand& integrand, if (!Legendre::GL(wg1,xg1,n1)) return false; - Vector N(p1); - Matrix dNdu, dNdX, Xnod, Jac; + FiniteElement fe(p1); + Matrix dNdu, Xnod, Jac; Vec4 X; @@ -110,29 +111,27 @@ bool ASMs1DSpec::integrate (Integrand& integrand, for (int i = 0; i < n1; i++) { - // Weight of current integration point - double weight = wg1[i]; - // Compute basis function derivatives at current integration point if (nGauss < 1) { - N.fill(0.0); - N(i+1) = 1.0; + fe.N.fill(0.0); + fe.N(i+1) = 1.0; dNdu.fillColumn(1,D1.getRow(i+1)); } else - if (!Lagrange::computeBasis(N,dNdu,points1,xg1[i])) + if (!Lagrange::computeBasis(fe.N,dNdu,points1,xg1[i])) return false; // Compute Jacobian inverse of coordinate mapping and derivatives - double detJ = utl::Jacobian(Jac,dNdX,Xnod,dNdu); + fe.detJxW = utl::Jacobian(Jac,fe.dNdX,Xnod,dNdu); // Cartesian coordinates of current integration point - X = Xnod*N; + X = Xnod*fe.N; X.t = time.t; // Evaluate the integrand and accumulate element contributions - if (!integrand.evalInt(elmInt,time,detJ*weight,N,dNdX,X)) + fe.detJxW *= wg1[i]; + if (!integrand.evalInt(elmInt,fe,time,X)) return false; } diff --git a/src/ASM/ASMs2D.C b/src/ASM/ASMs2D.C index abaaf499..8b122612 100644 --- a/src/ASM/ASMs2D.C +++ b/src/ASM/ASMs2D.C @@ -17,6 +17,7 @@ #include "ASMs2D.h" #include "TimeDomain.h" +#include "FiniteElement.h" #include "GlobalIntegral.h" #include "IntegrandBase.h" #include "CoordinateMapping.h" @@ -819,10 +820,9 @@ bool ASMs2D::integrate (Integrand& integrand, const int n2 = surf->numCoefs_v(); const int nel1 = n1 - p1 + 1; - Vector N(p1*p2), Navg; - Matrix dNdu, dNdX, Xnod, Jac; - Matrix3D d2Ndu2, d2NdX2, Hess; - double h = 0.0; + FiniteElement fe(p1*p2); + Matrix dNdu, Xnod, Jac; + Matrix3D d2Ndu2, Hess; Vec4 X; @@ -843,33 +843,33 @@ bool ASMs2D::integrate (Integrand& integrand, // Compute characteristic element length, if needed if (integrand.getIntegrandType() == 2) - h = getElmSize(p1,p2,Xnod); + fe.h = getElmSize(p1,p2,Xnod); else if (integrand.getIntegrandType() == 3) { // --- Compute average value of basis functions over the element ------- - Navg.resize(p1*p2,true); + fe.Navg.resize(p1*p2,true); double area = 0.0; int ip = ((i2-p2)*nGauss*nel1 + i1-p1)*nGauss; for (int j = 0; j < nGauss; j++, ip += nGauss*(nel1-1)) for (int i = 0; i < nGauss; i++, ip++) { // Fetch basis function derivatives at current integration point - extractBasis(spline[ip],N,dNdu); + extractBasis(spline[ip],fe.N,dNdu); // Compute Jacobian determinant of coordinate mapping // and multiply by weight of current integration point + double detJac = utl::Jacobian(Jac,fe.dNdX,Xnod,dNdu,false); double weight = 0.25*dA*wg[i]*wg[j]; - double detJxW = utl::Jacobian(Jac,dNdX,Xnod,dNdu,false)*weight; // Numerical quadrature - Navg += N*detJxW; - area += detJxW; + fe.Navg.add(fe.N,detJac*weight); + area += detJac*weight; } // Divide by element area - Navg /= area; + fe.Navg /= area; } else if (integrand.getIntegrandType() == 4) @@ -899,47 +899,38 @@ bool ASMs2D::integrate (Integrand& integrand, for (int j = 0; j < nGauss; j++, ip += nGauss*(nel1-1)) for (int i = 0; i < nGauss; i++, ip++) { - // Weight of current integration point - double weight = 0.25*dA*wg[i]*wg[j]; + // Parameter values of current integration point + fe.u = gpar[0](i+1,i1-p1+1); + fe.v = gpar[1](j+1,i2-p2+1); // Fetch basis function derivatives at current integration point if (integrand.getIntegrandType() == 2) - extractBasis(spline2[ip],N,dNdu,d2Ndu2); + extractBasis(spline2[ip],fe.N,dNdu,d2Ndu2); else - extractBasis(spline[ip],N,dNdu); + extractBasis(spline[ip],fe.N,dNdu); // Compute Jacobian inverse of coordinate mapping and derivatives - double detJ = utl::Jacobian(Jac,dNdX,Xnod,dNdu); - if (detJ == 0.0) continue; // skip singular points + fe.detJxW = utl::Jacobian(Jac,fe.dNdX,Xnod,dNdu); + if (fe.detJxW == 0.0) continue; // skip singular points // Compute Hessian of coordinate mapping and 2nd order derivatives if (integrand.getIntegrandType() == 2) - if (!utl::Hessian(Hess,d2NdX2,Jac,Xnod,d2Ndu2,dNdu)) + if (!utl::Hessian(Hess,fe.d2NdX2,Jac,Xnod,d2Ndu2,dNdu)) return false; #if SP_DEBUG > 4 std::cout <<"\niel, ip = "<< iel <<" "<< ip - <<"\nN ="<< N <<"dNdX ="<< dNdX << std::endl; + <<"\nN ="<< fe.N <<"dNdX ="<< fe.dNdX << std::endl; #endif // Cartesian coordinates of current integration point - X = Xnod * N; + X = Xnod * fe.N; X.t = time.t; // Evaluate the integrand and accumulate element contributions - bool ok = false; - switch (integrand.getIntegrandType()) - { - case 2: - ok = integrand.evalInt(elmInt,time,detJ*weight,N,dNdX,d2NdX2,X,h); - break; - case 3: - ok = integrand.evalInt(elmInt,time,detJ*weight,N,dNdX,Navg,X); - break; - default: - ok = integrand.evalInt(elmInt,time,detJ*weight,N,dNdX,X); - } - if (!ok) return false; + fe.detJxW *= 0.25*dA*wg[i]*wg[j]; + if (!integrand.evalInt(elmInt,fe,time,X)) + return false; } // Finalize the element quantities @@ -1013,8 +1004,8 @@ bool ASMs2D::integrate (Integrand& integrand, int lIndex, const int n1 = surf->numCoefs_u(); const int n2 = surf->numCoefs_v(); - Vector N(p1*p2); - Matrix dNdu, dNdX, Xnod, Jac; + FiniteElement fe(p1*p2); + Matrix dNdu, Xnod, Jac; Vec4 X; Vec3 normal; @@ -1060,24 +1051,22 @@ bool ASMs2D::integrate (Integrand& integrand, int lIndex, int ip = (abs(edgeDir) == 1 ? i2-p2 : i1-p1)*nGauss; for (int i = 0; i < nGauss; i++, ip++) { - // Weight of current integration point - double weight = 0.5*dS*wg[i]; - // Fetch basis function derivatives at current integration point - extractBasis(spline[ip],N,dNdu); + extractBasis(spline[ip],fe.N,dNdu); // Compute basis function derivatives and the edge normal - double dT = utl::Jacobian(Jac,normal,dNdX,Xnod,dNdu,t1,t2); - if (dT == 0.0) continue; // skip singular points + fe.detJxW = utl::Jacobian(Jac,normal,fe.dNdX,Xnod,dNdu,t1,t2); + if (fe.detJxW == 0.0) continue; // skip singular points if (edgeDir < 0) normal *= -1.0; // Cartesian coordinates of current integration point - X = Xnod * N; + X = Xnod * fe.N; X.t = time.t; // Evaluate the integrand and accumulate element contributions - if (!integrand.evalBou(elmInt,time,dT*weight,N,dNdX,X,normal)) + fe.detJxW *= 0.5*dS*wg[i]; + if (!integrand.evalBou(elmInt,fe,time,X,normal)) return false; } @@ -1102,12 +1091,14 @@ bool ASMs2D::getGridParameters (RealArray& prm, int dir, int nSegPerSpan) const } DoubleIter uit = surf->basis(dir).begin(); - double ucurr, uprev = *(uit++); + double ucurr = 0.0, uprev = *(uit++); while (uit != surf->basis(dir).end()) { ucurr = *(uit++); if (ucurr > uprev) - for (int i = 0; i < nSegPerSpan; i++) + if (nSegPerSpan == 1) + prm.push_back(uprev); + else for (int i = 0; i < nSegPerSpan; i++) { double xg = (double)(2*i-nSegPerSpan)/(double)nSegPerSpan; prm.push_back(0.5*(ucurr*(1.0+xg) + uprev*(1.0-xg))); @@ -1115,7 +1106,8 @@ bool ASMs2D::getGridParameters (RealArray& prm, int dir, int nSegPerSpan) const uprev = ucurr; } - prm.push_back(surf->basis(dir).endparam()); + if (ucurr > prm.back()) + prm.push_back(ucurr); return true; } @@ -1236,7 +1228,7 @@ bool ASMs2D::evalSolution (Matrix& sField, const Integrand& integrand, if (npe) { // Compute parameter values of the result sampling points - DoubleVec gpar[2]; + RealArray gpar[2]; for (int dir = 0; dir < 2 && retVal; dir++) retVal = this->getGridParameters(gpar[dir],dir,npe[dir]-1); diff --git a/src/ASM/ASMs2DLag.C b/src/ASM/ASMs2DLag.C index ce168916..6ed962df 100644 --- a/src/ASM/ASMs2DLag.C +++ b/src/ASM/ASMs2DLag.C @@ -16,6 +16,7 @@ #include "ASMs2DLag.h" #include "Lagrange.h" #include "TimeDomain.h" +#include "FiniteElement.h" #include "GlobalIntegral.h" #include "IntegrandBase.h" #include "CoordinateMapping.h" @@ -176,19 +177,29 @@ bool ASMs2DLag::integrate (Integrand& integrand, const double* wg = GaussQuadrature::getWeight(nGauss); if (!xg || !wg) return false; + // Get parametric coordinates of the elements + RealArray upar, vpar; + this->getGridParameters(upar,0,1); + this->getGridParameters(vpar,1,1); + + // Number of elements in each direction + const int nelx = upar.size() - 1; + const int nely = vpar.size() - 1; + // Order of basis in the two parametric directions (order = degree + 1) const int p1 = surf->order_u(); const int p2 = surf->order_v(); - Vector N(p1*p2); - Matrix dNdu, dNdX, Xnod, Jac; + FiniteElement fe(p1*p2); + Matrix dNdu, Xnod, Jac; Vec4 X; // === Assembly loop over all elements in the patch ========================== - const int nel = this->getNoElms(); - for (int iel = 1; iel <= nel; iel++) + int iel = 1; + for (int i2 = 0; i2 < nely; i2++) + for (int i1 = 0; i1 < nelx; i1++, iel++) { // Set up control point coordinates for current element if (!this->getElementCoordinates(Xnod,iel)) return false; @@ -219,24 +230,26 @@ bool ASMs2DLag::integrate (Integrand& integrand, for (int j = 0; j < nGauss; j++) for (int i = 0; i < nGauss; i++) { - // Weight of current integration point - double weight = wg[i]*wg[j]; + // Parameter value of current integration point + fe.u = 0.5*(upar[i1]*(1.0-xg[i]) + upar[i1+1]*(1.0+xg[i])); + fe.v = 0.5*(vpar[i2]*(1.0-xg[j]) + vpar[i2+1]*(1.0+xg[j])); // Compute basis function derivatives at current integration point // using tensor product of one-dimensional Lagrange polynomials - if (!Lagrange::computeBasis(N,dNdu,p1,xg[i],p2,xg[j])) + if (!Lagrange::computeBasis(fe.N,dNdu,p1,xg[i],p2,xg[j])) return false; // Compute Jacobian inverse of coordinate mapping and derivatives - double dA = utl::Jacobian(Jac,dNdX,Xnod,dNdu); - if (dA == 0.0) continue; // skip singular points + fe.detJxW = utl::Jacobian(Jac,fe.dNdX,Xnod,dNdu); + if (fe.detJxW == 0.0) continue; // skip singular points // Cartesian coordinates of current integration point - X = Xnod * N; + X = Xnod * fe.N; X.t = time.t; // Evaluate the integrand and accumulate element contributions - if (!integrand.evalInt(elmInt,time,dA*weight,N,dNdX,X)) + fe.detJxW *= wg[i]*wg[j]; + if (!integrand.evalInt(elmInt,fe,time,X)) return false; } @@ -279,8 +292,8 @@ bool ASMs2DLag::integrate (Integrand& integrand, int lIndex, const int nelx = (nx-1)/(p1-1); const int nely = (ny-1)/(p2-1); - Vector N(p1*p2); - Matrix dNdu, dNdX, Xnod, Jac; + FiniteElement fe(p1*p2); + Matrix dNdu, Xnod, Jac; Vec4 X; Vec3 normal; double xi[2]; @@ -326,21 +339,22 @@ bool ASMs2DLag::integrate (Integrand& integrand, int lIndex, // Compute the basis functions and their derivatives, using // tensor product of one-dimensional Lagrange polynomials - if (!Lagrange::computeBasis(N,dNdu,p1,xi[0],p2,xi[1])) + if (!Lagrange::computeBasis(fe.N,dNdu,p1,xi[0],p2,xi[1])) return false; // Compute basis function derivatives and the edge normal - double dS = utl::Jacobian(Jac,normal,dNdX,Xnod,dNdu,t1,t2); - if (dS == 0.0) continue; // skip singular points + fe.detJxW = utl::Jacobian(Jac,normal,fe.dNdX,Xnod,dNdu,t1,t2); + if (fe.detJxW == 0.0) continue; // skip singular points if (edgeDir < 0) normal *= -1.0; // Cartesian coordinates of current integration point - X = Xnod * N; + X = Xnod * fe.N; X.t = time.t; // Evaluate the integrand and accumulate element contributions - if (!integrand.evalBou(elmInt,time,dS*wg[i],N,dNdX,X,normal)) + fe.detJxW *= wg[i]; + if (!integrand.evalBou(elmInt,fe,time,X,normal)) return false; } diff --git a/src/ASM/ASMs2DSpec.C b/src/ASM/ASMs2DSpec.C index 6850674c..3b7c5556 100644 --- a/src/ASM/ASMs2DSpec.C +++ b/src/ASM/ASMs2DSpec.C @@ -15,6 +15,7 @@ #include "ASMs2DSpec.h" #include "TimeDomain.h" +#include "FiniteElement.h" #include "GlobalIntegral.h" #include "IntegrandBase.h" #include "CoordinateMapping.h" @@ -95,9 +96,8 @@ bool ASMs2DSpec::integrate (Integrand& integrand, if (!Legendre::basisDerivatives(p1,D1)) return false; if (!Legendre::basisDerivatives(p2,D2)) return false; - Vector N (p1*p2); - Matrix dNdu(p1*p2,2); - Matrix dNdX, Xnod, Jac; + FiniteElement fe(p1*p2); + Matrix dNdu(p1*p2,2), Xnod, Jac; Vec4 X; @@ -125,16 +125,13 @@ bool ASMs2DSpec::integrate (Integrand& integrand, for (int j = 1; j <= p2; j++) for (int i = 1; i <= p1; i++, count++) { - // Weight of current integration point - double weight = wg1(i)*wg2(j); - // Evaluate the basis functions and gradients using // tensor product of one-dimensional Lagrange polynomials - evalBasis(i,j,p1,p2,D1,D2,N,dNdu); + evalBasis(i,j,p1,p2,D1,D2,fe.N,dNdu); // Compute Jacobian inverse of coordinate mapping and derivatives - double dA = utl::Jacobian(Jac,dNdX,Xnod,dNdu); - if (dA == 0.0) continue; // skip singular points + fe.detJxW = utl::Jacobian(Jac,fe.dNdX,Xnod,dNdu); + if (fe.detJxW == 0.0) continue; // skip singular points // Cartesian coordinates of current integration point X.x = Xnod(1,count); @@ -142,7 +139,8 @@ bool ASMs2DSpec::integrate (Integrand& integrand, X.t = time.t; // Evaluate the integrand and accumulate element contributions - if (!integrand.evalInt(elmInt,time,dA*weight,N,dNdX,X)) + fe.detJxW *= wg1(i)*wg2(j); + if (!integrand.evalInt(elmInt,fe,time,X)) return false; } @@ -191,9 +189,8 @@ bool ASMs2DSpec::integrate (Integrand& integrand, int lIndex, int nen = p[0]*p[1]; - Vector N(nen); - Matrix dNdu(nen,2); - Matrix dNdX, Xnod, Jac; + FiniteElement fe(nen); + Matrix dNdu(nen,2), Xnod, Jac; Vec4 X; Vec3 normal; int xi[2]; @@ -233,29 +230,27 @@ bool ASMs2DSpec::integrate (Integrand& integrand, int lIndex, for (int i = 0; i < p[t2-1]; i++) { - // Weight of current integration point - double weight = wg[t2-1][i]; - // "Coordinates" along the edge xi[t1-1] = edgeDir < 0 ? 1 : p[t1-1]; xi[t2-1] = i+1; // Evaluate the basis functions and gradients using // tensor product of one-dimensional Lagrange polynomials - evalBasis(xi[0],xi[1],p[0],p[1],D1,D2,N,dNdu); + evalBasis(xi[0],xi[1],p[0],p[1],D1,D2,fe.N,dNdu); // Compute basis function derivatives and the edge normal - double dS = utl::Jacobian(Jac,normal,dNdX,Xnod,dNdu,t1,t2); - if (dS == 0.0) continue; // skip singular points + fe.detJxW = utl::Jacobian(Jac,normal,fe.dNdX,Xnod,dNdu,t1,t2); + if (fe.detJxW == 0.0) continue; // skip singular points if (edgeDir < 0) normal *= -1.0; // Cartesian coordinates of current integration point - X = Xnod * N; + X = Xnod * fe.N; X.t = time.t; // Evaluate the integrand and accumulate element contributions - if (!integrand.evalBou(elmInt,time,dS*weight,N,dNdX,X,normal)) + fe.detJxW *= wg[t2-1][i]; + if (!integrand.evalBou(elmInt,fe,time,X,normal)) return false; } diff --git a/src/ASM/ASMs2Dmx.C b/src/ASM/ASMs2Dmx.C index 77f8aa91..0642f75c 100644 --- a/src/ASM/ASMs2Dmx.C +++ b/src/ASM/ASMs2Dmx.C @@ -15,6 +15,7 @@ #include "ASMs2Dmx.h" #include "TimeDomain.h" +#include "FiniteElement.h" #include "GlobalIntegral.h" #include "IntegrandBase.h" #include "CoordinateMapping.h" @@ -400,11 +401,10 @@ bool ASMs2Dmx::integrate (Integrand& integrand, const int n2 = surf->numCoefs_v(); const int nel1 = n1 - p1 + 1; - Vector N1(basis1->order_u()*basis1->order_v()); - Vector N2(basis2->order_u()*basis2->order_v()); - Matrix dN1du, dN1dX, dN2du, dN2dX, Xnod, Jac; + MxFiniteElement fe(basis1->order_u()*basis1->order_v(), + basis2->order_u()*basis2->order_v()); + Matrix dN1du, dN2du, Xnod, Jac; Vec4 X; - double detJ; // === Assembly loop over all elements in the patch ========================== @@ -423,7 +423,7 @@ bool ASMs2Dmx::integrate (Integrand& integrand, if (!this->getElementCoordinates(Xnod,iel)) return false; // Initialize element quantities - IntVec::iterator f2start = MNPC[iel-1].begin() + N1.size(); + IntVec::iterator f2start = MNPC[iel-1].begin() + fe.N1.size(); if (!integrand.initElement(IntVec(MNPC[iel-1].begin(),f2start), IntVec(f2start,MNPC[iel-1].end()),nb1)) return false; @@ -441,33 +441,35 @@ bool ASMs2Dmx::integrate (Integrand& integrand, for (int j = 0; j < nGauss; j++, ip += nGauss*(nel1-1)) for (int i = 0; i < nGauss; i++, ip++) { - // Weight of current integration point - double weight = 0.25*dA*wg[i]*wg[j]; + // Parameter values of current integration point + fe.u = gpar[0](i+1,i1-p1+1); + fe.v = gpar[1](j+1,i2-p2+1); // Fetch basis function derivatives at current integration point - extractBasis(spline1[ip],N1,dN1du); - extractBasis(spline2[ip],N2,dN2du); + extractBasis(spline1[ip],fe.N1,dN1du); + extractBasis(spline2[ip],fe.N2,dN2du); // Compute Jacobian inverse of the coordinate mapping and // basis function derivatives w.r.t. Cartesian coordinates if (geoUsesBasis1) { - detJ = utl::Jacobian(Jac,dN1dX,Xnod,dN1du); - dN2dX.multiply(dN2du,Jac); // dN2dX = dN2du * J^-1 + fe.detJxW = utl::Jacobian(Jac,fe.dN1dX,Xnod,dN1du); + fe.dN2dX.multiply(dN2du,Jac); // dN2dX = dN2du * J^-1 } else { - detJ = utl::Jacobian(Jac,dN2dX,Xnod,dN2du); - dN1dX.multiply(dN1du,Jac); // dN1dX = dN1du * J^-1 + fe.detJxW = utl::Jacobian(Jac,fe.dN2dX,Xnod,dN2du); + fe.dN1dX.multiply(dN1du,Jac); // dN1dX = dN1du * J^-1 } - if (detJ == 0.0) continue; // skip singular points + if (fe.detJxW == 0.0) continue; // skip singular points // Cartesian coordinates of current integration point - X = Xnod * (geoUsesBasis1 ? N1 : N2); + X = Xnod * (geoUsesBasis1 ? fe.N1 : fe.N2); X.t = time.t; // Evaluate the integrand and accumulate element contributions - if (!integrand.evalInt(elmInt,time,detJ*weight,N1,N2,dN1dX,dN2dX,X)) + fe.detJxW *= 0.25*dA*wg[i]*wg[j]; + if (!integrand.evalIntMx(elmInt,fe,time,X)) return false; } @@ -540,12 +542,11 @@ bool ASMs2Dmx::integrate (Integrand& integrand, int lIndex, const int n1 = surf->numCoefs_u(); const int n2 = surf->numCoefs_v(); - Vector N1(basis1->order_u()*basis1->order_v()); - Vector N2(basis2->order_u()*basis2->order_v()); - Matrix dN1du, dN1dX, dN2du, dN2dX, Xnod, Jac; + MxFiniteElement fe(basis1->order_u()*basis1->order_v(), + basis2->order_u()*basis2->order_v()); + Matrix dN1du, dN2du, Xnod, Jac; Vec4 X; Vec3 normal; - double detJ; // === Assembly loop over all elements on the patch edge ===================== @@ -575,7 +576,7 @@ bool ASMs2Dmx::integrate (Integrand& integrand, int lIndex, if (!this->getElementCoordinates(Xnod,iel)) return false; // Initialize element quantities - IntVec::iterator f2start = MNPC[iel-1].begin() + N1.size(); + IntVec::iterator f2start = MNPC[iel-1].begin() + fe.N1.size(); if (!integrand.initElementBou(IntVec(MNPC[iel-1].begin(),f2start), IntVec(f2start,MNPC[iel-1].end()),nb1)) return false; @@ -592,36 +593,33 @@ bool ASMs2Dmx::integrate (Integrand& integrand, int lIndex, int ip = (abs(edgeDir) == 1 ? i2-p2 : i1-p1)*nGauss; for (int i = 0; i < nGauss; i++, ip++) { - // Weight of current integration point - double weight = 0.5*dS*wg[i]; - // Fetch basis function derivatives at current integration point - extractBasis(spline1[ip],N1,dN1du); - extractBasis(spline2[ip],N2,dN2du); + extractBasis(spline1[ip],fe.N1,dN1du); + extractBasis(spline2[ip],fe.N2,dN2du); // Compute Jacobian inverse of the coordinate mapping and // basis function derivatives w.r.t. Cartesian coordinates if (geoUsesBasis1) { - detJ = utl::Jacobian(Jac,normal,dN1dX,Xnod,dN1du,t1,t2); - dN2dX.multiply(dN2du,Jac); // dN2dX = dN2du * J^-1 + fe.detJxW = utl::Jacobian(Jac,normal,fe.dN1dX,Xnod,dN1du,t1,t2); + fe.dN2dX.multiply(dN2du,Jac); // dN2dX = dN2du * J^-1 } else { - detJ = utl::Jacobian(Jac,normal,dN2dX,Xnod,dN2du,t1,t2); - dN1dX.multiply(dN1du,Jac); // dN1dX = dN1du * J^-1 + fe.detJxW = utl::Jacobian(Jac,normal,fe.dN2dX,Xnod,dN2du,t1,t2); + fe.dN1dX.multiply(dN1du,Jac); // dN1dX = dN1du * J^-1 } - if (detJ == 0.0) continue; // skip singular points + if (fe.detJxW == 0.0) continue; // skip singular points if (edgeDir < 0) normal *= -1.0; // Cartesian coordinates of current integration point - X = Xnod * (geoUsesBasis1 ? N1 : N2); + X = Xnod * (geoUsesBasis1 ? fe.N1 : fe.N2); X.t = time.t; // Evaluate the integrand and accumulate element contributions - if (!integrand.evalBou(elmInt,time,detJ*weight, - N1,N2,dN1dX,dN2dX,X,normal)) + fe.detJxW *= 0.5*dS*wg[i]; + if (!integrand.evalBouMx(elmInt,fe,time,X,normal)) return false; } diff --git a/src/ASM/ASMs2DmxLag.C b/src/ASM/ASMs2DmxLag.C index c2c569b8..262ec07d 100644 --- a/src/ASM/ASMs2DmxLag.C +++ b/src/ASM/ASMs2DmxLag.C @@ -16,6 +16,7 @@ #include "ASMs2DmxLag.h" #include "Lagrange.h" #include "TimeDomain.h" +#include "FiniteElement.h" #include "GlobalIntegral.h" #include "IntegrandBase.h" #include "CoordinateMapping.h" @@ -175,6 +176,15 @@ bool ASMs2DmxLag::integrate (Integrand& integrand, const double* wg = GaussQuadrature::getWeight(nGauss); if (!xg || !wg) return false; + // Get parametric coordinates of the elements + RealArray upar, vpar; + this->getGridParameters(upar,0,1); + this->getGridParameters(vpar,1,1); + + // Number of elements in each direction + const int nelx = upar.size() - 1; + const int nely = vpar.size() - 1; + // Order of basis in the two parametric directions (order = degree + 1) const int p1 = surf->order_u(); const int p2 = surf->order_v(); @@ -182,21 +192,22 @@ bool ASMs2DmxLag::integrate (Integrand& integrand, const int q1 = p1 - 1; const int q2 = p2 - 1; - Vector N1(p1*p2), N2(q1*q2); - Matrix dN1du, dN2du, dN1dX, dN2dX, Xnod, Jac; + MxFiniteElement fe(p1*p2,q1*q2); + Matrix dN1du, dN2du, Xnod, Jac; Vec4 X; // === Assembly loop over all elements in the patch ========================== - const int nel = this->getNoElms(); - for (int iel = 1; iel <= nel; iel++) + int iel = 1; + for (int i2 = 0; i2 < nely; i2++) + for (int i1 = 0; i1 < nelx; i1++, iel++) { // Set up control point coordinates for current element if (!this->getElementCoordinates(Xnod,iel)) return false; // Initialize element quantities - IntVec::iterator f2start = MNPC[iel-1].begin() + N1.size(); + IntVec::iterator f2start = MNPC[iel-1].begin() + fe.N1.size(); if (!integrand.initElement(IntVec(MNPC[iel-1].begin(),f2start), IntVec(f2start,MNPC[iel-1].end()),nb1)) return false; @@ -213,28 +224,30 @@ bool ASMs2DmxLag::integrate (Integrand& integrand, for (int j = 0; j < nGauss; j++) for (int i = 0; i < nGauss; i++) { - // Weight of current integration point - double weight = wg[i]*wg[j]; + // Parameter value of current integration point + fe.u = 0.5*(upar[i1]*(1.0-xg[i]) + upar[i1+1]*(1.0+xg[i])); + fe.v = 0.5*(vpar[i2]*(1.0-xg[j]) + vpar[i2+1]*(1.0+xg[j])); // Compute basis function derivatives at current integration point // using tensor product of one-dimensional Lagrange polynomials - if (!Lagrange::computeBasis(N1,dN1du,p1,xg[i],p2,xg[j])) + if (!Lagrange::computeBasis(fe.N1,dN1du,p1,xg[i],p2,xg[j])) return false; - if (!Lagrange::computeBasis(N2,dN2du,q1,xg[i],q2,xg[j])) + if (!Lagrange::computeBasis(fe.N2,dN2du,q1,xg[i],q2,xg[j])) return false; // Compute Jacobian inverse of coordinate mapping and derivatives - double dA = utl::Jacobian(Jac,dN1dX,Xnod,dN1du); - if (dA == 0.0) continue; // skip singular points + fe.detJxW = utl::Jacobian(Jac,fe.dN1dX,Xnod,dN1du); + if (fe.detJxW == 0.0) continue; // skip singular points - dN2dX.multiply(dN2du,Jac); // dN2dX = dN2du * J^-1 + fe.dN2dX.multiply(dN2du,Jac); // dN2dX = dN2du * J^-1 // Cartesian coordinates of current integration point - X = Xnod * N1; + X = Xnod * fe.N1; X.t = time.t; // Evaluate the integrand and accumulate element contributions - if (!integrand.evalInt(elmInt,time,dA*weight,N1,N2,dN1dX,dN2dX,X)) + fe.detJxW *= wg[i]*wg[j]; + if (!integrand.evalInt(elmInt,fe,time,X)) return false; } @@ -276,8 +289,8 @@ bool ASMs2DmxLag::integrate (Integrand& integrand, int lIndex, const int nelx = (nx2-1)/(q1-1); const int nely = (ny2-1)/(q2-1); - Vector N1(p1*p2), N2(q1*q2); - Matrix dN1du, dN2du, dN1dX, dN2dX, Xnod, Jac; + MxFiniteElement fe(p1*p2,q1*q2); + Matrix dN1du, dN2du, Xnod, Jac; Vec4 X; Vec3 normal; double xi[2]; @@ -304,7 +317,7 @@ bool ASMs2DmxLag::integrate (Integrand& integrand, int lIndex, if (!this->getElementCoordinates(Xnod,iel)) return false; // Initialize element quantities - IntVec::iterator f2start = MNPC[iel-1].begin() + N1.size(); + IntVec::iterator f2start = MNPC[iel-1].begin() + fe.N1.size(); if (!integrand.initElementBou(IntVec(MNPC[iel-1].begin(),f2start), IntVec(f2start,MNPC[iel-1].end()),nb1)) return false; @@ -326,26 +339,26 @@ bool ASMs2DmxLag::integrate (Integrand& integrand, int lIndex, // Compute the basis functions and their derivatives, using // tensor product of one-dimensional Lagrange polynomials - if (!Lagrange::computeBasis(N1,dN1du,p1,xi[0],p2,xi[1])) + if (!Lagrange::computeBasis(fe.N1,dN1du,p1,xi[0],p2,xi[1])) return false; - if (!Lagrange::computeBasis(N2,dN2du,q1,xi[0],q2,xi[1])) + if (!Lagrange::computeBasis(fe.N2,dN2du,q1,xi[0],q2,xi[1])) return false; // Compute basis function derivatives and the edge normal - double dS = utl::Jacobian(Jac,normal,dN1dX,Xnod,dN1du,t1,t2); - if (dS == 0.0) continue; // skip singular points + fe.detJxW = utl::Jacobian(Jac,normal,fe.dN1dX,Xnod,dN1du,t1,t2); + if (fe.detJxW == 0.0) continue; // skip singular points - dN2dX.multiply(dN2du,Jac); // dN2dX = dN2du * J^-1 + fe.dN2dX.multiply(dN2du,Jac); // dN2dX = dN2du * J^-1 if (edgeDir < 0) normal *= -1.0; // Cartesian coordinates of current integration point - X = Xnod * N1; + X = Xnod * fe.N1; X.t = time.t; // Evaluate the integrand and accumulate element contributions - if (!integrand.evalBou(elmInt,time,dS*wg[i], - N1,N2,dN1dX,dN2dX,X,normal)) + fe.detJxW *= wg[i]; + if (!integrand.evalBou(elmInt,fe,time,X,normal)) return false; } diff --git a/src/ASM/ASMs3D.C b/src/ASM/ASMs3D.C index d783276d..734588f7 100644 --- a/src/ASM/ASMs3D.C +++ b/src/ASM/ASMs3D.C @@ -17,6 +17,7 @@ #include "ASMs3D.h" #include "TimeDomain.h" +#include "FiniteElement.h" #include "GlobalIntegral.h" #include "IntegrandBase.h" #include "CoordinateMapping.h" @@ -1132,10 +1133,9 @@ bool ASMs3D::integrate (Integrand& integrand, const int nel1 = n1 - p1 + 1; const int nel2 = n2 - p2 + 1; - Vector N(p1*p2*p3), Navg; - Matrix dNdu, dNdX, Xnod, Jac; - Matrix3D d2Ndu2, d2NdX2, Hess; - double h = 0.0; + FiniteElement fe(p1*p2*p3); + Matrix dNdu, Xnod, Jac; + Matrix3D d2Ndu2, Hess; Vec4 X; @@ -1157,13 +1157,13 @@ bool ASMs3D::integrate (Integrand& integrand, // Compute characteristic element length, if needed if (integrand.getIntegrandType() == 2) - h = getElmSize(p1,p2,p3,Xnod); + fe.h = getElmSize(p1,p2,p3,Xnod); else if (integrand.getIntegrandType() == 3) { // --- Compute average value of basis functions over the element ----- - Navg.resize(p1*p2*p3,true); + fe.Navg.resize(p1*p2*p3,true); double vol = 0.0; int ip = (((i3-p3)*nGauss*nel2 + i2-p2)*nGauss*nel1 + i1-p1)*nGauss; for (int k = 0; k < nGauss; k++, ip += nGauss*(nel2-1)*nGauss*nel1) @@ -1171,20 +1171,20 @@ bool ASMs3D::integrate (Integrand& integrand, for (int i = 0; i < nGauss; i++, ip++) { // Fetch basis function derivatives at current integration point - extractBasis(spline[ip],N,dNdu); + extractBasis(spline[ip],fe.N,dNdu); // Compute Jacobian determinant of coordinate mapping // and multiply by weight of current integration point + double detJac = utl::Jacobian(Jac,fe.dNdX,Xnod,dNdu,false); double weight = 0.125*dV*wg[i]*wg[j]*wg[k]; - double detJxW = utl::Jacobian(Jac,dNdX,Xnod,dNdu,false)*weight; // Numerical quadrature - Navg += N*detJxW; - vol += detJxW; + fe.Navg.add(fe.N,detJac*weight); + vol += detJac*weight; } // Divide by element volume - Navg /= vol; + fe.Navg /= vol; } else if (integrand.getIntegrandType() == 4) @@ -1218,50 +1218,39 @@ bool ASMs3D::integrate (Integrand& integrand, for (int j = 0; j < nGauss; j++, ip += nGauss*(nel1-1)) for (int i = 0; i < nGauss; i++, ip++) { - // Weight of current integration point - double weight = 0.125*dV*wg[i]*wg[j]*wg[k]; + // Parameter values of current integration point + fe.u = gpar[0](i+1,i1-p1+1); + fe.v = gpar[1](j+1,i2-p2+1); + fe.w = gpar[2](k+1,i3-p3+1); // Fetch basis function derivatives at current integration point if (integrand.getIntegrandType() == 2) - extractBasis(spline2[ip],N,dNdu,d2Ndu2); + extractBasis(spline2[ip],fe.N,dNdu,d2Ndu2); else - extractBasis(spline[ip],N,dNdu); + extractBasis(spline[ip],fe.N,dNdu); // Compute Jacobian inverse of coordinate mapping and derivatives - double detJ = utl::Jacobian(Jac,dNdX,Xnod,dNdu); - if (detJ == 0.0) continue; // skip singular points + fe.detJxW = utl::Jacobian(Jac,fe.dNdX,Xnod,dNdu); + if (fe.detJxW == 0.0) continue; // skip singular points // Compute Hessian of coordinate mapping and 2nd order derivatives if (integrand.getIntegrandType() == 2) - if (!utl::Hessian(Hess,d2NdX2,Jac,Xnod,d2Ndu2,dNdu)) + if (!utl::Hessian(Hess,fe.d2NdX2,Jac,Xnod,d2Ndu2,dNdu)) return false; #if SP_DEBUG > 4 std::cout <<"\niel, ip = "<< iel <<" "<< ip - <<"\nN ="<< N <<"dNdX ="<< dNdX << std::endl; + <<"\nN ="<< fe.N <<"dNdX ="<< fe.dNdX << std::endl; #endif // Cartesian coordinates of current integration point - X = Xnod * N; + X = Xnod * fe.N; X.t = time.t; // Evaluate the integrand and accumulate element contributions - bool ok = false; - switch (integrand.getIntegrandType()) - { - case 2: - ok = integrand.evalInt(elmInt,time,detJ*weight, - N,dNdX,d2NdX2,X,h); - break; - case 3: - ok = integrand.evalInt(elmInt,time,detJ*weight, - N,dNdX,Navg,X); - break; - default: - ok = integrand.evalInt(elmInt,time,detJ*weight, - N,dNdX,X); - } - if (!ok) return false; + fe.detJxW *= 0.125*dV*wg[i]*wg[j]*wg[k]; + if (!integrand.evalInt(elmInt,fe,time,X)) + return false; } // Finalize the element quantities @@ -1344,8 +1333,8 @@ bool ASMs3D::integrate (Integrand& integrand, int lIndex, const int nel1 = n1 - p1 + 1; const int nel2 = n2 - p2 + 1; - Vector N(p1*p2*p3); - Matrix dNdu, dNdX, Xnod, Jac; + FiniteElement fe(p1*p2*p3); + Matrix dNdu, Xnod, Jac; Vec4 X; Vec3 normal; @@ -1405,24 +1394,22 @@ bool ASMs3D::integrate (Integrand& integrand, int lIndex, for (int j = 0; j < nGauss; j++, ip += nGauss*(nf1-1)) for (int i = 0; i < nGauss; i++, ip++) { - // Weight of current integration point - double weight = 0.25*dA*wg[i]*wg[j]; - // Fetch basis function derivatives at current integration point - extractBasis(spline[ip],N,dNdu); + extractBasis(spline[ip],fe.N,dNdu); // Compute basis function derivatives and the face normal - double dS = utl::Jacobian(Jac,normal,dNdX,Xnod,dNdu,t1,t2); - if (dS == 0.0) continue; // skip singular points + fe.detJxW = utl::Jacobian(Jac,normal,fe.dNdX,Xnod,dNdu,t1,t2); + if (fe.detJxW == 0.0) continue; // skip singular points if (faceDir < 0) normal *= -1.0; // Cartesian coordinates of current integration point - X = Xnod * N; + X = Xnod * fe.N; X.t = time.t; // Evaluate the integrand and accumulate element contributions - if (!integrand.evalBou(elmInt,time,dS*weight,N,dNdX,X,normal)) + fe.detJxW *= 0.25*dA*wg[i]*wg[j]; + if (!integrand.evalBou(elmInt,fe,time,X,normal)) return false; } @@ -1496,10 +1483,10 @@ bool ASMs3D::integrateEdge (Integrand& integrand, int lEdge, const int p2 = svol->order(1); const int p3 = svol->order(2); - Vector N(p1*p2*p3); - Matrix dNdu, dNdX, Xnod, Jac; + FiniteElement fe(p1*p2*p3); + Matrix dNdu, Xnod, Jac; Vec4 X; - Vec3 tangent; + Vec3 tang; // === Assembly loop over all elements on the patch edge ===================== @@ -1564,22 +1551,20 @@ bool ASMs3D::integrateEdge (Integrand& integrand, int lEdge, LocalIntegral* elmInt = 0; for (int i = 0; i < nGauss; i++, ip++) { - // Weight of current integration point - double weight = 0.5*dS*wg[i]; - // Fetch basis function derivatives at current integration point - extractBasis(spline[ip],N,dNdu); + extractBasis(spline[ip],fe.N,dNdu); - // Compute basis function derivatives and the edge tangent - double dT = utl::Jacobian(Jac,tangent,dNdX,Xnod,dNdu,1+(lEdge-1)/4); - if (dT == 0.0) continue; // skip singular points + // Compute basis function derivatives and the edge tang + fe.detJxW = utl::Jacobian(Jac,tang,fe.dNdX,Xnod,dNdu,1+(lEdge-1)/4); + if (fe.detJxW == 0.0) continue; // skip singular points // Cartesian coordinates of current integration point - X = Xnod * N; + X = Xnod * fe.N; X.t = time.t; // Evaluate the integrand and accumulate element contributions - if (!integrand.evalBou(elmInt,time,dT*weight,N,dNdX,X,tangent)) + fe.detJxW *= 0.5*dS*wg[i]; + if (!integrand.evalBou(elmInt,fe,time,X,tang)) return false; } @@ -1604,12 +1589,14 @@ bool ASMs3D::getGridParameters (RealArray& prm, int dir, int nSegPerSpan) const } DoubleIter uit = svol->basis(dir).begin(); - double ucurr, uprev = *(uit++); + double ucurr = 0.0, uprev = *(uit++); while (uit != svol->basis(dir).end()) { ucurr = *(uit++); if (ucurr > uprev) - for (int i = 0; i < nSegPerSpan; i++) + if (nSegPerSpan == 1) + prm.push_back(uprev); + else for (int i = 0; i < nSegPerSpan; i++) { double xg = (double)(2*i-nSegPerSpan)/(double)nSegPerSpan; prm.push_back(0.5*(ucurr*(1.0+xg) + uprev*(1.0-xg))); @@ -1617,7 +1604,8 @@ bool ASMs3D::getGridParameters (RealArray& prm, int dir, int nSegPerSpan) const uprev = ucurr; } - prm.push_back(svol->basis(dir).endparam()); + if (ucurr > prm.back()) + prm.push_back(ucurr); return true; } @@ -1745,7 +1733,7 @@ bool ASMs3D::evalSolution (Matrix& sField, const Integrand& integrand, if (npe) { // Compute parameter values of the result sampling points - DoubleVec gpar[3]; + RealArray gpar[3]; for (int dir = 0; dir < 3 && retVal; dir++) retVal = this->getGridParameters(gpar[dir],dir,npe[dir]-1); diff --git a/src/ASM/ASMs3DLag.C b/src/ASM/ASMs3DLag.C index 2d3f975a..eabb6487 100644 --- a/src/ASM/ASMs3DLag.C +++ b/src/ASM/ASMs3DLag.C @@ -16,6 +16,7 @@ #include "ASMs3DLag.h" #include "Lagrange.h" #include "TimeDomain.h" +#include "FiniteElement.h" #include "GlobalIntegral.h" #include "IntegrandBase.h" #include "CoordinateMapping.h" @@ -104,7 +105,7 @@ bool ASMs3DLag::generateFEMTopology() int i, j, k, a, b, c, iel = 0; for (k = 0; k < nelz; k++) for (j = 0; j < nely; j++) - for ( i = 0; i < nelx; i++, iel++) + for (i = 0; i < nelx; i++, iel++) { MLGE[iel] = ++gEl; MNPC[iel].resize(nen); @@ -191,21 +192,34 @@ bool ASMs3DLag::integrate (Integrand& integrand, const double* wg = GaussQuadrature::getWeight(nGauss); if (!xg || !wg) return false; + // Get parametric coordinates of the elements + RealArray upar, vpar, wpar; + this->getGridParameters(upar,0,1); + this->getGridParameters(vpar,1,1); + this->getGridParameters(wpar,2,1); + + // Number of elements in each direction + const int nelx = upar.size() - 1; + const int nely = vpar.size() - 1; + const int nelz = wpar.size() - 1; + // Order of basis in the three parametric directions (order = degree + 1) const int p1 = svol->order(0); const int p2 = svol->order(1); const int p3 = svol->order(2); - Vector N(p1*p2*p3); - Matrix dNdu, dNdX, Xnod, Jac; + FiniteElement fe(p1*p2*p3); + Matrix dNdu, Xnod, Jac; Vec4 X; // === Assembly loop over all elements in the patch ========================== - const int nel = this->getNoElms(); - for (int iel = 1; iel <= nel; iel++) - { + int iel = 1; + for (int i3 = 0; i3 < nelz; i3++) + for (int i2 = 0; i2 < nely; i2++) + for (int i1 = 0; i1 < nelx; i1++, iel++) + { // Set up nodal point coordinates for current element if (!this->getElementCoordinates(Xnod,iel)) return false; @@ -237,24 +251,27 @@ bool ASMs3DLag::integrate (Integrand& integrand, for (int j = 0; j < nGauss; j++) for (int i = 0; i < nGauss; i++) { - // Weight of current integration point - double weight = wg[i]*wg[j]*wg[k]; + // Parameter value of current integration point + fe.u = 0.5*(upar[i1]*(1.0-xg[i]) + upar[i1+1]*(1.0+xg[i])); + fe.v = 0.5*(vpar[i2]*(1.0-xg[j]) + vpar[i2+1]*(1.0+xg[j])); + fe.w = 0.5*(wpar[i3]*(1.0-xg[k]) + wpar[i3+1]*(1.0+xg[k])); // Compute basis function derivatives at current integration point // using tensor product of one-dimensional Lagrange polynomials - if (!Lagrange::computeBasis(N,dNdu,p1,xg[i],p2,xg[j],p3,xg[k])) + if (!Lagrange::computeBasis(fe.N,dNdu,p1,xg[i],p2,xg[j],p3,xg[k])) return false; // Compute Jacobian inverse of coordinate mapping and derivatives - double detJ = utl::Jacobian(Jac,dNdX,Xnod,dNdu); - if (detJ == 0.0) continue; // skip singular points + fe.detJxW = utl::Jacobian(Jac,fe.dNdX,Xnod,dNdu); + if (fe.detJxW == 0.0) continue; // skip singular points // Cartesian coordinates of current integration point - X = Xnod * N; + X = Xnod * fe.N; X.t = time.t; // Evaluate the integrand and accumulate element contributions - if (!integrand.evalInt(elmInt,time,detJ*weight,N,dNdX,X)) + fe.detJxW *= wg[i]*wg[j]*wg[k]; + if (!integrand.evalInt(elmInt,fe,time,X)) return false; } @@ -265,7 +282,7 @@ bool ASMs3DLag::integrate (Integrand& integrand, // Assembly of global system integral if (!glInt.assemble(elmInt,MLGE[iel-1])) return false; - } + } return true; } @@ -300,8 +317,8 @@ bool ASMs3DLag::integrate (Integrand& integrand, int lIndex, const int nely = (ny-1)/(p2-1); const int nelz = (nz-1)/(p3-1); - Vector N(p1*p2*p3); - Matrix dNdu, dNdX, Xnod, Jac; + FiniteElement fe(p1*p2*p3); + Matrix dNdu, Xnod, Jac; Vec4 X; Vec3 normal; double xi[3]; @@ -345,9 +362,6 @@ bool ASMs3DLag::integrate (Integrand& integrand, int lIndex, for (int j = 0; j < nGauss; j++) for (int i = 0; i < nGauss; i++) { - // Weight of current integration point - double weight = wg[i]*wg[j]; - // Gauss point coordinates on the face xi[t0-1] = faceDir < 0 ? -1.0 : 1.0; xi[t1-1] = xg[i]; @@ -355,21 +369,22 @@ bool ASMs3DLag::integrate (Integrand& integrand, int lIndex, // Compute the basis functions and their derivatives, using // tensor product of one-dimensional Lagrange polynomials - if (!Lagrange::computeBasis(N,dNdu,p1,xi[0],p2,xi[1],p3,xi[2])) + if (!Lagrange::computeBasis(fe.N,dNdu,p1,xi[0],p2,xi[1],p3,xi[2])) return false; // Compute basis function derivatives and the face normal - double dS = utl::Jacobian(Jac,normal,dNdX,Xnod,dNdu,t1,t2); - if (dS == 0.0) continue; // skip singular points + fe.detJxW = utl::Jacobian(Jac,normal,fe.dNdX,Xnod,dNdu,t1,t2); + if (fe.detJxW == 0.0) continue; // skip singular points if (faceDir < 0) normal *= -1.0; // Cartesian coordinates of current integration point - X = Xnod * N; + X = Xnod * fe.N; X.t = time.t; // Evaluate the integrand and accumulate element contributions - if (!integrand.evalBou(elmInt,time,dS*weight,N,dNdX,X,normal)) + fe.detJxW *= wg[i]*wg[j]; + if (!integrand.evalBou(elmInt,fe,time,X,normal)) return false; } @@ -406,8 +421,8 @@ bool ASMs3DLag::integrateEdge (Integrand& integrand, int lEdge, const int nely = (ny-1)/(p2-1); const int nelz = (nz-1)/(p3-1); - Vector N(p1*p2*p3); - Matrix dNdu, dNdX, Xnod, Jac; + FiniteElement fe(p1*p2*p3); + Matrix dNdu, Xnod, Jac; Vec4 X; Vec3 tangent; double xi[3]; @@ -472,19 +487,19 @@ bool ASMs3DLag::integrateEdge (Integrand& integrand, int lEdge, // Compute the basis functions and their derivatives, using // tensor product of one-dimensional Lagrange polynomials - if (!Lagrange::computeBasis(N,dNdu,p1,xi[0],p2,xi[1],p3,xi[2])) + if (!Lagrange::computeBasis(fe.N,dNdu,p1,xi[0],p2,xi[1],p3,xi[2])) return false; // Compute basis function derivatives and the edge tangent - double dS = utl::Jacobian(Jac,tangent,dNdX,Xnod,dNdu,1+lDir); - if (dS == 0.0) continue; // skip singular points + fe.detJxW = utl::Jacobian(Jac,tangent,fe.dNdX,Xnod,dNdu,1+lDir); + if (fe.detJxW == 0.0) continue; // skip singular points // Cartesian coordinates of current integration point - X = Xnod * N; + X = Xnod * fe.N; X.t = time.t; // Evaluate the integrand and accumulate element contributions - if (!integrand.evalBou(elmInt,time,dS*wg[i],N,dNdX,X,tangent)) + if (!integrand.evalBou(elmInt,fe,time,X,tangent)) return false; } diff --git a/src/ASM/ASMs3DSpec.C b/src/ASM/ASMs3DSpec.C index 6951af0a..992841f4 100644 --- a/src/ASM/ASMs3DSpec.C +++ b/src/ASM/ASMs3DSpec.C @@ -15,6 +15,7 @@ #include "ASMs3DSpec.h" #include "TimeDomain.h" +#include "FiniteElement.h" #include "GlobalIntegral.h" #include "IntegrandBase.h" #include "CoordinateMapping.h" @@ -100,9 +101,8 @@ bool ASMs3DSpec::integrate (Integrand& integrand, if (!Legendre::basisDerivatives(p2,D2)) return false; if (!Legendre::basisDerivatives(p3,D3)) return false; - Vector N(p1*p2*p3); - Matrix dNdu(p1*p2*p3,3); - Matrix dNdX, Xnod, Jac; + FiniteElement fe(p1*p2*p3); + Matrix dNdu(p1*p2*p3,3), Xnod, Jac; Vec4 X; @@ -131,16 +131,13 @@ bool ASMs3DSpec::integrate (Integrand& integrand, for (int j = 1; j <= p2; j++) for (int i = 1; i <= p1; i++, count++) { - // Weight of current integration point - double weight = wg1(i)*wg2(j)*wg3(k); - // Evaluate the basis functions and gradients using tensor product // of the one-dimensional Lagrange polynomials - evalBasis(i,j,k,p1,p2,p3,D1,D2,D3,N,dNdu); + evalBasis(i,j,k,p1,p2,p3,D1,D2,D3,fe.N,dNdu); // Compute Jacobian inverse of coordinate mapping and derivatives - double detJ = utl::Jacobian(Jac,dNdX,Xnod,dNdu); - if (detJ == 0.0) continue; // skip singular points + fe.detJxW = utl::Jacobian(Jac,fe.dNdX,Xnod,dNdu); + if (fe.detJxW == 0.0) continue; // skip singular points // Cartesian coordinates of current integration point X.x = Xnod(1,count); @@ -149,7 +146,8 @@ bool ASMs3DSpec::integrate (Integrand& integrand, X.t = time.t; // Evaluate the integrand and accumulate element contributions - if (!integrand.evalInt(elmInt,time,detJ*weight,N,dNdX,X)) + fe.detJxW *= wg1(i)*wg2(j)*wg3(k); + if (!integrand.evalInt(elmInt,fe,time,X)) return false; } @@ -203,9 +201,8 @@ bool ASMs3DSpec::integrate (Integrand& integrand, int lIndex, int nen = p[0]*p[1]*p[2]; - Vector N(nen); - Matrix dNdu(nen,3); - Matrix dNdX, Xnod, Jac; + FiniteElement fe(nen); + Matrix dNdu(nen,3), Xnod, Jac; Vec4 X; Vec3 normal; int xi[3]; @@ -249,9 +246,6 @@ bool ASMs3DSpec::integrate (Integrand& integrand, int lIndex, for (int j = 0; j < p[t2-1]; j++) for (int i = 0; i < p[t1-1]; i++) { - // Weight of current integration point - double weight = wg[t1-1][i]*wg[t2-1][j]; - // "Coordinates" on the face xi[t0-1] = faceDir < 0 ? 1 : p[t0-1]; xi[t1-1] = i+1; @@ -259,20 +253,21 @@ bool ASMs3DSpec::integrate (Integrand& integrand, int lIndex, // Compute the basis functions and their derivatives, using // tensor product of one-dimensional Lagrange polynomials - evalBasis(xi[0],xi[1],xi[2],p[0],p[1],p[2],D1,D2,D3,N,dNdu); + evalBasis(xi[0],xi[1],xi[2],p[0],p[1],p[2],D1,D2,D3,fe.N,dNdu); // Compute basis function derivatives and the face normal - double dS = utl::Jacobian(Jac,normal,dNdX,Xnod,dNdu,t1,t2); - if (dS == 0.0) continue; // skip singular points + fe.detJxW = utl::Jacobian(Jac,normal,fe.dNdX,Xnod,dNdu,t1,t2); + if (fe.detJxW == 0.0) continue; // skip singular points if (faceDir < 0) normal *= -1.0; // Cartesian coordinates of current integration point - X = Xnod * N; + X = Xnod * fe.N; X.t = time.t; // Evaluate the integrand and accumulate element contributions - if (!integrand.evalBou(elmInt,time,dS*weight,N,dNdX,X,normal)) + fe.detJxW *= wg[t1-1][i]*wg[t2-1][j]; + if (!integrand.evalBou(elmInt,fe,time,X,normal)) return false; } @@ -321,9 +316,8 @@ bool ASMs3DSpec::integrateEdge (Integrand& integrand, int lEdge, const int nen = p1*p2*p3; - Vector N(nen); - Matrix dNdu(nen,3); - Matrix dNdX, Xnod, Jac; + FiniteElement fe(nen); + Matrix dNdu(nen,3), Xnod, Jac; Vec4 X; Vec3 tangent; int xi[3]; @@ -388,18 +382,19 @@ bool ASMs3DSpec::integrateEdge (Integrand& integrand, int lEdge, // Compute the basis functions and their derivatives, using // tensor product of one-dimensional Lagrange polynomials - evalBasis(xi[0],xi[1],xi[2],p1,p2,p3,D1,D2,D3,N,dNdu); + evalBasis(xi[0],xi[1],xi[2],p1,p2,p3,D1,D2,D3,fe.N,dNdu); // Compute basis function derivatives and the edge tangent - double dS = utl::Jacobian(Jac,tangent,dNdX,Xnod,dNdu,1+lDir); - if (dS == 0.0) continue; // skip singular points + fe.detJxW = utl::Jacobian(Jac,tangent,fe.dNdX,Xnod,dNdu,1+lDir); + if (fe.detJxW == 0.0) continue; // skip singular points // Cartesian coordinates of current integration point - X = Xnod * N; + X = Xnod * fe.N; X.t = time.t; // Evaluate the integrand and accumulate element contributions - if (!integrand.evalBou(elmInt,time,dS*wg[lDir][i],N,dNdX,X,tangent)) + fe.detJxW *= wg[lDir][i]; + if (!integrand.evalBou(elmInt,fe,time,X,tangent)) return false; } @@ -434,9 +429,9 @@ bool ASMs3DSpec::evalSolution (Matrix& sField, const Integrand& integrand, size_t nPoints = this->getNoNodes(); IntVec check(nPoints,0); - Vector N(p1*p2*p3), solPt; - std::vector globSolPt(nPoints); - Matrix dNdu(p1*p2*p3,3), dNdX, Xnod, Jac; + Vector N(p1*p2*p3), solPt; + Vectors globSolPt(nPoints); + Matrix dNdu(p1*p2*p3,3), dNdX, Xnod, Jac; // Evaluate the secondary solution field at each point const int nel = this->getNoElms(); diff --git a/src/ASM/ASMs3Dmx.C b/src/ASM/ASMs3Dmx.C index cd487851..bb60a8e3 100644 --- a/src/ASM/ASMs3Dmx.C +++ b/src/ASM/ASMs3Dmx.C @@ -15,6 +15,7 @@ #include "ASMs3Dmx.h" #include "TimeDomain.h" +#include "FiniteElement.h" #include "GlobalIntegral.h" #include "IntegrandBase.h" #include "CoordinateMapping.h" @@ -429,11 +430,10 @@ bool ASMs3Dmx::integrate (Integrand& integrand, const int nel1 = n1 - p1 + 1; const int nel2 = n2 - p2 + 1; - Vector N1(basis1->order(0)*basis1->order(1)*basis1->order(2)); - Vector N2(basis2->order(0)*basis2->order(1)*basis2->order(2)); - Matrix dN1du, dN1dX, dN2du, dN2dX, Xnod, Jac; + MxFiniteElement fe(basis1->order(0)*basis1->order(1)*basis1->order(2), + basis2->order(0)*basis2->order(1)*basis2->order(2)); + Matrix dN1du, dN2du, Xnod, Jac; Vec4 X; - double detJ; // === Assembly loop over all elements in the patch ========================== @@ -453,7 +453,7 @@ bool ASMs3Dmx::integrate (Integrand& integrand, if (!this->getElementCoordinates(Xnod,iel)) return false; // Initialize element quantities - IntVec::iterator f2start = MNPC[iel-1].begin() + N1.size(); + IntVec::iterator f2start = MNPC[iel-1].begin() + fe.N1.size(); if (!integrand.initElement(IntVec(MNPC[iel-1].begin(),f2start), IntVec(f2start,MNPC[iel-1].end()),nb1)) return false; @@ -472,34 +472,36 @@ bool ASMs3Dmx::integrate (Integrand& integrand, for (int j = 0; j < nGauss; j++, ip += nGauss*(nel1-1)) for (int i = 0; i < nGauss; i++, ip++) { - // Weight of current integration point - double weight = 0.125*dV*wg[i]*wg[j]*wg[k]; + // Parameter values of current integration point + fe.u = gpar[0](i+1,i1-p1+1); + fe.v = gpar[1](j+1,i2-p2+1); + fe.w = gpar[2](k+1,i3-p3+1); // Fetch basis function derivatives at current integration point - extractBasis(spline1[ip],N1,dN1du); - extractBasis(spline2[ip],N2,dN2du); + extractBasis(spline1[ip],fe.N1,dN1du); + extractBasis(spline2[ip],fe.N2,dN2du); // Compute Jacobian inverse of the coordinate mapping and // basis function derivatives w.r.t. Cartesian coordinates if (geoUsesBasis1) { - detJ = utl::Jacobian(Jac,dN1dX,Xnod,dN1du); - dN2dX.multiply(dN2du,Jac); // dN2dX = dN2du * J^-1 + fe.detJxW = utl::Jacobian(Jac,fe.dN1dX,Xnod,dN1du); + fe.dN2dX.multiply(dN2du,Jac); // dN2dX = dN2du * J^-1 } else { - detJ = utl::Jacobian(Jac,dN2dX,Xnod,dN2du); - dN1dX.multiply(dN1du,Jac); // dN1dX = dN1du * J^-1 + fe.detJxW = utl::Jacobian(Jac,fe.dN2dX,Xnod,dN2du); + fe.dN1dX.multiply(dN1du,Jac); // dN1dX = dN1du * J^-1 } - if (detJ == 0.0) continue; // skip singular points + if (fe.detJxW == 0.0) continue; // skip singular points // Cartesian coordinates of current integration point - X = Xnod * (geoUsesBasis1 ? N1 : N2); + X = Xnod * (geoUsesBasis1 ? fe.N1 : fe.N2); X.t = time.t; // Evaluate the integrand and accumulate element contributions - if (!integrand.evalInt(elmInt,time,detJ*weight, - N1,N2,dN1dX,dN2dX,X)) + fe.detJxW *= 0.125*dV*wg[i]*wg[j]*wg[k]; + if (!integrand.evalIntMx(elmInt,fe,time,X)) return false; } @@ -533,7 +535,7 @@ bool ASMs3Dmx::integrate (Integrand& integrand, int lIndex, const int t1 = 1 + abs(faceDir)%3; // first tangent direction const int t2 = 1 + t1%3; // second tangent direction - // Compute parameter values of the Gauss points along the whole patch edge + // Compute parameter values of the Gauss points over the whole patch face Matrix gpar[3]; for (int d = 0; d < 3; d++) if (-1-d == faceDir) @@ -578,12 +580,12 @@ bool ASMs3Dmx::integrate (Integrand& integrand, int lIndex, const int nel1 = n1 - p1 + 1; const int nel2 = n2 - p2 + 1; - Vector N1(basis1->order(0)*basis1->order(1)*basis1->order(2)); - Vector N2(basis2->order(0)*basis2->order(1)*basis2->order(2)); - Matrix dN1du, dN1dX, dN2du, dN2dX, Xnod, Jac; + MxFiniteElement fe(basis1->order(0)*basis1->order(1)*basis1->order(2), + basis2->order(0)*basis2->order(1)*basis2->order(2)); + + Matrix dN1du, dN2du, Xnod, Jac; Vec4 X; Vec3 normal; - double dS; // === Assembly loop over all elements on the patch face ===================== @@ -616,7 +618,7 @@ bool ASMs3Dmx::integrate (Integrand& integrand, int lIndex, if (!this->getElementCoordinates(Xnod,iel)) return false; // Initialize element quantities - IntVec::iterator f2start = MNPC[iel-1].begin() + N1.size(); + IntVec::iterator f2start = MNPC[iel-1].begin() + fe.N1.size(); if (!integrand.initElementBou(IntVec(MNPC[iel-1].begin(),f2start), IntVec(f2start,MNPC[iel-1].end()),nb1)) return false; @@ -644,36 +646,33 @@ bool ASMs3Dmx::integrate (Integrand& integrand, int lIndex, for (int j = 0; j < nGauss; j++, ip += nGauss*(nf1-1)) for (int i = 0; i < nGauss; i++, ip++) { - // Weight of current integration point - double weight = 0.25*dA*wg[i]*wg[j]; - // Fetch basis function derivatives at current integration point - extractBasis(spline1[ip],N1,dN1du); - extractBasis(spline2[ip],N2,dN2du); + extractBasis(spline1[ip],fe.N1,dN1du); + extractBasis(spline2[ip],fe.N2,dN2du); // Compute Jacobian inverse of the coordinate mapping and // basis function derivatives w.r.t. Cartesian coordinates if (geoUsesBasis1) { - dS = utl::Jacobian(Jac,normal,dN1dX,Xnod,dN1du,t1,t2); - dN2dX.multiply(dN2du,Jac); // dN2dX = dN2du * J^-1 + fe.detJxW = utl::Jacobian(Jac,normal,fe.dN1dX,Xnod,dN1du,t1,t2); + fe.dN2dX.multiply(dN2du,Jac); // dN2dX = dN2du * J^-1 } else { - dS = utl::Jacobian(Jac,normal,dN2dX,Xnod,dN2du,t1,t2); - dN1dX.multiply(dN1du,Jac); // dN1dX = dN1du * J^-1 + fe.detJxW = utl::Jacobian(Jac,normal,fe.dN2dX,Xnod,dN2du,t1,t2); + fe.dN1dX.multiply(dN1du,Jac); // dN1dX = dN1du * J^-1 } - if (dS == 0.0) continue; // skip singular points + if (fe.detJxW == 0.0) continue; // skip singular points if (faceDir < 0) normal *= -1.0; // Cartesian coordinates of current integration point - X = Xnod * (geoUsesBasis1 ? N1 : N2); + X = Xnod * (geoUsesBasis1 ? fe.N1 : fe.N2); X.t = time.t; // Evaluate the integrand and accumulate element contributions - if (!integrand.evalBou(elmInt,time,dS*weight, - N1,N2,dN1dX,dN2dX,X,normal)) + fe.detJxW *= 0.25*dA*wg[i]*wg[j]; + if (!integrand.evalBouMx(elmInt,fe,time,X,normal)) return false; } diff --git a/src/ASM/ASMs3DmxLag.C b/src/ASM/ASMs3DmxLag.C index 0f2e7695..2bf074c7 100644 --- a/src/ASM/ASMs3DmxLag.C +++ b/src/ASM/ASMs3DmxLag.C @@ -16,6 +16,7 @@ #include "ASMs3DmxLag.h" #include "Lagrange.h" #include "TimeDomain.h" +#include "FiniteElement.h" #include "GlobalIntegral.h" #include "IntegrandBase.h" #include "CoordinateMapping.h" @@ -186,9 +187,20 @@ bool ASMs3DmxLag::integrate (Integrand& integrand, if (!svol) return true; // silently ignore empty patches // Get Gaussian quadrature points and weights - const double* xg = GaussQuadrature::getCoord(nGauss); - const double* wg = GaussQuadrature::getWeight(nGauss); - if (!xg || !wg) return false; + const double* x = GaussQuadrature::getCoord(nGauss); + const double* w = GaussQuadrature::getWeight(nGauss); + if (!x || !w) return false; + + // Get parametric coordinates of the elements + RealArray upar, vpar, wpar; + this->getGridParameters(upar,0,1); + this->getGridParameters(vpar,1,1); + this->getGridParameters(wpar,2,1); + + // Number of elements in each direction + const int nelx = upar.size() - 1; + const int nely = vpar.size() - 1; + const int nelz = wpar.size() - 1; // Order of basis in the two parametric directions (order = degree + 1) const int p1 = svol->order(0); @@ -199,21 +211,23 @@ bool ASMs3DmxLag::integrate (Integrand& integrand, const int q2 = p2 - 1; const int q3 = p3 - 1; - Vector N1(p1*p2*p3), N2(q1*q2*q3); - Matrix dN1du, dN2du, dN1dX, dN2dX, Xnod, Jac; + MxFiniteElement fe(p1*p2*p3,q1*q2*q3); + Matrix dN1du, dN2du, Xnod, Jac; Vec4 X; // === Assembly loop over all elements in the patch ========================== - const int nel = this->getNoElms(); - for (int iel = 1; iel <= nel; iel++) - { + int iel = 1; + for (int i3 = 0; i3 < nelz; i3++) + for (int i2 = 0; i2 < nely; i2++) + for (int i1 = 0; i1 < nelx; i1++, iel++) + { // Set up control point coordinates for current element if (!this->getElementCoordinates(Xnod,iel)) return false; // Initialize element quantities - IntVec::iterator f2start = MNPC[iel-1].begin() + N1.size(); + IntVec::iterator f2start = MNPC[iel-1].begin() + fe.N1.size(); if (!integrand.initElement(IntVec(MNPC[iel-1].begin(),f2start), IntVec(f2start,MNPC[iel-1].end()),nb1)) return false; @@ -231,35 +245,38 @@ bool ASMs3DmxLag::integrate (Integrand& integrand, for (int j = 0; j < nGauss; j++) for (int i = 0; i < nGauss; i++) { - // Weight of current integration point - double weight = wg[i]*wg[j]*wg[k]; + // Parameter value of current integration point + fe.u = 0.5*(upar[i1]*(1.0-x[i]) + upar[i1+1]*(1.0+x[i])); + fe.v = 0.5*(vpar[i2]*(1.0-x[j]) + vpar[i2+1]*(1.0+x[j])); + fe.w = 0.5*(wpar[i3]*(1.0-x[k]) + wpar[i3+1]*(1.0+x[k])); // Compute basis function derivatives at current integration point // using tensor product of one-dimensional Lagrange polynomials - if (!Lagrange::computeBasis(N1,dN1du,p1,xg[i],p2,xg[j],p3,xg[k])) + if (!Lagrange::computeBasis(fe.N1,dN1du,p1,x[i],p2,x[j],p3,x[k])) return false; - if (!Lagrange::computeBasis(N2,dN2du,q1,xg[i],q2,xg[j],q3,xg[k])) + if (!Lagrange::computeBasis(fe.N2,dN2du,q1,x[i],q2,x[j],q3,x[k])) return false; // Compute Jacobian inverse of coordinate mapping and derivatives - double dA = utl::Jacobian(Jac,dN1dX,Xnod,dN1du); - if (dA == 0.0) continue; // skip singular points + fe.detJxW = utl::Jacobian(Jac,fe.dN1dX,Xnod,dN1du); + if (fe.detJxW == 0.0) continue; // skip singular points - dN2dX.multiply(dN2du,Jac); // dN2dX = dN2du * J^-1 + fe.dN2dX.multiply(dN2du,Jac); // dN2dX = dN2du * J^-1 // Cartesian coordinates of current integration point - X = Xnod * N1; + X = Xnod * fe.N1; X.t = time.t; // Evaluate the integrand and accumulate element contributions - if (!integrand.evalInt(elmInt,time,dA*weight,N1,N2,dN1dX,dN2dX,X)) + fe.detJxW *= w[i]*w[j]*w[k]; + if (!integrand.evalIntMx(elmInt,fe,time,X)) return false; } // Assembly of global system integral if (!glInt.assemble(elmInt,MLGE[iel-1])) return false; - } + } return true; } @@ -298,8 +315,8 @@ bool ASMs3DmxLag::integrate (Integrand& integrand, int lIndex, const int nely = (ny2-1)/(q2-1); const int nelz = (nz2-1)/(q3-1); - Vector N1(p1*p2*p3), N2(q1*q2*q3); - Matrix dN1du, dN2du, dN1dX, dN2dX, Xnod, Jac; + MxFiniteElement fe(p1*p2*p3,q1*q2*q3); + Matrix dN1du, dN2du, Xnod, Jac; Vec4 X; Vec3 normal; double xi[3]; @@ -329,7 +346,7 @@ bool ASMs3DmxLag::integrate (Integrand& integrand, int lIndex, if (!this->getElementCoordinates(Xnod,iel)) return false; // Initialize element quantities - IntVec::iterator f2start = MNPC[iel-1].begin() + N1.size(); + IntVec::iterator f2start = MNPC[iel-1].begin() + fe.N1.size(); if (!integrand.initElementBou(IntVec(MNPC[iel-1].begin(),f2start), IntVec(f2start,MNPC[iel-1].end()),nb1)) return false; @@ -346,9 +363,6 @@ bool ASMs3DmxLag::integrate (Integrand& integrand, int lIndex, for (int j = 0; j < nGauss; j++) for (int i = 0; i < nGauss; i++) { - // Weight of current integration point - double weight = wg[i]*wg[j]; - // Gauss point coordinates on the face xi[t0-1] = faceDir < 0 ? -1.0 : 1.0; xi[t1-1] = xg[i]; @@ -356,26 +370,26 @@ bool ASMs3DmxLag::integrate (Integrand& integrand, int lIndex, // Compute the basis functions and their derivatives, using // tensor product of one-dimensional Lagrange polynomials - if (!Lagrange::computeBasis(N1,dN1du,p1,xi[0],p2,xi[1],p3,xi[2])) + if (!Lagrange::computeBasis(fe.N1,dN1du,p1,xi[0],p2,xi[1],p3,xi[2])) return false; - if (!Lagrange::computeBasis(N2,dN2du,q1,xi[0],q2,xi[1],q3,xi[2])) + if (!Lagrange::computeBasis(fe.N2,dN2du,q1,xi[0],q2,xi[1],q3,xi[2])) return false; // Compute basis function derivatives and the edge normal - double dS = utl::Jacobian(Jac,normal,dN1dX,Xnod,dN1du,t1,t2); - if (dS == 0.0) continue; // skip singular points + fe.detJxW = utl::Jacobian(Jac,normal,fe.dN1dX,Xnod,dN1du,t1,t2); + if (fe.detJxW == 0.0) continue; // skip singular points - dN2dX.multiply(dN2du,Jac); // dN2dX = dN2du * J^-1 + fe.dN2dX.multiply(dN2du,Jac); // dN2dX = dN2du * J^-1 if (faceDir < 0) normal *= -1.0; // Cartesian coordinates of current integration point - X = Xnod * N1; + X = Xnod * fe.N1; X.t = time.t; // Evaluate the integrand and accumulate element contributions - if (!integrand.evalBou(elmInt,time,dS*weight, - N1,N2,dN1dX,dN2dX,X,normal)) + fe.detJxW *= wg[i]*wg[j]; + if (!integrand.evalBouMx(elmInt,fe,time,X,normal)) return false; } diff --git a/src/ASM/FiniteElement.h b/src/ASM/FiniteElement.h new file mode 100644 index 00000000..a1320a75 --- /dev/null +++ b/src/ASM/FiniteElement.h @@ -0,0 +1,59 @@ +// $Id$ +//============================================================================== +//! +//! \file FiniteElement.h +//! +//! \date Mar 30 2011 +//! +//! \author Knut Morten Okstad / SINTEF +//! +//! \brief Finite element quantities at an integration point. +//! +//============================================================================== + +#ifndef _FINITE_ELEMENT_H +#define _FINITE_ELEMENT_H + +#include "MatVec.h" + + +/*! + \brief Class representing a finite element. +*/ + +class FiniteElement +{ +public: + //! \brief Default constructor. + FiniteElement(size_t nb = 0) : h(0.0), N(nb), detJxW(1.0) {} + + double u; //!< First parameter of current point + double v; //!< Second parameter of current point + double w; //!< Third parameter of current point + double h; //!< Characteristic element size + Vector N; //!< Basis function values + Vector Navg; //!< Volume-averaged basis function values + Matrix dNdX; //!< First derivatives (gradient) of the basis functions + Matrix3D d2NdX2; //!< Second derivates of the basis functions + double detJxW; //!< Weighted determinant of the coordinate mapping +}; + + +/*! + \brief Class representing a mixed finite element. +*/ + +class MxFiniteElement : public FiniteElement +{ +public: + //! \brief Default constructor. + MxFiniteElement(size_t nb1 = 0, size_t nb2 = 0) : FiniteElement(nb1), + N1(N), N2(nb2), dN1dX(dNdX) {} + + Vector& N1; //!< Basis function values for the first basis + Vector N2; //!< Basis function values for the second basis + Matrix& dN1dX; //!< Basis function gradients for the first basis + Matrix dN2dX; //!< Basis function gradients for the second basis +}; + +#endif diff --git a/src/ASM/IntegrandBase.h b/src/ASM/IntegrandBase.h index b251d877..6c49eb73 100644 --- a/src/ASM/IntegrandBase.h +++ b/src/ASM/IntegrandBase.h @@ -19,6 +19,8 @@ #include "MatVec.h" struct TimeDomain; +class FiniteElement; +class MxFiniteElement; class LocalIntegral; class NormBase; class AnaSol; @@ -134,159 +136,80 @@ public: //! \brief Evaluates the integrand at an interior point. //! \param elmInt The local integral object to receive the contributions + //! \param[in] fe Finite element data of current integration point //! \param[in] time Parameters for nonlinear and time-dependent simulations - //! \param[in] detJW Jacobian determinant times integration point weight - //! \param[in] N Basis function values - //! \param[in] dNdX Basis function gradients //! \param[in] X Cartesian coordinates of current integration point //! - //! \details This interface is used when \a getIntegrandType returns 1. - //! The default implementation forwards to the stationary version. + //! \details The default implementation forwards to the stationary version. //! Reimplement this method for time-dependent or non-linear problems. - virtual bool evalInt(LocalIntegral*& elmInt, - const TimeDomain& time, double detJW, - const Vector& N, const Matrix& dNdX, - const Vec3& X) const + virtual bool evalInt(LocalIntegral*& elmInt, const FiniteElement& fe, + const TimeDomain& time, const Vec3& X) const { - return this->evalInt(elmInt,detJW,N,dNdX,X); + return this->evalInt(elmInt,fe,X); } //! \brief Evaluates the integrand at an interior point. //! \param elmInt The local integral object to receive the contributions - //! \param[in] detJW Jacobian determinant times integration point weight + //! \param[in] fe Mixed finite element data of current integration point //! \param[in] time Parameters for nonlinear and time-dependent simulations - //! \param[in] N Basis function values - //! \param[in] dNdX Basis function gradients - //! \param[in] d2NdX2 Basis function second derivatives - //! \param[in] X Cartesian coordinates of current integration point - //! \param[in] h Characteristic element length - //! - //! \details This interface is used when \a getIntegrandType returns 2. - //! The default implementation forwards to the stationary version. - //! Reimplement this method for time-dependent or non-linear problems. - virtual bool evalInt(LocalIntegral*& elmInt, - const TimeDomain& time, double detJW, - const Vector& N, const Matrix& dNdX, - const Matrix3D& d2NdX2, const Vec3& X, - double h = 0.0) const - { - return this->evalInt(elmInt,detJW,N,dNdX,d2NdX2,X,h); - } - - //! \brief Evaluates the integrand at an interior point. - //! \param elmInt The local integral object to receive the contributions - //! \param[in] detJW Jacobian determinant times integration point weight - //! \param[in] time Parameters for nonlinear and time-dependent simulations - //! \param[in] N Basis function values - //! \param[in] dNdX Basis function gradients - //! \param[in] Navg Volume-averaged basis function values over the element - //! \param[in] X Cartesian coordinates of current integration point - //! - //! \details This interface is used when \a getIntegrandType returns 3. - //! Use when the integrand requires second-derivatives of the basis functions. - //! The default implementation forwards to the stationary version. - //! Reimplement this method for time-dependent or non-linear problems. - virtual bool evalInt(LocalIntegral*& elmInt, - const TimeDomain& time, double detJW, - const Vector& N, const Matrix& dNdX, - const Vector& Navg, const Vec3& X) const - { - return this->evalInt(elmInt,detJW,N,dNdX,Navg,X); - } - - //! \brief Evaluates the integrand at an interior point. - //! \param elmInt The local integral object to receive the contributions - //! \param[in] time Parameters for nonlinear and time-dependent simulations - //! \param[in] detJW Jacobian determinant times integration point weight - //! \param[in] N1 Basis function values, field 1 - //! \param[in] N2 Basis function values, field 2 - //! \param[in] dN1dX Basis function gradients, field 1 - //! \param[in] dN2dX Basis function gradients, field 2 //! \param[in] X Cartesian coordinates of current integration point //! //! \details This interface is used for mixed formulations only. //! The default implementation forwards to the stationary version. //! Reimplement this method for time-dependent or non-linear problems. - virtual bool evalInt(LocalIntegral*& elmInt, - const TimeDomain& time, double detJW, - const Vector& N1, const Vector& N2, - const Matrix& dN1dX, const Matrix& dN2dX, - const Vec3& X) const + virtual bool evalIntMx(LocalIntegral*& elmInt, const MxFiniteElement& fe, + const TimeDomain& time, const Vec3& X) const { - return this->evalInt(elmInt,detJW,N1,N2,dN1dX,dN2dX,X); + return this->evalIntMx(elmInt,fe,X); } //! \brief Evaluates the integrand at a boundary point. //! \param elmInt The local integral object to receive the contributions - //! \param[in] detJW Jacobian determinant times integration point weight + //! \param[in] fe Finite element data of current integration point //! \param[in] time Parameters for nonlinear and time-dependent simulations - //! \param[in] N Basis function values - //! \param[in] dNdX Basis function gradients //! \param[in] X Cartesian coordinates of current integration point //! \param[in] normal Boundary normal vector at current integration point //! //! \details The default implementation forwards to the stationary version. //! Reimplement this method for time-dependent or non-linear problems. - virtual bool evalBou(LocalIntegral*& elmInt, - const TimeDomain& time, double detJW, - const Vector& N, const Matrix& dNdX, + virtual bool evalBou(LocalIntegral*& elmInt, const FiniteElement& fe, + const TimeDomain& time, const Vec3& X, const Vec3& normal) const { - return this->evalBou(elmInt,detJW,N,dNdX,X,normal); + return this->evalBou(elmInt,fe,X,normal); } //! \brief Evaluates the integrand at a boundary point. //! \param elmInt The local integral object to receive the contributions - //! \param[in] detJW Jacobian determinant times integration point weight + //! \param[in] fe Mixed finite element data of current integration point //! \param[in] time Parameters for nonlinear and time-dependent simulations - //! \param[in] N1 Basis function values, field 1 - //! \param[in] N2 Basis function values, field 2 - //! \param[in] dN1dX Basis function gradients, field 1 - //! \param[in] dN2dX Basis function gradients, field 2 //! \param[in] X Cartesian coordinates of current integration point //! \param[in] normal Boundary normal vector at current integration point //! //! \details This interface is used for mixed formulations. //! The default implementation forwards to the stationary version. //! Reimplement this method for time-dependent or non-linear problems. - virtual bool evalBou(LocalIntegral*& elmInt, - const TimeDomain& time, double detJW, - const Vector& N1, const Vector& N2, - const Matrix& dN1dX, const Matrix& dN2dX, - const Vec3& X, const Vec3& normal) const + virtual bool evalBouMx(LocalIntegral*& elmInt, const MxFiniteElement& fe, + const TimeDomain& time, + const Vec3& X, const Vec3& normal) const { - return this->evalBou(elmInt,detJW,N1,N2,dN1dX,dN2dX,X,normal); + return this->evalBouMx(elmInt,fe,X,normal); } protected: //! \brief Evaluates the integrand at interior points for stationary problems. - virtual bool evalInt(LocalIntegral*& elmInt, double detJW, - const Vector& N, const Matrix& dNdX, - const Vec3& X) const { return false; } + virtual bool evalInt(LocalIntegral*&, const FiniteElement& fe, + const Vec3&) const { return false; } //! \brief Evaluates the integrand at interior points for stationary problems. - virtual bool evalInt(LocalIntegral*& elmInt, double detJW, - const Vector& N, const Matrix& dNdX, - const Matrix3D& d2NdX2, const Vec3& X, - double h = 0.0) const { return false; } - //! \brief Evaluates the integrand at interior points for stationary problems. - virtual bool evalInt(LocalIntegral*& elmInt, double detJW, - const Vector& N, const Matrix& dNdX, const Vector& Navg, - const Vec3& X) const { return false; } - //! \brief Evaluates the integrand at interior points for stationary problems. - virtual bool evalInt(LocalIntegral*& elmInt, double detJW, - const Vector& N1, const Vector& N2, - const Matrix& dN1dX, const Matrix& dN2dX, - const Vec3& X) const { return false; } + virtual bool evalIntMx(LocalIntegral*&, const MxFiniteElement& fe, + const Vec3&) const { return false; } //! \brief Evaluates the integrand at boundary points for stationary problems. - virtual bool evalBou(LocalIntegral*& elmInt, double detJW, - const Vector& N, const Matrix& dNdX, const Vec3& X, - const Vec3& normal) const { return false; } + virtual bool evalBou(LocalIntegral*&, const FiniteElement&, + const Vec3&, const Vec3&) const { return false; } //! \brief Evaluates the integrand at boundary points for stationary problems. - virtual bool evalBou(LocalIntegral*& elmInt, double detJW, - const Vector& N1, const Vector& N2, - const Matrix& dN1dX, const Matrix& dN2dX, const Vec3& X, - const Vec3& normal) const { return false; } + virtual bool evalBouMx(LocalIntegral*&, const MxFiniteElement&, + const Vec3&, const Vec3&) const { return false; } public: //! \brief Evaluates the secondary solution at a result point. diff --git a/src/Integrands/Elasticity.C b/src/Integrands/Elasticity.C index 3a746666..5ed429eb 100644 --- a/src/Integrands/Elasticity.C +++ b/src/Integrands/Elasticity.C @@ -12,7 +12,8 @@ //============================================================================== #include "Elasticity.h" -#include "MaterialBase.h" +#include "LinIsotropic.h" +#include "FiniteElement.h" #include "Utilities.h" #include "ElmMats.h" #include "ElmNorm.h" @@ -51,7 +52,13 @@ void Elasticity::print (std::ostream& os) const for (unsigned short int d = 0; d < nsd; d++) std::cout <<" "<< grav[d]; std::cout << std::endl; - if (material) material->print(os); + + if (!material) + { + static LinIsotropic defaultMat; + const_cast(this)->material = &defaultMat; + } + material->print(os); } @@ -575,8 +582,7 @@ ElmNorm& ElasticityNorm::getElmNormBuffer (LocalIntegral*& elmInt) } -bool ElasticityNorm::evalInt (LocalIntegral*& elmInt, double detJW, - const Vector& N, const Matrix& dNdX, +bool ElasticityNorm::evalInt (LocalIntegral*& elmInt, const FiniteElement& fe, const Vec3& X) const { ElmNorm& pnorm = ElasticityNorm::getElmNormBuffer(elmInt); @@ -587,19 +593,19 @@ bool ElasticityNorm::evalInt (LocalIntegral*& elmInt, double detJW, // Evaluate the finite element stress field Vector sigmah, sigma; - if (!problem.evalSol(sigmah,dNdX,X)) return false; + if (!problem.evalSol(sigmah,fe.dNdX,X)) return false; // Integrate the energy norm a(u^h,u^h) - pnorm[0] += sigmah.dot(Cinv*sigmah)*detJW; + pnorm[0] += sigmah.dot(Cinv*sigmah)*fe.detJxW; if (problem.haveLoads()) { // Evaluate the body load Vec3 f = problem.getBodyforce(X); // Evaluate the displacement field - Vec3 u = problem.evalSol(N); + Vec3 u = problem.evalSol(fe.N); // Integrate the external energy (f,u^h) - pnorm[1] += f*u*detJW; + pnorm[1] += f*u*fe.detJxW; } if (anasol) @@ -607,18 +613,17 @@ bool ElasticityNorm::evalInt (LocalIntegral*& elmInt, double detJW, // Evaluate the analytical stress field sigma = (*anasol)(X); // Integrate the energy norm a(u,u) - pnorm[2] += sigma.dot(Cinv*sigma)*detJW; + pnorm[2] += sigma.dot(Cinv*sigma)*fe.detJxW; // Integrate the error in energy norm a(u-u^h,u-u^h) sigma -= sigmah; - pnorm[3] += sigma.dot(Cinv*sigma)*detJW; + pnorm[3] += sigma.dot(Cinv*sigma)*fe.detJxW; } return true; } -bool ElasticityNorm::evalBou (LocalIntegral*& elmInt, double detJW, - const Vector& N, const Matrix&, +bool ElasticityNorm::evalBou (LocalIntegral*& elmInt, const FiniteElement& fe, const Vec3& X, const Vec3& normal) const { if (!problem.haveLoads()) return true; @@ -628,9 +633,9 @@ bool ElasticityNorm::evalBou (LocalIntegral*& elmInt, double detJW, // Evaluate the surface traction Vec3 T = problem.getTraction(X,normal); // Evaluate the displacement field - Vec3 u = problem.evalSol(N); + Vec3 u = problem.evalSol(fe.N); // Integrate the external energy - pnorm[1] += T*u*detJW; + pnorm[1] += T*u*fe.detJxW; return true; } diff --git a/src/Integrands/Elasticity.h b/src/Integrands/Elasticity.h index 21f02339..828e4650 100644 --- a/src/Integrands/Elasticity.h +++ b/src/Integrands/Elasticity.h @@ -240,23 +240,17 @@ public: //! \brief Evaluates the integrand at an interior point. //! \param elmInt The local integral object to receive the contributions - //! \param[in] detJW Jacobian determinant times integration point weight - //! \param[in] N Basis function values - //! \param[in] dNdX Basis function gradients + //! \param[in] fe Finite element data of current integration point //! \param[in] X Cartesian coordinates of current integration point - virtual bool evalInt(LocalIntegral*& elmInt, double detJW, - const Vector& N, const Matrix& dNdX, + virtual bool evalInt(LocalIntegral*& elmInt, const FiniteElement& fe, const Vec3& X) const; //! \brief Evaluates the integrand at a boundary point. //! \param elmInt The local integral object to receive the contributions - //! \param[in] detJW Jacobian determinant times integration point weight - //! \param[in] N Basis function values - //! \param[in] dNdX Basis function gradients + //! \param[in] fe Finite element data of current integration point //! \param[in] X Cartesian coordinates of current integration point //! \param[in] normal Boundary normal vector at current integration point - virtual bool evalBou(LocalIntegral*& elmInt, double detJW, - const Vector& N, const Matrix& dNdX, + virtual bool evalBou(LocalIntegral*& elmInt, const FiniteElement& fe, const Vec3& X, const Vec3& normal) const; //! \brief Returns the number of norm quantities. diff --git a/src/Integrands/LinearElasticity.C b/src/Integrands/LinearElasticity.C index 265e40e9..80d21ff9 100644 --- a/src/Integrands/LinearElasticity.C +++ b/src/Integrands/LinearElasticity.C @@ -12,6 +12,7 @@ //============================================================================== #include "LinearElasticity.h" +#include "FiniteElement.h" #include "MaterialBase.h" #include "Tensor.h" #include "Vec3Oper.h" @@ -24,8 +25,7 @@ LinearElasticity::LinearElasticity (unsigned short int n) : Elasticity(n) } -bool LinearElasticity::evalInt (LocalIntegral*& elmInt, double detJW, - const Vector& N, const Matrix& dNdX, +bool LinearElasticity::evalInt (LocalIntegral*& elmInt, const FiniteElement& fe, const Vec3& X) const { bool lHaveStrains = false; @@ -35,7 +35,7 @@ bool LinearElasticity::evalInt (LocalIntegral*& elmInt, double detJW, { // Compute the strain-displacement matrix B from dNdX // and evaluate the symmetric strain tensor if displacements are available - if (!this->kinematics(dNdX,eps,eps)) return false; + if (!this->kinematics(fe.dNdX,eps,eps)) return false; if (!eps.isZero(1.0e-16)) lHaveStrains = true; // Evaluate the constitutive matrix and the stress tensor at this point @@ -47,36 +47,35 @@ bool LinearElasticity::evalInt (LocalIntegral*& elmInt, double detJW, if (eKm) { // Integrate the material stiffness matrix - CB.multiply(Cmat,Bmat).multiply(detJW); // CB = C*B*|J|*w + CB.multiply(Cmat,Bmat).multiply(fe.detJxW); // CB = C*B*|J|*w eKm->multiply(Bmat,CB,true,false,true); // EK += B^T * CB } if (eKg && lHaveStrains) // Integrate the geometric stiffness matrix - this->formKG(*eKg,dNdX,sigma,detJW); + this->formKG(*eKg,fe.dNdX,sigma,fe.detJxW); if (eM) // Integrate the mass matrix - this->formMassMatrix(*eM,N,X,detJW); + this->formMassMatrix(*eM,fe.N,X,fe.detJxW); if (iS && lHaveStrains) { // Integrate the internal forces - sigma *= -detJW; + sigma *= -fe.detJxW; if (!Bmat.multiply(sigma,*iS,true,true)) // ES -= B^T*sigma return false; } if (eS) // Integrate the load vector due to gravitation and other body forces - this->formBodyForce(*eS,N,X,detJW); + this->formBodyForce(*eS,fe.N,X,fe.detJxW); return this->getIntegralResult(elmInt); } -bool LinearElasticity::evalBou (LocalIntegral*& elmInt, double detJW, - const Vector& N, const Matrix& dNdX, +bool LinearElasticity::evalBou (LocalIntegral*& elmInt, const FiniteElement& fe, const Vec3& X, const Vec3& normal) const { if (!tracFld) @@ -98,9 +97,9 @@ bool LinearElasticity::evalBou (LocalIntegral*& elmInt, double detJW, // Integrate the force vector Vector& ES = *eS; - for (size_t a = 1; a <= N.size(); a++) + for (size_t a = 1; a <= fe.N.size(); a++) for (unsigned short int i = 1; i <= nsd; i++) - ES(nsd*(a-1)+i) += T[i-1]*N(a)*detJW; + ES(nsd*(a-1)+i) += T[i-1]*fe.N(a)*fe.detJxW; return this->getIntegralResult(elmInt); } diff --git a/src/Integrands/LinearElasticity.h b/src/Integrands/LinearElasticity.h index 5b767569..ec285957 100644 --- a/src/Integrands/LinearElasticity.h +++ b/src/Integrands/LinearElasticity.h @@ -32,23 +32,17 @@ public: //! \brief Evaluates the integrand at an interior point. //! \param elmInt The local integral object to receive the contributions - //! \param[in] detJW Jacobian determinant times integration point weight - //! \param[in] N Basis function values - //! \param[in] dNdX Basis function gradients + //! \param[in] fe Finite element data of current integration point //! \param[in] X Cartesian coordinates of current integration point - virtual bool evalInt(LocalIntegral*& elmInt, double detJW, - const Vector& N, const Matrix& dNdX, + virtual bool evalInt(LocalIntegral*& elmInt, const FiniteElement& fe, const Vec3& X) const; //! \brief Evaluates the integrand at a boundary point. //! \param elmInt The local integral object to receive the contributions - //! \param[in] detJW Jacobian determinant times integration point weight - //! \param[in] N Basis function values - //! \param[in] dNdX Basis function gradients + //! \param[in] fe Finite element data of current integration point //! \param[in] X Cartesian coordinates of current integration point //! \param[in] normal Boundary normal vector at current integration point - virtual bool evalBou(LocalIntegral*& elmInt, double detJW, - const Vector& N, const Matrix& dNdX, + virtual bool evalBou(LocalIntegral*& elmInt, const FiniteElement& fe, const Vec3& X, const Vec3& normal) const; private: diff --git a/src/Integrands/Poisson.C b/src/Integrands/Poisson.C index 093d952c..431c4ea0 100644 --- a/src/Integrands/Poisson.C +++ b/src/Integrands/Poisson.C @@ -12,6 +12,7 @@ //============================================================================== #include "Poisson.h" +#include "FiniteElement.h" #include "Utilities.h" #include "ElmNorm.h" #include "Tensor.h" @@ -105,8 +106,7 @@ bool Poisson::initElementBou (const std::vector& MNPC) } -bool Poisson::evalInt (LocalIntegral*& elmInt, double detJW, - const Vector& N, const Matrix& dNdX, +bool Poisson::evalInt (LocalIntegral*& elmInt, const FiniteElement& fe, const Vec3& X) const { elmInt = &myMats; @@ -116,20 +116,19 @@ bool Poisson::evalInt (LocalIntegral*& elmInt, double detJW, // Evaluate the constitutive matrix at this point if (!this->formCmatrix(C,X)) return false; - CB.multiply(C,dNdX,false,true).multiply(detJW); // CB = C*dNdX^T*|J|*w - eM->multiply(dNdX,CB,false,false,true); // EK += dNdX * CB + CB.multiply(C,fe.dNdX,false,true).multiply(fe.detJxW); // = C*dNdX^T*|J|*w + eM->multiply(fe.dNdX,CB,false,false,true); // EK += dNdX * CB } // Integrate heat source if defined if (eS && heatSrc) - eS->add(N,(*heatSrc)(X)*detJW); + eS->add(fe.N,(*heatSrc)(X)*fe.detJxW); return true; } -bool Poisson::evalBou (LocalIntegral*& elmInt, double detJW, - const Vector& N, const Matrix& dNdX, +bool Poisson::evalBou (LocalIntegral*& elmInt, const FiniteElement& fe, const Vec3& X, const Vec3& normal) const { elmInt = &myMats; @@ -155,7 +154,7 @@ bool Poisson::evalBou (LocalIntegral*& elmInt, double detJW, tracVal.insert(std::make_pair(X,trac*normal)); // Integrate the Neumann value - eS->add(N,trac*detJW); + eS->add(fe.N,trac*fe.detJxW); return true; } @@ -278,8 +277,7 @@ bool PoissonNorm::initElement (const std::vector& MNPC) } -bool PoissonNorm::evalInt (LocalIntegral*& elmInt, double detJW, - const Vector& N, const Matrix& dNdX, +bool PoissonNorm::evalInt (LocalIntegral*& elmInt, const FiniteElement& fe, const Vec3& X) const { ElmNorm* eNorm = dynamic_cast(elmInt); @@ -295,20 +293,20 @@ bool PoissonNorm::evalInt (LocalIntegral*& elmInt, double detJW, // Evaluate the finite element heat flux field Vector sigmah; - if (!problem.evalSol(sigmah,dNdX,X)) return false; + if (!problem.evalSol(sigmah,fe.dNdX,X)) return false; // Integrate the energy norm a(u^h,u^h) ElmNorm& pnorm = *eNorm; - pnorm[0] += sigmah.dot(Cinv*sigmah)*detJW; + pnorm[0] += sigmah.dot(Cinv*sigmah)*fe.detJxW; if (anasol) { // Evaluate the analytical heat flux Vector sigma((*anasol)(X).ptr(),sigmah.size()); // Integrate the energy norm a(u,u) - pnorm[1] += sigma.dot(Cinv*sigma)*detJW; + pnorm[1] += sigma.dot(Cinv*sigma)*fe.detJxW; // Integrate the error in energy norm a(u-u^h,u-u^h) sigma -= sigmah; - pnorm[2] += sigma.dot(Cinv*sigma)*detJW; + pnorm[2] += sigma.dot(Cinv*sigma)*fe.detJxW; } return true; diff --git a/src/Integrands/Poisson.h b/src/Integrands/Poisson.h index 40e2e472..745a3858 100644 --- a/src/Integrands/Poisson.h +++ b/src/Integrands/Poisson.h @@ -61,23 +61,17 @@ public: //! \brief Evaluates the integrand at an interior point. //! \param elmInt The local integral object to receive the contributions - //! \param[in] detJW Jacobian determinant times integration point weight - //! \param[in] N Basis function values - //! \param[in] dNdX Basis function gradients + //! \param[in] fe Finite element data of current integration point //! \param[in] X Cartesian coordinates of current integration point - virtual bool evalInt(LocalIntegral*& elmInt, double detJW, - const Vector& N, const Matrix& dNdX, + virtual bool evalInt(LocalIntegral*& elmInt, const FiniteElement& fe, const Vec3& X) const; //! \brief Evaluates the integrand at a boundary point. //! \param elmInt The local integral object to receive the contributions - //! \param[in] detJW Jacobian determinant times integration point weight - //! \param[in] N Basis function values - //! \param[in] dNdX Basis function gradients + //! \param[in] fe Finite element data of current integration point //! \param[in] X Cartesian coordinates of current integration point //! \param[in] normal Boundary normal vector at current integration point - virtual bool evalBou(LocalIntegral*& elmInt, double detJW, - const Vector& N, const Matrix& dNdX, + virtual bool evalBou(LocalIntegral*& elmInt, const FiniteElement& fe, const Vec3& X, const Vec3& normal) const; //! \brief Evaluates the secondary solution at a result point. @@ -178,12 +172,9 @@ public: //! \brief Evaluates the integrand at an interior point. //! \param elmInt The local integral object to receive the contributions - //! \param[in] detJW Jacobian determinant times integration point weight - //! \param[in] N Basis function values - //! \param[in] dNdX Basis function gradients + //! \param[in] fe Finite element data of current integration point //! \param[in] X Cartesian coordinates of current integration point - virtual bool evalInt(LocalIntegral*& elmInt, double detJW, - const Vector& N, const Matrix& dNdX, + virtual bool evalInt(LocalIntegral*& elmInt, const FiniteElement& fe, const Vec3& X) const; //! \brief Returns the number of norm quantities.