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

Kernel change: Added base class IntegrandBase with some generic functionality, Integrand is now just an abstract interface

Browse Source 
git-svn-id: http://svn.sintef.no/trondheim/IFEM/trunk@1074 e10b68d5-8a6e-419e-a041-bce267b0401d
This commit is contained in:
kmo
2011-06-23 18:13:16 +00:00
committed by Knut Morten Okstad
parent d652c21dad
commit 0de9f510bf
15 changed files with 650 additions and 528 deletions
Download Patch File Download Diff File Expand all files Collapse all files

11
Apps/FiniteDefElasticity/NonlinearElasticityTL.C
View File

@@ -47,13 +47,13 @@ void NonlinearElasticityTL::setMode (SIM::SolutionMode mode)
switch (mode)
{
case SIM::STATIC:
myMats->resize(1,2);
myMats->resize(1,1);
eKm = &myMats->A[0];
eKg = &myMats->A[0];
break;
case SIM::DYNAMIC:
myMats->resize(2,2);
myMats->resize(2,1);
eKm = &myMats->A[0];
eKg = &myMats->A[0];
eM = &myMats->A[1];
@@ -61,9 +61,9 @@ void NonlinearElasticityTL::setMode (SIM::SolutionMode mode)
case SIM::RHS_ONLY:
if (myMats->A.empty())
myMats->resize(1,2);
myMats->resize(1,1);
else
myMats->b.resize(2);
myMats->b.resize(1);
eKm = &myMats->A[0];
eKg = &myMats->A[0];
myMats->rhsOnly = true;
@@ -76,9 +76,10 @@ void NonlinearElasticityTL::setMode (SIM::SolutionMode mode)
}
// We always need the force+displacement vectors in nonlinear simulations
mySols.resize(1);
iS = &myMats->b[0];
eS = &myMats->b[0];
eV = &myMats->b[1];
eV = &mySols[0];
tracVal.clear();
}

25
Apps/FiniteDefElasticity/NonlinearElasticityUL.C
View File

