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

feat(polyMFEMUtils): changed slope constraint to look at all connected elements

Browse Source
This commit is contained in:
Emily Boudreaux
2025-03-14 09:07:51 -04:00
parent 7ea20369b2
commit e75f9ada09
1 changed files with 86 additions and 100 deletions
Download Patch File Download Diff File Expand all files Collapse all files

186
src/poly/utils/private/polyMFEMUtils.cpp
View File

@@ -375,7 +375,9 @@ namespace polyMFEMUtils {
}
ZeroSlopeNewtonSolver::ZeroSlopeNewtonSolver(double alpha_, std::vector<double> zeroSlopeCoordinate_)
: alpha(alpha_), zeroSlopeCoordinate(zeroSlopeCoordinate_) {}
: alpha(alpha_), zeroSlopeCoordinate(zeroSlopeCoordinate_) {
zeroIP.Set3w(zeroIPReferenceCoord);
}
ZeroSlopeNewtonSolver::~ZeroSlopeNewtonSolver() {}
@@ -421,70 +423,34 @@ namespace polyMFEMUtils {
mesh->FindPoints(zeroSlopeCoordinateMatrix, elementsIDs, ips);
zeroSlopeElemID = elementsIDs[0];
zeroSlopeIP = ips[0];
mfem::Array<int> elementVertexIDs;
mesh->GetElementVertices(zeroSlopeElemID, elementVertexIDs);
int centralVertexID = -1;
double smallestDistance = std::numeric_limits<double>::max();
for (int i = 0; i < elementVertexIDs.Size(); i++) {
const double *vertex = mesh->GetVertex(elementVertexIDs[i]);
double distance = 0.0;
for (int dimID = 0; dimID < mesh->SpaceDimension(); dimID++) {
distance += std::pow(vertex[dimID] - zeroSlopeCoordinate[dimID], 2);
}
distance = std::sqrt(distance);
if (distance < smallestDistance) {
smallestDistance = distance;
centralVertexID = elementVertexIDs[i];
}
LOG_INFO(logger, "Getting element dofs for zero slope constraint...");
fes->GetElementDofs(zeroSlopeElemID, zeroSlopeDofs);
LOG_INFO(logger, "Getting element dofs for zero slope constraint...done");
LOG_INFO(logger, "Building location of zero slope constraint...done");
}
mfem::Table* vertexToElementTable = mesh->GetVertexToElementTable();
vertexToElementTable->GetRow(centralVertexID, zeroSlopeConnectedElements);
mfem::Array<int> tempZeroSlopeDofs;
for (auto elemID: zeroSlopeConnectedElements) {
fes->GetElementDofs(elemID, tempZeroSlopeDofs);
zeroSlopeDofs.push_back(tempZeroSlopeDofs);
}
}
// void ZeroSlopeNewtonSolver::ProcessNewState(const mfem::Vector &x) const {
// LOG_INFO(logger, "Processing new state for zero slope constraint...");
// if (zeroSlopeElemID < 0) {
// LOG_ERROR(logger, "ZeroSlopeNewtonSolver::ProcessNewState: zero slope element ID is not set");
// MFEM_ABORT("ZeroSlopeNewtonSolver::ProcessNewState: zero slope element ID is not set");
// }
// mfem::NonlinearForm *nlf = dynamic_cast<mfem::NonlinearForm*>(const_cast<mfem::Operator*>(oper));
// if (!nlf) {
// LOG_ERROR(logger, "ZeroSlopeNewtonSolver::ProcessNewState: input operator is not a NonlinearForm");
// MFEM_ABORT("ZeroSlopeNewtonSolver::ProcessNewState: input operator is not a NonlinearForm");
// }
// mfem::FiniteElementSpace *fes = nlf->FESpace();
// if (!fes) {
// LOG_ERROR(logger, "ZeroSlopeNewtonSolver::ProcessNewState: input operator does not have a finite element space");
// MFEM_ABORT("ZeroSlopeNewtonSolver::ProcessNewState: input operator does not have a finite element space");
// }
// mfem::Mesh *mesh = fes->GetMesh();
// if (!mesh) {
// LOG_ERROR(logger, "ZeroSlopeNewtonSolver::ProcessNewState: input operator does not have a mesh");
// MFEM_ABORT("ZeroSlopeNewtonSolver::ProcessNewState: input operator does not have a mesh");
// }
// mfem::ElementTransformation *T = mesh->GetElementTransformation(zeroSlopeElemID);
// if (!T) {
// LOG_ERROR(logger, "ZeroSlopeNewtonSolver::ProcessNewState: element transformation is not found");
// MFEM_ABORT("ZeroSlopeNewtonSolver::ProcessNewState: element transformation is not found");
// }
// mfem::Vector grad_u(3);
// mfem::GridFunction u_gf(fes);
// DEPRECATION_WARNING_OFF
// u_gf.SetData(x.GetData());
// DEPRECATION_WARNING_ON
// T->SetIntPoint(&zeroSlopeIP);
// u_gf.GetGradient(*T, grad_u);
// int dof;
// LOG_DEBUG(logger, "Adjusting the residual to enforce the zero slope constraint by {:0.4E}...", -alpha*grad_u[0]);
// double rNorm = r.Norml2();
// LOG_INFO(logger, "||r_B|| = {:0.4E}", rNorm);
// for (int i = 0; i < zeroSlopeDofs.Size(); i++) {
// dof = zeroSlopeDofs[i];
// r[dof] -= alpha * grad_u[0];
// r[dof] -= alpha * grad_u[1];
// r[dof] -= alpha * grad_u[2];
// }
// rNorm = r.Norml2();
// LOG_INFO(logger, "||r_A|| = {:0.4E}", rNorm);
// // This still is not working; however, I think I am close. I also need to modify the jacobain.
// }
void ZeroSlopeNewtonSolver::Mult(const mfem::Vector &b, mfem::Vector &x) const {
using namespace mfem;
using namespace std;
@@ -501,14 +467,22 @@ namespace polyMFEMUtils {
x = 0.0;
}
if ( config.get<bool>("Debug", false) ) {
Probe::glVisView(x, *dynamic_cast<mfem::NonlinearForm*>(const_cast<mfem::Operator*>(oper))->FESpace(), "initial guess");
Probe::getRaySolution(x, *dynamic_cast<mfem::NonlinearForm*>(const_cast<mfem::Operator*>(oper))->FESpace(), {0.0, 0.0}, 100, "output/initial_guess.csv");
}
ProcessNewState(x);
DEPRECATION_WARNING_OFF
u_gf->SetData(x.GetData());
DEPRECATION_WARNING_ON
oper->Mult(x, r);
if (have_b)
{
r -= b;
}
// ComputeConstrainedResidual(x, r);
ComputeConstrainedResidual(x, r);
norm0 = norm = initial_norm = Norm(r);
if (print_options.first_and_last && !print_options.iterations)
@@ -579,15 +553,19 @@ namespace polyMFEMUtils {
ProcessNewState(x);
// Probe::glVisView(x, *dynamic_cast<mfem::NonlinearForm*>(const_cast<mfem::Operator*>(oper))->FESpace(), "solution " + it);
oper->Mult(x, r);
if (have_b)
{
r -= b;
}
// ComputeConstrainedResidual(x, r);
ComputeConstrainedResidual(x, r);
norm = Norm(r);
}
LOG_INFO(logger, "Final Computation of residual...");
ComputeConstrainedResidual(x, r);
final_iter = it;
final_norm = norm;
@@ -623,25 +601,29 @@ namespace polyMFEMUtils {
MFEM_ABORT("ZeroSlopeNewtonSolver::ProcessNewState: input operator does not have a mesh");
}
mfem::ElementTransformation *T = mesh->GetElementTransformation(zeroSlopeElemID);
if (!T) {
LOG_ERROR(logger, "ZeroSlopeNewtonSolver::ProcessNewState: element transformation is not found");
MFEM_ABORT("ZeroSlopeNewtonSolver::ProcessNewState: element transformation is not found");
}
int i = 0;
double grad_x_avg=0.0, grad_y_avg=0.0, grad_z_avg=0.0;
for (auto elemID : zeroSlopeConnectedElements) {
mfem::ElementTransformation *T = mesh->GetElementTransformation(elemID);
if (!T) {
LOG_ERROR(logger, "ZeroSlopeNewtonSolver::ProcessNewState: element transformation is not found");
MFEM_ABORT("ZeroSlopeNewtonSolver::ProcessNewState: element transformation is not found");
}
DEPRECATION_WARNING_OFF
u_gf->SetData(x.GetData());
DEPRECATION_WARNING_ON
T->SetIntPoint(&zeroIP);
mfem::Vector grad_u(3); // TODO make this a unique pointer so it can be dimensionally adaptive
u_gf->GetGradient(*T, grad_u);
grad_x_avg += grad_u[0];
grad_y_avg += grad_u[1];
grad_z_avg += grad_u[2];
T->SetIntPoint(&zeroSlopeIP);
mfem::Vector grad_u(3); // TODO make this a unique pointer so it can be dimensionally adaptive
u_gf->GetGradient(*T, grad_u);
for (int i = 0; i < zeroSlopeDofs.Size(); i++) {
int dof = zeroSlopeDofs[i];
residual[dof] -= alpha * grad_u[0];
residual[dof] -= alpha * grad_u[1];
residual[dof] -= alpha * grad_u[2];
for (int j = 0; j < zeroSlopeDofs[i].Size(); j++) {
int dof = zeroSlopeDofs[i][j];
residual[dof] -= alpha * grad_u[0];
residual[dof] -= alpha * grad_u[1];
residual[dof] -= alpha * grad_u[2];
}
i++;
}
}
@@ -665,31 +647,35 @@ namespace polyMFEMUtils {
MFEM_ABORT("ZeroSlopeNewtonSolver::ProcessNewState: input operator does not have a mesh");
}
mfem::ElementTransformation *T = mesh->GetElementTransformation(zeroSlopeElemID);
if (!T) {
LOG_ERROR(logger, "ZeroSlopeNewtonSolver::ProcessNewState: element transformation is not found");
MFEM_ABORT("ZeroSlopeNewtonSolver::ProcessNewState: element transformation is not found");
}
const mfem::FiniteElement* fe = fes->GetFE(zeroSlopeElemID); // Get FE *once*.
mfem::DenseMatrix dshape; // For shape function derivatives.
dshape.SetSize(fe->GetDof(), mesh->Dimension());
T->SetIntPoint(&zeroSlopeIP);
fe->CalcDShape(zeroSlopeIP, dshape);
if (!grad) {
LOG_ERROR(logger, "ZeroSlopeNewtonSolver::ComputeConstrainedGradient: Grad is not set");
MFEM_ABORT("ZeroSlopeNewtonSolver::ComputeConstrainedGradient: Grad is not set");
}
// --- Modify Jacobian ---
LOG_INFO(logger, "Adjusting the Jacobian to enforce the zero slope constraint...");
for (int i = 0; i < zeroSlopeDofs.Size(); i++) {
for (int j = 0; j < zeroSlopeDofs.Size(); j++) {
grad->Add(zeroSlopeDofs[i], zeroSlopeDofs[j], alpha * dshape(j, 0));
grad->Add(zeroSlopeDofs[i], zeroSlopeDofs[j], alpha * dshape(j, 1));
grad->Add(zeroSlopeDofs[i], zeroSlopeDofs[j], alpha * dshape(j, 2));
}
int dofID = 0;
for (auto elemID : zeroSlopeConnectedElements) {
mfem::ElementTransformation *T = mesh->GetElementTransformation(elemID);
if (!T) {
LOG_ERROR(logger, "ZeroSlopeNewtonSolver::ProcessNewState: element transformation is not found");
MFEM_ABORT("ZeroSlopeNewtonSolver::ProcessNewState: element transformation is not found");
}
const mfem::FiniteElement* fe = fes->GetFE(elemID); // Get FE *once*.
mfem::DenseMatrix dshape; // For shape function derivatives.
dshape.SetSize(fe->GetDof(), mesh->Dimension());
T->SetIntPoint(&zeroIP);
fe->CalcDShape(zeroIP, dshape);
// --- Modify Jacobian ---
for (int i = 0; i < zeroSlopeDofs[dofID].Size(); i++) {
for (int j = 0; j < zeroSlopeDofs[dofID].Size(); j++) {
grad->Add(zeroSlopeDofs[dofID][i], zeroSlopeDofs[dofID][j], alpha * dshape(j, 0));
grad->Add(zeroSlopeDofs[dofID][i], zeroSlopeDofs[dofID][j], alpha * dshape(j, 1));
grad->Add(zeroSlopeDofs[dofID][i], zeroSlopeDofs[dofID][j], alpha * dshape(j, 2));
}
}
}
LOG_INFO(logger, "Adjusting the Jacobian to enforce the zero slope constraint...done");
dofID++;
}
} // namespace polyMFEMUtils