add support for interior face integrals

cmake/SeracBasics.cmake

@@ -43,7 +43,9 @@ endif()
 
 option(SERAC_ENABLE_PROFILING "Enable profiling functionality" OFF)
 
-cmake_dependent_option(SERAC_ENABLE_BENCHMARKS "Enable benchmark executables" ON "ENABLE_BENCHMARKS" OFF)
+if (ENABLE_BENCHMARKS)
+    set(SERAC_ENABLE_BENCHMARKS ON)
+endif()
 
 # User turned on benchmarking but explicitly turned off profiling. Error out.
 if ((ENABLE_BENCHMARKS OR SERAC_ENABLE_BENCHMARKS) AND NOT SERAC_ENABLE_PROFILING)

examples/buckling/cylinder.cpp

@@ -124,9 +124,8 @@ int main(int argc, char* argv[])
 
   // Create and refine mesh
   std::string filename = SERAC_REPO_DIR "/data/meshes/hollow-cylinder.mesh";
-  auto        mesh     = serac::buildMeshFromFile(filename);
-  auto        pmesh    = mesh::refineAndDistribute(std::move(mesh), serial_refinement, parallel_refinement);
-  serac::StateManager::setMesh(std::move(pmesh), mesh_tag);
+  auto  mesh  = mesh::refineAndDistribute(serac::buildMeshFromFile(filename), serial_refinement, parallel_refinement);
+  auto& pmesh = serac::StateManager::setMesh(std::move(mesh), mesh_tag);
 
   // Surface attributes for boundary conditions
   std::set<int> xneg{2};
@@ -160,8 +159,8 @@ int main(int argc, char* argv[])
   auto                        lambda = 1.0;
   auto                        G      = 0.1;
   solid_mechanics::NeoHookean mat{.density = 1.0, .K = (3 * lambda + 2 * G) / 3, .G = G};
-
-  solid_solver->setMaterial(mat);
+  Domain                      whole_mesh = EntireDomain(pmesh);
+  solid_solver->setMaterial(mat, whole_mesh);
 
   // Set up essential boundary conditions
   // Bottom of cylinder is fixed

examples/contact/beam_bending.cpp

@@ -35,8 +35,8 @@ int main(int argc, char* argv[])
   // Construct the appropriate dimension mesh and give it to the data store
   std::string filename = SERAC_REPO_DIR "/data/meshes/beam-hex-with-contact-block.mesh";
 
-  auto mesh = serac::mesh::refineAndDistribute(serac::buildMeshFromFile(filename), 2, 0);
-  serac::StateManager::setMesh(std::move(mesh), "beam_mesh");
+  auto  mesh  = serac::mesh::refineAndDistribute(serac::buildMeshFromFile(filename), 2, 0);
+  auto& pmesh = serac::StateManager::setMesh(std::move(mesh), "beam_mesh");
 
   serac::LinearSolverOptions linear_options{.linear_solver = serac::LinearSolver::Strumpack, .print_level = 1};
 #ifndef MFEM_USE_STRUMPACK
@@ -78,7 +78,8 @@ int main(int argc, char* argv[])
   solid_solver.setParameter(1, G_field);
 
   serac::solid_mechanics::ParameterizedNeoHookeanSolid mat{1.0, 0.0, 0.0};
-  solid_solver.setMaterial(serac::DependsOn<0, 1>{}, mat);
+  serac::Domain                                        whole_mesh = serac::EntireDomain(pmesh);
+  solid_solver.setMaterial(serac::DependsOn<0, 1>{}, mat, whole_mesh);
 
   // Pass the BC information to the solver object
   solid_solver.setDisplacementBCs({1}, [](const mfem::Vector&, mfem::Vector& u) {
@@ -117,4 +118,4 @@ int main(int argc, char* argv[])
   serac::exitGracefully();
 
   return 0;
-}
+}

examples/contact/ironing.cpp

@@ -35,10 +35,11 @@ int main(int argc, char* argv[])
   // Construct the appropriate dimension mesh and give it to the data store
   std::string filename = SERAC_REPO_DIR "/data/meshes/ironing.mesh";
 
-  auto mesh = serac::mesh::refineAndDistribute(serac::buildMeshFromFile(filename), 2, 0);
-  serac::StateManager::setMesh(std::move(mesh), "ironing_mesh");
+  auto  mesh  = serac::mesh::refineAndDistribute(serac::buildMeshFromFile(filename), 2, 0);
+  auto& pmesh = serac::StateManager::setMesh(std::move(mesh), "ironing_mesh");
 
   serac::LinearSolverOptions linear_options{.linear_solver = serac::LinearSolver::Strumpack, .print_level = 1};
