Simplify displacement BC interface

examples/buckling/cylinder.cpp

@@ -32,6 +32,8 @@
 
 #include "mfem.hpp"
 
+#include "serac/numerics/functional/domain.hpp"
+#include "serac/physics/boundary_conditions/components.hpp"
 #include "serac/physics/solid_mechanics_contact.hpp"
 #include "serac/infrastructure/terminator.hpp"
 #include "serac/mesh/mesh_utils.hpp"
@@ -123,15 +125,18 @@ int main(int argc, char* argv[])
   serac::StateManager::initialize(datastore, name + "_data");
 
   // Create and refine mesh
-  std::string filename = SERAC_REPO_DIR "/data/meshes/hollow-cylinder.mesh";
-  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};
-  std::set<int> xpos{3};
-  std::set<int> bottom{1};
-  std::set<int> top{4};
+  std::string filename     = SERAC_REPO_DIR "/data/meshes/hollow-cylinder.mesh";
+  auto        mesh         = serac::buildMeshFromFile(filename);
+  auto        refined_mesh = mesh::refineAndDistribute(std::move(mesh), serial_refinement, parallel_refinement);
+  auto&       pmesh        = serac::StateManager::setMesh(std::move(refined_mesh), mesh_tag);
+
+  // Surfaces for boundary conditions
+  constexpr int xneg_attr{2};
+  constexpr int xpos_attr{3};
+  auto          xneg   = serac::Domain::ofBoundaryElements(pmesh, serac::by_attr<dim>(xneg_attr));
+  auto          xpos   = serac::Domain::ofBoundaryElements(pmesh, serac::by_attr<dim>(xpos_attr));
+  auto          bottom = serac::Domain::ofBoundaryElements(pmesh, serac::by_attr<dim>(1));
+  auto          top    = serac::Domain::ofBoundaryElements(pmesh, serac::by_attr<dim>(4));
 
   // Create solver, either with or without contact
   std::unique_ptr<SolidMechanics<p, dim>> solid_solver;
@@ -145,14 +150,12 @@ int main(int argc, char* argv[])
                                           .type        = serac::ContactType::Frictionless,
                                           .penalty     = penalty};
     auto                  contact_interaction_id = 0;
-    solid_contact_solver->addContactInteraction(contact_interaction_id, xpos, xneg, contact_options);
+    solid_contact_solver->addContactInteraction(contact_interaction_id, {xpos_attr}, {xneg_attr}, contact_options);
     solid_solver = std::move(solid_contact_solver);
   } else {
     solid_solver = std::make_unique<serac::SolidMechanics<p, dim>>(
         nonlinear_options, linear_options, serac::solid_mechanics::default_quasistatic_options, name, mesh_tag);
-    auto domain = serac::Domain::ofBoundaryElements(
-        StateManager::mesh(mesh_tag), [&](std::vector<vec3>, int attr) { return xpos.find(attr) != xpos.end(); });
-    solid_solver->setPressure([&](auto&, double t) { return 0.01 * t; }, domain);
+    solid_solver->setPressure([&](auto&, double t) { return 0.01 * t; }, xpos);
   }
 
   // Define a Neo-Hookean material
@@ -164,19 +167,17 @@ int main(int argc, char* argv[])
 
   // Set up essential boundary conditions
   // Bottom of cylinder is fixed
-  auto clamp = [](const mfem::Vector&, mfem::Vector& u) {
-    u.SetSize(dim);
-    u = 0.0;
-  };
-  solid_solver->setDisplacementBCs(bottom, clamp);
+  solid_solver->setFixedBCs(bottom);
 
   // Top of cylinder has prescribed displacement of magnitude in x-z direction
