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

feat(geometry/deformation): Added deformation subsystem

Browse Source 
The model must be able to deform to capture the shape induced by rotation or an external potential. Here we add the precription to descibe how we allow deformation. Note that the code to define how deformation is actually implemented is in the backend.
This commit is contained in:
Emily Boudreaux
2026-11-09 13:11:29 -05:00
parent 9760e46477
commit 4bbbf1c418
29 changed files with 1705 additions and 269 deletions
Download Patch File Download Diff File Expand all files Collapse all files

126
src/lib/serif/geometry/deformation/backend/mfem/radial_locator.cpp
View File

@@ -1,126 +0,0 @@
#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;
}
}

157
src/lib/serif/geometry/deformation/determinant_polynomial.cpp Normal file
View File

@@ -0,0 +1,157 @@
#include "serif/geometry/deformation/detail/determinant_polynomial.hpp"
#include <algorithm>
#include <bit>
#include <cmath>
#include <limits>
#include <numeric>
namespace serif::geometry::deformation::detail {
namespace {
/** @brief Use determinant multilinearity to separate powers of the displacement direction. */
[[nodiscard]] double selected_column_determinant(const JacobianEntries& base, const JacobianEntries& direction, const unsigned int direction_columns, const int dimension, const double maximum_step) noexcept {
const auto entry = [&](const int row, const int column) {
if ((direction_columns & (1U << static_cast<unsigned int>(column))) != 0U) {
return maximum_step * direction[row][column];
}
return base[row][column];
};
if (dimension == 1) {
return entry(0, 0);
}
if (dimension == 2) {
return entry(0, 0) * entry(1, 1) - entry(0, 1) * entry(1, 0);
}
return entry(0, 0) * (entry(1, 1) * entry(2, 2) - entry(1, 2) * entry(2, 1)) -
entry(0, 1) * (entry(1, 0) * entry(2, 2) - entry(1, 2) * entry(2, 0)) +
entry(0, 2) * (entry(1, 0) * entry(2, 1) - entry(1, 1) * entry(2, 0));
}
/** @brief Ignore numerically absent leading powers using the MeanField relative tolerance. */
[[nodiscard]] int polynomial_degree(const DeterminantPolynomial& coefficients) noexcept {
double scale = 0.0;
for (const double coefficient : coefficients) {
scale = std::max(scale, std::abs(coefficient));
}
const double tolerance = 64.0 * std::numeric_limits<double>::epsilon() * scale;
for (int degree = 3; degree > 0; --degree) {
if (std::abs(coefficients[degree]) > tolerance) {
return degree;
}
}
return 0;
}
/** @brief Keep distinct critical points strictly inside the interval being tested. */
void append_unit_interval_root(std::array<double, 2>& roots, int& root_count, const double root) noexcept {
if (!std::isfinite(root) || root <= 0.0 || root >= 1.0) {
return;
}
if (root_count > 0 && std::abs(root - roots[0]) <= 64.0 * std::numeric_limits<double>::epsilon()) {
return;
}
roots[root_count] = root;
++root_count;
}
/** @brief Split the cubic into monotone intervals without losing a repeated determinant root. */
[[nodiscard]] int derivative_critical_points(const DeterminantPolynomial& coefficients, const int degree, std::array<double, 2>& critical_points) noexcept {
int count = 0;
if (degree == 2) {
append_unit_interval_root(critical_points, count, -coefficients[1] / (2.0 * coefficients[2]));
} else if (degree == 3) {
const double quadratic = 3.0 * coefficients[3];
const double linear = 2.0 * coefficients[2];
const double constant = coefficients[1];
const double discriminant = std::fma(linear, linear, -4.0 * quadratic * constant);
const double discriminant_scale = linear * linear + std::abs(4.0 * quadratic * constant);
const double discriminant_tolerance = 64.0 * std::numeric_limits<double>::epsilon() * discriminant_scale;
if (discriminant >= -discriminant_tolerance) {
const double square_root = std::sqrt(std::max(0.0, discriminant));
if (square_root == 0.0) {
append_unit_interval_root(critical_points, count, -linear / (2.0 * quadratic));
} else {
const double q = -0.5 * (linear + std::copysign(square_root, linear));
append_unit_interval_root(critical_points, count, q / quadratic);
append_unit_interval_root(critical_points, count, constant / q);
}
}
}
std::sort(critical_points.begin(), critical_points.begin() + count);
return count;
}
/** @brief Retain the nonpositive side of a bracket, avoiding an overestimate of the safe interval. */
[[nodiscard]] double bisect_first_nonpositive_value(const DeterminantPolynomial& coefficients, double lower, double upper) noexcept {
for (int iteration = 0; iteration < 80; ++iteration) {
const double middle = std::midpoint(lower, upper);
if (evaluate_polynomial(coefficients, middle) > 0.0) {
lower = middle;
} else {
upper = middle;
}
}
return upper;
}
}
DeterminantPolynomial determinant_polynomial(const JacobianEntries& base, const JacobianEntries& direction, const int dimension, const double maximum_step, const double determinant_floor) noexcept {
DeterminantPolynomial coefficients{};
const unsigned int term_count = 1U << static_cast<unsigned int>(dimension);
for (unsigned int mask = 0; mask < term_count; ++mask) {
coefficients[std::popcount(mask)] += selected_column_determinant(base, direction, mask, dimension, maximum_step);
}
coefficients[0] -= determinant_floor;
return coefficients;
}
double evaluate_polynomial(const DeterminantPolynomial& coefficients, const double parameter) noexcept {
return std::fma(parameter, std::fma(parameter, std::fma(parameter, coefficients[3], coefficients[2]), coefficients[1]), coefficients[0]);
}
double first_boundary_parameter(const DeterminantPolynomial& input_coefficients) noexcept {
// Scaling leaves roots and the relative tolerances unchanged, and prevents overflow in the cubic discriminant.
double scale = 0.0;
for (const double coefficient : input_coefficients) {
scale = std::max(scale, std::abs(coefficient));
}
DeterminantPolynomial coefficients{};
for (std::size_t index = 0; index < coefficients.size(); ++index) {
coefficients[index] = input_coefficients[index] / scale;
}
const int degree = polynomial_degree(coefficients);
if (degree == 0) {
return std::numeric_limits<double>::infinity();
}
double coefficient_scale = 0.0;
for (const double coefficient : coefficients) {
coefficient_scale += std::abs(coefficient);
}
const double value_tolerance = 128.0 * std::numeric_limits<double>::epsilon() * coefficient_scale;
std::array<double, 2> critical_points{};
const int critical_point_count = derivative_critical_points(coefficients, degree, critical_points);
std::array<double, 4> interval_ends{};
for (int index = 0; index < critical_point_count; ++index) {
interval_ends[index + 1] = critical_points[index];
}
interval_ends[critical_point_count + 1] = 1.0;
for (int interval = 0; interval <= critical_point_count; ++interval) {
const double lower = interval_ends[interval];
const double upper = interval_ends[interval + 1];
const double upper_value = evaluate_polynomial(coefficients, upper);
if (upper_value <= 0.0) {
return bisect_first_nonpositive_value(coefficients, lower, upper);
}
if (upper_value <= value_tolerance) {
// A tangent zero may evaluate a few ulps above zero at its derivative root.
return upper;
}
}
return std::numeric_limits<double>::infinity();
}
}