+
 #ifndef MFEM_USE_STRUMPACK
   SLIC_INFO_ROOT("Contact requires MFEM built with strumpack.");
   return 1;
@@ -78,7 +79,8 @@ int main(int argc, char* argv[])
   solid_solver.setParameter(1, G_field);
 
   serac::solid_mechanics::ParameterizedNeoHookeanSolid mat{1.0, 0.0, 0.0};
-  solid_solver.setMaterial(serac::DependsOn<0, 1>{}, mat);
+  serac::Domain                                        whole_mesh = serac::EntireDomain(pmesh);
+  solid_solver.setMaterial(serac::DependsOn<0, 1>{}, mat, whole_mesh);
 
   // Pass the BC information to the solver object
   solid_solver.setDisplacementBCs({5}, [](const mfem::Vector&, mfem::Vector& u) {

examples/contact/sphere.cpp

@@ -51,8 +51,8 @@ int main(int argc, char* argv[])
   cube_mesh.SetCurvature(p);
 
   std::vector<mfem::Mesh*> mesh_ptrs{&ball_mesh, &cube_mesh};
-  auto mesh = serac::mesh::refineAndDistribute(mfem::Mesh(mesh_ptrs.data(), static_cast<int>(mesh_ptrs.size())), 0, 0);
-  serac::StateManager::setMesh(std::move(mesh), "sphere_mesh");
+  auto  mesh = serac::mesh::refineAndDistribute(mfem::Mesh(mesh_ptrs.data(), static_cast<int>(mesh_ptrs.size())), 0, 0);
+  auto& pmesh = serac::StateManager::setMesh(std::move(mesh), "sphere_mesh");
 
   serac::LinearSolverOptions linear_options{.linear_solver = serac::LinearSolver::Strumpack, .print_level = 1};
 #ifndef MFEM_USE_STRUMPACK
@@ -75,7 +75,8 @@ int main(int argc, char* argv[])
       nonlinear_options, linear_options, serac::solid_mechanics::default_quasistatic_options, name, "sphere_mesh");
 
   serac::solid_mechanics::NeoHookean mat{1.0, 10.0, 0.25};
-  solid_solver.setMaterial(mat);
+  serac::Domain                      whole_mesh = serac::EntireDomain(pmesh);
+  solid_solver.setMaterial(mat, whole_mesh);
 
   // Pass the BC information to the solver object
   solid_solver.setDisplacementBCs({3}, [](const mfem::Vector&, mfem::Vector& u) {
@@ -122,4 +123,4 @@ int main(int argc, char* argv[])
   serac::exitGracefully();
 
   return 0;
-}
+}

examples/contact/twist.cpp

@@ -37,8 +37,8 @@ int main(int argc, char* argv[])
   // Construct the appropriate dimension mesh and give it to the data store
   std::string filename = SERAC_REPO_DIR "/data/meshes/twohex_for_contact.mesh";
 
-  auto mesh = serac::mesh::refineAndDistribute(serac::buildMeshFromFile(filename), 3, 0);
-  serac::StateManager::setMesh(std::move(mesh), "twist_mesh");
+  auto  mesh  = serac::mesh::refineAndDistribute(serac::buildMeshFromFile(filename), 3, 0);
+  auto& pmesh = serac::StateManager::setMesh(std::move(mesh), "twist_mesh");
 
   serac::LinearSolverOptions linear_options{.linear_solver = serac::LinearSolver::Strumpack, .print_level = 1};
 #ifndef MFEM_USE_STRUMPACK
@@ -61,7 +61,8 @@ int main(int argc, char* argv[])
       nonlinear_options, linear_options, serac::solid_mechanics::default_quasistatic_options, name, "twist_mesh");
 
   serac::solid_mechanics::NeoHookean mat{1.0, 10.0, 10.0};
-  solid_solver.setMaterial(mat);
+  serac::Domain                      whole_mesh = serac::EntireDomain(pmesh);
+  solid_solver.setMaterial(mat, whole_mesh);
 
   // Pass the BC information to the solver object
   solid_solver.setDisplacementBCs({3}, [](const mfem::Vector&, mfem::Vector& u) {
@@ -108,4 +109,4 @@ int main(int argc, char* argv[])
   serac::exitGracefully();
 
   return 0;
-}
+}