-  auto compress = [&](const mfem::Vector&, double t, mfem::Vector& u) {
-    u.SetSize(dim);
-    u    = 0.0;
-    u(0) = u(2) = -1.5 / std::sqrt(2.0) * t;
+  auto compress = [&](const serac::tensor<double, dim>, double t) {
+    serac::tensor<double, dim> u{};
+    u[0] = u[2] = -1.5 / std::sqrt(2.0) * t;
+    return u;
   };
-  solid_solver->setDisplacementBCs(top, compress);
+  solid_solver->setDisplacementBCs(compress, top, Component::X + Component::Z);
+  solid_solver->setDisplacementBCs(compress, top,
+                                   Component::Y);  // BT: Would it be better to leave this component free?
 
   // Finalize the data structures
   solid_solver->completeSetup();

examples/contact/beam_bending.cpp

@@ -38,6 +38,10 @@ int main(int argc, char* argv[])
   auto  mesh  = serac::mesh::refineAndDistribute(serac::buildMeshFromFile(filename), 2, 0);
   auto& pmesh = serac::StateManager::setMesh(std::move(mesh), "beam_mesh");
 
+  // create boundary domains for boundary conditions
+  auto support                      = serac::Domain::ofBoundaryElements(pmesh, serac::by_attr<dim>(1));
+  auto applied_displacement_surface = serac::Domain::ofBoundaryElements(pmesh, serac::by_attr<dim>(6));
+
   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.");
@@ -82,15 +86,15 @@ int main(int argc, char* argv[])
   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) {
-    u.SetSize(dim);
-    u = 0.0;
-  });
-  solid_solver.setDisplacementBCs({6}, [](const mfem::Vector&, double t, mfem::Vector& u) {
-    u.SetSize(dim);
-    u    = 0.0;
+  solid_solver.setFixedBCs(support);
+
+  auto applied_displacement = [](serac::tensor<double, dim>, double t) {
+    serac::tensor<double, dim> u{};
     u[2] = -0.05 * t;
-  });
+    return u;
+  };
+
+  solid_solver.setDisplacementBCs(applied_displacement, applied_displacement_surface);
 
   // Add the contact interaction
   auto          contact_interaction_id = 0;

examples/contact/ironing.cpp

@@ -14,6 +14,7 @@
 #include "serac/physics/solid_mechanics_contact.hpp"
 #include "serac/infrastructure/terminator.hpp"
 #include "serac/mesh/mesh_utils.hpp"
+#include "serac/numerics/functional/domain.hpp"
 #include "serac/physics/state/state_manager.hpp"
 #include "serac/physics/materials/parameterized_solid_material.hpp"
 #include "serac/serac_config.hpp"
@@ -83,21 +84,22 @@ int main(int argc, char* argv[])
   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) {
-    u.SetSize(dim);
-    u = 0.0;
-  });
-  solid_solver.setDisplacementBCs({12}, [](const mfem::Vector&, double t, mfem::Vector& u) {
-    constexpr double init_steps = 2.0;
-    u.SetSize(dim);
-    u = 0.0;
+  serac::Domain bottom_of_substrate = serac::Domain::ofBoundaryElements(pmesh, serac::by_attr<dim>(5));
+  solid_solver.setFixedBCs(bottom_of_substrate);
+
+  serac::Domain top_of_indenter      = serac::Domain::ofBoundaryElements(pmesh, serac::by_attr<dim>(12));
+  auto          applied_displacement = [](serac::tensor<double, dim>, double t) {
+    constexpr double           init_steps = 2.0;
+    serac::tensor<double, dim> u{};
     if (t <= init_steps + 1.0e-12) {
       u[2] = -t * 0.3 / init_steps;
     } else {
       u[0] = -(t - init_steps) * 0.25;
       u[2] = -0.3;
     }
-  });
+    return u;
+  };
+  solid_solver.setDisplacementBCs(applied_displacement, top_of_indenter);
 
   // Add the contact interaction
   auto          contact_interaction_id = 0;

examples/contact/sphere.cpp

@@ -14,6 +14,7 @@
 
 #include "mfem.hpp"
 
+#include "serac/numerics/functional/domain.hpp"
 #include "serac/physics/solid_mechanics_contact.hpp"
 #include "serac/infrastructure/terminator.hpp"
 #include "serac/mesh/mesh_utils.hpp"
@@ -54,6 +55,9 @@ int main(int argc, char* argv[])
   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");
 
+  auto fixed_boundary = serac::Domain::ofBoundaryElements(pmesh, serac::by_attr<dim>(3));
+  auto driven_surface = serac::Domain::ofBoundaryElements(pmesh, serac::by_attr<dim>(12));
+
   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.");
@@ -79,13 +83,10 @@ int main(int argc, char* argv[])
   solid_solver.setMaterial(mat, whole_mesh);
 
   // Pass the BC information to the solver object
-  solid_solver.setDisplacementBCs({3}, [](const mfem::Vector&, mfem::Vector& u) {
-    u.SetSize(dim);
-    u = 0.0;
-  });
-  solid_solver.setDisplacementBCs({12}, [](const mfem::Vector& x, double t, mfem::Vector& u) {
-    u.SetSize(dim);
-    u = 0.0;
+  solid_solver.setFixedBCs(fixed_boundary);
+
+  auto applied_displacement = [](serac::tensor<double, dim> x, double t) {
+    serac::tensor<double, dim> u{};
     if (t <= 3.0 + 1.0e-12) {
       u[2] = -t * 0.02;
     } else {
@@ -95,7 +96,10 @@ int main(int argc, char* argv[])
           std::sin(M_PI / 40.0 * (t - 3.0)) * (x[0] - 0.5) + (std::cos(M_PI / 40.0 * (t - 3.0)) - 1.0) * (x[1] - 0.5);
       u[2] = -0.06;
     }
-  });
+    return u;
+  };
+
+  solid_solver.setDisplacementBCs(applied_displacement, driven_surface);
 
   // Add the contact interaction
   auto          contact_interaction_id = 0;

