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

changed: rip ElmMats and finite element vectors out of integrand

Browse Source 
this is necessary to allow multi-threaded assembly

git-svn-id: http://svn.sintef.no/trondheim/IFEM/trunk@1402 e10b68d5-8a6e-419e-a041-bce267b0401d
This commit is contained in:
akva 2012-01-18 13:43:28 +00:00 committed by Knut Morten Okstad
parent c6a3a640b2
commit 5e9a5d88e7
44 changed files with 660 additions and 760 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

8
src/ASM/ASMbase.h
View File

@ -256,20 +256,16 @@ public:
//! \param integrand Object with problem-specific data and methods
//! \param glbInt The integrated quantity
//! \param[in] time Parameters for nonlinear/time-dependent simulations
//! \param locInt Vector of element-wise contributions to \a glbInt
virtual bool integrate(Integrand& integrand,
GlobalIntegral& glbInt, const TimeDomain& time,
const LintegralVec& locInt = LintegralVec()) = 0;
GlobalIntegral& glbInt, const TimeDomain& time) = 0;
//! \brief Evaluates a boundary integral over a patch face/edge.
//! \param integrand Object with problem-specific data and methods
//! \param[in] lIndex Local index of the boundary face/edge
//! \param glbInt The integrated quantity
//! \param[in] time Parameters for nonlinear/time-dependent simulations
//! \param locInt Vector of element-wise contributions to \a glbInt
virtual bool integrate(Integrand& integrand, int lIndex,
GlobalIntegral& glbInt, const TimeDomain& time,
const LintegralVec& locInt = LintegralVec()) = 0;
GlobalIntegral& glbInt, const TimeDomain& time) = 0;
//! \brief Evaluates a boundary integral over a patch edge.
//! \param integrand Object with problem-specific data and methods

48
src/ASM/ASMs1D.C
View File

@ -19,6 +19,7 @@
#include "TimeDomain.h"
#include "FiniteElement.h"
#include "GlobalIntegral.h"
#include "LocalIntegral.h"
#include "IntegrandBase.h"
#include "CoordinateMapping.h"
#include "GaussQuadrature.h"
@ -112,7 +113,6 @@ void ASMs1D::clear (bool retainGeometry)
}
bool ASMs1D::refine (const RealArray& xi)
{
if (!curv || xi.empty()) return false;
@ -514,8 +514,7 @@ void ASMs1D::extractBasis (double u, Vector& N, Matrix& dNdu,
bool ASMs1D::integrate (Integrand& integrand,
GlobalIntegral& glInt,
const TimeDomain& time,
const LintegralVec& locInt)
const TimeDomain& time)
{
if (!curv) return true; // silently ignore empty patches
@ -563,14 +562,8 @@ bool ASMs1D::integrate (Integrand& integrand,
if (!this->getElementCoordinates(Xnod,iel)) return false;
// Initialize element matrices
if (!integrand.initElement(MNPC[iel-1],X,nRed)) return false;
// Caution: Unless locInt is empty, we assume it points to an array of
// LocalIntegral pointers, of length at least the number of elements in
// the model (as defined by the highest number in the MLGE array).
// If the array is shorter than this, expect a segmentation fault.
LocalIntegral* elmInt = locInt.empty() ? 0 : locInt[fe.iel-1];
LocalIntegral* A = integrand.getLocalIntegral(fe.N.size(),fe.iel);
if (!integrand.initElement(MNPC[iel-1],X,nRed,*A)) return false;
if (integrand.getIntegrandType() > 10)
{
@ -579,7 +572,7 @@ bool ASMs1D::integrate (Integrand& integrand,
for (int i = 0; i < nRed; i++)
{
// Local element coordinates of current integration point
fe.xi = xr[i];
fe.xi = xr[i];
// Parameter values of current integration point
fe.u = redpar(i+1,iel);
@ -595,7 +588,7 @@ bool ASMs1D::integrate (Integrand& integrand,
X.t = time.t;
// Compute the reduced integration terms of the integrand
if (!integrand.reducedInt(fe,X))
if (!integrand.reducedInt(*A,fe,X))
return false;
}
}
@ -632,17 +625,19 @@ bool ASMs1D::integrate (Integrand& integrand,
// Evaluate the integrand and accumulate element contributions
fe.detJxW *= 0.5*dL*wg[i];
if (!integrand.evalInt(elmInt,fe,time,X))
if (!integrand.evalInt(*A,fe,time,X))
return false;
}
// Finalize the element quantities
if (!integrand.finalizeElement(elmInt,time))
if (!integrand.finalizeElement(*A,time))
return false;
// Assembly of global system integral
if (!glInt.assemble(elmInt,fe.iel))
if (!glInt.assemble(A->ref(),fe.iel))
return false;
A->destruct();
}
return true;
@ -651,14 +646,13 @@ bool ASMs1D::integrate (Integrand& integrand,
bool ASMs1D::integrate (Integrand& integrand, int lIndex,
GlobalIntegral& glInt,
const TimeDomain& time,
const LintegralVec& locInt)
const TimeDomain& time)
{
if (!curv) return true; // silently ignore empty patches
// Integration of boundary point
FiniteElement fe;
FiniteElement fe(curv->order());
int iel;
switch (lIndex)
{
@ -686,13 +680,8 @@ bool ASMs1D::integrate (Integrand& integrand, int lIndex,
if (!this->getElementCoordinates(Xnod,iel)) return false;
// Initialize element matrices
if (!integrand.initElementBou(MNPC[iel-1])) return false;
// Caution: Unless locInt is empty, we assume it points to an array of
// LocalIntegral pointers, of length at least the number of elements in
// the model (as defined by the highest number in the MLGE array).
// If the array is shorter than this, expect a segmentation fault.
LocalIntegral* elmInt = locInt.empty() ? 0 : locInt[fe.iel-1];
LocalIntegral* A = integrand.getLocalIntegral(fe.N.size(),fe.iel,true);
if (!integrand.initElementBou(MNPC[iel-1],*A)) return false;
// Evaluate basis functions and corresponding derivatives
Matrix dNdu;
@ -713,11 +702,14 @@ bool ASMs1D::integrate (Integrand& integrand, int lIndex,
normal.x = copysign(1.0,Jac(1,1));
// Evaluate the integrand and accumulate element contributions
if (!integrand.evalBou(elmInt,fe,time,X,normal))
if (!integrand.evalBou(*A,fe,time,X,normal))
return false;
// Assembly of global system integral
return glInt.assemble(elmInt,fe.iel);
bool result = glInt.assemble(A->ref(),fe.iel);
A->destruct();
return result;
}

8
src/ASM/ASMs1D.h
View File

@ -120,20 +120,16 @@ public:
//! \param integrand Object with problem-specific data and methods
//! \param glbInt The integrated quantity
//! \param[in] time Parameters for nonlinear/time-dependent simulations
//! \param locInt Vector of element-wise contributions to \a glbInt
virtual bool integrate(Integrand& integrand,
GlobalIntegral& glbInt, const TimeDomain& time,
const LintegralVec& locInt = LintegralVec());
GlobalIntegral& glbInt, const TimeDomain& time);
//! \brief Evaluates a boundary integral over a patch end.
//! \param integrand Object with problem-specific data and methods
//! \param[in] lIndex Local index of the end point
//! \param glbInt The integrated quantity
//! \param[in] time Parameters for nonlinear/time-dependent simulations
//! \param locInt Vector of element-wise contributions to \a glbInt
virtual bool integrate(Integrand& integrand, int lIndex,
GlobalIntegral& glbInt, const TimeDomain& time,
const LintegralVec& locInt = LintegralVec());
GlobalIntegral& glbInt, const TimeDomain& time);
// Post-processing methods

55
src/ASM/ASMs1DLag.C
View File

@ -18,6 +18,7 @@
#include "TimeDomain.h"
#include "FiniteElement.h"
#include "GlobalIntegral.h"
#include "LocalIntegral.h"
#include "IntegrandBase.h"
#include "CoordinateMapping.h"
#include "GaussQuadrature.h"
@ -161,8 +162,7 @@ bool ASMs1DLag::updateCoords (const Vector& displ)
bool ASMs1DLag::integrate (Integrand& integrand,
GlobalIntegral& glInt,
const TimeDomain& time,
const LintegralVec& locInt)
const TimeDomain& time)
{
if (!curv) return true; // silently ignore empty patches
@ -196,20 +196,13 @@ bool ASMs1DLag::integrate (Integrand& integrand,
const int nel = this->getNoElms();
for (int iel = 1; iel <= nel; iel++)
{
fe.iel = MLGE[iel-1];
// Set up nodal point coordinates for current element
if (!this->getElementCoordinates(Xnod,iel)) return false;
// Initialize element quantities
if (!integrand.initElement(MNPC[iel-1],X,nRed)) return false;
// Caution: Unless locInt is empty, we assume it points to an array of
// LocalIntegral pointers, of length at least the number of elements in
// the model (as defined by the highest number in the MLGE array).
// If the array is shorter than this, expect a segmentation fault.
LocalIntegral* elmInt = locInt.empty() ? 0 : locInt[fe.iel-1];
fe.iel = MLGE[iel-1];
LocalIntegral* A = integrand.getLocalIntegral(fe.N.size(),fe.iel);
if (!integrand.initElement(MNPC[iel-1],X,nRed,*A)) return false;
if (integrand.getIntegrandType() > 10)
@ -235,7 +228,7 @@ bool ASMs1DLag::integrate (Integrand& integrand,
X.t = time.t;
// Compute the reduced integration terms of the integrand
if (!integrand.reducedInt(fe,X))
if (!integrand.reducedInt(*A,fe,X))
return false;
}
@ -262,13 +255,19 @@ bool ASMs1DLag::integrate (Integrand& integrand,
X.t = time.t;
// Evaluate the integrand and accumulate element contributions
if (!integrand.evalInt(elmInt,fe,time,X))
if (!integrand.evalInt(*A,fe,time,X))
return false;
}
// Assembly of global system integral
if (!glInt.assemble(elmInt,fe.iel))
// Finalize the element quantities
if (!integrand.finalizeElement(*A,time))
return false;
// Assembly of global system integral
if (!glInt.assemble(A->ref(),fe.iel))
return false;
A->destruct();
}
return true;
@ -277,14 +276,13 @@ bool ASMs1DLag::integrate (Integrand& integrand,
bool ASMs1DLag::integrate (Integrand& integrand, int lIndex,
GlobalIntegral& glInt,
const TimeDomain& time,
const LintegralVec& locInt)
const TimeDomain& time)
{
if (!curv) return true; // silently ignore empty patches
// Integration of boundary point
FiniteElement fe;
FiniteElement fe(curv->order());
int iel;
switch (lIndex)
{
@ -303,20 +301,14 @@ bool ASMs1DLag::integrate (Integrand& integrand, int lIndex,
return false;
}
fe.iel = MLGE[iel-1];
// Set up nodal point coordinates for current element
Matrix Xnod;
if (!this->getElementCoordinates(Xnod,iel)) return false;
// Initialize element quantities
if (!integrand.initElementBou(MNPC[iel-1])) return false;
// Caution: Unless locInt is empty, we assume it points to an array of
// LocalIntegral pointers, of length at least the number of elements in
// the model (as defined by the highest number in the MLGE array).
// If the array is shorter than this, expect a segmentation fault.
LocalIntegral* elmInt = locInt.empty() ? 0 : locInt[fe.iel-1];
fe.iel = MLGE[iel-1];
LocalIntegral* A = integrand.getLocalIntegral(fe.N.size(),fe.iel,true);
if (!integrand.initElementBou(MNPC[iel-1],*A)) return false;
// Evaluate basis functions and corresponding derivatives
Matrix dNdu, Jac;
@ -336,11 +328,14 @@ bool ASMs1DLag::integrate (Integrand& integrand, int lIndex,
normal.x = copysign(1.0,Jac(1,1));
// Evaluate the integrand and accumulate element contributions
if (!integrand.evalBou(elmInt,fe,time,X,normal))
if (!integrand.evalBou(*A,fe,time,X,normal))
return false;
// Assembly of global system integral
return glInt.assemble(elmInt,fe.iel);
bool result = glInt.assemble(A->ref(),fe.iel);
A->destruct();
return result;
}

8
src/ASM/ASMs1DLag.h
View File

@ -75,20 +75,16 @@ public:
//! \param integrand Object with problem-specific data and methods
//! \param glbInt The integrated quantity
//! \param[in] time Parameters for nonlinear/time-dependent simulations
//! \param locInt Vector of element-wise contributions to \a glbInt
virtual bool integrate(Integrand& integrand,
GlobalIntegral& glbInt, const TimeDomain& time,
const LintegralVec& locInt = LintegralVec());
GlobalIntegral& glbInt, const TimeDomain& time);
//! \brief Evaluates a boundary integral over a patch end.
//! \param integrand Object with problem-specific data and methods
//! \param[in] lIndex Local index of the end point
//! \param glbInt The integrated quantity
//! \param[in] time Parameters for nonlinear/time-dependent simulations
//! \param locInt Vector of element-wise contributions to \a glbInt
virtual bool integrate(Integrand& integrand, int lIndex,
GlobalIntegral& glbInt, const TimeDomain& time,
const LintegralVec& locInt = LintegralVec());
GlobalIntegral& glbInt, const TimeDomain& time);
// Post-processing methods

25
src/ASM/ASMs1DSpec.C
View File

@ -17,6 +17,7 @@
#include "TimeDomain.h"
#include "FiniteElement.h"
#include "GlobalIntegral.h"
#include "LocalIntegral.h"
#include "IntegrandBase.h"
#include "CoordinateMapping.h"
#include "Vec3Oper.h"
@ -57,8 +58,7 @@ bool ASMs1DSpec::getGridParameters (RealArray& prm, int nSegPerSpan) const
bool ASMs1DSpec::integrate (Integrand& integrand,
GlobalIntegral& glInt,
const TimeDomain& time,
const LintegralVec& locInt)
const TimeDomain& time)
{
if (!curv) return true; // silently ignore empty patches
@ -98,14 +98,9 @@ bool ASMs1DSpec::integrate (Integrand& integrand,
if (!this->getElementCoordinates(Xnod,iel)) return false;
// Initialize element quantities
if (!integrand.initElement(MNPC[iel-1])) return false;
// Caution: Unless locInt is empty, we assume it points to an array of
// LocalIntegral pointers, of length at least the number of elements in
// the model (as defined by the highest number in the MLGE array).
// If the array is shorter than this, expect a segmentation fault.
LocalIntegral* elmInt = locInt.empty() ? 0 : locInt[MLGE[iel-1]-1];
fe.iel = MLGE[iel-1];
LocalIntegral* A = integrand.getLocalIntegral(p1,fe.iel);
if (!integrand.initElement(MNPC[iel-1],*A)) return false;
// --- Integration loop over integration points ----------------------------
@ -131,13 +126,15 @@ bool ASMs1DSpec::integrate (Integrand& integrand,
// Evaluate the integrand and accumulate element contributions
fe.detJxW *= wg1[i];
if (!integrand.evalInt(elmInt,fe,time,X))
if (!integrand.evalInt(*A,fe,time,X))
return false;
}
// Assembly of global system integral
if (!glInt.assemble(elmInt,MLGE[iel-1]))
if (!glInt.assemble(A->ref(),fe.iel))
return false;
A->destruct();
}
return true;
@ -146,12 +143,10 @@ bool ASMs1DSpec::integrate (Integrand& integrand,
bool ASMs1DSpec::integrate (Integrand& integrand, int lIndex,
GlobalIntegral& glInt,
const TimeDomain& time,
const LintegralVec& locInt)
const TimeDomain& time)
{
if (!curv) return true; // silently ignore empty patches
return false; // not implemented
}

8
src/ASM/ASMs1DSpec.h
View File

@ -43,20 +43,16 @@ public:
//! \param integrand Object with problem-specific data and methods
//! \param glbInt The integrated quantity
//! \param[in] time Parameters for nonlinear/time-dependent simulations
//! \param locInt Vector of element-wise contributions to \a glbInt
virtual bool integrate(Integrand& integrand,
GlobalIntegral& glbInt, const TimeDomain& time,
const LintegralVec& locInt = LintegralVec());
GlobalIntegral& glbInt, const TimeDomain& time);
//! \brief Evaluates a boundary integral over a patch end.
//! \param integrand Object with problem-specific data and methods
//! \param[in] lIndex Local index of the end point
//! \param glbInt The integrated quantity
//! \param[in] time Parameters for nonlinear/time-dependent simulations
//! \param locInt Vector of element-wise contributions to \a glbInt
virtual bool integrate(Integrand& integrand, int lIndex,
GlobalIntegral& glbInt, const TimeDomain& time,
const LintegralVec& locInt = LintegralVec());
GlobalIntegral& glbInt, const TimeDomain& time);
// Post-processing methods

49
src/ASM/ASMs2D.C
View File

@ -30,6 +30,7 @@
#include "Tensor.h"
#include "MPC.h"
#include <algorithm>
#include "ElmNorm.h"
ASMs2D::ASMs2D (unsigned char n_s, unsigned char n_f)
@ -961,8 +962,7 @@ std::ostream& operator<< (std::ostream& os, const Go::BasisDerivsSf& bder)
bool ASMs2D::integrate (Integrand& integrand,
GlobalIntegral& glInt,
const TimeDomain& time,
const LintegralVec& locInt)
const TimeDomain& time)
{
if (!surf) return true; // silently ignore empty patches
@ -1075,15 +1075,12 @@ bool ASMs2D::integrate (Integrand& integrand,
}
// Initialize element quantities
if (!integrand.initElement(MNPC[iel-1],X,nRed*nRed)) return false;
Vectors vec;
LocalIntegral* A = integrand.getLocalIntegral(MNPC[iel-1].size(),iel-1,false);
if (!integrand.initElement(MNPC[iel-1],X,nRed*nRed,vec,*A))
return false;
// Caution: Unless locInt is empty, we assume it points to an array of
// LocalIntegral pointers, of length at least the number of elements in
// the model (as defined by the highest number in the MLGE array).
// If the array is shorter than this, expect a segmentation fault.
LocalIntegral* elmInt = locInt.empty() ? 0 : locInt[fe.iel-1];
if (integrand.getIntegrandType() > 10)
/* if (integrand.getIntegrandType() > 10)
{
// --- Selective reduced integration loop ------------------------------
@ -1114,7 +1111,7 @@ bool ASMs2D::integrate (Integrand& integrand,
return false;
}
}
*/
// --- Integration loop over all Gauss points in each direction ----------
@ -1164,17 +1161,19 @@ bool ASMs2D::integrate (Integrand& integrand,
// Evaluate the integrand and accumulate element contributions
fe.detJxW *= 0.25*dA*wg[i]*wg[j];
if (!integrand.evalInt(elmInt,fe,time,X))
if (!integrand.evalInt(*A,fe,time,X,vec))
return false;
}
// Finalize the element quantities
if (!integrand.finalizeElement(elmInt,time))
if (!integrand.finalizeElement(A,time))
return false;
// Assembly of global system integral
if (!glInt.assemble(elmInt,fe.iel))
return false;
if (!glInt.assemble(A,fe.iel))
return false;
if (!dynamic_cast<ElmNorm*>(A))
delete A;
}
return true;
@ -1183,8 +1182,7 @@ bool ASMs2D::integrate (Integrand& integrand,
bool ASMs2D::integrate (Integrand& integrand, int lIndex,
GlobalIntegral& glInt,
const TimeDomain& time,
const LintegralVec& locInt)
const TimeDomain& time)
{
if (!surf) return true; // silently ignore empty patches
@ -1264,14 +1262,9 @@ bool ASMs2D::integrate (Integrand& integrand, int lIndex,
if (!this->getElementCoordinates(Xnod,iel)) return false;
// Initialize element quantities
if (!integrand.initElementBou(MNPC[iel-1])) return false;
// Caution: Unless locInt is empty, we assume it points to an array of
// LocalIntegral pointers, of length at least the number of elements in
// the model (as defined by the highest number in the MLGE array).
// If the array is shorter than this, expect a segmentation fault.
LocalIntegral* elmInt = locInt.empty() ? 0 : locInt[fe.iel-1];
Vectors vec;
LocalIntegral* A = integrand.getLocalIntegral(MNPC[iel-1].size(),iel-1,true);
if (!integrand.initElementBou(MNPC[iel-1],vec,*A)) return false;
// --- Integration loop over all Gauss points along the edge -------------
@ -1306,13 +1299,15 @@ bool ASMs2D::integrate (Integrand& integrand, int lIndex,
// Evaluate the integrand and accumulate element contributions
fe.detJxW *= 0.5*dS*wg[i];
if (!integrand.evalBou(elmInt,fe,time,X,normal))
if (!integrand.evalBou(*A,fe,time,X,normal,vec))
return false;
}
// Assembly of global system integral
if (!glInt.assemble(elmInt,fe.iel))
if (!glInt.assemble(A,fe.iel))
return false;
if (!dynamic_cast<ElmNorm*>(A))
delete A;
}
return true;

6
src/ASM/ASMs2D.h
View File

@ -209,8 +209,7 @@ public:
//! \param[in] time Parameters for nonlinear/time-dependent simulations
//! \param locInt Vector of element-wise contributions to \a glbInt
virtual bool integrate(Integrand& integrand,
GlobalIntegral& glbInt, const TimeDomain& time,
const LintegralVec& locInt = LintegralVec());
GlobalIntegral& glbInt, const TimeDomain& time);
//! \brief Evaluates a boundary integral over a patch edge.
//! \param integrand Object with problem-specific data and methods
@ -219,8 +218,7 @@ public:
//! \param[in] time Parameters for nonlinear/time-dependent simulations
//! \param locInt Vector of element-wise contributions to \a glbInt
virtual bool integrate(Integrand& integrand, int lIndex,
GlobalIntegral& glbInt, const TimeDomain& time,
const LintegralVec& locInt = LintegralVec());
GlobalIntegral& glbInt, const TimeDomain& time);
// Post-processing methods

47
src/ASM/ASMs2DLag.C
View File

@ -24,6 +24,7 @@
#include "ElementBlock.h"
#include "Utilities.h"
#include "Vec3Oper.h"
#include "ElmNorm.h"
ASMs2DLag::ASMs2DLag (unsigned char n_s, unsigned char n_f)
@ -188,8 +189,7 @@ bool ASMs2DLag::getSize (int& n1, int& n2, int) const
bool ASMs2DLag::integrate (Integrand& integrand,
GlobalIntegral& glInt,
const TimeDomain& time,
const LintegralVec& locInt)
const TimeDomain& time)
{
if (!surf) return true; // silently ignore empty patches
@ -247,18 +247,13 @@ bool ASMs2DLag::integrate (Integrand& integrand,
}
// Initialize element quantities
if (!integrand.initElement(MNPC[iel-1],X,nRed*nRed)) return false;
// Caution: Unless locInt is empty, we assume it points to an array of
// LocalIntegral pointers, of length at least the number of elements in
// the model (as defined by the highest number in the MLGE array).
// If the array is shorter than this, expect a segmentation fault.
LocalIntegral* elmInt = locInt.empty() ? 0 : locInt[fe.iel-1];
Vectors vec;
LocalIntegral* A = integrand.getLocalIntegral(MNPC[iel-1].size(),iel-1,false);
if (!integrand.initElement(MNPC[iel-1],X,nRed*nRed,vec,*A)) return false;
// --- Selective reduced integration loop --------------------------------
if (integrand.getIntegrandType() > 10)
/* if (integrand.getIntegrandType() > 10)
for (int j = 0; j < nRed; j++)
for (int i = 0; i < nRed; i++)
{
@ -286,7 +281,7 @@ bool ASMs2DLag::integrate (Integrand& integrand,
if (!integrand.reducedInt(fe,X))
return false;
}
*/
// --- Integration loop over all Gauss points in each direction ----------
@ -316,17 +311,19 @@ bool ASMs2DLag::integrate (Integrand& integrand,
// Evaluate the integrand and accumulate element contributions
fe.detJxW *= wg[i]*wg[j];
if (!integrand.evalInt(elmInt,fe,time,X))
if (!integrand.evalInt(*A,fe,time,X,vec))
return false;
}
// Finalize the element quantities
if (!integrand.finalizeElement(elmInt,time))
if (!integrand.finalizeElement(A,time))
return false;
// Assembly of global system integral
if (!glInt.assemble(elmInt,fe.iel))
if (!glInt.assemble(A,fe.iel))
return false;
if (!dynamic_cast<ElmNorm*>(A))
delete A;
}
return true;
@ -335,8 +332,7 @@ bool ASMs2DLag::integrate (Integrand& integrand,
bool ASMs2DLag::integrate (Integrand& integrand, int lIndex,
GlobalIntegral& glInt,
const TimeDomain& time,
const LintegralVec& locInt)
const TimeDomain& time)
{
if (!surf) return true; // silently ignore empty patches
@ -402,14 +398,9 @@ bool ASMs2DLag::integrate (Integrand& integrand, int lIndex,
if (!this->getElementCoordinates(Xnod,iel)) return false;
// Initialize element quantities
if (!integrand.initElementBou(MNPC[iel-1])) return false;
// Caution: Unless locInt is empty, we assume it points to an array of
// LocalIntegral pointers, of length at least the number of elements in
// the model (as defined by the highest number in the MLGE array).
// If the array is shorter than this, expect a segmentation fault.
LocalIntegral* elmInt = locInt.empty() ? 0 : locInt[fe.iel-1];
Vectors vec;
LocalIntegral* A = integrand.getLocalIntegral(MNPC[iel-1].size(),iel-1,true);
if (!integrand.initElementBou(MNPC[iel-1],vec,*A)) return false;
// --- Integration loop over all Gauss points along the edge -------------
@ -444,13 +435,15 @@ bool ASMs2DLag::integrate (Integrand& integrand, int lIndex,
// Evaluate the integrand and accumulate element contributions
fe.detJxW *= wg[i];
if (!integrand.evalBou(elmInt,fe,time,X,normal))
if (!integrand.evalBou(*A,fe,time,X,normal,vec))
return false;
}
// Assembly of global system integral
if (!glInt.assemble(elmInt,fe.iel))
if (!glInt.assemble(A,fe.iel))
return false;
if (!dynamic_cast<ElmNorm*>(A))
delete A;
}
return true;

8
src/ASM/ASMs2DLag.h
View File

@ -75,20 +75,16 @@ public:
//! \param integrand Object with problem-specific data and methods
//! \param glbInt The integrated quantity
//! \param[in] time Parameters for nonlinear/time-dependent simulations
//! \param locInt Vector of element-wise contributions to \a glbInt
virtual bool integrate(Integrand& integrand,
GlobalIntegral& glbInt, const TimeDomain& time,
const LintegralVec& locInt = LintegralVec());
GlobalIntegral& glbInt, const TimeDomain& time);
//! \brief Evaluates a boundary integral over a patch edge.
//! \param integrand Object with problem-specific data and methods
//! \param[in] lIndex Local index of the boundary edge
//! \param glbInt The integrated quantity
//! \param[in] time Parameters for nonlinear/time-dependent simulations
//! \param locInt Vector of element-wise contributions to \a glbInt
virtual bool integrate(Integrand& integrand, int lIndex,
GlobalIntegral& glbInt, const TimeDomain& time,
const LintegralVec& locInt = LintegralVec());
GlobalIntegral& glbInt, const TimeDomain& time);
// Post-processing methods

