stress extrapolation

src/serac/infrastructure/variant.hpp

@@ -229,20 +229,20 @@ struct variant {
   /**
    * @brief "Parameterized" constructor with which a value can be assigned
    * @tparam T The type of the parameter to assign to one of the variant elements
-   * @param[in] t The parameter to assign to the variant's contents
+   * @param[in] t_ The parameter to assign to the variant's contents
    * @pre The parameter type @p T must be equal to or assignable to either @p T0 or @p T1
    * @note If the conversion is ambiguous, i.e., if @a t is equal or convertible to *both* @p T0 and @p T1,
    * the first element of the variant - of type @p T0 - will be assigned to
    */
   template <typename T, typename SFINAE = std::enable_if_t<detail::is_variant_assignable<T, T0, T1>::value>>
-  constexpr variant(T&& t)
+  constexpr variant(T&& t_)
   {
     if constexpr (std::is_same_v<std::decay_t<T>, T0> || std::is_assignable_v<T0, T>) {
       storage_.index_ = 0;
-      new (&storage_.t0_) T0(std::forward<T>(t));
+      new (&storage_.t0_) T0(std::forward<T>(t_));
     } else if constexpr (std::is_same_v<std::decay_t<T>, T1> || std::is_assignable_v<T1, T>) {
       storage_.index_ = 1;
-      new (&storage_.t1_) T1(std::forward<T>(t));
+      new (&storage_.t1_) T1(std::forward<T>(t_));
     } else {
       static_assert(sizeof(T) < 0, "Type not supported");
     }

src/serac/numerics/functional/tensor.hpp

@@ -1064,6 +1064,22 @@ SERAC_HOST_DEVICE auto normalize(const tensor<T, n...>& A)
   return A / norm(A);
 }
 
+/**
+ * @brief promotes a 2x2 matrix to a 3x3 matrix, by populating the upper left block, leaving zeroes in the third row /
+ * column
+ * @param[in] A the 2x2 matrix to be promoted to a 3x3 matrix
+ */
+template <typename T>
+SERAC_HOST_DEVICE tensor<T, 3, 3> to_3x3(const tensor<T, 2, 2>& A)
+{
+  tensor<T, 3, 3> output{};
+  output[0][0] = A[0][0];
+  output[0][1] = A[0][1];
+  output[1][0] = A[1][0];
+  output[1][1] = A[1][1];
+  return output;
+}
+
 /**
  * @brief Returns the trace of a square matrix
  * @param[in] A The matrix to compute the trace of
@@ -1207,6 +1223,17 @@ SERAC_HOST_DEVICE constexpr auto transpose(const tensor<T, m, n>& A)
   return AT;
 }
 
+/**
+ * @brief Returns the second invariant of a 3x3 matrix
+ * @param[in] A The matrix
+ */
+template <typename T>
+SERAC_HOST_DEVICE constexpr auto I2(const tensor<T, 3, 3>& A)
+{
+  return +A[0][0] * A[1][1] + A[1][1] * A[2][2] + A[2][2] * A[0][0] - A[0][1] * A[1][0] - A[1][2] * A[2][1] -
+         A[2][0] * A[0][2];
+}
+
 /**
  * @brief Returns the determinant of a matrix
  * @param[in] A The matrix to obtain the determinant of
@@ -1216,6 +1243,7 @@ SERAC_HOST_DEVICE constexpr auto det(const tensor<T, 2, 2>& A)
 {
   return A[0][0] * A[1][1] - A[0][1] * A[1][0];
 }
+
 /// @overload
 template <typename T>
 SERAC_HOST_DEVICE constexpr auto det(const tensor<T, 3, 3>& A)

src/serac/physics/CMakeLists.txt

@@ -20,6 +20,7 @@ set(physics_sources
 set(physics_headers
     base_physics.hpp
     common.hpp
+    fit.hpp
     heat_transfer.hpp
     heat_transfer_input.hpp
     solid_mechanics.hpp

src/serac/physics/fit.hpp

@@ -0,0 +1,72 @@
+// Copyright (c) 2019-2024, Lawrence Livermore National Security, LLC and
+// other Serac Project Developers. See the top-level LICENSE file for
+// details.
+//
+// SPDX-License-Identifier: (BSD-3-Clause)
+
+#include <fstream>
+#include <iostream>
+
+#include "serac/numerics/functional/functional.hpp"
+#include "serac/numerics/functional/tensor.hpp"
+#include "serac/physics/state/finite_element_state.hpp"
+
+namespace serac {
+
+namespace detail {
+
+/// @overload
+template <int dim, typename signature, int... i, typename func, typename... T>
+FiniteElementState fit(std::integer_sequence<int, i...>, func f, mfem::ParMesh& mesh, const T&... solution_fields)
+{
+  // signature looks like return_type(arg0_type, arg1_type);
+  // so this unpacks the return type
+  using output_space = typename FunctionSignature<signature>::return_type;
+
+  FiniteElementState fitted_field(mesh, output_space{});
+  fitted_field = 0.0;
+
+  // mass term
+  serac::Functional<output_space(output_space)> phi_phi(&fitted_field.space(), {&fitted_field.space()});
+  phi_phi.AddDomainIntegral(
+      Dimension<dim>{}, DependsOn<0>{},
+      [](double /*t*/, auto /*x*/, auto u) {
+        return tuple{get<0>(u), zero{}};
+      },
+      mesh);
+  auto M = get<1>(phi_phi(DifferentiateWRT<0>{}, 0.0 /* t */, fitted_field));
+
+  // rhs
+  std::array<const mfem::ParFiniteElementSpace*, sizeof...(T)> trial_spaces = {&solution_fields.space()...};
+  serac::Functional<signature>                                 phi_f(&fitted_field.space(), trial_spaces);
+  phi_f.AddDomainIntegral(Dimension<dim>{}, DependsOn<i...>{}, f, mesh);
+  mfem::Vector b = phi_f(0.0, solution_fields...);
+
+  mfem::CGSolver cg(MPI_COMM_WORLD);
+  cg.SetOperator(M);
+  cg.SetRelTol(1e-12);
+  cg.SetMaxIter(500);
+  cg.SetPrintLevel(2);
+  cg.Mult(b, fitted_field);
+
+  return fitted_field;
+}
+
+}  // namespace detail
+
+/**
+ * @brief determine field parameters to approximate the output of a user-provided q-function
+ * @param[in] f the user-provided function to approximate
+ * @param[in] mesh the region over which to approximate the function f
+ * @param[in] solution_fields [optional] any auxiliary field quantities needed to evaluate f
+ *
+ * @note: mesh is passed by non-const ref because mfem mutates the mesh when creating ParGridFunctions
+ */
+template <int dim, typename signature, int... n, typename func, typename... T>
+FiniteElementState fit(func f, mfem::ParMesh& mesh, const T&... solution_fields)
+{
+  auto iseq = std::make_integer_sequence<int, sizeof...(T)>{};
+  return detail::fit<dim, signature>(iseq, f, mesh, solution_fields...);
+}
+
+}  // namespace serac

