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

Changed: Some cosmetics in [extract|inject]NodVecMx to clarify logic.

Browse Source 
Added: Array elem_size as member in ASMmxBase so we don't have to
recalculate it over and over again.
This commit is contained in:
Knut Morten Okstad 2018-07-24 19:46:56 +02:00
parent 4b804a98de
commit 1f6c033861
4 changed files with 131 additions and 172 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

119
src/ASM/ASMmxBase.C
View File

@ -18,93 +18,88 @@
#include "GoTools/trivariate/VolumeInterpolator.h" #include "GoTools/trivariate/VolumeInterpolator.h"
#include <numeric> #include <numeric>
char ASMmxBase::geoBasis = 2; char ASMmxBase::geoBasis = 2;
ASMmxBase::MixedType ASMmxBase::Type = ASMmxBase::FULL_CONT_RAISE_BASIS1; ASMmxBase::MixedType ASMmxBase::Type = ASMmxBase::FULL_CONT_RAISE_BASIS1;
ASMmxBase::ASMmxBase (const std::vector<unsigned char>& n_f)
{
nfx = n_f;
}
void ASMmxBase::initMx (const std::vector<int>& MLGN, const int* sysMadof) void ASMmxBase::initMx (const std::vector<int>& MLGN, const int* sysMadof)
{ {
MADOF.resize(MLGN.size()); MADOF.clear();
for (size_t i = 0; i < MADOF.size(); i++) MADOF.reserve(MLGN.size());
MADOF[i] = sysMadof[MLGN[i]-1]-1; for (int n : MLGN)
MADOF.push_back(sysMadof[n-1]-1);
} }
void ASMmxBase::extractNodeVecMx (const Vector& globRes, Vector& nodeVec, void ASMmxBase::extractNodeVecMx (const Vector& glbVec, Vector& nodVec,
int basis) const int basis) const
{ {
if (basis > (int)nfx.size()) int b0 = basis-1, b1 = basis;
basis = 0; size_t ofs = 0, len = 0;
if (basis < 1 || basis > (int)nfx.size())
size_t len=0; {
if (basis == 0) b0 = 0;
for (size_t i=0;i<nfx.size();++i) b1 = nfx.size();
for (int i = b0; i < b1; i++)
len += nfx[i]*nb[i]; len += nfx[i]*nb[i];
else
len = nfx[basis-1]*nb[basis-1];
nodeVec.resize(len);
size_t i, j;
int idof, ldof = 0;
size_t k=basis==0?1:basis;
size_t ofs = std::accumulate(nb.begin(), nb.begin()+k-1, 0);
for (; k < (basis==0?nfx.size()+1:(size_t)basis+1); ++k) {
for (i = ofs; i < nb[k-1]+ofs; i++)
{
idof = MADOF[i];
for (j = 0; j < nfx[k-1]; j++, ldof++)
nodeVec[ldof] = globRes[idof++];
}
ofs += nb[k-1];
} }
else
{
len = nfx[basis-1]*nb[basis-1];
if (basis > 1)
ofs = std::accumulate(nb.begin(), nb.begin()+basis-1, 0);
}
nodVec.resize(len);
int ldof = 0;
for (int b = b0; b < b1; b++, ofs += nb[b])
for (size_t i = ofs; i < nb[b]+ofs; i++)
{
int idof = MADOF[i];
for (size_t j = 0; j < nfx[b]; j++)
nodVec[ldof++] = glbVec[idof++];
}
} }
void ASMmxBase::injectNodeVecMx (Vector& globRes, const Vector& nodeVec, void ASMmxBase::injectNodeVecMx (Vector& glbVec, const Vector& nodVec,
int basis) const int basis) const
{ {
if (basis > (int)nfx.size()) int b0 = basis-1, b1 = basis;
basis = 0; size_t ofs = 0;
if (basis < 1 || basis > (int)nfx.size())
size_t i, j; {
int idof, ldof = 0; b0 = 0;
size_t k=basis==0?1:basis; b1 = nfx.size();
size_t ofs = std::accumulate(nb.begin(), nb.begin()+k-1, 0);
for (; k < (basis==0?nfx.size()+1:(size_t)basis+1); ++k) {
for (i = ofs; i < nb[k-1]+ofs; i++)
{
idof = MADOF[i];
for (j = 0; j < nfx[k-1]; j++, ldof++)
globRes[idof++] = nodeVec[ldof];
}
} }
else if (basis > 1)
ofs = std::accumulate(nb.begin(), nb.begin()+basis-1, 0);
int ldof = 0;
for (int b = b0; b < b1; b++)
for (size_t i = ofs; i < nb[b]+ofs; i++)
{
int idof = MADOF[i];
for (size_t j = 0; j < nfx[b]; j++)
glbVec[idof++] = nodVec[ldof++];
}
} }
bool ASMmxBase::getSolutionMx (Matrix& sField, const Vector& locSol, bool ASMmxBase::getSolutionMx (Matrix& sField, const Vector& locSol,
const std::vector<int>& nodes) const const std::vector<int>& nodes) const
{ {
if (nodes.empty()) return true; if (nodes.empty()) return true;
int low, high, nvar; int nvar, low = 1, high = nb.front();
if ((size_t)nodes.front() <= nb[0]) if (nodes.front() <= high)
{
nvar = nfx[0]; nvar = nfx[0];
low = 1;
high = nb[0];
}
else else
{ {
nvar = nfx[1]; nvar = nfx[1];
low = nb[0]+1; low += nb[0];
high = nb[0]+nb[1]; high += nb[1];
} }
sField.resize(nvar,nodes.size()); sField.resize(nvar,nodes.size());
@ -112,15 +107,15 @@ bool ASMmxBase::getSolutionMx (Matrix& sField, const Vector& locSol,
if (nodes[i] < low || nodes[i] > high) if (nodes[i] < low || nodes[i] > high)
{ {
std::cerr <<" *** ASMmxBase::getSolutionMx: Node #"<< nodes[i] std::cerr <<" *** ASMmxBase::getSolutionMx: Node #"<< nodes[i]
<<" is out of range ["<< low <<","<< high <<"]."<< std::endl; <<" is out of range ["<< low <<","<< high <<"]."<< std::endl;
return false; return false;
} }
else else
{ {
int idof = nvar*(nodes[i]-1); int idof = nvar*(nodes[i]-1);
if (low > 1) idof += nfx[0]*nb[0]; if (low > 1) idof += nfx.front()*nb.front();
for (int j = 0; j < nvar; j++) for (int j = 0; j < nvar; j++)
sField(j+1,i+1) = locSol[idof++]; sField(j+1,i+1) = locSol[idof++];
} }
return true; return true;

24
src/ASM/ASMmxBase.h
View File

@ -32,7 +32,7 @@ class ASMmxBase
protected: protected:
//! \brief The constructor sets the number of field variables. //! \brief The constructor sets the number of field variables.
//! \param[in] n_f Number of nodal variables in each field //! \param[in] n_f Number of nodal variables in each field
explicit ASMmxBase(const std::vector<unsigned char>& n_f); explicit ASMmxBase(const std::vector<unsigned char>& n_f) : nfx(n_f) {}
//! \brief Initializes the patch level MADOF array. //! \brief Initializes the patch level MADOF array.
//! \param[in] MLGN Matrix of local-to-global node numbers //! \param[in] MLGN Matrix of local-to-global node numbers
@ -40,25 +40,23 @@ protected:
void initMx(const std::vector<int>& MLGN, const int* sysMadof); void initMx(const std::vector<int>& MLGN, const int* sysMadof);
//! \brief Extracts nodal results for this patch from the global vector. //! \brief Extracts nodal results for this patch from the global vector.
//! \param[in] globVec Global solution vector in DOF-order //! \param[in] glbVec Global solution vector in DOF-order
//! \param[out] nodeVec Nodal result vector for this patch //! \param[out] nodVec Nodal result vector for this patch
//! \param[in] basis Which basis to extract the nodal values for //! \param[in] basis Which basis to extract the nodal values for
void extractNodeVecMx(const Vector& globVec, Vector& nodeVec, void extractNodeVecMx(const Vector& glbVec, Vector& nodVec, int basis) const;
int basis = 0) const;
//! \brief Injects nodal results for this patch into a global vector. //! \brief Injects nodal results for this patch into a global vector.
//! \param[out] globVec Global solution vector in DOF-order //! \param[out] glbVec Global solution vector in DOF-order
//! \param[in] nodeVec Nodal result vector for this patch //! \param[in] nodVec Nodal result vector for this patch
//! \param[in] basis Which basis to inject the nodal values for //! \param[in] basis Which basis to inject the nodal values for
void injectNodeVecMx(Vector& globVec, const Vector& nodeVec, void injectNodeVecMx(Vector& glbVec, const Vector& nodVec, int basis) const;
int basis = 0) const;
//! \brief Extracts the primary solution field at the specified nodes. //! \brief Extracts the primary solution field at the specified nodes.
//! \param[out] sField Solution field //! \param[out] sField Solution field
//! \param[in] locSol Solution vector local to current patch //! \param[in] locSol Solution vector local to current patch
//! \param[in] nodes 1-based local node numbers to extract solution for //! \param[in] nodes 1-based local node numbers to extract solution for
bool getSolutionMx(Matrix& sField, const Vector& locSol, bool getSolutionMx(Matrix& sField, const Vector& locSol,
const std::vector<int>& nodes) const; const std::vector<int>& nodes) const;
public: public:
//! \brief Enum defining available mixed formulation types. //! \brief Enum defining available mixed formulation types.
@ -95,7 +93,11 @@ private:
std::vector<int> MADOF; //!< Matrix of accumulated DOFs for this patch std::vector<int> MADOF; //!< Matrix of accumulated DOFs for this patch
protected: protected:
std::vector<size_t> nb; //!< Number of basis functions in each basis //! Number of basis functions per element in each basis
std::vector<size_t> elem_size;
//! Total number of basis functions in each basis
std::vector<size_t> nb;
std::vector<unsigned char> nfx; //!< Number of fields on each basis std::vector<unsigned char> nfx; //!< Number of fields on each basis
}; };

79
src/ASM/ASMs2Dmx.C
View File

@ -35,15 +35,13 @@
#endif #endif
ASMs2Dmx::ASMs2Dmx (unsigned char n_s, ASMs2Dmx::ASMs2Dmx (unsigned char n_s, const CharVec& n_f)
const std::vector<unsigned char>& n_f)
: ASMs2D(n_s, *std::max_element(n_f.begin(),n_f.end())), ASMmxBase(n_f) : ASMs2D(n_s, *std::max_element(n_f.begin(),n_f.end())), ASMmxBase(n_f)
{ {
} }
ASMs2Dmx::ASMs2Dmx (const ASMs2Dmx& patch, ASMs2Dmx::ASMs2Dmx (const ASMs2Dmx& patch, const CharVec& n_f)
const std::vector<unsigned char>& n_f)
: ASMs2D(patch), ASMmxBase(n_f) : ASMs2D(patch), ASMmxBase(n_f)
{ {
nb = patch.nb; nb = patch.nb;
@ -191,8 +189,13 @@ bool ASMs2Dmx::generateFEMTopology ()
geo = surf = m_basis[geoBasis-1]->clone(); geo = surf = m_basis[geoBasis-1]->clone();
nb.clear(); nb.clear();
for (auto it : m_basis) nb.reserve(m_basis.size());
elem_size.clear();
elem_size.reserve(m_basis.size());
for (auto& it : m_basis) {
nb.push_back(it->numCoefs_u()*it->numCoefs_v()); nb.push_back(it->numCoefs_u()*it->numCoefs_v());
elem_size.push_back(it->order_u()*it->order_v());
}
if (!nodeInd.empty() && !shareFE) if (!nodeInd.empty() && !shareFE)
{ {
@ -440,12 +443,8 @@ double ASMs2Dmx::getParametricArea (int iel) const
if (MNPC[iel-1].empty()) if (MNPC[iel-1].empty())
return 0.0; return 0.0;
std::vector<size_t> elem_sizes; int inod1 = MNPC[iel-1][std::accumulate(elem_size.begin(),
for (auto& it : m_basis) elem_size.begin()+geoBasis, -1)];
elem_sizes.push_back(it->order_u()*it->order_v());
int inod1 = MNPC[iel-1][std::accumulate(elem_sizes.begin(),
elem_sizes.begin()+geoBasis, -1)];
#ifdef INDEX_CHECK #ifdef INDEX_CHECK
if (inod1 < 0 || (size_t)inod1 >= nnod) if (inod1 < 0 || (size_t)inod1 >= nnod)
{ {
@ -475,12 +474,8 @@ double ASMs2Dmx::getParametricLength (int iel, int dir) const
if (MNPC[iel-1].empty()) if (MNPC[iel-1].empty())
return 0.0; return 0.0;
std::vector<size_t> elem_sizes; int inod1 = MNPC[iel-1][std::accumulate(elem_size.begin(),
for (auto& it : m_basis) elem_size.begin()+geoBasis, -1)];
elem_sizes.push_back(it->order_u()*it->order_v());
int inod1 = MNPC[iel-1][std::accumulate(elem_sizes.begin(),
elem_sizes.begin()+geoBasis, -1)];
#ifdef INDEX_CHECK #ifdef INDEX_CHECK
if (inod1 < 0 || (size_t)inod1 >= nnod) if (inod1 < 0 || (size_t)inod1 >= nnod)
{ {
@ -537,9 +532,6 @@ bool ASMs2Dmx::integrate (Integrand& integrand,
const int n1 = surf->numCoefs_u(); const int n1 = surf->numCoefs_u();
const int nel1 = n1 - p1 + 1; const int nel1 = n1 - p1 + 1;
std::vector<size_t> elem_sizes;
for (auto& it : m_basis)
elem_sizes.push_back(it->order_u()*it->order_v());
// === Assembly loop over all elements in the patch ========================== // === Assembly loop over all elements in the patch ==========================
@ -547,7 +539,7 @@ bool ASMs2Dmx::integrate (Integrand& integrand,
for (size_t g=0;g<threadGroups.size() && ok;++g) { for (size_t g=0;g<threadGroups.size() && ok;++g) {
#pragma omp parallel for schedule(static) #pragma omp parallel for schedule(static)
for (size_t t=0;t<threadGroups[g].size();++t) { for (size_t t=0;t<threadGroups[g].size();++t) {
MxFiniteElement fe(elem_sizes); MxFiniteElement fe(elem_size);
std::vector<Matrix> dNxdu(m_basis.size()); std::vector<Matrix> dNxdu(m_basis.size());
std::vector<Matrix3D> d2Nxdu2(m_basis.size()); std::vector<Matrix3D> d2Nxdu2(m_basis.size());
Matrix3D Hess; Matrix3D Hess;
@ -590,8 +582,8 @@ bool ASMs2Dmx::integrate (Integrand& integrand,
} }
// Initialize element quantities // Initialize element quantities
LocalIntegral* A = integrand.getLocalIntegral(elem_sizes,fe.iel,false); LocalIntegral* A = integrand.getLocalIntegral(elem_size,fe.iel,false);
if (!integrand.initElement(MNPC[iel-1],elem_sizes,nb,*A)) if (!integrand.initElement(MNPC[iel-1],elem_size,nb,*A))
{ {
A->destruct(); A->destruct();
ok = false; ok = false;
@ -726,19 +718,15 @@ bool ASMs2Dmx::integrate (Integrand& integrand, int lIndex,
const int n1 = surf->numCoefs_u(); const int n1 = surf->numCoefs_u();
const int n2 = surf->numCoefs_v(); const int n2 = surf->numCoefs_v();
std::vector<size_t> elem_sizes;
for (auto& it : m_basis)
elem_sizes.push_back(it->order_u()*it->order_v());
std::map<char,size_t>::const_iterator iit = firstBp.find(lIndex%10); std::map<char,size_t>::const_iterator iit = firstBp.find(lIndex%10);
size_t firstp = iit == firstBp.end() ? 0 : iit->second; size_t firstp = iit == firstBp.end() ? 0 : iit->second;
MxFiniteElement fe(elem_sizes); MxFiniteElement fe(elem_size);
fe.xi = fe.eta = edgeDir < 0 ? -1.0 : 1.0; fe.xi = fe.eta = edgeDir < 0 ? -1.0 : 1.0;
fe.u = gpar[0](1,1); fe.u = gpar[0](1,1);
fe.v = gpar[1](1,1); fe.v = gpar[1](1,1);
std::vector<Matrix> dNxdu(m_basis.size()); Matrices dNxdu(m_basis.size());
Matrix Xnod, Jac; Matrix Xnod, Jac;
Vec4 X; Vec4 X;
Vec3 normal; Vec3 normal;
@ -775,8 +763,8 @@ bool ASMs2Dmx::integrate (Integrand& integrand, int lIndex,
fe.h = this->getElementCorners(i1-1,i2-1,fe.XC); fe.h = this->getElementCorners(i1-1,i2-1,fe.XC);
// Initialize element quantities // Initialize element quantities
LocalIntegral* A = integrand.getLocalIntegral(elem_sizes,fe.iel,true); LocalIntegral* A = integrand.getLocalIntegral(elem_size,fe.iel,true);
bool ok = integrand.initElementBou(MNPC[iel-1],elem_sizes,nb,*A); bool ok = integrand.initElementBou(MNPC[iel-1],elem_size,nb,*A);
// --- Integration loop over all Gauss points along the edge ------------- // --- Integration loop over all Gauss points along the edge -------------
@ -859,11 +847,8 @@ bool ASMs2Dmx::integrate (Integrand& integrand,
const int n1 = surf->numCoefs_u(); const int n1 = surf->numCoefs_u();
const int n2 = surf->numCoefs_v(); const int n2 = surf->numCoefs_v();
std::vector<size_t> elem_sizes; std::vector<size_t> elem_sizes2(elem_size);
for (auto& it : m_basis) std::copy(elem_size.begin(), elem_size.end(), std::back_inserter(elem_sizes2));
elem_sizes.push_back(it->order_u()*it->order_v());
std::vector<size_t> elem_sizes2(elem_sizes);
std::copy(elem_sizes.begin(), elem_sizes.end(), std::back_inserter(elem_sizes2));
MxFiniteElement fe(elem_sizes2); MxFiniteElement fe(elem_sizes2);
Matrix dNdu, Xnod, Jac; Matrix dNdu, Xnod, Jac;
@ -903,8 +888,8 @@ bool ASMs2Dmx::integrate (Integrand& integrand,
fe.h = this->getElementCorners(i1-1,i2-1,fe.XC); fe.h = this->getElementCorners(i1-1,i2-1,fe.XC);
// Initialize element quantities // Initialize element quantities
LocalIntegral* A = integrand.getLocalIntegral(elem_sizes,fe.iel); LocalIntegral* A = integrand.getLocalIntegral(elem_size,fe.iel);
bool ok = integrand.initElement(MNPC[iel],elem_sizes,nb,*A); bool ok = integrand.initElement(MNPC[iel],elem_size,nb,*A);
size_t origSize = A->vec.size(); size_t origSize = A->vec.size();
// Loop over the element edges with contributions // Loop over the element edges with contributions
@ -924,8 +909,8 @@ bool ASMs2Dmx::integrate (Integrand& integrand,
kel += iedge > 3 ? n1-p1+1 : -(n1-p1+1); kel += iedge > 3 ? n1-p1+1 : -(n1-p1+1);
// initialize neighbor element // initialize neighbor element
LocalIntegral* A_neigh = integrand.getLocalIntegral(elem_sizes,kel+1); LocalIntegral* A_neigh = integrand.getLocalIntegral(elem_size,kel+1);
ok &= integrand.initElement(MNPC[kel],elem_sizes,nb,*A_neigh); ok &= integrand.initElement(MNPC[kel],elem_size,nb,*A_neigh);
if (!A_neigh->vec.empty()) { if (!A_neigh->vec.empty()) {
A->vec.resize(origSize+A_neigh->vec.size()); A->vec.resize(origSize+A_neigh->vec.size());
std::copy(A_neigh->vec.begin(), A_neigh->vec.end(), A->vec.begin()+origSize); std::copy(A_neigh->vec.begin(), A_neigh->vec.end(), A->vec.begin()+origSize);
@ -961,7 +946,7 @@ bool ASMs2Dmx::integrate (Integrand& integrand,
} }
// Fetch basis function derivatives at current integration point // Fetch basis function derivatives at current integration point
std::vector<Matrix> dNxdu(m_basis.size()*2); Matrices dNxdu(m_basis.size()*2);
for (size_t b = 0; b < m_basis.size(); ++b) { for (size_t b = 0; b < m_basis.size(); ++b) {
Go::BasisDerivsSf spline; Go::BasisDerivsSf spline;
this->getBasis(b+1)->computeBasis(fe.u, fe.v, spline, edgeDir < 0); this->getBasis(b+1)->computeBasis(fe.u, fe.v, spline, edgeDir < 0);
@ -1120,17 +1105,13 @@ bool ASMs2Dmx::evalSolution (Matrix& sField, const IntegrandBase& integrand,
} }
} }
std::vector<size_t> elem_sizes;
for (auto& it : m_basis)
elem_sizes.push_back(it->order_u()*it->order_v());
// Fetch nodal (control point) coordinates // Fetch nodal (control point) coordinates
Matrix Xnod, Xtmp; Matrix Xnod, Xtmp;
this->getNodalCoordinates(Xnod); this->getNodalCoordinates(Xnod);
MxFiniteElement fe(elem_sizes,firstIp); MxFiniteElement fe(elem_size,firstIp);
Vector solPt; Vector solPt;
std::vector<Matrix> dNxdu(m_basis.size()); Matrices dNxdu(m_basis.size());
Matrix Jac; Matrix Jac;
Vec3 X; Vec3 X;
@ -1139,7 +1120,7 @@ bool ASMs2Dmx::evalSolution (Matrix& sField, const IntegrandBase& integrand,
for (size_t i = 0; i < nPoints; i++, fe.iGP++) for (size_t i = 0; i < nPoints; i++, fe.iGP++)
{ {
// Fetch indices of the non-zero basis functions at this point // Fetch indices of the non-zero basis functions at this point
std::vector<IntVec> ip(m_basis.size()); IntMat ip(m_basis.size());
IntVec ipa; IntVec ipa;
size_t ofs = 0; size_t ofs = 0;
for (size_t b = 0; b < m_basis.size(); ++b) { for (size_t b = 0; b < m_basis.size(); ++b) {
@ -1180,7 +1161,7 @@ bool ASMs2Dmx::evalSolution (Matrix& sField, const IntegrandBase& integrand,
X = Xtmp * fe.basis(geoBasis); X = Xtmp * fe.basis(geoBasis);
// Now evaluate the solution field // Now evaluate the solution field
if (!integrand.evalSol(solPt,fe,X,ipa,elem_sizes,nb)) if (!integrand.evalSol(solPt,fe,X,ipa,elem_size,nb))
return false; return false;
else if (sField.empty()) else if (sField.empty())
sField.resize(solPt.size(),nPoints,true); sField.resize(solPt.size(),nPoints,true);

81
src/ASM/ASMs3Dmx.C
View File

@ -35,14 +35,13 @@
#endif #endif
ASMs3Dmx::ASMs3Dmx (const std::vector<unsigned char>& n_f) ASMs3Dmx::ASMs3Dmx (const CharVec& n_f)
: ASMs3D(std::accumulate(n_f.begin(), n_f.end(), 0)), ASMmxBase(n_f) : ASMs3D(std::accumulate(n_f.begin(), n_f.end(), 0)), ASMmxBase(n_f)
{ {
} }
ASMs3Dmx::ASMs3Dmx (const ASMs3Dmx& patch, ASMs3Dmx::ASMs3Dmx (const ASMs3Dmx& patch, const CharVec& n_f)
const std::vector<unsigned char>& n_f)
: ASMs3D(patch), ASMmxBase(n_f) : ASMs3D(patch), ASMmxBase(n_f)
{ {
m_basis = patch.m_basis; m_basis = patch.m_basis;
@ -192,8 +191,13 @@ bool ASMs3Dmx::generateFEMTopology ()
geo = svol = m_basis[geoBasis-1]->clone(); geo = svol = m_basis[geoBasis-1]->clone();
nb.clear(); nb.clear();
for (auto it : m_basis) elem_size.clear();
nb.reserve(m_basis.size());
elem_size.reserve(m_basis.size());
for (auto& it : m_basis) {
nb.push_back(it->numCoefs(0)*it->numCoefs(1)*it->numCoefs(2)); nb.push_back(it->numCoefs(0)*it->numCoefs(1)*it->numCoefs(2));
elem_size.push_back(it->order(0)*it->order(1)*it->order(2));
}
if (!nodeInd.empty() && !shareFE) if (!nodeInd.empty() && !shareFE)
{ {
@ -474,10 +478,6 @@ bool ASMs3Dmx::integrate (Integrand& integrand,
for (size_t i=0;i<m_basis.size();++i) for (size_t i=0;i<m_basis.size();++i)
m_basis[i]->computeBasisGrid(gpar[0],gpar[1],gpar[2],splinex[i]); m_basis[i]->computeBasisGrid(gpar[0],gpar[1],gpar[2],splinex[i]);
std::vector<size_t> elem_sizes;
for (auto& it : m_basis)
elem_sizes.push_back(it->order(0)*it->order(1)*it->order(2));
const int p1 = svol->order(0); const int p1 = svol->order(0);
const int p2 = svol->order(1); const int p2 = svol->order(1);
const int p3 = svol->order(2); const int p3 = svol->order(2);
@ -493,7 +493,7 @@ bool ASMs3Dmx::integrate (Integrand& integrand,
for (size_t g=0;g<threadGroupsVol.size() && ok;++g) { for (size_t g=0;g<threadGroupsVol.size() && ok;++g) {
#pragma omp parallel for schedule(static) #pragma omp parallel for schedule(static)
for (size_t t=0;t<threadGroupsVol[g].size();++t) { for (size_t t=0;t<threadGroupsVol[g].size();++t) {
MxFiniteElement fe(elem_sizes); MxFiniteElement fe(elem_size);
std::vector<Matrix> dNxdu(m_basis.size()); std::vector<Matrix> dNxdu(m_basis.size());
std::vector<Matrix3D> d2Nxdu2(m_basis.size()); std::vector<Matrix3D> d2Nxdu2(m_basis.size());
Matrix3D Hess; Matrix3D Hess;
@ -538,8 +538,8 @@ bool ASMs3Dmx::integrate (Integrand& integrand,
} }
// Initialize element quantities // Initialize element quantities
LocalIntegral* A = integrand.getLocalIntegral(elem_sizes,fe.iel,false); LocalIntegral* A = integrand.getLocalIntegral(elem_size,fe.iel,false);
if (!integrand.initElement(MNPC[iel-1],elem_sizes,nb,*A)) if (!integrand.initElement(MNPC[iel-1],elem_size,nb,*A))
{ {
A->destruct(); A->destruct();
ok = false; ok = false;
@ -690,10 +690,6 @@ bool ASMs3Dmx::integrate (Integrand& integrand, int lIndex,
const int p2 = svol->order(1); const int p2 = svol->order(1);
const int p3 = svol->order(2); const int p3 = svol->order(2);
std::vector<size_t> elem_sizes;
for (auto& it : m_basis)
elem_sizes.push_back(it->order(0)*it->order(1)*it->order(2));
const int nel1 = n1 - p1 + 1; const int nel1 = n1 - p1 + 1;
const int nel2 = n2 - p2 + 1; const int nel2 = n2 - p2 + 1;
@ -707,15 +703,15 @@ bool ASMs3Dmx::integrate (Integrand& integrand, int lIndex,
for (size_t g = 0; g < threadGrp.size() && ok; ++g) { for (size_t g = 0; g < threadGrp.size() && ok; ++g) {
#pragma omp parallel for schedule(static) #pragma omp parallel for schedule(static)
for (size_t t = 0; t < threadGrp[g].size(); ++t) { for (size_t t = 0; t < threadGrp[g].size(); ++t) {
MxFiniteElement fe(elem_sizes); MxFiniteElement fe(elem_size);
fe.xi = fe.eta = fe.zeta = faceDir < 0 ? -1.0 : 1.0; fe.xi = fe.eta = fe.zeta = faceDir < 0 ? -1.0 : 1.0;
fe.u = gpar[0](1,1); fe.u = gpar[0](1,1);
fe.v = gpar[1](1,1); fe.v = gpar[1](1,1);
fe.w = gpar[2](1,1); fe.w = gpar[2](1,1);
std::vector<Matrix> dNxdu(m_basis.size()); Matrices dNxdu(m_basis.size());
Matrix Xnod, Jac; Matrix Xnod, Jac;
double param[3] = { fe.u, fe.v, fe.w }; double param[3] = { fe.u, fe.v, fe.w };
Vec4 X; Vec4 X;
Vec3 normal; Vec3 normal;
for (size_t l = 0; l < threadGrp[g][t].size() && ok; ++l) for (size_t l = 0; l < threadGrp[g][t].size() && ok; ++l)
@ -747,8 +743,8 @@ bool ASMs3Dmx::integrate (Integrand& integrand, int lIndex,
fe.h = this->getElementCorners(i1-1,i2-1,i3-1,fe.XC); fe.h = this->getElementCorners(i1-1,i2-1,i3-1,fe.XC);
// Initialize element quantities // Initialize element quantities
LocalIntegral* A = integrand.getLocalIntegral(elem_sizes,fe.iel,true); LocalIntegral* A = integrand.getLocalIntegral(elem_size,fe.iel,true);
if (!integrand.initElementBou(MNPC[iel-1],elem_sizes,nb,*A)) if (!integrand.initElementBou(MNPC[iel-1],elem_size,nb,*A))
{ {
A->destruct(); A->destruct();
ok = false; ok = false;
@ -867,11 +863,8 @@ bool ASMs3Dmx::integrate (Integrand& integrand,
const int nel1 = n1 - p1 + 1; const int nel1 = n1 - p1 + 1;
const int nel2 = n2 - p2 + 1; const int nel2 = n2 - p2 + 1;
std::vector<size_t> elem_sizes; std::vector<size_t> elem_sizes2(elem_size);
for (auto& it : m_basis) std::copy(elem_size.begin(), elem_size.end(), std::back_inserter(elem_sizes2));
elem_sizes.push_back(it->order(0)*it->order(1)*it->order(2));
std::vector<size_t> elem_sizes2(elem_sizes);
std::copy(elem_sizes.begin(), elem_sizes.end(), std::back_inserter(elem_sizes2));
MxFiniteElement fe(elem_sizes2); MxFiniteElement fe(elem_sizes2);
Matrix dNdu, Xnod, Jac; Matrix dNdu, Xnod, Jac;
@ -912,8 +905,8 @@ bool ASMs3Dmx::integrate (Integrand& integrand,
fe.h = this->getElementCorners(i1-1,i2-1,i3-1,fe.XC); fe.h = this->getElementCorners(i1-1,i2-1,i3-1,fe.XC);
// Initialize element quantities // Initialize element quantities
LocalIntegral* A = integrand.getLocalIntegral(elem_sizes,fe.iel); LocalIntegral* A = integrand.getLocalIntegral(elem_size,fe.iel);
bool ok = integrand.initElement(MNPC[iel],elem_sizes,nb,*A); bool ok = integrand.initElement(MNPC[iel],elem_size,nb,*A);
size_t origSize = A->vec.size(); size_t origSize = A->vec.size();
// Loop over the element edges with contributions // Loop over the element edges with contributions
@ -943,8 +936,8 @@ bool ASMs3Dmx::integrate (Integrand& integrand,
kel += iface == 6 ? (n2-p2+1)*(n1-p1+1) : -(n2-p2+1)*(n1-p1+1); kel += iface == 6 ? (n2-p2+1)*(n1-p1+1) : -(n2-p2+1)*(n1-p1+1);
// initialize neighbor element // initialize neighbor element
LocalIntegral* A_neigh = integrand.getLocalIntegral(elem_sizes,kel+1); LocalIntegral* A_neigh = integrand.getLocalIntegral(elem_size,kel+1);
ok &= integrand.initElement(MNPC[kel],elem_sizes,nb,*A_neigh); ok &= integrand.initElement(MNPC[kel],elem_size,nb,*A_neigh);
if (!A_neigh->vec.empty()) { if (!A_neigh->vec.empty()) {
A->vec.resize(origSize+A_neigh->vec.size()); A->vec.resize(origSize+A_neigh->vec.size());
std::copy(A_neigh->vec.begin(), A_neigh->vec.end(), A->vec.begin()+origSize); std::copy(A_neigh->vec.begin(), A_neigh->vec.end(), A->vec.begin()+origSize);
@ -1007,7 +1000,7 @@ bool ASMs3Dmx::integrate (Integrand& integrand,
} }
// Fetch basis function derivatives at current integration point // Fetch basis function derivatives at current integration point
std::vector<Matrix> dNxdu(m_basis.size()*2); Matrices dNxdu(m_basis.size()*2);
for (size_t b = 0; b < m_basis.size(); ++b) { for (size_t b = 0; b < m_basis.size(); ++b) {
Go::BasisDerivs spline; Go::BasisDerivs spline;
this->getBasis(b+1)->computeBasis(fe.u, fe.v, fe.w, spline, faceDir < 0); this->getBasis(b+1)->computeBasis(fe.u, fe.v, fe.w, spline, faceDir < 0);
@ -1172,17 +1165,13 @@ bool ASMs3Dmx::evalSolution (Matrix& sField, const IntegrandBase& integrand,
} }
} }
std::vector<size_t> elem_sizes;
for (auto& it : m_basis)
elem_sizes.push_back(it->order(0)*it->order(1)*it->order(2));
// Fetch nodal (control point) coordinates // Fetch nodal (control point) coordinates
Matrix Xnod, Xtmp; Matrix Xnod, Xtmp;
this->getNodalCoordinates(Xnod); this->getNodalCoordinates(Xnod);
MxFiniteElement fe(elem_sizes,firstIp); MxFiniteElement fe(elem_size,firstIp);
Vector solPt; Vector solPt;
std::vector<Matrix> dNxdu(m_basis.size()); Matrices dNxdu(m_basis.size());
Matrix Jac; Matrix Jac;
Vec3 X; Vec3 X;
@ -1191,7 +1180,7 @@ bool ASMs3Dmx::evalSolution (Matrix& sField, const IntegrandBase& integrand,
for (size_t i = 0; i < nPoints; i++, fe.iGP++) for (size_t i = 0; i < nPoints; i++, fe.iGP++)
{ {
// Fetch indices of the non-zero basis functions at this point // Fetch indices of the non-zero basis functions at this point
std::vector<IntVec> ip(m_basis.size()); IntMat ip(m_basis.size());
IntVec ipa; IntVec ipa;
size_t ofs = 0; size_t ofs = 0;
for (size_t b = 0; b < m_basis.size(); ++b) { for (size_t b = 0; b < m_basis.size(); ++b) {
@ -1234,7 +1223,7 @@ bool ASMs3Dmx::evalSolution (Matrix& sField, const IntegrandBase& integrand,
X = Xtmp * fe.basis(geoBasis); X = Xtmp * fe.basis(geoBasis);
// Now evaluate the solution field // Now evaluate the solution field
if (!integrand.evalSol(solPt,fe,X,ipa,elem_sizes,nb)) if (!integrand.evalSol(solPt,fe,X,ipa,elem_size,nb))
return false; return false;
else if (sField.empty()) else if (sField.empty())
sField.resize(solPt.size(),nPoints,true); sField.resize(solPt.size(),nPoints,true);
@ -1283,12 +1272,8 @@ double ASMs3Dmx::getParametricVolume (int iel) const
if (MNPC[iel-1].empty()) if (MNPC[iel-1].empty())
return 0.0; return 0.0;
std::vector<size_t> elem_sizes; int inod1 = MNPC[iel-1][std::accumulate(elem_size.begin(),
for (auto& it : m_basis) elem_size.begin()+geoBasis, -1)];
elem_sizes.push_back(it->order(0)*it->order(1)*it->order(2));
int inod1 = MNPC[iel-1][std::accumulate(elem_sizes.begin(),
elem_sizes.begin()+geoBasis, -1)];
#ifdef INDEX_CHECK #ifdef INDEX_CHECK
if (inod1 < 0 || (size_t)inod1 >= nnod) if (inod1 < 0 || (size_t)inod1 >= nnod)
{ {
@ -1318,12 +1303,8 @@ double ASMs3Dmx::getParametricArea (int iel, int dir) const
if (MNPC[iel-1].empty()) if (MNPC[iel-1].empty())
return 0.0; return 0.0;
std::vector<size_t> elem_sizes; int inod1 = MNPC[iel-1][std::accumulate(elem_size.begin(),
for (auto& it : m_basis) elem_size.begin()+geoBasis, -1)];
elem_sizes.push_back(it->order(0)*it->order(1)*it->order(2));
int inod1 = MNPC[iel-1][std::accumulate(elem_sizes.begin(),
elem_sizes.begin()+geoBasis, -1)];
#ifdef INDEX_CHECK #ifdef INDEX_CHECK
if (inod1 < 0 || (size_t)inod1 >= nnod) if (inod1 < 0 || (size_t)inod1 >= nnod)
{ {