41
src/ASM/ASMs2DSpec.C
View File

@ -21,6 +21,7 @@
#include "CoordinateMapping.h"
#include "Vec3Oper.h"
#include "Legendre.h"
#include "ElmNorm.h"
bool ASMs2DSpec::getGridParameters (RealArray& prm, int dir,
@ -78,8 +79,7 @@ static void evalBasis (int i, int j, int p1, int p2,
bool ASMs2DSpec::integrate (Integrand& integrand,
GlobalIntegral& glInt,
const TimeDomain& time,
const LintegralVec& locInt)
const TimeDomain& time)
{
if (!surf) return true; // silently ignore empty patches
@ -110,14 +110,9 @@ bool ASMs2DSpec::integrate (Integrand& integrand,
if (!this->getElementCoordinates(Xnod,iel)) return false;
// Initialize element quantities
if (!integrand.initElement(MNPC[iel-1])) return false;
// Caution: Unless locInt is empty, we assume it points to an array of
// LocalIntegral pointers, of length at least the number of elements in
// the model (as defined by the highest number in the MLGE array).
// If the array is shorter than this, expect a segmentation fault.
LocalIntegral* elmInt = locInt.empty() ? 0 : locInt[MLGE[iel-1]-1];
Vectors vec;
LocalIntegral* A = integrand.getLocalIntegral(MNPC[iel-1].size(),iel-1,false);
if (!integrand.initElement(MNPC[iel-1],vec,*A)) return false;
// --- Integration loop over all Gauss points in each direction ----------
@ -140,13 +135,15 @@ bool ASMs2DSpec::integrate (Integrand& integrand,
// Evaluate the integrand and accumulate element contributions
fe.detJxW *= wg1(i)*wg2(j);
if (!integrand.evalInt(elmInt,fe,time,X))
if (!integrand.evalInt(*A,fe,time,X,vec))
return false;
}
// Assembly of global system integral
if (!glInt.assemble(elmInt,MLGE[iel-1]))
if (!glInt.assemble(A,MLGE[iel-1]))
return false;
if (!dynamic_cast<ElmNorm*>(A))
delete A;
}
return true;
@ -155,8 +152,7 @@ bool ASMs2DSpec::integrate (Integrand& integrand,
bool ASMs2DSpec::integrate (Integrand& integrand, int lIndex,
GlobalIntegral& glInt,
const TimeDomain& time,
const LintegralVec& locInt)
const TimeDomain& time)
{
if (!surf) return true; // silently ignore empty patches
@ -217,14 +213,9 @@ bool ASMs2DSpec::integrate (Integrand& integrand, int lIndex,
if (!this->getElementCoordinates(Xnod,iel)) return false;
// Initialize element quantities
if (!integrand.initElementBou(MNPC[iel-1])) return false;
// Caution: Unless locInt is empty, we assume it points to an array of
// LocalIntegral pointers, of length at least the number of elements in
// the model (as defined by the highest number in the MLGE array).
// If the array is shorter than this, expect a segmentation fault.
LocalIntegral* elmInt = locInt.empty() ? 0 : locInt[MLGE[iel-1]-1];
Vectors vec;
LocalIntegral* A = integrand.getLocalIntegral(MNPC[iel-1].size(),iel-1,true);
if (!integrand.initElementBou(MNPC[iel-1],vec,*A)) return false;
// --- Integration loop over all Gauss points along the edge -------------
@ -250,13 +241,15 @@ bool ASMs2DSpec::integrate (Integrand& integrand, int lIndex,
// Evaluate the integrand and accumulate element contributions
fe.detJxW *= wg[t2-1][i];
if (!integrand.evalBou(elmInt,fe,time,X,normal))
if (!integrand.evalBou(*A,fe,time,X,normal,vec))
return false;
}
// Assembly of global system integral
if (!glInt.assemble(elmInt,MLGE[iel-1]))
if (!glInt.assemble(A,MLGE[iel-1]))
return false;
if (!dynamic_cast<ElmNorm*>(A))
delete A;
}
return true;

8
src/ASM/ASMs2DSpec.h
View File

@ -43,20 +43,16 @@ public:
//! \param integrand Object with problem-specific data and methods
//! \param glbInt The integrated quantity
//! \param[in] time Parameters for nonlinear/time-dependent simulations
//! \param locInt Vector of element-wise contributions to \a glbInt
virtual bool integrate(Integrand& integrand,
GlobalIntegral& glbInt, const TimeDomain& time,
const LintegralVec& locInt = LintegralVec());
GlobalIntegral& glbInt, const TimeDomain& time);
//! \brief Evaluates a boundary integral over a patch edge.
//! \param integrand Object with problem-specific data and methods
//! \param[in] lIndex Local index of the boundary edge
//! \param glbInt The integrated quantity
//! \param[in] time Parameters for nonlinear/time-dependent simulations
//! \param locInt Vector of element-wise contributions to \a glbInt
virtual bool integrate(Integrand& integrand, int lIndex,
GlobalIntegral& glbInt, const TimeDomain& time,
const LintegralVec& locInt = LintegralVec());
GlobalIntegral& glbInt, const TimeDomain& time);
// Post-processing methods

39
src/ASM/ASMs2Dmx.C
View File

@ -25,6 +25,7 @@
#include "Profiler.h"
#include "Vec3Oper.h"
#include "Vec3.h"
#include "ElmNorm.h"
ASMs2Dmx::ASMs2Dmx (unsigned char n_s,
@ -477,8 +478,7 @@ bool ASMs2Dmx::getSize (int& n1, int& n2, int basis) const
bool ASMs2Dmx::integrate (Integrand& integrand,
GlobalIntegral& glInt,
const TimeDomain& time,
const LintegralVec& locInt)
const TimeDomain& time)
{
if (!surf) return true; // silently ignore empty patches
if (!basis1 || !basis2) return false;
@ -529,17 +529,11 @@ bool ASMs2Dmx::integrate (Integrand& integrand,
// Initialize element quantities
IntVec::const_iterator f2start = MNPC[iel-1].begin() + fe.N1.size();
LocalIntegral* A = integrand.getLocalIntegral(MNPC[iel-1].size(),iel-1,false);
if (!integrand.initElement(IntVec(MNPC[iel-1].begin(),f2start),
IntVec(f2start,MNPC[iel-1].end()),nb1))
IntVec(f2start,MNPC[iel-1].end()),nb1,*A))
return false;
// Caution: Unless locInt is empty, we assume it points to an array of
// LocalIntegral pointers, of length at least the number of elements in
// the model (as defined by the highest number in the MLGE array).
// If the array is shorter than this, expect a segmentation fault.
LocalIntegral* elmInt = locInt.empty() ? 0 : locInt[fe.iel-1];
// --- Integration loop over all Gauss points in each direction ----------
int ip = ((i2-p2)*nGauss*nel1 + i1-p1)*nGauss;
@ -578,13 +572,15 @@ bool ASMs2Dmx::integrate (Integrand& integrand,
// Evaluate the integrand and accumulate element contributions
fe.detJxW *= 0.25*dA*wg[i]*wg[j];
if (!integrand.evalIntMx(elmInt,fe,time,X))
if (!integrand.evalIntMx(*A,fe,time,X))
return false;
}
// Assembly of global system integral
if (!glInt.assemble(elmInt,fe.iel))
if (!glInt.assemble(A,fe.iel))
return false;
if (!dynamic_cast<ElmNorm*>(A))
delete A;
}
return true;
@ -593,8 +589,7 @@ bool ASMs2Dmx::integrate (Integrand& integrand,
bool ASMs2Dmx::integrate (Integrand& integrand, int lIndex,
GlobalIntegral& glInt,
const TimeDomain& time,
const LintegralVec& locInt)
const TimeDomain& time)
{
if (!surf) return true; // silently ignore empty patches
if (!basis1 || !basis2) return false;
@ -678,17 +673,11 @@ bool ASMs2Dmx::integrate (Integrand& integrand, int lIndex,
// Initialize element quantities
IntVec::const_iterator f2start = MNPC[iel-1].begin() + fe.N1.size();
LocalIntegral* A = integrand.getLocalIntegral(MNPC[iel-1].size(),iel-1,true);
if (!integrand.initElementBou(IntVec(MNPC[iel-1].begin(),f2start),
IntVec(f2start,MNPC[iel-1].end()),nb1))
IntVec(f2start,MNPC[iel-1].end()),nb1,*A))
return false;
// Caution: Unless locInt is empty, we assume it points to an array of
// LocalIntegral pointers, of length at least the number of elements in
// the model (as defined by the highest number in the MLGE array).
// If the array is shorter than this, expect a segmentation fault.
LocalIntegral* elmInt = locInt.empty() ? 0 : locInt[fe.iel-1];
// --- Integration loop over all Gauss points along the edge -------------
int ip = (t1 == 1 ? i2-p2 : i1-p1)*nGauss;
@ -732,13 +721,15 @@ bool ASMs2Dmx::integrate (Integrand& integrand, int lIndex,
// Evaluate the integrand and accumulate element contributions
fe.detJxW *= 0.5*dS*wg[i];
if (!integrand.evalBouMx(elmInt,fe,time,X,normal))
if (!integrand.evalBouMx(A,fe,time,X,normal))
return false;
}
// Assembly of global system integral
if (!glInt.assemble(elmInt,fe.iel))
if (!glInt.assemble(A,fe.iel))
return false;
if (!dynamic_cast<ElmNorm*>(A))
delete A;
}
return true;

