Merge pull request #1310 from LLNL/talamini/feature/initial_condition…

…s_with_domain

Change displacement and velocity setters to use serac tensors

src/serac/physics/solid_mechanics.hpp

@@ -820,30 +820,42 @@ class SolidMechanics<order, dim, Parameters<parameter_space...>, std::integer_se
   }
 
   /**
-   * @brief Set the underlying finite element state to a prescribed displacement
+   * @brief Set the underlying finite element state to a prescribed displacement field
    *
-   * @param disp The function describing the displacement field
+   * @param applied_displacement The displacement field as a callable,
+   *   which must have the signature:
+   *   tensor<double, dim> applied_displacement(tensor<double, dim> X)
+   *
+   *   args:
+   *   X: coordinates of the material point in the reference configuration
+   *
+   *   returns: u, the displacement at X.
    */
-  void setDisplacement(std::function<void(const mfem::Vector& x, mfem::Vector& disp)> disp)
+  template <typename AppliedDisplacementFunction>
+  void setDisplacement(AppliedDisplacementFunction applied_displacement)
   {
-    // Project the coefficient onto the grid function
-    mfem::VectorFunctionCoefficient disp_coef(dim, disp);
-    displacement_.project(disp_coef);
+    displacement_.setFromFieldFunction(applied_displacement);
   }
 
   /// @overload
   void setDisplacement(const FiniteElementState& temp) { displacement_ = temp; }
 
   /**
-   * @brief Set the underlying finite element state to a prescribed velocity
+   * @brief Set the underlying finite element state to a prescribed velocity field
+   *
+   * @param applied_velocity The velocity field as a callable,
+   *   which must have the signature:
+   *   tensor<double, dim> applied_velocity(tensor<double, dim> X)
+   *
+   *   args:
+   *   X: coordinates of the material point in the reference configuration
    *
-   * @param vel The function describing the velocity field
+   *   returns: v, the velocity at X.
    */
-  void setVelocity(std::function<void(const mfem::Vector& x, mfem::Vector& vel)> vel)
+  template <typename AppliedVelocityFunction>
+  void setVelocity(AppliedVelocityFunction applied_velocity)
   {
-    // Project the coefficient onto the grid function
-    mfem::VectorFunctionCoefficient vel_coef(dim, vel);
-    velocity_.project(vel_coef);
+    velocity_.setFromFieldFunction(applied_velocity);
   }
 
   /// @overload

src/serac/physics/state/CMakeLists.txt

@@ -32,3 +32,7 @@ install(TARGETS              serac_state
         EXPORT               serac-targets
         DESTINATION          lib
         )
+
+if(SERAC_ENABLE_TESTS)
+    add_subdirectory(tests)
+endif()

src/serac/physics/state/finite_element_state.hpp

