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

move projection bases into ASMLRSpline

Browse Source
This commit is contained in:
Arne Morten Kvarving 2023-08-28 15:57:08 +02:00
parent dce989f1c6
commit 505d9d9a6c
12 changed files with 111 additions and 111 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
src/ASM/LR/ASMLRSpline.C
View File

@ -195,6 +195,8 @@ ASMLRSpline::ASMLRSpline (const ASMLRSpline& patch, unsigned char n_f)
: ASMbase(patch,n_f)
{
geomB = patch.geomB;
projB = patch.projB;
projB2 = patch.projB2;
}

4
src/ASM/LR/ASMLRSpline.h
View File

@ -203,7 +203,9 @@ protected:
//! \param[in] groups The generated thread groups
static void analyzeThreadGroups(const IntMat& groups);
std::shared_ptr<LR::LRSpline> geomB; //!< Pointer to the actual spline geometry object
std::shared_ptr<LR::LRSpline> geomB; //!< Pointer to spline object of the geometry basis
std::shared_ptr<LR::LRSpline> projB; //!< Pointer to spline object of the projection basis
std::shared_ptr<LR::LRSpline> projB2; //!< Pointer to spline object of the secondary projection basis
};
#endif

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

@ -56,7 +56,6 @@ ASMu2D::ASMu2D (const ASMu2D& patch, unsigned char n_f)
{
aMin = 0.0;
tensorspline = tensorPrjBas = nullptr;
projBasis = patch.projBasis;
is_rational = patch.is_rational;
// Need to set nnod here,
@ -401,8 +400,8 @@ bool ASMu2D::createProjectionBasis (bool init)
tensorPrjBas = tensorspline->clone();
std::swap(tensorspline,tensorPrjBas);
std::swap(lrspline,projBasis);
geomB = lrspline;
std::swap(geomB,projB);
lrspline = std::static_pointer_cast<LR::LRSplineSurface>(geomB);
return true;
}
@ -551,14 +550,14 @@ bool ASMu2D::generateFEMTopology ()
if (tensorPrjBas)
{
projBasis.reset(tensorPrjBas->rational() ? createLRNurbs(*tensorPrjBas)
: new LR::LRSplineSurface(tensorPrjBas));
projBasis->generateIDs();
projB.reset(tensorPrjBas->rational() ? createLRNurbs(*tensorPrjBas)
: new LR::LRSplineSurface(tensorPrjBas));
projB->generateIDs();
delete tensorPrjBas;
tensorPrjBas = nullptr;
}
else if (!projBasis)
projBasis = lrspline;
else if (!projB)
projB = lrspline;
if (!lrspline) return false;
@ -2395,20 +2394,20 @@ size_t ASMu2D::getNoNodes (int basis) const
size_t ASMu2D::getNoProjectionNodes () const
{
return projBasis->nBasisFunctions();
return projB->nBasisFunctions();
}
bool ASMu2D::separateProjectionBasis () const
{
return projBasis.get() != this->getBasis(1);
return projB.get() != this->getBasis(1);
}
Field* ASMu2D::getProjectedField (const Vector& coefs) const
{
if (coefs.size() == this->getNoProjectionNodes())
return new LRSplineField2D(projBasis.get(),coefs,is_rational);
return new LRSplineField2D(static_cast<const LR::LRSplineSurface*>(projB.get()),coefs,is_rational);
std::cerr <<" *** ASMu2D::getProjectedFields: Non-matching coefficent array,"
<<" size="<< coefs.size() <<" nnod="<< this->getNoProjectionNodes()
@ -2419,12 +2418,12 @@ Field* ASMu2D::getProjectedField (const Vector& coefs) const
Fields* ASMu2D::getProjectedFields (const Vector& coefs, size_t) const
{
if (projBasis.get() == this->getBasis(1))
if (!this->separateProjectionBasis())
return nullptr;
size_t ncmp = coefs.size() / this->getNoProjectionNodes();
if (ncmp*this->getNoProjectionNodes() == coefs.size())
return new LRSplineFields2D(projBasis.get(),coefs,ncmp,is_rational);
return new LRSplineFields2D(static_cast<const LR::LRSplineSurface*>(projB.get()),coefs,ncmp,is_rational);
std::cerr <<" *** ASMu2D::getProjectedFields: Non-matching coefficent array,"
<<" size="<< coefs.size() <<" nnod="<< this->getNoProjectionNodes()
@ -2585,8 +2584,8 @@ void ASMu2D::generateThreadGroups (const Integrand& integrand, bool silence,
bool ignoreGlobalLM)
{
LR::generateThreadGroups(threadGroups, this->getBasis(1));
if (projBasis != lrspline)
LR::generateThreadGroups(projThreadGroups, projBasis.get());
if (this->separateProjectionBasis())
LR::generateThreadGroups(projThreadGroups, projB.get());
if (silence || threadGroups[0].size() < 2) return;
IFEM::cout <<"\nMultiple threads are utilized during element assembly.";
@ -2741,21 +2740,22 @@ bool ASMu2D::refine (const LR::RefineData& prm, Vectors& sol)
prm.elements.size() + prm.errors.size() == 0)
return ok;
LR::LRSplineSurface* proj = static_cast<LR::LRSplineSurface*>(projB.get());
for (const LR::Meshline* line : lrspline->getAllMeshlines())
if (line->span_u_line_)
projBasis->insert_const_v_edge(line->const_par_,
line->start_, line->stop_,
line->multiplicity());
proj->insert_const_v_edge(line->const_par_,
line->start_, line->stop_,
line->multiplicity());
else
projBasis->insert_const_u_edge(line->const_par_,
line->start_, line->stop_,
line->multiplicity());
proj->insert_const_u_edge(line->const_par_,
line->start_, line->stop_,
line->multiplicity());
if (projBasis != lrspline)
projBasis->generateIDs();
projB->generateIDs();
IFEM::cout <<"Refined projection basis: "<< projBasis->nElements()
<<" elements "<< projBasis->nBasisFunctions() <<" nodes."
IFEM::cout <<"Refined projection basis: "<< projB->nElements()
<<" elements "<< projB->nBasisFunctions() <<" nodes."
<< std::endl;
return true;
}

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

@ -742,7 +742,6 @@ public:
protected:
std::shared_ptr<LR::LRSplineSurface> lrspline; //!< Pointer to the LR-spline surface object
std::shared_ptr<LR::LRSplineSurface> projBasis; //!< Basis to project onto
bool is_rational = false; //!< True if basis is rational

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

@ -95,7 +95,7 @@ bool ASMu2Dmx::readBasis (std::istream& is, size_t basis)
bool ASMu2Dmx::write (std::ostream& os, int basis) const
{
if (basis == -1)
os << *projBasis;
os << *projB;
else
os << *m_basis[basis-1];
@ -207,7 +207,7 @@ bool ASMu2Dmx::generateFEMTopology ()
if (tensorPrjBas)
{
projBasis.reset(createLR(*tensorPrjBas));
projB.reset(createLR(*tensorPrjBas));
delete tensorPrjBas;
tensorPrjBas = nullptr;
}
@ -224,33 +224,33 @@ bool ASMu2Dmx::generateFEMTopology ()
ASMmxBase::Type == ASMmxBase::DIV_COMPATIBLE ||
ASMmxBase::Type == ASMmxBase::SUBGRID) {
Go::SplineSurface* otherBasis = nullptr;
if (!projBasis)
if (!projB)
otherBasis = ASMmxBase::raiseBasis(tensorspline);
if (ASMmxBase::Type == ASMmxBase::SUBGRID) {
refBasis.reset(createLR(*otherBasis));
if (!projBasis)
projBasis = m_basis.front();
altProjBasis = refBasis;
if (!projB)
projB = m_basis.front();
projB2 = refBasis;
}
else {
if (!projBasis)
projBasis.reset(createLR(*otherBasis));
refBasis = projBasis;
if (!projB)
projB.reset(createLR(*otherBasis));
refBasis = std::static_pointer_cast<LR::LRSplineSurface>(projB);
}
delete otherBasis;
}
else {
if (!projBasis)
projBasis = m_basis[2-ASMmxBase::geoBasis];
refBasis = projBasis;
if (!projB)
projB = m_basis[2-ASMmxBase::geoBasis];
refBasis = std::static_pointer_cast<LR::LRSplineSurface>(projB);
}
is_rational = tensorspline->rational();
delete tensorspline;
tensorspline = nullptr;
}
projBasis->generateIDs();
projB->generateIDs();
refBasis->generateIDs();
lrspline = m_basis[geoBasis-1];
@ -988,15 +988,15 @@ bool ASMu2Dmx::refine (const LR::RefineData& prm, Vectors& sol)
// Uniformly refine to find basis 1
if (ASMmxBase::Type == ASMmxBase::SUBGRID) {
m_basis[0].reset(refBasis->copy());
projBasis = m_basis.front();
projB = m_basis.front();
size_t nFunc = refBasis->nBasisFunctions();
IntVec elems(nFunc);
std::iota(elems.begin(),elems.end(),0);
m_basis[0]->refineBasisFunction(elems);
}
if (altProjBasis)
altProjBasis->generateIDs();
if (projB2)
projB2->generateIDs();
size_t len = 0;
for (size_t j = 0; j< m_basis.size(); ++j) {
@ -1022,8 +1022,8 @@ bool ASMu2Dmx::refine (const LR::RefineData& prm, Vectors& sol)
std::cout << it->nBasisFunctions() <<" ";
std::cout <<"nodes."<< std::endl;
std::cout << "Projection basis: "
<< projBasis->nElements() << " elements "
<< projBasis->nBasisFunctions() << " nodes" << std::endl;
<< projB->nElements() << " elements "
<< projB->nBasisFunctions() << " nodes" << std::endl;
std::cout << "Refinement basis: "
<< refBasis->nElements() << " elements "
<< refBasis->nBasisFunctions() << " nodes" << std::endl;
@ -1078,9 +1078,9 @@ void ASMu2Dmx::generateThreadGroups (const Integrand& integrand, bool silence,
secConstraint = {this->getBasis(1), this->getBasis(2)};
LR::generateThreadGroups(threadGroups,threadBasis,secConstraint);
LR::generateThreadGroups(projThreadGroups,projBasis.get());
if (altProjBasis)
LR::generateThreadGroups(altProjThreadGroups,altProjBasis.get());
LR::generateThreadGroups(projThreadGroups,projB.get());
if (projB2)
LR::generateThreadGroups(altProjThreadGroups,projB2.get());
std::vector<const LR::LRSpline*> bases;
for (const SplinePtr& basis : m_basis)
@ -1192,7 +1192,7 @@ void ASMu2Dmx::storeMesh (const std::string& fName, int fType) const
writeBasis(patch, btag);
++btag.back();
}
writeBasis(projBasis, "proj");
writeBasis(std::static_pointer_cast<LR::LRSplineSurface>(projB), "proj");
writeBasis(refBasis, "ref");
}
@ -1214,9 +1214,9 @@ void ASMu2Dmx::copyRefinement (LR::LRSplineSurface* basis,
void ASMu2Dmx::swapProjectionBasis ()
{
if (altProjBasis) {
if (projB2) {
ASMmxBase::geoBasis = ASMmxBase::geoBasis == 1 ? 2 : 1;
std::swap(projBasis, altProjBasis);
std::swap(projB, projB2);
std::swap(projThreadGroups, altProjThreadGroups);
lrspline = m_basis[ASMmxBase::geoBasis-1];
geomB = lrspline;

1
src/ASM/LR/ASMu2Dmx.h
View File

@ -270,7 +270,6 @@ private:
std::vector<SplinePtr> m_basis; //!< All bases
LR::LRSplineSurface* threadBasis; //!< Basis for thread groups
SplinePtr refBasis; //!< Basis to refine based on
SplinePtr altProjBasis; //!< Alternative projection basis
ThreadGroups altProjThreadGroups; //!< Element groups for multi-threaded assembly - alternative projection basis
};

18
src/ASM/LR/ASMu2Drecovery.C
View File

@ -108,8 +108,8 @@ bool ASMu2D::assembleL2matrices (SparseMatrix& A, StdVector& B,
{
size_t nnod = this->getNoProjectionNodes();
const int p1 = projBasis->order(0);
const int p2 = projBasis->order(1);
const int p1 = projB->order(0);
const int p2 = projB->order(1);
const int pm = p1 > p2 ? p1 : p2;
// Get Gaussian quadrature points
@ -121,12 +121,12 @@ bool ASMu2D::assembleL2matrices (SparseMatrix& A, StdVector& B,
if (!xg || !yg) return false;
if (continuous && !wg) return false;
bool singleBasis = (this->getNoBasis() == 1 && projBasis == lrspline);
bool singleBasis = (this->getNoBasis() == 1 && projB == lrspline);
IntMat lmnpc;
const IntMat& gmnpc = singleBasis ? MNPC : lmnpc;
if (!singleBasis) {
lmnpc.resize(projBasis->nElements());
for (const LR::Element* elm : projBasis->getAllElements()) {
lmnpc.resize(projB->nElements());
for (const LR::Element* elm : projB->getAllElements()) {
lmnpc[elm->getId()].reserve(elm->nBasisFunctions());
for (const LR::Basisfunction* f : elm->support())
lmnpc[elm->getId()].push_back(f->getId());
@ -147,7 +147,7 @@ bool ASMu2D::assembleL2matrices (SparseMatrix& A, StdVector& B,
Go::BasisPtsSf spl0;
Go::BasisDerivsSf spl1, spl2;
int ielp = group[t][e];
const LR::Element* elm = projBasis->getElement(ielp);
const LR::Element* elm = projB->getElement(ielp);
int iel = lrspline->getElementContaining(elm->midpoint())+1;
if (continuous)
@ -189,14 +189,16 @@ bool ASMu2D::assembleL2matrices (SparseMatrix& A, StdVector& B,
{
if (continuous)
{
this->computeBasis(gpar[0][i],gpar[1][j],spl1,ielp,projBasis.get());
this->computeBasis(gpar[0][i],gpar[1][j],spl1,ielp,
static_cast<const LR::LRSplineSurface*>(projB.get()));
SplineUtils::extractBasis(spl1,phi,dNdu);
this->computeBasis(gpar[0][i],gpar[1][j],spl2,iel-1);
SplineUtils::extractBasis(spl2,phi2,dNdu);
}
else
{
this->computeBasis(gpar[0][i],gpar[1][j],spl0,ielp,projBasis.get());
this->computeBasis(gpar[0][i],gpar[1][j],spl0,ielp,
static_cast<const LR::LRSplineSurface*>(projB.get()));
phi = spl0.basisValues;
}

27
src/ASM/LR/ASMu3D.C
View File

@ -55,7 +55,6 @@ ASMu3D::ASMu3D (const ASMu3D& patch, unsigned char n_f)
{
vMin = 0.0;
tensorspline = tensorPrjBas = nullptr;
projBasis = patch.projBasis;
// Need to set nnod here,
// as hasXNodes might be invoked before the FE data is generated
@ -230,8 +229,8 @@ bool ASMu3D::createProjectionBasis (bool init)
tensorPrjBas = tensorspline->clone();
std::swap(tensorspline,tensorPrjBas);
std::swap(lrspline,projBasis);
geomB = lrspline;
std::swap(geomB,projB);
lrspline = std::static_pointer_cast<LR::LRSplineVolume>(geomB);
return true;
}
@ -255,12 +254,12 @@ bool ASMu3D::generateFEMTopology ()
if (tensorPrjBas)
{
projBasis.reset(new LR::LRSplineVolume(tensorPrjBas));
projB.reset(new LR::LRSplineVolume(tensorPrjBas));
delete tensorPrjBas;
tensorPrjBas = nullptr;
}
else if (!projBasis)
projBasis = lrspline;
else if (!projB)
projB = lrspline;
if (!lrspline) return false;
@ -1889,8 +1888,8 @@ void ASMu3D::generateThreadGroups (const Integrand& integrand, bool silence,
bool ignoreGlobalLM)
{
LR::generateThreadGroups(threadGroups, this->getBasis(1));
if (projBasis != lrspline)
LR::generateThreadGroups(projThreadGroups, projBasis.get());
if (this->separateProjectionBasis())
LR::generateThreadGroups(projThreadGroups, projB.get());
if (silence || threadGroups[0].size() < 2) return;
IFEM::cout <<"\nMultiple threads are utilized during element assembly.";
@ -1970,20 +1969,20 @@ size_t ASMu3D::getNoNodes (int basis) const
size_t ASMu3D::getNoProjectionNodes () const
{
return projBasis->nBasisFunctions();
return projB->nBasisFunctions();
}
bool ASMu3D::separateProjectionBasis () const
{
return projBasis.get() != this->getBasis(1);
return projB.get() != this->getBasis(1);
}
Field* ASMu3D::getProjectedField (const Vector& coefs) const
{
if (coefs.size() == this->getNoProjectionNodes())
return new LRSplineField3D(projBasis.get(),coefs);
return new LRSplineField3D(static_cast<const LR::LRSplineVolume*>(projB.get()),coefs);
std::cerr <<" *** ASMu3D::getProjectedFields: Non-matching coefficent array,"
<<" size="<< coefs.size() <<" nnod="<< this->getNoProjectionNodes()
@ -1994,12 +1993,12 @@ Field* ASMu3D::getProjectedField (const Vector& coefs) const
Fields* ASMu3D::getProjectedFields (const Vector& coefs, size_t) const
{
if (projBasis.get() == this->getBasis(1))
if (!this->separateProjectionBasis())
return nullptr;
size_t ncmp = coefs.size() / this->getNoProjectionNodes();
if (ncmp*this->getNoProjectionNodes() == coefs.size())
return new LRSplineFields3D(projBasis.get(),coefs,ncmp);
return new LRSplineFields3D(static_cast<const LR::LRSplineVolume*>(projB.get()),coefs,ncmp);
std::cerr <<" *** ASMu3D::getProjectedFields: Non-matching coefficent array,"
<<" size="<< coefs.size() <<" nnod="<< this->getNoProjectionNodes()
@ -2317,7 +2316,7 @@ bool ASMu3D::refine (const LR::RefineData& prm, Vectors& sol)
// TODO: check this
for (const LR::MeshRectangle* rect : lrspline->getAllMeshRectangles())
projBasis->insert_line(rect->copy());
std::static_pointer_cast<LR::LRSplineVolume>(projB)->insert_line(rect->copy());
return true;
}

1
src/ASM/LR/ASMu3D.h
View File

@ -693,7 +693,6 @@ public:
protected:
std::shared_ptr<LR::LRSplineVolume> lrspline; //!< Pointer to the LR-spline volume object
std::shared_ptr<LR::LRSplineVolume> projBasis; //!< Basis to project onto
Go::SplineVolume* tensorspline; //!< Pointer to original tensor spline object
Go::SplineVolume* tensorPrjBas; //!< Pointer to tensor spline projection base

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

@ -98,7 +98,7 @@ bool ASMu3Dmx::readBasis (std::istream& is, size_t basis)
bool ASMu3Dmx::write (std::ostream& os, int basis) const
{
if (basis == -1)
os << *projBasis;
os << *projB;
else
os << *m_basis[basis-1];
@ -204,9 +204,8 @@ bool ASMu3Dmx::generateFEMTopology ()
if (m_basis.empty()) {
VolumeVec vvec = ASMmxBase::establishBases(tensorspline, ASMmxBase::Type);
m_basis.resize(vvec.size());
for (size_t b = 0; b < vvec.size(); b++)
m_basis[b].reset(new LR::LRSplineVolume(vvec[b].get()));
m_basis.push_back(std::make_shared<LR::LRSplineVolume>(vvec[b].get()));
// we need to project on something that is not one of our bases
if (ASMmxBase::Type == ASMmxBase::REDUCED_CONT_RAISE_BASIS1 ||
@ -219,29 +218,29 @@ bool ASMu3Dmx::generateFEMTopology ()
ASMmxBase::FULL_CONT_RAISE_BASIS1).front();
geoBasis = geoB;
if (ASMmxBase::Type == ASMmxBase::SUBGRID) {
refBasis.reset(new LR::LRSplineVolume(otherBasis.get()));
projBasis = m_basis.front();
refBasis = std::make_shared<LR::LRSplineVolume>(otherBasis.get());
projB = m_basis.front();
refBasis->generateIDs();
altProjBasis = refBasis;
projB2 = refBasis;
} else {
projBasis.reset(new LR::LRSplineVolume(otherBasis.get()));
refBasis = projBasis;
projB = std::make_shared<LR::LRSplineVolume>(otherBasis.get());
refBasis = std::static_pointer_cast<LR::LRSplineVolume>(projB);
}
} else {
if (!projBasis)
projBasis = m_basis[2-ASMmxBase::geoBasis];
refBasis = projBasis;
if (!projB)
projB = m_basis[2-ASMmxBase::geoBasis];
refBasis = std::static_pointer_cast<LR::LRSplineVolume>(projB);
}
delete tensorspline;
tensorspline = nullptr;
}
lrspline = m_basis[geoBasis-1];
projBasis->generateIDs();
projBasis->getElementContaining(projBasis->getElement(0)->midpoint()); // to force cache generation
if (altProjBasis) {
altProjBasis->generateIDs();
altProjBasis->getElementContaining(projBasis->getElement(0)->midpoint()); // to force cache generation
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;
@ -866,15 +865,15 @@ bool ASMu3Dmx::refine (const LR::RefineData& prm, Vectors& sol)
// Uniformly refine to find basis 1
if (ASMmxBase::Type == ASMmxBase::SUBGRID) {
m_basis[0].reset(refBasis->copy());
projBasis = m_basis.front();
projB = m_basis.front();
size_t nFunc = refBasis->nBasisFunctions();
IntVec elems(nFunc);
std::iota(elems.begin(),elems.end(),0);
m_basis[0]->refineBasisFunction(elems);
}
if (altProjBasis)
altProjBasis->generateIDs();
if (projB2)
projB2->generateIDs();
size_t len = 0;
for (size_t j = 0; j< m_basis.size(); ++j) {
@ -900,8 +899,8 @@ bool ASMu3Dmx::refine (const LR::RefineData& prm, Vectors& sol)
std::cout << it->nBasisFunctions() <<" ";
std::cout <<"nodes."<< std::endl;
std::cout << "Projection basis: "
<< projBasis->nElements() << " elements "
<< projBasis->nBasisFunctions() << " nodes" << std::endl;
<< projB->nElements() << " elements "
<< projB->nBasisFunctions() << " nodes" << std::endl;
std::cout << "Refinement basis: "
<< refBasis->nElements() << " elements "
<< refBasis->nBasisFunctions() << " nodes" << std::endl;
@ -956,9 +955,9 @@ void ASMu3Dmx::generateThreadGroups (const Integrand& integrand, bool silence,
secConstraint = {this->getBasis(1),this->getBasis(2),this->getBasis(3)};
LR::generateThreadGroups(threadGroups,threadBasis,secConstraint);
LR::generateThreadGroups(projThreadGroups,projBasis.get());
if (altProjBasis)
LR::generateThreadGroups(altProjThreadGroups,altProjBasis.get());
LR::generateThreadGroups(projThreadGroups,projB.get());
if (projB2)
LR::generateThreadGroups(altProjThreadGroups,projB2.get());
std::vector<const LR::LRSpline*> bases;
for (const SplinePtr& basis : m_basis)
@ -1023,9 +1022,9 @@ void ASMu3Dmx::copyRefinement (LR::LRSplineVolume* basis,
void ASMu3Dmx::swapProjectionBasis ()
{
if (altProjBasis) {
if (projB2) {
ASMmxBase::geoBasis = ASMmxBase::geoBasis == 1 ? 2 : 1;
std::swap(projBasis, altProjBasis);
std::swap(projB, projB2);
std::swap(projThreadGroups, altProjThreadGroups);
}
}

1
src/ASM/LR/ASMu3Dmx.h
View File

@ -234,7 +234,6 @@ private:
std::vector<SplinePtr> m_basis; //!< All bases
LR::LRSplineVolume* threadBasis; //!< Basis for thread groups
SplinePtr refBasis; //!< Basis to refine based on
SplinePtr altProjBasis; //!< Alternative projection basis
const std::vector<Matrices>& bezierExtractmx; //!< Bezier extraction matrices
std::vector<Matrices> myBezierExtractmx; //!< Bezier extraction matrices

18
src/ASM/LR/ASMu3Drecovery.C
View File

@ -110,9 +110,9 @@ bool ASMu3D::assembleL2matrices (SparseMatrix& A, StdVector& B,
{
size_t nnod = this->getNoProjectionNodes();
const int p1 = projBasis->order(0);
const int p2 = projBasis->order(1);
const int p3 = projBasis->order(2);
const int p1 = projB->order(0);
const int p2 = projB->order(1);
const int p3 = projB->order(2);
const int pm = std::max(std::max(p1,p2),p3);
// Get Gaussian quadrature points
@ -126,12 +126,12 @@ bool ASMu3D::assembleL2matrices (SparseMatrix& A, StdVector& B,
if (!xg || !yg || !zg) return false;
if (continuous && !wg) return false;
bool singleBasis = (this->getNoBasis() == 1 && projBasis == lrspline);
bool singleBasis = (this->getNoBasis() == 1 && projB == lrspline);
IntMat lmnpc;
const IntMat& gmnpc = singleBasis ? MNPC : lmnpc;
if (!singleBasis) {
lmnpc.resize(projBasis->nElements());
for (const LR::Element* elm : projBasis->getAllElements()) {
lmnpc.resize(projB->nElements());
for (const LR::Element* elm : projB->getAllElements()) {
lmnpc[elm->getId()].reserve(elm->nBasisFunctions());
for (const LR::Basisfunction* f : elm->support())
lmnpc[elm->getId()].push_back(f->getId());
@ -152,7 +152,7 @@ bool ASMu3D::assembleL2matrices (SparseMatrix& A, StdVector& B,
Go::BasisPts spl0;
Go::BasisDerivs spl1, spl2;
int ielp = group[t][e];
const LR::Element* elm = projBasis->getElement(ielp);
const LR::Element* elm = projB->getElement(ielp);
int iel = lrspline->getElementContaining(elm->midpoint())+1;
if (continuous)
@ -196,14 +196,14 @@ bool ASMu3D::assembleL2matrices (SparseMatrix& A, StdVector& B,
{
if (continuous)
{
projBasis->computeBasis(gpar[0][i],gpar[1][j],gpar[2][k],spl1,ielp);
static_cast<const LR::LRSplineVolume*>(projB.get())->computeBasis(gpar[0][i],gpar[1][j],gpar[2][k],spl1,ielp);
SplineUtils::extractBasis(spl1,phi,dNdu);
lrspline->computeBasis(gpar[0][i],gpar[1][j],gpar[2][k],spl2,iel-1);
SplineUtils::extractBasis(spl2,phi2,dNdu);
}
else
{
projBasis->computeBasis(gpar[0][i],gpar[1][j],gpar[2][k],spl0,ielp);
static_cast<const LR::LRSplineVolume*>(projB.get())->computeBasis(gpar[0][i],gpar[1][j],gpar[2][k],spl0,ielp);
phi = spl0.basisValues;
}