From 8f2517d9b72e568efd9f209e5b648f4a44452b71 Mon Sep 17 00:00:00 2001 From: Michael Tupek Date: Thu, 27 Aug 2026 09:59:30 -0700 Subject: [PATCH 1/3] Add a capability to looad 1D tabulated data and have it be approximated by pchip splines (which respects local data bounds, total variation). Demo how to use within a user defined solid mechanics integral. --- src/smith/numerics/CMakeLists.txt | 1 + src/smith/numerics/pchip.hpp | 194 ++++++++++++++++++ src/smith/numerics/tests/CMakeLists.txt | 1 + src/smith/numerics/tests/test_pchip.cpp | 57 +++++ src/smith/physics/tests/CMakeLists.txt | 1 + .../tests/test_tabulated_solid_properties.cpp | 123 +++++++++++ 6 files changed, 377 insertions(+) create mode 100644 src/smith/numerics/pchip.hpp create mode 100644 src/smith/numerics/tests/test_pchip.cpp create mode 100644 src/smith/physics/tests/test_tabulated_solid_properties.cpp diff --git a/src/smith/numerics/CMakeLists.txt b/src/smith/numerics/CMakeLists.txt index cda8cb0bd6..69b8b392e6 100644 --- a/src/smith/numerics/CMakeLists.txt +++ b/src/smith/numerics/CMakeLists.txt @@ -9,6 +9,7 @@ set(numerics_headers equation_solver.hpp nonlinear_convergence.hpp odes.hpp + pchip.hpp solver_config.hpp solver_with_preconditioner.hpp stdfunction_operator.hpp diff --git a/src/smith/numerics/pchip.hpp b/src/smith/numerics/pchip.hpp new file mode 100644 index 0000000000..03b2307457 --- /dev/null +++ b/src/smith/numerics/pchip.hpp @@ -0,0 +1,194 @@ +// Copyright (c) Lawrence Livermore National Security, LLC and +// other Smith Project Developers. See the top-level LICENSE file for +// details. +// +// SPDX-License-Identifier: (BSD-3-Clause) + +#pragma once + +#include +#include +#include +#include +#include + +#include "axom/core/Array.hpp" +#include "smith/numerics/functional/dual.hpp" +#include "smith/numerics/functional/tensor.hpp" + +namespace smith { + +/** + * @brief Non-owning evaluator for shape-preserving piecewise cubic Hermite data. + * + * Each interval uses cubic Bernstein control points. Values outside + * the tabulated coordinate range clamp to endpoint values. + * + * @note Referenced PchipData must outlive every use of this view. + */ +class PchipView { + public: + /** + * @brief Evaluate interpolated value. + */ + template + SMITH_HOST_DEVICE auto operator()(T input) const + { + const double input_value = get_value(input); + if (input_value < coordinates_[0]) { + return control_points_[0][0] + 0.0 * input; + } + if (input_value > coordinates_[num_points_ - 1]) { + return control_points_[num_points_ - 2][3] + 0.0 * input; + } + + axom::IndexType lower = 0; + axom::IndexType upper = num_points_ - 1; + while (upper - lower > 1) { + const axom::IndexType middle = lower + (upper - lower) / 2; + if (input_value < coordinates_[middle]) { + upper = middle; + } else { + lower = middle; + } + } + + const auto local_coordinate = (input - coordinates_[lower]) / (coordinates_[lower + 1] - coordinates_[lower]); + const auto complement = 1.0 - local_coordinate; + const auto first_linear = complement * control_points_[lower][0] + local_coordinate * control_points_[lower][1]; + const auto second_linear = complement * control_points_[lower][1] + local_coordinate * control_points_[lower][2]; + const auto third_linear = complement * control_points_[lower][2] + local_coordinate * control_points_[lower][3]; + const auto first_quadratic = complement * first_linear + local_coordinate * second_linear; + const auto second_quadratic = complement * second_linear + local_coordinate * third_linear; + return complement * first_quadratic + local_coordinate * second_quadratic; + } + + private: + friend class PchipData; + + PchipView(const double* coordinates, const tensor* control_points, axom::IndexType num_points) + : coordinates_(coordinates), control_points_(control_points), num_points_(num_points) + { + } + + const double* coordinates_; + const tensor* control_points_; + axom::IndexType num_points_; +}; + +/** + * @brief Owns runtime-sized PCHIP data in unified memory. + */ +class PchipData { + public: + /** + * @brief Construct from strictly increasing coordinates and tabulated values. + */ + PchipData(std::span coordinates, std::span values) + : coordinates_(checkedSize(coordinates, values)), control_points_(checkedSize(coordinates, values) - 1) + { + std::vector widths(coordinates.size() - 1); + std::vector secant_slopes(coordinates.size() - 1); + std::vector nodal_slopes(coordinates.size()); + + for (std::size_t i = 0; i < coordinates.size(); ++i) { + if (!std::isfinite(coordinates[i]) || !std::isfinite(values[i])) { + throw std::invalid_argument("Pchip coordinates and values must be finite"); + } + if (i > 0 && coordinates[i] <= coordinates[i - 1]) { + throw std::invalid_argument("Pchip coordinates must be strictly increasing"); + } + coordinates_[static_cast(i)] = coordinates[i]; + } + + for (std::size_t i = 0; i + 1 < coordinates.size(); ++i) { + widths[i] = coordinates[i + 1] - coordinates[i]; + secant_slopes[i] = (values[i + 1] - values[i]) / widths[i]; + } + + if (coordinates.size() == 2) { + nodal_slopes[0] = secant_slopes[0]; + nodal_slopes[1] = secant_slopes[0]; + } else { + nodal_slopes[0] = endpointSlope(widths[0], widths[1], secant_slopes[0], secant_slopes[1]); + const std::size_t last = coordinates.size() - 1; + nodal_slopes[last] = + endpointSlope(widths[last - 1], widths[last - 2], secant_slopes[last - 1], secant_slopes[last - 2]); + + for (std::size_t i = 1; i < last; ++i) { + const double left_slope = secant_slopes[i - 1]; + const double right_slope = secant_slopes[i]; + if (left_slope == 0.0 || right_slope == 0.0 || std::signbit(left_slope) != std::signbit(right_slope)) { + nodal_slopes[i] = 0.0; + } else { + const double left_weight = 2.0 * widths[i] + widths[i - 1]; + const double right_weight = widths[i] + 2.0 * widths[i - 1]; + nodal_slopes[i] = (left_weight + right_weight) / (left_weight / left_slope + right_weight / right_slope); + } + } + } + + for (std::size_t i = 0; i + 1 < coordinates.size(); ++i) { + auto& control_points = control_points_[static_cast(i)]; + control_points[0] = values[i]; + control_points[1] = values[i] + widths[i] * nodal_slopes[i] / 3.0; + control_points[2] = values[i + 1] - widths[i] * nodal_slopes[i + 1] / 3.0; + control_points[3] = values[i + 1]; + } + } + + /** + * @brief Return lightweight callable view. + */ + PchipView view() const { return {coordinates_.data(), control_points_.data(), coordinates_.size()}; } + + private: + static axom::IndexType checkedSize(std::span coordinates, std::span values) + { + if (coordinates.size() != values.size()) { + throw std::invalid_argument("Pchip coordinates and values must have matching sizes"); + } + if (coordinates.size() < 2) { + throw std::invalid_argument("Pchip requires at least two points"); + } + if (coordinates.size() > static_cast(std::numeric_limits::max())) { + throw std::length_error("Pchip point count exceeds supported index range"); + } + return static_cast(coordinates.size()); + } + + /** + * @brief Compute shape-preserving derivative at one endpoint. + * + * PCHIP first forms a one-sided, three-point derivative estimate + * + * \f[ m = \frac{(2h_0 + h_1)d_0 - h_0d_1}{h_0 + h_1}, \f] + * + * where \f$h_i\f$ are interval widths and \f$d_i\f$ are secant slopes. + * If \f$m\f$ points opposite \f$d_0\f$, setting \f$m=0\f$ prevents the + * interpolant from initially moving away from the first interval data. When + * \f$d_0\f$ and \f$d_1\f$ have opposite signs, limiting \f$|m|\f$ to + * \f$3|d_0|\f$ prevents endpoint overshoot near the neighboring extremum. + * + * Explicit zero checks precede signbit comparisons. Comparing sign bits + * avoids multiplying slopes solely to determine whether signs differ. + * Reversing interval arguments applies the same rule at the right endpoint. + */ + static double endpointSlope(double first_width, double second_width, double first_secant, double second_secant) + { + double slope = ((2.0 * first_width + second_width) * first_secant - first_width * second_secant) / + (first_width + second_width); + if (slope == 0.0 || first_secant == 0.0 || std::signbit(slope) != std::signbit(first_secant)) { + return 0.0; + } + if (std::signbit(first_secant) != std::signbit(second_secant) && std::abs(slope) > 3.0 * std::abs(first_secant)) { + slope = 3.0 * first_secant; + } + return slope; + } + + axom::Array coordinates_; + axom::Array> control_points_; +}; + +} // namespace smith diff --git a/src/smith/numerics/tests/CMakeLists.txt b/src/smith/numerics/tests/CMakeLists.txt index 2fd4cedf6f..c359a0e305 100644 --- a/src/smith/numerics/tests/CMakeLists.txt +++ b/src/smith/numerics/tests/CMakeLists.txt @@ -10,6 +10,7 @@ set(numerics_serial_test_sources test_equationsolver.cpp test_operator.cpp test_odes.cpp + test_pchip.cpp test_steihaug_toint_cg.cpp test_block_preconditioner.cpp test_block_preconditioner_backend.cpp diff --git a/src/smith/numerics/tests/test_pchip.cpp b/src/smith/numerics/tests/test_pchip.cpp new file mode 100644 index 0000000000..609ca452a1 --- /dev/null +++ b/src/smith/numerics/tests/test_pchip.cpp @@ -0,0 +1,57 @@ +// Copyright (c) Lawrence Livermore National Security, LLC and +// other Smith Project Developers. See the top-level LICENSE file for +// details. +// +// SPDX-License-Identifier: (BSD-3-Clause) + +#include + +#include "gtest/gtest.h" + +#include "smith/infrastructure/application_manager.hpp" +#include "smith/numerics/pchip.hpp" + +namespace smith { + +TEST(Pchip, MonotoneAndClamped) +{ + const std::vector coordinates{0.0, 1.0, 2.0, 4.0}; + const std::vector values{10.0, 15.0, 18.0, 20.0}; + const PchipData data(coordinates, values); + const auto interpolant = data.view(); + + EXPECT_DOUBLE_EQ(interpolant(-1.0), values.front()); + EXPECT_DOUBLE_EQ(interpolant(5.0), values.back()); + + for (std::size_t i = 0; i + 1 < coordinates.size(); ++i) { + EXPECT_DOUBLE_EQ(interpolant(coordinates[i]), values[i]); + const double midpoint_value = interpolant(0.5 * (coordinates[i] + coordinates[i + 1])); + EXPECT_GE(midpoint_value, values[i]); + EXPECT_LE(midpoint_value, values[i + 1]); + } + EXPECT_DOUBLE_EQ(interpolant(coordinates.back()), values.back()); +} + +TEST(Pchip, SupportsRuntimeSizedData) +{ + constexpr std::size_t point_count = 64; + std::vector coordinates(point_count); + std::vector values(point_count); + for (std::size_t i = 0; i < point_count; ++i) { + coordinates[i] = static_cast(i); + values[i] = 2.0 * coordinates[i] + 1.0; + } + + const PchipData data(coordinates, values); + const auto interpolant = data.view(); + EXPECT_NEAR(interpolant(31.5), 64.0, 1.0e-14); +} + +} // namespace smith + +int main(int argc, char* argv[]) +{ + testing::InitGoogleTest(&argc, argv); + smith::ApplicationManager application_manager(argc, argv); + return RUN_ALL_TESTS(); +} diff --git a/src/smith/physics/tests/CMakeLists.txt b/src/smith/physics/tests/CMakeLists.txt index abfc9abac5..17fce5604d 100644 --- a/src/smith/physics/tests/CMakeLists.txt +++ b/src/smith/physics/tests/CMakeLists.txt @@ -34,6 +34,7 @@ set(physics_serial_test_sources test_functional_weak_form.cpp thermomech_statics_patch.cpp parameterized_thermomechanics_example.cpp + test_tabulated_solid_properties.cpp test_kinematic_objective.cpp ) diff --git a/src/smith/physics/tests/test_tabulated_solid_properties.cpp b/src/smith/physics/tests/test_tabulated_solid_properties.cpp new file mode 100644 index 0000000000..473b18385f --- /dev/null +++ b/src/smith/physics/tests/test_tabulated_solid_properties.cpp @@ -0,0 +1,123 @@ +// Copyright (c) Lawrence Livermore National Security, LLC and +// other Smith Project Developers. See the top-level LICENSE file for +// details. +// +// SPDX-License-Identifier: (BSD-3-Clause) + +#include +#include +#include + +#include "gtest/gtest.h" +#include "mfem.hpp" + +#include "smith/infrastructure/application_manager.hpp" +#include "smith/mesh_utils/mesh_utils.hpp" +#include "smith/numerics/functional/domain.hpp" +#include "smith/numerics/functional/finite_element.hpp" +#include "smith/numerics/functional/tensor.hpp" +#include "smith/numerics/pchip.hpp" +#include "smith/physics/materials/solid_material.hpp" +#include "smith/physics/mesh.hpp" +#include "smith/physics/solid_mechanics.hpp" +#include "smith/physics/state/state_manager.hpp" +#include "smith/smith_config.hpp" + +namespace smith { + +namespace { + +template +auto linearIsotropicStress(ModulusType youngs_modulus, double poisson_ratio, + const tensor& displacement_gradient) +{ + const auto shear_modulus = youngs_modulus / (2.0 * (1.0 + poisson_ratio)); + const auto bulk_modulus = youngs_modulus / (3.0 * (1.0 - 2.0 * poisson_ratio)); + const auto lambda = bulk_modulus - 2.0 * shear_modulus / 3.0; + const auto strain = 0.5 * (displacement_gradient + transpose(displacement_gradient)); + return lambda * tr(strain) * Identity() + 2.0 * shear_modulus * strain; +} + +struct TemperatureDependentLinearIsotropic { + using State = Empty; + + double density; + double poisson_ratio; + PchipView youngs_modulus; + + template + auto operator()(State&, const tensor& displacement_gradient, + TemperatureType temperature) const + { + const auto modulus = youngs_modulus(get<0>(temperature)); + return linearIsotropicStress(modulus, poisson_ratio, displacement_gradient); + } +}; + +} // namespace + +TEST(TabulatedSolidProperties, MaterialAndTractionCallbacks) +{ + constexpr int order = 1; + constexpr int dim = 2; + constexpr double poisson_ratio = 0.25; + const std::vector temperatures{0.0, 1.0, 2.0, 4.0}; + const std::vector moduli{10.0, 15.0, 18.0, 20.0}; + const PchipData youngs_modulus_data(temperatures, moduli); + const auto youngs_modulus = youngs_modulus_data.view(); + + axom::sidre::DataStore datastore; + StateManager::initialize(datastore, "test_tabulated_solid_properties"); + + auto mesh = std::make_shared(buildMeshFromFile(SMITH_REPO_DIR "/data/meshes/patch2D_quads.mesh"), "mesh"); + mesh->addDomainOfBoundaryElements("essential_boundary", by_attr(std::set{1, 4})); + + NonlinearSolverOptions nonlinear_options{.relative_tol = 1.0e-13, .absolute_tol = 1.0e-13}; + SolidMechanics>> solid(nonlinear_options, solid_mechanics::default_linear_options, + solid_mechanics::default_quasistatic_options, "solid", mesh, + {"temperature"}); + + FiniteElementState temperature(mesh->mfemParMesh(), H1{}, "temperature"); + temperature = 1.5; + solid.setParameter(0, temperature); + + TemperatureDependentLinearIsotropic material{ + .density = 1.0, .poisson_ratio = poisson_ratio, .youngs_modulus = youngs_modulus}; + solid.setMaterial(DependsOn<0>{}, material, mesh->entireBody()); + + constexpr tensor displacement_gradient{{{0.02, 0.01}, {-0.01, 0.03}}}; + constexpr tensor translation{{0.1, -0.2}}; + solid.setDisplacementBCs( + [=](tensor position, double) { return displacement_gradient * position + translation; }, + mesh->domain("essential_boundary")); + + solid.setTraction( + DependsOn<0>{}, + [=](auto, auto normal, double, auto temperature_value) { + const auto modulus = youngs_modulus(get<0>(temperature_value)); + const auto stress = linearIsotropicStress(modulus, poisson_ratio, displacement_gradient); + return stress * normal; + }, + mesh->entireBoundary()); + + solid.completeSetup(); + solid.advanceTimestep(1.0); + + auto exact_displacement = [=](const mfem::Vector& position, mfem::Vector& displacement) { + const tensor position_tensor{{position[0], position[1]}}; + const auto displacement_tensor = displacement_gradient * position_tensor + translation; + displacement[0] = displacement_tensor[0]; + displacement[1] = displacement_tensor[1]; + }; + mfem::VectorFunctionCoefficient exact_solution(dim, exact_displacement); + EXPECT_LT(computeL2Error(solid.displacement(), exact_solution), 1.0e-11); +} + +} // namespace smith + +int main(int argc, char* argv[]) +{ + testing::InitGoogleTest(&argc, argv); + smith::ApplicationManager application_manager(argc, argv); + return RUN_ALL_TESTS(); +} From 967a88443432b209c4c23fc8f6032c0d8c9b8679 Mon Sep 17 00:00:00 2001 From: Michael Tupek Date: Mon, 31 Aug 2026 09:49:43 -0700 Subject: [PATCH 2/3] Use array view to simplify a bit. --- src/smith/numerics/pchip.hpp | 17 ++++++++--------- 1 file changed, 8 insertions(+), 9 deletions(-) diff --git a/src/smith/numerics/pchip.hpp b/src/smith/numerics/pchip.hpp index 03b2307457..f5d5b616c4 100644 --- a/src/smith/numerics/pchip.hpp +++ b/src/smith/numerics/pchip.hpp @@ -38,12 +38,12 @@ class PchipView { if (input_value < coordinates_[0]) { return control_points_[0][0] + 0.0 * input; } - if (input_value > coordinates_[num_points_ - 1]) { - return control_points_[num_points_ - 2][3] + 0.0 * input; + if (input_value > coordinates_[coordinates_.size() - 1]) { + return control_points_[control_points_.size() - 1][3] + 0.0 * input; } axom::IndexType lower = 0; - axom::IndexType upper = num_points_ - 1; + axom::IndexType upper = coordinates_.size() - 1; while (upper - lower > 1) { const axom::IndexType middle = lower + (upper - lower) / 2; if (input_value < coordinates_[middle]) { @@ -66,14 +66,13 @@ class PchipView { private: friend class PchipData; - PchipView(const double* coordinates, const tensor* control_points, axom::IndexType num_points) - : coordinates_(coordinates), control_points_(control_points), num_points_(num_points) + PchipView(axom::ArrayView coordinates, axom::ArrayView> control_points) + : coordinates_(coordinates), control_points_(control_points) { } - const double* coordinates_; - const tensor* control_points_; - axom::IndexType num_points_; + axom::ArrayView coordinates_; + axom::ArrayView> control_points_; }; /** @@ -140,7 +139,7 @@ class PchipData { /** * @brief Return lightweight callable view. */ - PchipView view() const { return {coordinates_.data(), control_points_.data(), coordinates_.size()}; } + PchipView view() const { return {coordinates_.view(), control_points_.view()}; } private: static axom::IndexType checkedSize(std::span coordinates, std::span values) From 5ee3eb1dab624d9a1218df293187a1c424c071be Mon Sep 17 00:00:00 2001 From: Brandon Talamini <30813018+btalamini@users.noreply.github.com> Date: Mon, 14 Sep 2026 17:41:54 -0700 Subject: [PATCH 3/3] Add missing doxygen file description Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- src/smith/numerics/pchip.hpp | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/smith/numerics/pchip.hpp b/src/smith/numerics/pchip.hpp index f5d5b616c4..424fa5d418 100644 --- a/src/smith/numerics/pchip.hpp +++ b/src/smith/numerics/pchip.hpp @@ -4,6 +4,12 @@ // // SPDX-License-Identifier: (BSD-3-Clause) +/** + * @file pchip.hpp + * + * @brief Shape-preserving piecewise cubic Hermite (PCHIP) interpolation for 1D tabulated data. + */ + #pragma once #include