124
src/lib/serif/geometry/deformation/harmonics.cpp Normal file
View File

@@ -0,0 +1,124 @@
#include "serif/geometry/deformation/detail/harmonics.hpp"
#include <algorithm>
#include <cmath>
#include <cstdint>
#include <limits>
#include <numbers>
namespace serif::geometry::deformation::detail {
namespace {
/** @brief Keep accumulation in extended precision when checking sampled independence. */
[[nodiscard]] double column_product(std::span<const double> basis, std::size_t rows, std::size_t columns, std::size_t first, std::size_t second) {
long double result = 0.0L;
for (std::size_t row = 0; row < rows; ++row) {
result += static_cast<long double>(basis[row * columns + first]) * basis[row * columns + second];
}
return static_cast<double>(result);
}
/** @brief Build normalized associated Legendre functions without factorial overflow. */
void evaluate_legendre(const std::array<double, 3>& direction, int maximum_degree, int maximum_order, std::span<double> legendre) {
const std::size_t stride = static_cast<std::size_t>(maximum_order) + 1;
const double cylindrical_radius = std::hypot(direction[0], direction[1]);
/// The triangular recurrence uses a rectangular workspace to keep indexing explicit.
const auto polynomial = [&](int degree, int order) -> double& {
return legendre[static_cast<std::size_t>(degree) * stride + static_cast<std::size_t>(order)];
};
polynomial(0, 0) = 1.0 / std::sqrt(4.0 * std::numbers::pi);
for (int order = 0; order <= maximum_order; ++order) {
if (order > 0) {
polynomial(order, order) = -std::sqrt((2.0 * order + 1.0) / (2.0 * order)) * cylindrical_radius * polynomial(order - 1, order - 1);
}
if (order < maximum_degree) {
polynomial(order + 1, order) = std::sqrt(2.0 * order + 3.0) * direction[2] * polynomial(order, order);
}
for (int degree = order + 2; degree <= maximum_degree; ++degree) {
const double degree_value = degree;
const double order_value = order;
const double leading = std::sqrt((4.0 * degree_value * degree_value - 1.0) / (degree_value * degree_value - order_value * order_value));
const double trailing = std::sqrt(((2.0 * degree_value + 1.0) * ((degree_value - 1.0) * (degree_value - 1.0) - order_value * order_value)) / ((2.0 * degree_value - 3.0) * (degree_value * degree_value - order_value * order_value)));
polynomial(degree, order) = leading * direction[2] * polynomial(degree - 1, order) - trailing * polynomial(degree - 2, order);
}
}
}
}
std::expected<std::vector<SphericalHarmonicMode>, DeformationFailure> make_harmonic_modes(const SphericalHarmonicSurface& model, std::size_t node_count) {
const int maximum_order = model.maximum_azimuthal_order.value_or(model.maximum_degree);
if (model.maximum_degree < 0 || maximum_order < 0 || maximum_order > model.maximum_degree) {
return std::unexpected(DeformationFailure{.reason = DeformationFailureReason::invalid_configuration, .message = "The maximum harmonic degree must be nonnegative and the azimuthal cutoff must lie between zero and that degree."});
}
const std::int64_t degree = model.maximum_degree;
const std::int64_t order = maximum_order;
const std::int64_t count = (order + 1) * (order + 1) + (degree - order) * (2 * order + 1);
if (count > std::numeric_limits<int>::max()) {
return std::unexpected(DeformationFailure{.reason = DeformationFailureReason::invalid_configuration, .message = "The requested harmonic mode count exceeds supported vector indexing."});
}
if (static_cast<std::size_t>(count) > node_count) {
return std::unexpected(DeformationFailure{.reason = DeformationFailureReason::unresolved_surface, .message = "The requested harmonic mode count exceeds the number of surface nodes; reduce the degree or refine the surface."});
}
std::vector<SphericalHarmonicMode> modes;
modes.reserve(static_cast<std::size_t>(count));
for (int degree_index = 0; degree_index <= model.maximum_degree; ++degree_index) {
const int cutoff = std::min(degree_index, maximum_order);
for (int order_index = -cutoff; order_index <= cutoff; ++order_index) {
modes.push_back(SphericalHarmonicMode{.degree = degree_index, .order = order_index});
}
}
return modes;
}
std::expected<std::vector<double>, DeformationFailure> sample_harmonic_basis(const SphericalHarmonicSurface& model, std::span<const SphericalHarmonicMode> modes, std::span<const std::array<double, 3>> directions) {
if (modes.empty() || directions.size() > std::vector<double>{}.max_size() / modes.size()) {
return std::unexpected(DeformationFailure{.reason = DeformationFailureReason::invalid_size, .message = "The requested sampled harmonic basis exceeds supported storage indexing."});
}
const int maximum_order = model.maximum_azimuthal_order.value_or(model.maximum_degree);
const std::size_t stride = static_cast<std::size_t>(maximum_order) + 1;
std::vector<double> legendre((static_cast<std::size_t>(model.maximum_degree) + 1) * stride, 0.0);
std::vector<double> basis(directions.size() * modes.size());
for (std::size_t node = 0; node < directions.size(); ++node) {
evaluate_legendre(directions[node], model.maximum_degree, maximum_order, legendre);
const double azimuth = std::atan2(directions[node][1], directions[node][0]);
for (std::size_t column = 0; column < modes.size(); ++column) {
const SphericalHarmonicMode mode = modes[column];
double value = legendre[static_cast<std::size_t>(mode.degree) * stride + static_cast<std::size_t>(std::abs(mode.order))];
if (mode.order < 0) {
value *= std::sqrt(2.0) * std::sin(-mode.order * azimuth);
}
if (mode.order > 0) {
value *= std::sqrt(2.0) * std::cos(mode.order * azimuth);
}
if (!std::isfinite(value)) {
return std::unexpected(DeformationFailure{.reason = DeformationFailureReason::non_finite, .message = "Sampling the spherical harmonic basis produced a non-finite value."});
}
basis[node * modes.size() + column] = value;
}
}
return basis;
}
DeformationStatus check_independent_harmonics(std::span<const double> basis, std::size_t node_count, std::size_t mode_count) {
std::vector<double> orthogonal(basis.begin(), basis.end());
constexpr double minimum_relative_pivot = 1e-12;
for (std::size_t column = 0; column < mode_count; ++column) {
const double original = std::sqrt(column_product(basis, node_count, mode_count, column, column));
for (int pass = 0; pass < 2; ++pass) {
for (std::size_t previous = 0; previous < column; ++previous) {
const double projection = column_product(orthogonal, node_count, mode_count, previous, column);
for (std::size_t row = 0; row < node_count; ++row) {
orthogonal[row * mode_count + column] -= projection * orthogonal[row * mode_count + previous];
}
}
}
const double pivot = std::sqrt(column_product(orthogonal, node_count, mode_count, column, column));
if (!std::isfinite(original) || !(original > 0.0) || !std::isfinite(pivot) || !(pivot > minimum_relative_pivot * original)) {
return std::unexpected(DeformationFailure{.reason = DeformationFailureReason::unresolved_surface, .message = "The selected harmonics are numerically dependent on the surface nodes; reduce their degree or order, or refine the surface."});
}
for (std::size_t row = 0; row < node_count; ++row) {
orthogonal[row * mode_count + column] /= pivot;
}
}
return {};
}
}