8
src/ASM/ASMs2Dmx.h
View File

@ -105,20 +105,16 @@ public:
//! \param integrand Object with problem-specific data and methods
//! \param glbInt The integrated quantity
//! \param[in] time Parameters for nonlinear/time-dependent simulations
//! \param locInt Vector of element-wise contributions to \a glbInt
virtual bool integrate(Integrand& integrand,
GlobalIntegral& glbInt, const TimeDomain& time,
const LintegralVec& locInt = LintegralVec());
GlobalIntegral& glbInt, const TimeDomain& time);
//! \brief Evaluates a boundary integral over a patch edge.
//! \param integrand Object with problem-specific data and methods
//! \param[in] lIndex Local index of the boundary edge
//! \param glbInt The integrated quantity
//! \param[in] time Parameters for nonlinear/time-dependent simulations
//! \param locInt Vector of element-wise contributions to \a glbInt
virtual bool integrate(Integrand& integrand, int lIndex,
GlobalIntegral& glbInt, const TimeDomain& time,
const LintegralVec& locInt = LintegralVec());
GlobalIntegral& glbInt, const TimeDomain& time);
// Post-processing methods

39
src/ASM/ASMs2DmxLag.C
View File

@ -23,6 +23,7 @@
#include "GaussQuadrature.h"
#include "Utilities.h"
#include "Vec3Oper.h"
#include "ElmNorm.h"
ASMs2DmxLag::ASMs2DmxLag (unsigned char n_s,
@ -207,8 +208,7 @@ bool ASMs2DmxLag::getSize (int& n1, int& n2, int basis) const
bool ASMs2DmxLag::integrate (Integrand& integrand,
GlobalIntegral& glInt,
const TimeDomain& time,
const LintegralVec& locInt)
const TimeDomain& time)
{
if (!surf) return true; // silently ignore empty patches
@ -251,17 +251,11 @@ bool ASMs2DmxLag::integrate (Integrand& integrand,
// Initialize element quantities
IntVec::const_iterator f2start = MNPC[iel-1].begin() + fe.N1.size();
LocalIntegral* A = integrand.getLocalIntegral(MNPC[iel-1].size(),iel-1,false);
if (!integrand.initElement(IntVec(MNPC[iel-1].begin(),f2start),
IntVec(f2start,MNPC[iel-1].end()),nb1))
IntVec(f2start,MNPC[iel-1].end()),nb1,*A))
return false;
// Caution: Unless locInt is empty, we assume it points to an array of
// LocalIntegral pointers, of length at least the number of elements in
// the model (as defined by the highest number in the MLGE array).
// If the array is shorter than this, expect a segmentation fault.
LocalIntegral* elmInt = locInt.empty() ? 0 : locInt[fe.iel-1];
// --- Integration loop over all Gauss points in each direction ----------
for (int j = 0; j < nGauss; j++)
@ -294,13 +288,15 @@ bool ASMs2DmxLag::integrate (Integrand& integrand,
// Evaluate the integrand and accumulate element contributions
fe.detJxW *= wg[i]*wg[j];
if (!integrand.evalIntMx(elmInt,fe,time,X))
if (!integrand.evalIntMx(*A,fe,time,X))
return false;
}
// Assembly of global system integral
if (!glInt.assemble(elmInt,fe.iel))
if (!glInt.assemble(A,fe.iel))
return false;
if (!dynamic_cast<ElmNorm*>(A))
delete A;
}
return true;
@ -309,8 +305,7 @@ bool ASMs2DmxLag::integrate (Integrand& integrand,
bool ASMs2DmxLag::integrate (Integrand& integrand, int lIndex,
GlobalIntegral& glInt,
const TimeDomain& time,
const LintegralVec& locInt)
const TimeDomain& time)
{
if (!surf) return true; // silently ignore empty patches
@ -367,17 +362,11 @@ bool ASMs2DmxLag::integrate (Integrand& integrand, int lIndex,
// Initialize element quantities
IntVec::const_iterator f2start = MNPC[iel-1].begin() + fe.N1.size();
LocalIntegral* A = integrand.getLocalIntegral(MNPC[iel-1].size(),iel-1,true);
if (!integrand.initElementBou(IntVec(MNPC[iel-1].begin(),f2start),
IntVec(f2start,MNPC[iel-1].end()),nb1))
IntVec(f2start,MNPC[iel-1].end()),nb1,*A))
return false;
// Caution: Unless locInt is empty, we assume it points to an array of
// LocalIntegral pointers, of length at least the number of elements in
// the model (as defined by the highest number in the MLGE array).
// If the array is shorter than this, expect a segmentation fault.
LocalIntegral* elmInt = locInt.empty() ? 0 : locInt[fe.iel-1];
// --- Integration loop over all Gauss points along the edge -------------
for (int i = 0; i < nGauss; i++)
@ -407,13 +396,15 @@ bool ASMs2DmxLag::integrate (Integrand& integrand, int lIndex,
// Evaluate the integrand and accumulate element contributions
fe.detJxW *= wg[i];
if (!integrand.evalBouMx(elmInt,fe,time,X,normal))
if (!integrand.evalBouMx(A,fe,time,X,normal))
return false;
}
// Assembly of global system integral
if (!glInt.assemble(elmInt,fe.iel))
if (!glInt.assemble(A,fe.iel))
return false;
if (!dynamic_cast<ElmNorm*>(A))
delete A;
}
return true;

8
src/ASM/ASMs2DmxLag.h
View File

@ -85,20 +85,16 @@ public:
//! \param integrand Object with problem-specific data and methods
//! \param glbInt The integrated quantity
//! \param[in] time Parameters for nonlinear/time-dependent simulations
//! \param locInt Vector of element-wise contributions to \a glbInt
virtual bool integrate(Integrand& integrand,
GlobalIntegral& glbInt, const TimeDomain& time,
const LintegralVec& locInt = LintegralVec());
GlobalIntegral& glbInt, const TimeDomain& time);
//! \brief Evaluates a boundary integral over a patch edge.
//! \param integrand Object with problem-specific data and methods
//! \param[in] lIndex Local index of the boundary edge
//! \param glbInt The integrated quantity
//! \param[in] time Parameters for nonlinear/time-dependent simulations
//! \param locInt Vector of element-wise contributions to \a glbInt
virtual bool integrate(Integrand& integrand, int lIndex,
GlobalIntegral& glbInt, const TimeDomain& time,
const LintegralVec& locInt = LintegralVec());
GlobalIntegral& glbInt, const TimeDomain& time);
// Post-processing methods

57
src/ASM/ASMs3D.C
View File

@ -26,6 +26,8 @@
#include "Utilities.h"
#include "Profiler.h"
#include "Vec3Oper.h"
#include "ElmMats.h"
#include "ElmNorm.h"
ASMs3D::ASMs3D (unsigned char n_f)
@ -1160,8 +1162,7 @@ void ASMs3D::extractBasis (const Go::BasisDerivs2& spline,
bool ASMs3D::integrate (Integrand& integrand,
GlobalIntegral& glInt,
const TimeDomain& time,
const LintegralVec& locInt)
const TimeDomain& time)
{
if (!svol) return true; // silently ignore empty patches
@ -1283,17 +1284,12 @@ bool ASMs3D::integrate (Integrand& integrand,
}
// Initialize element quantities
if (!integrand.initElement(MNPC[iel-1],X,nRed*nRed*nRed))
Vectors vec;
LocalIntegral* A = integrand.getLocalIntegral(MNPC[iel-1].size(),iel-1,false);
if (!integrand.initElement(MNPC[iel-1],X,nRed*nRed*nRed,vec,*A))
return false;
// Caution: Unless locInt is empty, we assume it points to an array of
// LocalIntegral pointers, of length at least the number of elements in
// the model (as defined by the highest number in the MLGE array).
// If the array is shorter than this, expect a segmentation fault.
LocalIntegral* elmInt = locInt.empty() ? 0 : locInt[fe.iel-1];
if (integrand.getIntegrandType() > 10)
/*if (integrand.getIntegrandType() > 10)
{
// --- Selective reduced integration loop ----------------------------
@ -1326,7 +1322,7 @@ bool ASMs3D::integrate (Integrand& integrand,
if (!integrand.reducedInt(fe,X))
return false;
}
}
}*/
// --- Integration loop over all Gauss points in each direction --------
@ -1381,17 +1377,19 @@ bool ASMs3D::integrate (Integrand& integrand,
// Evaluate the integrand and accumulate element contributions
fe.detJxW *= 0.125*dV*wg[i]*wg[j]*wg[k];
if (!integrand.evalInt(elmInt,fe,time,X))
if (!integrand.evalInt(*A,fe,time,X,vec))
return false;
}
// Finalize the element quantities
if (!integrand.finalizeElement(elmInt,time))
if (!integrand.finalizeElement(A,time))
return false;
// Assembly of global system integral
if (!glInt.assemble(elmInt,fe.iel))
if (!glInt.assemble(A,fe.iel))
return false;
if (!dynamic_cast<ElmNorm*>(A))
delete A;
}
return true;
@ -1400,8 +1398,7 @@ bool ASMs3D::integrate (Integrand& integrand,
bool ASMs3D::integrate (Integrand& integrand, int lIndex,
GlobalIntegral& glInt,
const TimeDomain& time,
const LintegralVec& locInt)
const TimeDomain& time)
{
if (!svol) return true; // silently ignore empty patches
@ -1494,7 +1491,9 @@ bool ASMs3D::integrate (Integrand& integrand, int lIndex,
if (!this->getElementCoordinates(Xnod,iel)) return false;
// Initialize element quantities
if (!integrand.initElementBou(MNPC[iel-1])) return false;
Vectors vec;
LocalIntegral* A = integrand.getLocalIntegral(MNPC[iel-1].size(),iel-1,true);
if (!integrand.initElementBou(MNPC[iel-1],vec,*A)) return false;
// Define some loop control variables depending on which face we are on
int nf1, j1, j2;
@ -1506,13 +1505,6 @@ bool ASMs3D::integrate (Integrand& integrand, int lIndex,
default: nf1 = j1 = j2 = 0;
}
// Caution: Unless locInt is empty, we assume it points to an array of
// LocalIntegral pointers, of length at least the number of elements in
// the model (as defined by the highest number in the MLGE array).
// If the array is shorter than this, expect a segmentation fault.
LocalIntegral* elmInt = locInt.empty() ? 0 : locInt[fe.iel-1];
// --- Integration loop over all Gauss points in each direction --------
int k1, k2, k3;
@ -1559,13 +1551,15 @@ bool ASMs3D::integrate (Integrand& integrand, int lIndex,
// Evaluate the integrand and accumulate element contributions
fe.detJxW *= 0.25*dA*wg[i]*wg[j];
if (!integrand.evalBou(elmInt,fe,time,X,normal))
if (!integrand.evalBou(*A,fe,time,X,normal,vec))
return false;
}
// Assembly of global system integral
if (!glInt.assemble(elmInt,fe.iel))
if (!glInt.assemble(A,fe.iel))
return false;
if (!dynamic_cast<ElmNorm*>(A))
delete A;
}
return true;
@ -1701,12 +1695,13 @@ bool ASMs3D::integrateEdge (Integrand& integrand, int lEdge,
if (!this->getElementCoordinates(Xnod,iel)) return false;
// Initialize element quantities
if (!integrand.initElementBou(MNPC[iel-1])) return false;
Vectors vec;
LocalIntegral* A = integrand.getLocalIntegral(MNPC[iel-1].size(),iel-1,false);
if (!integrand.initElementBou(MNPC[iel-1],vec,*A)) return false;
// --- Integration loop over all Gauss points along the edge -----------
LocalIntegral* elmInt = 0;
for (int i = 0; i < nGauss; i++, ip++)
{
// Parameter values of current integration point
@ -1727,12 +1722,12 @@ bool ASMs3D::integrateEdge (Integrand& integrand, int lEdge,
// Evaluate the integrand and accumulate element contributions
fe.detJxW *= 0.5*dS*wg[i];
if (!integrand.evalBou(elmInt,fe,time,X,tang))
if (!integrand.evalBou(*A,fe,time,X,tang,vec))
return false;
}
// Assembly of global system integral
if (!glInt.assemble(elmInt,fe.iel))
if (!glInt.assemble(A,fe.iel))
return false;
}

8
src/ASM/ASMs3D.h
View File

@ -258,20 +258,16 @@ public:
//! \param integrand Object with problem-specific data and methods
//! \param glbInt The integrated quantity
//! \param[in] time Parameters for nonlinear/time-dependent simulations
//! \param locInt Vector of element-wise contributions to \a glbInt
virtual bool integrate(Integrand& integrand,
GlobalIntegral& glbInt, const TimeDomain& time,
const LintegralVec& locInt = LintegralVec());
GlobalIntegral& glbInt, const TimeDomain& time);
//! \brief Evaluates a boundary integral over a patch face.
//! \param integrand Object with problem-specific data and methods
//! \param[in] lIndex Local index [1,6] of the boundary face
//! \param glbInt The integrated quantity
//! \param[in] time Parameters for nonlinear/time-dependent simulations
//! \param locInt Vector of element-wise contributions to \a glbInt
virtual bool integrate(Integrand& integrand, int lIndex,
GlobalIntegral& glbInt, const TimeDomain& time,
const LintegralVec& locInt = LintegralVec());
GlobalIntegral& glbInt, const TimeDomain& time);
//! \brief Evaluates a boundary integral over a patch edge.
//! \param integrand Object with problem-specific data and methods

