diff --git a/src/include/serif/geometry/deformation/actions.hpp b/src/include/serif/geometry/deformation/actions.hpp new file mode 100644 index 00000000..f41a5d27 --- /dev/null +++ b/src/include/serif/geometry/deformation/actions.hpp @@ -0,0 +1,156 @@ +#pragma once + +#include +#include +#include + +#include "serif/geometry/deformation/failure.hpp" +#include "serif/geometry/deformation/detail/action_storage.hpp" +#include "serif/geometry/deformation/spaces.hpp" +#include "serif/geometry/deformation/workspace.hpp" + +namespace serif::geometry::deformation +{ + /** @brief Linearity is an explicit property of the mathematics, independent of its action interface. */ + template + inline constexpr bool is_linear_deformation = requires + { + requires Map::is_linear; + }; + + /** @brief Matching lengths never replace the physical meanings of the input and output spaces. */ + template + concept DeformationDimensions = requires(const Map& map) + { + typename Map::Workspace; + typename Map::InputSpace; + typename Map::OutputSpace; + { map.input_size() } -> std::same_as; + { map.output_size() } -> std::same_as; + }; + + /** @brief Fixed maps provide just forward and transpose actions with caller-owned scratch. */ + template + concept LinearDeformationMap = DeformationDimensions && is_linear_deformation && + (std::same_as, EmptyDeformationWorkspace> || requires(const Map& map) + { + { map.make_workspace() } -> std::same_as>; + }) && requires(const Map& map, std::span input, std::span output, DeformationWorkspace& workspace) + { + { map.apply(input, output, workspace) } -> std::same_as; + { map.apply_transpose(input, output, workspace) } -> std::same_as; + }; + + /** @brief Nonlinear maps explicitly provide every term needed by a Newton force pullback. */ + template + concept DifferentiableDeformationMap = DeformationDimensions && requires(const Map& map, std::span state, std::span direction, std::span dual, std::span output, DeformationWorkspace& workspace) + { + { map.make_workspace() } -> std::same_as>; + { map.build_displacement(state, output, workspace) } -> std::same_as; + { map.apply_jacobian(state, direction, output, workspace) } -> std::same_as; + { map.apply_jacobian_transpose(state, dual, output, workspace) } -> std::same_as; + { map.apply_pullback_derivative(state, direction, dual, output, workspace) } -> std::same_as; + }; + + template + concept DeformationMap = LinearDeformationMap || DifferentiableDeformationMap; + + /** @brief Allocate action scratch explicitly, once per independent caller. */ + template + [[nodiscard]] DeformationWorkspace make_deformation_workspace(const Map& map) + { + if constexpr (requires { map.make_workspace(); }) + { + return map.make_workspace(); + } + else + { + return {}; + } + } + + /** + * @brief Evaluate physical displacement with caller-owned output and scratch. + * All action outputs must be disjoint from their inputs. Workspace storage is reserved + * for the map and must not supply public input or output spans. + */ + template + [[nodiscard]] DeformationStatus build_displacement(const Map& map, std::span parameters, std::span displacement, DeformationWorkspace& workspace) + { + if (const auto status = detail::validate_disjoint_output(displacement, {parameters}); !status) + { + return status; + } + if constexpr (LinearDeformationMap) + { + return map.apply(parameters, displacement, workspace); + } + else + { + return map.build_displacement(parameters, displacement, workspace); + } + } + + /** @brief Linear maps reuse their primal action; nonlinear maps supply their exact tangent. */ + template + [[nodiscard]] DeformationStatus apply_jacobian(const Map& map, std::span parameters, std::span direction, std::span displacement_direction, DeformationWorkspace& workspace) + { + if (const DeformationStatus status = detail::validate_action_arguments(displacement_direction, map.output_size(), {{parameters, map.input_size()}}); !status) + { + return status; + } + if (const DeformationStatus status = detail::validate_disjoint_output(displacement_direction, {direction}); !status) + { + return status; + } + if constexpr (LinearDeformationMap) + { + return map.apply(direction, displacement_direction, workspace); + } + else + { + return map.apply_jacobian(parameters, direction, displacement_direction, workspace); + } + } + + /** @brief Transfer virtual work from the displacement field to its independent shape coefficients. */ + template + [[nodiscard]] DeformationStatus apply_jacobian_transpose(const Map& map, std::span parameters, std::span displacement_dual, std::span parameter_dual, DeformationWorkspace& workspace) + { + if (const DeformationStatus status = detail::validate_action_arguments(parameter_dual, map.input_size(), {{parameters, map.input_size()}}); !status) + { + return status; + } + if (const DeformationStatus status = detail::validate_disjoint_output(parameter_dual, {displacement_dual}); !status) + { + return status; + } + if constexpr (LinearDeformationMap) + { + return map.apply_transpose(displacement_dual, parameter_dual, workspace); + } + else + { + return map.apply_jacobian_transpose(parameters, displacement_dual, parameter_dual, workspace); + } + } + + /** @brief Differentiate the transpose at fixed dual; zero is justified only for a declared linear map. */ + template + [[nodiscard]] DeformationStatus apply_pullback_derivative(const Map& map, std::span parameters, std::span direction, std::span displacement_dual, std::span parameter_dual, DeformationWorkspace& workspace) + { + if (const DeformationStatus status = detail::validate_action_arguments(parameter_dual, map.input_size(), {{parameters, map.input_size()}, {direction, map.input_size()}, {displacement_dual, map.output_size()}}); !status) + { + return status; + } + if constexpr (is_linear_deformation) + { + std::ranges::fill(parameter_dual, 0.0); + return {}; + } + else + { + return map.apply_pullback_derivative(parameters, direction, displacement_dual, parameter_dual, workspace); + } + } +} // namespace serif::geometry::deformation diff --git a/src/include/serif/geometry/deformation/backend/mfem/detail/radial_locator.hpp b/src/include/serif/geometry/deformation/backend/mfem/detail/radial_locator.hpp deleted file mode 100644 index a82f0864..00000000 --- a/src/include/serif/geometry/deformation/backend/mfem/detail/radial_locator.hpp +++ /dev/null @@ -1,79 +0,0 @@ -#pragma once - -/** - * @file radial_locator.hpp - * @author Emily M. Boudreaux - * - * @brief File containing some utility functions for locating points on a radial shell. - */ - -#include "mfem/fem/intrules.hpp" -#include "mfem/linalg/vector.hpp" -#include "serif/geometry/deformation/failure.hpp" - -namespace mfem { - /** - * @brief Forward definition of the mfem::Mesh class to avoid having to parse the entire mfem mesh header file. - */ - class Mesh; -} - -namespace serif::geometry::deformation::backend::mfem::detail { - /** @struct LocatedBoundaryPoint - * - * @brief A structure representing a located boundary point, including the boundary element index and the integration point on that element. - */ - struct LocatedBoundaryPoint { - int boundary_element; ///< The index of the boundary element where the point is located. - ::mfem::IntegrationPoint point; ///< The integration point on the boundary element where the point is located. - }; - - /** - * @brief compute the cubical local radius of a given position vector. Alternative put: - * compute the \f$L\infty\f$ norm of the position vector (i.e. the maximum absolute value of its components) - * - * @param position The position vector for which to compute the logical radius. - * @returns The logical radius of the position vector, defined as the maximum absolute value of its components (i.e., the L\f$\infty\f$ norm). - */ - [[nodiscard]] double logical_radius(const ::mfem::Vector& position); - - - /** - * @brief Compute the logical radius of a boundary shell in a mesh, given the mesh and the boundary attribute. - * - * @param mesh The mesh containing the boundary shell. - * @param attribute The boundary attribute identifying the shell to compute the logical radius for. - * - * @returns The logical radius of the boundary shell if successful, or a DeformationFailure object if an error occurs (e.g., incompatible discretization, non-finite vertex, or invalid geometry). - * - * @par Explanation - * The algorithm here is as follows - * 1. Loop over all boundary elements. validate they they have the correct attribute id (if not skip them) - * 2. Get the geometry type of the element and validate that it has a valid transformation and is a square (which for mfem really means a quadrilateral). If not return an error - * 3. Get the vertices of the element and loop over them. For each vertex - * a. Transform the vertex into a position vector - * b. validate that each component of the position vector is finite. If not return an error - * c. Compute the logical radius of the position vector (L\f$\infty\f$ norm) and update the minimum and maximum logical radius values seen so far - * d. Increment the number of samples seen so far - * 4. validate that at least one sample was seen, that the minimum values is positive, and that are at least separated by 1e-12 * max(1, maximum). If not return an error - * 5. Return the average of the minimum and maximum logical radius values seen. - * - * Effectively then what this does is find an estimate for the radius of some enclosing surface over the boundary elements. - */ - [[nodiscard]] std::expected boundary_logical_radius(::mfem::Mesh& mesh, int attribute); - - /** - * @brief Locate a projected logical ray on a quadrilateral shell. - * - * @param mesh The mesh containing the boundary shell. - * @param attribute The boundary attribute identifying the shell to locate the point on. - * @param target The target position vector to locate on the boundary shell. - * - * @returns A LocatedBoundaryPoint object representing the located boundary point if successful, or a DeformationFailure object if an error occurs - */ - [[nodiscard]] std::expected locate_boundary_point( - ::mfem::Mesh& mesh, - int attribute, - const ::mfem::Vector& target - ); -} diff --git a/src/include/serif/geometry/deformation/backend/mfem/detail/reference_tracking.hpp b/src/include/serif/geometry/deformation/backend/mfem/detail/reference_tracking.hpp deleted file mode 100644 index dcf10da6..00000000 --- a/src/include/serif/geometry/deformation/backend/mfem/detail/reference_tracking.hpp +++ /dev/null @@ -1,64 +0,0 @@ -#pragma once - -#include - -namespace mfem -{ - //TODO: Really we should move the forward declarations of these mfem types to some unified location so that we are - // not managing haphazardly on a per-translation unit level (its really easy to imagining how this would lead - // to an ORD violation). Further, most of the time in headers we really don't need the implementation of the type - // so this does save a lot of compile time. - // I would suggest something like a utils/types/mfem/forward.hpp module which just forward defines these - // things and then that header can be included. - - - class Mesh; ///< Forward declaration of the mfem::Mesh class to avoid including the entire mfem mesh header file. - class FiniteElementSpace; ///< Forward declaration of the mfem::FiniteElementSpace class to avoid including the entire mfem finite element space header file. - class FiniteElementCollection; ///< Forward declaration of the mfem::FiniteElementCollection class to avoid including the entire mfem finite element collection header file. -} // namespace mfem - -namespace serif::geometry::deformation::backend::mfem::detail -{ - /** - * @brief Shared mesh metadata; a negative node revision and empty material list mean those values are irrelevant. - * - * @todo This is a holdover from MeanFields bespoke reference tracking system, consider if we can remove this with SERiF's more robust - * dependencies graph - */ - struct ReferenceMeshSnapshot - { - const ::mfem::Mesh* identity{}; ///< The identity of the mesh, which is a pointer to the mfem::Mesh object representing the mesh. - long sequence{-1}; ///< The sequence number of the mesh, which is a long integer representing the version or revision of the mesh. - long nodes_sequence{-1}; ///< The sequence number of the nodes in the mesh, which is a long integer representing the version or revision of the nodes. - int vertices{}; ///< The number of vertices in the mesh, which is an integer representing the count of vertices in the mesh. - int elements{}; ///< The number of elements in the mesh, which is an integer representing the count of elements in the mesh. - std::vector element_attributes; ///< A vector of integers representing the attributes of the elements in the mesh. - std::vector boundary_attributes; ///< A vector of integers representing the attributes of the boundaries in the mesh. - }; - - /** - * @brief Reference space metadata; by that we mean that this struct contains the information needed to uniquely identify - * some finite element space so that we can determine if two finite element spaces are equivalent or not. Used for cacheing. - * Like I say in the todo, this is a hold over from the older prototype code MeanField which had this bespoke tracking - * system per module. I think that at this point these could likely be removed in favor of using MeanField's - * more robust and extensible dependencies graph system. - * - * @todo See the todo on ReferenceMeshSnapshot. Same points apply here. - */ - struct ReferenceSpaceSnapshot - { - const ::mfem::FiniteElementSpace* identity; ///< The identity of the finite element space, which is a pointer to the mfem::FiniteElementSpace object representing the space. - const ::mfem::Mesh* mesh; ///< The mesh associated with the finite element space, which is a pointer to the mfem::Mesh object representing the mesh. - const ::mfem::FiniteElementCollection* collection; ///< The finite element collection associated with the finite element space, which is a pointer to the mfem::FiniteElementCollection object representing the collection. - long sequence; ///< The sequence number of the finite element space, which is a long integer representing the version or revision of the space. - int size; ///< The size of the finite element space, which is an integer representing the number of degrees of freedom in the space. - int vector_dimension; ///< The vector dimension of the finite element space, which is an integer representing the dimension of the vector space. - int ordering; ///< The ordering of the finite element space, which is an integer representing the ordering of the degrees of freedom. - - - /** - * @brief Compare coefficient interpretations exactly. - */ - bool operator==(const ReferenceSpaceSnapshot&) const = default; - }; -} // namespace serif::geometry::deformation::backend::mfem::detail diff --git a/src/include/serif/geometry/deformation/composition.hpp b/src/include/serif/geometry/deformation/composition.hpp new file mode 100644 index 00000000..50e06c47 --- /dev/null +++ b/src/include/serif/geometry/deformation/composition.hpp @@ -0,0 +1,217 @@ +#pragma once + +#include +#include +#include +#include + +#include "serif/geometry/deformation/actions.hpp" + +namespace serif::geometry::deformation +{ + /** + * @brief Compose shape-to-surface and surface-to-volume maps with the complete derivative chain rule. + * @note Numerical factors are immutable and owned. Workspaces belong to callers and must not overlap public inputs or outputs; this class has no preparation state. + */ + template + requires std::same_as + class ComposedDeformation + { + public: + using InputSpace = typename SurfaceMap::InputSpace; + using OutputSpace = typename ExtensionMap::OutputSpace; + + static constexpr bool is_linear = is_linear_deformation && is_linear_deformation; + static constexpr bool uses_linear_actions = LinearDeformationMap && LinearDeformationMap; + + /** @brief Nonlinear actions retain the tangent and the two contributions to the shape-force derivative. */ + struct CurvatureWorkspace + { + std::vector surface_direction; + std::vector surface_dual; + std::vector parameter_pullback; + }; + + struct Workspace + { + DeformationWorkspace surface_workspace; + DeformationWorkspace extension_workspace; + std::vector surface_displacement; + std::conditional_t curvature; + }; + + /** @brief Keep solver coordinates in the surface model's native parameter space. */ + [[nodiscard]] std::size_t input_size() const noexcept + { + return m_surface->input_size(); + } + + /** @brief Expose the displacement field size consumed by the existing geometry mapper. */ + [[nodiscard]] std::size_t output_size() const noexcept + { + return m_extension->output_size(); + } + + /** @brief Reserve every intermediate once, outside repeated Newton and Krylov actions. */ + [[nodiscard]] Workspace make_workspace() const + { + Workspace workspace{ + .surface_workspace = make_deformation_workspace(*m_surface), + .extension_workspace = make_deformation_workspace(*m_extension), + .surface_displacement = std::vector(m_surface->output_size()), + .curvature = {}}; + if constexpr (!uses_linear_actions) + { + workspace.curvature = { + .surface_direction = std::vector(m_surface->output_size()), + .surface_dual = std::vector(m_surface->output_size()), + .parameter_pullback = std::vector(input_size())}; + } + return workspace; + } + + /** @brief A fixed composition uses the same primal action as its Jacobian. */ + [[nodiscard]] DeformationStatus apply(std::span input, std::span output, Workspace& workspace) const + requires uses_linear_actions + { + return build_displacement(input, output, workspace); + } + + /** @brief Reverse a fixed composition using one intermediate and no artificial base state. */ + [[nodiscard]] DeformationStatus apply_transpose(std::span dual, std::span output, Workspace& workspace) const + requires uses_linear_actions + { + if (const DeformationStatus status = detail::validate_disjoint_output(output, {dual}); !status) + { + return status; + } + if (const DeformationStatus status = m_extension->apply_transpose(dual, workspace.surface_displacement, workspace.extension_workspace); !status) + { + return status; + } + return m_surface->apply_transpose(workspace.surface_displacement, output, workspace.surface_workspace); + } + + /** @brief Move the surface first, then extend that same vector trace into the volume. */ + [[nodiscard]] DeformationStatus build_displacement(std::span parameters, std::span displacement, Workspace& workspace) const + { + if (const auto status = detail::validate_disjoint_output(displacement, {parameters}); !status) + { + return status; + } + if (const auto status = deformation::build_displacement(*m_surface, parameters, workspace.surface_displacement, workspace.surface_workspace); !status) + { + return status; + } + return deformation::build_displacement(*m_extension, workspace.surface_displacement, displacement, workspace.extension_workspace); + } + + /** @brief Evaluate the tangent of each nonlinear factor at the accepted surface displacement. */ + [[nodiscard]] DeformationStatus apply_jacobian(std::span parameters, std::span direction, std::span displacement_direction, Workspace& workspace) const + requires(!uses_linear_actions) + { + if (const auto status = detail::validate_disjoint_output(displacement_direction, {parameters, direction}); !status) + { + return status; + } + if (const auto status = deformation::build_displacement(*m_surface, parameters, workspace.surface_displacement, workspace.surface_workspace); !status) + { + return status; + } + if (const auto status = deformation::apply_jacobian(*m_surface, parameters, direction, workspace.curvature.surface_direction, workspace.surface_workspace); !status) + { + return status; + } + return deformation::apply_jacobian(*m_extension, workspace.surface_displacement, workspace.curvature.surface_direction, displacement_direction, workspace.extension_workspace); + } + + /** @brief Reverse both factors in the forward composition, preserving their work pairing. */ + [[nodiscard]] DeformationStatus apply_jacobian_transpose(std::span parameters, std::span displacement_dual, std::span parameter_dual, Workspace& workspace) const + requires(!uses_linear_actions) + { + if (const auto status = detail::validate_disjoint_output(parameter_dual, {parameters, displacement_dual}); !status) + { + return status; + } + if (const auto status = deformation::build_displacement(*m_surface, parameters, workspace.surface_displacement, workspace.surface_workspace); !status) + { + return status; + } + if (const auto status = deformation::apply_jacobian_transpose(*m_extension, workspace.surface_displacement, displacement_dual, workspace.curvature.surface_dual, workspace.extension_workspace); !status) + { + return status; + } + return deformation::apply_jacobian_transpose(*m_surface, parameters, workspace.curvature.surface_dual, parameter_dual, workspace.surface_workspace); + } + + /** @brief Retain both the extension curvature and the derivative of the surface transpose. */ + [[nodiscard]] DeformationStatus apply_pullback_derivative(std::span parameters, std::span direction, std::span displacement_dual, std::span parameter_dual, Workspace& workspace) const + requires(!uses_linear_actions) + { + if (const auto status = detail::validate_disjoint_output(parameter_dual, {parameters, direction, displacement_dual}); !status) + { + return status; + } + if (const auto status = deformation::build_displacement(*m_surface, parameters, workspace.surface_displacement, workspace.surface_workspace); !status) + { + return status; + } + if (const auto status = deformation::apply_jacobian(*m_surface, parameters, direction, workspace.curvature.surface_direction, workspace.surface_workspace); !status) + { + return status; + } + if (const auto status = deformation::apply_pullback_derivative(*m_extension, workspace.surface_displacement, workspace.curvature.surface_direction, displacement_dual, workspace.curvature.surface_dual, workspace.extension_workspace); !status) + { + return status; + } + if (const auto status = deformation::apply_jacobian_transpose(*m_surface, parameters, workspace.curvature.surface_dual, parameter_dual, workspace.surface_workspace); !status) + { + return status; + } + if (const auto status = deformation::apply_jacobian_transpose(*m_extension, workspace.surface_displacement, displacement_dual, workspace.curvature.surface_dual, workspace.extension_workspace); !status) + { + return status; + } + if (const auto status = deformation::apply_pullback_derivative(*m_surface, parameters, direction, workspace.curvature.surface_dual, workspace.curvature.parameter_pullback, workspace.surface_workspace); !status) + { + return status; + } + for (std::size_t parameter = 0; parameter < parameter_dual.size(); ++parameter) + { + parameter_dual[parameter] += workspace.curvature.parameter_pullback[parameter]; + } + return validate_deformation_values(parameter_dual, input_size()); + } + + /** @brief Expose immutable surface mathematics for reference metrics and diagnostics. */ + [[nodiscard]] const SurfaceMap& surface_map() const noexcept + { + return *m_surface; + } + + private: + /** @brief Only the checked composition factory may establish factor compatibility. */ + ComposedDeformation(std::shared_ptr surface, std::shared_ptr extension) : m_surface(std::move(surface)), m_extension(std::move(extension)) + { + } + + template + requires std::same_as + friend std::expected, DeformationFailure> compose_deformation(std::shared_ptr, std::shared_ptr); + + std::shared_ptr m_surface; + std::shared_ptr m_extension; + }; + + /** @brief Validate factor sizes before constructing an owning mathematical composition. */ + template + requires std::same_as + [[nodiscard]] std::expected, DeformationFailure> compose_deformation(std::shared_ptr surface, std::shared_ptr extension) + { + if (!surface || !extension || surface->output_size() != extension->input_size()) + { + return std::unexpected(DeformationFailure{.reason = DeformationFailureReason::invalid_size, .message = "Surface and volume extension must share the same displacement trace size."}); + } + return ComposedDeformation(std::move(surface), std::move(extension)); + } +} // namespace serif::geometry::deformation diff --git a/src/include/serif/geometry/deformation/computations.hpp b/src/include/serif/geometry/deformation/computations.hpp new file mode 100644 index 00000000..3f711cfa --- /dev/null +++ b/src/include/serif/geometry/deformation/computations.hpp @@ -0,0 +1,111 @@ +#pragma once + +#include +#include +#include "serif/geometry/deformation/composition.hpp" +#include "serif/geometry/deformation/linearization.hpp" +#include "serif/utils/dependencies/computation_graph.hpp" + +namespace serif::geometry::deformation { + namespace keys { + struct ReferenceDiscretization { + }; + + struct SurfaceDiscretization { + }; + + struct ReferenceCoordinates { + }; + + struct SurfaceNodes { + }; + + struct SurfaceModel { + }; + + struct InteriorModel { + }; + + struct VacuumModel { + }; + + struct ShapeCoefficients { + }; + + struct ReferenceSurface { + }; + + struct RadialInterpolation { + }; + + struct SurfaceMap { + }; + + struct VolumeExtension { + }; + + struct DomainMap { + }; + + struct SurfaceDisplacement { + }; + + struct VolumeDisplacement { + }; + + struct Linearization { + }; + + struct SurfaceMetric { + }; + } // namespace keys + + /** + * @brief Declare the numerical composition as ordinary dependency computations. + * @note Shared immutable factors, coefficient snapshots, and results belong to the dependency cache. Analytical actions remain independent functions. + */ + template + requires std::same_as + [[nodiscard]] auto make_deformation_computations() { + using namespace utils::dependencies; + using DomainMap = ComposedDeformation; + + const auto domain_map = make_computation>( + [](const std::shared_ptr& surface, const std::shared_ptr& extension) -> std::expected, DeformationFailure> { + auto result = compose_deformation(surface, extension); + if (!result) { + return std::unexpected(result.error()); + } + return std::make_shared(std::move(*result)); + }); + const auto surface_displacement = make_computation>( + [](const std::shared_ptr& surface, std::span parameters) -> std::expected, DeformationFailure> { + if (!surface) { + return std::unexpected(DeformationFailure{.reason = DeformationFailureReason::invalid_configuration, .message = "Shape evaluation requires an owned surface map."}); + } + return evaluate_displacement(*surface, parameters); + }); + const auto volume_displacement = make_computation>( + [](const std::shared_ptr& extension, const std::vector& surface) -> std::expected, DeformationFailure> { + if (!extension) { + return std::unexpected(DeformationFailure{.reason = DeformationFailureReason::invalid_configuration, .message = "Volume deformation requires an owned extension map."}); + } + return evaluate_displacement(*extension, surface); + }); + const auto linearization = make_computation>( + [](const std::shared_ptr& map, std::span parameters, const std::vector&) { + return DeformationLinearization{.map = map, .parameters = {parameters.begin(), parameters.end()}}; + }); + return std::tuple{domain_map, surface_displacement, volume_displacement, linearization}; + } + + /** @brief Compose independently compiled mathematical factors without introducing a prepared-operator facade. */ + template + requires std::same_as + [[nodiscard]] auto make_deformation_graph() { + return std::apply([](auto... computations) { + return utils::dependencies::make_computation_graph(std::move(computations)...); + }, + make_deformation_computations()); + } +} // namespace serif::geometry::deformation diff --git a/src/include/serif/geometry/deformation/config_options.hpp b/src/include/serif/geometry/deformation/config_options.hpp deleted file mode 100644 index e69de29b..00000000 diff --git a/src/include/serif/geometry/deformation/deformation.hpp b/src/include/serif/geometry/deformation/deformation.hpp index e69de29b..5338e0c2 100644 --- a/src/include/serif/geometry/deformation/deformation.hpp +++ b/src/include/serif/geometry/deformation/deformation.hpp @@ -0,0 +1,8 @@ +#pragma once + +#include "serif/geometry/deformation/actions.hpp" +#include "serif/geometry/deformation/computations.hpp" +#include "serif/geometry/deformation/geometry_step.hpp" +#include "serif/geometry/deformation/models.hpp" +#include "serif/geometry/deformation/radial_extension.hpp" +#include "serif/geometry/deformation/surface_embedding.hpp" diff --git a/src/include/serif/geometry/deformation/descriptors.hpp b/src/include/serif/geometry/deformation/descriptors.hpp deleted file mode 100644 index e69de29b..00000000 diff --git a/src/include/serif/geometry/deformation/detail/action_storage.hpp b/src/include/serif/geometry/deformation/detail/action_storage.hpp new file mode 100644 index 00000000..fac61a1c --- /dev/null +++ b/src/include/serif/geometry/deformation/detail/action_storage.hpp @@ -0,0 +1,86 @@ +#pragma once + +#include +#include +#include +#include +#include +#include "serif/geometry/deformation/failure.hpp" + +namespace serif::geometry::deformation +{ + /** @brief Reject incompatible lengths and non-finite values before evaluating a deformation action. */ + [[nodiscard]] inline DeformationStatus validate_deformation_values(std::span values, std::size_t expected_size) + { + if (values.size() != expected_size) + { + return std::unexpected(DeformationFailure{.reason = DeformationFailureReason::invalid_size, .message = "Deformation vector size does not match its map."}); + } + if (!std::ranges::all_of(values, [](double value) + { + return std::isfinite(value); + })) + { + return std::unexpected(DeformationFailure{.reason = DeformationFailureReason::non_finite, .message = "Deformation coefficients, directions, forces and action results must be finite."}); + } + return {}; + } +} // namespace serif::geometry::deformation + +namespace serif::geometry::deformation::detail +{ + /** + * @brief Preserve every action input until the last use, including nonlinear transpose derivatives. + * The list contains only the fixed number of mathematical arguments to an action; checking + * address intervals requires no allocation and never scans numerical coefficients. + */ + [[nodiscard]] inline DeformationStatus validate_disjoint_output(std::span output, std::initializer_list> inputs) + { + if (output.empty()) + { + return {}; + } + const std::less before; + const double* output_end = output.data() + output.size(); + for (const std::span input : inputs) + { + if (input.empty()) + { + continue; + } + if (before(input.data(), output_end) && before(output.data(), input.data() + input.size())) + { + return std::unexpected(DeformationFailure{.reason = DeformationFailureReason::invalid_configuration, .message = "Deformation action output must not overlap its parameters, directions, or forces."}); + } + } + return {}; + } + + /** @brief Pair each action input with the dimension of its mathematical space. */ + struct ActionInput + { + std::span values; + std::size_t size; + }; + + /** @brief Validate caller buffers uniformly, retaining disjoint inputs until an action has finished. */ + [[nodiscard]] inline DeformationStatus validate_action_arguments(std::span output, std::size_t output_size, std::initializer_list inputs) + { + if (output.size() != output_size) + { + return std::unexpected(DeformationFailure{.reason = DeformationFailureReason::invalid_size, .message = "Deformation output size does not match its map."}); + } + for (const ActionInput input : inputs) + { + if (const DeformationStatus status = validate_disjoint_output(output, {input.values}); !status) + { + return status; + } + if (const DeformationStatus status = validate_deformation_values(input.values, input.size); !status) + { + return status; + } + } + return {}; + } +} // namespace serif::geometry::deformation::detail diff --git a/src/include/serif/geometry/deformation/detail/determinant_polynomial.hpp b/src/include/serif/geometry/deformation/detail/determinant_polynomial.hpp new file mode 100644 index 00000000..6a11914f --- /dev/null +++ b/src/include/serif/geometry/deformation/detail/determinant_polynomial.hpp @@ -0,0 +1,21 @@ +#pragma once + +#include + +namespace serif::geometry::deformation::detail { + using DeterminantPolynomial = std::array; + using JacobianEntries = std::array, 3>; + + /** @brief Expand det(base + t maximum_step direction) minus the floor for dimensions one through three. */ + [[nodiscard]] DeterminantPolynomial determinant_polynomial(const JacobianEntries& base, const JacobianEntries& direction, int dimension, double maximum_step, double determinant_floor) noexcept; + + /** @brief Evaluate the cubic with fused operations, retaining small values near repeated roots. */ + [[nodiscard]] double evaluate_polynomial(const DeterminantPolynomial& coefficients, double parameter) noexcept; + + /** + * @brief Find the first loss of strict positivity in (0, 1], including a tangent zero. + * @pre Coefficients are finite and the constant coefficient is strictly positive. + * @return The boundary parameter, or infinity when the entire interval remains admissible. + */ + [[nodiscard]] double first_boundary_parameter(const DeterminantPolynomial& coefficients) noexcept; +} diff --git a/src/include/serif/geometry/deformation/detail/harmonics.hpp b/src/include/serif/geometry/deformation/detail/harmonics.hpp new file mode 100644 index 00000000..670431a4 --- /dev/null +++ b/src/include/serif/geometry/deformation/detail/harmonics.hpp @@ -0,0 +1,20 @@ +#pragma once + +#include +#include +#include +#include +#include + +#include "serif/geometry/deformation/surface_embedding.hpp" + +namespace serif::geometry::deformation::detail { + /** @brief Validate harmonic cutoffs before allocating a basis or enumerating modes. */ + [[nodiscard]] std::expected, DeformationFailure> make_harmonic_modes(const SphericalHarmonicSurface& model, std::size_t node_count); + + /** @brief Use normalized Legendre recurrence to avoid factorial growth while preserving the Condon–Shortley phase. */ + [[nodiscard]] std::expected, DeformationFailure> sample_harmonic_basis(const SphericalHarmonicSurface& model, std::span modes, std::span> directions); + + /** @brief Check sampled rank on a private copy, leaving physical coefficients and basis normalization unchanged. */ + [[nodiscard]] DeformationStatus check_independent_harmonics(std::span basis, std::size_t node_count, std::size_t mode_count); +} diff --git a/src/include/serif/geometry/deformation/failure.hpp b/src/include/serif/geometry/deformation/failure.hpp new file mode 100644 index 00000000..3705d2ee --- /dev/null +++ b/src/include/serif/geometry/deformation/failure.hpp @@ -0,0 +1,25 @@ +#pragma once + +#include +#include + +namespace serif::geometry::deformation { + enum class DeformationFailureReason { + invalid_configuration, + incompatible_discretization, + invalid_size, + non_finite, + unresolved_surface, + invalid_geometry + }; + + /** @brief Describe a rejected geometry calculation; the dependency cache publishes results only after success. */ + struct DeformationFailure { + DeformationFailureReason reason; + std::string message; + int element{-1}; + int point{-1}; + }; + + using DeformationStatus = std::expected; +} diff --git a/src/include/serif/geometry/deformation/finite_elements.hpp b/src/include/serif/geometry/deformation/finite_elements.hpp new file mode 100644 index 00000000..d7acfad0 --- /dev/null +++ b/src/include/serif/geometry/deformation/finite_elements.hpp @@ -0,0 +1,9 @@ +#pragma once + +#include "serif/backend/mfem/geometry/deformation/computations.hpp" + +namespace serif::geometry::deformation { + using backend::mfem::make_deformation_inputs; + using backend::mfem::make_deformation_graph; + using backend::mfem::make_deformation_computations; +} diff --git a/src/include/serif/geometry/deformation/geometry_step.hpp b/src/include/serif/geometry/deformation/geometry_step.hpp new file mode 100644 index 00000000..b95338eb --- /dev/null +++ b/src/include/serif/geometry/deformation/geometry_step.hpp @@ -0,0 +1,33 @@ +#pragma once + +#include +#include + +namespace serif::geometry::deformation { + /** @brief Numerical controls for keeping a straight volume-displacement step above a Jacobian determinant floor. */ + struct GeometryStepOptions { + double maximum_step{1.0}; + double determinant_floor{0.0}; + double fraction_to_boundary{0.9}; + }; + + /** + * @brief Geometry-only bound at the supplied sample points, with a location for the first obstruction. + * + * The determinant is that of the physical map relative to the reference mesh. This does not bound + * unsampled points, thermodynamic admissibility, or a nonlinear trajectory in shape coefficients. + * A boundary at maximum_step still applies fraction_to_boundary. An unrestricted step does not. + */ + struct GeometryStepLimit { + double step{0.0}; + double boundary_step{0.0}; + double minimum_determinant_at_accepted{std::numeric_limits::infinity()}; + double minimum_determinant_at_maximum{std::numeric_limits::infinity()}; + double limiting_point_determinant_at_step{0.0}; + std::size_t sampled_point_count{0}; + bool limited_by_geometry{false}; + int limiting_element{-1}; + int limiting_rule{-1}; + int limiting_point{-1}; + }; +} diff --git a/src/include/serif/geometry/deformation/linearization.hpp b/src/include/serif/geometry/deformation/linearization.hpp new file mode 100644 index 00000000..6487925f --- /dev/null +++ b/src/include/serif/geometry/deformation/linearization.hpp @@ -0,0 +1,52 @@ +#pragma once + +#include +#include +#include "serif/geometry/deformation/actions.hpp" + +namespace serif::geometry::deformation { + /** @brief Freeze the shape coordinates associated with a derivative, independently of mutable caller inputs. */ + template + struct DeformationLinearization { + std::shared_ptr map; + std::vector parameters; + }; + + /** @brief Allocate a result for graph publication while leaving repeated low-level actions allocation-free. */ + template + [[nodiscard]] std::expected, DeformationFailure> evaluate_displacement(const Map& map, std::span parameters) { + std::vector result(map.output_size()); + auto workspace = make_deformation_workspace(map); + if (const auto status = build_displacement(map, parameters, result, workspace); !status) { + return std::unexpected(status.error()); + } + return result; + } + + /** @brief Evaluate a tangent about the cache's accepted shape, not the caller's next trial state. */ + template + [[nodiscard]] DeformationStatus apply_jacobian(const DeformationLinearization& state, std::span direction, std::span output, DeformationWorkspace& workspace) { + if (!state.map) { + return std::unexpected(DeformationFailure{.reason = DeformationFailureReason::invalid_configuration, .message = "A deformation linearization must own a valid mathematical map."}); + } + return apply_jacobian(*state.map, state.parameters, direction, output, workspace); + } + + /** @brief Pull displacement forces back to the accepted independent shape coordinates. */ + template + [[nodiscard]] DeformationStatus apply_jacobian_transpose(const DeformationLinearization& state, std::span dual, std::span output, DeformationWorkspace& workspace) { + if (!state.map) { + return std::unexpected(DeformationFailure{.reason = DeformationFailureReason::invalid_configuration, .message = "A deformation linearization must own a valid mathematical map."}); + } + return apply_jacobian_transpose(*state.map, state.parameters, dual, output, workspace); + } + + /** @brief Keep the geometric contribution to the Newton force derivative available for nonlinear models. */ + template + [[nodiscard]] DeformationStatus apply_pullback_derivative(const DeformationLinearization& state, std::span direction, std::span dual, std::span output, DeformationWorkspace& workspace) { + if (!state.map) { + return std::unexpected(DeformationFailure{.reason = DeformationFailureReason::invalid_configuration, .message = "A deformation linearization must own a valid mathematical map."}); + } + return apply_pullback_derivative(*state.map, state.parameters, direction, dual, output, workspace); + } +} diff --git a/src/include/serif/geometry/deformation/models.hpp b/src/include/serif/geometry/deformation/models.hpp new file mode 100644 index 00000000..7ded1b27 --- /dev/null +++ b/src/include/serif/geometry/deformation/models.hpp @@ -0,0 +1,49 @@ +#pragma once + +#include +#include + +namespace serif::geometry::deformation { + /** @brief Give each surface node an independent radial displacement, measured in physical length. */ + struct NodalRadialSurface { + std::array reference_center{}; + + /// Reference centers are physical inputs to dependency tracking. + [[nodiscard]] bool operator==(const NodalRadialSurface&) const = default; + }; + + /** @brief Represent radial displacement lengths with normalized real spherical harmonics on the reference surface. */ + struct SphericalHarmonicSurface { + std::array reference_center{}; + int maximum_degree{4}; + std::optional maximum_azimuthal_order{}; + + /// A changed center or harmonic selection changes the physical representation. + [[nodiscard]] bool operator==(const SphericalHarmonicSurface&) const = default; + }; + + /** @brief Extend surface motion toward a fixed stellar center using a power of logical reference radius. */ + struct PowerLawRadialExtension { + double radial_power{2.0}; + + /// The radial power changes the interior response even for unchanged surface motion. + [[nodiscard]] bool operator==(const PowerLawRadialExtension&) const = default; + }; + + /** @brief Reduce surface motion linearly to zero at reference infinity in logical reference radius. */ + struct FixedInfinityRadialExtension { + /// All fixed-infinity prescriptions have the same physical behavior. + [[nodiscard]] bool operator==(const FixedInfinityRadialExtension&) const = default; + }; + + /** @brief State how the stellar surface moves and how that motion extends into the two surrounding regions. */ + template + struct DomainDeformation { + SurfaceT surface; + InteriorT stellar_interior; + VacuumT vacuum; + + /// Compare physical choices without attaching preparation state to the model. + [[nodiscard]] bool operator==(const DomainDeformation&) const = default; + }; +} diff --git a/src/include/serif/geometry/deformation/radial_extension.hpp b/src/include/serif/geometry/deformation/radial_extension.hpp new file mode 100644 index 00000000..a819e8fb --- /dev/null +++ b/src/include/serif/geometry/deformation/radial_extension.hpp @@ -0,0 +1,77 @@ +#pragma once + +#include +#include +#include +#include "serif/geometry/deformation/radial_interpolation.hpp" +#include "serif/geometry/deformation/models.hpp" +#include "serif/geometry/deformation/spaces.hpp" +#include "serif/geometry/deformation/workspace.hpp" + +namespace serif::geometry::deformation +{ + /** + * @brief Extend one interleaved surface vector into the MFEM component-block volume layout. + * Stellar ownership wins at shared nodes. Both actions use the same scalar interpolation, + * so their pairing is Euclidean virtual work without an additional surface-area weight. + */ + class RadialExtension + { + public: + using Workspace = EmptyDeformationWorkspace; + using InputSpace = spaces::SurfaceDisplacement; + using OutputSpace = spaces::VolumeDisplacement; + static constexpr bool is_linear = true; + + /** @brief Expose shared interpolation read-only for diagnostics and independent geometry checks. */ + [[nodiscard]] const std::shared_ptr& interpolation() const noexcept + { + return m_interpolation; + } + + /** @brief Expose the ownership-selected radial response without permitting mutation. */ + [[nodiscard]] std::span radial_weights() const noexcept + { + return m_radial_weights; + } + + /** @brief Return the interleaved Cartesian surface-vector size. */ + [[nodiscard]] std::size_t input_size() const noexcept + { + return 3 * m_interpolation->surface_node_count(); + } + + /** @brief Return the Cartesian component-block volume-vector size. */ + [[nodiscard]] std::size_t output_size() const noexcept + { + return 3 * m_radial_weights.size(); + } + + /** + * @brief Apply the fused radial extension without allocating action scratch. + * Input and output must not overlap. Discard output if an action reports failure. + */ + [[nodiscard]] DeformationStatus apply(std::span surface_displacement, std::span volume_displacement, Workspace = {}) const; + + /** @brief Transfer volume forces to surface forces using exactly the forward ownership. */ + [[nodiscard]] DeformationStatus apply_transpose(std::span volume_force, std::span surface_force, Workspace = {}) const; + + private: + /** @brief Store already validated data without performing hidden numerical preparation. */ + RadialExtension(std::shared_ptr interpolation, std::vector radial_weights); + + friend std::expected make_radial_extension(std::shared_ptr, PowerLawRadialExtension, FixedInfinityRadialExtension); + + std::shared_ptr m_interpolation; + std::vector m_radial_weights; + }; + + /** + * @brief Compile regional weights while sharing the immutable surface interpolation. + * The interpolation factory has already established every reference-table invariant. + */ + [[nodiscard]] std::expected make_radial_extension( + std::shared_ptr interpolation, + PowerLawRadialExtension interior = {}, + FixedInfinityRadialExtension vacuum = {}); +} // namespace serif::geometry::deformation diff --git a/src/include/serif/geometry/deformation/radial_interpolation.hpp b/src/include/serif/geometry/deformation/radial_interpolation.hpp new file mode 100644 index 00000000..173ecb7c --- /dev/null +++ b/src/include/serif/geometry/deformation/radial_interpolation.hpp @@ -0,0 +1,52 @@ +#pragma once + +#include +#include +#include "serif/geometry/deformation/failure.hpp" + +namespace serif::geometry::deformation +{ + /** @brief Own the surface interpolation shared by the stellar and vacuum extensions. */ + struct RadialInterpolationData + { + std::vector surface_scalar_dofs; + std::vector row_offsets; + std::vector surface_columns; + std::vector interpolation_weights; + std::vector logical_radii; + std::vector stellar_support; + std::vector vacuum_support; + double stellar_radius{}; + double infinity_radius{}; + }; + + /** @brief A checked, immutable interpolation table; changing a radial law never revalidates its topology. */ + class RadialInterpolation + { + public: + /** @brief Inspect the canonical nodes, interpolation rows and logical shell geometry without mutating them. */ + [[nodiscard]] const RadialInterpolationData& data() const noexcept { + return m_data; + } + + /** @brief Return the volume scalar-node count, independent of vector layout. */ + [[nodiscard]] std::size_t scalar_dof_count() const noexcept { + return m_data.logical_radii.size(); + } + + /** @brief Return the number of canonical stellar surface nodes. */ + [[nodiscard]] std::size_t surface_node_count() const noexcept { + return m_data.surface_scalar_dofs.size(); + } + + private: + /** @brief Accept ownership only after the explicit factory has checked the table. */ + explicit RadialInterpolation(RadialInterpolationData data) : m_data(std::move(data)) {} + + friend std::expected make_radial_interpolation(RadialInterpolationData); + RadialInterpolationData m_data; + }; + + /** @brief Check imported or backend-generated reference rows once before sharing them with radial laws. */ + [[nodiscard]] std::expected make_radial_interpolation(RadialInterpolationData data); +} // namespace serif::geometry::deformation diff --git a/src/include/serif/geometry/deformation/reference_surface.hpp b/src/include/serif/geometry/deformation/reference_surface.hpp new file mode 100644 index 00000000..c9e1a593 --- /dev/null +++ b/src/include/serif/geometry/deformation/reference_surface.hpp @@ -0,0 +1,12 @@ +#pragma once + +#include +#include + +namespace serif::geometry::deformation { + /** @brief Own fixed surface positions and their canonical ascending scalar-node indices. */ + struct ReferenceSurface { + std::vector> positions; + std::vector scalar_dofs; + }; +} diff --git a/src/include/serif/geometry/deformation/spaces.hpp b/src/include/serif/geometry/deformation/spaces.hpp new file mode 100644 index 00000000..cf1a84a4 --- /dev/null +++ b/src/include/serif/geometry/deformation/spaces.hpp @@ -0,0 +1,12 @@ +#pragma once + +namespace serif::geometry::deformation::spaces { + /** @brief Independent coordinates belonging to the selected surface representation. */ + struct ShapeParameters {}; + + /** @brief Cartesian surface traces with components interleaved at each canonical node. */ + struct SurfaceDisplacement {}; + + /** @brief Displacement coefficients in the declared volume discretization. */ + struct VolumeDisplacement {}; +} diff --git a/src/include/serif/geometry/deformation/surface_embedding.hpp b/src/include/serif/geometry/deformation/surface_embedding.hpp new file mode 100644 index 00000000..efbd1fbd --- /dev/null +++ b/src/include/serif/geometry/deformation/surface_embedding.hpp @@ -0,0 +1,90 @@ +#pragma once + +#include +#include +#include +#include +#include + +#include "serif/geometry/deformation/failure.hpp" +#include "serif/geometry/deformation/models.hpp" +#include "serif/geometry/deformation/spaces.hpp" +#include "serif/geometry/deformation/workspace.hpp" + +namespace serif::geometry::deformation +{ + /** @brief Identify a real harmonic, with negative orders using sine and positive orders using cosine. */ + struct SphericalHarmonicMode + { + int degree; + int order; + + /// Mode identity includes the sign of the azimuthal order. + [[nodiscard]] bool operator==(const SphericalHarmonicMode&) const = default; + }; + + /** + * @brief Own a linear conversion from radial displacement lengths to Cartesian surface displacements. + * @note Directions and harmonics remain fixed on the reference surface. Finite-element interpolation of the resulting nodal vectors determines the between-node geometry. + */ + class SurfaceEmbedding + { + public: + using Workspace = EmptyDeformationWorkspace; + using InputSpace = spaces::ShapeParameters; + using OutputSpace = spaces::SurfaceDisplacement; + static constexpr bool is_linear = true; + + /// Nodal coordinates are independent lengths; harmonic coordinates retain their physical normalization. + [[nodiscard]] std::size_t input_size() const noexcept; + + /// Surface vectors interleave the three Cartesian components at each node. + [[nodiscard]] std::size_t output_size() const noexcept; + + /// Expose the ordered surface-node count for independent metric and backend calculations. + [[nodiscard]] std::size_t surface_node_count() const noexcept; + + /// An empty harmonic basis represents the identity between parameters and nodal radial lengths. + [[nodiscard]] bool has_nodal_parameters() const noexcept; + + /// Reference directions determine radial motion and remain immutable when coefficients change. + [[nodiscard]] std::span> radial_directions() const noexcept; + + /// Reference radii retain physical lengths independently of the logical radius used by volume extensions. + [[nodiscard]] std::span reference_radii() const noexcept; + + /// Preserve the exact ordered reference points so boundary metrics can reject an unrelated trace. + [[nodiscard]] std::span> reference_positions() const noexcept; + + /// Harmonic coordinates use increasing degree, then increasing signed azimuthal order. + [[nodiscard]] std::span modes() const noexcept; + + /// Return the physical row-major harmonic basis; nodal parameterizations have an implicit identity basis. + [[nodiscard]] std::span nodal_basis() const noexcept; + + /** @brief Apply radial motion without allocation. Input and output storage must not overlap; rejected output must be discarded. */ + [[nodiscard]] DeformationStatus apply(std::span parameters, std::span surface_displacement, Workspace = {}) const; + + /** @brief Transfer Cartesian surface forces to shape coordinates while preserving virtual work, without allocation or metric weighting. */ + [[nodiscard]] DeformationStatus apply_transpose(std::span surface_force, std::span parameter_force, Workspace = {}) const; + + private: + /// Factories validate all reference geometry before publishing an embedding. + SurfaceEmbedding() = default; + + friend std::expected build_surface_embedding(const NodalRadialSurface&, std::span>); + friend std::expected build_surface_embedding(const SphericalHarmonicSurface&, std::span>); + + std::vector> radial_directions_; + std::vector reference_radii_; + std::vector> reference_positions_; + std::vector modes_; + std::vector nodal_basis_; + }; + + /** @brief Compile fixed radial directions from ordered physical reference positions without retaining backend objects. */ + [[nodiscard]] std::expected build_surface_embedding(const NodalRadialSurface& model, std::span> positions); + + /** @brief Sample physically normalized real harmonics and reject unresolved modes without renormalizing or dropping coordinates. */ + [[nodiscard]] std::expected build_surface_embedding(const SphericalHarmonicSurface& model, std::span> positions); +} // namespace serif::geometry::deformation diff --git a/src/include/serif/geometry/deformation/vacuum_extension.hpp b/src/include/serif/geometry/deformation/vacuum_extension.hpp deleted file mode 100644 index e69de29b..00000000 diff --git a/src/include/serif/geometry/deformation/workspace.hpp b/src/include/serif/geometry/deformation/workspace.hpp new file mode 100644 index 00000000..3408565e --- /dev/null +++ b/src/include/serif/geometry/deformation/workspace.hpp @@ -0,0 +1,12 @@ +#pragma once + +namespace serif::geometry::deformation +{ + /** @brief Stateless maps require no numerical scratch, while retaining the same explicit action interface. */ + struct EmptyDeformationWorkspace + { + }; + + template + using DeformationWorkspace = typename Map::Workspace; +} // namespace serif::geometry::deformation diff --git a/src/lib/serif/geometry/deformation/backend/mfem/radial_locator.cpp b/src/lib/serif/geometry/deformation/backend/mfem/radial_locator.cpp deleted file mode 100644 index c3c4128b..00000000 --- a/src/lib/serif/geometry/deformation/backend/mfem/radial_locator.cpp +++ /dev/null @@ -1,126 +0,0 @@ -#include -#include -#include -#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 boundary_logical_radius(::mfem::Mesh& mesh, const int attribute) { - double minimum = std::numeric_limits::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 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::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::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; - } -} diff --git a/src/lib/serif/geometry/deformation/determinant_polynomial.cpp b/src/lib/serif/geometry/deformation/determinant_polynomial.cpp new file mode 100644 index 00000000..cbfcf5e2 --- /dev/null +++ b/src/lib/serif/geometry/deformation/determinant_polynomial.cpp @@ -0,0 +1,157 @@ +#include "serif/geometry/deformation/detail/determinant_polynomial.hpp" + +#include +#include +#include +#include +#include + +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(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::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& 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::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& 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::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(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::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::epsilon() * coefficient_scale; + + std::array critical_points{}; + const int critical_point_count = derivative_critical_points(coefficients, degree, critical_points); + std::array 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::infinity(); + } +} diff --git a/src/lib/serif/geometry/deformation/harmonics.cpp b/src/lib/serif/geometry/deformation/harmonics.cpp new file mode 100644 index 00000000..2cf188e3 --- /dev/null +++ b/src/lib/serif/geometry/deformation/harmonics.cpp @@ -0,0 +1,124 @@ +#include "serif/geometry/deformation/detail/harmonics.hpp" + +#include +#include +#include +#include +#include + +namespace serif::geometry::deformation::detail { + namespace { + /** @brief Keep accumulation in extended precision when checking sampled independence. */ + [[nodiscard]] double column_product(std::span 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(basis[row * columns + first]) * basis[row * columns + second]; + } + return static_cast(result); + } + + /** @brief Build normalized associated Legendre functions without factorial overflow. */ + void evaluate_legendre(const std::array& direction, int maximum_degree, int maximum_order, std::span legendre) { + const std::size_t stride = static_cast(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(degree) * stride + static_cast(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, 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::max()) { + return std::unexpected(DeformationFailure{.reason = DeformationFailureReason::invalid_configuration, .message = "The requested harmonic mode count exceeds supported vector indexing."}); + } + if (static_cast(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 modes; + modes.reserve(static_cast(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, DeformationFailure> sample_harmonic_basis(const SphericalHarmonicSurface& model, std::span modes, std::span> directions) { + if (modes.empty() || directions.size() > std::vector{}.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(maximum_order) + 1; + std::vector legendre((static_cast(model.maximum_degree) + 1) * stride, 0.0); + std::vector 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(mode.degree) * stride + static_cast(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 basis, std::size_t node_count, std::size_t mode_count) { + std::vector 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 {}; + } +} diff --git a/src/lib/serif/geometry/deformation/radial_extension.cpp b/src/lib/serif/geometry/deformation/radial_extension.cpp new file mode 100644 index 00000000..f55bbffe --- /dev/null +++ b/src/lib/serif/geometry/deformation/radial_extension.cpp @@ -0,0 +1,104 @@ +#include +#include +#include +#include +#include "serif/geometry/deformation/radial_extension.hpp" +#include "serif/geometry/deformation/detail/action_storage.hpp" + +namespace serif::geometry::deformation +{ + RadialExtension::RadialExtension(std::shared_ptr interpolation, std::vector radial_weights) + : m_interpolation(std::move(interpolation)), m_radial_weights(std::move(radial_weights)) + { + } + + DeformationStatus RadialExtension::apply(std::span surface_displacement, std::span 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 volume_force, std::span 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 make_radial_extension(std::shared_ptr 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(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::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 diff --git a/src/lib/serif/geometry/deformation/radial_interpolation.cpp b/src/lib/serif/geometry/deformation/radial_interpolation.cpp new file mode 100644 index 00000000..626df3f4 --- /dev/null +++ b/src/lib/serif/geometry/deformation/radial_interpolation.cpp @@ -0,0 +1,68 @@ +#include +#include +#include "serif/geometry/deformation/radial_interpolation.hpp" + +namespace serif::geometry::deformation +{ + std::expected 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(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(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::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(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::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 diff --git a/src/lib/serif/geometry/deformation/surface_embedding.cpp b/src/lib/serif/geometry/deformation/surface_embedding.cpp new file mode 100644 index 00000000..09651c57 --- /dev/null +++ b/src/lib/serif/geometry/deformation/surface_embedding.cpp @@ -0,0 +1,210 @@ +#include "serif/geometry/deformation/surface_embedding.hpp" + +#include +#include +#include +#include + +#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 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& direction, std::span force) + { + long double result = 0.0L; + for (std::size_t component = 0; component < 3; ++component) + { + result += static_cast(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> SurfaceEmbedding::radial_directions() const noexcept + { + return radial_directions_; + } + + std::span SurfaceEmbedding::reference_radii() const noexcept + { + return reference_radii_; + } + + std::span> SurfaceEmbedding::reference_positions() const noexcept + { + return reference_positions_; + } + + std::span SurfaceEmbedding::modes() const noexcept + { + return modes_; + } + + std::span SurfaceEmbedding::nodal_basis() const noexcept + { + return nodal_basis_; + } + + DeformationStatus SurfaceEmbedding::apply(std::span parameters, std::span 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(nodal_basis_[node * input_size() + mode]) * parameters[mode]; + } + } + const double radial_amplitude = static_cast(amplitude); + if (!std::isfinite(radial_amplitude)) + { + return action_overflow(); + } + for (std::size_t component = 0; component < 3; ++component) + { + const double displacement = static_cast(static_cast(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 surface_force, std::span 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(nodal_basis_[node * input_size() + parameter]) * radial_work(radial_directions_[node], surface_force.subspan(3 * node, 3)); + } + } + const double parameter_work = static_cast(work); + if (!std::isfinite(parameter_work)) + { + return action_overflow(); + } + parameter_force[parameter] = parameter_work; + } + return {}; + } + + std::expected build_surface_embedding(const NodalRadialSurface& model, std::span> positions) + { + if (positions.empty() || positions.size() > std::numeric_limits::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& position : positions) + { + std::array 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 build_surface_embedding(const SphericalHarmonicSurface& model, std::span> 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