104
src/lib/serif/geometry/deformation/radial_extension.cpp Normal file
View File

@@ -0,0 +1,104 @@
#include <algorithm>
#include <cmath>
#include <limits>
#include <utility>
#include "serif/geometry/deformation/radial_extension.hpp"
#include "serif/geometry/deformation/detail/action_storage.hpp"
namespace serif::geometry::deformation
{
RadialExtension::RadialExtension(std::shared_ptr<const RadialInterpolation> interpolation, std::vector<double> radial_weights)
: m_interpolation(std::move(interpolation)), m_radial_weights(std::move(radial_weights))
{
}
DeformationStatus RadialExtension::apply(std::span<const double> surface_displacement, std::span<double> volume_displacement, Workspace) const
{
const DeformationStatus status = detail::validate_action_arguments(volume_displacement, output_size(), {{surface_displacement, input_size()}});
if (!status)
{
return status;
}
std::fill(volume_displacement.begin(), volume_displacement.end(), 0.0);
const std::size_t count = m_radial_weights.size();
const RadialInterpolationData& reference = m_interpolation->data();
for (std::size_t node = 0; node < count; ++node)
{
if (m_radial_weights[node] == 0.0)
{
continue;
}
for (std::size_t component = 0; component < 3; ++component)
{
double value = 0.0;
for (int entry = reference.row_offsets[node]; entry < reference.row_offsets[node + 1]; ++entry)
{
value += reference.interpolation_weights[entry] * surface_displacement[3 * reference.surface_columns[entry] + component];
}
volume_displacement[node + component * count] = m_radial_weights[node] * value;
}
}
return validate_deformation_values(volume_displacement, output_size());
}
DeformationStatus RadialExtension::apply_transpose(std::span<const double> volume_force, std::span<double> surface_force, Workspace) const
{
const DeformationStatus status = detail::validate_action_arguments(surface_force, input_size(), {{volume_force, output_size()}});
if (!status)
{
return status;
}
std::fill(surface_force.begin(), surface_force.end(), 0.0);
const std::size_t count = m_radial_weights.size();
const RadialInterpolationData& reference = m_interpolation->data();
for (std::size_t node = 0; node < count; ++node)
{
if (m_radial_weights[node] == 0.0)
{
continue;
}
for (int entry = reference.row_offsets[node]; entry < reference.row_offsets[node + 1]; ++entry)
{
for (std::size_t component = 0; component < 3; ++component)
{
surface_force[3 * reference.surface_columns[entry] + component] += m_radial_weights[node] * reference.interpolation_weights[entry] * volume_force[node + component * count];
}
}
}
return validate_deformation_values(surface_force, input_size());
}
std::expected<RadialExtension, DeformationFailure> make_radial_extension(std::shared_ptr<const RadialInterpolation> interpolation, PowerLawRadialExtension interior, FixedInfinityRadialExtension)
{
if (!interpolation || !std::isfinite(interior.radial_power) || interior.radial_power < 1.0)
{
return std::unexpected(DeformationFailure{.reason = DeformationFailureReason::invalid_configuration, .message = "Radial extension requires interpolation data and a finite stellar power of at least one."});
}
const RadialInterpolationData& reference = interpolation->data();
const std::size_t count = interpolation->scalar_dof_count();
RadialExtension result(std::move(interpolation), std::vector<double>(count, 0.0));
for (std::size_t node = 0; node < count; ++node)
{
const double radius = reference.logical_radii[node];
if (reference.stellar_support[node])
{
result.m_radial_weights[node] = std::pow(radius / reference.stellar_radius, interior.radial_power);
}
else if (reference.vacuum_support[node])
{
result.m_radial_weights[node] = (reference.infinity_radius - radius) / (reference.infinity_radius - reference.stellar_radius);
}
if (!std::isfinite(result.m_radial_weights[node]))
{
return std::unexpected(DeformationFailure{.reason = DeformationFailureReason::non_finite, .message = "A radial extension weight is non-finite."});
}
// Compilation omits interpolation at the fixed center using this same roundoff threshold.
const double center_tolerance = 64.0 * std::numeric_limits<double>::epsilon() * reference.infinity_radius;
if (radius <= center_tolerance && reference.stellar_support[node])
{
result.m_radial_weights[node] = 0.0;
}
}
return result;
}
} // namespace serif::geometry::deformation