examples/simple_conduction/without_input_file.cpp

@@ -38,7 +38,7 @@ int main(int argc, char* argv[])
 
   std::string mesh_tag{"mesh"};
 
-  serac::StateManager::setMesh(std::move(mesh), mesh_tag);
+  auto& pmesh = serac::StateManager::setMesh(std::move(mesh), mesh_tag);
   // _create_mesh_end
 
   // _create_module_start
@@ -54,7 +54,9 @@ int main(int argc, char* argv[])
   // _conductivity_start
   constexpr double                               kappa = 0.5;
   serac::heat_transfer::LinearIsotropicConductor mat(1.0, 1.0, kappa);
-  heat_transfer.setMaterial(mat);
+
+  serac::Domain whole_domain = serac::EntireDomain(pmesh);
+  heat_transfer.setMaterial(mat, whole_domain);
 
   // _conductivity_end
   // _bc_start

src/drivers/CMakeLists.txt

@@ -4,33 +4,33 @@
 #
 # SPDX-License-Identifier: (BSD-3-Clause)
 
-blt_add_executable( NAME        serac_driver
-                    SOURCES     serac.cpp
-                    DEPENDS_ON  serac_physics serac_mesh
-                    OUTPUT_NAME serac
-                    )
-
-if (SERAC_ENABLE_TESTS)
-    set(input_files_dir ${CMAKE_CURRENT_SOURCE_DIR}/../../data/input_files/tests)
-
-    # Run basic test for the Serac driver
-    blt_add_test(NAME          serac_driver_solid
-                 COMMAND       ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}/serac -o driver_solid -i ${input_files_dir}/solid/dyn_solve.lua
-                 NUM_MPI_TASKS 1 )
-
-    blt_add_test(NAME          serac_driver_heat_transfer
-                 COMMAND       ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}/serac -o driver_heat_transfer -i ${input_files_dir}/heat_transfer/static_solve.lua
-                 NUM_MPI_TASKS 1 )
-
-    blt_add_test(NAME          serac_driver_help
-                 COMMAND       ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}/serac --help 
-                 NUM_MPI_TASKS 1 )
-
-    blt_add_test(NAME          serac_driver_docs
-                 COMMAND       ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}/serac -o docs -d -i ${input_files_dir}/solid/qs_linear.lua
-                 NUM_MPI_TASKS 1 )
-endif()
-
-install( TARGETS serac_driver
-         RUNTIME DESTINATION bin
-         )
+#blt_add_executable( NAME        serac_driver
+#                    SOURCES     serac.cpp
+#                    DEPENDS_ON  serac_physics serac_mesh
+#                    OUTPUT_NAME serac
+#                    )
+#
+#if (SERAC_ENABLE_TESTS)
+#    set(input_files_dir ${CMAKE_CURRENT_SOURCE_DIR}/../../data/input_files/tests)
+#
+#    # Run basic test for the Serac driver
+#    blt_add_test(NAME          serac_driver_solid
+#                 COMMAND       ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}/serac -o driver_solid -i ${input_files_dir}/solid/dyn_solve.lua
+#                 NUM_MPI_TASKS 1 )
+#
+#    blt_add_test(NAME          serac_driver_heat_transfer
+#                 COMMAND       ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}/serac -o driver_heat_transfer -i ${input_files_dir}/heat_transfer/static_solve.lua
+#                 NUM_MPI_TASKS 1 )
+#
+#    blt_add_test(NAME          serac_driver_help
+#                 COMMAND       ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}/serac --help 
+#                 NUM_MPI_TASKS 1 )
+#
+#    blt_add_test(NAME          serac_driver_docs
+#                 COMMAND       ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}/serac -o docs -d -i ${input_files_dir}/solid/qs_linear.lua
+#                 NUM_MPI_TASKS 1 )
+#endif()
+#
+#install( TARGETS serac_driver
+#         RUNTIME DESTINATION bin
+#         )

src/serac/infrastructure/accelerator.hpp

@@ -118,6 +118,13 @@ void zero_out(axom::Array<T, dim, space>& arr)
 {
   std::memset(arr.data(), 0, static_cast<std::size_t>(arr.size()) * sizeof(T));
 }
+
+/// @brief set the contents of an array to zero, byte-wise
+template <typename T, int dim>
+void zero_out(axom::ArrayView<T, dim, axom::MemorySpace::Host>& arr)
+{
+  std::memset(arr.data(), 0, static_cast<std::size_t>(arr.size()) * sizeof(T));
+}
 #ifdef __CUDACC__
 /// @overload
 template <typename T, int dim>