@@ -25,6 +25,71 @@
 
 namespace serac {
 
+namespace detail {
+/**
+ * @brief Helper function to copy a tensor into an mfem Vector
+ */
+template <int dim>
+void setMfemVectorFromTensorOrDouble(mfem::Vector& v_mfem, const tensor<double, dim>& v)
+{
+  SLIC_ERROR_IF(v_mfem.Size() != dim, "Cannot copy tensor into an MFEM Vector with incompatible size.");
+  for (int i = 0; i < dim; i++) v_mfem(i) = v[i];
+}
+
+/**
+ * @brief Helper function to copy a double into a singleton mfem Vector
+ */
+inline void setMfemVectorFromTensorOrDouble(mfem::Vector& v_mfem, double v)
+{
+  SLIC_ERROR_IF(v_mfem.Size() != 1, "Mfem Vector is not a singleton.");
+  v_mfem = v;
+}
+
+/**
+ * @brief Helper for extracting type of first argument of a free function
+ */
+template <typename Ret, typename Arg, typename... Rest>
+Arg first_argument_helper(Ret (*)(Arg, Rest...));
+
+/**
+ * @brief Helper for extracting type of first argument of a class method
+ */
+template <typename Ret, typename F, typename Arg, typename... Rest>
+Arg first_argument_helper(Ret (F::*)(Arg, Rest...));
+
+/**
+ * @brief Helper for extracting type of first argument of a const class method
+ */
+template <typename Ret, typename F, typename Arg, typename... Rest>
+Arg first_argument_helper(Ret (F::*)(Arg, Rest...) const);
+
+/**
+ * Extract type of first argument of a class method
+ */
+template <typename F>
+decltype(first_argument_helper(&F::operator())) first_argument_helper(F);
+
+/**
+ * Extract type of first argument of a free callable
+ */
+template <typename T>
+using first_argument = decltype(first_argument_helper(std::declval<T>()));
+
+/**
+ * @brief Evaluate a function of a tensor with an mfem Vector object
+ */
+template <typename Callable>
+auto evaluateTensorFunctionOnMfemVector(const mfem::Vector& X_mfem, Callable&& f)
+{
+  first_argument<Callable> X;
+  SLIC_ERROR_IF(X_mfem.Size() != size(X),
+                "Size of tensor in callable does not match spatial dimension of MFEM Vector.");
+  for (int i = 0; i < X_mfem.Size(); i++) X[i] = X_mfem[i];
+  return f(X);
+}
+
+}  // namespace detail
+
 /**
  * @brief convenience function for querying the type stored in a GeneralCoefficient
  */
@@ -206,6 +271,40 @@ class FiniteElementState : public FiniteElementVector {
   /// \overload
   void project(mfem::VectorCoefficient& coef, const Domain& d);
 
+  /**
+   * @brief Set state as interpolant of an analytical function
+   *
+   * In other words, this sets the dofs by direct evaluation of the analytical
+   * field at the nodal points.
+   *
+   * @tparam FieldFunction A callable type
+   * @param field_function A function that should have the signature:
+   *   tensor<double, field_dim> field_function(tensor<double, spatial_dim> X)
+   *
+   *   template params:
+   *     spatial_dim: number of components in the mesh coordinates
+   *     field_dim: number of components of the FiniteElementState field
+   *
+   *   args:
+   *     X: coordinates of the material point
+   *
+   *   returns: the value of the field at X.
+   *
+   *   If the field is scalar-valued, then alternatively the signature may be:
+   *   double field_function(tensor<double, spatial_dim> X)
+   */
+  template <typename FieldFunction>
+  void setFromFieldFunction(FieldFunction&& field_function)
+  {
+    auto evaluate_mfem = [&field_function](const mfem::Vector& X_mfem, mfem::Vector& u_mfem) {
+      auto u = detail::evaluateTensorFunctionOnMfemVector(X_mfem, field_function);
+      detail::setMfemVectorFromTensorOrDouble(u_mfem, u);
+    };
+
+    mfem::VectorFunctionCoefficient coef(space_->GetVDim(), evaluate_mfem);
+    project(coef);
+  }
+
   /**
    * @brief Construct a grid function from the finite element state true vector
    *

src/serac/physics/state/tests/CMakeLists.txt

@@ -0,0 +1,13 @@
+# 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)
+
+set(state_test_depends serac_mesh serac_state gtest)
+
+set(state_test_sources
+    test_finite_element_state.cpp)
+
+serac_add_tests(SOURCES    ${state_test_sources}
+                DEPENDS_ON ${state_test_depends})

src/serac/physics/state/tests/test_finite_element_state.cpp

@@ -0,0 +1,138 @@
+// 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)
+
+/**
+ * @file test_finite_element_staet.cpp
+ */
+
+#include "serac/physics/state/finite_element_state.hpp"
+
+#include "axom/slic/core/SimpleLogger.hpp"
+#include <gtest/gtest.h>
+#include "mfem.hpp"
+
+#include "serac/infrastructure/initialize.hpp"
+#include "serac/infrastructure/terminator.hpp"
+#include "serac/mesh/mesh_utils.hpp"
+#include "serac/numerics/functional/tensor.hpp"
+
+namespace serac {
+
+class TestFiniteElementState : public testing::Test {
+ protected:
+  void SetUp() override
+  {
+    int serial_refinement = 0;
+    int parallel_refinement = 0;
+
+    std::string filename = SERAC_REPO_DIR "/data/meshes/beam-hex.mesh";
+    mesh = mesh::refineAndDistribute(buildMeshFromFile(filename), serial_refinement, parallel_refinement);
+    ASSERT_EQ(spatial_dim, mesh->SpaceDimension())
+        << "Test configured incorrectly. The variable spatial_dim must match the spatial dimension of the mesh.";
+  }
+
+  static constexpr int spatial_dim{3};
+  std::unique_ptr<mfem::ParMesh> mesh;
+};
+
+TEST_F(TestFiniteElementState, SetScalarStateFromFieldFunction)
+{
+  // make a scalar-valued state
+  constexpr int p = 1;
+  FiniteElementState scalar_state(*mesh, H1<p>{}, "scalar_field");
+
+  // Set state with field function.
+  // Check that lambda captures work with this.
+  double c = 2.0;
+  auto scalar_field = [c](tensor<double, spatial_dim> X) -> double { return c * X[0]; };
+  scalar_state.setFromFieldFunction(scalar_field);
+
+  // Get the nodal positions corresponding to state dofs in a grid function
+  auto [coords_fe_space, coords_fe_coll] = serac::generateParFiniteElementSpace<H1<p, spatial_dim>>(mesh.get());
+  mfem::ParGridFunction nodal_coords_gf(coords_fe_space.get());
+  mesh->GetNodes(nodal_coords_gf);
+
+  for (int node = 0; node < scalar_state.space().GetNDofs(); node++) {
+    // Fill a tensor with the coordinates of the node
+    tensor<double, spatial_dim> Xn;
+    for (int i = 0; i < spatial_dim; i++) {
+      int dof_index = nodal_coords_gf.FESpace()->DofToVDof(node, i);
+      Xn[i] = nodal_coords_gf(dof_index);
+    }
+
+    // check that value set in the state for this node matches the field function
+    EXPECT_DOUBLE_EQ(scalar_field(Xn), scalar_state(node));
+  }
+}
+
+TEST_F(TestFiniteElementState, SetVectorStateFromFieldFunction)
+{
+  constexpr int p = 2;
+
+  // Choose vector dimension for state field that is different from spatial dimension
+  // to test the field indexing more thoroughly.
+  constexpr int vdim = 2;
+  FiniteElementState state(*mesh, H1<p, vdim>{}, "vector_field");
+
+  // set the field with an arbitrarily chosen field function
+  auto vector_field = [](tensor<double, spatial_dim> X) {
+    return tensor<double, vdim>{norm(X), 1.0 / (1.0 + norm(X))};
+  };
+  state.setFromFieldFunction(vector_field);
+
+  // Get the nodal positions for the state in a grid function
+  auto [coords_fe_space, coords_fe_coll] = serac::generateParFiniteElementSpace<H1<p, spatial_dim>>(mesh.get());
+  mfem::ParGridFunction nodal_coords_gf(coords_fe_space.get());
+  mesh->GetNodes(nodal_coords_gf);
+
+  // we need the state values and the nodal coordinates in the same kind of container,
+  // so we will get the grid function view of the state
+  mfem::ParGridFunction& state_gf = state.gridFunction();
+
+  for (int node = 0; node < state_gf.FESpace()->GetNDofs(); node++) {
+    // Fill a tensor with the coordinates of the node
+    tensor<double, spatial_dim> Xn;
+    for (int i = 0; i < spatial_dim; i++) {
+      int dof_index = nodal_coords_gf.FESpace()->DofToVDof(node, i);
+      Xn[i] = nodal_coords_gf(dof_index);
+    }
+
+    // apply the field function to the node coords
+    auto v = vector_field(Xn);
+
+    // check that value set in the state matches the field function
+    for (int j = 0; j < vdim; j++) {
+      int dof_index = state_gf.FESpace()->DofToVDof(node, j);
+      EXPECT_DOUBLE_EQ(v[j], state_gf(dof_index));
+    }
+  }
+}
+
+TEST_F(TestFiniteElementState, DISABLED_ErrorsIfFieldFunctionDimensionMismatchedToState)
+{
+  constexpr int p = 2;
+
+  // Choose vector dimension for state field that is different from spatial dimension
+  constexpr int vdim = 2;
+  FiniteElementState state(*mesh, H1<p, vdim>{}, "vector_field");
+
+  // Set the field with a field function with the wrong vector dimension.
+  // Should return tensor of size vdim!
+  auto vector_field = [](tensor<double, spatial_dim> X) { return X; };
+
+  EXPECT_DEATH(state.setFromFieldFunction(vector_field),
+               "Cannot copy tensor into an MFEM Vector with incompatible size.");
+}
+
+}  // namespace serac
+
+int main(int argc, char* argv[])
+{
+  ::testing::InitGoogleTest(&argc, argv);
+  serac::initialize(argc, argv);
+  int result = RUN_ALL_TESTS();
+  serac::exitGracefully(result);
+}

src/serac/physics/tests/beam_bending.cpp

@@ -76,9 +76,6 @@ TEST(BeamBending, TwoDimensional)
   Domain support = Domain::ofBoundaryElements(pmesh, by_attr<dim>(1));
   solid_solver.setFixedBCs(support);
 
-  // initial displacement
-  solid_solver.setDisplacement([](const mfem::Vector&, mfem::Vector& bc_vec) { bc_vec = 0.0; });
-
   Domain top_face = Domain::ofBoundaryElements(
       pmesh, [](std::vector<vec2> vertices, int /*attr*/) { return (average(vertices)[1] > 0.99); });
 

src/serac/physics/tests/fit_test.cpp

@@ -83,9 +83,6 @@ void stress_extrapolation_test()
   auto down = [up](tensor<double, dim> X, double time) { return -up(X, time); };
   solid_solver.setDisplacementBCs(down, bottom_hole);
 
-  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();
 

src/serac/physics/tests/lce_Bertoldi_lattice.cpp

@@ -124,9 +124,9 @@ TEST(LiquidCrystalElastomer, Bertoldi)
   auto support = Domain::ofBoundaryElements(pmesh, by_attr<dim>(2));
   solid_solver.setFixedBCs(support);
 
-  double iniDispVal = 5.0e-6;
-  auto ini_displacement = [iniDispVal](const mfem::Vector&, mfem::Vector& u) -> void { u = iniDispVal; };
-  solid_solver.setDisplacement(ini_displacement);
+  constexpr double iniDispVal = 5.0e-6;
+  auto initial_displacement = [lx](vec3 X) { return vec3{0.0, (X[0] / lx) * iniDispVal, 0.0}; };
+  solid_solver.setDisplacement(initial_displacement);
 
   // Finalize the data structures
   solid_solver.completeSetup();

src/serac/physics/tests/lce_Brighenti_tensile.cpp

@@ -135,16 +135,11 @@ TEST(LiquidCrystalElastomer, Brighenti)
   solid_solver.setFixedBCs(ymin_face, Component::Y);
   solid_solver.setFixedBCs(zmin_face, Component::Z);
 
-  // set initila displacement different than zero to help solver
-  auto ini_displacement = [](const mfem::Vector&, mfem::Vector& u) -> void { u = 1.0e-5; };
-
   double iniLoadVal = 1.0e0;
   double maxLoadVal = 4 * 1.3e0 / lx / lz;
   double loadVal = iniLoadVal + 0.0 * maxLoadVal;
   solid_solver.setTraction([&loadVal](auto, auto n, auto) { return loadVal * n; }, ymax_face);
 
-  solid_solver.setDisplacement(ini_displacement);
-
   // Finalize the data structures
   solid_solver.completeSetup();
 

src/serac/physics/tests/parameterized_thermomechanics_example.cpp

@@ -134,10 +134,6 @@ TEST(Thermomechanics, ParameterizedMaterial)
   simulation.setFixedBCs(y_equals_0, Component::Y);
   simulation.setFixedBCs(z_equals_0, Component::Z);
 
-  // set up initial conditions
-  auto zero_vector = [](const mfem::Vector&, mfem::Vector& u) -> void { u = 0.0; };
-  simulation.setDisplacement(zero_vector);
-
   // Finalize the data structures
   simulation.completeSetup();
 

src/serac/physics/tests/solid.cpp

@@ -196,7 +196,6 @@ void functional_parameterized_solid_test(double expected_disp_norm)
   Domain essential_boundary = Domain::ofBoundaryElements(pmesh, by_attr<dim>(1));
 
   solid_solver.setFixedBCs(essential_boundary);
-  solid_solver.setDisplacement([](const mfem::Vector&, mfem::Vector& bc_vec) -> void { bc_vec = 0.0; });
 
   tensor<double, dim> constant_force;