@@ -42,7 +42,6 @@ void NonlinearElasticityUL::setMode (SIM::SolutionMode mode)
{
if (!myMats) return;
size_t nvec = 1 + primsol.size();
myMats->rhsOnly = false;
eM = eKm = eKg = 0;
iS = eS = eV = 0;
@@ -50,13 +49,13 @@ void NonlinearElasticityUL::setMode (SIM::SolutionMode mode)
switch (mode)
{
case SIM::STATIC:
myMats->resize(1,nvec);
myMats->resize(1,1);
eKm = &myMats->A[0];
eKg = &myMats->A[0];
break;
case SIM::DYNAMIC:
myMats->resize(2,nvec);
myMats->resize(2,1);
eKm = &myMats->A[0];
eKg = &myMats->A[0];
eM = &myMats->A[1];
@@ -64,9 +63,9 @@ void NonlinearElasticityUL::setMode (SIM::SolutionMode mode)
case SIM::RHS_ONLY:
if (myMats->A.empty())
myMats->resize(1,nvec);
myMats->resize(1,1);
else
myMats->b.resize(nvec);
myMats->b.resize(1);
eKm = &myMats->A[0];
eKg = &myMats->A[0];
myMats->rhsOnly = true;
@@ -80,7 +79,8 @@ void NonlinearElasticityUL::setMode (SIM::SolutionMode mode)
// We always need the force+displacement vectors in nonlinear simulations
iS = &myMats->b[0];
eS = &myMats->b[0];
eV = &myMats->b[1];
mySols.resize(1);
eV = &mySols[0];
tracVal.clear();
}
@@ -330,9 +330,9 @@ bool ElasticityNormUL::evalInt (LocalIntegral*& elmInt,
const TimeDomain& prm,
const Vec3& X) const
{
ElmNorm& pnorm = ElasticityNorm::getElmNormBuffer(elmInt);
NonlinearElasticityUL& ulp = static_cast<NonlinearElasticityUL&>(myProblem);
NonlinearElasticityUL& ulp = static_cast<NonlinearElasticityUL&>(problem);
ElmNorm& pnorm = NormBase::getElmNormBuffer(elmInt);
// Evaluate the deformation gradient, F, and the Green-Lagrange strains, E
Tensor F(fe.dNdX.cols());
@@ -357,14 +357,15 @@ bool ElasticityNormUL::evalBou (LocalIntegral*& elmInt,
const FiniteElement& fe,
const Vec3& X, const Vec3& normal) const
{
if (!problem.haveLoads()) return true;
NonlinearElasticityUL& ulp = static_cast<NonlinearElasticityUL&>(myProblem);
if (!ulp.haveLoads()) return true;
ElmNorm& pnorm = ElasticityNorm::getElmNormBuffer(elmInt);
ElmNorm& pnorm = NormBase::getElmNormBuffer(elmInt);
// Evaluate the current surface traction
Vec3 t = problem.getTraction(X,normal);
Vec3 t = ulp.getTraction(X,normal);
// Evaluate the current displacement field
Vec3 u = problem.evalSol(fe.N);
Vec3 u = ulp.evalSol(fe.N);
// Integrate the external energy (path integral)
if (iP < Ux.size())

73
Apps/FiniteDefElasticity/NonlinearElasticityULMX.C
View File

@@ -73,14 +73,14 @@ void NonlinearElasticityULMX::setMode (SIM::SolutionMode mode)
switch (mode)
{
case SIM::STATIC:
myMats->resize(2,3);
myMats->resize(2,1);
eKm = &myMats->A[0];
eKg = &myMats->A[0];
Hh = &myMats->A[1];
break;
case SIM::DYNAMIC:
myMats->resize(3,3);
myMats->resize(3,1);
eKm = &myMats->A[0];
eKg = &myMats->A[0];
eM = &myMats->A[1];
@@ -89,9 +89,9 @@ void NonlinearElasticityULMX::setMode (SIM::SolutionMode mode)
case SIM::RHS_ONLY:
if (myMats->A.size() < 2)
myMats->resize(2,3);
myMats->resize(2,1);
else
myMats->b.resize(3);
myMats->b.resize(1);
eKm = &myMats->A[0];
eKg = &myMats->A[0];
Hh = &myMats->A.back();
@@ -99,12 +99,12 @@ void NonlinearElasticityULMX::setMode (SIM::SolutionMode mode)
break;
case SIM::RECOVERY:
if (myMats->A.size() < 1)
myMats->resize(1,3);
else
myMats->b.resize(3);
if (myMats->A.empty())
myMats->A.resize(1);
if (mySols.empty())
mySols.resize(1);
Hh = &myMats->A.back();
eV = &myMats->b[1];
eV = &mySols[0];
myMats->rhsOnly = true;
return;
@@ -116,7 +116,8 @@ void NonlinearElasticityULMX::setMode (SIM::SolutionMode mode)
// We always need the force+displacement vectors in nonlinear simulations
iS = &myMats->b[0];
eS = &myMats->b[0];
eV = &myMats->b[1];
mySols.resize(2);
eV = &mySols[0];
tracVal.clear();
}
@@ -144,21 +145,22 @@ bool NonlinearElasticityULMX::evalInt (LocalIntegral*& elmInt,
const FiniteElement& fe,
const TimeDomain&, const Vec3& X) const
{
if (myMats->b.size() < 3) return false;
if (mySols.size() < 2) return false;
ItgPtData& ptData = myData[iP++];
Vectors& eVs = const_cast<Vectors&>(mySols);
#if INT_DEBUG > 1
std::cout <<"NonlinearElasticityUL::dNdX ="<< fe.dNdX;
#endif
// Evaluate the deformation gradient, Fp, at previous configuration
const_cast<NonlinearElasticityULMX*>(this)->eV = &myMats->b[2];
const_cast<NonlinearElasticityULMX*>(this)->eV = &eVs[1];
if (!this->formDefGradient(fe.dNdX,ptData.Fp))
return false;
// Evaluate the deformation gradient, F, at current configuration
const_cast<NonlinearElasticityULMX*>(this)->eV = &myMats->b[1];
const_cast<NonlinearElasticityULMX*>(this)->eV = &eVs[0];
if (!this->formDefGradient(fe.dNdX,ptData.F))
return false;
@@ -456,16 +458,16 @@ NormBase* NonlinearElasticityULMX::getNormIntegrand (AnaSol*) const
bool ElasticityNormULMX::initElement (const std::vector<int>& MNPC,
const Vec3& Xc, size_t nPt)
{
NonlinearElasticityULMX& elp = static_cast<NonlinearElasticityULMX&>(problem);
NonlinearElasticityULMX& p = static_cast<NonlinearElasticityULMX&>(myProblem);
elp.iP = 0;
elp.X0 = Xc;
elp.myData.resize(nPt);
size_t nPm = utl::Pascal(elp.p,elp.nsd);
if (elp.Hh) elp.Hh->resize(nPm,nPm,true);
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 elp.NonlinearElasticityUL::initElement(MNPC);
return p.NonlinearElasticityUL::initElement(MNPC);
}
@@ -473,38 +475,39 @@ bool ElasticityNormULMX::evalInt (LocalIntegral*& elmInt,
const FiniteElement& fe,
const TimeDomain& td, const Vec3& X) const
{
return static_cast<NonlinearElasticityULMX&>(problem).evalInt(elmInt,fe,td,X);
return static_cast<NonlinearElasticityULMX&>(myProblem).evalInt(elmInt,
fe,td,X);
}
bool ElasticityNormULMX::finalizeElement (LocalIntegral*& elmInt,
const TimeDomain& prm)
{
NonlinearElasticityULMX& elp = static_cast<NonlinearElasticityULMX&>(problem);
if (!elp.Hh) return false;
NonlinearElasticityULMX& p = static_cast<NonlinearElasticityULMX&>(myProblem);
if (!p.Hh) return false;
size_t iP;
#if INT_DEBUG > 0
std::cout <<"\n\n *** Entering ElasticityNormULMX::finalizeElement\n";
for (iP = 1; iP <= elp.myData.size(); iP++)
for (iP = 1; iP <= p.myData.size(); iP++)
{
const NonlinearElasticityULMX::ItgPtData& pt = elp.myData[iP-1];
const NonlinearElasticityULMX::ItgPtData& pt = p.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 ="<< *(elp.Hh) << std::endl;
std::cout <<"\n H ="<< *(p.Hh) << std::endl;
#endif
// 1. Eliminate the internal pressure DOFs by static condensation.
size_t nPM = elp.Hh->rows();
size_t nGP = elp.myData.size();
size_t nPM = p.Hh->rows();
size_t nGP = p.myData.size();
// Invert the H-matrix
Matrix& Hi = *elp.Hh;
Matrix& Hi = *p.Hh;
if (!utl::invert(Hi)) return false;
Vector Theta(nGP);
@@ -515,7 +518,7 @@ bool ElasticityNormULMX::finalizeElement (LocalIntegral*& elmInt,
double h1 = 0.0;
for (iP = 0; iP < nGP; iP++)
h1 += elp.myData[iP].F.det() * elp.myData[iP].detJW;
h1 += p.myData[iP].F.det() * p.myData[iP].detJW;
// All gauss point values are identical
Theta.front() = h1 * Hi(1,1);
@@ -529,7 +532,7 @@ bool ElasticityNormULMX::finalizeElement (LocalIntegral*& elmInt,
// Integrate the Ji-matrix
for (iP = 0; iP < nGP; iP++)
{
const NonlinearElasticityULMX::ItgPtData& pt = elp.myData[iP];
const NonlinearElasticityULMX::ItgPtData& pt = p.myData[iP];
for (size_t i = 0; i < nPM; i++)
Ji[i] += pt.Phi[i] * pt.F.det() * pt.detJW;
}
@@ -544,7 +547,7 @@ bool ElasticityNormULMX::finalizeElement (LocalIntegral*& elmInt,
// Calculate Theta = Hj*Phi
for (iP = 0; iP < nGP; iP++)
Theta[iP] = Hj.dot(elp.myData[iP].Phi);
Theta[iP] = Hj.dot(p.myData[iP].Phi);
}
#if INT_DEBUG > 1
@@ -553,7 +556,7 @@ bool ElasticityNormULMX::finalizeElement (LocalIntegral*& elmInt,
// 2. Evaluate the constitutive relation and integrate energy.
ElmNorm& pnorm = ElasticityNorm::getElmNormBuffer(elmInt);
ElmNorm& pnorm = NormBase::getElmNormBuffer(elmInt);
SymmTensor Sig(3), Eps(3);
for (iP = 0; iP < nGP; iP++)
@@ -561,7 +564,7 @@ bool ElasticityNormULMX::finalizeElement (LocalIntegral*& elmInt,
#if INT_DEBUG > 0
std::cout <<"\n iGP = "<< iP+1 << std::endl;
#endif
NonlinearElasticityULMX::ItgPtData& pt = elp.myData[iP];
NonlinearElasticityULMX::ItgPtData& pt = p.myData[iP];
Eps.rightCauchyGreen(pt.F); // Green Lagrange strain tensor
Eps -= 1.0;
Eps *= 0.5;
@@ -571,7 +574,7 @@ bool ElasticityNormULMX::finalizeElement (LocalIntegral*& elmInt,
// Compute the strain energy density, U(Eps) = Int_Eps (Sig:E) dE
double U = 0.0;
if (!elp.material->evaluate(elp.Cmat,Sig,U,pt.X,pt.F,Eps,3,&prm))
if (!p.material->evaluate(p.Cmat,Sig,U,pt.X,pt.F,Eps,3,&prm))
return false;
// Integrate energy norm a(u^h,u^h) = Int_Omega0 U(Eps) dV0

53
Apps/FiniteDefElasticity/NonlinearElasticityULMixed.C
View File

@@ -46,16 +46,17 @@ extern "C" {
enum SolutionVectors {
U = 0, // displacement
T = 1, // volumetric change (theta)
P = 2 // pressure
P = 2, // pressure
NSOL = 3
};
//! \brief Enum for element-level right-hand-side vectors.
enum ResidualVectors {
Ru = 3,
Rt = 4,
Rp = 5,
Rres = 6,
NVEC = 7
Ru = 0,
Rt = 1,
Rp = 2,
Rres = 3,
NVEC = 4
};
//! \brief Enum for element-level left-hand-side matrices.
@@ -198,13 +199,14 @@ NonlinearElasticityULMixed::NonlinearElasticityULMixed (unsigned short int n)
{
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 = &myMats->b[U];
eV = &mySols[U];
}
@@ -267,8 +269,8 @@ bool NonlinearElasticityULMixed::initElement (const std::vector<int>& MNPC1,
// 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(),myMats->b[T],nsd*n1,n1) +
utl::gather(MNPC2,1,2,primsol.front(),myMats->b[P],nsd*n1,n1)))
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;
@@ -291,7 +293,7 @@ bool NonlinearElasticityULMixed::initElementBou (const std::vector<int>& MNPC1,
// Extract the element level displacement vector
int ierr = 0;
if (!primsol.front().empty())
if ((ierr = utl::gather(MNPC1,nsd,primsol.front(),myMats->b[U])))
if ((ierr = utl::gather(MNPC1,nsd,primsol.front(),mySols[U])))
std::cerr <<" *** NonlinearElasticityULMixed::initElementBou: Detected "
<< ierr <<" node numbers out of range."<< std::endl;
@@ -319,11 +321,11 @@ bool NonlinearElasticityULMixed::evalIntMx (LocalIntegral*& elmInt,
bool lHaveStrains = !E.isZero(1.0e-16);
// Evaluate the volumetric change and pressure fields
double Theta = fe.N2.dot(myMats->b[T]) + 1.0;
double Press = fe.N2.dot(myMats->b[P]);
double Theta = fe.N2.dot(mySols[T]) + 1.0;
double Press = fe.N2.dot(mySols[P]);
#if INT_DEBUG > 0
std::cout <<"NonlinearElasticityULMixed::b_theta ="<< myMats->b[T];
std::cout <<"NonlinearElasticityULMixed::b_p ="<< myMats->b[P];
std::cout <<"NonlinearElasticityULMixed::b_theta ="<< mySols[T];
std::cout <<"NonlinearElasticityULMixed::b_p ="<< mySols[P];
std::cout <<"NonlinearElasticityULMixed::Theta = "<< Theta
<<", Press = "<< Press << std::endl;
#endif
@@ -526,32 +528,15 @@ NormBase* NonlinearElasticityULMixed::getNormIntegrand (AnaSol*) const
}
bool ElasticityNormULMixed::initElement (const std::vector<int>& MNPC1,
const std::vector<int>& MNPC2,
size_t n1)
{
return static_cast<NonlinearElasticityULMixed&>(problem).initElement(MNPC1,
MNPC2,
n1);
}
bool ElasticityNormULMixed::initElementBou (const std::vector<int>& MNPC1,
const std::vector<int>&, size_t)
{
return problem.initElementBou(MNPC1);
}
bool ElasticityNormULMixed::evalIntMx (LocalIntegral*& elmInt,
const MxFiniteElement& fe,
const TimeDomain& prm,
const Vec3& X) const
{
ElmNorm& pnorm = ElasticityNorm::getElmNormBuffer(elmInt);
ElmNorm& pnorm = NormBase::getElmNormBuffer(elmInt);
NonlinearElasticityULMixed* ulp;
ulp = static_cast<NonlinearElasticityULMixed*>(&problem);
ulp = static_cast<NonlinearElasticityULMixed*>(&myProblem);
// Evaluate the deformation gradient, F, and the Green-Lagrange strains, E
Tensor F(fe.dN1dX.cols());
@@ -560,7 +545,7 @@ bool ElasticityNormULMixed::evalIntMx (LocalIntegral*& elmInt,
return false;
// Evaluate the volumetric change field
double Theta = fe.N2.dot(ulp->myMats->b[T]) + 1.0;
double Theta = fe.N2.dot(ulp->mySols[T]) + 1.0;
// Compute the mixed model deformation gradient, F_bar
ulp->Fbar = F; // notice that F_bar always has dimension 3

13
Apps/FiniteDefElasticity/NonlinearElasticityULMixed.h
View File

@@ -122,19 +122,6 @@ public:
//! \brief Empty destructor.
virtual ~ElasticityNormULMixed() {}
//! \brief Initializes current element for numerical integration (mixed).
//! \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
virtual bool initElement(const std::vector<int>& MNPC1,
const std::vector<int>& MNPC2, size_t n1);
//! \brief Initializes current element for boundary integration (mixed).
//! \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
virtual bool initElementBou(const std::vector<int>& MNPC1,
const std::vector<int>& MNPC2, size_t n1);
//! \brief Evaluates the integrand at an interior point (mixed).
//! \param elmInt The local integral object to receive the contributions
//! \param[in] fe Mixed finite element data of current integration point

243
src/ASM/Integrand.h Normal file
View File

@@ -0,0 +1,243 @@
// $Id$
//==============================================================================
//!
//! \file Integrand.h
//!
//! \date Nov 11 2009
//!
//! \author Knut Morten Okstad / SINTEF
//!
//! \brief Abstract interface for classes representing FEM integrands.
//!
//==============================================================================
#ifndef _INTEGRAND_H
#define _INTEGRAND_H
#include "MatVec.h"
struct TimeDomain;
class LocalIntegral;
class FiniteElement;
class MxFiniteElement;
class Vec3;
/*!
\brief Abstract base class representing a system level integrated quantity.
\details This class defines the interface between the finite element (FE)
assembly drivers of the ASM-hierarchy and the problem-dependent classes
containing all physical properties for the problem to be solved.
The interface consists of methods for evaluating the integrand at interior
integration points (\a evalInt), and at boundary integration points
(\a evalBou). The latter are used for Neumann boundary conditions, typically.
The integrand evaluation methods have access to the FE basis function values
and derivatives through the FiniteElement argument. There are also a set
of methods dedicated for mixed field interpolation problems, which take
a MxFiniteElement object as argument instead.
*/
class Integrand
{
protected:
//! \brief The default constructor is protected to allow sub-classes only.
Integrand() {}
public:
//! \brief Empty destructor.
virtual ~Integrand() {}
// Global initialization interface
// ===============================
//! \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.
virtual void initIntegration(const TimeDomain&) = 0;
// Element-level initialization interface
// ======================================
//! \brief Initializes current element for numerical integration.
//! \param[in] MNPC Matrix of nodal point correspondance for current element
//! \param[in] X0 Cartesian coordinates of the element center
//! \param[in] nPt Number of integration points in this element
//!
//! \details This method is invoked once before starting the numerical
//! integration loop over the Gaussian quadrature points over an element.
//! It is supposed to perform all the necessary internal initializations
//! needed before the numerical integration is started for current element.
//! The default implementation forwards to an overloaded method not taking
//! \a X0 and \a nPt as arguments.
//! Reimplement this method for problems requiring the element center and/or
//! the number of integration points during/before the integrand evaluations.
virtual bool initElement(const std::vector<int>& MNPC,
const Vec3& X0, size_t nPt)
{
return this->initElement(MNPC);
}
//! \brief Initializes current element for numerical integration.
//! \param[in] MNPC Matrix of nodal point correspondance for current element
//!
//! \details This method is invoked once before starting the numerical
//! integration loop over the Gaussian quadrature points over an element.
//! It is supposed to perform all the necessary internal initializations
//! needed before the numerical integration is started for current element.
//! Reimplement this method for problems \e not requiring the
//! the element center nor the number of integration points before the
//! integration loop is started.
virtual bool initElement(const std::vector<int>& MNPC) = 0;
//! \brief Initializes current element for numerical integration (mixed).
//! \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
virtual bool initElement(const std::vector<int>& MNPC1,
const std::vector<int>& MNPC2, size_t n1) = 0;
//! \brief Initializes current element for boundary integration.
//! \param[in] MNPC Matrix of nodal point correspondance for current element
virtual bool initElementBou(const std::vector<int>& MNPC) = 0;
//! \brief Initializes current element for boundary integration (mixed).
//! \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
virtual bool initElementBou(const std::vector<int>& MNPC1,
const std::vector<int>& MNPC2, size_t n1) = 0;
// Integrand evaluation interface
// ==============================
//! \brief Defines which FE quantities are needed by the integrand.
//! \return 1: The basis functions and their gradients are needed
//! \return 2: As 1, but in addition the second derivatives of the basis
//! functions and the characteristic element size is needed
//! \return 3: As 1, but in addition the volume-averaged basis functions
//! are needed
//! \return 4: As 1, but in addition the element center coordinates are needed
virtual int getIntegrandType() const { return 1; }
//! \brief Evaluates the integrand at an interior point.
//! \param elmInt The local integral object to receive the contributions
//! \param[in] fe Finite element data of current integration point
//! \param[in] time Parameters for nonlinear and time-dependent simulations
//! \param[in] X Cartesian coordinates of current integration point
//!
//! \details The default implementation forwards to the stationary version.
//! Reimplement this method for time-dependent or non-linear problems.
virtual bool evalInt(LocalIntegral*& elmInt, const FiniteElement& fe,
const TimeDomain& time, const Vec3& X) const
{
return this->evalInt(elmInt,fe,X);
}
//! \brief Evaluates the integrand at an interior point.
//! \param elmInt The local integral object to receive the contributions
//! \param[in] fe Mixed finite element data of current integration point
//! \param[in] time Parameters for nonlinear and time-dependent simulations
//! \param[in] X Cartesian coordinates of current integration point
//!
//! \details This interface is used for mixed formulations only.
//! The default implementation forwards to the stationary version.
//! Reimplement this method for time-dependent or non-linear problems.
virtual bool evalIntMx(LocalIntegral*& elmInt, const MxFiniteElement& fe,
const TimeDomain& time, const Vec3& X) const
{
return this->evalIntMx(elmInt,fe,X);
}
//! \brief Finalizes the element quantities after the numerical integration.
//! \details This method is invoked once for each element, after the numerical
//! integration loop over interior points is finished and before the resulting
//! element quantities are assembled into their system level equivalents.
//! It can also be used to implement multiple integration point loops within
//! the same element, provided the necessary integration point values are
//! stored internally in the object during the first integration loop.
virtual bool finalizeElement(LocalIntegral*&, const TimeDomain&)
{
return true;
}
//! \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] time Parameters for nonlinear and time-dependent simulations
//! \param[in] X Cartesian coordinates of current integration point
//! \param[in] normal Boundary normal vector at current integration point
//!
//! \details The default implementation forwards to the stationary version.
//! Reimplement this method for time-dependent or non-linear problems.
virtual bool evalBou(LocalIntegral*& elmInt, const FiniteElement& fe,
const TimeDomain& time,
const Vec3& X, const Vec3& normal) const
{
return this->evalBou(elmInt,fe,X,normal);
}
//! \brief Evaluates the integrand at a boundary point.
//! \param elmInt The local integral object to receive the contributions
//! \param[in] fe Mixed finite element data of current integration point
//! \param[in] time Parameters for nonlinear and time-dependent simulations
//! \param[in] X Cartesian coordinates of current integration point
//! \param[in] normal Boundary normal vector at current integration point
//!
//! \details This interface is used for mixed formulations.
//! The default implementation forwards to the stationary version.
//! Reimplement this method for time-dependent or non-linear problems.
virtual bool evalBouMx(LocalIntegral*& elmInt, const MxFiniteElement& fe,
const TimeDomain& time,
const Vec3& X, const Vec3& normal) const
{
return this->evalBouMx(elmInt,fe,X,normal);
}
protected:
//! \brief Evaluates the integrand at interior points for stationary problems.
virtual bool evalInt(LocalIntegral*&, const FiniteElement& fe,
const Vec3&) const { return false; }
//! \brief Evaluates the integrand at interior points for stationary problems.
virtual bool evalIntMx(LocalIntegral*&, const MxFiniteElement& fe,
const Vec3&) const { return false; }
//! \brief Evaluates the integrand at boundary points for stationary problems.
virtual bool evalBou(LocalIntegral*&, const FiniteElement&,
const Vec3&, const Vec3&) const { return false; }
//! \brief Evaluates the integrand at boundary points for stationary problems.
virtual bool evalBouMx(LocalIntegral*&, const MxFiniteElement&,
const Vec3&, const Vec3&) const { return false; }
public:
// Solution field evaluation interface
// ===================================
//! \brief Evaluates the secondary solution at a result point.
//! \param[out] s The solution field values at current point
//! \param[in] N Basis function values at current point
//! \param[in] dNdX Basis function gradients at current point
//! \param[in] X Cartesian coordinates of current point
//! \param[in] MNPC Nodal point correspondance for the basis function values
virtual bool evalSol(Vector& s,
const Vector& N, const Matrix& dNdX,
const Vec3& X, const std::vector<int>& MNPC) const;
//! \brief Evaluates the secondary solution at a result point (mixed problem).
//! \param[out] s The solution field values at current point
//! \param[in] N1 Basis function values at current point, field 1
//! \param[in] N2 Basis function values at current point, field 2
//! \param[in] dN1dX Basis function gradients at current point, field 1
//! \param[in] dN2dX Basis function gradients at current point, field 2
//! \param[in] X Cartesian coordinates of current point
//! \param[in] MNPC1 Nodal point correspondance for the basis 1
//! \param[in] MNPC2 Nodal point correspondance for the basis 2
virtual bool evalSol(Vector& s,
const Vector& N1, const Vector& N2,
const Matrix& dN1dX, const Matrix& dN2dX, const Vec3& X,
const std::vector<int>& MNPC1,
const std::vector<int>& MNPC2) const;
};
#endif

178
src/ASM/IntegrandBase.C Normal file
View File

@@ -0,0 +1,178 @@
// $Id$
//==============================================================================
//!
//! \file IntegrandBase.C
//!
//! \date Nov 11 2009
//!
//! \author Knut Morten Okstad / SINTEF
//!
//! \brief Base classes representing FEM integrands.
//!
//==============================================================================
#include "IntegrandBase.h"
#include "ElmMats.h"
#include "ElmNorm.h"
#include "Utilities.h"
IntegrandBase::~IntegrandBase ()
{
if (myMats) delete myMats;
}
bool IntegrandBase::initElement (const std::vector<int>& MNPC)
{
if (myMats)
myMats->withLHS = true;
// Extract the primary solution vectors for this element
int ierr = 0;
for (size_t i = 0; i < mySols.size() && i < primsol.size() && ierr == 0; i++)
if (!primsol[i].empty())
ierr = utl::gather(MNPC,npv,primsol[i],mySols[i]);
if (ierr == 0) return true;
std::cerr <<" *** IntegrandBase::initElement: Detected "
<< ierr <<" node numbers out of range."<< std::endl;
return false;
}
bool IntegrandBase::initElement (const std::vector<int>&,
const std::vector<int>&, size_t)
{
std::cerr <<" *** IntegrandBase::initElement must be reimplemented"
<<" for mixed problems."<< std::endl;
return false;
}
bool IntegrandBase::initElementBou (const std::vector<int>& MNPC)
{
if (myMats)
myMats->withLHS = false;
// Extract current primary solution vector for this element
int ierr = 0;
if (!mySols.empty() && !primsol.empty() && !primsol.front().empty())
ierr = utl::gather(MNPC,npv,primsol.front(),mySols.front());
if (ierr == 0) return true;
std::cerr <<" *** IntegrandBase::initElementBou: Detected "
<< ierr <<" node numbers out of range."<< std::endl;
return false;
}
bool IntegrandBase::initElementBou(const std::vector<int>&,
const std::vector<int>&, size_t)
{
std::cerr <<" *** IntegrandBase::initElementBou must be reimplemented"
<<" for mixed problems."<< std::endl;
return false;
}
bool Integrand::evalSol (Vector&,
const Vector&, const Matrix&,
const Vec3&, const std::vector<int>&) const
{
std::cerr <<" *** Integrand::evalSol not implemented."<< std::endl;
return false;
}
bool Integrand::evalSol (Vector& s,
const Vector& N1, const Vector&,
const Matrix& dN1dX, const Matrix&, const Vec3& X,
const std::vector<int>& MNPC1,
const std::vector<int>&) const
{
return this->evalSol(s,N1,dN1dX,X,MNPC1);
}
bool IntegrandBase::evalPrimSol (Vector&, const VecFunc&, const Vec3&) const
{
std::cerr <<" *** IntegrandBase::evalPrimSol not implemented."<< std::endl;
return false;
}
bool IntegrandBase::evalSol (Vector&, const TensorFunc&, const Vec3&) const
{
std::cerr <<" *** IntegrandBase::evalSol not implemented."<< std::endl;
return false;
}
bool IntegrandBase::evalSol (Vector&, const STensorFunc&, const Vec3&) const
{
std::cerr <<" *** IntegrandBase::evalSol not implemented."<< std::endl;
return false;
}
bool IntegrandBase::evalPrimSol (double&, const RealFunc&, const Vec3&) const
{
std::cerr <<" *** IntegrandBase::evalPrimSol not implemented."<< std::endl;
return false;
}
bool IntegrandBase::evalSol (Vector&, const VecFunc&, const Vec3&) const
{
std::cerr <<" *** IntegrandBase::evalSol not implemented."<< std::endl;
return false;
}
void NormBase::initIntegration (const TimeDomain& time)
{
myProblem.initIntegration(time);
}
bool NormBase::initElement (const std::vector<int>& MNPC)
{
return myProblem.initElement(MNPC);
}
bool NormBase::initElement (const std::vector<int>& MNPC1,
const std::vector<int>& MNPC2, size_t n1)
{
return myProblem.initElement(MNPC1,MNPC2,n1);
}
bool NormBase::initElementBou (const std::vector<int>& MNPC)
{
return myProblem.initElementBou(MNPC);
}
bool NormBase::initElementBou (const std::vector<int>& MNPC1,
const std::vector<int>& MNPC2, size_t n1)
{
return myProblem.initElementBou(MNPC1,MNPC2,n1);
}
ElmNorm& NormBase::getElmNormBuffer (LocalIntegral*& elmInt, const size_t nn)
{
ElmNorm* eNorm = dynamic_cast<ElmNorm*>(elmInt);
if (eNorm) return *eNorm;
static double* data = 0;
if (!data) data = new double[nn];
static ElmNorm buf(data);
memset(data,0,nn*sizeof(double));
elmInt = &buf;
return buf;
}

285
src/ASM/IntegrandBase.h
View File

@@ -7,51 +7,37 @@
//!
//! \author Knut Morten Okstad / SINTEF
//!
//! \brief Abstract interface for classes representing FEM integrands.
//! \brief Base classes representing FEM integrands.
//!
//==============================================================================
#ifndef _INTEGRAND_BASE_H
#define _INTEGRAND_BASE_H
#include "Integrand.h"
#include "SIMenums.h"
#include "Function.h"
#include "MatVec.h"
struct TimeDomain;
class FiniteElement;
class MxFiniteElement;
class LocalIntegral;
class NormBase;
class ElmMats;
class ElmNorm;
class AnaSol;
class Vec3;
class VTF;
/*!
\brief Abstract base class representing a system level integrated quantity.
\details This class defines the interface between the finite element (FE)
assembly drivers of the ASM-hierarchy and the problem-dependent classes
containing all physical properties for the problem to be solved.
The interface consists of methods for evaluating the integrand at interior
integration points (\a evalInt), and at boundary integration points
(\a evalBou). The latter are used for Neumann boundary conditions, typically.
The integrand evaluation methods have access to the FE basis function values
and derivatives through the FiniteElement argument. There are also a set
of methods dedicated for mixed field interpolation problems, which take
and MxFiniteElement object as argument instead.
\brief Base class representing a system level integrated quantity.
*/
class Integrand
class IntegrandBase : public Integrand
{
protected:
//! \brief The default constructor is protected to allow sub-classes only.
Integrand() {}
IntegrandBase() : npv(0) {}
public:
//! \brief Empty destructor.
virtual ~Integrand() {}
//! \brief The destructor frees the dynamically allocated data objects.
virtual ~IntegrandBase();
//! \brief Prints out the problem definition to the given output stream.
virtual void print(std::ostream&) const {}
@@ -78,250 +64,67 @@ public:
//! \brief Initializes current element for numerical integration.
//! \param[in] MNPC Matrix of nodal point correspondance for current element
//! \param[in] X0 Cartesian coordinates of the element center
//! \param[in] nPt Number of integration points in this element
//!
//! \details This method is invoked once before starting the numerical
//! integration loop over the Gaussian quadrature points over an element.
//! It is supposed to perform all the necessary internal initializations
//! needed before the numerical integration is started for current element.
//! The default implementation forwards to an overloaded method not taking
//! \a X0 and \a nPt as arguments.
//! Reimplement this method for problems requiring the element center and/or
//! the number of integration points during/before the integrand evaluations.
virtual bool initElement(const std::vector<int>& MNPC,
const Vec3& X0, size_t nPt)
{
return this->initElement(MNPC);
}
//! \brief Initializes current element for numerical integration.
//! \param[in] MNPC Matrix of nodal point correspondance for current element
//!
//! \details Reimplement this method for problems \e not requiring the
//! the element center nor the number of integration points before the
//! integration loop is started.
virtual bool initElement(const std::vector<int>& MNPC)
{
std::cerr <<" *** Integrand::initElement not implemented."<< std::endl;
return false;
}
virtual bool initElement(const std::vector<int>& MNPC);
//! \brief Initializes current element for numerical integration (mixed).
//! \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
virtual bool initElement(const std::vector<int>& MNPC1,
const std::vector<int>& MNPC2, size_t n1)
{
std::cerr <<" *** Integrand::initElement not implemented."<< std::endl;
return false;
}
const std::vector<int>& MNPC2, size_t n1);
//! \brief Initializes current element for boundary integration.
//! \param[in] MNPC Matrix of nodal point correspondance for current element
virtual bool initElementBou(const std::vector<int>& MNPC)
{
std::cerr <<" *** Integrand::initElementBou not implemented."<< std::endl;
return false;
}
virtual bool initElementBou(const std::vector<int>& MNPC);
//! \brief Initializes current element for boundary integration (mixed).
//! \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
virtual bool initElementBou(const std::vector<int>& MNPC1,
const std::vector<int>& MNPC2, size_t n1)
{
std::cerr <<" *** Integrand::initElementBou not implemented."<< std::endl;
return false;
}
const std::vector<int>& MNPC2, size_t n1);
// Integrand evaluation interface
// ==============================
//! \brief Defines which FE quantities are needed by the integrand.
//! \return 1: The basis functions and their gradients are needed
//! \return 2: As 1, but in addition the second derivatives of the basis
//! functions and the characteristic element size is needed
//! \return 3: As 1, but in addition the volume-averaged basis functions
//! are needed
//! \return 4: As 1, but in addition the element center coordinates are needed
virtual int getIntegrandType() const { return 1; }
//! \brief Evaluates the integrand at an interior point.
//! \param elmInt The local integral object to receive the contributions
//! \param[in] fe Finite element data of current integration point
//! \param[in] time Parameters for nonlinear and time-dependent simulations
//! \param[in] X Cartesian coordinates of current integration point
//!
//! \details The default implementation forwards to the stationary version.
//! Reimplement this method for time-dependent or non-linear problems.
virtual bool evalInt(LocalIntegral*& elmInt, const FiniteElement& fe,
const TimeDomain& time, const Vec3& X) const
{
return this->evalInt(elmInt,fe,X);
}
//! \brief Evaluates the integrand at an interior point.
//! \param elmInt The local integral object to receive the contributions
//! \param[in] fe Mixed finite element data of current integration point
//! \param[in] time Parameters for nonlinear and time-dependent simulations
//! \param[in] X Cartesian coordinates of current integration point
//!
//! \details This interface is used for mixed formulations only.
//! The default implementation forwards to the stationary version.
//! Reimplement this method for time-dependent or non-linear problems.
virtual bool evalIntMx(LocalIntegral*& elmInt, const MxFiniteElement& fe,
const TimeDomain& time, const Vec3& X) const
{
return this->evalIntMx(elmInt,fe,X);
}
//! \brief Finalizes the element quantities after the numerical integration.
//! \details This method is invoked once for each element, after the numerical
//! integration loop over interior points is finished and before the resulting
//! element quantities are assembled into their system level equivalents.
//! It can also be used to implement multiple integration point loops within
//! the same element, provided the necessary integration point values are
//! stored internally in the object during the first integration loop.
virtual bool finalizeElement(LocalIntegral*&, const TimeDomain&)
{
return true;
}
//! \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] time Parameters for nonlinear and time-dependent simulations
//! \param[in] X Cartesian coordinates of current integration point
//! \param[in] normal Boundary normal vector at current integration point
//!
//! \details The default implementation forwards to the stationary version.
//! Reimplement this method for time-dependent or non-linear problems.
virtual bool evalBou(LocalIntegral*& elmInt, const FiniteElement& fe,
const TimeDomain& time,
const Vec3& X, const Vec3& normal) const
{
return this->evalBou(elmInt,fe,X,normal);
}
//! \brief Evaluates the integrand at a boundary point.
//! \param elmInt The local integral object to receive the contributions
//! \param[in] fe Mixed finite element data of current integration point
//! \param[in] time Parameters for nonlinear and time-dependent simulations
//! \param[in] X Cartesian coordinates of current integration point
//! \param[in] normal Boundary normal vector at current integration point
//!
//! \details This interface is used for mixed formulations.
//! The default implementation forwards to the stationary version.
//! Reimplement this method for time-dependent or non-linear problems.
virtual bool evalBouMx(LocalIntegral*& elmInt, const MxFiniteElement& fe,
const TimeDomain& time,
const Vec3& X, const Vec3& normal) const
{
return this->evalBouMx(elmInt,fe,X,normal);
}
protected:
//! \brief Evaluates the integrand at interior points for stationary problems.
virtual bool evalInt(LocalIntegral*&, const FiniteElement& fe,
const Vec3&) const { return false; }
//! \brief Evaluates the integrand at interior points for stationary problems.
virtual bool evalIntMx(LocalIntegral*&, const MxFiniteElement& fe,
const Vec3&) const { return false; }
//! \brief Evaluates the integrand at boundary points for stationary problems.
virtual bool evalBou(LocalIntegral*&, const FiniteElement&,
const Vec3&, const Vec3&) const { return false; }
//! \brief Evaluates the integrand at boundary points for stationary problems.
virtual bool evalBouMx(LocalIntegral*&, const MxFiniteElement&,
const Vec3&, const Vec3&) const { return false; }
public:
// Solution field evaluation interface
// ===================================
//! \brief Evaluates the secondary solution at a result point.
//! \param[out] s The solution field values at current point
//! \param[in] N Basis function values at current point
//! \param[in] dNdX Basis function gradients at current point
//! \param[in] X Cartesian coordinates of current point
//! \param[in] MNPC Nodal point correspondance for the basis function values
virtual bool evalSol(Vector& s,
const Vector& N, const Matrix& dNdX,
const Vec3& X, const std::vector<int>& MNPC) const
{
std::cerr <<" *** Integrand::evalSol not implemented"<< std::endl;
return false;
}
//! \brief Evaluates the secondary solution at a result point (mixed problem).
//! \param[out] s The solution field values at current point
//! \param[in] N1 Basis function values at current point, field 1
//! \param[in] N2 Basis function values at current point, field 2
//! \param[in] dN1dX Basis function gradients at current point, field 1
//! \param[in] dN2dX Basis function gradients at current point, field 2
//! \param[in] X Cartesian coordinates of current point
//! \param[in] MNPC1 Nodal point correspondance for the basis 1
//! \param[in] MNPC2 Nodal point correspondance for the basis 2
virtual bool evalSol(Vector& s,
const Vector& N1, const Vector& N2,
const Matrix& dN1dX, const Matrix& dN2dX, const Vec3& X,
const std::vector<int>& MNPC1,
const std::vector<int>& MNPC2) const
{
return this->evalSol(s,N1,dN1dX,X,MNPC1);
}
//! \brief Evaluates the analytical primary solution at a result point.
//! \param[out] s The solution field values at current point
//! \param[in] asol The analytical solution field (vector field)
//! \param[in] X Cartesian coordinates of current point
virtual bool evalPrimSol(Vector& s, const VecFunc& asol, const Vec3& X) const
{
std::cerr <<" *** Integrand::evalPrimSol not implemented"<< std::endl;
return false;
}
virtual bool evalPrimSol(Vector& s,
const VecFunc& asol, const Vec3& X) const;
//! \brief Evaluates the analytical secondary solution at a result point.
//! \param[out] s The solution field values at current point
//! \param[in] asol The analytical solution field (tensor field)
//! \param[in] X Cartesian coordinates of current point
virtual bool evalSol(Vector& s, const TensorFunc& asol, const Vec3& X) const
{
std::cerr <<" *** Integrand::evalSol (exact) not implemented"<< std::endl;
return false;
}
virtual bool evalSol(Vector& s,
const TensorFunc& asol, const Vec3& X) const;
//! \brief Evaluates the analytical secondary solution at a result point.
//! \param[out] s The solution field values at current point
//! \param[in] asol The analytical solution field (symmetric tensor field)
//! \param[in] X Cartesian coordinates of current point
virtual bool evalSol(Vector& s, const STensorFunc& asol, const Vec3& X) const
{
std::cerr <<" *** Integrand::evalSol (exact) not implemented"<< std::endl;
return false;
}
virtual bool evalSol(Vector& s,
const STensorFunc& asol, const Vec3& X) const;
//! \brief Evaluates the analytical primary solution at a result point.
//! \param[out] s The solution field value at current point
//! \param[in] asol The analytical solution field (scalar field)
//! \param[in] X Cartesian coordinates of current point
virtual bool evalPrimSol(double& s, const RealFunc& asol, const Vec3& X) const
{
std::cerr <<" *** Integrand::evalPrimSol not implemented"<< std::endl;
return false;
}
virtual bool evalPrimSol(double& s,
const RealFunc& asol, const Vec3& X) const;
//! \brief Evaluates the analytical secondary solution at a result point.
//! \param[out] s The solution field values at current point
//! \param[in] asol The analytical solution field (vector field)
//! \param[in] X Cartesian coordinates of current point
virtual bool evalSol(Vector& s, const VecFunc& asol, const Vec3& X) const
{
std::cerr <<" *** Integrand::evalSol (exact) not implemented"<< std::endl;
return false;
}
virtual bool evalSol(Vector& s,
const VecFunc& asol, const Vec3& X) const;
// Various service methods
@@ -353,29 +156,65 @@ public:
Vector& getSolution(size_t n = 0) { return primsol[n]; }
protected:
Vectors primsol; //!< Primary solution vectors for current patch
Vectors primsol; //!< Primary solution vectors for current patch
ElmMats* myMats; //!< Local element matrices
Vectors mySols; //!< Local element solution vectors
unsigned short int npv; //!< Number of primary solution variables per node
};
/*!
\brief Abstract base class representing a system level norm quantity.
\brief Base class representing a system level norm quantity.
*/
class NormBase : public Integrand
{
protected:
//! \brief The default constructor is protected to allow sub-classes only.
NormBase() {}
NormBase(IntegrandBase& p) : myProblem(p) {}
public:
//! \brief Empty destructor.
virtual ~NormBase() {}
//! \brief Initializes the integrand for a new integration loop.
virtual void initIntegration(const TimeDomain& time);
//! \brief Initializes current element for numerical integration.
virtual bool initElement(const std::vector<int>& MNPC);
//! \brief Initializes current element for numerical integration (mixed).
virtual bool initElement(const std::vector<int>& MNPC1,
const std::vector<int>& MNPC2, size_t n1);
//! \brief Initializes current element for boundary integration.
virtual bool initElementBou(const std::vector<int>& MNPC);
//! \brief Initializes current element for boundary integration (mixed).
virtual bool initElementBou(const std::vector<int>& MNPC1,
const std::vector<int>& MNPC2, size_t n1);
//! \brief Returns whether this norm has explicit boundary contributions.
virtual bool hasBoundaryTerms() const { return false; }
//! \brief Returns a 1-based index of the external energy norm.
virtual size_t indExt() const { return 0; }
//! \brief Returns the number of field components.
virtual size_t getNoFields() const { return 0; }
protected:
//! \brief Returns the element norm object to use in the integration.
//! \param elmInt The local integral object to receive norm contributions
//! \param[in] nn Size of static norm buffer
//!
//! \details If \a elmInt is NULL or cannot be casted to an ElmNorm pointer,
//! a local static buffer is used instead.
static ElmNorm& getElmNormBuffer(LocalIntegral*& elmInt, const size_t nn = 4);
protected:
IntegrandBase& myProblem; //!< The problem-specific data
};
#endif

83
src/Integrands/Elasticity.C
View File

@@ -25,6 +25,8 @@
Elasticity::Elasticity (unsigned short int n) : nsd(n)
{
npv = n; // Number of primary unknowns per node
// Default is zero gravity
grav[0] = grav[1] = grav[2] = 0.0;
@@ -41,7 +43,6 @@ Elasticity::Elasticity (unsigned short int n) : nsd(n)
Elasticity::~Elasticity ()
{
if (myMats) delete myMats;
if (locSys) delete locSys;
}
@@ -94,10 +95,11 @@ void Elasticity::setMode (SIM::SolutionMode mode)
break;
case SIM::BUCKLING:
myMats->resize(2,1);
myMats->resize(2,0);
mySols.resize(1);
eKm = &myMats->A[0];
eKg = &myMats->A[1];
eV = &myMats->b[0];
eV = &mySols[0];
break;
case SIM::STIFF_ONLY:
@@ -107,7 +109,7 @@ void Elasticity::setMode (SIM::SolutionMode mode)
case SIM::MASS_ONLY:
myMats->resize(1,0);
eM = &myMats->A[0];
eM = &myMats->A[0];
break;
case SIM::RHS_ONLY:
@@ -120,13 +122,13 @@ void Elasticity::setMode (SIM::SolutionMode mode)
case SIM::RECOVERY:
myMats->rhsOnly = true;
if (myMats->b.empty())
myMats->b.resize(1);
eV = &myMats->b[0];
mySols.resize(1);
eV = &mySols[0];
break;
default:
myMats->resize(0,0);
mySols.clear();
tracVal.clear();
}
}
@@ -160,23 +162,7 @@ bool Elasticity::initElement (const std::vector<int>& MNPC)
if (eM) eM->resize(nsd*nen,nsd*nen,true);
if (eS) eS->resize(nsd*nen,true);
// Extract the current primary solution and store in the array *eV
int ierr = 0;
if (eV && !primsol.front().empty())
ierr = utl::gather(MNPC,nsd,primsol.front(),*eV);
// If previous solutions are required, they are stored in the
// (primsol.size()-1) last entries in myMats->b
int j = 1+myMats->b.size()-primsol.size();
for (size_t i = 1; i < primsol.size() && ierr == 0; i++, j++)
if (!primsol[i].empty())
ierr = utl::gather(MNPC,nsd,primsol[i],myMats->b[j]);
if (ierr == 0) return true;
std::cerr <<" *** Elasticity::initElement: Detected "
<< ierr <<" node numbers out of range."<< std::endl;
return false;
return this->IntegrandBase::initElement(MNPC);
}
@@ -187,15 +173,7 @@ bool Elasticity::initElementBou (const std::vector<int>& MNPC)
if (eS) eS->resize(nsd*MNPC.size(),true);
int ierr = 0;
if (eV && !primsol.front().empty())
ierr = utl::gather(MNPC,nsd,primsol.front(),*eV);
if (ierr == 0) return true;
std::cerr <<" *** Elasticity::initElement: Detected "
<< ierr <<" node numbers out of range."<< std::endl;
return false;
return this->IntegrandBase::initElementBou(MNPC);
}
@@ -572,47 +550,17 @@ NormBase* Elasticity::getNormIntegrand (AnaSol* asol) const
}
void ElasticityNorm::initIntegration (const TimeDomain& time)
{
problem.initIntegration(time);
}
bool ElasticityNorm::initElement (const std::vector<int>& MNPC)
{
return problem.initElement(MNPC);
}
bool ElasticityNorm::initElementBou (const std::vector<int>& MNPC)
{
return problem.initElementBou(MNPC);
}
size_t ElasticityNorm::getNoFields (int) const
size_t ElasticityNorm::getNoFields () const
{
return anasol ? 4 : 2;
}
ElmNorm& ElasticityNorm::getElmNormBuffer (LocalIntegral*& elmInt)
{
ElmNorm* eNorm = dynamic_cast<ElmNorm*>(elmInt);
if (eNorm) return *eNorm;
static double data[4];
static ElmNorm buf(data);
memset(data,0,4*sizeof(double));
elmInt = &buf;
return buf;
}
bool ElasticityNorm::evalInt (LocalIntegral*& elmInt, const FiniteElement& fe,
const Vec3& X) const
{
ElmNorm& pnorm = ElasticityNorm::getElmNormBuffer(elmInt);
Elasticity& problem = static_cast<Elasticity&>(myProblem);
ElmNorm& pnorm = NormBase::getElmNormBuffer(elmInt);
// Evaluate the inverse constitutive matrix at this point
Matrix Cinv;
@@ -663,9 +611,10 @@ bool ElasticityNorm::evalInt (LocalIntegral*& elmInt, const FiniteElement& fe,
bool ElasticityNorm::evalBou (LocalIntegral*& elmInt, const FiniteElement& fe,
const Vec3& X, const Vec3& normal) const
{
Elasticity& problem = static_cast<Elasticity&>(myProblem);
if (!problem.haveLoads()) return true;
ElmNorm& pnorm = ElasticityNorm::getElmNormBuffer(elmInt);
ElmNorm& pnorm = NormBase::getElmNormBuffer(elmInt);
// Evaluate the surface traction
Vec3 T = problem.getTraction(X,normal);

31
src/Integrands/Elasticity.h
View File

@@ -20,9 +20,6 @@
class LocalSystem;
class Material;
class ElmMats;
class ElmNorm;
class VTF;
/*!
@@ -32,7 +29,7 @@ class VTF;
Thus, it is regarded as an abstract base class with a protected constructor.
*/
class Elasticity : public Integrand
class Elasticity : public IntegrandBase
{
protected:
//! \brief The default constructor initializes all pointers to zero.
@@ -199,8 +196,6 @@ protected:
Vector* iS; //!< Pointer to element internal force vector
Vector* eV; //!< Pointer to element displacement vector
ElmMats* myMats; //!< Local element matrices, result of numerical integration
// Physical properties
Material* material; //!< Material data and constitutive relation
double grav[3]; //!< Gravitation vector
@@ -230,7 +225,7 @@ public:
//! \brief The only constructor initializes its data members.
//! \param[in] p The linear elasticity problem to evaluate norms for
//! \param[in] a The analytical stress field (optional)
ElasticityNorm(Elasticity& p, STensorFunc* a = 0) : problem(p), anasol(a) {}
ElasticityNorm(Elasticity& p, STensorFunc* a = 0) : NormBase(p), anasol(a) {}
//! \brief Empty destructor.
virtual ~ElasticityNorm() {}
@@ -239,16 +234,6 @@ public:
//! \brief Returns a 1-based index of the external energy norm.
virtual size_t indExt() const { return 2; }
//! \brief Initializes the integrand for a new integration loop.
//! \param[in] time Parameters for nonlinear and time-dependent simulations
virtual void initIntegration(const TimeDomain& time);
//! \brief Initializes current element for numerical integration.
//! \param[in] MNPC Matrix of nodal point correspondance for current element
virtual bool initElement(const std::vector<int>& MNPC);
//! \brief Initializes current element for boundary numerical integration.
//! \param[in] MNPC Matrix of nodal point correspondance for current element
virtual bool initElementBou(const std::vector<int>& MNPC);
//! \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
@@ -265,17 +250,7 @@ public:
const Vec3& X, const Vec3& normal) const;
//! \brief Returns the number of norm quantities.
virtual size_t getNoFields(int = 0) const;
protected:
//! \brief Returns the element norm object to use in the integration.
//! \param elmInt The local integral object to receive norm contributions
//!
//! \details If \a elmInt is NULL or cannot be casted to an ElmNorm pointer,
//! a local static buffer is used instead.
static ElmNorm& getElmNormBuffer(LocalIntegral*& elmInt);
Elasticity& problem; //!< The problem-specific data
virtual size_t getNoFields() const;
private:
STensorFunc* anasol; //!< Analytical stress field

111
src/Integrands/Poisson.C
View File

@@ -14,6 +14,7 @@
#include "Poisson.h"
#include "FiniteElement.h"
#include "Utilities.h"
#include "ElmMats.h"
#include "ElmNorm.h"
#include "Tensor.h"
#include "Vec3Oper.h"
@@ -23,6 +24,8 @@
Poisson::Poisson (unsigned short int n) : nsd(n)
{
npv = 1; // One primary unknown per node (scalar equation)
kappa = 1.0;
tracFld = 0;
@@ -32,48 +35,49 @@ Poisson::Poisson (unsigned short int n) : nsd(n)
// Only the current solution is needed
primsol.resize(1);
mySols.resize(1);
myMats = new ElmMats();
}
void Poisson::setMode (SIM::SolutionMode mode)
{
myMats.rhsOnly = false;
myMats->rhsOnly = false;
eM = 0;
eS = eV = 0;
switch (mode)
{
case SIM::STATIC:
myMats.A.resize(1);
myMats.b.resize(1);
eM = &myMats.A[0];
eS = &myMats.b[0];
myMats->A.resize(1);
myMats->b.resize(1);
eM = &myMats->A[0];
eS = &myMats->b[0];
tracVal.clear();
break;
case SIM::VIBRATION:
myMats.A.resize(1);
eM = &myMats.A[0];
myMats->A.resize(1);
eM = &myMats->A[0];
break;
case SIM::RHS_ONLY:
myMats.rhsOnly = true;
if (myMats.b.empty())
myMats.b.resize(1);
eS = &myMats.b[0];
myMats->rhsOnly = true;
if (myMats->b.empty())
myMats->b.resize(1);
eS = &myMats->b[0];
tracVal.clear();
break;
case SIM::RECOVERY:
myMats.rhsOnly = true;
if (myMats.b.empty())
myMats.b.resize(1);
eV = &myMats.b[0];
myMats->rhsOnly = true;
eV = &mySols[0];
break;
default:
myMats.A.clear();
myMats.b.clear();
myMats->A.clear();
myMats->b.clear();
tracVal.clear();
}
}
@@ -87,6 +91,12 @@ double Poisson::getTraction (const Vec3& X, const Vec3& n) const
}
void Poisson::setTraction (VecFunc* tf)
{
tracFld = tf;
myMats->rhsOnly = true;
}
bool Poisson::initElement (const std::vector<int>& MNPC)
{
@@ -95,14 +105,7 @@ bool Poisson::initElement (const std::vector<int>& MNPC)
if (eM) eM->resize(nen,nen,true);
if (eS) eS->resize(nen,true);
int ierr = 0;
if (eV && !primsol.front().empty())
if ((ierr = utl::gather(MNPC,1,primsol.front(),*eV)))
std::cerr <<" *** Poisson::initElement: Detected "
<< ierr <<" node numbers out of range."<< std::endl;
myMats.withLHS = true;
return ierr == 0;
return this->IntegrandBase::initElement(MNPC);
}
@@ -110,7 +113,7 @@ bool Poisson::initElementBou (const std::vector<int>& MNPC)
{
if (eS) eS->resize(MNPC.size(),true);
myMats.withLHS = false;
myMats->withLHS = false;
return true;
}
@@ -118,7 +121,7 @@ bool Poisson::initElementBou (const std::vector<int>& MNPC)
bool Poisson::evalInt (LocalIntegral*& elmInt, const FiniteElement& fe,
const Vec3& X) const
{
elmInt = &myMats;
elmInt = myMats;
if (eM)
{
@@ -140,7 +143,7 @@ bool Poisson::evalInt (LocalIntegral*& elmInt, const FiniteElement& fe,
bool Poisson::evalBou (LocalIntegral*& elmInt, const FiniteElement& fe,
const Vec3& X, const Vec3& normal) const
{
elmInt = &myMats;
elmInt = myMats;
if (!tracFld)
{
std::cerr <<" *** Poisson::evalBou: No tractions."<< std::endl;
@@ -183,15 +186,6 @@ bool Poisson::writeGlvT (VTF* vtf, int iStep, int& nBlock) const
return vtf->writeVblk(nBlock,"Tractions",1,iStep);
}
double Poisson::evalSol (const Vector& N) const
{
double u = 0;
if (eV && eV->size() == N.size()*nsd)
u = eV->dot(N);
return u;
}
bool Poisson::formCmatrix (Matrix& C, const Vec3&, bool invers) const
{
@@ -200,6 +194,12 @@ bool Poisson::formCmatrix (Matrix& C, const Vec3&, bool invers) const
}
double Poisson::evalSol (const Vector& N) const
{
return eV ? eV->dot(N) : 0.0;
}
bool Poisson::evalSol (Vector& q, const Vector&,
const Matrix& dNdX, const Vec3& X,
const std::vector<int>& MNPC) const
@@ -255,18 +255,6 @@ bool Poisson::evalSol (Vector& q, const Matrix& dNdX, const Vec3& X) const
return true;
}
ElmNorm* PoissonNorm::getElmNormBuffer (LocalIntegral*& elmInt)
{
ElmNorm* eNorm = dynamic_cast<ElmNorm*>(elmInt);
if (eNorm) return eNorm;
static double data[4];
static ElmNorm buf(data);
memset(data,0,4*sizeof(double));
elmInt = eNorm = &buf;
return eNorm;
}
bool Poisson::evalSol (Vector& s, const VecFunc& asol, const Vec3& X) const
{
@@ -310,21 +298,11 @@ NormBase* Poisson::getNormIntegrand (AnaSol* asol) const
}
bool PoissonNorm::initElement (const std::vector<int>& MNPC)
{
return problem.initElement(MNPC);
}
bool PoissonNorm::evalInt (LocalIntegral*& elmInt, const FiniteElement& fe,
const Vec3& X) const
{
ElmNorm* eNorm = dynamic_cast<ElmNorm*>(elmInt);
if (!eNorm)
{
std::cerr <<" *** PoissonNorm::evalInt: No norm array."<< std::endl;
return false;
}
Poisson& problem = static_cast<Poisson&>(myProblem);
ElmNorm& pnorm = NormBase::getElmNormBuffer(elmInt);
// Evaluate the inverse constitutive matrix at this point
Matrix Cinv;
@@ -335,7 +313,6 @@ bool PoissonNorm::evalInt (LocalIntegral*& elmInt, const FiniteElement& fe,
if (!problem.evalSol(sigmah,fe.dNdX,X)) return false;
// Integrate the energy norm a(u^h,u^h)
ElmNorm& pnorm = *eNorm;
pnorm[0] += sigmah.dot(Cinv*sigmah)*fe.detJxW;
if (anasol)
{
@@ -351,10 +328,12 @@ bool PoissonNorm::evalInt (LocalIntegral*& elmInt, const FiniteElement& fe,
return true;
}
bool PoissonNorm::evalBou(LocalIntegral*& elmInt, const FiniteElement& fe,
const Vec3& X, const Vec3& normal) const
bool PoissonNorm::evalBou (LocalIntegral*& elmInt, const FiniteElement& fe,
const Vec3& X, const Vec3& normal) const
{
ElmNorm* eNorm = PoissonNorm::getElmNormBuffer(elmInt);
Poisson& problem = static_cast<Poisson&>(myProblem);
ElmNorm& pnorm = NormBase::getElmNormBuffer(elmInt);
// Evaluate the surface flux
double T = problem.getTraction(X,normal);
@@ -362,6 +341,6 @@ bool PoissonNorm::evalBou(LocalIntegral*& elmInt, const FiniteElement& fe,
double u = problem.evalSol(fe.N);
// Integrate the external energy
(*eNorm)[1] += T*u*fe.detJxW;
pnorm[1] += T*u*fe.detJxW;
return true;
}

40
src/Integrands/Poisson.h
View File

@@ -15,21 +15,18 @@
#define _POISSON_H
#include "IntegrandBase.h"
#include "ElmMats.h"
#include "Vec3.h"
#include <map>
class VTF;
class ElmNorm;
/*!
\brief Class representing the integrand of the Poisson problem.
\details This class supports constant isotropic constitutive properties only.
Properties with spatial variation has to be implemented as sub-classes.
*/
class Poisson : public Integrand
class Poisson : public IntegrandBase
{
public:
//! \brief The default constructor initializes all pointers to zero.
@@ -43,7 +40,7 @@ public:
virtual void setMode(SIM::SolutionMode mode);
//! \brief Defines the traction field to use in Neumann boundary conditions.
void setTraction(VecFunc* tf) { tracFld = tf; myMats.rhsOnly = true; }
void setTraction(VecFunc* tf);
//! \brief Clears the integration point traction values.
void clearTracVal() { tracVal.clear(); }
@@ -78,6 +75,11 @@ public:
virtual bool evalBou(LocalIntegral*& elmInt, const FiniteElement& fe,
const Vec3& X, const Vec3& normal) const;
//! \brief Evaluates the primary solution at a result point.
//! \param[in] N Basis function values at current point
//! \return Primary solution value at current point
double evalSol(const Vector& N) const;
//! \brief Evaluates the secondary solution at a result point.
//! \param[out] s Array of solution field values at current point
//! \param[in] N Basis function values at current point
@@ -88,11 +90,6 @@ public:
const Vector& N, const Matrix& dNdX,
const Vec3& X, const std::vector<int>& MNPC) const;
//! \brief Evaluates the primary solution at a result point.
//! \param[in] N Basis function values at current point
//! \return Primary solution vector at current point
double evalSol(const Vector& N) const;
//! \brief Evaluates the finite element (FE) solution at an integration point.
//! \param[out] s The FE solution values at current point
//! \param[in] dNdX Basis function gradients at current point
@@ -148,8 +145,6 @@ protected:
Vector* eS; //!< Pointer to element right-hand-side vector
Vector* eV; //!< Pointer to element solution vector
mutable ElmMats myMats; //!< Local element matrices
VecFunc* tracFld; //!< Pointer to boundary traction field
RealFunc* heatSrc; //!< Pointer to interior heat source
@@ -174,14 +169,10 @@ public:
//! \brief The only constructor initializes its data members.
//! \param[in] p The Poisson problem to evaluate norms for
//! \param[in] a The analytical heat flux (optional)
PoissonNorm(Poisson& p, VecFunc* a = 0) : problem(p), anasol(a) {}
PoissonNorm(Poisson& p, VecFunc* a = 0) : NormBase(p), anasol(a) {}
//! \brief Empty destructor.
virtual ~PoissonNorm() {}
//! \brief Initializes current element for numerical integration.
//! \param[in] MNPC Matrix of nodal point correspondance for current element
virtual bool initElement(const std::vector<int>& MNPC);
//! \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
@@ -189,10 +180,6 @@ public:
virtual bool evalInt(LocalIntegral*& elmInt, const FiniteElement& fe,
const Vec3& X) const;
//! \brief Returns the number of norm quantities.
virtual size_t getNoFields(int = 0) const { return anasol ? 3 : 1; }
protected:
//! \brief Evaluates the integrand at a boundary point.
//! \param elmInt The local integral object to receive the contributions
//! \param[fe] Finite Element quantities
@@ -201,16 +188,11 @@ protected:
virtual bool evalBou(LocalIntegral*& elmInt, const FiniteElement& fe,
const Vec3& X, const Vec3& normal) const;
//! \brief Get a pointer to the element norm object to use in the integration.
//! \param elmInt The local integral object to receive norm contributions
//!
//! \details If \a elmInt is NULL or cannot be casted to a ElmNorm pointer,
//! a local static buffer is used instead.
static ElmNorm* getElmNormBuffer(LocalIntegral*& elmInt);
//! \brief Returns the number of norm quantities.
virtual size_t getNoFields() const { return anasol ? 3 : 1; }
private:
Poisson& problem; //!< The problem-specific data
VecFunc* anasol; //!< Analytical heat flux
VecFunc* anasol; //!< Analytical heat flux
};
#endif

22
src/SIM/SIMbase.h
View File

@@ -23,7 +23,7 @@
#include <map>
class ASMbase;
class Integrand;
class IntegrandBase;
class AnaSol;
class VTF;
class SAMpatch;
@@ -123,7 +123,7 @@ public:
void printProblem(std::ostream& os) const;
//! \brief Returns a pointer to the problem-specific data object.
const Integrand* getProblem() const { return myProblem; }
const IntegrandBase* getProblem() const { return myProblem; }
virtual std::string getName() const { return "SIMbase"; }
@@ -455,15 +455,15 @@ protected:
typedef std::map<int,TractionFunc*> TracFuncMap;
// Model attributes
FEModelVec myModel; //!< The actual NURBS/spline model
PropertyVec myProps; //!< Physical property mapping
SclFuncMap myScalars; //!< Scalar property fields
VecFuncMap myVectors; //!< Vector property fields
TracFuncMap myTracs; //!< Traction property fields
Integrand* myProblem; //!< Problem-specific data and methods
AnaSol* mySol; //!< Analytical/Exact solution
VTF* myVtf; //!< VTF-file for result visualization
ResPointVec myPoints; //!< User-defined result sampling points
FEModelVec myModel; //!< The actual NURBS/spline model
PropertyVec myProps; //!< Physical property mapping
SclFuncMap myScalars; //!< Scalar property fields
VecFuncMap myVectors; //!< Vector property fields
TracFuncMap myTracs; //!< Traction property fields
IntegrandBase* myProblem; //!< Problem-specific data and methods
AnaSol* mySol; //!< Analytical/Exact solution
VTF* myVtf; //!< VTF-file for result visualization
ResPointVec myPoints; //!< User-defined result sampling points
// Parallel computing attributes
int nGlPatches; //!< Number of global patches

6
src/Utility/HDF5Writer.C
View File

@@ -186,7 +186,7 @@ bool HDF5Writer::readSIM (int level, const DataEntry& entry)
SIMbase* sim = static_cast<SIMbase*>(entry.second.data);
Vector* sol = static_cast<Vector*>(entry.second.data2);
if (!sol) return false;
const Integrand* prob = sim->getProblem();
const IntegrandBase* prob = sim->getProblem();
for (int i = 0; i < sim->getNoPatches() && ok; ++i) {
std::stringstream str;
@@ -267,7 +267,7 @@ void HDF5Writer::writeSIM (int level, const DataEntry& entry)
Vector* sol = static_cast<Vector*>(entry.second.data2);
if (!sol) return;
const Integrand* prob = sim->getProblem();
const IntegrandBase* prob = sim->getProblem();
if (level == 0) { // TODO time dependent geometries
writeBasis(sim,sim->getName()+"-1",1,level);
if (prob->mixedFormulation())
@@ -288,7 +288,7 @@ void HDF5Writer::writeSIM (int level, const DataEntry& entry)
if (loc > 0) // we own the patch
{
size_t ndof1 = sim->extractPatchSolution(*sol,loc-1);
Vector& psol = const_cast<Integrand*>(prob)->getSolution();
Vector& psol = const_cast<IntegrandBase*>(prob)->getSolution();
if (prob->mixedFormulation())
{
// Mixed methods: The primary solution vector is referring to two bases

4
src/Utility/XMLWriter.C
View File

@@ -127,8 +127,8 @@ void XMLWriter::writeSIM (int level, const DataEntry& entry)
if (m_rank != 0)
return;
const SIMbase* sim = static_cast<const SIMbase*>(entry.second.data);
const Integrand* prob = sim->getProblem();
const SIMbase* sim = static_cast<const SIMbase*>(entry.second.data);
const IntegrandBase* prob = sim->getProblem();
std::string g2file;
if (prob->mixedFormulation())