56
src/ASM/ASMs3DLag.C
View File

@ -24,6 +24,7 @@
#include "ElementBlock.h"
#include "Utilities.h"
#include "Vec3Oper.h"
#include "ElmNorm.h"
ASMs3DLag::ASMs3DLag (unsigned char n_f) : ASMs3D(n_f), coord(myCoord)
@ -202,8 +203,7 @@ bool ASMs3DLag::getSize (int& n1, int& n2, int& n3, int) const
bool ASMs3DLag::integrate (Integrand& integrand,
GlobalIntegral& glInt,
const TimeDomain& time,
const LintegralVec& locInt)
const TimeDomain& time)
{
if (!svol) return true; // silently ignore empty patches
@ -265,19 +265,14 @@ bool ASMs3DLag::integrate (Integrand& integrand,
}
// Initialize element quantities
if (!integrand.initElement(MNPC[iel-1],X,nRed*nRed*nRed))
Vectors vec;
LocalIntegral* A = integrand.getLocalIntegral(MNPC[iel-1].size(),iel-1,false);
if (!integrand.initElement(MNPC[iel-1],X,nRed*nRed*nRed,vec,*A))
return false;
// Caution: Unless locInt is empty, we assume it points to an array of
// LocalIntegral pointers, of length at least the number of elements in
// the model (as defined by the highest number in the MLGE array).
// If the array is shorter than this, expect a segmentation fault.
LocalIntegral* elmInt = locInt.empty() ? 0 : locInt[fe.iel-1];
// --- Selective reduced integration loop ------------------------------
if (integrand.getIntegrandType() > 10)
/*if (integrand.getIntegrandType() > 10)
for (int k = 0; k < nRed; k++)
for (int j = 0; j < nRed; j++)
for (int i = 0; i < nRed; i++)
@ -310,7 +305,7 @@ bool ASMs3DLag::integrate (Integrand& integrand,
return false;
}
*/
// --- Integration loop over all Gauss points in each direction --------
for (int k = 0; k < nGauss; k++)
@ -342,17 +337,19 @@ bool ASMs3DLag::integrate (Integrand& integrand,
// Evaluate the integrand and accumulate element contributions
fe.detJxW *= wg[i]*wg[j]*wg[k];
if (!integrand.evalInt(elmInt,fe,time,X))
if (!integrand.evalInt(*A,fe,time,X,vec))
return false;
}
// Finalize the element quantities
if (!integrand.finalizeElement(elmInt,time))
if (!integrand.finalizeElement(A,time))
return false;
// Assembly of global system integral
if (!glInt.assemble(elmInt,fe.iel))
if (!glInt.assemble(A,fe.iel))
return false;
if (!dynamic_cast<ElmNorm*>(A))
delete A;
}
return true;
@ -361,8 +358,7 @@ bool ASMs3DLag::integrate (Integrand& integrand,
bool ASMs3DLag::integrate (Integrand& integrand, int lIndex,
GlobalIntegral& glInt,
const TimeDomain& time,
const LintegralVec& locInt)
const TimeDomain& time)
{
if (!svol) return true; // silently ignore empty patches
@ -438,14 +434,9 @@ bool ASMs3DLag::integrate (Integrand& integrand, int lIndex,
if (!this->getElementCoordinates(Xnod,iel)) return false;
// Initialize element quantities
if (!integrand.initElementBou(MNPC[iel-1])) return false;
// Caution: Unless locInt is empty, we assume it points to an array of
// LocalIntegral pointers, of length at least the number of elements in
// the model (as defined by the highest number in the MLGE array).
// If the array is shorter than this, expect a segmentation fault.
LocalIntegral* elmInt = locInt.empty() ? 0 : locInt[fe.iel-1];
Vectors vec;
LocalIntegral* A = integrand.getLocalIntegral(MNPC[iel-1].size(),iel-1,true);
if (!integrand.initElementBou(MNPC[iel-1],vec,*A)) return false;
// --- Integration loop over all Gauss points in each direction --------
@ -493,13 +484,15 @@ bool ASMs3DLag::integrate (Integrand& integrand, int lIndex,
// Evaluate the integrand and accumulate element contributions
fe.detJxW *= wg[i]*wg[j];
if (!integrand.evalBou(elmInt,fe,time,X,normal))
if (!integrand.evalBou(*A,fe,time,X,normal,vec))
return false;
}
// Assembly of global system integral
if (!glInt.assemble(elmInt,fe.iel))
if (!glInt.assemble(A,fe.iel))
return false;
if (!dynamic_cast<ElmNorm*>(A))
delete A;
}
return true;
@ -585,12 +578,13 @@ bool ASMs3DLag::integrateEdge (Integrand& integrand, int lEdge,
if (!this->getElementCoordinates(Xnod,iel)) return false;
// Initialize element quantities
if (!integrand.initElementBou(MNPC[iel-1])) return false;
Vectors vec;
LocalIntegral* A = integrand.getLocalIntegral(MNPC[iel-1].size(),iel-1,false);
if (!integrand.initElementBou(MNPC[iel-1],vec,*A)) return false;
// --- Integration loop over all Gauss points along the edge -----------
LocalIntegral* elmInt = 0;
for (int i = 0; i < nGauss; i++)
{
// Gauss point coordinates on the edge
@ -610,12 +604,12 @@ bool ASMs3DLag::integrateEdge (Integrand& integrand, int lEdge,
X.t = time.t;
// Evaluate the integrand and accumulate element contributions
if (!integrand.evalBou(elmInt,fe,time,X,tangent))
if (!integrand.evalBou(*A,fe,time,X,tangent,vec))
return false;
}
// Assembly of global system integral
if (!glInt.assemble(elmInt,fe.iel))
if (!glInt.assemble(A,fe.iel))
return false;
}

8
src/ASM/ASMs3DLag.h
View File

@ -64,20 +64,16 @@ public:
//! \param integrand Object with problem-specific data and methods
//! \param glbInt The integrated quantity
//! \param[in] time Parameters for nonlinear/time-dependent simulations
//! \param locInt Vector of element-wise contributions to \a glbInt
virtual bool integrate(Integrand& integrand,
GlobalIntegral& glbInt, const TimeDomain& time,
const LintegralVec& locInt = LintegralVec());
GlobalIntegral& glbInt, const TimeDomain& time);
//! \brief Evaluates a boundary integral over a patch face.
//! \param integrand Object with problem-specific data and methods
//! \param[in] lIndex Local index [1,6] of the boundary face
//! \param glbInt The integrated quantity
//! \param[in] time Parameters for nonlinear/time-dependent simulations
//! \param locInt Vector of element-wise contributions to \a glbInt
virtual bool integrate(Integrand& integrand, int lIndex,
GlobalIntegral& glbInt, const TimeDomain& time,
const LintegralVec& locInt = LintegralVec());
GlobalIntegral& glbInt, const TimeDomain& time);
//! \brief Evaluates a boundary integral over a patch edge.
//! \param integrand Object with problem-specific data and methods

51
src/ASM/ASMs3DSpec.C
View File

@ -21,6 +21,7 @@
#include "CoordinateMapping.h"
#include "Vec3Oper.h"
#include "Legendre.h"
#include "ElmNorm.h"
bool ASMs3DSpec::getGridParameters (RealArray& prm, int dir,
@ -79,8 +80,7 @@ static void evalBasis (int i, int j, int k, int p1, int p2, int p3,
bool ASMs3DSpec::integrate (Integrand& integrand,
GlobalIntegral& glInt,
const TimeDomain& time,
const LintegralVec& locInt)
const TimeDomain& time)
{
if (!svol) return true; // silently ignore empty patches
@ -115,14 +115,9 @@ bool ASMs3DSpec::integrate (Integrand& integrand,
if (!this->getElementCoordinates(Xnod,iel)) return false;
// Initialize element quantities
if (!integrand.initElement(MNPC[iel-1])) return false;
// Caution: Unless locInt is empty, we assume it points to an array of
// LocalIntegral pointers, of length at least the number of elements in
// the model (as defined by the highest number in the MLGE array).
// If the array is shorter than this, expect a segmentation fault.
LocalIntegral* elmInt = locInt.empty() ? 0 : locInt[MLGE[iel-1]-1];
Vectors vec;
LocalIntegral* A = integrand.getLocalIntegral(MNPC[iel-1].size(),iel-1,false);
if (!integrand.initElement(MNPC[iel-1],vec,*A)) return false;
// --- Integration loop over all Gauss points in each direction --------
@ -147,13 +142,15 @@ bool ASMs3DSpec::integrate (Integrand& integrand,
// Evaluate the integrand and accumulate element contributions
fe.detJxW *= wg1(i)*wg2(j)*wg3(k);
if (!integrand.evalInt(elmInt,fe,time,X))
if (!integrand.evalInt(*A,fe,time,X,vec))
return false;
}
// Assembly of global system integral
if (!glInt.assemble(elmInt,MLGE[iel-1]))
if (!glInt.assemble(A,MLGE[iel-1]))
return false;
if (!dynamic_cast<ElmNorm*>(A))
delete A;
}
return true;
@ -162,8 +159,7 @@ bool ASMs3DSpec::integrate (Integrand& integrand,
bool ASMs3DSpec::integrate (Integrand& integrand, int lIndex,
GlobalIntegral& glInt,
const TimeDomain& time,
const LintegralVec& locInt)
const TimeDomain& time)
{
if (!svol) return true; // silently ignore empty patches
@ -232,14 +228,9 @@ bool ASMs3DSpec::integrate (Integrand& integrand, int lIndex,
if (!this->getElementCoordinates(Xnod,iel)) return false;
// Initialize element quantities
if (!integrand.initElementBou(MNPC[iel-1])) return false;
// Caution: Unless locInt is empty, we assume it points to an array of
// LocalIntegral pointers, of length at least the number of elements in
// the model (as defined by the highest number in the MLGE array).
// If the array is shorter than this, expect a segmentation fault.
LocalIntegral* elmInt = locInt.empty() ? 0 : locInt[MLGE[iel-1]-1];
Vectors vec;
LocalIntegral* A = integrand.getLocalIntegral(MNPC[iel-1].size(),iel-1,true);
if (!integrand.initElementBou(MNPC[iel-1],vec,*A)) return false;
// --- Integration loop over all Gauss points in each direction --------
@ -267,13 +258,15 @@ bool ASMs3DSpec::integrate (Integrand& integrand, int lIndex,
// Evaluate the integrand and accumulate element contributions
fe.detJxW *= wg[t1-1][i]*wg[t2-1][j];
if (!integrand.evalBou(elmInt,fe,time,X,normal))
if (!integrand.evalBou(*A,fe,time,X,normal,vec))
return false;
}
// Assembly of global system integral
if (!glInt.assemble(elmInt,MLGE[iel-1]))
if (!glInt.assemble(A,MLGE[iel-1]))
return false;
if (!dynamic_cast<ElmNorm*>(A))
delete A;
}
return true;
@ -369,12 +362,12 @@ bool ASMs3DSpec::integrateEdge (Integrand& integrand, int lEdge,
if (!this->getElementCoordinates(Xnod,iel)) return false;
// Initialize element quantities
if (!integrand.initElementBou(MNPC[iel-1])) return false;
Vectors vec;
LocalIntegral* A = integrand.getLocalIntegral(MNPC[iel-1].size(),iel-1,false);
if (!integrand.initElementBou(MNPC[iel-1],vec,*A)) return false;
// --- Integration loop over all Gauss points along the edge -----------
LocalIntegral* elmInt = 0;
for (int i = 0; i < pe; i++)
{
// "Coordinate" on the edge
@ -394,12 +387,12 @@ bool ASMs3DSpec::integrateEdge (Integrand& integrand, int lEdge,
// Evaluate the integrand and accumulate element contributions
fe.detJxW *= wg[lDir][i];
if (!integrand.evalBou(elmInt,fe,time,X,tangent))
if (!integrand.evalBou(*A,fe,time,X,tangent,vec))
return false;
}
// Assembly of global system integral
if (!glInt.assemble(elmInt,MLGE[iel-1]))
if (!glInt.assemble(A,MLGE[iel-1]))
return false;
}

8
src/ASM/ASMs3DSpec.h
View File

@ -42,20 +42,16 @@ public:
//! \param integrand Object with problem-specific data and methods
//! \param glbInt The integrated quantity
//! \param[in] time Parameters for nonlinear/time-dependent simulations
//! \param locInt Vector of element-wise contributions to \a glbInt
virtual bool integrate(Integrand& integrand,
GlobalIntegral& glbInt, const TimeDomain& time,
const LintegralVec& locInt = LintegralVec());
GlobalIntegral& glbInt, const TimeDomain& time);
//! \brief Evaluates a boundary integral over a patch face.
//! \param integrand Object with problem-specific data and methods
//! \param[in] lIndex Local index [1,6] of the boundary face
//! \param glbInt The integrated quantity
//! \param[in] time Parameters for nonlinear/time-dependent simulations
//! \param locInt Vector of element-wise contributions to \a glbInt
virtual bool integrate(Integrand& integrand, int lIndex,
GlobalIntegral& glbInt, const TimeDomain& time,
const LintegralVec& locInt = LintegralVec());
GlobalIntegral& glbInt, const TimeDomain& time);
//! \brief Evaluates a boundary integral over a patch edge.
//! \param integrand Object with problem-specific data and methods

39
src/ASM/ASMs3Dmx.C
View File

@ -25,6 +25,7 @@
#include "Profiler.h"
#include "Vec3Oper.h"
#include "Vec3.h"
#include "ElmNorm.h"
ASMs3Dmx::ASMs3Dmx (unsigned char n_f1, unsigned char n_f2)
@ -503,8 +504,7 @@ bool ASMs3Dmx::getSize (int& n1, int& n2, int& n3, int basis) const
bool ASMs3Dmx::integrate (Integrand& integrand,
GlobalIntegral& glInt,
const TimeDomain& time,
const LintegralVec& locInt)
const TimeDomain& time)
{
if (!svol) return true; // silently ignore empty patches
if (!basis1 || !basis2) return false;
@ -562,17 +562,11 @@ bool ASMs3Dmx::integrate (Integrand& integrand,
// Initialize element quantities
IntVec::const_iterator f2start = MNPC[iel-1].begin() + fe.N1.size();
LocalIntegral* A = integrand.getLocalIntegral(MNPC[iel-1].size(),iel-1,false);
if (!integrand.initElement(IntVec(MNPC[iel-1].begin(),f2start),
IntVec(f2start,MNPC[iel-1].end()),nb1))
IntVec(f2start,MNPC[iel-1].end()),nb1,*A))
return false;
// Caution: Unless locInt is empty, we assume it points to an array of
// LocalIntegral pointers, of length at least the number of elements in
// the model (as defined by the highest number in the MLGE array).
// If the array is shorter than this, expect a segmentation fault.
LocalIntegral* elmInt = locInt.empty() ? 0 : locInt[fe.iel-1];
// --- Integration loop over all Gauss points in each direction --------
int ip = (((i3-p3)*nGauss*nel2 + i2-p2)*nGauss*nel1 + i1-p1)*nGauss;
@ -614,13 +608,15 @@ bool ASMs3Dmx::integrate (Integrand& integrand,
// Evaluate the integrand and accumulate element contributions
fe.detJxW *= 0.125*dV*wg[i]*wg[j]*wg[k];
if (!integrand.evalIntMx(elmInt,fe,time,X))
if (!integrand.evalIntMx(*A,fe,time,X))
return false;
}
// Assembly of global system integral
if (!glInt.assemble(elmInt,fe.iel))
if (!glInt.assemble(A,fe.iel))
return false;
if (!dynamic_cast<ElmNorm*>(A))
delete A;
}
return true;
@ -629,8 +625,7 @@ bool ASMs3Dmx::integrate (Integrand& integrand,
bool ASMs3Dmx::integrate (Integrand& integrand, int lIndex,
GlobalIntegral& glInt,
const TimeDomain& time,
const LintegralVec& locInt)
const TimeDomain& time)
{
if (!svol) return true; // silently ignore empty patches
if (!basis1 || !basis2) return false;
@ -724,8 +719,9 @@ bool ASMs3Dmx::integrate (Integrand& integrand, int lIndex,
// Initialize element quantities
IntVec::const_iterator f2start = MNPC[iel-1].begin() + fe.N1.size();
LocalIntegral* A = integrand.getLocalIntegral(MNPC[iel-1].size(),iel-1,true);
if (!integrand.initElementBou(IntVec(MNPC[iel-1].begin(),f2start),
IntVec(f2start,MNPC[iel-1].end()),nb1))
IntVec(f2start,MNPC[iel-1].end()),nb1,*A))
return false;
// Define some loop control variables depending on which face we are on
@ -738,13 +734,6 @@ bool ASMs3Dmx::integrate (Integrand& integrand, int lIndex,
default: nf1 = j1 = j2 = 0;
}
// Caution: Unless locInt is empty, we assume it points to an array of
// LocalIntegral pointers, of length at least the number of elements in
// the model (as defined by the highest number in the MLGE array).
// If the array is shorter than this, expect a segmentation fault.
LocalIntegral* elmInt = locInt.empty() ? 0 : locInt[fe.iel-1];
// --- Integration loop over all Gauss points in each direction --------
int k1, k2, k3;
@ -802,13 +791,15 @@ bool ASMs3Dmx::integrate (Integrand& integrand, int lIndex,
// Evaluate the integrand and accumulate element contributions
fe.detJxW *= 0.25*dA*wg[i]*wg[j];
if (!integrand.evalBouMx(elmInt,fe,time,X,normal))
if (!integrand.evalBouMx(A,fe,time,X,normal))
return false;
}
// Assembly of global system integral
if (!glInt.assemble(elmInt,fe.iel))
if (!glInt.assemble(A,fe.iel))
return false;
if (!dynamic_cast<ElmNorm*>(A))
delete A;
}
return true;

