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

added: COMPATIBLE integrand

Browse Source
This commit is contained in:
Timo van Opstal
2015-12-17 16:30:51 +01:00
committed by Knut Morten Okstad
parent 0d08357cb7
commit 6b3954c502
14 changed files with 75 additions and 40 deletions
Download Patch File Download Diff File Expand all files Collapse all files

18
Apps/Common/WeakOperators.C
View File

@@ -140,10 +140,10 @@ namespace WeakOperators {
void Laplacian(Matrix& EM, const FiniteElement& fe,
double scale, size_t cmp, size_t nf,
bool stress, size_t scmp)
bool stress, size_t scmp, unsigned char basis)
{
Matrix A;
A.multiply(fe.dNdX,fe.dNdX,false,true);
A.multiply(fe.grad(basis),fe.grad(basis),false,true);
A *= scale*fe.detJxW;
addComponents(EM, A, cmp, nf, scmp);
if (stress)
@@ -151,7 +151,7 @@ namespace WeakOperators {
for (size_t j = 1; j <= fe.N.size(); j++)
for (size_t k = 1; k <= cmp; k++)
for (size_t l = 1; l <= cmp; l++)
EM(nf*(j-1)+k+scmp,nf*(i-1)+l+scmp) += scale*fe.dNdX(i,k)*fe.dNdX(j,l)*fe.detJxW;
EM(nf*(j-1)+k+scmp,nf*(i-1)+l+scmp) += scale*fe.grad(basis)(i,k)*fe.grad(basis)(j,l)*fe.detJxW;
}
@@ -166,24 +166,26 @@ namespace WeakOperators {
void Mass(Matrix& EM, const FiniteElement& fe,
double scale, size_t cmp, size_t nf, size_t scmp)
double scale, size_t cmp, size_t nf, size_t scmp,
unsigned char basis)
{
Matrix A;
A.outer_product(fe.N,fe.N,false);
A.outer_product(fe.basis(basis),fe.basis(basis),false);
A *= scale*fe.detJxW;
addComponents(EM, A, cmp, nf, scmp);
}
void Source(Vector& EV, const FiniteElement& fe,
double scale, size_t cmp, size_t nf, size_t scmp)
double scale, size_t cmp, size_t nf, size_t scmp,
unsigned char basis)
{
if (cmp == 1 && nf == 1)
EV.add(fe.N, scale*fe.detJxW);
EV.add(fe.basis(basis), scale*fe.detJxW);
else {
for (size_t i = 1; i <= fe.N.size(); ++i)
for (size_t k = 1; k <= cmp; ++k)
EV(nf*(i-1)+k+scmp) += scale*fe.N(i)*fe.detJxW;
EV(nf*(i-1)+k+scmp) += scale*fe.basis(basis)(i)*fe.detJxW;
}
}

11
Apps/Common/WeakOperators.h
View File

@@ -95,9 +95,10 @@ namespace WeakOperators
//! \param[in] nf Number of fields in basis.
//! \param[in] stress Whether to add extra stress formulation terms.
//! \param[in] scmp Starting component.
//! \param[in] basis Basis to use.
void Laplacian(Matrix& EM, const FiniteElement& fe,
double scale=1.0, size_t cmp=1, size_t nf=1,
bool stress=false, size_t scmp=0);
bool stress=false, size_t scmp=0, unsigned char basis=1);
//! \brief Compute a heteregenous coefficient laplacian.
//! \param[out] EM The element matrix to add contribution to.
@@ -114,8 +115,10 @@ namespace WeakOperators
//! \param[in] cmp Number of components to add.
//! \param[in] nf Number of fields in basis.
//! \param[in] scmp Starting component.
//! \param[in] basis Basis to use.
void Mass(Matrix& EM, const FiniteElement& fe,
double scale=1.0, size_t cmp=1, size_t nf=1, size_t scmp=0);
double scale=1.0, size_t cmp=1, size_t nf=1, size_t scmp=0,
unsigned char basis=1);
//! \brief Compute a source term.
//! \param[out] EV The element vector to add contribution to.
@@ -124,8 +127,10 @@ namespace WeakOperators
//! \param[in] cmp Number of components to add.
//! \param[in] nf Number of fields in basis.
//! \param[in] scmp Starting component.
//! \param[in] basis Basis to use.
void Source(Vector& EV, const FiniteElement& fe,
double scale=1.0, size_t cmp=1, size_t nf=1, size_t scmp=0);
double scale=1.0, size_t cmp=1, size_t nf=1, size_t scmp=0,
unsigned char basis=1);
//! \brief Compute a vector-source term.
//! \param[out] EV The element vector to add contribution to.