examples/contact/twist.cpp

@@ -13,6 +13,7 @@
 
 #include "mfem.hpp"
 
+#include "serac/numerics/functional/domain.hpp"
 #include "serac/physics/solid_mechanics_contact.hpp"
 #include "serac/infrastructure/terminator.hpp"
 #include "serac/mesh/mesh_utils.hpp"
@@ -40,6 +41,9 @@ int main(int argc, char* argv[])
   auto  mesh  = serac::mesh::refineAndDistribute(serac::buildMeshFromFile(filename), 3, 0);
   auto& pmesh = serac::StateManager::setMesh(std::move(mesh), "twist_mesh");
 
+  auto fixed_surface  = serac::Domain::ofBoundaryElements(pmesh, serac::by_attr<dim>(3));
+  auto driven_surface = serac::Domain::ofBoundaryElements(pmesh, serac::by_attr<dim>(6));
+
   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.");
@@ -65,13 +69,10 @@ int main(int argc, char* argv[])
   solid_solver.setMaterial(mat, whole_mesh);
 
   // Pass the BC information to the solver object
-  solid_solver.setDisplacementBCs({3}, [](const mfem::Vector&, mfem::Vector& u) {
-    u.SetSize(dim);
-    u = 0.0;
-  });
-  solid_solver.setDisplacementBCs({6}, [](const mfem::Vector& x, double t, mfem::Vector& u) {
-    u.SetSize(dim);
-    u = 0.0;
+  solid_solver.setFixedBCs(fixed_surface);
+
+  auto applied_displacement = [](serac::tensor<double, dim> x, double t) {
+    serac::tensor<double, dim> u{};
     if (t <= 3.0 + 1.0e-12) {
       u[2] = -t * 0.02;
     } else {
@@ -81,7 +82,10 @@ int main(int argc, char* argv[])
           std::sin(M_PI / 40.0 * (t - 3.0)) * (x[0] - 0.5) + (std::cos(M_PI / 40.0 * (t - 3.0)) - 1.0) * (x[1] - 0.5);
       u[2] = -0.06;
     }
-  });
+    return u;
+  };
+
+  solid_solver.setDisplacementBCs(applied_displacement, driven_surface);
 
   // Add the contact interaction
   auto          contact_interaction_id = 0;

src/serac/numerics/functional/domain.cpp