8
src/ASM/ASMs3Dmx.h
View File

@ -104,20 +104,16 @@ public:
//! \param integrand Object with problem-specific data and methods
//! \param glbInt The integrated quantity
//! \param[in] time Parameters for nonlinear/time-dependent simulations
//! \param locInt Vector of element-wise contributions to \a glbInt
virtual bool integrate(Integrand& integrand,
GlobalIntegral& glbInt, const TimeDomain& time,
const LintegralVec& locInt = LintegralVec());
GlobalIntegral& glbInt, const TimeDomain& time);
//! \brief Evaluates a boundary integral over a patch edge.
//! \param integrand Object with problem-specific data and methods
//! \param[in] lIndex Local index of the boundary edge
//! \param glbInt The integrated quantity
//! \param[in] time Parameters for nonlinear/time-dependent simulations
//! \param locInt Vector of element-wise contributions to \a glbInt
virtual bool integrate(Integrand& integrand, int lIndex,
GlobalIntegral& glbInt, const TimeDomain& time,
const LintegralVec& locInt = LintegralVec());
GlobalIntegral& glbInt, const TimeDomain& time);
// Post-processing methods

39
src/ASM/ASMs3DmxLag.C
View File

@ -23,6 +23,7 @@
#include "GaussQuadrature.h"
#include "Utilities.h"
#include "Vec3Oper.h"
#include "ElmNorm.h"
ASMs3DmxLag::ASMs3DmxLag (unsigned char n_f1, unsigned char n_f2)
@ -227,8 +228,7 @@ bool ASMs3DmxLag::getSize (int& n1, int& n2, int& n3, int basis) const
bool ASMs3DmxLag::integrate (Integrand& integrand,
GlobalIntegral& glInt,
const TimeDomain& time,
const LintegralVec& locInt)
const TimeDomain& time)
{
if (!svol) return true; // silently ignore empty patches
@ -276,17 +276,11 @@ bool ASMs3DmxLag::integrate (Integrand& integrand,
// Initialize element quantities
IntVec::const_iterator f2start = MNPC[iel-1].begin() + fe.N1.size();
LocalIntegral* A = integrand.getLocalIntegral(MNPC[iel-1].size(),iel-1,false);
if (!integrand.initElement(IntVec(MNPC[iel-1].begin(),f2start),
IntVec(f2start,MNPC[iel-1].end()),nb1))
IntVec(f2start,MNPC[iel-1].end()),nb1,*A))
return false;
// Caution: Unless locInt is empty, we assume it points to an array of
// LocalIntegral pointers, of length at least the number of elements in
// the model (as defined by the highest number in the MLGE array).
// If the array is shorter than this, expect a segmentation fault.
LocalIntegral* elmInt = locInt.empty() ? 0 : locInt[fe.iel-1];
// --- Integration loop over all Gauss points in each direction --------
for (int k = 0; k < nGauss; k++)
@ -322,13 +316,15 @@ bool ASMs3DmxLag::integrate (Integrand& integrand,
// Evaluate the integrand and accumulate element contributions
fe.detJxW *= w[i]*w[j]*w[k];
if (!integrand.evalIntMx(elmInt,fe,time,X))
if (!integrand.evalIntMx(*A,fe,time,X))
return false;
}
// Assembly of global system integral
if (!glInt.assemble(elmInt,fe.iel))
if (!glInt.assemble(A,fe.iel))
return false;
if (!dynamic_cast<ElmNorm*>(A))
delete A;
}
return true;
@ -337,8 +333,7 @@ bool ASMs3DmxLag::integrate (Integrand& integrand,
bool ASMs3DmxLag::integrate (Integrand& integrand, int lIndex,
GlobalIntegral& glInt,
const TimeDomain& time,
const LintegralVec& locInt)
const TimeDomain& time)
{
if (!svol) return true; // silently ignore empty patches
@ -402,17 +397,11 @@ bool ASMs3DmxLag::integrate (Integrand& integrand, int lIndex,
// Initialize element quantities
IntVec::const_iterator f2start = MNPC[iel-1].begin() + fe.N1.size();
LocalIntegral* A = integrand.getLocalIntegral(MNPC[iel-1].size(),iel-1,true);
if (!integrand.initElementBou(IntVec(MNPC[iel-1].begin(),f2start),
IntVec(f2start,MNPC[iel-1].end()),nb1))
IntVec(f2start,MNPC[iel-1].end()),nb1,*A))
return false;
// Caution: Unless locInt is empty, we assume it points to an array of
// LocalIntegral pointers, of length at least the number of elements in
// the model (as defined by the highest number in the MLGE array).
// If the array is shorter than this, expect a segmentation fault.
LocalIntegral* elmInt = locInt.empty() ? 0 : locInt[fe.iel-1];
// --- Integration loop over all Gauss points in each direction --------
for (int j = 0; j < nGauss; j++)
@ -444,13 +433,15 @@ bool ASMs3DmxLag::integrate (Integrand& integrand, int lIndex,
// Evaluate the integrand and accumulate element contributions
fe.detJxW *= wg[i]*wg[j];
if (!integrand.evalBouMx(elmInt,fe,time,X,normal))
if (!integrand.evalBouMx(A,fe,time,X,normal))
return false;
}
// Assembly of global system integral
if (!glInt.assemble(elmInt,fe.iel))
if (!glInt.assemble(A,fe.iel))
return false;
if (!dynamic_cast<ElmNorm*>(A))
delete A;
}
return true;

8
src/ASM/ASMs3DmxLag.h
View File

@ -84,20 +84,16 @@ public:
//! \param integrand Object with problem-specific data and methods
//! \param glbInt The integrated quantity
//! \param[in] time Parameters for nonlinear/time-dependent simulations
//! \param locInt Vector of element-wise contributions to \a glbInt
virtual bool integrate(Integrand& integrand,
GlobalIntegral& glbInt, const TimeDomain& time,
const LintegralVec& locInt = LintegralVec());
GlobalIntegral& glbInt, const TimeDomain& time);
//! \brief Evaluates a boundary integral over a patch face.
//! \param integrand Object with problem-specific data and methods
//! \param[in] lIndex Local index [1,6] of the boundary face
//! \param glbInt The integrated quantity
//! \param[in] time Parameters for nonlinear/time-dependent simulations
//! \param locInt Vector of element-wise contributions to \a glbInt
virtual bool integrate(Integrand& integrand, int lIndex,
GlobalIntegral& glbInt, const TimeDomain& time,
const LintegralVec& locInt = LintegralVec());
GlobalIntegral& glbInt, const TimeDomain& time);
// Post-processing methods

3
src/ASM/ElmMats.h
View File

@ -70,6 +70,9 @@ public:
return b.front();
}
//! \brief Virtual destruction method to clean up after numerical integration.
virtual void destruct() { delete this; }
std::vector<Matrix> A; //!< The element coefficient matrices
std::vector<Vector> b; //!< The element right-hand-side vectors

6
src/ASM/ElmNorm.h
View File

@ -41,6 +41,12 @@ public:
//! \brief Returns the number of norm values.
size_t size() const { return nnv; }
//! \brief Virtual destruction method to clean up after numerical integration.
//! \details For this class we only clear the two solution vector containers
//! (to save memory). We do NOT delete the object itself, since it might be
//! needed in a second integration loop over element boundaries, for instance.
virtual void destruct() { vec.clear(); psol.clear(); }
private:
double* ptr; //!< Pointer to the actual norm values
size_t nnv; //!< Number of norm values

65
src/ASM/Integrand.h
View File

