Merge pull request #1249 from LLNL/talamini/feature/multi_material

Multi-material solid mechanics

src/serac/numerics/functional/functional.hpp

@@ -301,7 +301,7 @@ class Functional<test(trials...), exec> {
 
   /// @overload
   template <int dim, int... args, typename lambda, typename qpt_data_type = Nothing>
-  void AddDomainIntegral(Dimension<dim>, DependsOn<args...>, const lambda& integrand, Domain& domain,
+  void AddDomainIntegral(Dimension<dim>, DependsOn<args...>, const lambda& integrand, const Domain& domain,
                          std::shared_ptr<QuadratureData<qpt_data_type>> qdata = NoQData)
   {
     if (domain.mesh_.GetNE() == 0) return;

src/serac/numerics/functional/geometry.hpp

@@ -2,6 +2,8 @@
 
 #include "mfem.hpp"
 
+#include "serac/numerics/functional/domain.hpp"
+
 namespace serac {
 
 /**
@@ -76,6 +78,23 @@ inline std::array<uint32_t, mfem::Geometry::NUM_GEOMETRIES> geometry_counts(cons
   return counts;
 }
 
+/**
+ * @brief count the number of elements of each geometry in a domain
+ * @param domain the domain to count
+ */
+inline std::array<uint32_t, mfem::Geometry::NUM_GEOMETRIES> geometry_counts(const Domain& domain)
+{
+  std::array<uint32_t, mfem::Geometry::NUM_GEOMETRIES> counts{};
+
+  constexpr std::array<mfem::Geometry::Type, 5> geometries = {mfem::Geometry::SEGMENT, mfem::Geometry::TRIANGLE,
+                                                              mfem::Geometry::SQUARE, mfem::Geometry::TETRAHEDRON,
+                                                              mfem::Geometry::CUBE};
+  for (auto geom : geometries) {
+    counts[uint32_t(geom)] = uint32_t(domain.get(geom).size());
+  }
+  return counts;
+}
+
 /**
  * @brief count the number of boundary elements of each geometry in a mesh
  * @param mesh the mesh to count

src/serac/physics/solid_mechanics.hpp

@@ -416,12 +416,14 @@ class SolidMechanics<order, dim, Parameters<parameter_space...>, std::integer_se
    *
    * @tparam T the type to be created at each quadrature point
    * @param initial_state the value to be broadcast to each quadrature point
+   * @param optional_domain restricts the quadrature buffer to the given domain
    * @return std::shared_ptr< QuadratureData<T> >
    */
   template <typename T>
-  qdata_type<T> createQuadratureDataBuffer(T initial_state)
+  qdata_type<T> createQuadratureDataBuffer(T initial_state, const std::optional<Domain>& optional_domain = std::nullopt)
   {
-    return StateManager::newQuadratureDataBuffer(mesh_tag_, order, dim, initial_state);
+    Domain domain = (optional_domain) ? *optional_domain : EntireDomain(mesh_);
+    return StateManager::newQuadratureDataBuffer(domain, order, dim, initial_state);
   }
 
   /**
@@ -895,6 +897,7 @@ class SolidMechanics<order, dim, Parameters<parameter_space...>, std::integer_se
    *    values will change to `dual` numbers rather than `double`. (e.g. `tensor<double,3>` becomes `tensor<dual<...>,
    * 3>`)
    *
+   * @param domain The subdomain which will use this material. Must be a domain of elements.
    * @param qdata the buffer of material internal variables at each quadrature point
    *
    * @pre MaterialType must have a public member variable `density`
@@ -903,7 +906,7 @@ class SolidMechanics<order, dim, Parameters<parameter_space...>, std::integer_se
    * @note This method must be called prior to completeSetup()
    */
   template <int... active_parameters, typename MaterialType, typename StateType = Empty>
-  void setMaterial(DependsOn<active_parameters...>, const MaterialType& material,
+  void setMaterial(DependsOn<active_parameters...>, const MaterialType& material, const Domain& domain,
                    qdata_type<StateType> qdata = EmptyQData)
   {
     static_assert(std::is_same_v<StateType, Empty> || std::is_same_v<StateType, typename MaterialType::State>,
@@ -916,14 +919,30 @@ class SolidMechanics<order, dim, Parameters<parameter_space...>, std::integer_se
                                                              // fact that the displacement, acceleration, and shape
                                                              // fields are always-on and come first, so the `n`th
                                                              // parameter will actually be argument `n + NUM_STATE_VARS`
-        std::move(material_functor), mesh_, qdata);
+        std::move(material_functor), domain, qdata);
+  }
+
+  /// @overload
+  template <int... active_parameters, typename MaterialType, typename StateType = Empty>
+  void setMaterial(DependsOn<active_parameters...>, const MaterialType& material,
+                   qdata_type<StateType> qdata = EmptyQData)
+  {
+    setMaterial(DependsOn<active_parameters...>{}, material, EntireDomain(mesh_), qdata);
+  }
+
+  /// @overload
+  template <typename MaterialType, typename StateType = Empty>
+  void setMaterial(const MaterialType& material, const Domain& domain,
+                   std::shared_ptr<QuadratureData<StateType>> qdata = EmptyQData)
+  {
+    setMaterial(DependsOn<>{}, material, domain, qdata);
   }
 
   /// @overload
   template <typename MaterialType, typename StateType = Empty>
   void setMaterial(const MaterialType& material, std::shared_ptr<QuadratureData<StateType>> qdata = EmptyQData)
   {
-    setMaterial(DependsOn<>{}, material, qdata);
+    setMaterial(DependsOn<>{}, material, EntireDomain(mesh_), qdata);
   }
 
   /**

src/serac/physics/state/state_manager.hpp

@@ -99,21 +99,19 @@ class StateManager {
    * @brief Create a shared ptr to a quadrature data buffer for the given material type
    *
    * @tparam T the type to be created at each quadrature point
-   * @param mesh_tag The tag for the stored mesh used to construct the finite element state
-   * @param order The order of the discretization of the displacement and velocity fields
+   * @param domain The spatial domain over which to allocate the quadrature data
+   * @param order The order of the discretization of the primal fields
    * @param dim The spatial dimension of the mesh
    * @param initial_state the value to be broadcast to each quadrature point
    * @return shared pointer to quadrature data buffer
    */
   template <typename T>
-  static std::shared_ptr<QuadratureData<T>> newQuadratureDataBuffer(const std::string& mesh_tag, int order, int dim,
+  static std::shared_ptr<QuadratureData<T>> newQuadratureDataBuffer(const Domain& domain, int order, int dim,
                                                                     T initial_state)
   {
-    SLIC_ERROR_ROOT_IF(!hasMesh(mesh_tag), axom::fmt::format("Mesh tag '{}' not found in the data store", mesh_tag));
-
     int Q = order + 1;
 
-    std::array<uint32_t, mfem::Geometry::NUM_GEOMETRIES> elems = geometry_counts(mesh(mesh_tag));
+    std::array<uint32_t, mfem::Geometry::NUM_GEOMETRIES> elems = geometry_counts(domain);
     std::array<uint32_t, mfem::Geometry::NUM_GEOMETRIES> qpts_per_elem{};
 
     std::vector<mfem::Geometry::Type> geometries;

src/serac/physics/tests/CMakeLists.txt

@@ -18,6 +18,7 @@ set(physics_serial_test_sources
     dynamic_solid_adjoint.cpp
     quasistatic_solid_adjoint.cpp
     finite_element_vector_set_over_domain.cpp
+    solid_multi_material.cpp
     )
 
 serac_add_tests(SOURCES       ${physics_serial_test_sources}