68
src/lib/serif/geometry/deformation/radial_interpolation.cpp Normal file
View File

@@ -0,0 +1,68 @@
#include <cmath>
#include <limits>
#include "serif/geometry/deformation/radial_interpolation.hpp"
namespace serif::geometry::deformation
{
std::expected<RadialInterpolation, DeformationFailure> make_radial_interpolation(RadialInterpolationData data)
{
const std::size_t count = data.logical_radii.size();
if (count == 0 || data.surface_scalar_dofs.size() == 0 || data.stellar_support.size() != count || data.vacuum_support.size() != count || data.row_offsets.size() != count + 1 || !(data.stellar_radius > 0.0) || !(data.infinity_radius > data.stellar_radius) || !std::isfinite(data.infinity_radius))
{
return std::unexpected(DeformationFailure{.reason = DeformationFailureReason::invalid_configuration, .message = "Radial interpolation supports, rows and shell radii are inconsistent."});
}
if (data.surface_columns.size() != data.interpolation_weights.size() || data.row_offsets.front() != 0 || data.row_offsets.back() < 0 || static_cast<std::size_t>(data.row_offsets.back()) != data.surface_columns.size())
{
return std::unexpected(DeformationFailure{.reason = DeformationFailureReason::invalid_configuration, .message = "Radial interpolation row offsets do not span the stored entries."});
}
int previous_surface_dof = -1;
for (const int surface_dof : data.surface_scalar_dofs)
{
if (surface_dof <= previous_surface_dof || static_cast<std::size_t>(surface_dof) >= count)
{
return std::unexpected(DeformationFailure{.reason = DeformationFailureReason::invalid_configuration, .message = "Radial surface nodes must be distinct, ascending and inside the volume scalar space."});
}
previous_surface_dof = surface_dof;
}
const double radial_tolerance = 128.0 * std::numeric_limits<double>::epsilon() * data.infinity_radius;
for (std::size_t node = 0; node < count; ++node)
{
const int stellar = data.stellar_support[node];
const int vacuum = data.vacuum_support[node];
const double radius = data.logical_radii[node];
if ((stellar != 0 && stellar != 1) || (vacuum != 0 && vacuum != 1) || (stellar && radius > data.stellar_radius + radial_tolerance) || (vacuum && (radius < data.stellar_radius - radial_tolerance || radius > data.infinity_radius + radial_tolerance)))
{
return std::unexpected(DeformationFailure{.reason = DeformationFailureReason::invalid_configuration, .message = "Radial node support and logical shell radii are inconsistent."});
}
if (data.row_offsets[node] > data.row_offsets[node + 1] || !std::isfinite(data.logical_radii[node]) || data.logical_radii[node] < 0.0)
{
return std::unexpected(DeformationFailure{.reason = DeformationFailureReason::invalid_configuration, .message = "Radial interpolation rows must increase and logical radii must be finite and nonnegative."});
}
}
for (std::size_t entry = 0; entry < data.surface_columns.size(); ++entry)
{
if (data.surface_columns[entry] < 0 || static_cast<std::size_t>(data.surface_columns[entry]) >= data.surface_scalar_dofs.size() || !std::isfinite(data.interpolation_weights[entry]))
{
return std::unexpected(DeformationFailure{.reason = DeformationFailureReason::invalid_configuration, .message = "A radial interpolation entry has an invalid column or weight."});
}
}
const double center_tolerance = 64.0 * std::numeric_limits<double>::epsilon() * data.infinity_radius;
for (std::size_t node = 0; node < count; ++node)
{
const bool moving = data.stellar_support[node] ? data.logical_radii[node] > center_tolerance : data.vacuum_support[node] && data.logical_radii[node] != data.infinity_radius;
if (moving)
{
double partition = 0.0;
for (int entry = data.row_offsets[node]; entry < data.row_offsets[node + 1]; ++entry)
{
partition += data.interpolation_weights[entry];
}
if (!std::isfinite(partition) || std::abs(partition - 1.0) > 2.0e-12)
{
return std::unexpected(DeformationFailure{.reason = DeformationFailureReason::invalid_configuration, .message = "A moving radial interpolation row must preserve a constant surface trace."});
}
}
}
return RadialInterpolation(std::move(data));
}
} // namespace serif::geometry::deformation