src/serac/physics/solid_mechanics.hpp

@@ -1625,8 +1625,8 @@ class SolidMechanics<order, dim, Parameters<parameter_space...>, std::integer_se
   std::array<std::function<decltype((*residual_)(DifferentiateWRT<1>{}, 0.0, shape_displacement_, displacement_,
                                                  acceleration_, *parameters_[parameter_indices].state...))(double)>,
              sizeof...(parameter_indices)>
-      d_residual_d_ = {[&](double t) {
-        return (*residual_)(DifferentiateWRT<NUM_STATE_VARS + 1 + parameter_indices>{}, t, shape_displacement_,
+      d_residual_d_ = {[&](double _t) {
+        return (*residual_)(DifferentiateWRT<NUM_STATE_VARS + 1 + parameter_indices>{}, _t, shape_displacement_,
                             displacement_, acceleration_, *parameters_[parameter_indices].state...);
       }...};
 

src/serac/physics/tests/CMakeLists.txt

@@ -8,6 +8,7 @@ set(physics_test_depends serac_physics serac_mesh gtest)
 
 set(physics_serial_test_sources
     beam_bending.cpp
+    fit_test.cpp
     thermal_finite_diff.cpp
     thermal_statics_patch.cpp
     thermal_dynamics_patch.cpp

src/serac/physics/tests/fit_test.cpp

@@ -0,0 +1,125 @@
+// Copyright (c) 2019-2024, Lawrence Livermore National Security, LLC and
+// other Serac Project Developers. See the top-level LICENSE file for
+// details.
+//
+// SPDX-License-Identifier: (BSD-3-Clause)
+
+#include "serac/physics/fit.hpp"
+#include "serac/numerics/functional/functional.hpp"
+
+#include "serac/mesh/mesh_utils.hpp"
+#include "serac/physics/solid_mechanics.hpp"
+#include "serac/physics/materials/solid_material.hpp"
+
+#include <gtest/gtest.h>
+
+using namespace serac;
+using namespace serac::profiling;
+
+int num_procs, myid;
+int nsamples = 1;  // because mfem doesn't take in unsigned int
+
+int    n = 0;
+double t = 0.0;
+
+template <typename output_space>
+void stress_extrapolation_test()
+{
+  int serial_refinement   = 2;
+  int parallel_refinement = 0;
+
+  std::string filename = SERAC_REPO_DIR "/data/meshes/notched_plate.mesh";
+
+  auto mesh_ = mesh::refineAndDistribute(buildMeshFromFile(filename), serial_refinement, parallel_refinement);
+
+  constexpr int p   = 2;
+  constexpr int dim = 2;
+
+  using input_space = H1<2, dim>;
+
+  // Create DataStore
+  axom::sidre::DataStore datastore;
+  serac::StateManager::initialize(datastore, "solid_mechanics_J2_test");
+
+  // Construct the appropriate dimension mesh and give it to the data store
+
+  std::string mesh_tag{"mesh"};
+
+  auto& pmesh = StateManager::setMesh(std::move(mesh_), mesh_tag);
+
+  LinearSolverOptions linear_options{.linear_solver = LinearSolver::SuperLU};
+
+  NonlinearSolverOptions nonlinear_options{.nonlin_solver  = NonlinearSolver::Newton,
+                                           .relative_tol   = 1.0e-12,
+                                           .absolute_tol   = 1.0e-12,
+                                           .max_iterations = 5000,
+                                           .print_level    = 1};
+
+  FiniteElementState sigma_J2(pmesh, output_space{}, "sigma_J2");
+
+  SolidMechanics<p, dim, serac::Parameters<output_space> > solid_solver(
+      nonlinear_options, linear_options, solid_mechanics::default_quasistatic_options, GeometricNonlinearities::Off,
+      "solid_mechanics", mesh_tag, {"sigma_J2"});
+
+  solid_mechanics::NeoHookean mat{
+      1.0,    // density
+      100.0,  // bulk modulus
+      50.0    // shear modulus
+  };
+
+  solid_solver.setMaterial(mat);
+
+  // prescribe small displacement at each hole, pulling the plate apart
+  std::set<int> top_hole = {2};
+  auto          up       = [](const mfem::Vector&, mfem::Vector& u) -> void { u[1] = 0.01; };
+  solid_solver.setDisplacementBCs(top_hole, up);
+
+  std::set<int> bottom_hole = {3};
+  auto          down        = [](const mfem::Vector&, mfem::Vector& u) -> void { u[1] = -0.01; };
+  solid_solver.setDisplacementBCs(bottom_hole, down);
+
+  auto zero_displacement = [](const mfem::Vector&, mfem::Vector& u) -> void { u = 0.0; };
+  solid_solver.setDisplacement(zero_displacement);
+
+  // Finalize the data structures
+  solid_solver.completeSetup();
+
+  solid_solver.outputStateToDisk("paraview" + std::to_string(n));
+
+  double dt = 1.0;
+  solid_solver.advanceTimestep(dt);
+
+  auto u = solid_solver.displacement();
+
+  Empty internal_variables{};
+
+  sigma_J2 = fit<dim, output_space(input_space)>(
+      [&](double /*t*/, [[maybe_unused]] auto position, [[maybe_unused]] auto displacement_) {
+        mat3 du_dx  = to_3x3(get_value(get<1>(displacement_)));
+        auto stress = mat(internal_variables, du_dx);
+        return tuple{I2(dev(stress)), zero{}};
+      },
+      pmesh, u);
+
+  solid_solver.setParameter(0, sigma_J2);
+
+  solid_solver.outputStateToDisk("paraview" + std::to_string(n));
+  n++;
+}
+
+TEST(StressExtrapolation, PiecewiseConstant2D) { stress_extrapolation_test<L2<0> >(); }
+TEST(StressExtrapolation, PiecewiseLinear2D) { stress_extrapolation_test<H1<1> >(); }
+
+int main(int argc, char* argv[])
+{
+  ::testing::InitGoogleTest(&argc, argv);
+  MPI_Init(&argc, &argv);
+  MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
+  MPI_Comm_rank(MPI_COMM_WORLD, &myid);
+
+  axom::slic::SimpleLogger logger;
+
+  int result = RUN_ALL_TESTS();
+  MPI_Finalize();
+  return result;
+}