@ -60,11 +60,35 @@ public:
// Element-level initialization interface
// ======================================
//! \brief Get a local integral contribution container for a given element
//! \param[in] nen Number of DOFs on element
//! \param[in] iEl The element number
//! \param[in] neumann Whether or not we are assembling Neumann b.c's
virtual LocalIntegral* getLocalIntegral(size_t nen, size_t iEl,
bool neumann = false) const = 0;
//! \brief Returns a local integral contribution object for the given element.
//! \param[in] nen1 Number of nodes on element for basis 1
//! \param[in] nen2 Number of nodes on element for basis 2
//! \param[in] iEl Global element number (1-based)
//! \param[in] neumann Whether or not we are assembling Neumann BCs
//!
//! \details This form is used for mixed formulations only.
//! The default implementation just forwards to the single-basis version.
//! Reimplement this method if your mixed formulation requires specialized
//! local integral objects.
virtual LocalIntegral* getLocalIntegral(size_t nen1, size_t nen2,
size_t iEl,
bool neumann = false) const
{
return this->getLocalIntegral(nen1,iEl,neumann);
}
//! \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
//! \param elmInt Local integral for element
//!
//! \details This method is invoked once before starting the numerical
//! integration loop over the Gaussian quadrature points over an element.
@ -73,9 +97,11 @@ public:
//! 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) = 0;
const Vec3& X0, size_t nPt,
LocalIntegral& elmInt) = 0;
//! \brief Initializes current element for numerical integration.
//! \param[in] MNPC Matrix of nodal point correspondance for current element
//! \param elmInt Local integral for element
//!
//! \details This method is invoked once before starting the numerical
//! integration loop over the Gaussian quadrature points over an element.
@ -84,23 +110,30 @@ public:
//! 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;
virtual bool initElement(const std::vector<int>& MNPC,
LocalIntegral& elmInt) = 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
//! \param elmInt Local integral for element
virtual bool initElement(const std::vector<int>& MNPC1,
const std::vector<int>& MNPC2, size_t n1) = 0;
const std::vector<int>& MNPC2, size_t n1,
LocalIntegral& elmInt) = 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;
//! \param elmInt Local integral for element
virtual bool initElementBou(const std::vector<int>& MNPC,
LocalIntegral& elmInt) = 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
//! \param elmInt Local integral for element
virtual bool initElementBou(const std::vector<int>& MNPC1,
const std::vector<int>& MNPC2, size_t n1) = 0;
const std::vector<int>& MNPC2, size_t n1,
LocalIntegral& elmInt) = 0;
// Integrand evaluation interface
@ -117,13 +150,15 @@ public:
virtual int getIntegrandType() const { return 1; }
//! \brief Evaluates reduced integration terms 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] X Cartesian coordinates of current integration point
//!
//! \details Reimplement this method to evaluate terms at other points than
//! the regular integration points. This method is invoked in a separate loop
//! prior to the main integration point loop.
virtual bool reducedInt(const FiniteElement& fe, const Vec3& X) const
virtual bool reducedInt(LocalIntegral& elmInt,
const FiniteElement& fe, const Vec3& X) const
{
return false;
}
@ -136,7 +171,7 @@ public:
//!
//! \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,
virtual bool evalInt(LocalIntegral& elmInt, const FiniteElement& fe,
const TimeDomain& time, const Vec3& X) const
{
return this->evalInt(elmInt,fe,X);
@ -151,7 +186,7 @@ public:
//! \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,
virtual bool evalIntMx(LocalIntegral& elmInt, const MxFiniteElement& fe,
const TimeDomain& time, const Vec3& X) const
{
return this->evalIntMx(elmInt,fe,X);
@ -164,7 +199,7 @@ public:
//! 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&)
virtual bool finalizeElement(LocalIntegral&, const TimeDomain&)
{
return true;
}
@ -178,7 +213,7 @@ public:
//!
//! \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,
virtual bool evalBou(LocalIntegral& elmInt, const FiniteElement& fe,
const TimeDomain& time,
const Vec3& X, const Vec3& normal) const
{
@ -195,7 +230,7 @@ public:
//! \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,
virtual bool evalBouMx(LocalIntegral& elmInt, const MxFiniteElement& fe,
const TimeDomain& time,
const Vec3& X, const Vec3& normal) const
{
@ -204,17 +239,17 @@ public:
protected:
//! \brief Evaluates the integrand at interior points for stationary problems.
virtual bool evalInt(LocalIntegral*&, const FiniteElement& fe,
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,
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&,
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&,
virtual bool evalBouMx(LocalIntegral&, const MxFiniteElement&,
const Vec3&, const Vec3&) const { return false; }
public:

110
src/ASM/IntegrandBase.C
View File

@ -17,12 +17,6 @@
#include "Utilities.h"
IntegrandBase::~IntegrandBase ()
{
if (myMats) delete myMats;
}
void IntegrandBase::resetSolution ()
{
for (size_t i = 0; i < primsol.size(); i++)
@ -31,22 +25,21 @@ void IntegrandBase::resetSolution ()
bool IntegrandBase::initElement (const std::vector<int>& MNPC,
const Vec3&, size_t)
const Vec3&, size_t, LocalIntegral& elmInt)
{
return this->initElement(MNPC);
return this->initElement(MNPC,elmInt);
}
bool IntegrandBase::initElement (const std::vector<int>& MNPC)
bool IntegrandBase::initElement (const std::vector<int>& MNPC,
LocalIntegral& elmInt)
{
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++)
elmInt.vec.resize(primsol.size());
for (size_t i = 0; i < primsol.size() && ierr == 0; i++)
if (!primsol[i].empty())
ierr = utl::gather(MNPC,npv,primsol[i],mySols[i]);
ierr = utl::gather(MNPC,npv,primsol[i],elmInt.vec[i]);
if (ierr == 0) return true;
@ -57,7 +50,8 @@ bool IntegrandBase::initElement (const std::vector<int>& MNPC)
bool IntegrandBase::initElement (const std::vector<int>&,
const std::vector<int>&, size_t)
const std::vector<int>&, size_t,
LocalIntegral&)
{
std::cerr <<" *** IntegrandBase::initElement must be reimplemented"
<<" for mixed problems."<< std::endl;
@ -65,15 +59,14 @@ bool IntegrandBase::initElement (const std::vector<int>&,
}
bool IntegrandBase::initElementBou (const std::vector<int>& MNPC)
bool IntegrandBase::initElementBou (const std::vector<int>& MNPC,
LocalIntegral& elmInt)
{
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());
elmInt.vec.resize(1);
if (!primsol.empty() && !primsol.front().empty())
ierr = utl::gather(MNPC,npv,primsol.front(),elmInt.vec.front());
if (ierr == 0) return true;
@ -84,7 +77,8 @@ bool IntegrandBase::initElementBou (const std::vector<int>& MNPC)
bool IntegrandBase::initElementBou (const std::vector<int>&,
const std::vector<int>&, size_t)
const std::vector<int>&, size_t,
LocalIntegral&)
{
std::cerr <<" *** IntegrandBase::initElementBou must be reimplemented"
<<" for mixed problems."<< std::endl;
@ -168,75 +162,81 @@ Vector& NormBase::getProjection (size_t i)
return dummy;
}
else if (prjsol.size() < i)
{
prjsol.resize(i);
mySols.resize(i);
}
return prjsol[i-1];
}
bool NormBase::initProjection (const std::vector<int>& MNPC)
bool NormBase::initProjection (const std::vector<int>& MNPC,
LocalIntegral& elmInt)
{
// Extract projected solution vectors for this element
Vectors& psol = static_cast<ElmNorm&>(elmInt).psol;
psol.resize(prjsol.size());
int ierr = 0;
for (size_t i = 0; i < mySols.size() && ierr == 0; i++)
for (size_t i = 0; i < prjsol.size() && ierr == 0; i++)
if (!prjsol[i].empty())
ierr = utl::gather(MNPC,nrcmp,prjsol[i],mySols[i]);
ierr = utl::gather(MNPC,nrcmp,prjsol[i],psol[i]);
if (ierr == 0) return true;
std::cerr <<" *** NormBase::initProjection: Detected "
<< ierr <<" node numbers out of range."<< std::endl;
std::cerr <<" *** NormBase::initProjection: Detected "<< ierr
<<" node numbers out of range."<< std::endl;
return false;
}
bool NormBase::initElement (const std::vector<int>& MNPC,
const Vec3& Xc, size_t nPt)
LocalIntegral* NormBase::getLocalIntegral(size_t nen, size_t iEl,
bool neumann) const
{
return this->initProjection(MNPC) && myProblem.initElement(MNPC,Xc,nPt);
if (lints && iEl > 0 && iEl <= lints->size()) return (*lints)[iEl-1];
// Element norms are not requested, so allocate one internally instead that
// will delete itself when invoking the destruct method.
return new ElmNorm(this->getNoFields());
}
bool NormBase::initElement (const std::vector<int>& MNPC)
bool NormBase::initElement (const std::vector<int>& MNPC,
const Vec3& Xc, size_t nPt,
LocalIntegral& elmInt)
{
return this->initProjection(MNPC) && myProblem.initElement(MNPC);
return this->initProjection(MNPC,elmInt) &&
myProblem.initElement(MNPC,Xc,nPt,elmInt);
}
bool NormBase::initElement (const std::vector<int>& MNPC,
LocalIntegral& elmInt)
{
return this->initProjection(MNPC,elmInt) &&
myProblem.initElement(MNPC,elmInt);
}
bool NormBase::initElement (const std::vector<int>& MNPC1,
const std::vector<int>& MNPC2, size_t n1)
const std::vector<int>& MNPC2, size_t n1,
LocalIntegral& elmInt)
{
return this->initProjection(MNPC1) && myProblem.initElement(MNPC1,MNPC2,n1);
return this->initProjection(MNPC1,elmInt) &&
myProblem.initElement(MNPC1,MNPC2,n1,elmInt);
}
bool NormBase::initElementBou (const std::vector<int>& MNPC)
bool NormBase::initElementBou (const std::vector<int>& MNPC,
LocalIntegral& elmInt)
{
return myProblem.initElementBou(MNPC);
return myProblem.initElementBou(MNPC,elmInt);
}
bool NormBase::initElementBou (const std::vector<int>& MNPC1,
const std::vector<int>& MNPC2, size_t n1)
const std::vector<int>& MNPC2, size_t n1,
LocalIntegral& elmInt)
{
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,nn);
memset(data,0,buf.size()*sizeof(double));
elmInt = &buf;
return buf;
return myProblem.initElementBou(MNPC1,MNPC2,n1,elmInt);
}

80
src/ASM/IntegrandBase.h
View File

@ -19,12 +19,9 @@
#include "Function.h"
class NormBase;
class ElmMats;
class ElmNorm;
class AnaSol;
class VTF;
/*!
\brief Base class representing a system level integrated quantity.
*/
@ -33,11 +30,11 @@ class IntegrandBase : public Integrand
{
protected:
//! \brief The default constructor is protected to allow sub-classes only.
IntegrandBase() : m_mode(SIM::STATIC), npv(0){}
IntegrandBase() : npv(1), m_mode(SIM::INIT) {}
public:
//! \brief The destructor frees the dynamically allocated data objects.
virtual ~IntegrandBase();
//! \brief Empty destructor.
virtual ~IntegrandBase() {}
//! \brief Prints out the problem definition to the given output stream.
virtual void print(std::ostream&) const {}
@ -66,6 +63,7 @@ public:
//! \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
//! \param[in] elmInt Local integral for element
//!
//! \details This method is invoked once before starting the numerical
//! integration loop over the Gaussian quadrature points over an element.
@ -74,31 +72,38 @@ public:
//! The default implementation forwards to an overloaded method not taking
//! \a X0 and \a nPt as arguments.
virtual bool initElement(const std::vector<int>& MNPC,
const Vec3& X0, size_t nPt);
const Vec3& X0, size_t nPt, LocalIntegral& elmInt);
//! \brief Initializes current element for numerical integration.
//! \param[in] MNPC Matrix of nodal point correspondance for current element
//! \param[in] elmInt Local integral for 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.
virtual bool initElement(const std::vector<int>& MNPC);
virtual bool initElement(const std::vector<int>& MNPC, LocalIntegral& elmInt);
//! \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
//! \param[in] elmInt Local integral for element
virtual bool initElement(const std::vector<int>& MNPC1,
const std::vector<int>& MNPC2, size_t n1);
const std::vector<int>& MNPC2, size_t n1,
LocalIntegral& elmInt);
//! \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);
//! \param[in] elmInt Local integral for element
virtual bool initElementBou(const std::vector<int>& MNPC,
LocalIntegral& elmInt);
//! \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
//! \param[in] elmInt Local integral for element
virtual bool initElementBou(const std::vector<int>& MNPC1,
const std::vector<int>& MNPC2, size_t n1);
const std::vector<int>& MNPC2, size_t n1,
LocalIntegral& elmInt);
// Solution field evaluation interface
@ -172,12 +177,9 @@ public:
void resetSolution();
protected:
Vectors primsol; //!< Primary solution vectors for current patch
ElmMats* myMats; //!< Local element matrices
Vectors mySols; //!< Local element vectors of the primary solution
SIM::SolutionMode m_mode; //!< Current solution mode
unsigned short int npv; //!< Number of primary solution variables per node
Vectors primsol; //!< Primary solution vectors for current patch
unsigned short int npv; //!< Number of primary solution variables per node
SIM::SolutionMode m_mode; //!< Current solution mode
};
@ -189,7 +191,7 @@ class NormBase : public Integrand
{
protected:
//! \brief The default constructor is protected to allow sub-classes only.
NormBase(IntegrandBase& p) : myProblem(p), nrcmp(0) {}
NormBase(IntegrandBase& p) : myProblem(p), nrcmp(0), lints(0) {}
public:
//! \brief Empty destructor.
@ -198,23 +200,39 @@ public:
//! \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,
const Vec3& X0, size_t nPt);
//! \brief Set a vector of LocalIntegrals to be used during norm integration
void setLocalIntegrals(LintegralVec* elementNorms)
{
lints = elementNorms;
}
//! \brief Get a local integral contribution container for a given element
//! \param[in] nen Number of DOFs on element
//! \param[in] iEl The element number
//! \param[in] neumann Whether or not we are assembling Neumann b.c's
virtual LocalIntegral* getLocalIntegral(size_t nen, size_t iEl,
bool neumann) const;
//! \brief Initializes current element for numerical integration.
virtual bool initElement(const std::vector<int>& MNPC);
virtual bool initElement(const std::vector<int>& MNPC,
const Vec3& X0, size_t nPt, LocalIntegral& elmInt);
//! \brief Initializes current element for numerical integration.
virtual bool initElement(const std::vector<int>& MNPC, LocalIntegral& elmInt);
//! \brief Initializes current element for numerical integration (mixed).
virtual bool initElement(const std::vector<int>& MNPC1,
const std::vector<int>& MNPC2, size_t n1);
const std::vector<int>& MNPC2, size_t n1,
LocalIntegral& elmInt);
//! \brief Initializes current element for boundary integration.
virtual bool initElementBou(const std::vector<int>& MNPC);
virtual bool initElementBou(const std::vector<int>& MNPC,
LocalIntegral& elmInt);
//! \brief Initializes current element for boundary integration (mixed).
virtual bool initElementBou(const std::vector<int>& MNPC1,
const std::vector<int>& MNPC2, size_t n1);
const std::vector<int>& MNPC2, size_t n1,
LocalIntegral& elmInt);
//! \brief Returns whether this norm has explicit boundary contributions.
virtual bool hasBoundaryTerms() const { return false; }
@ -231,24 +249,14 @@ public:
//! \brief Accesses a projected secondary solution vector of current patch.
Vector& getProjection(size_t i);
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);
private:
//! \brief Initializes projected field for current element.
bool initProjection(const std::vector<int>& MNPC);
bool initProjection(const std::vector<int>& MNPC, LocalIntegral& elmInt);
protected:
IntegrandBase& myProblem; //!< The problem-specific data
Vectors prjsol; //!< Projected secondary solution vectors for current patch
Vectors mySols; //!< Local element vectors of projected secondary solutions
unsigned short int nrcmp; //!< Number of projected solution components
};

81
src/ASM/LR/ASMu2D.C
View File

