Skip to content

Adoption of new domain iterators in G+Smo #40

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 2 commits into from
Feb 17, 2025
Merged

Adoption of new domain iterators in G+Smo #40

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
31d8eb8
adoption of new domain iterators
hverhelst Feb 14, 2025
489c6ac
fix test
hverhelst Feb 17, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
78 changes: 48 additions & 30 deletions gsGeoUtils.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -214,16 +214,20 @@ index_t checkGeometry(gsMultiPatch<T> const & domain)
numNodes.at(i) = domain.basis(p).degree(i)+1;
gsQuadRule<T> quRule = gsQuadrature::get<T>(gsQuadrature::GaussLegendre,numNodes);

typename gsBasis<T>::domainIter domIt = domain.basis(p).makeDomainIterator(boundary::none);
#ifdef _OPENMP
// Initialize domain element iterator
typename gsBasis<T>::domainIter domIt = domain.basis(p).domain()->beginAll();
typename gsBasis<T>::domainIter domItEnd = domain.basis(p).domain()->endAll();

# ifdef _OPENMP
const int tid = omp_get_thread_num();
const int nt = omp_get_num_threads();
for ( domIt->next(tid); domIt->good() && continueIt; domIt->next(nt) )
#else
for (; domIt->good() && continueIt; domIt->next() )
#endif
domIt += tid;
for ( ; domIt<domItEnd; domIt+=nt )
# else
for (; domIt<domItEnd && continueIt; ++domIt )
# endif
{
genSamplingPoints(domIt->lowerCorner(),domIt->upperCorner(),quRule,points);
genSamplingPoints(domIt.lowerCorner(),domIt.upperCorner(),quRule,points);
md.points = points;
domain.patch(p).computeMap(md);
for (index_t q = 0; q < points.cols() && continueIt; ++q)
Expand Down Expand Up @@ -260,16 +264,20 @@ index_t checkDisplacement(gsMultiPatch<T> const & domain, gsMultiPatch<T> const
numNodes.at(i) = displacement.basis(p).degree(i)+1;
gsQuadRule<T> quRule = gsQuadrature::get<T>(gsQuadrature::GaussLegendre,numNodes);

typename gsBasis<T>::domainIter domIt = displacement.basis(p).makeDomainIterator(boundary::none);
#ifdef _OPENMP
// Initialize domain element iterator
typename gsBasis<T>::domainIter domIt = displacement.basis(p).domain()->beginAll();
typename gsBasis<T>::domainIter domItEnd = displacement.basis(p).domain()->endAll();

# ifdef _OPENMP
const int tid = omp_get_thread_num();
const int nt = omp_get_num_threads();
for ( domIt->next(tid); domIt->good() && continueIt; domIt->next(nt) )
#else
for (; domIt->good() && continueIt; domIt->next() )
#endif
domIt += tid;
for ( ; domIt<domItEnd; domIt+=nt )
# else
for (; domIt<domItEnd && continueIt; ++domIt )
# endif
{
genSamplingPoints(domIt->lowerCorner(),domIt->upperCorner(),quRule,points);
genSamplingPoints(domIt.lowerCorner(),domIt.upperCorner(),quRule,points);
mdG.points = points;
mdU.points = points;
domain.patch(p).computeMap(mdG);
Expand Down Expand Up @@ -308,16 +316,20 @@ T normL2(gsMultiPatch<T> const & domain, gsMultiPatch<T> const & solution)
numNodes.at(i) = solution.basis(p).degree(i)+1;
gsQuadRule<T> quRule = gsQuadrature::get<T>(gsQuadrature::GaussLegendre,numNodes);

typename gsBasis<T>::domainIter domIt = solution.basis(p).makeDomainIterator(boundary::none);
#ifdef _OPENMP
// Initialize domain element iterator
typename gsBasis<T>::domainIter domIt = solution.basis(p).domain()->beginAll();
typename gsBasis<T>::domainIter domItEnd = solution.basis(p).domain()->endAll();

# ifdef _OPENMP
const int tid = omp_get_thread_num();
const int nt = omp_get_num_threads();
for ( domIt->next(tid); domIt->good(); domIt->next(nt) )
#else
for (; domIt->good(); domIt->next() )
#endif
domIt += tid;
for ( ; domIt<domItEnd; domIt+=nt )
# else
for (; domIt<domItEnd; ++domIt )
# endif
{
quRule.mapTo( domIt->lowerCorner(), domIt->upperCorner(), quNodes, quWeights );
quRule.mapTo( domIt.lowerCorner(), domIt.upperCorner(), quNodes, quWeights );
mdGeo.points = quNodes;
domain.patch(p).computeMap(mdGeo);
solution.patch(p).eval_into(quNodes,values);
Expand Down Expand Up @@ -347,10 +359,12 @@ T geometryJacRatio(gsMultiPatch<T> const & domain)
numNodes.at(i) = domain.basis(p).degree(i)+1;
gsQuadRule<T> quRule = gsQuadrature::get<T>(gsQuadrature::GaussLegendre,numNodes);

typename gsBasis<T>::domainIter domIt = domain.basis(p).makeDomainIterator(boundary::none);
for (; domIt->good(); domIt->next())
// Initialize domain element iterator
typename gsBasis<T>::domainIter domIt = domain.basis(p).domain()->beginAll();
typename gsBasis<T>::domainIter domItEnd = domain.basis(p).domain()->endAll();
for (; domIt<domItEnd; ++domIt )
{
genSamplingPoints(domIt->lowerCorner(),domIt->upperCorner(),quRule,points);
genSamplingPoints(domIt.lowerCorner(),domIt.upperCorner(),quRule,points);
md.points = points;
domain.patch(p).computeMap(md);

Expand Down Expand Up @@ -390,10 +404,12 @@ T displacementJacRatio(const gsMultiPatch<T> & domain,const gsMultiPatch<T> & di
numNodes.at(i) = displacement.basis(p).degree(i)+1;
gsQuadRule<T> quRule = gsQuadrature::get<T>(gsQuadrature::GaussLegendre,numNodes);

typename gsBasis<T>::domainIter domIt = domain.basis(p).makeDomainIterator(boundary::none);
for (; domIt->good(); domIt->next())
// Initialize domain element iterator
typename gsBasis<T>::domainIter domIt = domain.basis(p).domain()->beginAll();
typename gsBasis<T>::domainIter domItEnd = domain.basis(p).domain()->endAll();
for (; domIt<domItEnd; ++domIt )
{
genSamplingPoints(domIt->lowerCorner(),domIt->upperCorner(),quRule,points);
genSamplingPoints(domIt.lowerCorner(),domIt.upperCorner(),quRule,points);
mdG.points = points;
mdU.points = points;
domain.patch(p).computeMap(mdG);
Expand Down Expand Up @@ -520,10 +536,12 @@ T curveLength(const gsGeometry<T> & geo)
gsMapData<T> md;
md.flags = NEED_DERIV;

typename gsBasis<>::domainIter domIt = geo.basis().makeDomainIterator(boundary::none);
for (; domIt->good(); domIt->next())
// Initialize domain element iterator
typename gsBasis<T>::domainIter domIt = geo.basis().domain()->beginAll();
typename gsBasis<T>::domainIter domItEnd = geo.basis().domain()->endAll();
for (; domIt<domItEnd; ++domIt )
{
quRule.mapTo(domIt->lowerCorner(),domIt->upperCorner(),qPoints,qWeights);
quRule.mapTo(domIt.lowerCorner(),domIt.upperCorner(),qPoints,qWeights);
md.points = qPoints;
geo.computeMap(md);
for (index_t q = 0; q < qWeights.rows(); ++q)
Expand Down
9 changes: 6 additions & 3 deletions gsNsAssembler.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -199,11 +199,14 @@ gsMatrix<T> gsNsAssembler<T>::computeForce(const gsMultiPatch<T> & velocity, con
// setting quadrature rule for the boundary side
gsGaussRule<T> bdQuRule(basis,1.0,1,it->second.direction());
// loop over elements of the side
typename gsBasis<T>::domainIter elem = basis.makeDomainIterator(it->second);
for (; elem->good(); elem->next())

// Initialize domain element iterator
typename gsBasis<T>::domainIter domIt = basis.domain()->beginBdr(it->second);
typename gsBasis<T>::domainIter domItEnd = basis.domain()->endBdr(it->second);
for (; domIt<domItEnd; ++domIt )
{
// mapping quadrature rule to the element
bdQuRule.mapTo(elem->lowerCorner(),elem->upperCorner(),quNodes,quWeights);
bdQuRule.mapTo(domIt.lowerCorner(),domIt.upperCorner(),quNodes,quWeights);
// evaluate geoemtry mapping at the quad points
mdGeo.points = quNodes;
m_pde_ptr->patches().patch(it->first).computeMap(mdGeo);
Expand Down
2 changes: 1 addition & 1 deletion gsVisitorBiharmonic.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -71,7 +71,7 @@ class gsVisitorBiharmonic
pde_ptr->rhs()->eval_into(md.values[0],forceValues);
}

inline void assemble(gsDomainIterator<T> & element,
inline void assemble(gsDomainIteratorWrapper<T> & element,
const gsVector<T> & quWeights)
{
GISMO_UNUSED(element);
Expand Down
2 changes: 1 addition & 1 deletion gsVisitorBiharmonicMixed.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -71,7 +71,7 @@ class gsVisitorBiharmonicMixed
pde_ptr->rhs()->eval_into(md.values[0],forceValues);
}

inline void assemble(gsDomainIterator<T> & element,
inline void assemble(gsDomainIteratorWrapper<T> & element,
const gsVector<T> & quWeights)
{
GISMO_UNUSED(element);
Expand Down
2 changes: 1 addition & 1 deletion gsVisitorElPoisson.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -63,7 +63,7 @@ class gsVisitorElPoisson
pde_ptr->rhs()->eval_into(md.values[0],forceValues);
}

inline void assemble(gsDomainIterator<T> & element,
inline void assemble(gsDomainIteratorWrapper<T> & element,
const gsVector<T> & quWeights)
{
GISMO_UNUSED(element);
Expand Down
2 changes: 1 addition & 1 deletion gsVisitorElasticityNeumann.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -68,7 +68,7 @@ class gsVisitorElasticityNeumann
basisRefs.front().eval_into(quNodes,basisValuesDisp);
}

inline void assemble(gsDomainIterator<T> & element,
inline void assemble(gsDomainIteratorWrapper<T> & element,
const gsVector<T> & quWeights)
{
GISMO_UNUSED(element);
Expand Down
2 changes: 1 addition & 1 deletion gsVisitorLinearElasticity.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -83,7 +83,7 @@ class gsVisitorLinearElasticity
pde_ptr->rhs()->eval_into(md.values[0],forceValues);
}

inline void assemble(gsDomainIterator<T> & element,
inline void assemble(gsDomainIteratorWrapper<T> & element,
const gsVector<T> & quWeights)
{
GISMO_UNUSED(element);
Expand Down
2 changes: 1 addition & 1 deletion gsVisitorMass.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -64,7 +64,7 @@ class gsVisitorMassElasticity

}

inline void assemble(gsDomainIterator<T> & element,
inline void assemble(gsDomainIteratorWrapper<T> & element,
const gsVector<T> & quWeights)
{
GISMO_UNUSED(element);
Expand Down
2 changes: 1 addition & 1 deletion gsVisitorMassElasticity.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -64,7 +64,7 @@ class gsVisitorMassElasticity

}

inline void assemble(gsDomainIterator<T> & element,
inline void assemble(gsDomainIteratorWrapper<T> & element,
const gsVector<T> & quWeights)
{
GISMO_UNUSED(element);
Expand Down
2 changes: 1 addition & 1 deletion gsVisitorMixedLinearElasticity.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -79,7 +79,7 @@ class gsVisitorMixedLinearElasticity
pde_ptr->rhs()->eval_into(md.values[0],forceValues);
}

inline void assemble(gsDomainIterator<T> & element,
inline void assemble(gsDomainIteratorWrapper<T> & element,
const gsVector<T> & quWeights)
{
GISMO_UNUSED(element);
Expand Down
2 changes: 1 addition & 1 deletion gsVisitorMixedNonLinearElasticity.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -92,7 +92,7 @@ class gsVisitorMixedNonLinearElasticity
pressure.patch(patch).eval_into(quNodes,pressureValues);
}

inline void assemble(gsDomainIterator<T> & element,
inline void assemble(gsDomainIteratorWrapper<T> & element,
const gsVector<T> & quWeights)
{
GISMO_UNUSED(element);
Expand Down
2 changes: 1 addition & 1 deletion gsVisitorMuscle.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -108,7 +108,7 @@ class gsVisitorMuscle
muscleTendon.piece(patch).eval_into(quNodes,muscleTendonValues);
}

inline void assemble(gsDomainIterator<T> & element,
inline void assemble(gsDomainIteratorWrapper<T> & element,
const gsVector<T> & quWeights)
{
GISMO_UNUSED(element);
Expand Down
8 changes: 4 additions & 4 deletions gsVisitorNavierStokes.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -95,7 +95,7 @@ class gsVisitorNavierStokes
pressure.patch(patch).eval_into(quNodes,pressureValues);
}

inline void assemble(gsDomainIterator<T> & element,
inline void assemble(gsDomainIteratorWrapper<T> & element,
const gsVector<T> & quWeights)
{
GISMO_UNUSED(element);
Expand Down Expand Up @@ -127,7 +127,7 @@ class gsVisitorNavierStokes

protected:

void assembleNewtonUpdate(gsDomainIterator<T> & element,
void assembleNewtonUpdate(gsDomainIteratorWrapper<T> & element,
const gsVector<T> & quWeights)
{
GISMO_UNUSED(element);
Expand Down Expand Up @@ -186,7 +186,7 @@ class gsVisitorNavierStokes
}
}

void assembleNewtonFull(gsDomainIterator<T> & element,
void assembleNewtonFull(gsDomainIteratorWrapper<T> & element,
const gsVector<T> & quWeights)
{
GISMO_UNUSED(element);
Expand Down Expand Up @@ -234,7 +234,7 @@ class gsVisitorNavierStokes
}
}

void assembleOseen(gsDomainIterator<T> & element,
void assembleOseen(gsDomainIteratorWrapper<T> & element,
const gsVector<T> & quWeights)
{
GISMO_UNUSED(element);
Expand Down
2 changes: 1 addition & 1 deletion gsVisitorNonLinearElasticity.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -92,7 +92,7 @@ class gsVisitorNonLinearElasticity
displacement.patch(patch).computeMap(mdDisplacement);
}

inline void assemble(gsDomainIterator<T> & element,
inline void assemble(gsDomainIteratorWrapper<T> & element,
const gsVector<T> & quWeights)
{
GISMO_UNUSED(element);
Expand Down
2 changes: 1 addition & 1 deletion gsVisitorStokes.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -75,7 +75,7 @@ class gsVisitorStokes
pde_ptr->rhs()->eval_into(md.values[0],forceValues);
}

inline void assemble(gsDomainIterator<T> & element,
inline void assemble(gsDomainIteratorWrapper<T> & element,
const gsVector<T> & quWeights)
{
GISMO_UNUSED(element);
Expand Down
2 changes: 1 addition & 1 deletion gsVisitorThermo.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -78,7 +78,7 @@ class gsVisitorThermo
basisRefs.front().eval_into(quNodes,basisValuesDisp);
}

inline void assemble(gsDomainIterator<T> & element,
inline void assemble(gsDomainIteratorWrapper<T> & element,
const gsVector<T> & quWeights)
{
GISMO_UNUSED(element);
Expand Down
2 changes: 1 addition & 1 deletion gsVisitorThermoBoundary.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -80,7 +80,7 @@ class gsVisitorThermoBoundary
basisRefs.front().eval_into(quNodes,basisValuesDisp);
}

inline void assemble(gsDomainIterator<T> & element,
inline void assemble(gsDomainIteratorWrapper<T> & element,
const gsVector<T> & quWeights)
{
GISMO_UNUSED(element);
Expand Down
1 change: 0 additions & 1 deletion gsWriteParaviewMultiPhysics.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -77,7 +77,6 @@ void gsWriteParaviewMultiPhysics(std::map<std::string, const gsField<T>*> fields
std::string const & fn,
unsigned npts, bool mesh, bool ctrlNet)
{
gsDebugVar(fn);
const unsigned numP = fields.begin()->second->patches().nPatches();
gsParaviewCollection collection(fn);
std::string baseName = gsFileManager::getFilename(fn); // file name without a path
Expand Down