diff --git a/Apps/FiniteDefElasticity/NonlinearElasticityUL.C b/Apps/FiniteDefElasticity/NonlinearElasticityUL.C index f8d0b332..c2620b60 100644 --- a/Apps/FiniteDefElasticity/NonlinearElasticityUL.C +++ b/Apps/FiniteDefElasticity/NonlinearElasticityUL.C @@ -126,6 +126,15 @@ LocalIntegral* NonlinearElasticityUL::getLocalIntegral (size_t nen, size_t, } +void NonlinearElasticityUL::initIntegration (size_t nGp, size_t nBp) +{ + if (material) + material->initIntegration(nGp); + + this->Elasticity::initIntegration(nGp,nBp); +} + + void NonlinearElasticityUL::initIntegration (const TimeDomain& prm) { if (material) @@ -400,10 +409,19 @@ NormBase* NonlinearElasticityUL::getNormIntegrand (AnaSol*) const } +void ElasticityNormUL::initIntegration (size_t nGp, size_t nBp) +{ + Ux.resize(nBp,0.0); + up.resize(nBp); + tp.resize(nBp); + + this->ElasticityNorm::initIntegration(nGp,nBp); +} + + void ElasticityNormUL::initIntegration (const TimeDomain& prm) { this->ElasticityNorm::initIntegration(prm); - iP = 0; } @@ -474,24 +492,23 @@ bool ElasticityNormUL::evalBou (LocalIntegral& elmInt, Vec3 u = ulp.evalSol(elmInt.vec.front(),fe.N); // Integrate the external energy (path integral) - if (iP < Ux.size()) + size_t iP = fe.iGP; +#ifdef INDEX_CHECK + if (iP > Ux.size()) { - // New load step, update integration point values - Ux[iP] += 0.5*(t+tp[iP])*(u-up[iP]); - tp[iP] = t; - up[iP] = u; - } - else - { - // This is the first load step at this integration point - Ux.push_back(0.5*t*u); - tp.push_back(t); - up.push_back(u); + std::cerr <<" *** ElasticityNormUL::evalBou: Integration point "<< iP + <<" out of range [0,"<< Ux.size() ">"<< std::endl; + return false; } +#endif + + Ux[iP] += 0.5*(t+tp[iP])*(u-up[iP]); + tp[iP] = t; + up[iP] = u; // Axi-symmetric integration point volume; 2*pi*r*|J|*w double detJW = ulp.isAxiSymmetric() ? 2.0*M_PI*X.x*fe.detJxW : fe.detJxW; - static_cast(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 df29e325..93e608d4 100644 --- a/Apps/FiniteDefElasticity/NonlinearElasticityUL.h +++ b/Apps/FiniteDefElasticity/NonlinearElasticityUL.h @@ -53,6 +53,10 @@ public: virtual LocalIntegral* getLocalIntegral(size_t nen, size_t, bool neumann) const; + //! \brief Initializes the integrand with the number of integration points. + //! \param[in] nGp Total number of interior integration points + //! \param[in] nBp Total number of boundary integration points + virtual void initIntegration(size_t nGp, size_t nBp); //! \brief Initializes the integrand for a new integration loop. //! \param[in] prm Nonlinear solution algorithm parameters virtual void initIntegration(const TimeDomain& prm); @@ -130,10 +134,14 @@ class ElasticityNormUL : public ElasticityNorm public: //! \brief The only constructor initializes its data members. //! \param[in] p The linear elasticity problem to evaluate norms for - ElasticityNormUL(NonlinearElasticityUL& p) : ElasticityNorm(p), iP(0) {} + ElasticityNormUL(NonlinearElasticityUL& p) : ElasticityNorm(p) {} //! \brief Empty destructor. virtual ~ElasticityNormUL() {} + //! \brief Initializes the integrand with the number of integration points. + //! \param[in] nGp Total number of interior integration points + //! \param[in] nBp Total number of boundary integration points + virtual void initIntegration(size_t nGp, size_t nBp); //! \brief Initializes the integrand for a new integration loop. //! \param[in] prm Nonlinear solution algorithm parameters virtual void initIntegration(const TimeDomain& prm); @@ -172,7 +180,6 @@ protected: private: // Data for path-integral of the external energy due to boundary tractions - mutable size_t iP; //!< Global integration point counter mutable RealArray Ux; //!< External energy densities mutable std::vector up; //!< Previous displacements mutable std::vector tp; //!< Previous tractions diff --git a/src/ASM/ASMbase.h b/src/ASM/ASMbase.h index c764e1c6..b4674c1d 100644 --- a/src/ASM/ASMbase.h +++ b/src/ASM/ASMbase.h @@ -247,6 +247,11 @@ public: //! \brief Initializes the patch level MADOF array for mixed problems. virtual void initMADOF(const int*) {} + //! \brief Computes the total number of integration points in this patch. + virtual void getNoIntPoints(size_t& nPt) = 0; + //! \brief Computes the number of boundary integration points in this patch. + virtual void getNoBouPoints(size_t& nPt, int ldim, int lindx) = 0; + // Methods for integration of finite element quantities. // These are the main computational methods of the ASM class hierarchy. diff --git a/src/ASM/ASMs1D.C b/src/ASM/ASMs1D.C index 48693e31..5267ac2d 100644 --- a/src/ASM/ASMs1D.C +++ b/src/ASM/ASMs1D.C @@ -598,6 +598,8 @@ bool ASMs1D::integrate (Integrand& integrand, for (int i = 0; i < nGauss; i++) { + fe.iGP = firstIp + (iel-1)*nGauss + i; // Global integration point counter + // Local element coordinate of current integration point fe.xi = xg[i]; @@ -672,6 +674,7 @@ bool ASMs1D::integrate (Integrand& integrand, int lIndex, return false; } + fe.iGP = firstBp; fe.iel = MLGE[iel-1]; if (fe.iel < 1) return true; // zero-length element diff --git a/src/ASM/ASMs2D.C b/src/ASM/ASMs2D.C index 87dab177..979736f1 100644 --- a/src/ASM/ASMs2D.C +++ b/src/ASM/ASMs2D.C @@ -1139,6 +1139,8 @@ bool ASMs2D::integrate (Integrand& integrand, for (int j = 0; j < nGauss; j++, ip += nGauss*(nel1-1)) for (int i = 0; i < nGauss; i++, ip++) { + fe.iGP = firstIp + ip; // Global integration point counter + // Local element coordinates of current integration point fe.xi = xg[i]; fe.eta = xg[j]; @@ -1306,6 +1308,8 @@ bool ASMs2D::integrate (Integrand& integrand, int lIndex, int ip = (t1 == 1 ? i2-p2 : i1-p1)*nGauss; for (int i = 0; i < nGauss; i++, ip++) { + fe.iGP = firstBp + ip; // Global integration point counter + // Local element coordinates and parameter values // of current integration point if (gpar[0].size() > 1) diff --git a/src/ASM/ASMs3D.C b/src/ASM/ASMs3D.C index 87e9397f..691a2968 100644 --- a/src/ASM/ASMs3D.C +++ b/src/ASM/ASMs3D.C @@ -1345,6 +1345,8 @@ bool ASMs3D::integrate (Integrand& integrand, for (int j = 0; j < nGauss; j++, ip += nGauss*(nel1-1)) for (int i = 0; i < nGauss; i++, ip++) { + fe.iGP = firstIp + ip; // Global integration point counter + // Local element coordinates of current integration point fe.xi = xg[i]; fe.eta = xg[j]; @@ -1534,6 +1536,8 @@ bool ASMs3D::integrate (Integrand& integrand, int lIndex, for (int j = 0; j < nGauss; j++, ip += nGauss*(nf1-1)) for (int i = 0; i < nGauss; i++, ip++) { + fe.iGP = firstBp + ip; // Global integration point counter + // Local element coordinates and parameter values // of current integration point switch (abs(faceDir)) { @@ -1730,6 +1734,8 @@ bool ASMs3D::integrateEdge (Integrand& integrand, int lEdge, for (int i = 0; i < nGauss; i++, ip++) { + fe.iGP = firstBp + ip; // Global integration point counter + // Parameter values of current integration point if (gpar[0].size() > 1) fe.u = gpar[0](i+1,i1-p1+1); if (gpar[1].size() > 1) fe.v = gpar[1](i+1,i2-p2+1); diff --git a/src/ASM/ASMstruct.C b/src/ASM/ASMstruct.C index 541c0f02..6820c73a 100644 --- a/src/ASM/ASMstruct.C +++ b/src/ASM/ASMstruct.C @@ -38,3 +38,23 @@ ASMstruct::~ASMstruct () { if (geo && !shareFE) delete geo; } + + +void ASMstruct::getNoIntPoints (size_t& nPt) +{ + size_t nGp = 1; + for (unsigned char d = 0; d < ndim; d++) + nGp *= nGauss; + + firstIp = nPt; + + nPt += this->getNoElms(true)*nGp; // Note: Includes also the 0-span elements +} + + +void ASMstruct::getNoBouPoints (size_t& nPt, int ldim, int lindx) +{ + firstBp = nPt; + if (ldim == 0) nPt ++; // This is correct for 1D only, 2D and 3D must take + // care of how many element do we have on the boundary defined by lindx +} diff --git a/src/ASM/ASMstruct.h b/src/ASM/ASMstruct.h index 0c0e6ed4..58bfe6aa 100644 --- a/src/ASM/ASMstruct.h +++ b/src/ASM/ASMstruct.h @@ -56,12 +56,19 @@ public: //! \param[in] integrand Object with problem-specific data and methods virtual Go::GeomObject* evalSolution(const Integrand& integrand) const = 0; + //! \brief Computes the total number of integration points in this patch. + virtual void getNoIntPoints(size_t& nPt); + //! \brief Computes the number of boundary integration points in this patch. + virtual void getNoBouPoints(size_t& nPt, int ldim, int lindx); + protected: Go::GeomObject* geo; //!< Pointer to the actual spline geometry object - int nGauss; //!< Numerical integration scheme - static int gEl; //!< Global element counter - static int gNod; //!< Global node counter + size_t firstIp; //!< Global index to first interior integration point + size_t firstBp; //!< Global index to first boundary integration point + int nGauss; //!< Numerical integration scheme + static int gEl; //!< Global element counter + static int gNod; //!< Global node counter }; #endif diff --git a/src/ASM/FiniteElement.h b/src/ASM/FiniteElement.h index c8909a27..cbf67331 100644 --- a/src/ASM/FiniteElement.h +++ b/src/ASM/FiniteElement.h @@ -25,9 +25,10 @@ class FiniteElement { public: //! \brief Default constructor. - FiniteElement(size_t nb = 0) : iel(0), h(0.0), N(nb), detJxW(1.0) {} + FiniteElement(size_t nb = 0) : iel(0), iGP(0), h(0.0), N(nb), detJxW(1.0) {} int iel; //!< Element identifier + size_t iGP; //!< Global integration point counter double u; //!< First parameter of current point double v; //!< Second parameter of current point double w; //!< Third parameter of current point @@ -39,8 +40,8 @@ public: Vector Navg; //!< Volume-averaged basis function values Matrix dNdX; //!< First derivatives (gradient) of the basis functions Matrix3D d2NdX2; //!< Second derivatives of the basis functions - double detJxW; //!< Weighted determinant of the coordinate mapping Matrix G; //!< Matrix used for stabilized methods + double detJxW; //!< Weighted determinant of the coordinate mapping }; diff --git a/src/ASM/Integrand.h b/src/ASM/Integrand.h index 2bf658be..3e6cdbe1 100644 --- a/src/ASM/Integrand.h +++ b/src/ASM/Integrand.h @@ -52,6 +52,9 @@ public: // Global initialization interface // =============================== + //! \brief Initializes the integrand with the number of integration points. + //! \details This method is invoked only once during the preprocessing stage. + virtual void initIntegration(size_t, size_t) = 0; //! \brief Initializes the integrand for a new integration loop. //! \details This method is invoked once before starting the numerical //! integration over the entire spatial domain. diff --git a/src/ASM/IntegrandBase.h b/src/ASM/IntegrandBase.h index 92ee9746..75fa3ad1 100644 --- a/src/ASM/IntegrandBase.h +++ b/src/ASM/IntegrandBase.h @@ -47,6 +47,9 @@ public: //! \brief Defines the solution mode before the element assembly is started. virtual void setMode(SIM::SolutionMode mode) { m_mode = mode; } + //! \brief Initializes the integrand with the number of integration points. + //! \details This method is invoked only once during the preprocessing stage. + virtual void initIntegration(size_t, size_t) {} //! \brief Initializes the integrand for a new integration loop. //! \details This method is invoked once before starting the numerical //! integration over the entire spatial domain. @@ -199,6 +202,8 @@ public: //! \brief Empty destructor. virtual ~NormBase() {} + //! \brief Initializes the integrand with the number of integration points. + virtual void initIntegration(size_t, size_t) {} //! \brief Initializes the integrand for a new integration loop. virtual void initIntegration(const TimeDomain& time); diff --git a/src/ASM/LR/ASMunstruct.h b/src/ASM/LR/ASMunstruct.h index 423575ae..0cbc8dde 100644 --- a/src/ASM/LR/ASMunstruct.h +++ b/src/ASM/LR/ASMunstruct.h @@ -18,12 +18,12 @@ namespace Go { class GeomObject; - class BoundingBox; } namespace LR { class LRSplineSurface; } + /*! \brief Base class for unstructured spline-based FE assembly drivers. \details This class contains methods common for unstructured spline patches. @@ -53,15 +53,20 @@ public: //! \param[in] ng Number of Gauss points in each parameter direction void setGauss(int ng) { nGauss = ng; } - //! \brief Resets global element and node counters + //! \brief Resets global element and node counters. static void resetNumbering() { gEl = gNod = 0; } + //! \brief Computes the total number of integration points in this patch. + virtual void getNoIntPoints(size_t& nPt) { nPt = 0; } // later... + //! \brief Computes the number of boundary integration points in this patch. + virtual void getNoBouPoints(size_t& nPt, int, int) { nPt = 0; } // later... + protected: LR::LRSplineSurface* geo; //!< Pointer to the actual spline geometry object - int nGauss; //!< Numerical integration scheme - static int gEl; //!< Global element counter - static int gNod; //!< Global node counter + int nGauss; //!< Numerical integration scheme + static int gEl; //!< Global element counter + static int gNod; //!< Global node counter }; diff --git a/src/Integrands/MaterialBase.h b/src/Integrands/MaterialBase.h index 38e1fe2b..350557dd 100644 --- a/src/Integrands/MaterialBase.h +++ b/src/Integrands/MaterialBase.h @@ -42,6 +42,8 @@ public: //! \brief Prints out material parameters to the given output stream. virtual void print(std::ostream&) const {} + //! \brief Initializes the material with the number of integration points. + virtual void initIntegration(size_t) {} //! \brief Initializes the material model for a new integration loop. virtual void initIntegration(const TimeDomain&) {} //! \brief Initializes the material model for a new result point loop. diff --git a/src/SIM/SIMbase.C b/src/SIM/SIMbase.C index 0587452d..7678d6f8 100644 --- a/src/SIM/SIMbase.C +++ b/src/SIM/SIMbase.C @@ -34,12 +34,10 @@ #include "VTF.h" #include "Functions.h" #include "Utilities.h" +#include "tinyxml.h" #include #include #include -#include - -#include "tinyxml.h" ASM::Discretization SIMbase::discretization = ASM::Spline; @@ -660,9 +658,13 @@ bool SIMbase::preprocess (const std::vector& ignored, bool fixDup) if (renum > 0) std::cout <<"\nRenumbered "<< renum <<" nodes."<< std::endl; + size_t nIntP = 0, nBouP = 0; ASMs2DC1::renumberNodes(*g2l); for (mit = myModel.begin(); mit != myModel.end(); mit++) + { (*mit)->renumberNodes(*g2l); + (*mit)->getNoIntPoints(nIntP); // count number of global integration points + } // Process the specified Dirichlet boundary conditions std::cout <<"\nResolving Dirichlet boundary conditions"<< std::endl; @@ -685,6 +687,12 @@ bool SIMbase::preprocess (const std::vector& ignored, bool fixDup) ok = false; break; + case Property::NEUMANN: + // Count number of global boundary integration points + if (p->patch > 0 && p->patch <= myModel.size()) + myModel[p->patch-1]->getNoBouPoints(nBouP,p->ldim,p->lindx); + break; + default: break; } @@ -701,6 +709,9 @@ bool SIMbase::preprocess (const std::vector& ignored, bool fixDup) // Resolve possibly chaining of the MPC equations if (!allMPCs.empty()) ASMbase::resolveMPCchains(allMPCs); + // Let the integrand know how many integration points in total do we have + if (myProblem) myProblem->initIntegration(nIntP,nBouP); + // Preprocess the result points for (ResPointVec::iterator p = myPoints.begin(); p != myPoints.end();) {