@ -23,6 +23,7 @@
#include "TimeDomain.h"
#include "FiniteElement.h"
#include "GlobalIntegral.h"
#include "LocalIntegral.h"
#include "IntegrandBase.h"
#include "CoordinateMapping.h"
#include "GaussQuadrature.h"
@ -800,8 +801,7 @@ void ASMu2D::extractBasis (const Go::BasisDerivsSf2& spline,
bool ASMu2D::integrate (Integrand& integrand,
GlobalIntegral& glInt,
const TimeDomain& time,
const LintegralVec& locInt)
const TimeDomain& time)
{
if (!lrspline) return true; // silently ignore empty patches
@ -827,12 +827,11 @@ bool ASMu2D::integrate (Integrand& integrand,
// === Assembly loop over all elements in the patch ==========================
std::vector<LR::Element*>::iterator el;
int iel = 1;
for(el = lrspline->elementBegin(); el!=lrspline->elementEnd(); el++, iel++)
std::vector<LR::Element*>::iterator el = lrspline->elementBegin();
for (int iel = 1; el != lrspline->elementEnd(); el++, iel++)
{
FiniteElement fe(myMNPC[iel-1].size());
fe.iel = myMLGE[iel-1];
fe.iel = MLGE[iel-1];
// Get element area in the parameter space
double dA = this->getParametricArea(iel);
@ -896,16 +895,10 @@ bool ASMu2D::integrate (Integrand& integrand,
#endif
// Initialize element quantities
if (!integrand.initElement(myMNPC[iel-1],X,nRed*nRed)) return false;
LocalIntegral* A = integrand.getLocalIntegral(fe.N.size(),fe.iel);
if (!integrand.initElement(MNPC[iel-1],X,nRed*nRed,*A)) return false;
// Caution: Unless locInt is empty, we assume it points to an array of
// LocalIntegral pointers, of length at least the number of elements in
// the model (as defined by the highest number in the MLGE array).
// If the array is shorter than this, expect a segmentation fault.
LocalIntegral* elmInt = locInt.empty() ? 0 : locInt[fe.iel-1];
if (integrand.getIntegrandType() > 10)
/* if (integrand.getIntegrandType() > 10)
{
// --- Selective reduced integration loop ------------------------------
@ -937,7 +930,7 @@ bool ASMu2D::integrate (Integrand& integrand,
return false;
}
}
*/
// --- Integration loop over all Gauss points in each direction ----------
@ -991,17 +984,19 @@ bool ASMu2D::integrate (Integrand& integrand,
// Evaluate the integrand and accumulate element contributions
fe.detJxW *= 0.25*dA*wg[i]*wg[j];
if (!integrand.evalInt(elmInt,fe,time,X))
if (!integrand.evalInt(*A,fe,time,X))
return false;
}
// Finalize the element quantities
if (!integrand.finalizeElement(elmInt,time))
if (!integrand.finalizeElement(*A,time))
return false;
// Assembly of global system integral
if (!glInt.assemble(elmInt,fe.iel))
if (!glInt.assemble(A->ref(),fe.iel))
return false;
A->destruct();
}
return true;
@ -1010,8 +1005,7 @@ bool ASMu2D::integrate (Integrand& integrand,
bool ASMu2D::integrate (Integrand& integrand, int lIndex,
GlobalIntegral& glInt,
const TimeDomain& time,
const LintegralVec& locInt)
const TimeDomain& time)
{
if (!lrspline) return true; // silently ignore empty patches
@ -1048,14 +1042,9 @@ bool ASMu2D::integrate (Integrand& integrand, int lIndex,
// === Assembly loop over all elements on the patch edge =====================
std::vector<LR::Element*>::iterator el;
int iel = 1;
for(el = lrspline->elementBegin(); el!=lrspline->elementEnd(); el++, iel++)
std::vector<LR::Element*>::iterator el = lrspline->elementBegin();
for (int iel = 1; el != lrspline->elementEnd(); el++, iel++)
{
FiniteElement fe((**el).nBasisFunctions());
fe.xi = fe.eta = edgeDir < 0 ? -1.0 : 1.0;
fe.iel = myMLGE[iel-1];
// Skip elements that are not on current boundary edge
bool skipMe = false;
switch (edgeDir)
@ -1075,13 +1064,11 @@ bool ASMu2D::integrate (Integrand& integrand, int lIndex,
if (!this->getElementCoordinates(Xnod,iel)) return false;
// Initialize element quantities
if (!integrand.initElementBou(myMNPC[iel-1])) return false;
// Caution: Unless locInt is empty, we assume it points to an array of
// LocalIntegral pointers, of length at least the number of elements in
// the model (as defined by the highest number in the MLGE array).
// If the array is shorter than this, expect a segmentation fault.
LocalIntegral* elmInt = locInt.empty() ? 0 : locInt[fe.iel-1];
FiniteElement fe((**el).nBasisFunctions());
fe.iel = MLGE[iel-1];
fe.xi = fe.eta = edgeDir < 0 ? -1.0 : 1.0;
LocalIntegral* A = integrand.getLocalIntegral(fe.N.size(),fe.iel,true);
if (!integrand.initElementBou(MNPC[iel-1],*A)) return false;
// get integration gauss points over this element
int dim = 2-abs(edgeDir);
@ -1099,33 +1086,35 @@ bool ASMu2D::integrate (Integrand& integrand, int lIndex,
fe.eta = xg[i];
fe.u = gpar[0][i];
fe.v = gpar[1][i];
// Evaluate basis function derivatives at current integration points
Go::BasisDerivsSf spline;
lrspline->computeBasis(fe.u, fe.v, spline, iel-1);
// Fetch basis function derivatives at current integration point
extractBasis(spline,fe.N,dNdu);
// Compute basis function derivatives and the edge normal
fe.detJxW = utl::Jacobian(Jac,normal,fe.dNdX,Xnod,dNdu,t1,t2);
if (fe.detJxW == 0.0) continue; // skip singular points
if (edgeDir < 0) normal *= -1.0;
// Cartesian coordinates of current integration point
X = Xnod * fe.N;
X.t = time.t;
// Evaluate the integrand and accumulate element contributions
fe.detJxW *= 0.5*dS*wg[i];
if (!integrand.evalBou(elmInt,fe,time,X,normal))
if (!integrand.evalBou(*A,fe,time,X,normal))
return false;
}
// Assembly of global system integral
if (!glInt.assemble(elmInt,fe.iel))
return false;
}
// Assembly of global system integral
if (!glInt.assemble(A,fe.iel))
return false;
A->destruct();
}
return true;

8
src/ASM/LR/ASMu2D.h
View File

@ -176,20 +176,16 @@ public:
//! \param integrand Object with problem-specific data and methods
//! \param glbInt The integrated quantity
//! \param[in] time Parameters for nonlinear/time-dependent simulations
//! \param locInt Vector of element-wise contributions to \a glbInt
virtual bool integrate(Integrand& integrand,
GlobalIntegral& glbInt, const TimeDomain& time,
const LintegralVec& locInt = LintegralVec());
GlobalIntegral& glbInt, const TimeDomain& time);
//! \brief Evaluates a boundary integral over a patch edge.
//! \param integrand Object with problem-specific data and methods
//! \param[in] lIndex Local index of the boundary edge
//! \param glbInt The integrated quantity
//! \param[in] time Parameters for nonlinear/time-dependent simulations
//! \param locInt Vector of element-wise contributions to \a glbInt
virtual bool integrate(Integrand& integrand, int lIndex,
GlobalIntegral& glbInt, const TimeDomain& time,
const LintegralVec& locInt = LintegralVec());
GlobalIntegral& glbInt, const TimeDomain& time);
// Post-processing methods

9
src/ASM/LocalIntegral.h
View File

@ -14,6 +14,8 @@
#ifndef _LOCAL_INTEGRAL_H
#define _LOCAL_INTEGRAL_H
#include "MatVec.h"
/*!
\brief Abstract base class representing an element level integrated quantity.
@ -24,9 +26,16 @@ class LocalIntegral
protected:
//! \brief The default constructor is protected to allow sub-classes only.
LocalIntegral() {}
public:
//! \brief Empty destructor.
virtual ~LocalIntegral() {}
//! \brief Virtual destruction method to clean up after numerical integration.
virtual void destruct() { delete this; }
//! \brief Returns the LocalIntegral object to assemble into the global one.
virtual const LocalIntegral* ref() const { return this; }
Vectors vec; //!< Element-level primary solution vectors
};
#endif

36
src/Integrands/Elasticity.C
View File

@ -106,10 +106,10 @@ void Elasticity::setMode (SIM::SolutionMode mode)
case SIM::BUCKLING:
myMats->resize(2,0);
mySols.resize(1);
// mySols.resize(1);
eKm = &myMats->A[0];
eKg = &myMats->A[1];
eV = &mySols[0];
// eV = &mySols[0];
break;
case SIM::STIFF_ONLY:
@ -132,13 +132,13 @@ void Elasticity::setMode (SIM::SolutionMode mode)
case SIM::RECOVERY:
myMats->rhsOnly = true;
mySols.resize(1);
eV = &mySols[0];
// mySols.resize(1);
// eV = &mySols[0];
break;
default:
myMats->resize(0,0);
mySols.clear();
// mySols.clear();
tracVal.clear();
}
}
@ -160,7 +160,8 @@ void Elasticity::setTraction (TractionFunc* tf)
}
bool Elasticity::initElement (const std::vector<int>& MNPC)
bool Elasticity::initElement (const std::vector<int>& MNPC,
LocalIntegral& elmInt)
{
if (myMats)
myMats->withLHS = true;
@ -172,18 +173,19 @@ bool Elasticity::initElement (const std::vector<int>& MNPC)
if (eM) eM->resize(nsd*nen,nsd*nen,true);
if (eS) eS->resize(nsd*nen,true);
return this->IntegrandBase::initElement(MNPC);
return this->IntegrandBase::initElement(MNPC,elmInt);
}
bool Elasticity::initElementBou (const std::vector<int>& MNPC)
bool Elasticity::initElementBou (const std::vector<int>& MNPC,
LocalIntegral& elmInt)
{
if (myMats)
myMats->withLHS = false;
if (eS) eS->resize(nsd*MNPC.size(),true);
return this->IntegrandBase::initElementBou(MNPC);
return this->IntegrandBase::initElementBou(MNPC,elmInt);
}
@ -662,15 +664,17 @@ size_t ElasticityNorm::getNoFields () const
}
bool ElasticityNorm::initElement (const std::vector<int>& MNPC)
bool ElasticityNorm::initElement (const std::vector<int>& MNPC,
LocalIntegral& elmInt)
{
// Extract projected solution vectors for this element
int ierr = 0;
for (size_t i = 0; i < mySols.size() && ierr == 0; i++)
elmInt.vec.resize(prjsol.size());
for (size_t i = 0; i < prjsol.size() && ierr == 0; i++)
if (!prjsol[i].empty())
ierr = utl::gather(MNPC,nrcmp,prjsol[i],mySols[i]);
ierr = utl::gather(MNPC,nrcmp,prjsol[i],elmInt.vec[i]);
if (ierr == 0) return myProblem.initElement(MNPC);
// if (ierr == 0) return myProblem.initElement(MNPC);
std::cerr <<" *** ElasticityNorm::initElement: Detected "
<< ierr <<" node numbers out of range."<< std::endl;
@ -729,14 +733,14 @@ bool ElasticityNorm::evalInt (LocalIntegral*& elmInt, const FiniteElement& fe,
}
size_t i, j, k;
for (i = 0; i < mySols.size(); i++)
if (!mySols[i].empty())
for (i = 0; i < 0; i++) // mySols.size(); i++)
if (0) //!mySols[i].empty()) // TODO!
{
// Evaluate projected stress field
Vector sigmar(sigmah.size());
for (j = k = 0; j < nrcmp && k < sigmar.size(); j++)
if (!planeStrain || j != 2)
sigmar[k++] = mySols[i].dot(fe.N,j,nrcmp);
sigmar[k++] = 0;//mySols[i].dot(fe.N,j,nrcmp);
// Integrate the energy norm a(u^r,u^r)
pnorm[ip++] += sigmar.dot(Cinv*sigmar)*detJW;

11
src/Integrands/Elasticity.h
View File

@ -68,10 +68,15 @@ public:
//! \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);
//! \param[in] elmInt LocalIntegral for element
virtual bool initElement(const std::vector<int>& MNPC,
LocalIntegral& elmInt);
//! \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);
//! \param[in] elmInt LocalIntegral for element
virtual bool initElementBou(const std::vector<int>& MNPC,
LocalIntegral& elmInt);
//! \brief Evaluates the secondary solution at a result point.
//! \param[out] s Array of solution field values at current point
@ -249,7 +254,7 @@ public:
//! \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);
virtual bool initElement(const std::vector<int>& MNPC, LocalIntegral& elmInt);
//! \brief Evaluates the integrand at an interior point.
//! \param elmInt The local integral object to receive the contributions

35
src/Integrands/KirchhoffLovePlate.C
View File

@ -115,19 +115,20 @@ void KirchhoffLovePlate::setMode (SIM::SolutionMode mode)
case SIM::RECOVERY:
myMats->rhsOnly = true;
mySols.resize(1);
eV = &mySols[0];
// mySols.resize(1);
// eV = &mySols[0];
break;
default:
myMats->resize(0,0);
mySols.clear();
// mySols.clear();
presVal.clear();
}
}
bool KirchhoffLovePlate::initElement (const std::vector<int>& MNPC)
bool KirchhoffLovePlate::initElement (const std::vector<int>& MNPC,
LocalIntegral& elmInt)
{
if (myMats)
myMats->withLHS = true;
@ -138,18 +139,21 @@ bool KirchhoffLovePlate::initElement (const std::vector<int>& MNPC)
if (eM) eM->resize(nen,nen,true);
if (eS) eS->resize(nen,true);
return this->IntegrandBase::initElement(MNPC);
return false;
// return this->IntegrandBase::initElement(MNPC);
}
bool KirchhoffLovePlate::initElementBou (const std::vector<int>& MNPC)
bool KirchhoffLovePlate::initElementBou (const std::vector<int>& MNPC,
LocalIntegral& elmInt)
{
if (myMats)
myMats->withLHS = false;
if (eS) eS->resize(MNPC.size(),true);
return this->IntegrandBase::initElementBou(MNPC);
// return this->IntegrandBase::initElementBou(MNPC);
return false;
}
@ -464,15 +468,17 @@ size_t KirchhoffLovePlateNorm::getNoFields () const
}
bool KirchhoffLovePlateNorm::initElement (const std::vector<int>& MNPC)
bool KirchhoffLovePlateNorm::initElement (const std::vector<int>& MNPC,
LocalIntegral& elmInt)
{
// Extract projected solution vectors for this element
int ierr = 0;
for (size_t i = 0; i < mySols.size() && ierr == 0; i++)
elmInt.vec.resize(prjsol.size());
for (size_t i = 0; i < prjsol.size() && ierr == 0; i++)
if (!prjsol[i].empty())
ierr = utl::gather(MNPC,nrcmp,prjsol[i],mySols[i]);
ierr = utl::gather(MNPC,nrcmp,prjsol[i],elmInt.vec[i]);
if (ierr == 0) return myProblem.initElement(MNPC);
// if (ierr == 0) return myProblem.initElement(MNPC);
std::cerr <<" *** KirchhoffLovePlateNorm::initElement: Detected "
<< ierr <<" node numbers out of range."<< std::endl;
@ -522,13 +528,14 @@ bool KirchhoffLovePlateNorm::evalInt (LocalIntegral*& elmInt,
}
size_t i, j;
for (i = 0; i < mySols.size(); i++)
if (!mySols[i].empty())
// TODO!
for (i = 0; i < 0; i++) //mySols.size(); i++)
if (0) //!mySols[i].empty())
{
// Evaluate projected stress field
Vector mr(mh.size());
for (j = 0; j < nrcmp; j++)
mr[j] = mySols[i].dot(fe.N,j,nrcmp);
mr[j] = 0;//mySols[i].dot(fe.N,j,nrcmp);
// Integrate the energy norm a(w^r,w^r)
pnorm[ip++] += mr.dot(Cinv*mr)*fe.detJxW;

9
src/Integrands/KirchhoffLovePlate.h
View File

@ -61,10 +61,12 @@ public:
//! \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);
virtual bool initElement(const std::vector<int>& MNPC, LocalIntegral& elmInt);
//! \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);
virtual bool initElementBou(const std::vector<int>& MNPC,
LocalIntegral& elmInt);
//! \brief Defines which FE quantities are needed by the integrand.
virtual int getIntegrandType() const { return 2; }
@ -219,7 +221,8 @@ public:
//! \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);
virtual bool initElement(const std::vector<int>& MNPC,
LocalIntegral& elmInt);
//! \brief Defines which FE quantities are needed by the integrand.
virtual int getIntegrandType() const { return 2; }

206
src/Integrands/Poisson.C
View File

@ -21,6 +21,10 @@
#include "AnaSol.h"
#include "VTF.h"
#ifdef USE_OPENMP
#include <omp.h>
#endif
Poisson::Poisson (unsigned short int n) : nsd(n)
{
@ -30,56 +34,18 @@ Poisson::Poisson (unsigned short int n) : nsd(n)
tracFld = 0;
heatSrc = 0;
eM = 0;
eS = eV = 0;
// Only the current solution is needed
primsol.resize(1);
mySols.resize(1);
myMats = new ElmMats();
}
void Poisson::setMode (SIM::SolutionMode mode)
{
myMats->rhsOnly = false;
eM = 0;
eS = eV = 0;
m_mode = mode;
switch (mode)
{
case SIM::STATIC:
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];
break;
case SIM::RHS_ONLY:
myMats->rhsOnly = true;
if (myMats->b.empty())
myMats->b.resize(1);
eS = &myMats->b[0];
tracVal.clear();
break;
case SIM::RECOVERY:
myMats->rhsOnly = true;
eV = &mySols[0];
break;
default:
myMats->A.clear();
myMats->b.clear();
tracVal.clear();
}
if (mode != SIM::RECOVERY)
tracVal.clear();
}
@ -95,65 +61,81 @@ double Poisson::getTraction (const Vec3& X, const Vec3& n) const
}
void Poisson::setTraction (VecFunc* tf)
LocalIntegral* Poisson::getLocalIntegral (size_t nen, size_t,
bool neumann) const
{
tracFld = tf;
myMats->rhsOnly = true;
ElmMats* result = new ElmMats;
switch (m_mode)
{
case SIM::STATIC:
result->rhsOnly = neumann;
result->withLHS = !neumann;
result->resize(neumann?0:1,1);
break;
case SIM::VIBRATION:
result->withLHS = true;
result->resize(1,0);
break;
case SIM::RHS_ONLY:
result->rhsOnly = true;
result->resize(0,1);
break;
case SIM::RECOVERY:
result->rhsOnly = true;
break;
default:
break;
}
if (result->A.size())
result->A.front().resize(nen,nen);
if (result->b.size())
result->b.front().resize(nen);
return result;
}
bool Poisson::initElement (const std::vector<int>& MNPC)
{
const size_t nen = MNPC.size();
if (eM) eM->resize(nen,nen,true);
if (eS) eS->resize(nen,true);
return this->IntegrandBase::initElement(MNPC);
}
bool Poisson::initElementBou (const std::vector<int>& MNPC)
{
if (eS) eS->resize(MNPC.size(),true);
myMats->withLHS = false;
return true;
}
bool Poisson::evalInt (LocalIntegral*& elmInt, const FiniteElement& fe,
bool Poisson::evalInt (LocalIntegral& elmInt, const FiniteElement& fe,
const Vec3& X) const
{
elmInt = myMats;
ElmMats& elMat = static_cast<ElmMats&>(elmInt);
if (eM)
if (!elMat.A.empty())
{
// Evaluate the constitutive matrix at this point
if (!this->formCmatrix(C,X)) return false;
Matrix C, CB;
this->formCmatrix(C,X);
// Integrate the coefficient matrix
CB.multiply(C,fe.dNdX,false,true).multiply(fe.detJxW); // = C*dNdX^T*|J|*w
eM->multiply(fe.dNdX,CB,false,false,true); // EK += dNdX * CB
elMat.A.front().multiply(fe.dNdX,CB,false,false,true); // EK += dNdX * CB
}
// Integrate heat source, if defined
if (eS && heatSrc)
eS->add(fe.N,(*heatSrc)(X)*fe.detJxW);
if (heatSrc && !elMat.b.empty())
elMat.b.front().add(fe.N,(*heatSrc)(X)*fe.detJxW);
return true;
}
bool Poisson::evalBou (LocalIntegral*& elmInt, const FiniteElement& fe,
bool Poisson::evalBou (LocalIntegral& elmInt, const FiniteElement& fe,
const Vec3& X, const Vec3& normal) const
{
elmInt = myMats;
if (!tracFld)
{
std::cerr <<" *** Poisson::evalBou: No tractions."<< std::endl;
return false;
}
else if (!eS)
ElmMats& elMat = static_cast<ElmMats&>(elmInt);
if (elMat.b.empty())
{
std::cerr <<" *** Poisson::evalBou: No load vector."<< std::endl;
return false;
@ -167,10 +149,13 @@ bool Poisson::evalBou (LocalIntegral*& elmInt, const FiniteElement& fe,
// Store traction value for visualization
if (abs(trac) > 1.0e8)
tracVal.insert(std::make_pair(X,trac*normal));
#ifdef USE_OPENMP
if (omp_get_max_threads() == 1)
#endif
tracVal.insert(std::make_pair(X,trac*normal));
// Integrate the Neumann value
eS->add(fe.N,trac*fe.detJxW);
elMat.b.front().add(fe.N,trac*fe.detJxW);
return true;
}
@ -198,62 +183,46 @@ 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
{
if (primsol.front().empty())
{
std::cerr <<" *** Poisson::evalSol: No primary solution."
<< std::endl;
std::cerr <<" *** Poisson::evalSol: No primary solution."<< std::endl;
return false;
}
else if (!this->formCmatrix(C,X))
return false;
Vector Dtmp;
int ierr = utl::gather(MNPC,1,primsol.front(),Dtmp);
Vector eV;
int ierr = utl::gather(MNPC,1,primsol.front(),eV);
if (ierr > 0)
{
std::cerr <<" *** Poisson::evalSol: Detected "
<< ierr <<" node numbers out of range."<< std::endl;
std::cerr <<" *** Poisson::evalSol: Detected "<< ierr
<<" node numbers out of range."<< std::endl;
return false;
}
// Evaluate the heat flux vector
CB.multiply(C,dNdX,false,true).multiply(Dtmp,q); // q = C*dNdX^T*Dtmp
q *= -1.0;
return true;
return this->evalSol(q,eV,dNdX,X);
}
bool Poisson::evalSol (Vector& q, const Matrix& dNdX, const Vec3& X) const
bool Poisson::evalSol (Vector& q, const Vector& eV,
const Matrix& dNdX, const Vec3& X) const
{
if (!eV || eV->empty())
{
std::cerr <<" *** Poisson::evalSol: No solution vector."
<< std::endl;
return false;
}
else if (eV->size() != dNdX.rows())
if (eV.size() != dNdX.rows())
{
std::cerr <<" *** Poisson::evalSol: Invalid solution vector."
<<"\n size(eV) = "<< eV->size() <<" size(dNdX) = "
<< dNdX.rows() <<","<< dNdX.cols() << std::endl;
<<"\n size(eV) = "<< eV.size() <<" size(dNdX) = "
<< dNdX.rows() <<","<< dNdX.cols() << std::endl;
return false;
}
else if (!this->formCmatrix(C,X))
return false;
// Evaluate the constitutive matrix at this point
Matrix C, CB;
this->formCmatrix(C,X);
// Evaluate the heat flux vector
CB.multiply(C,dNdX,false,true).multiply(*eV,q); // q = C*dNdX^T*eV
CB.multiply(C,dNdX,false,true).multiply(eV,q); // q = C*dNdX^T*eV
q *= -1.0;
return true;
@ -302,24 +271,26 @@ NormBase* Poisson::getNormIntegrand (AnaSol* asol) const
}
bool PoissonNorm::evalInt (LocalIntegral*& elmInt, const FiniteElement& fe,
bool PoissonNorm::evalInt (LocalIntegral& elmInt, const FiniteElement& fe,
const Vec3& X) const
{
Poisson& problem = static_cast<Poisson&>(myProblem);
ElmNorm& pnorm = NormBase::getElmNormBuffer(elmInt);
ElmNorm& pnorm = static_cast<ElmNorm&>(elmInt);
// Evaluate the constitutive matrix at this point
// Evaluate the inverse constitutive matrix at this point
Matrix C;
if (!problem.formCmatrix(C,X,true)) return false;
problem.formCmatrix(C,X,true);
// Evaluate the finite element heat flux field
Vector sigmah;
if (!problem.evalSol(sigmah,fe.dNdX,X)) return false;
if (!problem.evalSol(sigmah,pnorm.vec.front(),fe.dNdX,X)) return false;
// Integrate the energy norm a(u^h,u^h)
pnorm[0] += sigmah.dot(C*sigmah)*fe.detJxW;
// Evaluate the temperature field
double u = pnorm.vec.front().dot(fe.N);
// Integrate the external energy (h,u^h)
pnorm[1] += problem.getHeat(X)*problem.evalSol(fe.N)*fe.detJxW;
pnorm[1] += problem.getHeat(X)*u*fe.detJxW;
if (anasol)
{
@ -336,18 +307,19 @@ bool PoissonNorm::evalInt (LocalIntegral*& elmInt, const FiniteElement& fe,
}
bool PoissonNorm::evalBou (LocalIntegral*& elmInt, const FiniteElement& fe,
bool PoissonNorm::evalBou (LocalIntegral& elmInt, const FiniteElement& fe,
const Vec3& X, const Vec3& normal) const
{
Poisson& problem = static_cast<Poisson&>(myProblem);
ElmNorm& pnorm = NormBase::getElmNormBuffer(elmInt);
ElmNorm& pnorm = static_cast<ElmNorm&>(elmInt);
// Evaluate the surface heat flux
double t = problem.getTraction(X,normal);
// Evaluate the temperature field
double u = problem.evalSol(fe.N);
double u = pnorm.vec.front().dot(fe.N);
// Integrate the external energy (t,u^h)
pnorm[1] += t*u*fe.detJxW;
return true;
}

41
src/Integrands/Poisson.h
View File

@ -20,6 +20,7 @@
class ElmNorm;
/*!
\brief Class representing the integrand of the Poisson problem.
\details This class supports constant isotropic constitutive properties only.
@ -40,7 +41,7 @@ public:
virtual void setMode(SIM::SolutionMode mode);
//! \brief Defines the traction field to use in Neumann boundary conditions.
void setTraction(VecFunc* tf);
void setTraction(VecFunc* tf) { tracFld = tf; }
//! \brief Clears the integration point traction values.
void clearTracVal() { tracVal.clear(); }
@ -55,18 +56,17 @@ public:
//! \brief Evaluates the heat source (if any) at specified point.
double getHeat(const Vec3& X) const;
//! \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 Returns a local integral container for the given element.
//! \param[in] nen Number of nodes on element
//! \param[in] neumann Whether or not we are assembling Neumann BC's
virtual LocalIntegral* getLocalIntegral(size_t nen, size_t,
bool neumann) const;
//! \brief 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] X Cartesian coordinates of current integration point
virtual bool evalInt(LocalIntegral*& elmInt, const FiniteElement& fe,
virtual bool evalInt(LocalIntegral& elmInt, const FiniteElement& fe,
const Vec3& X) const;
//! \brief Evaluates the integrand at a boundary point.
@ -74,14 +74,9 @@ public:
//! \param[in] fe Finite element data of current integration point
//! \param[in] X Cartesian coordinates of current integration point
//! \param[in] normal Boundary normal vector at current integration point
virtual bool evalBou(LocalIntegral*& elmInt, const FiniteElement& fe,
virtual bool evalBou(LocalIntegral& elmInt, const FiniteElement& fe,
const Vec3& X, const Vec3& normal) const;
//! \brief 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
@ -94,9 +89,11 @@ public:
//! \brief Evaluates the finite element (FE) solution at an integration point.
//! \param[out] s The FE solution values at current point
//! \param[in] eV Element solution vector
//! \param[in] dNdX Basis function gradients at current point
//! \param[in] X Cartesian coordinates of current point
virtual bool evalSol(Vector& s, const Matrix& dNdX, const Vec3& X) const;
virtual bool evalSol(Vector& s, const Vector& eV,
const Matrix& dNdX, const Vec3& X) const;
//! \brief Evaluates the analytical solution at an integration point.
//! \param[out] s The analytical solution values at current point
@ -142,22 +139,12 @@ private:
double kappa; //!< Conductivity
protected:
// Finite element quantities
Matrix* eM; //!< Pointer to element matrix
Vector* eS; //!< Pointer to element right-hand-side vector
Vector* eV; //!< Pointer to element solution vector
VecFunc* tracFld; //!< Pointer to boundary traction field
RealFunc* heatSrc; //!< Pointer to interior heat source
mutable std::map<Vec3,Vec3> tracVal; //!< Traction field point values
unsigned short int nsd; //!< Number of space dimensions (1, 2 or, 3)
// Work arrays declared as members to avoid frequent re-allocation
// within the numerical integration loop (for reduced overhead)
mutable Matrix C; //!< Constitutive matrix
mutable Matrix CB; //!< Result of the matrix-matrix product C*dNdX^T
};
@ -185,7 +172,7 @@ public:
//! \param elmInt The local integral object to receive the contributions
//! \param[in] fe Finite element data of current integration point
//! \param[in] X Cartesian coordinates of current integration point
virtual bool evalInt(LocalIntegral*& elmInt, const FiniteElement& fe,
virtual bool evalInt(LocalIntegral& elmInt, const FiniteElement& fe,
const Vec3& X) const;
//! \brief Evaluates the integrand at a boundary point.
@ -193,7 +180,7 @@ public:
//! \param[in] fe Finite Element quantities
//! \param[in] X Cartesian coordinates of current integration point
//! \param[in] normal Boundary normal vector at current integration point
virtual bool evalBou(LocalIntegral*& elmInt, const FiniteElement& fe,
virtual bool evalBou(LocalIntegral& elmInt, const FiniteElement& fe,
const Vec3& X, const Vec3& normal) const;
private:

12
src/SIM/SIM2D.C
View File

@ -18,6 +18,10 @@
#include "tinyxml.h"
#ifdef USE_OPENMP
#include <omp.h>
#endif
/*!
A struct defining a patch interface for C1-continuous models.
@ -199,6 +203,10 @@ bool SIM2D::parseGeometryTag(const TiXmlElement* elem)
std::cout <<"\tPeriodic "<< char('H'+pedir) <<"-direction P"<< patch
<< std::endl;
static_cast<ASMs2D*>(myModel[patch-1])->closeEdges(pedir);
// cannot do multi-threaded assembly with periodicities
#ifdef USE_OPENMP
omp_set_num_threads(1);
#endif
}
// TODO: constraints? fixpoints?
@ -382,6 +390,10 @@ bool SIM2D::parse (char* keyWord, std::istream& is)
std::cout <<"\tPeriodic "<< char('H'+pedir) <<"-direction P"<< patch
<< std::endl;
static_cast<ASMs2D*>(myModel[patch-1])->closeEdges(pedir);
// cannot do multi-threaded assembly with periodicities
#ifdef USE_OPENMP
omp_set_num_threads(1);
#endif
}
}

8
src/SIM/SIM3D.C
View File

@ -211,6 +211,10 @@ bool SIM3D::parse (char* keyWord, std::istream& is)
std::cout <<"\tPeriodic "<< char('H'+pfdir) <<"-direction P"<< patch
<< std::endl;
static_cast<ASMs3D*>(myModel[patch-1])->closeFaces(pfdir);
// cannot do multi-threaded assembly with periodicities
#ifdef USE_OPENMP
omp_set_num_threads(1);
#endif
}
}
@ -412,6 +416,10 @@ bool SIM3D::parseGeometryTag(const TiXmlElement* elem)
std::cout <<"\tPeriodic "<< char('H'+pfdir) <<"-direction P"<< patch
<< std::endl;
static_cast<ASMs3D*>(myModel[patch-1])->closeFaces(pfdir);
// cannot do multi-threaded assembly with periodicities
#ifdef USE_OPENMP
omp_set_num_threads(1);
#endif
}
// TODO: topology? constraints? fixpoints?

5
src/SIM/SIMbase.C
View File

@ -1187,7 +1187,7 @@ bool SIMbase::solutionNorms (const TimeDomain& time,
for (k = 0; k < ssol.size(); k++)
if (!ssol[k].empty())
myModel[j-1]->extractNodeVec(ssol[k],norm->getProjection(k+1),nCmp);
ok = myModel[j-1]->integrate(*norm,globalNorm,time,elementNorms);
ok = myModel[j-1]->integrate(*norm,globalNorm,time);
lp = j;
}
else
@ -1202,7 +1202,7 @@ bool SIMbase::solutionNorms (const TimeDomain& time,
for (k = 0; k < ssol.size(); k++)
if (!ssol[k].empty())
myModel[i]->extractNodeVec(ssol[k],norm->getProjection(k+1),nCmp);
ok = myModel[i]->integrate(*norm,globalNorm,time,elementNorms);
ok = myModel[i]->integrate(*norm,globalNorm,time);
lp = i+1;
}
@ -1211,6 +1211,7 @@ bool SIMbase::solutionNorms (const TimeDomain& time,
// The corresponding element-level norms are not stored. This is mainly
// because the current design only allows one loop over the elements
// in the element-norm calculations. Consider rivising this later.
if (norm->hasBoundaryTerms())
for (i = 0; i < myProps.size() && ok; i++)
if (myProps[i].pcode == Property::NEUMANN)