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

changed: rename ASMmxBase::geoBasis to ASMmxBase::itgBasis

Browse Source 
more precise as this represents the basis holding the integration elements
This commit is contained in:
Arne Morten Kvarving 2023-08-21 12:24:24 +02:00
parent 6aeefcd0f0
commit 3ad0b872c8
18 changed files with 181 additions and 181 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

22
src/ASM/ASMmxBase.C
View File

@ -21,7 +21,7 @@
#include <numeric>
char ASMmxBase::geoBasis = 2;
char ASMmxBase::itgBasis = 2;
ASMmxBase::MixedType ASMmxBase::Type = ASMmxBase::FULL_CONT_RAISE_BASIS1;
@ -135,7 +135,7 @@ ASMmxBase::SurfaceVec ASMmxBase::establishBases (Go::SplineSurface* surf,
// basis1 should be one degree higher than basis2 and C^p-1 continuous
result[0].reset(ASMmxBase::raiseBasis(surf));
result[1].reset(new Go::SplineSurface(*surf));
geoBasis = 2;
itgBasis = 2;
}
else if (type == REDUCED_CONT_RAISE_BASIS1 || type == REDUCED_CONT_RAISE_BASIS2)
{
@ -144,7 +144,7 @@ ASMmxBase::SurfaceVec ASMmxBase::establishBases (Go::SplineSurface* surf,
result[0].reset(new Go::SplineSurface(*surf));
result[0]->raiseOrder(1,1);
result[1].reset(new Go::SplineSurface(*surf));
geoBasis = 2;
itgBasis = 2;
}
else if (ASMmxBase::Type == ASMmxBase::DIV_COMPATIBLE)
{
@ -182,7 +182,7 @@ ASMmxBase::SurfaceVec ASMmxBase::establishBases (Go::SplineSurface* surf,
result[1].reset(Go::SurfaceInterpolator::regularInterpolation(a1,b2,
u0,vg,XYZ1,ndim,
false,XYZ1));
geoBasis = 3;
itgBasis = 3;
} else if (type == SUBGRID) {
// basis1 should be one degree higher than basis2 and C^p-1 continuous
int ndim = surf->dimension();
@ -242,11 +242,11 @@ ASMmxBase::SurfaceVec ASMmxBase::establishBases (Go::SplineSurface* surf,
ug,vg,XYZ,ndim,
false,XYZ));
}
geoBasis = 1;
itgBasis = 1;
}
if (type == FULL_CONT_RAISE_BASIS2 || type == REDUCED_CONT_RAISE_BASIS2)
std::swap(result[0], result[1]), geoBasis = 1;
std::swap(result[0], result[1]), itgBasis = 1;
return result;
}
@ -262,7 +262,7 @@ ASMmxBase::VolumeVec ASMmxBase::establishBases(Go::SplineVolume* svol,
// basis1 should be one degree higher than basis2 and C^p-1 continuous
result[0].reset(ASMmxBase::raiseBasis(svol));
result[1].reset(new Go::SplineVolume(*svol));
geoBasis = 2;
itgBasis = 2;
}
else if (type == REDUCED_CONT_RAISE_BASIS1 || type == REDUCED_CONT_RAISE_BASIS2)
{
@ -271,7 +271,7 @@ ASMmxBase::VolumeVec ASMmxBase::establishBases(Go::SplineVolume* svol,
result[0].reset(new Go::SplineVolume(*svol));
result[0]->raiseOrder(1,1,1);
result[1].reset(new Go::SplineVolume(*svol));
geoBasis = 2;
itgBasis = 2;
}
else if (ASMmxBase::Type == ASMmxBase::DIV_COMPATIBLE)
{
@ -319,7 +319,7 @@ ASMmxBase::VolumeVec ASMmxBase::establishBases(Go::SplineVolume* svol,
u0,v0,wg,XYZ2,ndim,
false,XYZ2));
result[3].reset(new Go::SplineVolume(*svol));
geoBasis = 4;
itgBasis = 4;
} else if (type == SUBGRID) {
// basis1 should be one degree higher than basis2 and C^p-1 continuous
int ndim = svol->dimension();
@ -380,11 +380,11 @@ ASMmxBase::VolumeVec ASMmxBase::establishBases(Go::SplineVolume* svol,
XYZ,ndim,
false,XYZ));
}
geoBasis = 1;
itgBasis = 1;
}
if (type == FULL_CONT_RAISE_BASIS2 || type == REDUCED_CONT_RAISE_BASIS2)
std::swap(result[0], result[1]), geoBasis = 1;
std::swap(result[0], result[1]), itgBasis = 1;
return result;
}

2
src/ASM/ASMmxBase.h
View File

@ -73,7 +73,7 @@ public:
};
static MixedType Type; //!< Type of mixed formulation used
static char geoBasis; //!< 1-based index of basis representing the geometry
static char itgBasis; //!< 1-based index of basis representing the integration elements
protected:
typedef std::vector<std::shared_ptr<Go::SplineSurface>> SurfaceVec; //!< Convenience type

50
src/ASM/ASMs2Dmx.C
View File