210
src/lib/serif/geometry/deformation/surface_embedding.cpp Normal file
View File

@@ -0,0 +1,210 @@
#include "serif/geometry/deformation/surface_embedding.hpp"
#include <algorithm>
#include <cmath>
#include <limits>
#include <utility>
#include "serif/geometry/deformation/detail/harmonics.hpp"
#include "serif/geometry/deformation/detail/action_storage.hpp"
namespace serif::geometry::deformation
{
namespace
{
/** @brief A finite input may overflow during summation; keep that failure distinct from an accepted deformation. */
[[nodiscard]] std::unexpected<DeformationFailure> action_overflow()
{
return std::unexpected(DeformationFailure{.reason = DeformationFailureReason::non_finite, .message = "Surface action produced a non-finite result; reduce the supplied parameter or force magnitude."});
}
/** @brief Accumulate radial virtual work without introducing surface-area or mass-matrix weights. */
[[nodiscard]] long double radial_work(const std::array<double, 3>& direction, std::span<const double> force)
{
long double result = 0.0L;
for (std::size_t component = 0; component < 3; ++component)
{
result += static_cast<long double>(direction[component]) * force[component];
}
return result;
}
} // namespace
std::size_t SurfaceEmbedding::input_size() const noexcept
{
return has_nodal_parameters() ? surface_node_count() : modes_.size();
}
std::size_t SurfaceEmbedding::output_size() const noexcept
{
return 3 * surface_node_count();
}
std::size_t SurfaceEmbedding::surface_node_count() const noexcept
{
return radial_directions_.size();
}
bool SurfaceEmbedding::has_nodal_parameters() const noexcept
{
return modes_.empty();
}
std::span<const std::array<double, 3>> SurfaceEmbedding::radial_directions() const noexcept
{
return radial_directions_;
}
std::span<const double> SurfaceEmbedding::reference_radii() const noexcept
{
return reference_radii_;
}
std::span<const std::array<double, 3>> SurfaceEmbedding::reference_positions() const noexcept
{
return reference_positions_;
}
std::span<const SphericalHarmonicMode> SurfaceEmbedding::modes() const noexcept
{
return modes_;
}
std::span<const double> SurfaceEmbedding::nodal_basis() const noexcept
{
return nodal_basis_;
}
DeformationStatus SurfaceEmbedding::apply(std::span<const double> parameters, std::span<double> surface_displacement, Workspace) const
{
if (const DeformationStatus status = detail::validate_action_arguments(surface_displacement, output_size(), {{parameters, input_size()}}); !status)
{
return status;
}
for (std::size_t node = 0; node < surface_node_count(); ++node)
{
long double amplitude = 0.0L;
if (has_nodal_parameters())
{
amplitude = parameters[node];
}
else
{
for (std::size_t mode = 0; mode < input_size(); ++mode)
{
amplitude += static_cast<long double>(nodal_basis_[node * input_size() + mode]) * parameters[mode];
}
}
const double radial_amplitude = static_cast<double>(amplitude);
if (!std::isfinite(radial_amplitude))
{
return action_overflow();
}
for (std::size_t component = 0; component < 3; ++component)
{
const double displacement = static_cast<double>(static_cast<long double>(radial_amplitude) * radial_directions_[node][component]);
if (!std::isfinite(displacement))
{
return action_overflow();
}
surface_displacement[3 * node + component] = displacement;
}
}
return {};
}
DeformationStatus SurfaceEmbedding::apply_transpose(std::span<const double> surface_force, std::span<double> parameter_force, Workspace) const
{
if (const DeformationStatus status = detail::validate_action_arguments(parameter_force, input_size(), {{surface_force, output_size()}}); !status)
{
return status;
}
for (std::size_t parameter = 0; parameter < input_size(); ++parameter)
{
long double work = 0.0L;
if (has_nodal_parameters())
{
work = radial_work(radial_directions_[parameter], surface_force.subspan(3 * parameter, 3));
}
else
{
for (std::size_t node = 0; node < surface_node_count(); ++node)
{
work += static_cast<long double>(nodal_basis_[node * input_size() + parameter]) * radial_work(radial_directions_[node], surface_force.subspan(3 * node, 3));
}
}
const double parameter_work = static_cast<double>(work);
if (!std::isfinite(parameter_work))
{
return action_overflow();
}
parameter_force[parameter] = parameter_work;
}
return {};
}
std::expected<SurfaceEmbedding, DeformationFailure> build_surface_embedding(const NodalRadialSurface& model, std::span<const std::array<double, 3>> positions)
{
if (positions.empty() || positions.size() > std::numeric_limits<std::size_t>::max() / 3)
{
return std::unexpected(DeformationFailure{.reason = DeformationFailureReason::invalid_size, .message = "A radial surface requires a nonempty, representable list of reference nodes."});
}
if (!std::ranges::all_of(model.reference_center, [](double value)
{
return std::isfinite(value);
}))
{
return std::unexpected(DeformationFailure{.reason = DeformationFailureReason::non_finite, .message = "The reference center must contain three finite coordinates."});
}
SurfaceEmbedding result;
result.reference_positions_.assign(positions.begin(), positions.end());
result.radial_directions_.reserve(positions.size());
result.reference_radii_.reserve(positions.size());
for (const std::array<double, 3>& position : positions)
{
std::array<double, 3> direction{};
for (std::size_t component = 0; component < 3; ++component)
{
direction[component] = position[component] - model.reference_center[component];
}
const double radius = std::hypot(direction[0], direction[1], direction[2]);
if (!std::isfinite(radius) || radius <= 0.0)
{
return std::unexpected(DeformationFailure{.reason = DeformationFailureReason::invalid_geometry, .message = "Every reference surface node must have a finite positive distance from the reference center."});
}
for (double& component : direction)
{
component /= radius;
}
result.radial_directions_.push_back(direction);
result.reference_radii_.push_back(radius);
}
return result;
}
std::expected<SurfaceEmbedding, DeformationFailure> build_surface_embedding(const SphericalHarmonicSurface& model, std::span<const std::array<double, 3>> positions)
{
auto modes = detail::make_harmonic_modes(model, positions.size());
if (!modes)
{
return std::unexpected(modes.error());
}
auto result = build_surface_embedding(NodalRadialSurface{.reference_center = model.reference_center}, positions);
if (!result)
{
return result;
}
auto basis = detail::sample_harmonic_basis(model, *modes, result->radial_directions());
if (!basis)
{
return std::unexpected(basis.error());
}
if (const DeformationStatus status = detail::check_independent_harmonics(*basis, positions.size(), modes->size()); !status)
{
return std::unexpected(status.error());
}
result->modes_ = std::move(*modes);
result->nodal_basis_ = std::move(*basis);
return result;
}
} // namespace serif::geometry::deformation