tboudreaux/SERiF
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

refactor(serif): fixed typos and updated names to reflect 4DSSE->SERiF

Browse Source
This commit is contained in:
Emily Boudreaux
2025-06-12 09:04:03 -04:00
parent 620f3d495c
commit 2eca802d01
9 changed files with 302 additions and 308 deletions
Download Patch File Download Diff File Expand all files Collapse all files

161
src/polytrope/utils/private/integrators.cpp
View File

@@ -25,9 +25,7 @@
#include "config.h"
#include <string>
namespace serif {
namespace polytrope {
namespace polyMFEMUtils {
namespace serif::polytrope::polyMFEMUtils {
NonlinearPowerIntegrator::NonlinearPowerIntegrator(const double n) :
m_polytropicIndex(n),
m_epsilon(serif::config::Config::getInstance().get<double>("Poly:Solver:Epsilon", 1.0e-8)) {
@@ -48,99 +46,94 @@ namespace polyMFEMUtils {
mfem::ElementTransformation &Trans,
const mfem::Vector &elfun,
mfem::Vector &elvect) {
const mfem::IntegrationRule *ir = &mfem::IntRules.Get(el.GetGeomType(), 2 * el.GetOrder() + 3);
int dof = el.GetDof();
elvect.SetSize(dof);
elvect = 0.0;
mfem::Vector shape(dof);
mfem::Vector physCoord;
for (int iqp = 0; iqp < ir->GetNPoints(); iqp++) {
mfem::IntegrationPoint ip = ir->IntPoint(iqp);
Trans.SetIntPoint(&ip);
const double weight = ip.weight * Trans.Weight();
el.CalcShape(ip, shape);
double u_val = 0.0;
for (int j = 0; j < dof; j++) {
u_val += elfun(j) * shape(j);
}
const mfem::IntegrationRule *ir = &mfem::IntRules.Get(el.GetGeomType(), 2 * el.GetOrder() + 3);
const int dof = el.GetDof();
elvect.SetSize(dof);
elvect = 0.0;
double u_nl;
Trans.Transform(ip, physCoord);
const double r = physCoord.Norml2();
std::ofstream outFile("r.dat", std::ios::app);
outFile << r << '\n';
outFile.close();
if (r > m_regularizationRadius) {
if (u_val < m_epsilon) {
u_nl = fmod(u_val, m_polytropicIndex, m_epsilon);
} else {
u_nl = std::pow(u_val, m_polytropicIndex);
}
mfem::Vector shape(dof);
mfem::Vector physCoord;
for (int iqp = 0; iqp < ir->GetNPoints(); iqp++) {
mfem::IntegrationPoint ip = ir->IntPoint(iqp);
Trans.SetIntPoint(&ip);
const double weight = ip.weight * Trans.Weight();
el.CalcShape(ip, shape);
double u_val = 0.0;
for (int j = 0; j < dof; j++) {
u_val += elfun(j) * shape(j);
}
double u_nl;
Trans.Transform(ip, physCoord);
const double r = physCoord.Norml2();
if (r > m_regularizationRadius) {
if (u_val < m_epsilon) {
u_nl = fmod(u_val, m_polytropicIndex, m_epsilon);
} else {
u_nl = 1.0 - m_polytropicIndex * m_regularizationCoeff * std::pow(r, 2);
}
for (int i = 0; i < dof; i++){
elvect(i) += shape(i) * u_nl * weight;
u_nl = std::pow(u_val, m_polytropicIndex);
}
} else {
u_nl = 1.0 - m_polytropicIndex * m_regularizationCoeff * std::pow(r, 2);
}
for (int i = 0; i < dof; i++){
elvect(i) += shape(i) * u_nl * weight;
}
}
}
void NonlinearPowerIntegrator::AssembleElementGrad (
const mfem::FiniteElement &el,
mfem::ElementTransformation &Trans,
const mfem::Vector &elfun,
mfem::DenseMatrix &elmat) {
const mfem::IntegrationRule *ir = &mfem::IntRules.Get(el.GetGeomType(), 2 * el.GetOrder() + 3);
const int dof = el.GetDof();
elmat.SetSize(dof);
elmat = 0.0;
mfem::Vector shape(dof);
mfem::DenseMatrix dshape(dof, 3);
mfem::DenseMatrix invJ(3, 3);
mfem::Vector physCoord;
for (int iqp = 0; iqp < ir->GetNPoints(); iqp++) {
const mfem::IntegrationPoint &ip = ir->IntPoint(iqp);
Trans.SetIntPoint(&ip);
const double weight = ip.weight * Trans.Weight();
Trans.Transform(ip, physCoord);
double r = physCoord.Norml2();
el.CalcShape(ip, shape);
void NonlinearPowerIntegrator::AssembleElementGrad (
const mfem::FiniteElement &el,
mfem::ElementTransformation &Trans,
const mfem::Vector &elfun,
mfem::DenseMatrix &elmat) {
double u_val = 0.0;
for (int j = 0; j < dof; j++) {
u_val += elfun(j) * shape(j);
}
const mfem::IntegrationRule *ir = &mfem::IntRules.Get(el.GetGeomType(), 2 * el.GetOrder() + 3);
const int dof = el.GetDof();
elmat.SetSize(dof);
elmat = 0.0;
mfem::Vector shape(dof);
mfem::DenseMatrix dshape(dof, 3);
mfem::DenseMatrix invJ(3, 3);
mfem::Vector physCoord;
double d_u_nl;
if (r > m_regularizationRadius) {
if (u_val < m_epsilon) {
d_u_nl = dfmod(m_epsilon, m_polytropicIndex);
} else {
d_u_nl = m_polytropicIndex * std::pow(u_val, m_polytropicIndex - 1.0);
}
} else {
d_u_nl = 0.0;
}
for (int iqp = 0; iqp < ir->GetNPoints(); iqp++) {
const mfem::IntegrationPoint &ip = ir->IntPoint(iqp);
Trans.SetIntPoint(&ip);
const double weight = ip.weight * Trans.Weight();
Trans.Transform(ip, physCoord);
double r = physCoord.Norml2();
for (int i = 0; i < dof; i++) {
for (int j = 0; j < dof; j++) {
elmat(i, j) += shape(i) * d_u_nl * shape(j) * weight;
}
}
el.CalcShape(ip, shape);
double u_val = 0.0;
for (int j = 0; j < dof; j++) {
u_val += elfun(j) * shape(j);
}
double d_u_nl;
if (r > m_regularizationRadius) {
if (u_val < m_epsilon) {
d_u_nl = dfmod(m_epsilon, m_polytropicIndex);
} else {
d_u_nl = m_polytropicIndex * std::pow(u_val, m_polytropicIndex - 1.0);
}
} else {
d_u_nl = 0.0;
}
for (int i = 0; i < dof; i++) {
for (int j = 0; j < dof; j++) {
elmat(i, j) += shape(i) * d_u_nl * shape(j) * weight;
}
}
}
}
} // namespace polyMFEMUtils
} // namespace polytrope
} // namespace serif
}