@ -217,7 +217,7 @@ bool ASMs2Dmx::generateFEMTopology ()
else if (ASMmxBase::Type == ASMmxBase::SUBGRID)
projB = m_basis.front().get();
else // FULL_CONT_RAISE_BASISx
projB = m_basis[2-geoBasis].get();
projB = m_basis[2-itgBasis].get();
}
if (ASMmxBase::Type == ASMmxBase::SUBGRID)
@ -225,7 +225,7 @@ bool ASMs2Dmx::generateFEMTopology ()
delete surf;
}
geomB = surf = m_basis[geoBasis-1].get();
geomB = surf = m_basis[itgBasis-1].get();
nb.clear();
nb.reserve(m_basis.size());
@ -299,7 +299,7 @@ bool ASMs2Dmx::generateFEMTopology ()
myMNPC[iel][elm_node++] = last_node - spline.numCoefs_u()*j - i;
};
inod = std::accumulate(nb.begin(),nb.begin()+geoBasis-1,0u);
inod = std::accumulate(nb.begin(),nb.begin()+itgBasis-1,0u);
auto knotv = itg->basis_v().begin();
for (int i2 = 1; i2 <= itg->numCoefs_v(); i2++, ++knotv) {
auto knotu = itg->basis_u().begin();
@ -316,7 +316,7 @@ bool ASMs2Dmx::generateFEMTopology ()
size_t basis_ofs = 0;
int basis = 0;
for (const std::shared_ptr<Go::SplineSurface>& spline : m_basis) {
if (basis != geoBasis-1) {
if (basis != itgBasis-1) {
double ku = *knotu;
double kv = *knotv;
int iu = spline->basis_u().knotIntervalFuzzy(ku);
@ -541,11 +541,11 @@ bool ASMs2Dmx::integrate (Integrand& integrand,
// Compute Jacobian inverse of the coordinate mapping and
// basis function derivatives w.r.t. Cartesian coordinates
if (!fe.Jacobian(Jac,Xnod,geoBasis,&bfs))
if (!fe.Jacobian(Jac,Xnod,itgBasis,&bfs))
continue; // skip singular points
// Compute Hessian of coordinate mapping and 2nd order derivatives
if (use2ndDer && !fe.Hessian(Hess,Jac,Xnod,geoBasis,&bfs))
if (use2ndDer && !fe.Hessian(Hess,Jac,Xnod,itgBasis,&bfs))
ok = false;
// Compute G-matrix
@ -553,7 +553,7 @@ bool ASMs2Dmx::integrate (Integrand& integrand,
utl::getGmat(Jac,dXidu,fe.G);
// Cartesian coordinates of current integration point
X.assign(Xnod * fe.basis(geoBasis));
X.assign(Xnod * fe.basis(itgBasis));
// Evaluate the integrand and accumulate element contributions
fe.detJxW *= dA*wg[0][i]*wg[1][j];
@ -710,18 +710,18 @@ bool ASMs2Dmx::integrate (Integrand& integrand, int lIndex,
// Compute Jacobian inverse of the coordinate mapping and
// basis function derivatives w.r.t. Cartesian coordinates
fe.detJxW = utl::Jacobian(Jac,normal,fe.grad(geoBasis),Xnod,
dNxdu[geoBasis-1],t1,t2);
fe.detJxW = utl::Jacobian(Jac,normal,fe.grad(itgBasis),Xnod,
dNxdu[itgBasis-1],t1,t2);
if (fe.detJxW == 0.0) continue; // skip singular points
for (size_t b = 0; b < m_basis.size(); ++b)
if (b != (size_t)geoBasis-1)
if (b != (size_t)itgBasis-1)
fe.grad(b+1).multiply(dNxdu[b],Jac);
if (edgeDir < 0) normal *= -1.0;
// Cartesian coordinates of current integration point
X.assign(Xnod * fe.basis(geoBasis));
X.assign(Xnod * fe.basis(itgBasis));
// Evaluate the integrand and accumulate element contributions
fe.detJxW *= dS*wg[i];
@ -888,13 +888,13 @@ bool ASMs2Dmx::integrate (Integrand& integrand,
}
// Compute basis function derivatives and the edge normal
fe.detJxW = utl::Jacobian(Jac,normal,fe.grad(geoBasis+m_basis.size()),Xnod,
dNxdu[geoBasis-1+m_basis.size()],t1,t2);
fe.detJxW = utl::Jacobian(Jac,normal,fe.grad(geoBasis),Xnod,
dNxdu[geoBasis-1],t1,t2);
fe.detJxW = utl::Jacobian(Jac,normal,fe.grad(itgBasis+m_basis.size()),Xnod,
dNxdu[itgBasis-1+m_basis.size()],t1,t2);
fe.detJxW = utl::Jacobian(Jac,normal,fe.grad(itgBasis),Xnod,
dNxdu[itgBasis-1],t1,t2);
if (fe.detJxW == 0.0) continue; // skip singular points
for (size_t b = 0; b < m_basis.size(); ++b)
if (b != (size_t)geoBasis-1) {
if (b != (size_t)itgBasis-1) {
fe.grad(b+1).multiply(dNxdu[b],Jac);
fe.grad(b+1+m_basis.size()).multiply(dNxdu[b+m_basis.size()],Jac);
}
@ -902,7 +902,7 @@ bool ASMs2Dmx::integrate (Integrand& integrand,
if (edgeDir < 0) normal *= -1.0;
// Cartesian coordinates of current integration point
X.assign(Xnod * fe.basis(geoBasis));
X.assign(Xnod * fe.basis(itgBasis));
// Evaluate the integrand and accumulate element contributions
fe.detJxW *= dS*wg[i];
@ -1053,8 +1053,8 @@ bool ASMs2Dmx::evalSolution (Matrix& sField, const IntegrandBase& integrand,
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);
if (b == (size_t)itgBasis-1)
utl::gather(ip[itgBasis-1], nsd, Xnod, Xtmp);
for (int& c : ip[b]) c += ofs;
ipa.insert(ipa.end(), ip[b].begin(), ip[b].end());
@ -1070,11 +1070,11 @@ bool ASMs2Dmx::evalSolution (Matrix& sField, const IntegrandBase& integrand,
// Compute Jacobian inverse of the coordinate mapping and
// basis function derivatives w.r.t. Cartesian coordinates
if (!fe.Jacobian(Jac,Xtmp,geoBasis,nullptr,&dNxdu))
if (!fe.Jacobian(Jac,Xtmp,itgBasis,nullptr,&dNxdu))
continue; // skip singular points
// Cartesian coordinates of current integration point
utl::Point X4(Xtmp*fe.basis(geoBasis),{fe.u,fe.v});
utl::Point X4(Xtmp*fe.basis(itgBasis),{fe.u,fe.v});
// Now evaluate the solution field
if (!integrand.evalSol(solPt,fe,X4,ipa,elem_size,nb))
@ -1118,22 +1118,22 @@ void ASMs2Dmx::getBoundaryNodes (int lIndex, IntVec& nodes, int basis,
void ASMs2Dmx::swapProjectionBasis ()
{
if (projB2) {
ASMmxBase::geoBasis = ASMmxBase::geoBasis == 1 ? 2 : 1;
ASMmxBase::itgBasis = ASMmxBase::itgBasis == 1 ? 2 : 1;
std::swap(projB, projB2);
surf = this->getBasis(ASMmxBase::geoBasis);
surf = this->getBasis(ASMmxBase::itgBasis);
}
}
int ASMs2Dmx::getFirstItgElmNode () const
{
return std::accumulate(elem_size.begin(), elem_size.begin() + geoBasis-1, 0);
return std::accumulate(elem_size.begin(), elem_size.begin() + itgBasis-1, 0);
}
int ASMs2Dmx::getLastItgElmNode () const
{
return std::accumulate(elem_size.begin(), elem_size.begin() + geoBasis, -1);
return std::accumulate(elem_size.begin(), elem_size.begin() + itgBasis, -1);
}

30
src/ASM/ASMs2DmxLag.C
View File

@ -111,7 +111,7 @@ bool ASMs2DmxLag::generateFEMTopology ()
// Generate/check FE data for the geometry/basis1
bool haveFEdata = !MLGN.empty();
bool basis1IsOK = this->ASMs2DLag::generateFEMTopology();
geoBasis = 1;
itgBasis = 1;
if ((haveFEdata && !shareFE) || !basis1IsOK) return basis1IsOK;
// Order of 2nd basis in the two parametric directions (order = degree + 1)
@ -319,11 +319,11 @@ bool ASMs2DmxLag::integrate (Integrand& integrand,
}
// Compute Jacobian inverse of coordinate mapping and derivatives
if (!fe.Jacobian(Jac,Xnod,geoBasis,&bfs))
if (!fe.Jacobian(Jac,Xnod,itgBasis,&bfs))
continue; // skip singular points
// Cartesian coordinates of current integration point
X.assign(Xnod * fe.basis(geoBasis));
X.assign(Xnod * fe.basis(itgBasis));
// Evaluate the integrand and accumulate element contributions
fe.detJxW *= wg[0][i]*wg[1][j];
@ -370,8 +370,8 @@ bool ASMs2DmxLag::integrate (Integrand& integrand, int lIndex,
// Number of elements in each direction
const int nelx = (nxx[geoBasis-1]-1)/(elem_sizes[geoBasis-1][0]-1);
const int nely = (nyx[geoBasis-1]-1)/(elem_sizes[geoBasis-1][1]-1);
const int nelx = (nxx[itgBasis-1]-1)/(elem_sizes[itgBasis-1][0]-1);
const int nely = (nyx[itgBasis-1]-1)/(elem_sizes[itgBasis-1][1]-1);
std::map<char,size_t>::const_iterator iit = firstBp.find(lIndex%10);
size_t firstp = iit == firstBp.end() ? 0 : iit->second;
@ -428,18 +428,18 @@ bool ASMs2DmxLag::integrate (Integrand& integrand, int lIndex,
ok = false;
// Compute basis function derivatives and the edge normal
fe.detJxW = utl::Jacobian(Jac,normal,fe.grad(geoBasis),Xnod,
dNxdu[geoBasis-1],t1,t2);
fe.detJxW = utl::Jacobian(Jac,normal,fe.grad(itgBasis),Xnod,
dNxdu[itgBasis-1],t1,t2);
if (fe.detJxW == 0.0) continue; // skip singular points
for (size_t b = 0; b < nxx.size(); ++b)
if (b != (size_t)geoBasis-1)
if (b != (size_t)itgBasis-1)
fe.grad(b+1).multiply(dNxdu[b],Jac);
if (edgeDir < 0) normal *= -1.0;
// Cartesian coordinates of current integration point
X.assign(Xnod * fe.basis(geoBasis));
X.assign(Xnod * fe.basis(itgBasis));
// Evaluate the integrand and accumulate element contributions
fe.detJxW *= wg[i];
@ -510,11 +510,11 @@ bool ASMs2DmxLag::evalSolution (Matrix& sField, const IntegrandBase& integrand,
// Evaluate the secondary solution field at each point
for (size_t iel = 1; iel <= nel; iel++)
{
IntVec::const_iterator f2start = geoBasis == 1? MNPC[iel-1].begin() :
IntVec::const_iterator f2start = itgBasis == 1? MNPC[iel-1].begin() :
MNPC[iel-1].begin() +
std::accumulate(elem_size.begin()+geoBasis-2,
elem_size.begin()+geoBasis-1, 0);
IntVec::const_iterator f2end = f2start + elem_size[geoBasis-1];
std::accumulate(elem_size.begin()+itgBasis-2,
elem_size.begin()+itgBasis-1, 0);
IntVec::const_iterator f2end = f2start + elem_size[itgBasis-1];
IntVec mnpc1(f2start,f2end);
this->getElementCoordinates(Xnod,iel);
@ -533,11 +533,11 @@ bool ASMs2DmxLag::evalSolution (Matrix& sField, const IntegrandBase& integrand,
// Compute Jacobian inverse of the coordinate mapping and
// basis function derivatives w.r.t. Cartesian coordinates
if (!fe.Jacobian(Jac,Xnod,geoBasis,nullptr,&dNxdu))
if (!fe.Jacobian(Jac,Xnod,itgBasis,nullptr,&dNxdu))
continue; // skip singular points
// Now evaluate the solution field
if (!integrand.evalSol(solPt,fe,Xnod*fe.basis(geoBasis),
if (!integrand.evalSol(solPt,fe,Xnod*fe.basis(itgBasis),
MNPC[iel-1],elem_size,nb))
return false;
else if (sField.empty())

52
src/ASM/ASMs3Dmx.C
View File

@ -215,7 +215,7 @@ bool ASMs3Dmx::generateFEMTopology ()
else if (ASMmxBase::Type == ASMmxBase::SUBGRID)
projB = m_basis.front().get();
else // FULL_CONT_RAISE_BASISx
projB = m_basis[2-geoBasis].get();
projB = m_basis[2-itgBasis].get();
}
if (ASMmxBase::Type == ASMmxBase::SUBGRID)
@ -223,7 +223,7 @@ bool ASMs3Dmx::generateFEMTopology ()
delete svol;
}
geomB = svol = m_basis[geoBasis-1].get();
geomB = svol = m_basis[itgBasis-1].get();
nb.clear();
nb.reserve(m_basis.size());
@ -307,7 +307,7 @@ bool ASMs3Dmx::generateFEMTopology ()
myMNPC[iel][elm_node++] = last_node - (k*nv + j)*nu - i;
};
inod = std::accumulate(nb.begin(),nb.begin()+geoBasis-1,0u);
inod = std::accumulate(nb.begin(),nb.begin()+itgBasis-1,0u);
auto knotw = itg->basis(2).begin();
for (int k = 1; k <= itg->numCoefs(2); k++, ++knotw) {
auto knotv = itg->basis(1).begin();
@ -328,7 +328,7 @@ bool ASMs3Dmx::generateFEMTopology ()
size_t basis_ofs = 0;
int basis = 0;
for (const std::shared_ptr<Go::SplineVolume>& spline : m_basis) {
if (basis != geoBasis-1) {
if (basis != itgBasis-1) {
double ku = *knotu;
double kv = *knotv;
double kw = *knotw;
@ -404,7 +404,7 @@ bool ASMs3Dmx::getElementCoordinates (Matrix& X, int iel) const
size_t nenod = svol->order(0)*svol->order(1)*svol->order(2);
size_t lnod0 = 0;
for (int i = 1; i < geoBasis; ++i)
for (int i = 1; i < itgBasis; ++i)
lnod0 += m_basis[i-1]->order(0)*m_basis[i-1]->order(1)*m_basis[i-1]->order(2);
X.resize(3,nenod);
@ -608,11 +608,11 @@ bool ASMs3Dmx::integrate (Integrand& integrand,
// Compute Jacobian inverse of the coordinate mapping and
// basis function derivatives w.r.t. Cartesian coordinates
if (!fe.Jacobian(Jac,Xnod,geoBasis,&bfs))
if (!fe.Jacobian(Jac,Xnod,itgBasis,&bfs))
continue; // skip singular points
// Compute Hessian of coordinate mapping and 2nd order derivatives
if (use2ndDer && !fe.Hessian(Hess,Jac,Xnod,geoBasis,&bfs))
if (use2ndDer && !fe.Hessian(Hess,Jac,Xnod,itgBasis,&bfs))
ok = false;
// Compute G-matrix
@ -620,7 +620,7 @@ bool ASMs3Dmx::integrate (Integrand& integrand,
utl::getGmat(Jac,dXidu,fe.G);
// Cartesian coordinates of current integration point
X.assign(Xnod * fe.basis(geoBasis));
X.assign(Xnod * fe.basis(itgBasis));
// Evaluate the integrand and accumulate element contributions
fe.detJxW *= dV*wg[0][i]*wg[1][j]*wg[2][k];
@ -822,18 +822,18 @@ bool ASMs3Dmx::integrate (Integrand& integrand, int lIndex,
// Compute Jacobian inverse of the coordinate mapping and
// basis function derivatives w.r.t. Cartesian coordinates
fe.detJxW = utl::Jacobian(Jac,normal,fe.grad(geoBasis),Xnod,
dNxdu[geoBasis-1],t1,t2);
fe.detJxW = utl::Jacobian(Jac,normal,fe.grad(itgBasis),Xnod,
dNxdu[itgBasis-1],t1,t2);
if (fe.detJxW == 0.0) continue; // skip singular points
for (size_t b = 0; b < m_basis.size(); ++b)
if (b != (size_t)geoBasis-1)
if (b != (size_t)itgBasis-1)
fe.grad(b+1).multiply(dNxdu[b],Jac);
if (faceDir < 0) normal *= -1.0;
// Cartesian coordinates of current integration point
X.assign(Xnod * fe.basis(geoBasis));
X.assign(Xnod * fe.basis(itgBasis));
// Evaluate the integrand and accumulate element contributions
fe.detJxW *= dA*wg[i]*wg[j];
@ -1039,13 +1039,13 @@ bool ASMs3Dmx::integrate (Integrand& integrand,
}
// Compute basis function derivatives and the edge normal
fe.detJxW = utl::Jacobian(Jac,normal,fe.grad(geoBasis+m_basis.size()),
Xnod,dNxdu[geoBasis-1+m_basis.size()],t1,t2);
fe.detJxW = utl::Jacobian(Jac,normal,fe.grad(geoBasis),Xnod,
dNxdu[geoBasis-1],t1,t2);
fe.detJxW = utl::Jacobian(Jac,normal,fe.grad(itgBasis+m_basis.size()),
Xnod,dNxdu[itgBasis-1+m_basis.size()],t1,t2);
fe.detJxW = utl::Jacobian(Jac,normal,fe.grad(itgBasis),Xnod,
dNxdu[itgBasis-1],t1,t2);
if (fe.detJxW == 0.0) continue; // skip singular points
for (size_t b = 0; b < m_basis.size(); ++b)
if (b != (size_t)geoBasis-1) {
if (b != (size_t)itgBasis-1) {
fe.grad(b+1).multiply(dNxdu[b],Jac);
fe.grad(b+1+m_basis.size()).multiply(dNxdu[b+m_basis.size()],Jac);
}
@ -1053,7 +1053,7 @@ bool ASMs3Dmx::integrate (Integrand& integrand,
if (faceDir < 0) normal *= -1.0;
// Cartesian coordinates of current integration point
X.assign(Xnod * fe.basis(geoBasis));
X.assign(Xnod * fe.basis(itgBasis));
// Evaluate the integrand and accumulate element contributions
fe.detJxW *= 0.25*dA*wg[i]*wg[j];
@ -1206,8 +1206,8 @@ bool ASMs3Dmx::evalSolution (Matrix& sField, const IntegrandBase& integrand,
splinex[b][i].left_idx,ip[b]);
// Fetch associated control point coordinates
if (b == (size_t)geoBasis-1)
utl::gather(ip[geoBasis-1], 3, Xnod, Xtmp);
if (b == (size_t)itgBasis-1)
utl::gather(ip[itgBasis-1], 3, Xnod, Xtmp);
for (int& c : ip[b]) c += ofs;
ipa.insert(ipa.end(), ip[b].begin(), ip[b].end());
@ -1224,11 +1224,11 @@ bool ASMs3Dmx::evalSolution (Matrix& sField, const IntegrandBase& integrand,
// Compute Jacobian inverse of the coordinate mapping and
// basis function derivatives w.r.t. Cartesian coordinate
if (!fe.Jacobian(Jac,Xtmp,geoBasis,nullptr,&dNxdu))
if (!fe.Jacobian(Jac,Xtmp,itgBasis,nullptr,&dNxdu))
continue; // skip singular points
// Cartesian coordinates of current integration point
utl::Point X4(Xtmp * fe.basis(geoBasis),{fe.u,fe.v,fe.w});
utl::Point X4(Xtmp * fe.basis(itgBasis),{fe.u,fe.v,fe.w});
// Now evaluate the solution field
if (!integrand.evalSol(solPt,fe,X4,ipa,elem_size,nb))
@ -1271,22 +1271,22 @@ void ASMs3Dmx::getBoundaryNodes (int lIndex, IntVec& nodes, int basis,
void ASMs3Dmx::swapProjectionBasis ()
{
if (projB2) {
ASMmxBase::geoBasis = ASMmxBase::geoBasis == 1 ? 2 : 1;
ASMmxBase::itgBasis = ASMmxBase::itgBasis == 1 ? 2 : 1;
std::swap(projB, projB2);
svol = this->getBasis(ASMmxBase::geoBasis);
svol = this->getBasis(ASMmxBase::itgBasis);
}
}
int ASMs3Dmx::getFirstItgElmNode () const
{
return std::accumulate(elem_size.begin(), elem_size.begin() + geoBasis-1, 0);
return std::accumulate(elem_size.begin(), elem_size.begin() + itgBasis-1, 0);
}
int ASMs3Dmx::getLastItgElmNode () const
{
return std::accumulate(elem_size.begin(), elem_size.begin() + geoBasis, -1);
return std::accumulate(elem_size.begin(), elem_size.begin() + itgBasis, -1);
}

30
src/ASM/ASMs3DmxLag.C
View File

@ -112,7 +112,7 @@ bool ASMs3DmxLag::generateFEMTopology ()
// Generate/check FE data for the geometry/basis1
bool haveFEdata = !MLGN.empty();
bool basis1IsOK = this->ASMs3DLag::generateFEMTopology();
geoBasis = 1;
itgBasis = 1;
if ((haveFEdata && !shareFE) || !basis1IsOK) return basis1IsOK;
// Order of 2nd basis in the three parametric directions (order = degree + 1)
@ -355,11 +355,11 @@ bool ASMs3DmxLag::integrate (Integrand& integrand,
}
// Compute Jacobian inverse of coordinate mapping and derivatives
if (!fe.Jacobian(Jac,Xnod,geoBasis,&bfs))
if (!fe.Jacobian(Jac,Xnod,itgBasis,&bfs))
continue; // skip singular points
// Cartesian coordinates of current integration point
X.assign(Xnod * fe.basis(geoBasis));
X.assign(Xnod * fe.basis(itgBasis));
// Evaluate the integrand and accumulate element contributions
fe.detJxW *= wg[0][i]*wg[1][j]*wg[2][k];
@ -415,8 +415,8 @@ bool ASMs3DmxLag::integrate (Integrand& integrand, int lIndex,
integrand.setNeumannOrder(1 + lIndex/10);
// Number of elements in each direction
const int nel1 = (nxx[geoBasis-1]-1)/(elem_sizes[geoBasis-1][0]-1);
const int nel2 = (nyx[geoBasis-1]-1)/(elem_sizes[geoBasis-1][1]-1);
const int nel1 = (nxx[itgBasis-1]-1)/(elem_sizes[itgBasis-1][0]-1);
const int nel2 = (nyx[itgBasis-1]-1)/(elem_sizes[itgBasis-1][1]-1);
std::map<char,size_t>::const_iterator iit = firstBp.find(lIndex%10);
size_t firstp = iit == firstBp.end() ? 0 : iit->second;
@ -494,18 +494,18 @@ bool ASMs3DmxLag::integrate (Integrand& integrand, int lIndex,
ok = false;
// Compute basis function derivatives and the edge normal
fe.detJxW = utl::Jacobian(Jac,normal,fe.grad(geoBasis),Xnod,
dNxdu[geoBasis-1],t1,t2);
fe.detJxW = utl::Jacobian(Jac,normal,fe.grad(itgBasis),Xnod,
dNxdu[itgBasis-1],t1,t2);
if (fe.detJxW == 0.0) continue; // skip singular points
for (size_t b = 0; b < nxx.size(); ++b)
if (b != (size_t)geoBasis-1)
if (b != (size_t)itgBasis-1)
fe.grad(b+1).multiply(dNxdu[b],Jac);
if (faceDir < 0) normal *= -1.0;
// Cartesian coordinates of current integration point
X.assign(Xnod * fe.basis(geoBasis));
X.assign(Xnod * fe.basis(itgBasis));
// Evaluate the integrand and accumulate element contributions
fe.detJxW *= wg[i]*wg[j];
@ -577,11 +577,11 @@ bool ASMs3DmxLag::evalSolution (Matrix& sField, const IntegrandBase& integrand,
// Evaluate the secondary solution field at each point
for (size_t iel = 1; iel <= nel; iel++)
{
IntVec::const_iterator f2start = geoBasis == 1? MNPC[iel-1].begin() :
IntVec::const_iterator f2start = itgBasis == 1? MNPC[iel-1].begin() :
MNPC[iel-1].begin() +
std::accumulate(elem_size.begin()+geoBasis-2,
elem_size.begin()+geoBasis-1, 0);
IntVec::const_iterator f2end = f2start + elem_size[geoBasis-1];
std::accumulate(elem_size.begin()+itgBasis-2,
elem_size.begin()+itgBasis-1, 0);
IntVec::const_iterator f2end = f2start + elem_size[itgBasis-1];
IntVec mnpc1(f2start,f2end);
this->getElementCoordinates(Xnod,iel);
@ -603,11 +603,11 @@ bool ASMs3DmxLag::evalSolution (Matrix& sField, const IntegrandBase& integrand,
// Compute Jacobian inverse of the coordinate mapping and
// basis function derivatives w.r.t. Cartesian coordinates
if (!fe.Jacobian(Jac,Xnod,geoBasis,nullptr,&dNxdu))
if (!fe.Jacobian(Jac,Xnod,itgBasis,nullptr,&dNxdu))
continue; // skip singular points
// Now evaluate the solution field
if (!integrand.evalSol(solPt,fe,Xnod*fe.basis(geoBasis),
if (!integrand.evalSol(solPt,fe,Xnod*fe.basis(itgBasis),
MNPC[iel-1],elem_size,nb))
return false;
else if (sField.empty())

80
src/ASM/LR/ASMu2Dmx.C
View File

@ -220,7 +220,7 @@ bool ASMu2Dmx::generateFEMTopology ()
else if (ASMmxBase::Type == ASMmxBase::SUBGRID)
projB = m_basis.front();
else // FULL_CONT_RAISE_BASISx
projB = m_basis[2-ASMmxBase::geoBasis];
projB = m_basis[2-ASMmxBase::itgBasis];
}
if (ASMmxBase::Type == ASMmxBase::SUBGRID)
@ -234,7 +234,7 @@ bool ASMu2Dmx::generateFEMTopology ()
}
projB->generateIDs();
refB->generateIDs();
lrspline = m_basis[geoBasis-1];
lrspline = m_basis[itgBasis-1];
nb.clear();
nb.reserve(m_basis.size());
@ -249,7 +249,7 @@ bool ASMu2Dmx::generateFEMTopology ()
if (shareFE == 'F') return true;
nel = m_basis[geoBasis-1]->nElements();
nel = m_basis[itgBasis-1]->nElements();
nnod = std::accumulate(nb.begin(), nb.end(), 0);
myMLGE.resize(nel,0);
@ -259,7 +259,7 @@ bool ASMu2Dmx::generateFEMTopology ()
it->generateIDs();
size_t iel = 0;
for (const LR::Element* el1 : m_basis[geoBasis-1]->getAllElements())
for (const LR::Element* el1 : m_basis[itgBasis-1]->getAllElements())
{
size_t nfunc = 0;
for (const SplinePtr& it : m_basis)
@ -284,7 +284,7 @@ bool ASMu2Dmx::generateFEMTopology ()
std::cout <<"NEL = "<< nel <<" NNOD = "<< nnod << std::endl;
#endif
geomB = m_basis[geoBasis-1];
geomB = m_basis[itgBasis-1];
this->generateBezierBasis();
this->generateBezierExtraction();
@ -355,7 +355,7 @@ bool ASMu2Dmx::integrate (Integrand& integrand,
double param[3] = { 0.0, 0.0, 0.0 };
Vec4 X(param,time.t);
int geoEl = els[geoBasis-1];
int geoEl = els[itgBasis-1];
fe.iel = MLGE[geoEl-1];
// Get element area in the parameter space
@ -417,11 +417,11 @@ bool ASMu2Dmx::integrate (Integrand& integrand,
// Compute Jacobian inverse of the coordinate mapping and
// basis function derivatives w.r.t. Cartesian coordinates
if (!fe.Jacobian(Jac,Xnod,geoBasis,&bfs))
if (!fe.Jacobian(Jac,Xnod,itgBasis,&bfs))
continue; // skip singular points
// Compute Hessian of coordinate mapping and 2nd order derivatives
if (use2ndDer && !fe.Hessian(Hess,Jac,Xnod,geoBasis,&bfs))
if (use2ndDer && !fe.Hessian(Hess,Jac,Xnod,itgBasis,&bfs))
ok = false;
// Compute G-matrix
@ -429,7 +429,7 @@ bool ASMu2Dmx::integrate (Integrand& integrand,
utl::getGmat(Jac,dXidu,fe.G);
// Cartesian coordinates of current integration point
X.assign(Xnod * fe.basis(geoBasis));
X.assign(Xnod * fe.basis(itgBasis));
// Evaluate the integrand and accumulate element contributions
fe.detJxW *= dA*wg[0][i]*wg[1][j];
@ -504,7 +504,7 @@ bool ASMu2Dmx::integrate (Integrand& integrand, int lIndex,
// === Assembly loop over all elements on the patch edge =====================
int iel = 0;
for (const LR::Element* el1 : m_basis[geoBasis-1]->getAllElements())
for (const LR::Element* el1 : m_basis[itgBasis-1]->getAllElements())
{
++iel;
@ -512,10 +512,10 @@ bool ASMu2Dmx::integrate (Integrand& integrand, int lIndex,
bool skipMe = false;
switch (edgeDir)
{
case -1: if (el1->umin() != m_basis[geoBasis-1]->startparam(0)) skipMe = true; break;
case 1: if (el1->umax() != m_basis[geoBasis-1]->endparam(0) ) skipMe = true; break;
case -2: if (el1->vmin() != m_basis[geoBasis-1]->startparam(1)) skipMe = true; break;
case 2: if (el1->vmax() != m_basis[geoBasis-1]->endparam(1) ) skipMe = true; break;
case -1: if (el1->umin() != m_basis[itgBasis-1]->startparam(0)) skipMe = true; break;
case 1: if (el1->umax() != m_basis[itgBasis-1]->endparam(0) ) skipMe = true; break;
case -2: if (el1->vmin() != m_basis[itgBasis-1]->startparam(1)) skipMe = true; break;
case 2: if (el1->vmax() != m_basis[itgBasis-1]->endparam(1) ) skipMe = true; break;
}
if (skipMe) continue;
@ -528,7 +528,7 @@ bool ASMu2Dmx::integrate (Integrand& integrand, int lIndex,
this->getElementsAt(el1->midpoint(),els,elem_sizes);
MxFiniteElement fe(elem_sizes,firstp);
int geoEl = els[geoBasis-1];
int geoEl = els[itgBasis-1];
fe.iel = MLGE[geoEl-1];
fe.xi = fe.eta = edgeDir < 0 ? -1.0 : 1.0;
firstp += nGP; // Global integration point counter
@ -574,19 +574,19 @@ bool ASMu2Dmx::integrate (Integrand& integrand, int lIndex,
// Compute Jacobian inverse of the coordinate mapping and
// basis function derivatives w.r.t. Cartesian coordinates
fe.detJxW = utl::Jacobian(Jac,normal,fe.grad(geoBasis),Xnod,
dNxdu[geoBasis-1],t1,t2);
fe.detJxW = utl::Jacobian(Jac,normal,fe.grad(itgBasis),Xnod,
dNxdu[itgBasis-1],t1,t2);
if (fe.detJxW == 0.0) continue; // skip singular points
for (size_t b = 1; b <= m_basis.size(); ++b)
if ((int)b != geoBasis)
if ((int)b != itgBasis)
fe.grad(b).multiply(dNxdu[b-1],Jac);
if (edgeDir < 0)
normal *= -1.0;
// Cartesian coordinates of current integration point
X.assign(Xnod * fe.basis(geoBasis));
X.assign(Xnod * fe.basis(itgBasis));
// Evaluate the integrand and accumulate element contributions
fe.detJxW *= dS*wg[i];
@ -654,11 +654,11 @@ bool ASMu2Dmx::integrate (Integrand& integrand,
elem_sizes.front() = el1->nBasisFunctions();
// Set up control point coordinates for current element
if (!this->getElementCoordinates(Xnod,els[geoBasis-1]))
if (!this->getElementCoordinates(Xnod,els[itgBasis-1]))
return false;
LocalIntegral* A = integrand.getLocalIntegral(elem_sizes, iel);
bool ok = integrand.initElement(MNPC[els[geoBasis-1]-1],elem_sizes,nb,*A);
bool ok = integrand.initElement(MNPC[els[itgBasis-1]-1],elem_sizes,nb,*A);
size_t origSize = A->vec.size();
int bit = 8;
@ -700,7 +700,7 @@ bool ASMu2Dmx::integrate (Integrand& integrand,
elem_sizes2.front() = el2->nBasisFunctions();
LocalIntegral* A_neigh = integrand.getLocalIntegral(elem_sizes2, el_neigh);
ok = integrand.initElement(MNPC[els2[geoBasis-1]-1],
ok = integrand.initElement(MNPC[els2[itgBasis-1]-1],
elem_sizes2,nb,*A_neigh);
// Element sizes for both elements
@ -709,7 +709,7 @@ bool ASMu2Dmx::integrate (Integrand& integrand,
std::back_inserter(elem_sizes3));
MxFiniteElement fe(elem_sizes3);
fe.h = this->getElementCorners(els2[geoBasis-1], fe.XC);
fe.h = this->getElementCorners(els2[itgBasis-1], fe.XC);
if (!A_neigh->vec.empty()) {
A->vec.resize(origSize+A_neigh->vec.size());
@ -735,7 +735,7 @@ bool ASMu2Dmx::integrate (Integrand& integrand,
gpar[1].fill(v1);
}
Matrix Xnod2, Jac2;
ok &= this->getElementCoordinates(Xnod2,els2[geoBasis-1]);
ok &= this->getElementCoordinates(Xnod2,els2[itgBasis-1]);
for (int g = 0; g < nGP && ok; g++, ++fe.iGP)
{
@ -760,15 +760,15 @@ bool ASMu2Dmx::integrate (Integrand& integrand,
// Compute Jacobian inverse of the coordinate mapping and
// basis function derivatives w.r.t. Cartesian coordinates
fe.detJxW = utl::Jacobian(Jac2,normal,
fe.grad(geoBasis+m_basis.size()),Xnod2,
dNxdu[geoBasis-1+m_basis.size()],t1,t2);
fe.grad(itgBasis+m_basis.size()),Xnod2,
dNxdu[itgBasis-1+m_basis.size()],t1,t2);
fe.detJxW = utl::Jacobian(Jac,normal,
fe.grad(geoBasis),Xnod,
dNxdu[geoBasis-1],t1,t2);
fe.grad(itgBasis),Xnod,
dNxdu[itgBasis-1],t1,t2);
if (fe.detJxW == 0.0) continue; // skip singular points
for (size_t b = 1; b <= m_basis.size(); ++b)
if ((int)b != geoBasis) {
if ((int)b != itgBasis) {
fe.grad(b).multiply(dNxdu[b-1],Jac);
fe.grad(b+m_basis.size()).multiply(dNxdu[b-1+m_basis.size()],Jac);
}
@ -777,7 +777,7 @@ bool ASMu2Dmx::integrate (Integrand& integrand,
normal *= -1.0;
// Cartesian coordinates of current integration point
X.assign(Xnod * fe.basis(geoBasis));
X.assign(Xnod * fe.basis(itgBasis));
// Evaluate the integrand and accumulate element contributions
fe.detJxW *= 0.5*dS*wg[g];
@ -797,7 +797,7 @@ bool ASMu2Dmx::integrate (Integrand& integrand,
ok = false;
// Assembly of global system integral
if (ok && !glInt.assemble(A,els[geoBasis-1]))
if (ok && !glInt.assemble(A,els[itgBasis-1]))
ok = false;
A->destruct();
@ -911,26 +911,26 @@ bool ASMu2Dmx::evalSolution (Matrix& sField, const IntegrandBase& integrand,
}
// Set up control point (nodal) coordinates for current element
if (!this->getElementCoordinates(Xnod,els[geoBasis-1])) return false;
if (!this->getElementCoordinates(Xnod,els[itgBasis-1])) return false;
// Compute Jacobian inverse of the coordinate mapping and
// basis function derivatives w.r.t. Cartesian coordinates
if (!fe.Jacobian(Jac,Xnod,geoBasis,nullptr,&dNxdu))
if (!fe.Jacobian(Jac,Xnod,itgBasis,nullptr,&dNxdu))
continue; // skip singular points
// Compute Hessian of coordinate mapping and 2nd order derivatives
if (use2ndDer && !fe.Hessian(Hess,Jac,Xnod,geoBasis,nullptr,&d2Nxdu2))
if (use2ndDer && !fe.Hessian(Hess,Jac,Xnod,itgBasis,nullptr,&d2Nxdu2))
return false;
// Cartesian coordinates of current integration point
fe.u = gpar[0][i];
fe.v = gpar[1][i];
utl::Point X4(Xnod*fe.basis(geoBasis),{fe.u,fe.v});
utl::Point X4(Xnod*fe.basis(itgBasis),{fe.u,fe.v});
// Now evaluate the solution field
Vector solPt;
if (!integrand.evalSol(solPt,fe,X4,
MNPC[els[geoBasis-1]-1],elem_sizes,nb))
MNPC[els[itgBasis-1]-1],elem_sizes,nb))
return false;
else if (sField.empty())
sField.resize(solPt.size(),nPoints,true);
@ -1114,7 +1114,7 @@ void ASMu2Dmx::getBoundaryNodes (int lIndex, IntVec& nodes, int basis,
void ASMu2Dmx::remapErrors (RealArray& errors,
const RealArray& origErr, bool elemErrors) const
{
const LR::LRSplineSurface* geo = this->getBasis(ASMmxBase::geoBasis);
const LR::LRSplineSurface* geo = this->getBasis(ASMmxBase::itgBasis);
for (const LR::Element* elm : geo->getAllElements()) {
int rEl = refB->getElementContaining(elm->midpoint());
if (elemErrors)
@ -1198,10 +1198,10 @@ void ASMu2Dmx::copyRefinement (LR::LRSplineSurface* basis,
void ASMu2Dmx::swapProjectionBasis ()
{
if (projB2) {
ASMmxBase::geoBasis = ASMmxBase::geoBasis == 1 ? 2 : 1;
ASMmxBase::itgBasis = ASMmxBase::itgBasis == 1 ? 2 : 1;
std::swap(projB, projB2);
std::swap(projThreadGroups, proj2ThreadGroups);
lrspline = m_basis[ASMmxBase::geoBasis-1];
lrspline = m_basis[ASMmxBase::itgBasis-1];
geomB = lrspline;
this->generateBezierBasis();
this->generateBezierExtraction();
@ -1236,7 +1236,7 @@ BasisFunctionVals ASMu2Dmx::BasisFunctionCache::calculatePt (size_t el,
const ASMu2Dmx& pch = static_cast<const ASMu2Dmx&>(patch);
const LR::Element* el1 = pch.getBasis(ASMmxBase::geoBasis)->getElement(el);
const LR::Element* el1 = pch.getBasis(ASMmxBase::itgBasis)->getElement(el);
size_t el_b = patch.getBasis(basis)->getElementContaining(el1->midpoint());
BasisFunctionVals result;

60
src/ASM/LR/ASMu3Dmx.C
View File

@ -42,7 +42,7 @@ ASMu3Dmx::ASMu3Dmx (const CharVec& n_f)
bezierExtractmx(myBezierExtractmx)
{
threadBasis = nullptr;
myGeoBasis = ASMmxBase::geoBasis;
myGeoBasis = ASMmxBase::itgBasis;
}
@ -54,7 +54,7 @@ ASMu3Dmx::ASMu3Dmx (const ASMu3Dmx& patch, const CharVec& n_f)
threadBasis = patch.threadBasis;
nfx = patch.nfx;
nb = patch.nb;
myGeoBasis = ASMmxBase::geoBasis;
myGeoBasis = ASMmxBase::itgBasis;
}
@ -198,11 +198,11 @@ bool ASMu3Dmx::generateFEMTopology ()
m_basis.push_back(std::make_shared<LR::LRSplineVolume>(vvec[b].get()));
// make a backup as establishBases resets it
int geoB = ASMmxBase::geoBasis;
int geoB = ASMmxBase::itgBasis;
std::shared_ptr<Go::SplineVolume> otherBasis =
ASMmxBase::establishBases(tensorspline,
ASMmxBase::FULL_CONT_RAISE_BASIS1).front();
geoBasis = geoB;
itgBasis = geoB;
// we need to project on something that is not one of our bases
if (!projB) {
@ -213,7 +213,7 @@ bool ASMu3Dmx::generateFEMTopology ()
else if (ASMmxBase::Type == ASMmxBase::SUBGRID)
projB = m_basis.front();
else // FULL_CONT_RAISE_BASISx
projB = m_basis[2-ASMmxBase::geoBasis];
projB = m_basis[2-ASMmxBase::itgBasis];
}
if (ASMmxBase::Type == ASMmxBase::SUBGRID)
@ -224,14 +224,14 @@ bool ASMu3Dmx::generateFEMTopology ()
delete tensorspline;
tensorspline = nullptr;
}
lrspline = m_basis[geoBasis-1];
lrspline = m_basis[itgBasis-1];
projB->generateIDs();
projB->getElementContaining(projB->getElement(0)->midpoint()); // to force cache generation
if (projB2) {
projB2->generateIDs();
projB2->getElementContaining(projB2->getElement(0)->midpoint()); // to force cache generation
}
myGeoBasis = ASMmxBase::geoBasis;
myGeoBasis = ASMmxBase::itgBasis;
nb.clear();
nb.reserve(m_basis.size());
@ -247,7 +247,7 @@ bool ASMu3Dmx::generateFEMTopology ()
if (shareFE == 'F') return true;
nel = m_basis[geoBasis-1]->nElements();
nel = m_basis[itgBasis-1]->nElements();
nnod = std::accumulate(nb.begin(), nb.end(), 0);
myMLGE.resize(nel,0);
@ -257,7 +257,7 @@ bool ASMu3Dmx::generateFEMTopology ()
it->generateIDs();
size_t iel = 0;
for (const LR::Element* el1 : m_basis[geoBasis-1]->getAllElements())
for (const LR::Element* el1 : m_basis[itgBasis-1]->getAllElements())
{
size_t nfunc = 0;
for (const SplinePtr& it : m_basis)
@ -299,7 +299,7 @@ bool ASMu3Dmx::generateFEMTopology ()
std::cout <<"NEL = "<< nel <<" NNOD = "<< nnod << std::endl;
#endif
geomB = m_basis[geoBasis-1];
geomB = m_basis[itgBasis-1];
return true;
}
@ -403,11 +403,11 @@ bool ASMu3Dmx::integrate (Integrand& integrand,
for (int i = 0; i < ng[0]; i++, jp++)
{
// Fetch basis function derivatives at current integration point
const BasisFunctionVals& bfs = myCache[geoBasis-1]->getVals(iEl,jp);
const BasisFunctionVals& bfs = myCache[itgBasis-1]->getVals(iEl,jp);
// Compute Jacobian determinant of coordinate mapping
// and multiply by weight of current integration point
double detJac = utl::Jacobian(Jac,fe.grad(geoBasis),
double detJac = utl::Jacobian(Jac,fe.grad(itgBasis),
Xnod,bfs.dNdu,false);
double weight = dV*wg[0][i]*wg[1][j]*wg[2][k];
@ -427,7 +427,7 @@ bool ASMu3Dmx::integrate (Integrand& integrand,
double u0 = 0.5*(el->getParmin(0) + el->getParmax(0));
double v0 = 0.5*(el->getParmin(1) + el->getParmax(1));
double w0 = 0.5*(el->getParmin(2) + el->getParmax(2));
this->getBasis(geoBasis)->point(X0,u0,v0,w0);
this->getBasis(itgBasis)->point(X0,u0,v0,w0);
X.assign(SplineUtils::toVec3(X0));
}
@ -468,11 +468,11 @@ bool ASMu3Dmx::integrate (Integrand& integrand,
// Compute Jacobian inverse of the coordinate mapping and
// basis function derivatives w.r.t. Cartesian coordinates
if (!fe.Jacobian(Jac,Xnod,geoBasis,&bfs))
if (!fe.Jacobian(Jac,Xnod,itgBasis,&bfs))
continue; // skip singular points
// Compute Hessian of coordinate mapping and 2nd order derivatives
if (use2ndDer && !fe.Hessian(Hess,Jac,Xnod,geoBasis,&bfs))
if (use2ndDer && !fe.Hessian(Hess,Jac,Xnod,itgBasis,&bfs))
ok = false;
// Compute G-matrix
@ -480,7 +480,7 @@ bool ASMu3Dmx::integrate (Integrand& integrand,
utl::getGmat(Jac,dXidu,fe.G);
// Cartesian coordinates of current integration point
X.assign(Xnod * fe.basis(geoBasis));
X.assign(Xnod * fe.basis(itgBasis));
// Evaluate the integrand and accumulate element contributions
fe.detJxW *= dV*wg[0][i]*wg[1][j]*wg[2][k];
@ -545,7 +545,7 @@ bool ASMu3Dmx::integrate (Integrand& integrand, int lIndex,
// Fetch all elements on the chosen face
std::vector<LR::Element*> edgeElms;
this->getBasis(geoBasis)->getEdgeElements(edgeElms,getFaceEnum(lIndex));
this->getBasis(itgBasis)->getEdgeElements(edgeElms,getFaceEnum(lIndex));
// === Assembly loop over all elements on the patch face =====================
@ -570,12 +570,12 @@ bool ASMu3Dmx::integrate (Integrand& integrand, int lIndex,
if (-1-d == faceDir)
{
gpar[d].resize(1);
gpar[d].fill(this->getBasis(geoBasis)->startparam(d));
gpar[d].fill(this->getBasis(itgBasis)->startparam(d));
}
else if (1+d == faceDir)
{
gpar[d].resize(1);
gpar[d].fill(this->getBasis(geoBasis)->endparam(d));
gpar[d].fill(this->getBasis(itgBasis)->endparam(d));
}
else
this->getGaussPointParameters(gpar[d],d,nGP,iEl+1,xg);
@ -653,12 +653,12 @@ bool ASMu3Dmx::integrate (Integrand& integrand, int lIndex,
this->evaluateBasis(iEl, fe, dNxdu[b-1], b);
// Compute basis function derivatives and the face normal
fe.detJxW = utl::Jacobian(Jac, normal, fe.grad(geoBasis), Xnod,
dNxdu[geoBasis-1], t1, t2);
fe.detJxW = utl::Jacobian(Jac, normal, fe.grad(itgBasis), Xnod,
dNxdu[itgBasis-1], t1, t2);
if (fe.detJxW == 0.0) continue; // skip singular points
for (size_t b = 1; b <= m_basis.size(); ++b)
if ((int)b != geoBasis)
if ((int)b != itgBasis)
fe.grad(b).multiply(dNxdu[b-1],Jac);
if (faceDir < 0) normal *= -1.0;
@ -668,7 +668,7 @@ bool ASMu3Dmx::integrate (Integrand& integrand, int lIndex,
utl::getGmat(Jac,dXidu,fe.G);
// Cartesian coordinates of current integration point
X.assign(Xnod * fe.basis(geoBasis));
X.assign(Xnod * fe.basis(itgBasis));
// Evaluate the integrand and accumulate element contributions
fe.detJxW *= dA*wg[i]*wg[j];
@ -794,27 +794,27 @@ bool ASMu3Dmx::evalSolution (Matrix& sField, const IntegrandBase& integrand,
}
// Set up control point (nodal) coordinates for current element
if (!this->getElementCoordinates(Xnod,els[geoBasis-1])) return false;
if (!this->getElementCoordinates(Xnod,els[itgBasis-1])) return false;
// Compute Jacobian inverse of the coordinate mapping and
// basis function derivatives w.r.t. Cartesian coordinates
if (!fe.Jacobian(Jac,Xnod,geoBasis,nullptr,&dNxdu))
if (!fe.Jacobian(Jac,Xnod,itgBasis,nullptr,&dNxdu))
continue; // skip singular points
// Compute Hessian of coordinate mapping and 2nd order derivatives
if (use2ndDer && !fe.Hessian(Hess,Jac,Xnod,geoBasis,nullptr,&d2Nxdu2))
if (use2ndDer && !fe.Hessian(Hess,Jac,Xnod,itgBasis,nullptr,&d2Nxdu2))
return false;
// Cartesian coordinates of current integration point
fe.u = gpar[0][i];
fe.v = gpar[1][i];
fe.w = gpar[2][i];
utl::Point X4(Xnod*fe.basis(geoBasis), {fe.u, fe.v, fe.w});
utl::Point X4(Xnod*fe.basis(itgBasis), {fe.u, fe.v, fe.w});
// Now evaluate the solution field
Vector solPt;
if (!integrand.evalSol(solPt,fe,X4,
MNPC[els[geoBasis-1]-1],elem_sizes,nb))
MNPC[els[itgBasis-1]-1],elem_sizes,nb))
return false;
else if (sField.empty())
sField.resize(solPt.size(),nPoints,true);
@ -971,7 +971,7 @@ void ASMu3Dmx::generateThreadGroups (const Integrand& integrand, bool silence,
void ASMu3Dmx::remapErrors (RealArray& errors,
const RealArray& origErr, bool elemErrors) const
{
const LR::LRSplineVolume* geo = this->getBasis(ASMmxBase::geoBasis);
const LR::LRSplineVolume* geo = this->getBasis(ASMmxBase::itgBasis);
for (const LR::Element* elm : geo->getAllElements()) {
int rEl = refB->getElementContaining(elm->midpoint());
if (elemErrors)
@ -1013,7 +1013,7 @@ void ASMu3Dmx::copyRefinement (LR::LRSplineVolume* basis,
void ASMu3Dmx::swapProjectionBasis ()
{
if (projB2) {
ASMmxBase::geoBasis = ASMmxBase::geoBasis == 1 ? 2 : 1;
ASMmxBase::itgBasis = ASMmxBase::itgBasis == 1 ? 2 : 1;
std::swap(projB, projB2);
std::swap(projThreadGroups, proj2ThreadGroups);
}

4
src/ASM/LR/LRSplineFields2Dmx.C
View File

@ -103,7 +103,7 @@ bool LRSplineFields2Dmx::gradFE (const ItgPoint& x, Matrix& grad) const
// Evaluate the basis functions at the given point
Matrix Xnod, Jac, dNdX;
const LR::Element* elm;
const LR::LRSplineSurface* geo = surf->getBasis(ASMmxBase::geoBasis);
const LR::LRSplineSurface* geo = surf->getBasis(ASMmxBase::itgBasis);
if (!geo)
geo = surf->getBasis(1);
if (!LRSplineField::evalMapping(*geo,x,elm,Xnod,Jac,dNdX,surf->rational()))
@ -142,7 +142,7 @@ bool LRSplineFields2Dmx::hessianFE (const ItgPoint& x, Matrix3D& H) const
Matrix Xnod, Jac, dNdX;
Matrix3D d2NdX2, Hess;
const LR::Element* elm;
const LR::LRSplineSurface* geo = surf->getBasis(ASMmxBase::geoBasis);
const LR::LRSplineSurface* geo = surf->getBasis(ASMmxBase::itgBasis);
if (!geo)
geo = surf->getBasis(1);
if (!LRSplineField::evalMapping(*geo,x,elm,Xnod,Jac,

4
src/ASM/LR/LRSplineFields3Dmx.C
View File

@ -101,7 +101,7 @@ bool LRSplineFields3Dmx::gradFE (const ItgPoint& x, Matrix& grad) const
// Evaluate the basis functions at the given point
Matrix Xnod, Jac, dNdX;
const LR::LRSplineVolume* geo = vol->getBasis(ASMmxBase::geoBasis);
const LR::LRSplineVolume* geo = vol->getBasis(ASMmxBase::itgBasis);
if (!geo)
geo = vol->getBasis(1);
const LR::Element* elm;
@ -141,7 +141,7 @@ bool LRSplineFields3Dmx::hessianFE (const ItgPoint& x, Matrix3D& H) const
// Evaluate the basis functions at the given point
Matrix Xnod, Jac, dNdX;
Matrix3D d2NdX2, Hess;
const LR::LRSplineVolume* geo = vol->getBasis(ASMmxBase::geoBasis);
const LR::LRSplineVolume* geo = vol->getBasis(ASMmxBase::itgBasis);
if (!geo)
geo = vol->getBasis(1);
const LR::Element* elm;

6
src/ASM/SplineFields2Dmx.C
View File

@ -59,7 +59,7 @@ bool SplineFields2Dmx::valueCoor (const Vec4& x, Vector& vals) const
// Use with caution, very slow!
Go::Point pt(x.x,x.y,x.z), clopt(3);
double clo_u, clo_v, dist;
const Go::SplineSurface* geo = surf->getBasis(ASMmxBase::geoBasis);
const Go::SplineSurface* geo = surf->getBasis(ASMmxBase::itgBasis);
if (!geo)
geo = surf->getBasis(1);
#pragma omp critical
@ -109,7 +109,7 @@ bool SplineFields2Dmx::gradFE (const ItgPoint& x, Matrix& grad) const
// Evaluate the basis functions at the given point
Matrix Xnod, Jac, dNdX;
std::vector<int> ip;
const Go::SplineSurface* geo = surf->getBasis(ASMmxBase::geoBasis);
const Go::SplineSurface* geo = surf->getBasis(ASMmxBase::itgBasis);
if (!geo)
geo = surf->getBasis(1);
if (!SplineField::evalMapping(*geo,surf->getNoSpaceDim(),x,ip,Xnod,Jac,dNdX))
@ -146,7 +146,7 @@ bool SplineFields2Dmx::hessianFE (const ItgPoint& x, Matrix3D& H) const
Matrix Xnod, Jac, dNdX;
Matrix3D d2NdX2, Hess;
std::vector<int> ip;
const Go::SplineSurface* geo = surf->getBasis(ASMmxBase::geoBasis);
const Go::SplineSurface* geo = surf->getBasis(ASMmxBase::itgBasis);
if (!geo)
geo = surf->getBasis(1);
if (!SplineField::evalMapping(*geo,surf->getNoSpaceDim(),x,ip,Xnod,Jac,dNdX,&d2NdX2,&Hess))

6
src/ASM/SplineFields3Dmx.C
View File

@ -59,7 +59,7 @@ bool SplineFields3Dmx::valueCoor (const Vec4& x, Vector& vals) const
// Use with caution, very slow!
Go::Point pt(x.x,x.y,x.z), clopt(3);
double clo_u, clo_v, clo_w, dist;
const Go::SplineVolume* geo = svol->getBasis(ASMmxBase::geoBasis);
const Go::SplineVolume* geo = svol->getBasis(ASMmxBase::itgBasis);
if (!geo)
geo = svol->getBasis(1);
#pragma omp critical
@ -110,7 +110,7 @@ bool SplineFields3Dmx::gradFE (const ItgPoint& x, Matrix& grad) const
// Evaluate the basis functions at the given point
Matrix Xnod, Jac, dNdX;
std::vector<int> ip;
const Go::SplineVolume* geo = svol->getBasis(ASMmxBase::geoBasis);
const Go::SplineVolume* geo = svol->getBasis(ASMmxBase::itgBasis);
if (!geo)
geo = svol->getBasis(1);
if (!SplineField::evalMapping(*geo,x,ip,Xnod,Jac,dNdX))
@ -146,7 +146,7 @@ bool SplineFields3Dmx::hessianFE (const ItgPoint& x, Matrix3D& H) const
Matrix Xnod, Jac, dNdX;
Matrix3D d2NdX2, Hess;
std::vector<int> ip;
const Go::SplineVolume* geo = svol->getBasis(ASMmxBase::geoBasis);
const Go::SplineVolume* geo = svol->getBasis(ASMmxBase::itgBasis);
if (!geo)
geo = svol->getBasis(1);
if (!SplineField::evalMapping(*geo,x,ip,Xnod,Jac,dNdX,&d2NdX2,&Hess))

2
src/LinAlg/Test/TestISTLMatrix.C
View File

@ -69,7 +69,7 @@ TEST(TestISTLMatrix, Assemble)
TEST(TestISTLMatrix, AssembleBasisBlocks)
{
ASMmxBase::Type = ASMmxBase::FULL_CONT_RAISE_BASIS1;
ASMmxBase::geoBasis = 2;
ASMmxBase::itgBasis = 2;
SIM2D sim({1,1});
sim.read("src/LinAlg/Test/refdata/petsc_test_blocks_basis.xinp");
sim.opt.solver = LinAlg::ISTL;

2
src/LinAlg/Test/TestPETScMatrix.C
View File

@ -81,7 +81,7 @@ TEST(TestPETScMatrix, Assemble)
TEST(TestPETScMatrix, AssembleBasisBlocks)
{
ASMmxBase::Type = ASMmxBase::FULL_CONT_RAISE_BASIS1;
ASMmxBase::geoBasis = 2;
ASMmxBase::itgBasis = 2;
SIM2D sim({1,1});
sim.read("src/LinAlg/Test/refdata/petsc_test_blocks_basis.xinp");
sim.opt.solver = LinAlg::PETSC;

2
src/LinAlg/Test/TestSAM.C
View File

@ -94,7 +94,7 @@ TEST(TestSAM, SingleBasisDirichlet1P)
TEST(TestSAM, MixedBasis1P)
{
ASMmxBase::Type = ASMmxBase::FULL_CONT_RAISE_BASIS1;
ASMmxBase::geoBasis = 2;
ASMmxBase::itgBasis = 2;
SIM2D sim({1,1});
sim.read("src/LinAlg/Test/refdata/sam_2D_1P.xinp");

2
src/SIM/SIMbase.C
View File

@ -845,7 +845,7 @@ bool SIMbase::updateGrid (const RealArray& displ)
Vector locdisp;
if (this->mixedProblem())
this->extractPatchSolution(displ,locdisp,pch,
pch->getNoSpaceDim(),ASMmxBase::geoBasis);
pch->getNoSpaceDim(),ASMmxBase::itgBasis);
else
pch->extractNodeVec(displ,locdisp,pch->getNoSpaceDim(),-1);

2
src/SIM/Test/TestSIM.C
View File

@ -159,7 +159,7 @@ TEST(TestSIM3D, ProjectSolutionMixed)
TEST(TestSIM2D, InjectPatchSolution)
{
ASMmxBase::Type = ASMmxBase::REDUCED_CONT_RAISE_BASIS1;
ASMmxBase::geoBasis = 2;
ASMmxBase::itgBasis = 2;
TestProjectSIM<SIM2D> sim({1,1});
ASMbase* pch = sim.getPatch(1);

6
src/Utility/Test/TestCoordinateMapping.C
View File

@ -305,14 +305,14 @@ TEST(TestCoordinateMapping, Hessian2D_mixed)
Matrix Jac;
std::array<Matrix, 2> grad;
int geoBasis = ASMmxBase::geoBasis - 1;
int itgBasis = ASMmxBase::itgBasis - 1;
Matrix Xnod;
p.getElementCoordinates(Xnod,1);
utl::Jacobian(Jac, grad[geoBasis], Xnod, dNxdu[geoBasis]);
utl::Jacobian(Jac, grad[itgBasis], Xnod, dNxdu[itgBasis]);
grad[1-geoBasis].multiply(dNxdu[1-geoBasis],Jac);
grad[1-itgBasis].multiply(dNxdu[1-itgBasis],Jac);
const double Hess_basis1[] = {
0.5,-1.0, 0.5, 1.0,-2.0, 1.0, 0.5,-1.0, 0.5,