11
src/ASM/ASMbase.h
View File

@@ -556,16 +556,17 @@ public:
//! \param[in] code Identifier for inhomogeneous Dirichlet condition field
bool add2PC(int slave, int dir, int master, int code = 0);
protected:
// Internal methods for preprocessing of boundary conditions
// =========================================================
//! \brief Adds a general multi-point-constraint (MPC) equation to this patch.
//! \param mpc Pointer to an MPC-object
//! \param[in] code Identifier for inhomogeneous Dirichlet condition field
//! \param[in] silence If \e true, suppress debug print
bool addMPC(MPC*& mpc, int code = 0, bool silence = false);
protected:
// Internal methods for preprocessing of boundary conditions
// =========================================================
//! \brief Creates and adds a three-point constraint to this patch.
//! \param[in] slave Global node number of the node to constrain
//! \param[in] dir Which local DOF to constrain (1, 2, 3)

14
src/ASM/ASMs2Dmx.C
View File

@@ -929,18 +929,24 @@ bool ASMs2Dmx::evalSolution (Matrix& sField, const IntegrandBase& integrand,
// Fetch indices of the non-zero basis functions at this point
std::vector<IntVec> ip(m_basis.size());
IntVec ipa;
size_t ofs = 0;
for (size_t b = 0; b < m_basis.size(); ++b) {
scatterInd(m_basis[b]->numCoefs_u(),m_basis[b]->numCoefs_v(),
m_basis[b]->order_u(),m_basis[b]->order_v(),splinex[b][i].left_idx,ip[b]);
// Fetch associated control point coordinates
if (b == (size_t)geoBasis-1)
utl::gather(ip[geoBasis-1], nsd, Xnod, Xtmp);
for (auto& it : ip[b])
it += ofs;
ipa.insert(ipa.end(), ip[b].begin(), ip[b].end());
ofs += nb[b];
}
fe.u = splinex[0][i].param[0];
fe.v = splinex[0][i].param[1];
// Fetch associated control point coordinates
utl::gather(ip[geoBasis-1], nsd, Xnod, Xtmp);
// Fetch basis function derivatives at current integration point
for (size_t b = 0; b < m_basis.size(); ++b)
SplineUtils::extractBasis(splinex[b][i],fe.basis(b+1),dNxdu[b]);
@@ -962,7 +968,7 @@ bool ASMs2Dmx::evalSolution (Matrix& sField, const IntegrandBase& integrand,
X = Xtmp * fe.basis(geoBasis);
// Now evaluate the solution field
if (!integrand.evalSol(solPt,fe,X,ipa,elem_sizes))
if (!integrand.evalSol(solPt,fe,X,ipa,elem_sizes,nb))
return false;
else if (sField.empty())
sField.resize(solPt.size(),nPoints,true);

2
src/ASM/ASMs2DmxLag.C
View File

@@ -552,7 +552,7 @@ bool ASMs2DmxLag::evalSolution (Matrix& sField, const IntegrandBase& integrand,
}
// Now evaluate the solution field
if (!integrand.evalSol(solPt,fe,Xnod*fe.basis(geoBasis),MNPC[iel-1],elem_size))
if (!integrand.evalSol(solPt,fe,Xnod*fe.basis(geoBasis),MNPC[iel-1],elem_size,nb))
return false;
else if (sField.empty())
sField.resize(solPt.size(),nPoints,true);

13
src/ASM/ASMs3Dmx.C
View File

@@ -286,7 +286,8 @@ bool ASMs3Dmx::generateFEMTopology ()
iel = 0;
if ((int)b != geoBasis-1)
addBasis(m_basis[b],false);
inod += nb[b];
else
inod += nb[b];
lnod2 += m_basis[b]->order(0)*m_basis[b]->order(1)*m_basis[b]->order(2);
}
@@ -391,10 +392,10 @@ Vec3 ASMs3Dmx::getCoord (size_t inod) const
const int J = nodeInd[inod-1].J;
const int K = nodeInd[inod-1].K;
int b = 0;
size_t nbb = 0;
while (inod < nbb)
nbb += nb[b++];
size_t b = 0;
size_t nbb = nb[0];
while (nbb < inod)
nbb += nb[++b];
cit = m_basis[b]->coefs_begin()
+ ((K*m_basis[b]->numCoefs(1)+J)*m_basis[b]->numCoefs(0)+I) * m_basis[b]->dimension();
@@ -930,7 +931,7 @@ bool ASMs3Dmx::evalSolution (Matrix& sField, const IntegrandBase& integrand,
X = Xtmp * fe.basis(geoBasis);
// Now evaluate the solution field
if (!integrand.evalSol(solPt,fe,X,ipa,elem_sizes))
if (!integrand.evalSol(solPt,fe,X,ipa,elem_sizes,nb))
return false;
else if (sField.empty())
sField.resize(solPt.size(),nPoints,true);

2
src/ASM/ASMs3DmxLag.C
View File

@@ -622,7 +622,7 @@ bool ASMs3DmxLag::evalSolution (Matrix& sField, const IntegrandBase& integrand,
}
// Now evaluate the solution field
if (!integrand.evalSol(solPt,fe,Xnod*fe.basis(geoBasis),MNPC[iel-1],elem_size))
if (!integrand.evalSol(solPt,fe,Xnod*fe.basis(geoBasis),MNPC[iel-1],elem_size,nb))
return false;
else if (sField.empty())
sField.resize(solPt.size(),nPoints,true);

6
src/ASM/GlbForceVec.C
View File

@@ -53,19 +53,19 @@ bool GlbForceVec::assemble (const LocalIntegral* elmObj, int elmId)
const ElmMats* elm = dynamic_cast<const ElmMats*>(elmObj);
if (!elm || elm->b.size() < 1) return false;
const Vector& ES = elm->b.front();
const Vector& ES = elm->getRHSVector();
const size_t nfc = F.rows();
std::vector<int> mnpc;
if (!sam.getElmNodes(mnpc,elmId))
return false;
else if (ES.size() < nfc*mnpc.size())
/* else if (ES.size() < nfc*mnpc.size())
{
std::cerr <<" *** GlbForceVec::assemble: Invalid element force vector,"
<<" size="<< ES.size() <<" should be (at least) "
<< nfc*mnpc.size() << std::endl;
return false;
}
*/
// Assemble the nodal forces into the Matrix F
size_t i, j, k, ninod = 0;
std::map<int,size_t>::const_iterator nit;

8
src/ASM/GlbL2projector.C
View File

@@ -40,7 +40,8 @@ public:
GlbL2& gl2Int; //!< The global L2 projection integrand
LocalIntegral* elmData; //!< Element data associated with problem integrand
IntVec mnpc; //!< Matrix of element nodal correspondance for bases
std::vector<size_t> sizes; //!< Size of each basis
std::vector<size_t> elem_sizes; //!< Size of each basis on the element
std::vector<size_t> basis_sizes; //!< Size of each basis on the patch
};
@@ -85,7 +86,8 @@ bool GlbL2::initElement (const IntVec& MNPC1,
L2Mats& gl2 = static_cast<L2Mats&>(elmInt);
gl2.mnpc = MNPC1;
gl2.sizes = elem_sizes;
gl2.elem_sizes = elem_sizes;
gl2.basis_sizes = basis_sizes;
return problem.initElement(MNPC1,elem_sizes,basis_sizes,*gl2.elmData);
}
@@ -127,7 +129,7 @@ bool GlbL2::evalIntMx (LocalIntegral& elmInt,
L2Mats& gl2 = static_cast<L2Mats&>(elmInt);
Vector solPt;
if (!problem.evalSol(solPt,fe,X,gl2.mnpc,gl2.sizes))
if (!problem.evalSol(solPt,fe,X,gl2.mnpc,gl2.elem_sizes,gl2.basis_sizes))
if (!problem.diverged(fe.iGP+1))
return false;

3
src/ASM/IntegrandBase.C
View File

@@ -140,7 +140,8 @@ bool IntegrandBase::evalSol (Vector& s, const FiniteElement& fe,
bool IntegrandBase::evalSol (Vector& s, const MxFiniteElement& fe,
const Vec3& X, const std::vector<int>& MNPC,
const std::vector<size_t>& elem_sizes) const
const std::vector<size_t>& elem_sizes,
const std::vector<size_t>& basis_sizes) const
{
std::vector<int> MNPC1(MNPC.begin(), MNPC.begin()+elem_sizes.front());
return this->evalSol(s,fe,X,MNPC1);

6
src/ASM/IntegrandBase.h
View File

@@ -59,6 +59,8 @@ public:
//! \brief Defines the solution mode before the element assembly is started.
virtual void setMode(SIM::SolutionMode mode) { m_mode = mode; }
//! \brief Obtain current solution mode
SIM::SolutionMode getMode() const { return m_mode; }
//! \brief Initializes an integration parameter for the integrand.
virtual void setIntegrationPrm(unsigned short int, double) {}
//! \brief Returns an integration parameter for the integrand.
@@ -158,9 +160,11 @@ public:
//! \param[in] X Cartesian coordinates of current point
//! \param[in] MNPC Nodal point correspondance for the bases
//! \param[in] elem_sizes Size of each basis on the element
//! \param[in] basis_sizes Size of each basis on the patch
virtual bool evalSol(Vector& s, const MxFiniteElement& fe, const Vec3& X,
const std::vector<int>& MNPC,
const std::vector<size_t>& elem_sizes) const;
const std::vector<size_t>& elem_sizes,
const std::vector<size_t>& basis_sizes) const;
//! \brief Evaluates the analytical secondary solution at a result point.
//! \param[out] s The solution field values at current point

11
src/SIM/ForceIntegrator.C
View File

@@ -103,6 +103,17 @@ bool SIM::integrate(const Vectors& solution, SIMbase* model, int code,
size_t prevPatch = 0;
GlobalIntegral dummy;
GlobalIntegral& frc = force ? *force : dummy;
if (code == 0) { // special case - volume integral over the entire model
for (int p = 1; p <= model->getNoPatches() && ok; ++p) {
ASMbase* patch = model->getPatch(p);
ok = model->extractPatchSolution(solution,p-1);
model->setPatchMaterial(p);
ok &= patch->integrate(*forceInt,frc,time);
}
return ok;
}
PropertyVec::const_iterator p;
for (p = model->begin_prop(); p != model->end_prop() && ok; p++)
if (abs(p->pindx) == code)

9
src/SIM/SIMbase.C
View File

@@ -309,6 +309,8 @@ bool SIMbase::parseBCTag (const TiXmlElement* elem)
std::string set, type;
utl::getAttribute(elem,"set",set);
utl::getAttribute(elem,"type",type,true);
int ndir = 0;
utl::getAttribute(elem,"direction",ndir);
int code = this->getUniquePropertyCode(set);
if (code == 0) utl::getAttribute(elem,"code",code);
if (type == "anasol") {
@@ -319,14 +321,12 @@ bool SIMbase::parseBCTag (const TiXmlElement* elem)
IFEM::cout <<"\tNeumann code "<< code <<" (generic)";
if (elem->FirstChild()) {
this->setPropertyType(code,Property::ROBIN);
this->setNeumann(elem->FirstChild()->Value(),"expression",0,code);
this->setNeumann(elem->FirstChild()->Value(),"expression",ndir,code);
}
else
this->setPropertyType(code,Property::NEUMANN_GENERIC);
}
else {
int ndir = 0;
utl::getAttribute(elem,"direction",ndir);
IFEM::cout <<"\tNeumann code "<< code <<" direction "<< ndir;
if (!type.empty()) IFEM::cout <<" ("<< type <<")";
if (elem->FirstChild())
@@ -1044,7 +1044,8 @@ bool SIMbase::preprocess (const IntVec& ignored, bool fixDup)
(*mit)->generateThreadGroups(*myProblem,silence);
for (q = myProps.begin(); q != myProps.end(); q++)
if (q->pcode == Property::NEUMANN || q->pcode == Property::NEUMANN_GENERIC)
if (q->pcode == Property::NEUMANN || q->pcode == Property::NEUMANN_GENERIC ||
q->pcode == Property::ROBIN)
this->generateThreadGroups(*q,silence);
// Preprocess the result points

1
src/Utility/HDF5Writer.C
View File

@@ -465,6 +465,7 @@ void HDF5Writer::writeSIM (int level, const DataEntry& entry,
if (abs(results) & DataExporter::SECONDARY) {
Matrix field;
if (prefix.empty()) {
sim->setMode(SIM::RECOVERY);
sim->extractPatchSolution(Vectors(1,*sol), loc-1);
sim->evalSecondarySolution(field,loc-1);
}