[TODO][DRY] Migrate the elements loops to the new format
With the recent implementation of FEUtilities::LoopOverElements
we can replace many loops, specially in the LinearProblem, to avoid code repetition.
Example:
void LinearProblem::assembleOnlyMassMatrix(PetscMatrix& mass)
{
// Allocation
mass.duplicateFrom(matrix(0),MAT_DO_NOT_COPY_VALUES);
mass.zeroEntries();
// Assemblyloop
ElementType eltType; // Geometric element type in the mesh.
int numPointPerElt = 0; // Number of points per geometric element.
felInt numEltPerLabel = 0; // Number of element for one label and one eltType.
// Use to define a "global" numbering of element in the mesh.
felInt numElement[ GeometricMeshRegion::m_numTypesOfElement ];
for (int ityp=0; ityp<GeometricMeshRegion::m_numTypesOfElement; ityp++ ) {
numElement[static_cast<ElementType>(ityp)] = 0;
}
std::vector<Point*> elemPoint;
std::vector<felInt> elemIdPoint;
felInt ielSupportDof;
// Assembly loop.
// First loop on geometric type.
const std::vector<ElementType>& bagElementTypeDomain = m_meshLocal.bagElementTypeDomain();
for (std::size_t i = 0; i < bagElementTypeDomain.size(); ++i) {
eltType = bagElementTypeDomain[i];
// Resize array.
numPointPerElt = GeometricMeshRegion::m_numPointsPerElt[eltType];
elemPoint.resize(numPointPerElt, nullptr);
elemIdPoint.resize(numPointPerElt, 0);
defineFiniteElement(eltType);
allocateArrayForAssembleMatrixRHS(FlagMatrixRHS::only_matrix);
// Second loop on region of the mesh.
for(auto itRef = m_meshLocal.intRefToBegEndMaps[eltType].begin(); itRef != m_meshLocal.intRefToBegEndMaps[eltType].end(); itRef++) {
numEltPerLabel = itRef->second.second;
// Third loop on element in the region with the type: eltType. ("real" loop on elements).
for ( felInt iel = 0; iel < numEltPerLabel; iel++) {
setElemPoint(eltType, numElement[eltType], elemPoint, elemIdPoint, &ielSupportDof);
m_elementMat[0]->zero();
userComputeMass(elemPoint, elemIdPoint, ielSupportDof);
setValueCustomMatrix(ielSupportDof, ielSupportDof, mass);
numElement[eltType]++;
}
}
// Desallocate array use for assemble with petsc.
desallocateArrayForAssembleMatrixRHS(FlagMatrixRHS::only_matrix);
}
mass.assembly(MAT_FINAL_ASSEMBLY);
}
Becomes:
void LinearProblem::assembleOnlyMassMatrix(PetscMatrix& mass)
{
// Allocation
mass.duplicateFrom(matrix(0),MAT_DO_NOT_COPY_VALUES);
mass.zeroEntries();
// The following it's the equivalent // TODO: Remove after check it works properly
/* Loop function */
auto function = [this,&mass](ElementType& eltType, felInt iel, std::vector<Point*>& elemPoint, std::vector<felInt>& elemIdPoint, felInt& ielSupportDof) {
setElemPoint(eltType, iel, elemPoint, elemIdPoint, &ielSupportDof);
m_elementMat[0]->zero();
userComputeMass(elemPoint, elemIdPoint, ielSupportDof);
setValueCustomMatrix(ielSupportDof, ielSupportDof, mass);
return 0.0;
};
/* Pre-function */
auto functionpre = [this](ElementType& eltType) {
defineFiniteElement(eltType);
allocateArrayForAssembleMatrixRHS(FlagMatrixRHS::only_matrix);
};
/* Post-function */
auto functionpost = [this](ElementType& eltType) {
desallocateArrayForAssembleMatrixRHS(FlagMatrixRHS::only_matrix);
};
FEUtilities::LoopOverElements(m_meshLocal, m_meshLocal.bagElementTypeDomain(), &function, &functionpre, &functionpost);
// Final assembly
mass.assembly(MAT_FINAL_ASSEMBLY);
}
NOTE: The pre and post functions cann be left blank, as it default argument is nullptr
.