src/serac/infrastructure/debug_print.hpp

@@ -61,12 +61,6 @@ std::ostream& operator<<(std::ostream& out, DoF dof)
   return out;
 }
 
-std::ostream& operator<<(std::ostream& out, serac::SignedIndex i)
-{
-  out << "{" << i.index_ << ", " << i.sign_ << "}";
-  return out;
-}
-
 /**
  * @brief write a 2D array of values out to file, in a space-separated format
  * @tparam T the type of each value in the array

src/serac/numerics/functional/CMakeLists.txt

@@ -10,7 +10,6 @@ set(functional_depends serac_mesh ${serac_device_depends})
 set(functional_headers
     differentiate_wrt.hpp
     boundary_integral_kernels.hpp
-    dof_numbering.hpp
     element_restriction.hpp
     geometry.hpp
     geometric_factors.hpp
@@ -22,6 +21,7 @@ set(functional_headers
     function_signature.hpp
     functional_qoi.inl
     integral.hpp
+    interior_face_integral_kernels.hpp
     isotropic_tensor.hpp
     polynomials.hpp
     quadrature.hpp
@@ -30,6 +30,7 @@ set(functional_headers
     tensor.hpp
     tuple.hpp
     tuple_tensor_dual_functions.hpp
+    typedefs.hpp
     )
 
 set(functional_sources 

src/serac/numerics/functional/boundary_integral_kernels.hpp

@@ -161,7 +161,7 @@ template <uint32_t differentiation_index, int Q, mfem::Geometry::Type geom, type
 void evaluation_kernel_impl(trial_element_type trial_elements, test_element, double t,
                             const std::vector<const double*>& inputs, double* outputs, const double* positions,
                             const double* jacobians, lambda_type qf, [[maybe_unused]] derivative_type* qf_derivatives,
-                            const int* elements, uint32_t num_elements, camp::int_seq<int, indices...>)
+                            uint32_t num_elements, camp::int_seq<int, indices...>)
 {
   // mfem provides this information as opaque arrays of doubles,
   // so we reinterpret the pointer with
@@ -185,7 +185,7 @@ void evaluation_kernel_impl(trial_element_type trial_elements, test_element, dou
 
     // batch-calculate values / derivatives of each trial space, at each quadrature point
     [[maybe_unused]] tuple qf_inputs = {promote_each_to_dual_when<indices == differentiation_index>(
-        get<indices>(trial_elements).interpolate(get<indices>(u)[elements[e]], rule))...};
+        get<indices>(trial_elements).interpolate(get<indices>(u)[e], rule))...};
 
     // (batch) evalute the q-function at each quadrature point
     auto qf_outputs = batch_apply_qf(qf, t, x_e, J_e, get<indices>(qf_inputs)...);
@@ -200,7 +200,7 @@ void evaluation_kernel_impl(trial_element_type trial_elements, test_element, dou
     }
 
     // (batch) integrate the material response against the test-space basis functions
-    test_element::integrate(get_value(qf_outputs), rule, &r[elements[e]]);
+    test_element::integrate(get_value(qf_outputs), rule, &r[e]);
   }
 }
 
@@ -250,8 +250,7 @@ SERAC_HOST_DEVICE auto batch_apply_chain_rule(derivative_type* qf_derivatives, c
  * @param[in] num_elements The number of elements in the mesh
  */
 template <int Q, mfem::Geometry::Type geom, typename test, typename trial, typename derivatives_type>
