diff --git a/Apps/FiniteDefElasticity/CMakeLists.txt b/Apps/FiniteDefElasticity/CMakeLists.txt index cf7feb0f..d88be656 100644 --- a/Apps/FiniteDefElasticity/CMakeLists.txt +++ b/Apps/FiniteDefElasticity/CMakeLists.txt @@ -78,6 +78,14 @@ IF(NOT "${DISABLE_HDF5}") FIND_PACKAGE(HDF5) ENDIF(NOT "${DISABLE_HDF5}") +IF(NOT "${DISABLE_OPENMP}") + FIND_PACKAGE(OpenMP) +ENDIF(NOT "${DISABLE_OPENMP}") + +IF(OPENMP_FOUND) + SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS} -DUSE_OPENMP") +ENDIF(OPENMP_FOUND) + # Required libraries SET(DEPLIBS ${IFEM_LIBRARIES} ${GoTrivariate_LIBRARIES} ${GoTools_LIBRARIES} diff --git a/Apps/FiniteDefElasticity/NonlinearElasticity.C b/Apps/FiniteDefElasticity/NonlinearElasticity.C index b136b8e4..73664942 100644 --- a/Apps/FiniteDefElasticity/NonlinearElasticity.C +++ b/Apps/FiniteDefElasticity/NonlinearElasticity.C @@ -14,6 +14,7 @@ #include "NonlinearElasticity.h" #include "MaterialBase.h" #include "FiniteElement.h" +#include "ElmMats.h" #include "Utilities.h" #include "Profiler.h" @@ -54,7 +55,7 @@ extern "C" { NonlinearElasticity::NonlinearElasticity (unsigned short int n) - : NonlinearElasticityTL(n), E(n) + : NonlinearElasticityTL(n) { fullCmat = false; } @@ -76,22 +77,27 @@ void NonlinearElasticity::setMode (SIM::SolutionMode mode) } -bool NonlinearElasticity::evalInt (LocalIntegral*& elmInt, +bool NonlinearElasticity::evalInt (LocalIntegral& elmInt, const FiniteElement& fe, const Vec3& X) const { PROFILE3("NonlinearEl::evalInt"); + ElmMats& elMat = static_cast(elmInt); + // Evaluate the kinematic quantities, F and E, at this point + Matrix Bmat; Tensor F(nsd); - if (!this->kinematics(fe.N,fe.dNdX,X.x,F,E)) + SymmTensor E(nsd); + if (!this->kinematics(elMat.vec.front(),fe.N,fe.dNdX,X.x,Bmat,F,E)) return false; // Evaluate current tangent at this point, that is // the incremental constitutive matrix, Cmat, and // the 2nd Piola-Kirchhoff stress tensor, S + Matrix Cmat; SymmTensor S(nsd); - if (!this->formTangent(Cmat,S,X,F)) + if (!this->formTangent(Cmat,S,X,F,E)) return false; bool haveStrains = !E.isZero(1.0e-16); @@ -102,33 +108,33 @@ bool NonlinearElasticity::evalInt (LocalIntegral*& elmInt, if (eKm) { // Integrate the material stiffness matrix + Matrix& EM = elMat.A[eKm-1]; #ifndef NO_FORTRAN // Using Fortran routines optimized for symmetric constitutive tensors PROFILE4("stiff_TL_"); if (nsd == 3) if (fullCmat) stiff_tl3d_(fe.N.size(),fe.detJxW,fe.dNdX.ptr(),F.ptr(), - Cmat.ptr(),eKm->ptr()); + Cmat.ptr(),EM.ptr()); else stiff_tl3d_isoel_(fe.N.size(),fe.detJxW,fe.dNdX.ptr(),F.ptr(), - Cmat(1,1),Cmat(1,2),Cmat(4,4),eKm->ptr()); + Cmat(1,1),Cmat(1,2),Cmat(4,4),EM.ptr()); else if (nsd == 2) if (fullCmat) stiff_tl2d_(fe.N.size(),fe.detJxW,fe.dNdX.ptr(),F.ptr(), - Cmat.ptr(),eKm->ptr()); + Cmat.ptr(),EM.ptr()); else stiff_tl2d_isoel_(fe.N.size(),fe.detJxW,fe.dNdX.ptr(),F.ptr(), - Cmat(1,1),Cmat(1,2),Cmat(3,3),eKm->ptr()); + Cmat(1,1),Cmat(1,2),Cmat(3,3),EM.ptr()); else if (nsd == 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); + EM(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 <= fe.N.size(); a++) for (b = 1; b <= fe.N.size(); b++) for (m = 1; m <= nsd; m++) @@ -148,16 +154,16 @@ bool NonlinearElasticity::evalInt (LocalIntegral*& elmInt, if (eKg && haveStrains) // Integrate the geometric stiffness matrix - this->formKG(*eKg,fe.N,fe.dNdX,X.x,S,fe.detJxW); + this->formKG(elMat.A[eKg-1],fe.N,fe.dNdX,X.x,S,fe.detJxW); if (eM) // Integrate the mass matrix - this->formMassMatrix(*eM,fe.N,X,fe.detJxW); + this->formMassMatrix(elMat.A[eM-1],fe.N,X,fe.detJxW); if (iS && haveStrains) { // Integrate the internal forces - Vector& ES = *iS; + Vector& ES = elMat.b[iS-1]; for (a = 1; a <= fe.N.size(); a++) for (k = 1; k <= nsd; k++) { @@ -171,9 +177,9 @@ bool NonlinearElasticity::evalInt (LocalIntegral*& elmInt, if (eS) // Integrate the load vector due to gravitation and other body forces - this->formBodyForce(*eS,fe.N,X,fe.detJxW); + this->formBodyForce(elMat.b[eS-1],fe.N,X,fe.detJxW); - return this->getIntegralResult(elmInt); + return true; } @@ -183,18 +189,20 @@ bool NonlinearElasticity::evalSol (Vector& s, const Vector&, { PROFILE3("NonlinearEl::evalSol"); + // Extract element displacements + Vector eV; int ierr = 0; - if (eV && !primsol.front().empty()) - if ((ierr = utl::gather(MNPC,nsd,primsol.front(),*eV))) + if (!primsol.front().empty()) + if ((ierr = utl::gather(MNPC,nsd,primsol.front(),eV))) { - std::cerr <<" *** NonlinearElasticity::evalSol: Detected " - << ierr <<" node numbers out of range."<< std::endl; + std::cerr <<" *** NonlinearElasticity::evalSol: Detected "<< ierr + <<" node numbers out of range."<< std::endl; return false; } // Evaluate the stress state at this point SymmTensor Sigma(nsd); - if (!this->formStressTensor(dNdX,X,Sigma)) + if (!this->formStressTensor(eV,dNdX,X,Sigma)) return false; // Congruence transformation to local coordinate system at current point @@ -206,13 +214,13 @@ bool NonlinearElasticity::evalSol (Vector& s, const Vector&, } -bool NonlinearElasticity::evalSol (Vector& s, +bool NonlinearElasticity::evalSol (Vector& s, const Vector& eV, const Matrix& dNdX, const Vec3& X) const { PROFILE3("NonlinearEl::evalSol"); SymmTensor Sigma(nsd); - if (!this->formStressTensor(dNdX,X,Sigma)) + if (!this->formStressTensor(eV,dNdX,X,Sigma)) return false; s = Sigma; @@ -220,10 +228,11 @@ bool NonlinearElasticity::evalSol (Vector& s, } -bool NonlinearElasticity::formStressTensor (const Matrix& dNdX, const Vec3& X, +bool NonlinearElasticity::formStressTensor (const Vector& eV, + const Matrix& dNdX, const Vec3& X, SymmTensor& S) const { - if (!eV || eV->empty()) + if (eV.empty()) { // Initial state (zero stresses) S.zero(); @@ -231,18 +240,22 @@ bool NonlinearElasticity::formStressTensor (const Matrix& dNdX, const Vec3& X, } // Evaluate the kinematic quantities, F and E, at this point + Matrix B; Tensor F(nsd); - if (!this->kinematics(Vector(),dNdX,X.x,F,E)) + SymmTensor E(nsd); + if (!this->kinematics(eV,Vector(),dNdX,X.x,B,F,E)) return false; // Evaluate the 2nd Piola-Kirchhoff stress tensor, S, at this point + Matrix Cmat; double U; return material->evaluate(Cmat,S,U,X,F,E,2); } bool NonlinearElasticity::formTangent (Matrix& Ctan, SymmTensor& S, - const Vec3& X, const Tensor& F) const + const Vec3& X, const Tensor& F, + const SymmTensor& E) const { double U; return material->evaluate(Ctan,S,U,X,F,E, E.isZero(1.0e-16) ? 0 : 2); diff --git a/Apps/FiniteDefElasticity/NonlinearElasticity.h b/Apps/FiniteDefElasticity/NonlinearElasticity.h index 84d5dfe0..08fb25bf 100644 --- a/Apps/FiniteDefElasticity/NonlinearElasticity.h +++ b/Apps/FiniteDefElasticity/NonlinearElasticity.h @@ -51,7 +51,7 @@ public: //! \param elmInt The local integral object to receive the contributions //! \param[in] fe Finite element data of current integration point //! \param[in] X Cartesian coordinates of current integration point - virtual bool evalInt(LocalIntegral*& elmInt, const FiniteElement& fe, + virtual bool evalInt(LocalIntegral& elmInt, const FiniteElement& fe, const Vec3& X) const; //! \brief Evaluates the secondary solution at a result point. @@ -66,9 +66,11 @@ public: //! \brief Evaluates the finite element (FE) solution at an integration point. //! \param[out] s The FE stress values at current point + //! \param[in] eV Element solution vector //! \param[in] dNdX Basis function gradients at current point //! \param[in] X Cartesian coordinates of current point - virtual bool evalSol(Vector& s, const Matrix& dNdX, const Vec3& X) const; + virtual bool evalSol(Vector& s, const Vector& eV, + const Matrix& dNdX, const Vec3& X) const; protected: //! \brief Forms tangential tensorial quantities needed by the evalInt method. @@ -76,19 +78,22 @@ protected: //! \param[out] S 2nd Piola-Kirchhoff stress tensor at current point //! \param[in] X Cartesian coordinates of current integration point //! \param[in] F Deformation gradient at current integration point + //! \param[in] E Green-Lagrange strain tensor at current integration point virtual bool formTangent(Matrix& Ctan, SymmTensor& S, - const Vec3& X, const Tensor& F) const; + const Vec3& X, const Tensor& F, + const SymmTensor& E) const; //! \brief Forms the 2nd Piola-Kirchhoff stress tensor. + //! \param[in] eV Element solution vector //! \param[in] dNdX Basis function gradients at current integration point //! \param[in] X Cartesian coordinates of current integration point //! \param[out] S 2nd Piola-Kirchhoff stress tensor at current point - virtual bool formStressTensor(const Matrix& dNdX, const Vec3& X, + virtual bool formStressTensor(const Vector& eV, + const Matrix& dNdX, const Vec3& X, SymmTensor& S) const; protected: - bool fullCmat; //!< If \e true, assume a full (but symmetric) C-matrix - mutable SymmTensor E; //!< Green-Lagrange strain tensor + bool fullCmat; //!< If \e true, assume a full (but symmetric) C-matrix }; #endif diff --git a/Apps/FiniteDefElasticity/NonlinearElasticityFbar.C b/Apps/FiniteDefElasticity/NonlinearElasticityFbar.C index 2d6b68cb..f762ab01 100644 --- a/Apps/FiniteDefElasticity/NonlinearElasticityFbar.C +++ b/Apps/FiniteDefElasticity/NonlinearElasticityFbar.C @@ -15,18 +15,97 @@ #include "MaterialBase.h" #include "FiniteElement.h" #include "GaussQuadrature.h" +#include "ElmMats.h" +#include "ElmNorm.h" #include "Lagrange.h" #include "TimeDomain.h" #include "Tensor.h" +/*! + \brief A struct with volumetric sampling point data. +*/ + +struct VolPtData +{ + double J; //!< Determinant of current deformation gradient + Vector Nr; //!< Basis function values (for axisymmetric problems) + Matrix dNdx; //!< Basis function gradients at current configuration + //! \brief Default constructor. + VolPtData() { J = 1.0; } +}; + + +/*! + \brief Class containing internal data for an Fbar element. +*/ + +class FbarElmData +{ +public: + //! \brief Default constructor. + FbarElmData() : pbar(0), scale(0.0), iP(0) {} + //! \brief Empty destructor. + virtual ~FbarElmData() {} + + //! \brief Initializes the element data + bool init(unsigned short int nsd, size_t nPt) + { + pbar = ceil(pow((double)nPt,1.0/(double)nsd)-0.5); + if (pow((double)pbar,(double)nsd) == (double)nPt) + myVolData.resize(nPt); + else + { + pbar = 0; + std::cerr <<" *** FbarElmData::init: Invalid element, "<< nPt + <<" volumetric sampling points specified."<< std::endl; + return false; + } + scale = pbar > 1 ? 1.0/GaussQuadrature::getCoord(pbar)[pbar-1] : 1.0; + iP = 0; + return true; + } + + int pbar; //!< Polynomial order of the internal volumetric data field + double scale; //!< Scaling factor for extrapolation from sampling points + size_t iP; //!< Volumetric sampling point counter + + std::vector myVolData; //!< Volumetric sampling point data +}; + + +/*! + \brief Class collecting element matrices associated with a Fbar FEM problem. +*/ + +class FbarMats : public FbarElmData, public ElmMats {}; + + +/*! + \brief Class representing integrated norm quantities for an Fbar element. +*/ + +class FbarNorm : public FbarElmData, public LocalIntegral +{ +public: + //! \brief The constructor initializes the ElmNorm pointer. + FbarNorm(LocalIntegral& n) { myNorm = static_cast(&n); } + //! \brief Alternative constructor allocating an internal ElmNorm. + FbarNorm(size_t n) { myNorm = new ElmNorm(n); } + //! \brief The destructor destroys the ElmNorm object. + virtual ~FbarNorm() { myNorm->destruct(); } + + //! \brief Returns the LocalIntegral object to assemble into the global one. + virtual const LocalIntegral* ref() const { return myNorm; } + + ElmNorm* myNorm; //!< Pointer to the actual element norms object +}; + + NonlinearElasticityFbar::NonlinearElasticityFbar (unsigned short int n, bool axS, int nvp) : NonlinearElasticityUL(n,axS), npt1(nvp) { - iP = 0; - pbar = 0; - scale = 1.0; } @@ -39,48 +118,104 @@ void NonlinearElasticityFbar::print (std::ostream& os) const } -bool NonlinearElasticityFbar::initElement (const std::vector& MNPC, - const Vec3&, size_t nPt) +LocalIntegral* NonlinearElasticityFbar::getLocalIntegral (size_t nen, size_t, + bool neumann) const { - iP = 0; - pbar = ceil(pow((double)nPt,1.0/(double)nsd)-0.5); - if (pow((double)pbar,(double)nsd) == (double)nPt) - myVolData.resize(nPt); - else + ElmMats* result = new FbarMats; + switch (m_mode) { - pbar = 0; - std::cerr <<" *** NonlinearElasticityFbar::initElement: Invalid element, " - << nPt <<" volumetric sampling points specified."<< std::endl; - return false; - } - scale = pbar > 1 ? 1.0/GaussQuadrature::getCoord(pbar)[pbar-1] : 1.0; -#if INT_DEBUG > 0 - std::cout <<"\n\n *** Entering NonlinearElasticityFbar::initElement"; - std::cout <<", nPt = "<< nPt <<", pbar = "<< pbar - <<", scale = "<< scale << std::endl; -#endif + case SIM::STATIC: + result->rhsOnly = neumann; + result->withLHS = !neumann; + result->resize(neumann?0:1,1); + break; - return this->NonlinearElasticityUL::initElement(MNPC); + case SIM::DYNAMIC: + result->rhsOnly = neumann; + result->withLHS = !neumann; + result->resize(neumann?0:2,1); + break; + + case SIM::VIBRATION: + case SIM::BUCKLING: + result->withLHS = true; + result->resize(2,0); + break; + + case SIM::STIFF_ONLY: + case SIM::MASS_ONLY: + result->withLHS = true; + result->resize(1,0); + break; + + case SIM::RHS_ONLY: + result->rhsOnly = true; + result->resize(neumann?0:1,1); + break; + + case SIM::RECOVERY: + result->rhsOnly = true; + break; + + default: + ; + } + + for (size_t i = 0; i < result->A.size(); i++) + result->A[i].resize(nsd*nen,nsd*nen); + + if (result->b.size()) + result->b.front().resize(nsd*nen); + + return result; } -bool NonlinearElasticityFbar::reducedInt (const FiniteElement& fe, +bool NonlinearElasticityFbar::initElement (const std::vector& MNPC, + const Vec3&, size_t nPt, + LocalIntegral& elmInt) +{ + FbarMats& fbar = static_cast(elmInt); + +#if INT_DEBUG > 0 + std::cout <<"\n\n *** Entering NonlinearElasticityFbar::initElement"; + std::cout <<", nPt = "<< nPt <<", pbar = "<< fbar.pbar + <<", scale = "<< fbar.scale << std::endl; +#endif + return fbar.init(nsd,nPt) && this->IntegrandBase::initElement(MNPC,elmInt); +} + + +bool NonlinearElasticityFbar::reducedInt (LocalIntegral& elmInt, + const FiniteElement& fe, const Vec3& X) const { + FbarNorm* fbNorm = 0; + FbarElmData* fbar = dynamic_cast(&elmInt); + if (!fbar) + { + if ((fbNorm = dynamic_cast(&elmInt))) + fbar = fbNorm; + else + return false; + } + #if INT_DEBUG > 1 std::cout <<"NonlinearElasticityFbar::u(red) = "<< fe.u; if (nsd > 1) std::cout <<" "<< fe.v; if (nsd > 2) std::cout <<" "<< fe.w; - std::cout <<" (iP="<< iP+1 <<")\n" + std::cout <<" (iP="<< fbar->iP+1 <<")\n" <<"NonlinearElasticityFbar::dNdX ="<< fe.dNdX; #endif - VolPtData& ptData = myVolData[iP++]; + const Vector& eV = fbNorm ? fbNorm->myNorm->vec.front() : elmInt.vec.front(); + VolPtData& ptData = fbar->myVolData[fbar->iP++]; // Evaluate the deformation gradient, F, at current configuration + Matrix B; Tensor F(nDF); SymmTensor E(nsd,axiSymmetry); - if (!this->kinematics(fe.N,fe.dNdX,X.x,F,E)) + if (!this->kinematics(eV,fe.N,fe.dNdX,X.x,B,F,E)) return false; if (E.isZero(1.0e-16)) @@ -110,7 +245,7 @@ bool NonlinearElasticityFbar::reducedInt (const FiniteElement& fe, // Push-forward the basis function gradients to current configuration ptData.dNdx.multiply(fe.dNdX,Fi); // dNdx = dNdX * F^-1 if (axiSymmetry && X.x > 0.0) - ptData.Nr = fe.N * (1.0/(X.x + eV->dot(fe.N,0,nsd))); + ptData.Nr = fe.N * (1.0/(X.x + eV.dot(fe.N,0,nsd))); #ifdef INT_DEBUG std::cout <<"NonlinearElasticityFbar::J = "<< ptData.J @@ -120,7 +255,7 @@ bool NonlinearElasticityFbar::reducedInt (const FiniteElement& fe, #endif } #ifdef INT_DEBUG - if (iP == myVolData.size()) + if (fbar->iP == fbar->myVolData.size()) std::cout <<"NonlinearElasticityFbar: Volumetric sampling points completed." <<"\n"<< std::endl; #endif @@ -326,7 +461,7 @@ static void getAmat3D (const Matrix& C, const SymmTensor& Sig, Matrix& A) } -bool NonlinearElasticityFbar::evalInt (LocalIntegral*& elmInt, +bool NonlinearElasticityFbar::evalInt (LocalIntegral& elmInt, const FiniteElement& fe, const TimeDomain& prm, const Vec3& X) const @@ -338,10 +473,13 @@ bool NonlinearElasticityFbar::evalInt (LocalIntegral*& elmInt, std::cout <<"\nNonlinearElasticityFbar::dNdX ="<< fe.dNdX; #endif + FbarMats& fbar = static_cast(elmInt); + // Evaluate the deformation gradient, F, at current configuration + Matrix B, dNdx; Tensor F(nDF); SymmTensor E(nsd,axiSymmetry); - if (!this->kinematics(fe.N,fe.dNdX,X.x,F,E)) + if (!this->kinematics(fbar.vec.front(),fe.N,fe.dNdX,X.x,B,F,E)) return false; double J, Jbar = 0.0; @@ -380,23 +518,27 @@ bool NonlinearElasticityFbar::evalInt (LocalIntegral*& elmInt, // Axi-symmetric integration point volume; 2*pi*r*|J|*w double detJW = (axiSymmetry ? 2.0*M_PI*X.x : 1.0)*fe.detJxW*J; - double r = axiSymmetry ? X.x + eV->dot(fe.N,0,nsd) : 0.0; + double r = axiSymmetry ? X.x + fbar.vec.front().dot(fe.N,0,nsd) : 0.0; - if (myVolData.size() == 1) + Vector M; + Matrix dMdx; + if (fbar.myVolData.size() == 1) { // Only one volumetric sampling point (linear element) - Jbar = myVolData.front().J; - dMdx = myVolData.front().dNdx; - if (axiSymmetry) M = myVolData.front().Nr; + Jbar = fbar.myVolData.front().J; + dMdx = fbar.myVolData.front().dNdx; + if (axiSymmetry) M = fbar.myVolData.front().Nr; } - else if (!myVolData.empty()) + else if (!fbar.myVolData.empty()) { // Evaluate the Lagrange polynomials extrapolating the volumetric // sampling points (assuming a regular distribution over the element) Vector Nbar; - int pbar3 = nsd > 2 ? pbar : 0; - if (!Lagrange::computeBasis(Nbar,pbar,fe.xi*scale,pbar,fe.eta*scale, - pbar3,fe.zeta*scale)) return false; + int pbar3 = nsd > 2 ? fbar.pbar : 0; + if (!Lagrange::computeBasis(Nbar, + fbar.pbar,fe.xi*fbar.scale, + fbar.pbar,fe.eta*fbar.scale, + pbar3,fe.zeta*fbar.scale)) return false; #ifdef INT_DEBUG std::cout <<"NonlinearElasticityFbar::Nbar ="<< Nbar; #endif @@ -408,10 +550,10 @@ bool NonlinearElasticityFbar::evalInt (LocalIntegral*& elmInt, if (axiSymmetry) M.resize(dNdx.rows(),true); for (size_t a = 0; a < Nbar.size(); a++) { - Jbar += myVolData[a].J*Nbar[a]; - dMdx.add(myVolData[a].dNdx,myVolData[a].J*Nbar[a]); + Jbar += fbar.myVolData[a].J*Nbar[a]; + dMdx.add(fbar.myVolData[a].dNdx,fbar.myVolData[a].J*Nbar[a]); if (axiSymmetry) - M.add(myVolData[a].Nr,myVolData[a].J*Nbar[a]); + M.add(fbar.myVolData[a].Nr,fbar.myVolData[a].J*Nbar[a]); } dMdx.multiply(1.0/Jbar); if (axiSymmetry) M /= Jbar; @@ -434,6 +576,7 @@ bool NonlinearElasticityFbar::evalInt (LocalIntegral*& elmInt, #endif // Evaluate the constitutive relation (Jbar is dummy here) + Matrix Cmat; SymmTensor Sig(nsd,axiSymmetry); if (!material->evaluate(Cmat,Sig,Jbar,X,F,E,lHaveStrains,&prm)) return false; @@ -445,7 +588,7 @@ bool NonlinearElasticityFbar::evalInt (LocalIntegral*& elmInt, if (iS && lHaveStrains) { // Compute and accumulate contribution to internal force vector - Vector& ES = *iS; + Vector& ES = fbar.b[iS-1]; size_t a, d; unsigned short int i, j; for (a = d = 1; a <= dNdx.rows(); a++) @@ -465,6 +608,7 @@ bool NonlinearElasticityFbar::evalInt (LocalIntegral*& elmInt, if (eKm) { // Compute standard and modified discrete spatial gradient operators + Matrix G, Gbar; getGradOperator(axiSymmetry ? r : -1.0, fe.N, dNdx, G); getGradOperator(axiSymmetry ? 1.0 : -1.0, M, dMdx, Gbar); @@ -496,78 +640,113 @@ bool NonlinearElasticityFbar::evalInt (LocalIntegral*& elmInt, } #ifdef INT_DEBUG - std::cout <<"NonlinearElasticityFbar::G ="<< G - <<"NonlinearElasticityFbar::Gbar ="<< Gbar - <<"NonlinearElasticityFbar::A ="<< A - <<"NonlinearElasticityFbar::Q ="<< Q; + std::cout <<"NonlinearElasticityFbar::G ="<< G + <<"NonlinearElasticityFbar::Gbar ="<< Gbar + <<"NonlinearElasticityFbar::A ="<< A + <<"NonlinearElasticityFbar::Q ="<< Q; #endif // Compute and accumulate contribution to tangent stiffness matrix. // First, standard (material and geometric) tangent stiffness terms - eKm->multiply(G,CB.multiply(A,G),true,false,true); // EK += G^T * A*G + Matrix CB; CB.multiply(A,G); + fbar.A[eKm-1].multiply(G,CB,true,false,true); // EK += G^T * A*G // Then, modify the tangent stiffness for the F-bar terms - Gbar -= G; - eKm->multiply(G,CB.multiply(Q,Gbar),true,false,true); // EK += G^T * Q*Gbar + Gbar -= G; CB.multiply(Q,Gbar); + fbar.A[eKm-1].multiply(G,CB,true,false,true); // EK += G^T * Q*Gbar } if (eM) // Integrate the mass matrix - this->formMassMatrix(*eM,fe.N,X,detJW); + this->formMassMatrix(fbar.A[eM-1],fe.N,X,detJW); if (eS) // Integrate the load vector due to gravitation and other body forces - this->formBodyForce(*eS,fe.N,X,detJW); + this->formBodyForce(fbar.b[eS-1],fe.N,X,detJW); - return this->getIntegralResult(elmInt); + return true; } - NormBase* NonlinearElasticityFbar::getNormIntegrand (AnaSol*) const { return new ElasticityNormFbar(*const_cast(this)); } -bool ElasticityNormFbar::reducedInt (const FiniteElement& fe, - const Vec3& X) const +LocalIntegral* ElasticityNormFbar::getLocalIntegral (size_t nen, size_t iEl, + bool neumann) const { - return myProblem.reducedInt(fe,X); + LocalIntegral* result = this->NormBase::getLocalIntegral(nen,iEl,neumann); + if (result) return new FbarNorm(static_cast(*result)); + + // Element norms are not requested, so allocate an internal object instead + // that will delete itself when invoking its destruct method. + return new FbarNorm(this->getNoFields()); } -bool ElasticityNormFbar::evalInt (LocalIntegral*& elmInt, +bool ElasticityNormFbar::initElement (const std::vector& MNPC, + const Vec3&, size_t nPt, + LocalIntegral& elmInt) +{ + NonlinearElasticityFbar& p = static_cast(myProblem); + FbarNorm& fbar = static_cast(elmInt); + + return fbar.init(p.nsd,nPt) && this->NormBase::initElement(MNPC,*fbar.myNorm); +} + + +bool ElasticityNormFbar::initElementBou (const std::vector& MNPC, + LocalIntegral& elmInt) +{ + return myProblem.initElementBou(MNPC,*static_cast(elmInt).myNorm); +} + + +bool ElasticityNormFbar::reducedInt (LocalIntegral& elmInt, + const FiniteElement& fe, + const Vec3& X) const +{ + return myProblem.reducedInt(elmInt,fe,X); +} + + +bool ElasticityNormFbar::evalInt (LocalIntegral& elmInt, const FiniteElement& fe, const TimeDomain& prm, const Vec3& X) const { size_t nsd = fe.dNdX.cols(); NonlinearElasticityFbar& p = static_cast(myProblem); + FbarNorm& fbar = static_cast(elmInt); // Evaluate the deformation gradient, F, and the Green-Lagrange strains, E + Matrix B; Tensor F(p.nDF); SymmTensor E(p.nDF); - if (!p.kinematics(fe.N,fe.dNdX,X.x,F,E)) + if (!p.kinematics(fbar.myNorm->vec.front(),fe.N,fe.dNdX,X.x,B,F,E)) return false; double Jbar = 0.0; - if (p.myVolData.size() == 1) - Jbar = p.myVolData.front().J; - else if (!p.myVolData.empty()) + if (fbar.myVolData.size() == 1) + Jbar = fbar.myVolData.front().J; + else if (!fbar.myVolData.empty()) { // Evaluate the Lagrange polynomials extrapolating the volumetric // sampling points (assuming a regular distribution over the element) RealArray Nbar; - int pbar3 = nsd > 2 ? p.pbar : 0; - if (!Lagrange::computeBasis(Nbar,p.pbar,fe.xi*p.scale,p.pbar,fe.eta*p.scale, - pbar3,fe.zeta*p.scale)) return false; + int pbar3 = nsd > 2 ? fbar.pbar : 0; + if (!Lagrange::computeBasis(Nbar, + fbar.pbar,fe.xi*fbar.scale, + fbar.pbar,fe.eta*fbar.scale, + pbar3,fe.zeta*fbar.scale)) return false; // Compute modified deformation gradient determinant // by extrapolating the volume sampling points for (size_t a = 0; a < Nbar.size(); a++) - Jbar += p.myVolData[a].J*Nbar[a]; + Jbar += fbar.myVolData[a].J*Nbar[a]; } else Jbar = F.det(); @@ -578,15 +757,23 @@ bool ElasticityNormFbar::evalInt (LocalIntegral*& elmInt, // Compute the strain energy density, U(E) = Int_E (S:Eps) dEps // and the Cauchy stress tensor, sigma - double U = 0.0; + Matrix Cmat; double U = 0.0; SymmTensor sigma(nsd, p.isAxiSymmetric() || p.material->isPlaneStrain()); - if (!p.material->evaluate(p.Cmat,sigma,U,X,Fbar,E,3,&prm,&F)) + if (!p.material->evaluate(Cmat,sigma,U,X,Fbar,E,3,&prm,&F)) return false; // Axi-symmetric integration point volume; 2*pi*r*|J|*w double detJW = p.isAxiSymmetric() ? 2.0*M_PI*X.x*fe.detJxW : fe.detJxW; // Integrate the norms - return ElasticityNormUL::evalInt(getElmNormBuffer(elmInt,6), - sigma,U,F.det(),detJW); + return ElasticityNormUL::evalInt(*fbar.myNorm,sigma,U,F.det(),detJW); +} + + +bool ElasticityNormFbar::evalBou (LocalIntegral& elmInt, + const FiniteElement& fe, + const Vec3& X, const Vec3& normal) const +{ + return this->ElasticityNormUL::evalBou(*static_cast(elmInt).myNorm, + fe,X,normal); } diff --git a/Apps/FiniteDefElasticity/NonlinearElasticityFbar.h b/Apps/FiniteDefElasticity/NonlinearElasticityFbar.h index 2e7dc482..7ee1656e 100644 --- a/Apps/FiniteDefElasticity/NonlinearElasticityFbar.h +++ b/Apps/FiniteDefElasticity/NonlinearElasticityFbar.h @@ -15,6 +15,8 @@ #define _NONLINEAR_ELASTICITY_FBAR_H #include "NonlinearElasticityUL.h" +#include "ElmMats.h" +#include "ElmNorm.h" /*! @@ -24,16 +26,6 @@ class NonlinearElasticityFbar : public NonlinearElasticityUL { - //! \brief A struct with volumetric sampling point data. - struct VolPtData - { - double J; //!< Determinant of current deformation gradient - Vector Nr; //!< Basis function values (for axisymmetric problems) - Matrix dNdx; //!< Basis function gradients at current configuration - //! \brief Default constructor. - VolPtData() { J = 1.0; } - }; - public: //! \brief The default constructor invokes the parent class constructor. //! \param[in] n Number of spatial dimensions @@ -47,23 +39,33 @@ public: //! \brief Prints out problem definition to the given output stream. virtual void print(std::ostream& os) const; + //! \brief Returns a local integral contribution object for the given element. + //! \param[in] nen Number of nodes on element + //! \param[in] neumann Whether or not we are assembling Neumann BC's + virtual LocalIntegral* getLocalIntegral(size_t nen, size_t, + bool neumann) const; + //! \brief Initializes current element for numerical integration. //! \param[in] MNPC Matrix of nodal point correspondance for current element //! \param[in] nPt Number of volumetric sampling points in this element + //! \param elmInt The local integral object for current element virtual bool initElement(const std::vector& MNPC, - const Vec3&, size_t nPt); + const Vec3&, size_t nPt, + LocalIntegral& elmInt); //! \brief Evaluates volumetric sampling point data at an interior point. + //! \param elmInt The local integral object to receive the contributions //! \param[in] fe Finite element data of current point //! \param[in] X Cartesian coordinates of current integration point - virtual bool reducedInt(const FiniteElement& fe, const Vec3& X) const; + virtual bool reducedInt(LocalIntegral& elmInt, const FiniteElement& fe, + const Vec3& X) const; //! \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] X Cartesian coordinates of current integration point - virtual bool evalInt(LocalIntegral*& elmInt, const FiniteElement& fe, + virtual bool evalInt(LocalIntegral& elmInt, const FiniteElement& fe, const TimeDomain& prm, const Vec3& X) const; //! \brief Returns a pointer to an Integrand for solution norm evaluation. @@ -76,17 +78,7 @@ public: virtual int getIntegrandType() const { return 10+npt1; } private: - int npt1; //!< Number of volumetric sampling points in each direction - int pbar; //!< Polynomial order of the internal volumetric data field - double scale; //!< Scaling factor for extrapolation from sampling points - - mutable std::vector myVolData; //!< Volumetric sampling point data - - mutable size_t iP; //!< Volumetric sampling point counter - mutable Vector M; //!< Modified basis function values - mutable Matrix dMdx; //!< Modified basis function gradients - mutable Matrix G; //!< Discrete gradient operator - mutable Matrix Gbar; //!< Modified discrete gradient operator + int npt1; //!< Number of volumetric sampling points in each direction friend class ElasticityNormFbar; }; @@ -105,19 +97,50 @@ public: //! \brief Empty destructor. virtual ~ElasticityNormFbar() {} + //! \brief Returns a local integral contribution object for the given element. + //! \param[in] nen Number of nodes on element + //! \param[in] iEl Global element number + //! \param[in] neumann Whether or not we are assembling Neumann BC's + virtual LocalIntegral* getLocalIntegral(size_t nen, size_t iEl, + bool neumann) const; + + //! \brief Initializes current element for numerical integration. + //! \param[in] MNPC Matrix of nodal point correspondance for current element + //! \param[in] nPt Number of volumetric sampling points in this element + //! \param elmInt The local integral object for current element + virtual bool initElement(const std::vector& MNPC, + const Vec3&, size_t nPt, + LocalIntegral& elmInt); + + //! \brief Initializes current element for boundary integration. + //! \param[in] MNPC Matrix of nodal point correspondance for current element + //! \param elmInt The local integral object for current element + virtual bool initElementBou(const std::vector& MNPC, + LocalIntegral& elmInt); + //! \brief Evaluates volumetric sampling point data at an interior point. + //! \param elmInt The local integral object to receive the contributions //! \param[in] fe Finite element data of current point //! \param[in] X Cartesian coordinates of current integration point - virtual bool reducedInt(const FiniteElement& fe, const Vec3& X) const; + virtual bool reducedInt(LocalIntegral& elmInt, const FiniteElement& fe, + const Vec3& X) const; //! \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] X Cartesian coordinates of current integration point - virtual bool evalInt(LocalIntegral*& elmInt, const FiniteElement& fe, + 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] 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, const FiniteElement& fe, + const Vec3& X, const Vec3& normal) const; + //! \brief Returns which integrand to be used. virtual int getIntegrandType() const { return myProblem.getIntegrandType(); } }; diff --git a/Apps/FiniteDefElasticity/NonlinearElasticityTL.C b/Apps/FiniteDefElasticity/NonlinearElasticityTL.C index f040542a..3ed5d70c 100644 --- a/Apps/FiniteDefElasticity/NonlinearElasticityTL.C +++ b/Apps/FiniteDefElasticity/NonlinearElasticityTL.C @@ -23,6 +23,10 @@ #define epsR 1.0e-16 #endif +#ifdef USE_OPENMP +#include +#endif + NonlinearElasticityTL::NonlinearElasticityTL (unsigned short int n, bool axS) : Elasticity(n,axS) @@ -135,7 +139,7 @@ bool NonlinearElasticityTL::evalInt (LocalIntegral& elmInt, Matrix Bmat; Tensor F(nDF); SymmTensor E(nsd,axiSymmetry), S(nsd,axiSymmetry); - if (!this->kinematics(elMat.vec,fe.N,fe.dNdX,X.x,Bmat,F,E)) + if (!this->kinematics(elMat.vec.front(),fe.N,fe.dNdX,X.x,Bmat,F,E)) return false; // Evaluate the constitutive relation @@ -206,7 +210,10 @@ bool NonlinearElasticityTL::evalBou (LocalIntegral& elmInt, Vec3 T = (*tracFld)(X,normal); // Store the traction value for vizualization - if (!T.isZero()) tracVal[X] = T; +#ifdef USE_OPENMP + if (omp_get_max_threads() == 1) +#endif + if (!T.isZero()) tracVal[X] = T; // Check for with-rotated pressure load unsigned short int i, j; @@ -216,7 +223,8 @@ bool NonlinearElasticityTL::evalBou (LocalIntegral& elmInt, Matrix B; Tensor F(nDF); SymmTensor dummy(0); - if (!this->kinematics(elMat.vec,fe.N,fe.dNdX,X.x,B,F,dummy)) return false; + if (!this->kinematics(elMat.vec.front(),fe.N,fe.dNdX,X.x,B,F,dummy)) + return false; // Compute its inverse and determinant, J double J = F.inverse(); @@ -244,12 +252,12 @@ bool NonlinearElasticityTL::evalBou (LocalIntegral& elmInt, } -bool NonlinearElasticityTL::kinematics (const Vectors& eV, +bool NonlinearElasticityTL::kinematics (const Vector& eV, const Vector& N, const Matrix& dNdX, double r, Matrix& Bmat, Tensor& F, SymmTensor& E) const { - if (eV.empty() || eV.front().empty()) + if (eV.empty()) { // Initial state, unit deformation gradient and linear B-matrix F = 1.0; @@ -262,7 +270,7 @@ bool NonlinearElasticityTL::kinematics (const Vectors& eV, const size_t nenod = dNdX.rows(); const size_t nstrc = axiSymmetry ? 4 : nsd*(nsd+1)/2; - if (eV.front().size() != nenod*nsd || dNdX.cols() < nsd) + if (eV.size() != nenod*nsd || dNdX.cols() < nsd) { std::cerr <<" *** NonlinearElasticityTL::kinematics: Invalid dimension," <<" dNdX("<< nenod <<","<< dNdX.cols() <<")"<< std::endl; @@ -273,10 +281,10 @@ bool NonlinearElasticityTL::kinematics (const Vectors& eV, // and the Green-Lagrange strains, E_ij = 0.5(F_ij+F_ji+F_ki*F_kj). // Notice that the matrix multiplication method used here treats the - // element displacement vector, *eV, as a matrix whose number of columns + // element displacement vector, eV, as a matrix whose number of columns // equals the number of rows in the matrix dNdX. Matrix dUdX; - if (!dUdX.multiplyMat(eV.front(),dNdX)) // dUdX = Grad{u} = eV*dNd + if (!dUdX.multiplyMat(eV,dNdX)) // dUdX = Grad{u} = eV*dNd return false; unsigned short int i, j, k; @@ -304,7 +312,7 @@ bool NonlinearElasticityTL::kinematics (const Vectors& eV, // Add the unit tensor to F to form the deformation gradient F += 1.0; // Add the dU/r term to the F(3,3)-term for axisymmetric problems - if (axiSymmetry && r > epsR) F(3,3) += eV.front().dot(N,0,nsd)/r; + if (axiSymmetry && r > epsR) F(3,3) += eV.dot(N,0,nsd)/r; #ifdef INT_DEBUG std::cout <<"NonlinearElasticityTL::F =\n"<< F; diff --git a/Apps/FiniteDefElasticity/NonlinearElasticityTL.h b/Apps/FiniteDefElasticity/NonlinearElasticityTL.h index 673f2734..1056b300 100644 --- a/Apps/FiniteDefElasticity/NonlinearElasticityTL.h +++ b/Apps/FiniteDefElasticity/NonlinearElasticityTL.h @@ -69,7 +69,7 @@ public: protected: //! \brief Calculates some kinematic quantities at current point. - //! \param[in] eV Element solution vectors + //! \param[in] eV Element solution vector //! \param[in] N Basis function values at current point //! \param[in] dNdX Basis function gradients at current point //! \param[in] r Radial coordinate of current point @@ -80,7 +80,7 @@ protected: //! \details The deformation gradient \b F and the nonlinear //! strain-displacement matrix \b B are established. //! The B-matrix is formed only when the variable \a formB is true. - virtual bool kinematics(const Vectors& eV, + virtual bool kinematics(const Vector& eV, const Vector& N, const Matrix& dNdX, double r, Matrix& Bmat, Tensor& F, SymmTensor& E) const; diff --git a/Apps/FiniteDefElasticity/NonlinearElasticityUL.C b/Apps/FiniteDefElasticity/NonlinearElasticityUL.C index 971cba40..f8d0b332 100644 --- a/Apps/FiniteDefElasticity/NonlinearElasticityUL.C +++ b/Apps/FiniteDefElasticity/NonlinearElasticityUL.C @@ -25,6 +25,10 @@ #define epsR 1.0e-16 #endif +#ifdef USE_OPENMP +#include +#endif + NonlinearElasticityUL::NonlinearElasticityUL (unsigned short int n, bool axS, char lop) @@ -46,35 +50,20 @@ void NonlinearElasticityUL::print (std::ostream& os) const void NonlinearElasticityUL::setMode (SIM::SolutionMode mode) { - if (!myMats) return; - - myMats->rhsOnly = false; + m_mode = mode; eM = eKm = eKg = 0; - iS = eS = eV = 0; + eS = iS = 0; switch (mode) { case SIM::STATIC: - myMats->resize(1,1); - eKm = &myMats->A[0]; - eKg = &myMats->A[0]; + case SIM::RHS_ONLY: + eKm = eKg = 1; break; case SIM::DYNAMIC: - myMats->resize(2,1); - eKm = &myMats->A[0]; - eKg = &myMats->A[0]; - eM = &myMats->A[1]; - break; - - case SIM::RHS_ONLY: - if (myMats->A.empty()) - myMats->resize(1,1); - else - myMats->b.resize(1); - eKm = &myMats->A[0]; - eKg = &myMats->A[0]; - myMats->rhsOnly = true; + eKm = eKg = 1; + eM = 2; break; default: @@ -82,15 +71,61 @@ void NonlinearElasticityUL::setMode (SIM::SolutionMode mode) return; } - // We always need the force+displacement vectors in nonlinear simulations - iS = &myMats->b[0]; - eS = &myMats->b[0]; - mySols.resize(1); - eV = &mySols[0]; + // We always need the force vectors in nonlinear simulations + iS = eS = 1; tracVal.clear(); } +LocalIntegral* NonlinearElasticityUL::getLocalIntegral (size_t nen, size_t, + bool neumann) const +{ + ElmMats* result = new ElmMats; + switch (m_mode) + { + case SIM::STATIC: + result->withLHS = !neumann; + case SIM::RHS_ONLY: + result->rhsOnly = neumann; + result->resize(neumann?0:1,1); + break; + + case SIM::DYNAMIC: + result->withLHS = !neumann; + result->rhsOnly = neumann; + result->resize(neumann?0:2,1); + break; + + case SIM::VIBRATION: + case SIM::BUCKLING: + result->withLHS = true; + result->resize(2,0); + break; + + case SIM::STIFF_ONLY: + case SIM::MASS_ONLY: + result->withLHS = true; + result->resize(1,0); + break; + + case SIM::RECOVERY: + result->rhsOnly = true; + break; + + default: + ; + } + + for (size_t i = 0; i < result->A.size(); i++) + result->A[i].resize(nsd*nen,nsd*nen); + + if (result->b.size()) + result->b.front().resize(nsd*nen); + + return result; +} + + void NonlinearElasticityUL::initIntegration (const TimeDomain& prm) { if (material) @@ -108,15 +143,18 @@ void NonlinearElasticityUL::initResultPoints (double lambda) } -bool NonlinearElasticityUL::evalInt (LocalIntegral*& elmInt, +bool NonlinearElasticityUL::evalInt (LocalIntegral& elmInt, const FiniteElement& fe, const TimeDomain& prm, const Vec3& X) const { + ElmMats& elMat = static_cast(elmInt); + // Evaluate the deformation gradient, F, and the Green-Lagrange strains, E + Matrix Bmat, dNdx; Tensor F(nDF); SymmTensor E(nsd,axiSymmetry); - if (!this->kinematics(fe.N,fe.dNdX,X.x,F,E)) + if (!this->kinematics(elMat.vec.front(),fe.N,fe.dNdX,X.x,Bmat,F,E)) return false; // Axi-symmetric integration point volume; 2*pi*r*|J|*w @@ -149,11 +187,11 @@ bool NonlinearElasticityUL::evalInt (LocalIntegral*& elmInt, // Compute the small-deformation strain-displacement matrix B from dNdx if (axiSymmetry) { - r += eV->dot(fe.N,0,nsd); - this->formBmatrix(fe.N,dNdx,r); + r += elMat.vec.front().dot(fe.N,0,nsd); + this->formBmatrix(Bmat,fe.N,dNdx,r); } else - this->formBmatrix(dNdx); + this->formBmatrix(Bmat,dNdx); #ifdef INT_DEBUG std::cout <<"NonlinearElasticityUL::dNdx ="<< dNdx; @@ -165,12 +203,13 @@ bool NonlinearElasticityUL::evalInt (LocalIntegral*& elmInt, { // Initial state, no deformation yet if (axiSymmetry) - this->formBmatrix(fe.N,fe.dNdX,r); + this->formBmatrix(Bmat,fe.N,fe.dNdX,r); else - this->formBmatrix(fe.dNdX); + this->formBmatrix(Bmat,fe.dNdX); } // Evaluate the constitutive relation + Matrix Cmat; SymmTensor sigma(nsd,axiSymmetry); if (eKm || eKg || iS) { @@ -182,35 +221,36 @@ bool NonlinearElasticityUL::evalInt (LocalIntegral*& elmInt, if (eKm) { // Integrate the material stiffness matrix + Matrix CB; CB.multiply(Cmat,Bmat).multiply(detJW); // CB = C*B*|J|*w - eKm->multiply(Bmat,CB,true,false,true); // EK += B^T * CB + elMat.A[eKm-1].multiply(Bmat,CB,true,false,true); // EK += B^T * CB } if (eKg && lHaveStrains) // Integrate the geometric stiffness matrix - this->formKG(*eKg,fe.N,dNdx,r,sigma,detJW); + this->formKG(elMat.A[eKg-1],fe.N,dNdx,r,sigma,detJW); if (eM) // Integrate the mass matrix - this->formMassMatrix(*eM,fe.N,X,detJW); + this->formMassMatrix(elMat.A[eM-1],fe.N,X,detJW); if (iS && lHaveStrains) { // Integrate the internal forces sigma *= -detJW; - if (!Bmat.multiply(sigma,*iS,true,true)) // ES -= B^T*sigma + if (!Bmat.multiply(sigma,elMat.b[iS-1],true,true)) // ES -= B^T*sigma return false; } if (eS) // Integrate the load vector due to gravitation and other body forces - this->formBodyForce(*eS,fe.N,X,detJW); + this->formBodyForce(elMat.b[eS-1],fe.N,X,detJW); - return this->getIntegralResult(elmInt); + return true; } -bool NonlinearElasticityUL::evalBou (LocalIntegral*& elmInt, +bool NonlinearElasticityUL::evalBou (LocalIntegral& elmInt, const FiniteElement& fe, const Vec3& X, const Vec3& normal) const { @@ -231,7 +271,10 @@ bool NonlinearElasticityUL::evalBou (LocalIntegral*& elmInt, Vec3 T = (*tracFld)(X,normal); // Store the traction value for vizualization - if (!T.isZero()) tracVal[X] = T; +#ifdef USE_OPENMP + if (omp_get_max_threads() == 1) +#endif + if (!T.isZero()) tracVal[X] = T; // Axi-symmetric integration point volume; 2*pi*r*|J|*w double detJW = axiSymmetry ? 2.0*M_PI*X.x*fe.detJxW : fe.detJxW; @@ -241,7 +284,8 @@ bool NonlinearElasticityUL::evalBou (LocalIntegral*& elmInt, { // Compute the deformation gradient, F Tensor F(nDF); - if (!this->formDefGradient(fe.N,fe.dNdX,X.x,F)) return false; + if (!this->formDefGradient(elmInt.vec.front(),fe.N,fe.dNdX,X.x,F)) + return false; // Check for with-rotated pressure load if (tracFld->isNormalPressure()) @@ -264,29 +308,31 @@ bool NonlinearElasticityUL::evalBou (LocalIntegral*& elmInt, } // Integrate the force vector - Vector& ES = *eS; + Vector& ES = static_cast(elmInt).b[eS-1]; for (size_t a = 1; a <= fe.N.size(); a++) for (i = 1; i <= nsd; i++) ES(nsd*(a-1)+i) += T[i-1]*fe.N(a)*detJW; - return this->getIntegralResult(elmInt); + return true; } -bool NonlinearElasticityUL::formDefGradient (const Vector& N, +bool NonlinearElasticityUL::formDefGradient (const Vector& eV, const Vector& N, const Matrix& dNdX, double r, Tensor& F) const { - static SymmTensor dummy(0); - return this->kinematics(N,dNdX,r,F,dummy); + Matrix B; + SymmTensor dummy(0); + return this->kinematics(eV,N,dNdX,r,B,F,dummy); } -bool NonlinearElasticityUL::kinematics (const Vector& N, const Matrix& dNdX, - double r, Tensor& F, +bool NonlinearElasticityUL::kinematics (const Vector& eV, + const Vector& N, const Matrix& dNdX, + double r, Matrix&, Tensor& F, SymmTensor& E) const { - if (!eV || eV->empty()) + if (eV.empty()) { // Initial state, unit deformation gradient and zero strains F = 1.0; @@ -295,7 +341,7 @@ bool NonlinearElasticityUL::kinematics (const Vector& N, const Matrix& dNdX, } const size_t nenod = dNdX.rows(); - if (eV->size() != nenod*nsd || dNdX.cols() < nsd) + if (eV.size() != nenod*nsd || dNdX.cols() < nsd) { std::cerr <<" *** NonlinearElasticityUL::kinematics: Invalid dimension," <<" dNdX("<< nenod <<","<< dNdX.cols() <<")"<< std::endl; @@ -309,7 +355,7 @@ bool NonlinearElasticityUL::kinematics (const Vector& N, const Matrix& dNdX, // element displacement vector, *eV, as a matrix whose number of columns // equals the number of rows in the matrix dNdX. Matrix dUdX; - if (!dUdX.multiplyMat(*eV,dNdX)) // dUdX = Grad{u} = eV*dNdX + if (!dUdX.multiplyMat(eV,dNdX)) // dUdX = Grad{u} = eV*dNdX return false; unsigned short int i, j, k; @@ -337,10 +383,10 @@ bool NonlinearElasticityUL::kinematics (const Vector& N, const Matrix& dNdX, // Add the unit tensor to F to form the deformation gradient F += 1.0; // Add the dU/r term to the F(3,3)-term for axisymmetric problems - if (axiSymmetry && r > epsR) F(3,3) += eV->dot(N,0,nsd)/r; + if (axiSymmetry && r > epsR) F(3,3) += eV.dot(N,0,nsd)/r; #ifdef INT_DEBUG - std::cout <<"NonlinearElasticityUL::eV ="<< *eV; + std::cout <<"NonlinearElasticityUL::eV ="<< eV.front(); std::cout <<"NonlinearElasticityUL::F =\n"<< F; #endif @@ -361,7 +407,7 @@ void ElasticityNormUL::initIntegration (const TimeDomain& prm) } -bool ElasticityNormUL::evalInt (LocalIntegral*& elmInt, +bool ElasticityNormUL::evalInt (LocalIntegral& elmInt, const FiniteElement& fe, const TimeDomain& prm, const Vec3& X) const @@ -369,23 +415,24 @@ bool ElasticityNormUL::evalInt (LocalIntegral*& elmInt, NonlinearElasticityUL& ulp = static_cast(myProblem); // Evaluate the deformation gradient, F, and the Green-Lagrange strains, E + Matrix B; Tensor F(ulp.nDF); SymmTensor E(ulp.nDF); - if (!ulp.kinematics(fe.N,fe.dNdX,X.x,F,E)) + if (!ulp.kinematics(elmInt.vec.front(),fe.N,fe.dNdX,X.x,B,F,E)) return false; // Compute the strain energy density, U(E) = Int_E (S:Eps) dEps // and the Cauchy stress tensor, sigma - double U = 0.0; + Matrix Cmat; double U = 0.0; SymmTensor sigma(E.dim(),ulp.isAxiSymmetric()||ulp.material->isPlaneStrain()); - if (!ulp.material->evaluate(ulp.Cmat,sigma,U,X,F,E,3,&prm)) + if (!ulp.material->evaluate(Cmat,sigma,U,X,F,E,3,&prm)) return false; // Axi-symmetric integration point volume; 2*pi*r*|J|*w double detJW = ulp.isAxiSymmetric() ? 2.0*M_PI*X.x*fe.detJxW : fe.detJxW; // Integrate the norms - return evalInt(getElmNormBuffer(elmInt,6),sigma,U,F.det(),detJW); + return evalInt(static_cast(elmInt),sigma,U,F.det(),detJW); } @@ -414,7 +461,7 @@ bool ElasticityNormUL::evalInt (ElmNorm& pnorm, const SymmTensor& S, } -bool ElasticityNormUL::evalBou (LocalIntegral*& elmInt, +bool ElasticityNormUL::evalBou (LocalIntegral& elmInt, const FiniteElement& fe, const Vec3& X, const Vec3& normal) const { @@ -424,7 +471,7 @@ bool ElasticityNormUL::evalBou (LocalIntegral*& elmInt, // Evaluate the current surface traction Vec3 t = ulp.getTraction(X,normal); // Evaluate the current displacement field - Vec3 u = ulp.evalSol(fe.N); + Vec3 u = ulp.evalSol(elmInt.vec.front(),fe.N); // Integrate the external energy (path integral) if (iP < Ux.size()) @@ -445,6 +492,6 @@ bool ElasticityNormUL::evalBou (LocalIntegral*& elmInt, // Axi-symmetric integration point volume; 2*pi*r*|J|*w double detJW = ulp.isAxiSymmetric() ? 2.0*M_PI*X.x*fe.detJxW : fe.detJxW; - getElmNormBuffer(elmInt)[1] += Ux[iP++]*detJW; + static_cast(elmInt)[1] += Ux[iP++]*detJW; return true; } diff --git a/Apps/FiniteDefElasticity/NonlinearElasticityUL.h b/Apps/FiniteDefElasticity/NonlinearElasticityUL.h index 743bd705..df29e325 100644 --- a/Apps/FiniteDefElasticity/NonlinearElasticityUL.h +++ b/Apps/FiniteDefElasticity/NonlinearElasticityUL.h @@ -16,6 +16,8 @@ #include "Elasticity.h" +class ElmNorm; + /*! \brief Class representing the integrand of the nonlinear elasticity problem. @@ -45,6 +47,12 @@ public: //! \param[in] mode The solution mode to use virtual void setMode(SIM::SolutionMode mode); + //! \brief Returns a local integral container for the given element. + //! \param[in] nen Number of DOFs on element + //! \param[in] neumann Whether or not we are assembling Neumann BC's + virtual LocalIntegral* getLocalIntegral(size_t nen, size_t, + bool neumann) const; + //! \brief Initializes the integrand for a new integration loop. //! \param[in] prm Nonlinear solution algorithm parameters virtual void initIntegration(const TimeDomain& prm); @@ -57,7 +65,7 @@ public: //! \param[in] fe Finite element data of current integration point //! \param[in] prm Nonlinear solution algorithm parameters //! \param[in] X Cartesian coordinates of current integration point - virtual bool evalInt(LocalIntegral*& elmInt, const FiniteElement& fe, + virtual bool evalInt(LocalIntegral& elmInt, const FiniteElement& fe, const TimeDomain& prm, const Vec3& X) const; //! \brief Evaluates the integrand at a boundary point. @@ -70,7 +78,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, const FiniteElement& fe, + 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. @@ -80,26 +88,26 @@ public: virtual NormBase* getNormIntegrand(AnaSol* = 0) const; //! \brief Calculates some kinematic quantities at current point. + //! \param[in] eV Element solution vector //! \param[in] N Basis function values at current point //! \param[in] dNdX Basis function gradients at current point //! \param[in] r Radial coordinate of current point //! \param[out] F Deformation gradient at current point + //! \param[out] B The strain-displacement matrix at current point //! \param[out] E Green-Lagrange strain tensor at current point - virtual bool kinematics(const Vector& N, const Matrix& dNdX, double r, - Tensor& F, SymmTensor& E) const; + virtual bool kinematics(const Vector& eV, + const Vector& N, const Matrix& dNdX, double r, + Matrix& B, Tensor& F, SymmTensor& E) const; protected: //! \brief Calculates the deformation gradient at current point. + //! \param[in] eV Element solution vector //! \param[in] N Basis function values at current point //! \param[in] dNdX Basis function gradients at current point //! \param[in] r Radial coordinate of current point //! \param[out] F Deformation gradient at current point - bool formDefGradient(const Vector& N, const Matrix& dNdX, double r, - Tensor& F) const; - -protected: - mutable Matrix dNdx; //!< Basis function gradients in current configuration - mutable Matrix CB; //!< Result of the matrix-matrix product C*B + bool formDefGradient(const Vector& eV, const Vector& N, + const Matrix& dNdX, double r, Tensor& F) const; private: char loadOp; //!< Load option @@ -135,7 +143,7 @@ public: //! \param[in] fe Finite element data of current integration point //! \param[in] prm Nonlinear solution algorithm parameters //! \param[in] X Cartesian coordinates of current integration point - virtual bool evalInt(LocalIntegral*& elmInt, const FiniteElement& fe, + virtual bool evalInt(LocalIntegral& elmInt, const FiniteElement& fe, const TimeDomain& prm, const Vec3& X) const; //! \brief Evaluates the integrand at a boundary point. @@ -143,7 +151,7 @@ public: //! \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, const FiniteElement& fe, + 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/Apps/FiniteDefElasticity/NonlinearElasticityULMX.C b/Apps/FiniteDefElasticity/NonlinearElasticityULMX.C index e85b5400..44de5eb2 100644 --- a/Apps/FiniteDefElasticity/NonlinearElasticityULMX.C +++ b/Apps/FiniteDefElasticity/NonlinearElasticityULMX.C @@ -15,8 +15,10 @@ #include "MaterialBase.h" #include "FiniteElement.h" #include "ElmMats.h" +#include "ElmNorm.h" #include "TimeDomain.h" #include "Utilities.h" +#include "Tensor.h" #include "Vec3Oper.h" #ifdef USE_FTNMAT @@ -42,13 +44,102 @@ extern "C" { #endif +/*! + \brief A struct with integration point data needed by \a finalizeElement. +*/ + +struct ItgPtData +{ + Tensor F; //!< Deformation gradient at current step/iteration + Tensor Fp; //!< Deformation gradient at previous (converged) step + Vector Nr; //!< Basis function values (for axisymmetric problems) + Matrix dNdx; //!< Basis function gradients at current configuration + Vector Phi; //!< Internal modes for the pressure/volumetric-change fields + Vec4 X; //!< Cartesian coordinates of current integration point + double detJW; //!< Jacobian determinant times integration point weight + //! \brief Default constructor. + ItgPtData() : F(3), Fp(3) { detJW = 0.0; } +}; + + +/*! + \brief Class containing internal data for a mixed element. +*/ + +class MxElmData +{ +public: + //! \brief Default constructor. + MxElmData() : Hh(0), iP(0) {} + //! \brief Empty destructor. + virtual ~MxElmData() {} + + //! \brief Initializes the element data. + bool init(const Vec3& Xc, size_t nPt, int p, unsigned short int nsd) + { + iP = 0; + X0 = Xc; + myData.resize(nPt); + if (Hh) + { + size_t nPm = utl::Pascal(p,nsd); + Hh->resize(nPm,nPm,true); + return true; + } + std::cerr <<" *** MxElmData::init: No Hh matrix defined."<< std::endl; + return false; + } + + Vec3 X0; //!< Cartesian coordinates of the element center + Matrix* Hh; //!< Pointer to element Hh-matrix associated with pressure modes + size_t iP; //!< Local integration point counter + + std::vector myData; //!< Local integration point data +}; + + +/*! + \brief Class collecting element matrices associated with a mixed FEM problem. +*/ + +class MxMats : public MxElmData, public ElmMats {}; + + +/*! + \brief Class representing integrated norm quantities for a mixed element. +*/ + +class MxNorm : public MxElmData, public LocalIntegral +{ +public: + //! \brief The constructor initializes the ElmNorm pointer. + MxNorm(LocalIntegral& n) + { + myNorm = static_cast(&n); + Hh = new Matrix(); + } + + //! \brief Alternative constructor allocating an internal ElmNorm. + MxNorm(size_t n) + { + myNorm = new ElmNorm(n); + Hh = new Matrix(); + } + + //! \brief The destructor destroys the ElmNorm object. + virtual ~MxNorm() { myNorm->destruct(); } + + //! \brief Returns the LocalIntegral object to assemble into the global one. + virtual const LocalIntegral* ref() const { return myNorm; } + + ElmNorm* myNorm; //!< Pointer to the actual element norms object +}; + + NonlinearElasticityULMX::NonlinearElasticityULMX (unsigned short int n, bool axS, int pp) - : NonlinearElasticityUL(n,axS) + : NonlinearElasticityUL(n,axS), p(pp) { - p = pp; - iP = 0; - // Both the current and previous solutions are needed primsol.resize(2); } @@ -63,106 +154,90 @@ void NonlinearElasticityULMX::print (std::ostream& os) const } -void NonlinearElasticityULMX::setMode (SIM::SolutionMode mode) +LocalIntegral* NonlinearElasticityULMX::getLocalIntegral (size_t nen, size_t, + bool neumann) const { - if (!myMats) return; + MxMats* result = new MxMats; - myMats->rhsOnly = false; - Hh = eM = eKm = eKg = 0; - iS = eS = eV = 0; - - switch (mode) + switch (m_mode) { case SIM::STATIC: - myMats->resize(2,1); - eKm = &myMats->A[0]; - eKg = &myMats->A[0]; - Hh = &myMats->A[1]; + result->withLHS = !neumann; + case SIM::RHS_ONLY: + result->rhsOnly = neumann; + result->resize(neumann?1:2,1); break; case SIM::DYNAMIC: - myMats->resize(3,1); - eKm = &myMats->A[0]; - eKg = &myMats->A[0]; - eM = &myMats->A[1]; - Hh = &myMats->A[2]; - break; - - case SIM::RHS_ONLY: - if (myMats->A.size() < 2) - myMats->resize(2,1); - else - myMats->b.resize(1); - eKm = &myMats->A[0]; - eKg = &myMats->A[0]; - Hh = &myMats->A.back(); - myMats->rhsOnly = true; + result->withLHS = !neumann; + result->rhsOnly = neumann; + result->resize(neumann?1:3,1); break; case SIM::RECOVERY: - if (myMats->A.empty()) - myMats->A.resize(1); - if (mySols.empty()) - mySols.resize(1); - Hh = &myMats->A.back(); - eV = &mySols[0]; - myMats->rhsOnly = true; - return; + result->rhsOnly = true; + result->resize(1,0); + break; default: - this->Elasticity::setMode(mode); - return; + return result; } - // We always need the force+displacement vectors in nonlinear simulations - iS = &myMats->b[0]; - eS = &myMats->b[0]; - mySols.resize(2); - eV = &mySols[0]; - tracVal.clear(); + result->Hh = &result->A.back(); + + for (size_t i = 1; i < result->A.size(); i++) + result->A[i-1].resize(nsd*nen,nsd*nen); + + if (result->b.size()) + result->b.front().resize(nsd*nen); + + return result; } bool NonlinearElasticityULMX::initElement (const std::vector& MNPC, - const Vec3& Xc, size_t nPt) + const Vec3& Xc, size_t nPt, + LocalIntegral& elmInt) { #if INT_DEBUG > 0 std::cout <<"\n\n *** Entering NonlinearElasticityULMX::initElement"; std::cout <<", Xc = "<< Xc << std::endl; #endif - iP = 0; - X0 = Xc; - myData.resize(nPt); - size_t nPm = utl::Pascal(p,nsd); - if (Hh) Hh->resize(nPm,nPm,true); - - // The other element matrices are initialized by the parent class method - return this->NonlinearElasticityUL::initElement(MNPC); + return static_cast(elmInt).init(Xc,nPt,p,nsd) && + this->IntegrandBase::initElement(MNPC,elmInt); } -bool NonlinearElasticityULMX::evalInt (LocalIntegral*& elmInt, +bool NonlinearElasticityULMX::evalInt (LocalIntegral& elmInt, const FiniteElement& fe, const TimeDomain&, const Vec3& X) const { - if (mySols.size() < 2) return false; + MxElmData* mx = 0; + MxNorm* mxNrm = 0; + MxMats* mxMat = dynamic_cast(&elmInt); + if (mxMat) + mx = mxMat; + else if ((mxNrm = dynamic_cast(&elmInt))) + mx = mxNrm; + else + return false; + + ItgPtData& ptData = mx->myData[mx->iP++]; + const Vectors& eV = mxNrm ? mxNrm->myNorm->vec : elmInt.vec; + if (eV.size() < 2) return false; - ItgPtData& ptData = myData[iP++]; - Vectors& eVs = const_cast(mySols); #if INT_DEBUG > 1 std::cout <<"NonlinearElasticityULMX::dNdX ="<< fe.dNdX; #endif // Evaluate the deformation gradient, Fp, at previous configuration - const_cast(this)->eV = &eVs[1]; - if (!this->formDefGradient(fe.N,fe.dNdX,X.x,ptData.Fp)) + if (!this->formDefGradient(eV[1],fe.N,fe.dNdX,X.x,ptData.Fp)) return false; // Evaluate the deformation gradient, F, at current configuration - const_cast(this)->eV = &eVs[0]; - if (!this->formDefGradient(fe.N,fe.dNdX,X.x,ptData.F)) + if (!this->formDefGradient(eV[0],fe.N,fe.dNdX,X.x,ptData.F)) return false; if (nDF == 2) // In 2D we always assume plane strain so set F(3,3)=1 @@ -184,13 +259,13 @@ bool NonlinearElasticityULMX::evalInt (LocalIntegral*& elmInt, // Push-forward the basis function gradients to current configuration ptData.dNdx.multiply(fe.dNdX,Fi); // dNdx = dNdX * F^-1 if (axiSymmetry) - ptData.Nr = fe.N * (1.0/(X.x + eV->dot(fe.N,0,nsd))); + ptData.Nr = fe.N * (1.0/(X.x + eV.front().dot(fe.N,0,nsd))); ptData.X.assign(X); ptData.detJW = axiSymmetry ? 2.0*M_PI*X.x*fe.detJxW : fe.detJxW; // Evaluate the pressure modes (generalized coordinates) - Vec3 Xg = X-X0; + Vec3 Xg = X - mx->X0; if (nsd == 3) utl::Pascal(p,Xg.x,Xg.y,Xg.z,ptData.Phi); else if (nsd == 2) @@ -198,33 +273,36 @@ bool NonlinearElasticityULMX::evalInt (LocalIntegral*& elmInt, else return false; - if (Hh) + if (mx->Hh) // Integrate the Hh-matrix - Hh->outer_product(ptData.Phi,ptData.Phi*ptData.detJW,true); + mx->Hh->outer_product(ptData.Phi,ptData.Phi*ptData.detJW,true); if (eM) // Integrate the mass matrix - this->formMassMatrix(*eM,fe.N,X,J*ptData.detJW); + this->formMassMatrix(mxMat->A[eM-1],fe.N,X,J*ptData.detJW); if (eS) // Integrate the load vector due to gravitation and other body forces - this->formBodyForce(*eS,fe.N,X,J*ptData.detJW); + this->formBodyForce(mxMat->b[eS-1],fe.N,X,J*ptData.detJW); - return this->getIntegralResult(elmInt); + return true; } -bool NonlinearElasticityULMX::finalizeElement (LocalIntegral*& elmInt, +bool NonlinearElasticityULMX::finalizeElement (LocalIntegral& elmInt, const TimeDomain& prm) { if (!iS && !eKm && !eKg) return true; - if (!Hh) return false; + MxMats& mx = static_cast(elmInt); + if (!mx.Hh) return false; + + size_t iP; #if INT_DEBUG > 0 std::cout <<"\n\n *** Entering NonlinearElasticityULMX::finalizeElement\n"; - for (iP = 1; iP <= myData.size(); iP++) + for (iP = 1; iP <= mx.myData.size(); iP++) { - const ItgPtData& pt = myData[iP-1]; + const ItgPtData& pt = mx.myData[iP-1]; std::cout <<"\n X" << iP <<" = "<< pt.X; std::cout <<"\n detJ"<< iP <<"W = "<< pt.detJW; std::cout <<"\n F" << iP <<" =\n"<< pt.F; @@ -237,14 +315,14 @@ bool NonlinearElasticityULMX::finalizeElement (LocalIntegral*& elmInt, // 1. Eliminate the internal pressure DOFs by static condensation. size_t i, j; - size_t nPM = Hh->rows(); - size_t nEN = myData.front().dNdx.rows(); - size_t nGP = myData.size(); + size_t nPM = mx.Hh->rows(); + size_t nEN = mx.myData.front().dNdx.rows(); + size_t nGP = mx.myData.size(); // Invert the H-matrix - if (!utl::invert(*Hh)) return false; + if (!utl::invert(*mx.Hh)) return false; - const Matrix& Hi = *Hh; + const Matrix& Hi = *mx.Hh; Matrix Theta(2,nGP), dNdxBar[nGP]; @@ -258,7 +336,7 @@ bool NonlinearElasticityULMX::finalizeElement (LocalIntegral*& elmInt, for (iP = 0; iP < nGP; iP++) { - const ItgPtData& pt = myData[iP]; + const ItgPtData& pt = mx.myData[iP]; double dVol = pt.F.det()*pt.detJW; h1 += dVol; h2 += pt.Fp.det()*pt.detJW; @@ -292,7 +370,7 @@ bool NonlinearElasticityULMX::finalizeElement (LocalIntegral*& elmInt, // Integrate the Ji- and G-matrices for (iP = 0; iP < nGP; iP++) { - const ItgPtData& pt = myData[iP]; + const ItgPtData& pt = mx.myData[iP]; for (i = 0; i < nPM; i++) { double h1 = pt.Phi[i] * pt.F.det() * pt.detJW; @@ -332,7 +410,7 @@ bool NonlinearElasticityULMX::finalizeElement (LocalIntegral*& elmInt, for (iP = 0; iP < nGP; iP++) { i = iP+1; - ItgPtData& pt = myData[iP]; + ItgPtData& pt = mx.myData[iP]; // Calculate Theta = Hj*Phi Theta(1,i) = Hj1.dot(pt.Phi); @@ -349,7 +427,7 @@ bool NonlinearElasticityULMX::finalizeElement (LocalIntegral*& elmInt, Vector detF(nGP); for (iP = 1; iP <= nGP; iP++) { - ItgPtData& pt = myData[iP-1]; + ItgPtData& pt = mx.myData[iP-1]; detF(iP) = pt.F.det(); pt.F *= pow(fabs(Theta(1,iP)/detF(iP)),1.0/3.0); pt.Fp *= pow(fabs(Theta(2,iP)/pt.Fp.det()),1.0/3.0); @@ -364,7 +442,7 @@ bool NonlinearElasticityULMX::finalizeElement (LocalIntegral*& elmInt, // 2. Evaluate the constitutive relation and calculate mixed pressure state. - Matrix D[nGP]; + Matrix Cmat, D[nGP]; Vector Sigm(nPM), Hsig; std::vector Sig; Sig.reserve(nGP); @@ -375,7 +453,7 @@ bool NonlinearElasticityULMX::finalizeElement (LocalIntegral*& elmInt, #if INT_DEBUG > 0 std::cout <<"\n iGP = "<< iP+1 << std::endl; #endif - const ItgPtData& pt = myData[iP]; + const ItgPtData& pt = mx.myData[iP]; // Evaluate the constitutive relation double U = 0.0; @@ -405,7 +483,7 @@ bool NonlinearElasticityULMX::finalizeElement (LocalIntegral*& elmInt, // Divide pressure by reference volume; Hsig = Hi*Sigm if (lHaveStress) - if (!Hh->multiply(Sigm,Hsig)) return false; + if (!mx.Hh->multiply(Sigm,Hsig)) return false; // 3. Evaluate tangent matrices. @@ -413,7 +491,7 @@ bool NonlinearElasticityULMX::finalizeElement (LocalIntegral*& elmInt, SymmTensor Sigma(nsd,axiSymmetry); for (iP = 0; iP < nGP; iP++) { - const ItgPtData& pt = myData[iP]; + const ItgPtData& pt = mx.myData[iP]; double dFmix = Theta(1,iP+1); if (lHaveStress) { @@ -425,7 +503,7 @@ bool NonlinearElasticityULMX::finalizeElement (LocalIntegral*& elmInt, #endif if (eKg) // Integrate the geometric stiffness matrix - this->formKG(*eKg,pt.Nr,pt.dNdx,1.0,Sigma,dFmix*pt.detJW); + this->formKG(mx.A[eKg-1],pt.Nr,pt.dNdx,1.0,Sigma,dFmix*pt.detJW); } if (eKm) @@ -437,31 +515,33 @@ bool NonlinearElasticityULMX::finalizeElement (LocalIntegral*& elmInt, D[iP] *= dFmix*pt.detJW; if (nsd == 2) acckmx2d_(axiSymmetry,nEN,pt.Nr.ptr(),pt.dNdx.ptr(),dNdxBar[iP].ptr(), - D[iP].ptr(),eKm->ptr()); + D[iP].ptr(),mx.A[eKm-1].ptr()); else - acckmx3d_(nEN,pt.dNdx.ptr(),dNdxBar[iP].ptr(),D[iP].ptr(),eKm->ptr()); + acckmx3d_(nEN,pt.dNdx.ptr(),dNdxBar[iP].ptr(), + D[iP].ptr(),mx.A[eKm-1].ptr()); #endif } if (iS && lHaveStress) { // Compute the small-deformation strain-displacement matrix B from dNdx + Matrix Bmat; if (axiSymmetry) - this->formBmatrix(pt.Nr,pt.dNdx,1.0); + this->formBmatrix(Bmat,pt.Nr,pt.dNdx,1.0); else - this->formBmatrix(pt.dNdx); + this->formBmatrix(Bmat,pt.dNdx); // Integrate the internal forces Sigma *= -dFmix*pt.detJW; - if (!Bmat.multiply(Sigma,*iS,true,true)) // ES -= B^T*sigma + if (!Bmat.multiply(Sigma,mx.b[iS-1],true,true)) // ES -= B^T*sigma return false; } } #if INT_DEBUG > 0 std::cout <<"\n *** Leaving NonlinearElasticityULMX::finalizeElement"; - if (eKm) std::cout <<", eKm ="<< *eKm; - if (iS) std::cout <<" iS ="<< *iS; + if (eKm) std::cout <<", eKm ="<< mx.A[eKm-1]; + if (iS) std::cout <<" iS ="<< mx.b[iS-1]; std::cout << std::endl; #endif return true; @@ -474,59 +554,86 @@ NormBase* NonlinearElasticityULMX::getNormIntegrand (AnaSol*) const } -bool ElasticityNormULMX::initElement (const std::vector& MNPC, - const Vec3& Xc, size_t nPt) +LocalIntegral* ElasticityNormULMX::getLocalIntegral (size_t nen, size_t iEl, + bool neumann) const { - NonlinearElasticityULMX& p = static_cast(myProblem); + LocalIntegral* result = this->NormBase::getLocalIntegral(nen,iEl,neumann); + if (result) return new MxNorm(static_cast(*result)); - p.iP = 0; - p.X0 = Xc; - p.myData.resize(nPt); - size_t nPm = utl::Pascal(p.p,p.nsd); - if (p.Hh) p.Hh->resize(nPm,nPm,true); - - // The other element matrices are initialized by the parent class method - return p.NonlinearElasticityUL::initElement(MNPC); + // Element norms are not requested, so allocate an inter object instead that + // will delete itself when invoking the destruct method. + return new MxNorm(this->getNoFields()); } -bool ElasticityNormULMX::evalInt (LocalIntegral*& elmInt, +bool ElasticityNormULMX::initElement (const std::vector& MNPC, + const Vec3& Xc, size_t nPt, + LocalIntegral& elmInt) +{ + NonlinearElasticityULMX& p = static_cast(myProblem); + MxNorm& mx = static_cast(elmInt); + + return mx.init(Xc,nPt,p.p,p.nsd) && + this->NormBase::initElement(MNPC,*mx.myNorm); +} + + +bool ElasticityNormULMX::initElementBou (const std::vector& MNPC, + LocalIntegral& elmInt) +{ + return myProblem.initElementBou(MNPC,*static_cast(elmInt).myNorm); +} + + +bool ElasticityNormULMX::evalInt (LocalIntegral& elmInt, const FiniteElement& fe, const TimeDomain& td, const Vec3& X) const { - return static_cast(myProblem).evalInt(elmInt, - fe,td,X); + NonlinearElasticityULMX& p = static_cast(myProblem); + + return p.evalInt(elmInt,fe,td,X); } -bool ElasticityNormULMX::finalizeElement (LocalIntegral*& elmInt, +bool ElasticityNormULMX::evalBou (LocalIntegral& elmInt, + const FiniteElement& fe, + const Vec3& X, const Vec3& normal) const +{ + return this->ElasticityNormUL::evalBou(*static_cast(elmInt).myNorm, + fe,X,normal); +} + + +bool ElasticityNormULMX::finalizeElement (LocalIntegral& elmInt, const TimeDomain& prm) { NonlinearElasticityULMX& p = static_cast(myProblem); - if (!p.Hh) return false; + MxNorm& mx = static_cast(elmInt); + + if (!mx.Hh) return false; size_t iP; #if INT_DEBUG > 0 std::cout <<"\n\n *** Entering ElasticityNormULMX::finalizeElement\n"; - for (iP = 1; iP <= p.myData.size(); iP++) + for (iP = 1; iP <= mx.myData.size(); iP++) { - const NonlinearElasticityULMX::ItgPtData& pt = p.myData[iP-1]; + const NonlinearElasticityULMX::ItgPtData& pt = mx.myData[iP-1]; std::cout <<"\n X" << iP <<" = "<< pt.X; std::cout <<"\n detJ"<< iP <<"W = "<< pt.detJW; std::cout <<"\n F" << iP <<" =\n"<< pt.F; std::cout <<" dNdx"<< iP <<" ="<< pt.dNdx; std::cout <<" Phi" << iP <<" ="<< pt.Phi; } - std::cout <<"\n H ="<< *(p.Hh) << std::endl; + std::cout <<"\n H ="<< *(mx.Hh) << std::endl; #endif // 1. Eliminate the internal pressure DOFs by static condensation. - size_t nPM = p.Hh->rows(); - size_t nGP = p.myData.size(); + size_t nPM = mx.Hh->rows(); + size_t nGP = mx.myData.size(); // Invert the H-matrix - Matrix& Hi = *p.Hh; + Matrix& Hi = *mx.Hh; if (!utl::invert(Hi)) return false; Vector Theta(nGP); @@ -537,7 +644,7 @@ bool ElasticityNormULMX::finalizeElement (LocalIntegral*& elmInt, double h1 = 0.0; for (iP = 0; iP < nGP; iP++) - h1 += p.myData[iP].F.det() * p.myData[iP].detJW; + h1 += mx.myData[iP].F.det() * mx.myData[iP].detJW; // All gauss point values are identical Theta.front() = h1 * Hi(1,1); @@ -551,7 +658,7 @@ bool ElasticityNormULMX::finalizeElement (LocalIntegral*& elmInt, // Integrate the Ji-matrix for (iP = 0; iP < nGP; iP++) { - const NonlinearElasticityULMX::ItgPtData& pt = p.myData[iP]; + const ItgPtData& pt = mx.myData[iP]; for (size_t i = 0; i < nPM; i++) Ji[i] += pt.Phi[i] * pt.F.det() * pt.detJW; } @@ -566,7 +673,7 @@ bool ElasticityNormULMX::finalizeElement (LocalIntegral*& elmInt, // Calculate Theta = Hj*Phi for (iP = 0; iP < nGP; iP++) - Theta[iP] = Hj.dot(p.myData[iP].Phi); + Theta[iP] = Hj.dot(mx.myData[iP].Phi); } #if INT_DEBUG > 1 @@ -575,15 +682,14 @@ bool ElasticityNormULMX::finalizeElement (LocalIntegral*& elmInt, // 2. Evaluate the constitutive relation and integrate energy. - ElmNorm& pnorm = NormBase::getElmNormBuffer(elmInt,6); - + Matrix Cmat; SymmTensor Sig(3), Eps(3); for (iP = 0; iP < nGP; iP++) { #if INT_DEBUG > 0 std::cout <<"\n iGP = "<< iP+1 << std::endl; #endif - NonlinearElasticityULMX::ItgPtData& pt = p.myData[iP]; + ItgPtData& pt = mx.myData[iP]; Eps.rightCauchyGreen(pt.F); // Green Lagrange strain tensor Eps -= 1.0; Eps *= 0.5; @@ -595,11 +701,11 @@ bool ElasticityNormULMX::finalizeElement (LocalIntegral*& elmInt, // Compute the strain energy density, U(Eps) = Int_Eps (S:E) dE // and the Cauchy stress tensor, Sig double U = 0.0; - if (!p.material->evaluate(p.Cmat,Sig,U,pt.X,Fbar,Eps,3,&prm,&pt.F)) + if (!p.material->evaluate(Cmat,Sig,U,pt.X,Fbar,Eps,3,&prm,&pt.F)) return false; // Integrate the norms - if (!ElasticityNormUL::evalInt(pnorm,Sig,U,Theta[iP],pt.detJW)) + if (!ElasticityNormUL::evalInt(*mx.myNorm,Sig,U,Theta[iP],pt.detJW)) return false; } diff --git a/Apps/FiniteDefElasticity/NonlinearElasticityULMX.h b/Apps/FiniteDefElasticity/NonlinearElasticityULMX.h index d267b660..700fd4b3 100644 --- a/Apps/FiniteDefElasticity/NonlinearElasticityULMX.h +++ b/Apps/FiniteDefElasticity/NonlinearElasticityULMX.h @@ -15,7 +15,6 @@ #define _NONLINEAR_ELASTICITY_UL_MX_H #include "NonlinearElasticityUL.h" -#include "Tensor.h" /*! @@ -26,20 +25,6 @@ class NonlinearElasticityULMX : public NonlinearElasticityUL { - //! \brief A struct with integration point data needed by \a finalizeElement. - struct ItgPtData - { - Tensor F; //!< Deformation gradient at current step/iteration - Tensor Fp; //!< Deformation gradient at previous (converged) step - Vector Nr; //!< Basis function values (for axisymmetric problems) - Matrix dNdx; //!< Basis function gradients at current configuration - Vector Phi; //!< Internal modes for the pressure/volumetric-change fields - Vec4 X; //!< Cartesian coordinates of current integration point - double detJW; //!< Jacobian determinant times integration point weight - //! \brief Default constructor. - ItgPtData() : F(3), Fp(3) { detJW = 0.0; } - }; - public: //! \brief The default constructor invokes the parent class constructor. //! \param[in] n Number of spatial dimensions @@ -53,16 +38,20 @@ public: //! \brief Prints out problem definition to the given output stream. virtual void print(std::ostream& os) const; - //! \brief Defines the solution mode before the element assembly is started. - //! \param[in] mode The solution mode to use - virtual void setMode(SIM::SolutionMode mode); + //! \brief Returns a local integral container for the given element. + //! \param[in] nen Number of nodes on element + //! \param[in] neumann Whether or not we are assembling Neumann BC's + virtual LocalIntegral* getLocalIntegral(size_t nen, size_t, + bool neumann) const; //! \brief Initializes current element for numerical integration. //! \param[in] MNPC Matrix of nodal point correspondance for current element //! \param[in] Xc Cartesian coordinates of the element center //! \param[in] nPt Number of integration points in this element + //! \param elmInt The local integral object for current element virtual bool initElement(const std::vector& MNPC, - const Vec3& Xc, size_t nPt); + const Vec3& Xc, size_t nPt, + LocalIntegral& elmInt); //! \brief Evaluates the integrand at an interior point. //! \param elmInt The local integral object to receive the contributions @@ -74,17 +63,17 @@ public: //! 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 FiniteElement& fe, + 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 //! loops within an element, required by the mixed formulation with internal //! modes. All data needed during the integration has been stored internally - //! in the \a myData member by the \a evalInt method. + //! in the \a elmInt object by the \a evalInt method. //! \param elmInt The local integral object to receive the contributions //! \param[in] prm Nonlinear solution algorithm parameters - virtual bool finalizeElement(LocalIntegral*& elmInt, const TimeDomain& prm); + virtual bool finalizeElement(LocalIntegral& elmInt, const TimeDomain& prm); //! \brief Returns a pointer to an Integrand for solution norm evaluation. //! \note The Integrand object is allocated dynamically and has to be deleted @@ -96,13 +85,7 @@ public: virtual int getIntegrandType() const { return 4; } private: - int p; //!< Polynomial order of the internal pressure field - Vec3 X0; //!< Cartesian coordinates of the element center - - Matrix* Hh; //!< Pointer to element Hh-matrix associated with pressure modes - - mutable size_t iP; //!< Local integration point counter - mutable std::vector myData; //!< Local integration point data + int p; //!< Polynomial order of the internal pressure field friend class ElasticityNormULMX; }; @@ -121,12 +104,27 @@ public: //! \brief Empty destructor. virtual ~ElasticityNormULMX() {} + //! \brief Returns a local integral contribution object for the given element. + //! \param[in] nen Number of nodes on element + //! \param[in] iEl The element number + //! \param[in] neumann Whether or not we are assembling Neumann BC's + virtual LocalIntegral* getLocalIntegral(size_t nen, size_t iEl, + bool neumann) const; + //! \brief Initializes current element for numerical integration. //! \param[in] MNPC Matrix of nodal point correspondance for current element //! \param[in] Xc Cartesian coordinates of the element center //! \param[in] nPt Number of integration points in this element + //! \param elmInt The local integral object for current element virtual bool initElement(const std::vector& MNPC, - const Vec3& Xc, size_t nPt); + const Vec3& Xc, size_t nPt, + LocalIntegral& elmInt); + + //! \brief Initializes current element for boundary integration. + //! \param[in] MNPC Matrix of nodal point correspondance for current element + //! \param elmInt The local integral object for current element + virtual bool initElementBou(const std::vector& MNPC, + LocalIntegral& elmInt); //! \brief Evaluates the integrand at an interior point. //! \param elmInt The local integral object to receive the contributions @@ -138,16 +136,24 @@ public: //! 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 FiniteElement& fe, + 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] 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, const FiniteElement& fe, + const Vec3& X, const Vec3& normal) const; + //! \brief Finalizes the element norms after the numerical integration. //! \details This method is used to implement the multiple integration point //! loops within an element required by the mixed formulation. It performs //! a subset of the tasks done by NonlinearElasticityULMX::finalizeElement. //! \param elmInt The local integral object to receive the contributions //! \param[in] prm Nonlinear solution algorithm parameters - virtual bool finalizeElement(LocalIntegral*& elmInt, const TimeDomain& prm); + virtual bool finalizeElement(LocalIntegral& elmInt, const TimeDomain& prm); //! \brief Returns which integrand to be used. virtual int getIntegrandType() const { return 4; } diff --git a/Apps/FiniteDefElasticity/NonlinearElasticityULMixed.C b/Apps/FiniteDefElasticity/NonlinearElasticityULMixed.C index df2ef46a..2988b5f5 100644 --- a/Apps/FiniteDefElasticity/NonlinearElasticityULMixed.C +++ b/Apps/FiniteDefElasticity/NonlinearElasticityULMixed.C @@ -15,7 +15,7 @@ #include "MaterialBase.h" #include "FiniteElement.h" #include "TimeDomain.h" -#include "ElmMats.h" +#include "ElmNorm.h" #include "Utilities.h" #include "Vec3Oper.h" @@ -35,6 +35,11 @@ extern "C" { void acckm3d_(const int& nEN, const double* dNdx, const double* D, double* eKt); } +#ifdef USE_OPENMP +#if INT_DEBUG > 0 +#undef INT_DEBUG +#endif +#endif #ifndef INT_DEBUG #define INT_DEBUG 0 #endif @@ -195,18 +200,10 @@ const Vector& NonlinearElasticityULMixed::MixedElmMats::getRHSVector () const NonlinearElasticityULMixed::NonlinearElasticityULMixed (unsigned short int n, bool axS) - : NonlinearElasticityUL(n,axS), Fbar(3), Dmat(7,7) + : NonlinearElasticityUL(n,axS) { - if (myMats) delete myMats; - myMats = new MixedElmMats(); - mySols.resize(NSOL); - - eKm = &myMats->A[Kuu]; - eKg = &myMats->A[Kuu]; - - iS = &myMats->b[Ru]; - eS = &myMats->b[Ru]; - eV = &mySols[U]; + eKm = eKg = Kuu+1; + iS = eS = Ru+1; } @@ -221,18 +218,14 @@ void NonlinearElasticityULMixed::print (std::ostream& os) const void NonlinearElasticityULMixed::setMode (SIM::SolutionMode mode) { + m_mode = mode; switch (mode) { case SIM::INIT: case SIM::STATIC: - tracVal.clear(); - myMats->rhsOnly = false; - break; - case SIM::RHS_ONLY: tracVal.clear(); case SIM::RECOVERY: - myMats->rhsOnly = true; break; default: @@ -242,100 +235,111 @@ void NonlinearElasticityULMixed::setMode (SIM::SolutionMode mode) } +LocalIntegral* NonlinearElasticityULMixed::getLocalIntegral (size_t nen1, + size_t nen2, + size_t, + bool neumann) const +{ + const size_t nedof1 = nsd*nen1; + const size_t nedof = nedof1 + 2*nen2; + + ElmMats* result = new MixedElmMats; + + result->rhsOnly = neumann || m_mode == SIM::RECOVERY; + result->b[Rres].resize(nedof); + result->b[Ru].resize(nedof1); + if (!neumann) + { + result->A[Kuu].resize(nedof1,nedof1); + result->A[Kut].resize(nedof1,nen2); + result->A[Kup].resize(nedof1,nen2); + result->A[Ktt].resize(nen2,nen2); + result->A[Ktp].resize(nen2,nen2); + result->A[Ktan].resize(nedof,nedof); + result->b[Rt].resize(nen2); + result->b[Rp].resize(nen2); + } + if (m_mode == SIM::STATIC) + result->withLHS = !neumann; + + return result; +} + + +#ifndef USE_OPENMP static int iP = 0; //!< Local integration point counter for debug output std::vector mixedDbgEl; //!< List of elements for additional output +#endif bool NonlinearElasticityULMixed::initElement (const std::vector& MNPC1, const std::vector& MNPC2, - size_t n1) + size_t n1, LocalIntegral& elmInt) { +#ifndef USE_OPENMP iP = 0; +#endif + if (primsol.front().empty()) return true; - const size_t nen1 = MNPC1.size(); - const size_t nen2 = MNPC2.size(); - const size_t nedof = nsd*nen1 + 2*nen2; + // Extract the element level solution vectors + elmInt.vec.resize(NSOL); + int ierr = utl::gather(MNPC1,nsd,primsol.front(),elmInt.vec[U]) + + utl::gather(MNPC2,0,2,primsol.front(),elmInt.vec[T],nsd*n1,n1) + + utl::gather(MNPC2,1,2,primsol.front(),elmInt.vec[P],nsd*n1,n1); + if (ierr == 0) return true; - myMats->A[Kut].resize(nsd*nen1,nen2,true); - myMats->A[Kup].resize(nsd*nen1,nen2,true); - myMats->A[Ktt].resize(nen2,nen2,true); - myMats->A[Ktp].resize(nen2,nen2,true); - myMats->A[Ktan].resize(nedof,nedof,true); - - myMats->b[Rt].resize(nen2,true); - myMats->b[Rp].resize(nen2,true); - myMats->b[Rres].resize(nedof,true); - - // Extract the element level volumetric change and pressure vectors - int ierr = 0; - if (!primsol.front().empty()) - if ((ierr = utl::gather(MNPC2,0,2,primsol.front(),mySols[T],nsd*n1,n1) + - utl::gather(MNPC2,1,2,primsol.front(),mySols[P],nsd*n1,n1))) - std::cerr <<" *** NonlinearElasticityULMixed::initElement: Detected " - << ierr/2 <<" node numbers out of range."<< std::endl; - - // The other element matrices are initialized by the parent class method - return this->NonlinearElasticityUL::initElement(MNPC1) && ierr == 0; + std::cerr <<" *** NonlinearElasticityULMixed::initElement: Detected " + << ierr/3 <<" node numbers out of range."<< std::endl; + return false; } bool NonlinearElasticityULMixed::initElementBou (const std::vector& MNPC1, - const std::vector& MNPC2, - size_t) + const std::vector&, + size_t, LocalIntegral& elmInt) { - const size_t nen1 = MNPC1.size(); - const size_t nen2 = MNPC2.size(); - const size_t nedof = nsd*nen1 + 2*nen2; - - myMats->b[Ru].resize(nsd*nen1,true); - myMats->b[Rres].resize(nedof,true); - - // Extract the element level displacement vector - int ierr = 0; - if (!primsol.front().empty()) - if ((ierr = utl::gather(MNPC1,nsd,primsol.front(),mySols[U]))) - std::cerr <<" *** NonlinearElasticityULMixed::initElementBou: Detected " - << ierr <<" node numbers out of range."<< std::endl; - - myMats->withLHS = false; - return ierr == 0; + return this->IntegrandBase::initElementBou(MNPC1,elmInt); } -bool NonlinearElasticityULMixed::evalIntMx (LocalIntegral*& elmInt, +bool NonlinearElasticityULMixed::evalIntMx (LocalIntegral& elmInt, const MxFiniteElement& fe, const TimeDomain& prm, const Vec3& X) const { + ElmMats& elMat = static_cast(elmInt); + #if INT_DEBUG > 0 std::cout <<"\n\n *** Entering NonlinearElasticityULMixed::evalIntMx: iP = " << ++iP <<", (X,t) = "<< X << std::endl; #endif // Evaluate the deformation gradient, F, and the Green-Lagrange strains, E + Matrix B; Tensor F(nDF); SymmTensor E(nsd,axiSymmetry); - if (!this->kinematics(fe.N1,fe.dN1dX,X.x,F,E)) + if (!this->kinematics(elmInt.vec[U],fe.N1,fe.dN1dX,X.x,B,F,E)) return false; bool lHaveStrains = !E.isZero(1.0e-16); // Evaluate the volumetric change and pressure fields - double Theta = fe.N2.dot(mySols[T]) + 1.0; - double Press = fe.N2.dot(mySols[P]); + double Theta = fe.N2.dot(elmInt.vec[T]) + 1.0; + double Press = fe.N2.dot(elmInt.vec[P]); #if INT_DEBUG > 0 - std::cout <<"NonlinearElasticityULMixed::b_theta ="<< mySols[T]; - std::cout <<"NonlinearElasticityULMixed::b_p ="<< mySols[P]; + std::cout <<"NonlinearElasticityULMixed::b_theta ="<< elmInt.vec[T]; + std::cout <<"NonlinearElasticityULMixed::b_p ="<< elmInt.vec[P]; std::cout <<"NonlinearElasticityULMixed::Theta = "<< Theta <<", Press = "<< Press << std::endl; #endif // Axi-symmetric integration point volume; 2*pi*r*|J|*w double detJW = axiSymmetry ? 2.0*M_PI*X.x*fe.detJxW : fe.detJxW; - double r = axiSymmetry ? X.x + eV->dot(fe.N1,0,nsd) : 0.0; + double r = axiSymmetry ? X.x + elmInt.vec[U].dot(fe.N1,0,nsd) : 0.0; // Compute the mixed model deformation gradient, F_bar - Fbar = F; // notice that F_bar always has dimension 3 + Tensor Fbar(3); // notice that F_bar always has dimension 3 + Fbar = F; double r1 = pow(fabs(Theta/F.det()),1.0/3.0); Fbar *= r1; @@ -356,6 +360,7 @@ bool NonlinearElasticityULMixed::evalIntMx (LocalIntegral*& elmInt, if (J == 0.0) return false; // Push-forward the basis function gradients to current configuration + Matrix dNdx; dNdx.multiply(fe.dN1dX,Fi); // dNdx = dNdX * F^-1 // Compute the mixed integration point volume @@ -373,18 +378,20 @@ bool NonlinearElasticityULMixed::evalIntMx (LocalIntegral*& elmInt, if (eM) // Integrate the mass matrix - this->formMassMatrix(*eM,fe.N1,X,J*detJW); + this->formMassMatrix(elMat.A[eM-1],fe.N1,X,J*detJW); if (eS) // Integrate the load vector due to gravitation and other body forces - this->formBodyForce(*eS,fe.N1,X,J*detJW); + this->formBodyForce(elMat.b[eS-1],fe.N1,X,J*detJW); // Evaluate the constitutive relation + Matrix Cmat; SymmTensor Sig(3), Sigma(nsd,axiSymmetry); double U, Bpres = 0.0, Mpres = 0.0; if (!material->evaluate(Cmat,Sig,U,X,Fbar,E,lHaveStrains,&prm)) return false; + Matrix Dmat(7,7); #ifdef USE_FTNMAT // Project the constitutive matrix for the mixed model pcst3d_(Cmat.ptr(),Dmat.ptr(),INT_DEBUG,6); @@ -397,6 +404,7 @@ bool NonlinearElasticityULMixed::evalIntMx (LocalIntegral*& elmInt, if (lHaveStrains) { +#ifndef USE_OPENMP if (prm.it == 0 && find(mixedDbgEl.begin(),mixedDbgEl.end(),fe.iel) != mixedDbgEl.end()) { @@ -413,6 +421,7 @@ bool NonlinearElasticityULMixed::evalIntMx (LocalIntegral*& elmInt, std::cout <<" "<< sigma[i]; std::cout <<" "<< Sig.vonMises() << std::endl; } +#endif // Compute mixed strees Bpres = Sig.trace()/3.0; @@ -420,7 +429,7 @@ bool NonlinearElasticityULMixed::evalIntMx (LocalIntegral*& elmInt, Sigma = Sig + (Mpres-Bpres); // Integrate the geometric stiffness matrix - this->formKG(*eKg,fe.N1,dNdx,r,Sigma,dVol); + this->formKG(elMat.A[eKg-1],fe.N1,dNdx,r,Sigma,dVol); } // Radial shape function for axi-symmetric problems @@ -434,9 +443,11 @@ bool NonlinearElasticityULMixed::evalIntMx (LocalIntegral*& elmInt, // Integrate the material stiffness matrix Dmat *= dVol; if (nsd == 2) - acckm2d_(axiSymmetry,Nr.size(),Nr.ptr(),dNdx.ptr(),Dmat.ptr(),eKm->ptr()); + acckm2d_(axiSymmetry,Nr.size(),Nr.ptr(),dNdx.ptr(),Dmat.ptr(), + elMat.A[eKm-1].ptr()); else - acckm3d_(fe.N1.size(),dNdx.ptr(),Dmat.ptr(),eKm->ptr()); + acckm3d_(fe.N1.size(),dNdx.ptr(),Dmat.ptr(), + elMat.A[eKm-1].ptr()); #endif // Integrate the volumetric change and pressure tangent terms @@ -450,55 +461,55 @@ bool NonlinearElasticityULMixed::evalIntMx (LocalIntegral*& elmInt, { for (i = 1; i <= nsd; i++) { - myMats->A[Kut](nsd*(a-1)+i,b) += dNdx(a,i)*fe.N2(b) * Dmat(i,7); + elMat.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)*fe.N2(b) * Dmat(4,7); + elMat.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)*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); + elMat.A[Kut](nsd*(a-1)+i,b) += dNdx(a,3-i)*fe.N2(b) * Dmat(4,7); + elMat.A[Kut](nsd*(a-1)+j,b) += dNdx(a,5-j)*fe.N2(b) * Dmat(5,7); + elMat.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)*fe.N2(b) * dVup; - } + elMat.A[Kup](nsd*(a-1)+i,b) += dNdx(a,i)*fe.N2(b) * dVup; + } if (axiSymmetry) { double Nr = fe.N1(a)/r; - myMats->A[Kut](nsd*(a-1)+1,b) += Nr*fe.N2(b) * Dmat(3,7); - myMats->A[Kup](nsd*(a-1)+1,b) += Nr*fe.N2(b) * dVup; + elMat.A[Kut](nsd*(a-1)+1,b) += Nr*fe.N2(b) * Dmat(3,7); + elMat.A[Kup](nsd*(a-1)+1,b) += Nr*fe.N2(b) * dVup; } } - 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*(-detJW),true); // -= N2*N2^T*|J| + elMat.A[Ktt].outer_product(fe.N2,fe.N2*Dmat(7,7),true); // += N2*N2^T*D77 + elMat.A[Ktp].outer_product(fe.N2,fe.N2*(-detJW),true); // -= N2*N2^T*|J| if (lHaveStrains) { // Compute the small-deformation strain-displacement matrix B from dNdx if (axiSymmetry) - this->formBmatrix(fe.N1,dNdx,r); + this->formBmatrix(B,fe.N1,dNdx,r); else - this->formBmatrix(dNdx); + this->formBmatrix(B,dNdx); // Integrate the internal forces Sigma *= -dVol; - if (!Bmat.multiply(Sigma,myMats->b[Ru],true,true)) // -= B^T*sigma + if (!B.multiply(Sigma,elMat.b[Ru],true,true)) // -= B^T*sigma return false; // Integrate the volumetric change and pressure forces - myMats->b[Rt].add(fe.N2,(Press - Bpres)*detJW); // += N2*(p-pBar)*|J| - myMats->b[Rp].add(fe.N2,(Theta - J)*detJW); // += N2*(Theta-J)*|J| + elMat.b[Rt].add(fe.N2,(Press - Bpres)*detJW); // += N2*(p-pBar)*|J| + elMat.b[Rp].add(fe.N2,(Theta - J)*detJW); // += N2*(Theta-J)*|J| #if INT_DEBUG > 4 std::cout <<"NonlinearElasticityULMixed::Sigma*dVol =\n"<< Sigma; - std::cout <<"NonlinearElasticityULMixed::Bmat ="<< Bmat; - std::cout <<"NonlinearElasticityULMixed::Ru("<< iP <<") ="<< myMats->b[Ru]; + std::cout <<"NonlinearElasticityULMixed::Bmat ="<< B; + std::cout <<"NonlinearElasticityULMixed::Ru("<< iP <<") ="<< elMat.b[Ru]; #endif } - return this->getIntegralResult(elmInt); + return true; } @@ -508,7 +519,7 @@ bool NonlinearElasticityULMixed::evalIntMx (LocalIntegral*& elmInt, forwarded to the corresponding single-field method of the parent class. */ -bool NonlinearElasticityULMixed::evalBouMx (LocalIntegral*& elmInt, +bool NonlinearElasticityULMixed::evalBouMx (LocalIntegral& elmInt, const MxFiniteElement& fe, const Vec3& X, const Vec3& normal) const @@ -557,7 +568,7 @@ NormBase* NonlinearElasticityULMixed::getNormIntegrand (AnaSol*) const } -bool ElasticityNormULMixed::evalIntMx (LocalIntegral*& elmInt, +bool ElasticityNormULMixed::evalIntMx (LocalIntegral& elmInt, const MxFiniteElement& fe, const TimeDomain& prm, const Vec3& X) const @@ -566,40 +577,44 @@ bool ElasticityNormULMixed::evalIntMx (LocalIntegral*& elmInt, ulp = static_cast(&myProblem); // Evaluate the deformation gradient, F, and the Green-Lagrange strains, E + Matrix B; Tensor F(ulp->nDF); SymmTensor E(ulp->nDF); - if (!ulp->kinematics(fe.N1,fe.dN1dX,X.x,F,E)) + if (!ulp->kinematics(elmInt.vec[U],fe.N1,fe.dN1dX,X.x,B,F,E)) return false; // Evaluate the volumetric change field - double Theta = fe.N2.dot(ulp->mySols[T]) + 1.0; + double Theta = fe.N2.dot(elmInt.vec[T]) + 1.0; // Compute the mixed model deformation gradient, F_bar - ulp->Fbar = F; // notice that F_bar always has dimension 3 + Tensor Fbar(3); // notice that F_bar always has dimension 3 + Fbar = F; // notice that F_bar always has dimension 3 double r1 = pow(fabs(Theta/F.det()),1.0/3.0); - ulp->Fbar *= r1; + Fbar *= r1; if (F.dim() == 2) // In 2D we always assume plane strain so set F(3,3)=1 - ulp->Fbar(3,3) = r1; + Fbar(3,3) = r1; // Compute the strain energy density, U(E) = Int_E (Sig:Eps) dEps // and the Cauchy stress tensor, Sig + Matrix Cmat; double U = 0.0; SymmTensor Sig(3); - if (!ulp->material->evaluate(ulp->Cmat,Sig,U,X,ulp->Fbar,E,3,&prm,&F)) + if (!ulp->material->evaluate(Cmat,Sig,U,X,Fbar,E,3,&prm,&F)) return false; // Axi-symmetric integration point volume; 2*pi*r*|J|*w double detJW = ulp->isAxiSymmetric() ? 2.0*M_PI*X.x*fe.detJxW : fe.detJxW; // Integrate the norms - return evalInt(getElmNormBuffer(elmInt,6),Sig,U,Theta,detJW); + return evalInt(static_cast(elmInt),Sig,U,Theta,detJW); } -bool ElasticityNormULMixed::evalBouMx (LocalIntegral*& elmInt, +bool ElasticityNormULMixed::evalBouMx (LocalIntegral& elmInt, const MxFiniteElement& fe, const Vec3& X, const Vec3& normal) const { return this->ElasticityNormUL::evalBou(elmInt,fe,X,normal); } + diff --git a/Apps/FiniteDefElasticity/NonlinearElasticityULMixed.h b/Apps/FiniteDefElasticity/NonlinearElasticityULMixed.h index dd6f18af..d38325f3 100644 --- a/Apps/FiniteDefElasticity/NonlinearElasticityULMixed.h +++ b/Apps/FiniteDefElasticity/NonlinearElasticityULMixed.h @@ -57,24 +57,35 @@ public: //! \param[in] mode The solution mode to use virtual void setMode(SIM::SolutionMode mode); + //! \brief Returns a local integral container for the given element. + //! \param[in] nen1 Number of nodes on element for basis 1 + //! \param[in] nen2 Number of nodes on element for basis 2 + //! \param[in] neumann Whether or not we are assembling Neumann BC's + virtual LocalIntegral* getLocalIntegral(size_t nen1, size_t nen2, size_t, + bool neumann) const; + //! \brief Initializes current element for numerical integration. //! \param[in] MNPC1 Nodal point correspondance for the basis 1 //! \param[in] MNPC2 Nodal point correspondance for the basis 2 //! \param[in] n1 Number of nodes in basis 1 on this patch + //! \param elmInt The local integral object for current element virtual bool initElement(const std::vector& MNPC1, - const std::vector& MNPC2, size_t n1); + const std::vector& MNPC2, size_t n1, + LocalIntegral& elmInt); //! \brief Initializes current element for numerical boundary integration. //! \param[in] MNPC1 Nodal point correspondance for the basis 1 //! \param[in] MNPC2 Nodal point correspondance for the basis 2 + //! \param elmInt The local integral object for current element virtual bool initElementBou(const std::vector& MNPC1, - const std::vector& MNPC2, size_t); + const std::vector& MNPC2, size_t, + LocalIntegral& elmInt); //! \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] X Cartesian coordinates of current integration point - virtual bool evalIntMx(LocalIntegral*& elmInt, const MxFiniteElement& fe, + virtual bool evalIntMx(LocalIntegral& elmInt, const MxFiniteElement& fe, const TimeDomain& prm, const Vec3& X) const; //! \brief Evaluates the integrand at a boundary point. @@ -82,7 +93,7 @@ public: //! \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 evalBouMx(LocalIntegral*& elmInt, const MxFiniteElement& fe, + virtual bool evalBouMx(LocalIntegral& elmInt, const MxFiniteElement& fe, const Vec3& X, const Vec3& normal) const; //! \brief Returns a pointer to an Integrand for solution norm evaluation. @@ -102,10 +113,6 @@ public: //! \param[in] prefix Name prefix for all components virtual const char* getField1Name(size_t i, const char* prefix = 0) const; -protected: - mutable Tensor Fbar; //!< Mixed model deformation gradient - mutable Matrix Dmat; //!< Projected mixed constitutive matrix - friend class ElasticityNormULMixed; }; @@ -128,7 +135,7 @@ public: //! \param[in] fe Mixed finite element data of current integration point //! \param[in] prm Nonlinear solution algorithm parameters //! \param[in] X Cartesian coordinates of current integration point - virtual bool evalIntMx(LocalIntegral*& elmInt, const MxFiniteElement& fe, + virtual bool evalIntMx(LocalIntegral& elmInt, const MxFiniteElement& fe, const TimeDomain& prm, const Vec3& X) const; //! \brief Evaluates the integrand at a boundary point (mixed). @@ -136,7 +143,7 @@ public: //! \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 evalBouMx(LocalIntegral*& elmInt, const MxFiniteElement& fe, + virtual bool evalBouMx(LocalIntegral& elmInt, const MxFiniteElement& fe, const Vec3& X, const Vec3& normal) const; }; diff --git a/Apps/FiniteDefElasticity/main_NonLinEl.C b/Apps/FiniteDefElasticity/main_NonLinEl.C index ea7972de..ee4901c5 100644 --- a/Apps/FiniteDefElasticity/main_NonLinEl.C +++ b/Apps/FiniteDefElasticity/main_NonLinEl.C @@ -25,7 +25,9 @@ #include #include +#ifndef USE_OPENMP extern std::vector mixedDbgEl; //!< List of elements for additional output +#endif /*! @@ -146,9 +148,11 @@ int main (int argc, char** argv) else --i; } +#ifndef USE_OPENMP else if (!strcmp(argv[i],"-dbgElm")) while (i < argc-1 && isdigit(argv[i+1][0])) utl::parseIntegers(mixedDbgEl,argv[++i]); +#endif else if (!strcmp(argv[i],"-ignore")) while (i < argc-1 && isdigit(argv[i+1][0])) utl::parseIntegers(ignoredPatches,argv[++i]);