OilfieldToolsNet/IFEM
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Added a global integration point counter needed by some integrands. It needs to be set up outside the integrands when multi-threading.

Browse Source 
git-svn-id: http://svn.sintef.no/trondheim/IFEM/trunk@1423 e10b68d5-8a6e-419e-a041-bce267b0401d
This commit is contained in:
kmo 2012-01-19 16:59:50 +00:00 committed by Knut Morten Okstad
parent 5ad7f3afda
commit 01f4bc5fb3
14 changed files with 125 additions and 29 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

41
Apps/FiniteDefElasticity/NonlinearElasticityUL.C
View File

@ -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
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;
}
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);
}
// 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<ElmNorm&>(elmInt)[1] += Ux[iP++]*detJW;
static_cast<ElmNorm&>(elmInt)[1] += Ux[iP]*detJW;
return true;
}

11
Apps/FiniteDefElasticity/NonlinearElasticityUL.h
View File

@ -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<Vec3> up; //!< Previous displacements
mutable std::vector<Vec3> tp; //!< Previous tractions

5
src/ASM/ASMbase.h
View File

@ -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.

3
src/ASM/ASMs1D.C
View File

@ -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

4
src/ASM/ASMs2D.C
View File

@ -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)

6
src/ASM/ASMs3D.C
View File

@ -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);

20
src/ASM/ASMstruct.C
View File

@ -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
}

7
src/ASM/ASMstruct.h
View File

@ -56,9 +56,16 @@ 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
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

5
src/ASM/FiniteElement.h
View File

@ -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
};

3
src/ASM/Integrand.h
View File

@ -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.

5
src/ASM/IntegrandBase.h
View File

@ -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);

9
src/ASM/LR/ASMunstruct.h
View File

@ -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,9 +53,14 @@ 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

2
src/Integrands/MaterialBase.h
View File

@ -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.

17
src/SIM/SIMbase.C
View File

@ -34,12 +34,10 @@
#include "VTF.h"
#include "Functions.h"
#include "Utilities.h"
#include "tinyxml.h"
#include <fstream>
#include <sstream>
#include <iomanip>
#include <fstream>
#include "tinyxml.h"
ASM::Discretization SIMbase::discretization = ASM::Spline;
@ -660,9 +658,13 @@ bool SIMbase::preprocess (const std::vector<int>& 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<int>& 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<int>& 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();)
{