-void action_of_gradient_kernel(const double* dU, double* dR, derivatives_type* qf_derivatives, const int* elements,
-                               std::size_t num_elements)
+void action_of_gradient_kernel(const double* dU, double* dR, derivatives_type* qf_derivatives, std::size_t num_elements)
 {
   using test_element  = finite_element<geom, test>;
   using trial_element = finite_element<geom, trial>;
@@ -266,13 +265,13 @@ void action_of_gradient_kernel(const double* dU, double* dR, derivatives_type* q
   // for each element in the domain
   for (uint32_t e = 0; e < num_elements; e++) {
     // (batch) interpolate each quadrature point's value
-    auto qf_inputs = trial_element::interpolate(du[elements[e]], rule);
+    auto qf_inputs = trial_element::interpolate(du[e], rule);
 
     // (batch) evalute the q-function at each quadrature point
     auto qf_outputs = batch_apply_chain_rule(qf_derivatives + e * nqp, qf_inputs);
 
     // (batch) integrate the material response against the test-space basis functions
-    test_element::integrate(qf_outputs, rule, &dr[elements[e]]);
+    test_element::integrate(qf_outputs, rule, &dr[e]);
   }
 }
 
@@ -299,7 +298,7 @@ void action_of_gradient_kernel(const double* dU, double* dR, derivatives_type* q
  */
 template <mfem::Geometry::Type g, typename test, typename trial, int Q, typename derivatives_type>
 void element_gradient_kernel(ExecArrayView<double, 3, ExecutionSpace::CPU> dK, derivatives_type* qf_derivatives,
-                             const int* elements, std::size_t num_elements)
+                             std::size_t num_elements)
 {
   using test_element  = finite_element<g, test>;
   using trial_element = finite_element<g, trial>;
@@ -310,7 +309,7 @@ void element_gradient_kernel(ExecArrayView<double, 3, ExecutionSpace::CPU> dK, d
 
   // for each element in the domain
   for (uint32_t e = 0; e < num_elements; e++) {
-    auto* output_ptr = reinterpret_cast<typename test_element::dof_type*>(&dK(elements[e], 0, 0));
+    auto* output_ptr = reinterpret_cast<typename test_element::dof_type*>(&dK(e, 0, 0));
 
     tensor<derivatives_type, nquad> derivatives{};
     for (int q = 0; q < nquad; q++) {
@@ -327,35 +326,35 @@ void element_gradient_kernel(ExecArrayView<double, 3, ExecutionSpace::CPU> dK, d
 template <uint32_t wrt, int Q, mfem::Geometry::Type geom, typename signature, typename lambda_type,
           typename derivative_type>
 auto evaluation_kernel(signature s, lambda_type qf, const double* positions, const double* jacobians,
-                       std::shared_ptr<derivative_type> qf_derivatives, const int* elements, uint32_t num_elements)
+                       std::shared_ptr<derivative_type> qf_derivatives, uint32_t num_elements)
 {
   auto trial_elements = trial_elements_tuple<geom>(s);
   auto test_element   = get_test_element<geom>(s);
   return [=](double time, const std::vector<const double*>& inputs, double* outputs, bool /* update state */) {
     evaluation_kernel_impl<wrt, Q, geom>(trial_elements, test_element, time, inputs, outputs, positions, jacobians, qf,
-                                         qf_derivatives.get(), elements, num_elements, s.index_seq);
+                                         qf_derivatives.get(), num_elements, s.index_seq);
   };
 }
 
 template <int wrt, int Q, mfem::Geometry::Type geom, typename signature, typename derivative_type>
 std::function<void(const double*, double*)> jacobian_vector_product_kernel(
-    signature, std::shared_ptr<derivative_type> qf_derivatives, const int* elements, uint32_t num_elements)
+    signature, std::shared_ptr<derivative_type> qf_derivatives, uint32_t num_elements)
 {
   return [=](const double* du, double* dr) {
     using test_space  = typename signature::return_type;
     using trial_space = typename std::tuple_element<wrt, typename signature::parameter_types>::type;
-    action_of_gradient_kernel<Q, geom, test_space, trial_space>(du, dr, qf_derivatives.get(), elements, num_elements);
+    action_of_gradient_kernel<Q, geom, test_space, trial_space>(du, dr, qf_derivatives.get(), num_elements);
   };
 }
 
 template <int wrt, int Q, mfem::Geometry::Type geom, typename signature, typename derivative_type>
 std::function<void(ExecArrayView<double, 3, ExecutionSpace::CPU>)> element_gradient_kernel(
-    signature, std::shared_ptr<derivative_type> qf_derivatives, const int* elements, uint32_t num_elements)
+    signature, std::shared_ptr<derivative_type> qf_derivatives, uint32_t num_elements)
 {
   return [=](ExecArrayView<double, 3, ExecutionSpace::CPU> K_elem) {
     using test_space  = typename signature::return_type;
     using trial_space = typename std::tuple_element<wrt, typename signature::parameter_types>::type;
-    element_gradient_kernel<geom, test_space, trial_space, Q>(K_elem, qf_derivatives.get(), elements, num_elements);
+    element_gradient_kernel<geom, test_space, trial_space, Q>(K_elem, qf_derivatives.get(), num_elements);
   };
 }