From 0a1fb8d8db1601a1a290d58ba9fc6e305bd8d54c Mon Sep 17 00:00:00 2001
From: Emily Boudreaux <thomas@boudreauxmail.com>
Date: Mon, 12 Oct 2026 14:10:52 -0400
Subject: [PATCH] feat(geometry/deformation/backend/mfem): added
 radial_location implementation

---
 .../backend/mfem/radial_locator.cpp           | 126 ++++++++++++++++++
 1 file changed, 126 insertions(+)
 create mode 100644 src/lib/serif/geometry/deformation/backend/mfem/radial_locator.cpp

diff --git a/src/lib/serif/geometry/deformation/backend/mfem/radial_locator.cpp b/src/lib/serif/geometry/deformation/backend/mfem/radial_locator.cpp
new file mode 100644
index 00000000..c3c4128b
--- /dev/null
+++ b/src/lib/serif/geometry/deformation/backend/mfem/radial_locator.cpp
@@ -0,0 +1,126 @@
+#include <algorithm>
+#include <cmath>
+#include <limits>
+#include "mfem/mesh/mesh.hpp"
+#include "serif/geometry/deformation/backend/mfem/detail/radial_locator.hpp"
+
+namespace serif::geometry::deformation::backend::mfem::detail {
+    double logical_radius(const ::mfem::Vector& position) {
+        double radius = 0.0;
+        for (int component = 0; component < position.Size(); ++component) {
+            radius = std::max(radius, std::abs(position[component]));
+        }
+        return radius;
+    }
+
+    std::expected<double, DeformationFailure> boundary_logical_radius(::mfem::Mesh& mesh, const int attribute) {
+        double minimum = std::numeric_limits<double>::infinity();
+        double maximum = 0.0;
+        int samples = 0;
+        ::mfem::Vector position(3);
+
+        for (int boundary = 0; boundary < mesh.GetNBE(); ++boundary) {
+            if (mesh.GetBdrAttribute(boundary) != attribute) {
+                continue;
+            }
+
+            ::mfem::ElementTransformation* transformation = mesh.GetBdrElementTransformation(boundary);
+
+            const int geometry = mesh.GetBdrElementGeometry(boundary);
+
+            if (geometry != ::mfem::Geometry::SQUARE || transformation == nullptr) {
+                return std::unexpected(DeformationFailure{.reason = DeformationFailureReason::incompatible_discretization, .message = "The logical radial foliation requires quadrilateral shell elements.", .element = boundary});
+            }
+
+            const ::mfem::IntegrationRule* vertices = ::mfem::Geometries.GetVertices(geometry);
+
+            for (int vertex = 0; vertex < vertices->GetNPoints(); ++vertex) {
+                transformation->Transform(vertices->IntPoint(vertex), position);
+                for (int component = 0; component < 3; ++component) {
+                    if (!std::isfinite(position[component])) {
+                        return std::unexpected(DeformationFailure{.reason = DeformationFailureReason::non_finite, .message = "A logical boundary vertex is non-finite.", .element = boundary, .point = vertex});
+                    }
+                }
+                const double radius = logical_radius(position);
+
+                minimum = std::min(minimum, radius);
+                maximum = std::max(maximum, radius);
+
+                ++samples;
+            }
+        }
+
+        if (samples == 0 || !(minimum > 0.0) || maximum - minimum > 1.0e-12 * std::max(1.0, maximum)) {
+            return std::unexpected(DeformationFailure{.reason = DeformationFailureReason::invalid_geometry, .message = "A required logical boundary is absent or is not a positive constant L-infinity-radius shell."});
+        }
+
+        return 0.5 * (minimum + maximum);
+    }
+
+    std::expected<LocatedBoundaryPoint, DeformationFailure> locate_boundary_point(::mfem::Mesh& mesh, const int attribute, const ::mfem::Vector& target) {
+        LocatedBoundaryPoint best{.boundary_element = -1, .point = {}};
+
+        double best_residual = std::numeric_limits<double>::infinity();
+        const double tolerance = 2.0e-11 * std::max(1.0, logical_radius(target));
+
+        ::mfem::Vector transformed(3), residual(3);
+
+        for (int boundary = 0; boundary < mesh.GetNBE(); ++boundary) {
+            if (mesh.GetBdrAttribute(boundary) != attribute) {
+                continue;
+            }
+
+            ::mfem::ElementTransformation* transformation = mesh.GetBdrElementTransformation(boundary);
+            ::mfem::IntegrationPoint point = ::mfem::Geometries.GetCenter(::mfem::Geometry::SQUARE);
+
+            bool nonsingular = true;
+
+            for (int iteration = 0; iteration < 8; ++iteration) {
+                transformation->Transform(point, transformed);
+                residual = target;
+                residual -= transformed;
+                if (residual.Norml2() <= tolerance) {
+                    break;
+                }
+
+                transformation->SetIntPoint(&point);
+                const ::mfem::DenseMatrix& jacobian = transformation->Jacobian();
+                double normal00 = 0.0, normal01 = 0.0, normal11 = 0.0;
+                double right0 = 0.0, right1 = 0.0;
+                for (int component = 0; component < 3; ++component) {
+                    normal00 += jacobian(component, 0) * jacobian(component, 0);
+                    normal01 += jacobian(component, 0) * jacobian(component, 1);
+                    normal11 += jacobian(component, 1) * jacobian(component, 1);
+                    right0 += jacobian(component, 0) * residual(component);
+                    right1 += jacobian(component, 1) * residual(component);
+                }
+
+                const double determinant = normal00 * normal11 - normal01 * normal01;
+                if (!(std::abs(determinant) > std::numeric_limits<double>::epsilon())) {
+                    nonsingular = false;
+                    break;
+                }
+                point.x += (normal11 * right0 - normal01 * right1) / determinant;
+                point.y += (normal00 * right1 - normal01 * right0) / determinant;
+            }
+
+            if (!nonsingular || !::mfem::Geometry::CheckPoint(::mfem::Geometry::SQUARE, point, 1.0e-10)) {
+                continue;
+            }
+
+            transformation->Transform(point, transformed);
+            transformed -= target;
+            const double norm = transformed.Norml2();
+            if (norm <= tolerance && norm < best_residual) {
+                best = {.boundary_element = boundary, .point = point};
+                best_residual = norm;
+            }
+        }
+
+        if (best.boundary_element < 0) {
+            return std::unexpected(DeformationFailure{.reason = DeformationFailureReason::unresolved_surface, .message = "No logical quadrilateral contains a radial projection onto the required shell."});
+        }
+
+        return best;
+    }
+}