@@ -77,7 +77,7 @@ static Domain domain_of_edges(const mesh_t& mesh, std::function<T> predicate)
 
     if constexpr (d == 2) {
       int bdr_id = edge_id_to_bdr_id[i];
-      int attr   = (bdr_id > 0) ? mesh.GetBdrAttribute(bdr_id) : -1;
+      int attr   = (bdr_id >= 0) ? mesh.GetBdrAttribute(bdr_id) : -1;
       if (predicate(x, attr)) {
         output.addElement(i, i, mfem::Geometry::SEGMENT);
       }
@@ -367,6 +367,94 @@ Domain Domain::ofBoundaryElements(const mesh_t& mesh, std::function<bool(std::ve
   return domain_of_boundary_elems<3>(mesh, func);
 }
 
+/**
+ * @brief  Get local dofs that are part of a domain, but are owned by a neighboring MPI rank
+ *
+ *  This is necessary for situations like this:
+ *  Mesh before parallel partition:
+ *      3 *--------* 2
+ *        |       /|
+ *        |      / |
+ *        |     /  |
+ *        |    /   |     <-- Edge we want in Domain
+ *        |   /    |
+ *        |  /     |
+ *        | /      |
+ *        |/       |
+ * Node 0 *--------* 1
+ *
+ * Possible mesh after partition into two ranks:
+ *
+ *    RANK 0           RANK 1
+ *
+ * 3 *--------* 2         2 o
+ *   |       /             /|
+ *   |      /             / |
+ *   |     /             /  |
+ *   |    /             /   |  <-- Edge we want in Domain
+ *   |   /             /    |
+ *   |  /             /     |
+ *   | /             /      |
+ *   |/             /       |
+ * 0 *           0 o--------* 1
+ *
+ *   *: locally owned node
+ *   o: node owned by a neighbor rank
+ *
+ * We create a domain containing the right vertical edge, and then ask for its local dofs.
+ * The dof list returned for Rank 1 will be correct, containing the local indices for nodes
+ * 1 and 2. However, the dof list on rank 0 will not be correct without parallel communication.
+ * It will see that it doesn't own the edge in question, so when it then goes to fetch the
+ * local dofs on the domain, it will be an empty list.
+ *
+ * This function corrects for that, flagging the dofs we want on the domain on each rank
+ * (using the local_dof_ids list), and then exchanging this info with
+ * neighboring ranks, so that rank 0 will be told that its local dof for node 2 should be
+ * added to the list of dofs on the domain.
+ *
+ * Before findDomainDofsOnNeighborRanks():
+ * dof list on Rank 0: {}
+ * dof list on Rank 1: {1, 2}
+ *
+ * After:
+ * dof list on Rank 0: {2}
+ * dof list on Rank 1: {1, 2}
+ *
+ * Note: the sets will actually contain the ldof indices corresponding to the global (tdof)
+ * indices in the sets above.
+ *
+ * This function operates on the local_dof_ids data in place.
+ */
+void findDomainDofsOnNeighborRanks(const serac::fes_t* fes, mfem::Array<int>& local_dof_ids)
+{
+  auto par_fes = dynamic_cast<const mfem::ParFiniteElementSpace*>(fes);
+  // There's no work to do unless the finite element space really is parallel
+  if (par_fes) {
+    // As far as I can tell, the parallel communication in mfem only works with
+    // vector field dof indexing. So we need to get the parallel-correct scalar
+    // dof ids, we do the following:
+    // (1) transform scalar ldof ids to vector ldof ids,
+    // (2) transform the vector ldof ids into a boolean "marker" ldof field
+    // (3) do our parallel sync, which applies an OR logic operator to the
+    //     boolean fields from all ranks at each dof
+    // (4) get the ldof indices of the TRUE values of the parallel-correct
+    //     boolean ldof field
+    // (5) transform the parallel-correct vector ldof ids back to scalar dof ids.
+    fes->DofsToVDofs(0, local_dof_ids);
+
+    mfem::Array<int> local_dof_markers;
+    mfem::FiniteElementSpace::ListToMarker(local_dof_ids, par_fes->GetVSize(), local_dof_markers, 1);
+
+    par_fes->Synchronize(local_dof_markers);
+
+    mfem::FiniteElementSpace::MarkerToList(local_dof_markers, local_dof_ids);
+
+    for (int i = 0; i < local_dof_ids.Size(); i++) {
+      local_dof_ids[i] = par_fes->VDofToDof(local_dof_ids[i]);
+    }
+  }
+}
+
 mfem::Array<int> Domain::dof_list(const serac::fes_t* fes) const
 {
   std::set<int>    dof_ids;
@@ -432,6 +520,8 @@ mfem::Array<int> Domain::dof_list(const serac::fes_t* fes) const
     uniq_dof_ids[i++] = id;
   }
 
+  findDomainDofsOnNeighborRanks(fes, uniq_dof_ids);
+
   return uniq_dof_ids;
 }
 

src/serac/numerics/functional/domain.hpp

@@ -268,7 +268,14 @@ Domain operator-(const Domain& a, const Domain& b);
 template <int dim>
 inline auto by_attr(int value)
 {
-  return [value](std::vector<tensor<double, dim> >, int attr) { return attr == value; };
+  return [value](std::vector<tensor<double, dim>>, int attr) { return value == attr; };
+}
+
+/// @brief convenience predicate for creating domains by a set of attributes
+template <int dim>
+inline auto by_attr(std::set<int> values)
+{
+  return [values](std::vector<tensor<double, dim>>, int attr) { return values.find(attr) != values.end(); };
 }
 
 /**
@@ -292,7 +299,7 @@ inline std::array<uint32_t, mfem::Geometry::NUM_GEOMETRIES> geometry_counts(cons
  * @brief convenience function for computing the arithmetic mean of some list of vectors
  */
 template <int dim>
-inline tensor<double, dim> average(std::vector<tensor<double, dim> >& positions)
+inline tensor<double, dim> average(std::vector<tensor<double, dim>>& positions)
 {
   tensor<double, dim> total{};
   for (auto x : positions) {

src/serac/numerics/functional/tests/domain_tests.cpp

@@ -103,6 +103,7 @@ TEST(domain, of_edges)
 
     // check that by_attr compiles
     Domain d4 = Domain::ofEdges(mesh, by_attr<dim>(3));
+    EXPECT_EQ(d4.mfem_edge_ids_.size(), 16);
 
     Domain d5 = Domain::ofBoundaryElements(mesh, [](std::vector<vec2>, int) { return true; });
     EXPECT_EQ(d5.edge_ids_.size(), 18);  // 1x8 row of quads has 18 boundary edges