From e75f9ada09daa87a899313bd543a02ffd597db00 Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Fri, 14 Mar 2025 09:07:51 -0400 Subject: [PATCH] feat(polyMFEMUtils): changed slope constraint to look at all connected elements --- src/poly/utils/private/polyMFEMUtils.cpp | 186 +++++++++++------------ 1 file changed, 86 insertions(+), 100 deletions(-) diff --git a/src/poly/utils/private/polyMFEMUtils.cpp b/src/poly/utils/private/polyMFEMUtils.cpp index 40231f7..b4c4755 100644 --- a/src/poly/utils/private/polyMFEMUtils.cpp +++ b/src/poly/utils/private/polyMFEMUtils.cpp @@ -375,7 +375,9 @@ namespace polyMFEMUtils { } ZeroSlopeNewtonSolver::ZeroSlopeNewtonSolver(double alpha_, std::vector 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 elementVertexIDs; + mesh->GetElementVertices(zeroSlopeElemID, elementVertexIDs); + int centralVertexID = -1; + double smallestDistance = std::numeric_limits::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 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(const_cast(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("Debug", false) ) { + Probe::glVisView(x, *dynamic_cast(const_cast(oper))->FESpace(), "initial guess"); + Probe::getRaySolution(x, *dynamic_cast(const_cast(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(const_cast(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 \ No newline at end of file