Merge branch 'main' into dokken/fix-em-demos

cpp/dolfinx/fem/utils.cpp

@@ -135,11 +135,12 @@ std::vector<std::int32_t> fem::compute_integration_domains(
  const int tdim = topology.dim();
  if ((integral_type == IntegralType::cell ? tdim : tdim - 1) != dim)
  {
- throw std::runtime_error("Invalid MeshTags dimension: "
+ throw std::runtime_error("Invalid mesh entity dimension: "
  + std::to_string(dim));
  }
 
  {
+ // Create span of the owned entities (leaves off any ghosts)
  assert(topology.index_map(dim));
  auto it1 = std::ranges::lower_bound(entities,
  topology.index_map(dim)->size_local());
@@ -150,15 +151,19 @@ std::vector<std::int32_t> fem::compute_integration_domains(
  switch (integral_type)
  {
  case IntegralType::cell:
+ {
  entity_data.insert(entity_data.begin(), entities.begin(), entities.end());
  break;
+ }
  default:
+ {
  auto f_to_c = topology.connectivity(tdim - 1, tdim);
  if (!f_to_c)
  {
  throw std::runtime_error(
  "Topology facet-to-cell connectivity has not been computed.");
  }
+
  auto c_to_f = topology.connectivity(tdim, tdim - 1);
  if (!c_to_f)
  {
@@ -186,9 +191,18 @@ std::vector<std::int32_t> fem::compute_integration_domains(
  break;
  case IntegralType::interior_facet:
  {
- for (std::size_t j = 0; j < entities.size(); ++j)
+ // Create indicator for interprocess facets
+ assert(topology.index_map(tdim - 1));
+ const std::vector<std::int32_t>& interprocess_facets
+ = topology.interprocess_facets();
+ std::vector<std::int8_t> interprocess_marker(
+ topology.index_map(tdim - 1)->size_local()
+ + topology.index_map(tdim - 1)->num_ghosts(),
+ 0);
+ std::ranges::for_each(interprocess_facets, [&interprocess_marker](auto f)
+ { interprocess_marker[f] = 1; });
+ for (auto f : entities)
  {
- const std::int32_t f = entities[j];
  if (f_to_c->num_links(f) == 2)
  {
  // Get the facet as a pair of (cell, local facet) pairs, one
@@ -197,6 +211,12 @@ std::vector<std::int32_t> fem::compute_integration_domains(
  = impl::get_cell_facet_pairs<2>(f, f_to_c->links(f), *c_to_f);
  entity_data.insert(entity_data.end(), facets.begin(), facets.end());
  }
+ else if (interprocess_marker[f])
+ {
+ throw std::runtime_error(
+ "Cannot compute interior facet integral over interprocess facet. "
+ "Use \"shared facet\" ghost mode when creating the mesh.");
+ }
  }
  }
  break;
@@ -205,6 +225,8 @@ std::vector<std::int32_t> fem::compute_integration_domains(
  "Cannot compute integration domains. Integral type not supported.");
  }
  }
+ }
+
  return entity_data;
 }
 //-----------------------------------------------------------------------------

cpp/dolfinx/fem/utils.h

@@ -573,6 +573,21 @@ Form<T, U> create_form_factory(
  num_integrals_type[interior_facet]);
  auto itg = integrals.insert({IntegralType::interior_facet, {}});
  auto sd = subdomains.find(IntegralType::interior_facet);
+
+ // Create indicator for interprocess facets
+ std::vector<std::int8_t> interprocess_marker;
+ if (num_integrals_type[interior_facet] > 0)
+ {
+ assert(topology->index_map(tdim - 1));
+ const std::vector<std::int32_t>& interprocess_facets
+ = topology->interprocess_facets();
+ std::int32_t num_facets = topology->index_map(tdim - 1)->size_local()
+ + topology->index_map(tdim - 1)->num_ghosts();
+ interprocess_marker.resize(num_facets, 0);
+ std::ranges::for_each(interprocess_facets, [&interprocess_marker](auto f)
+ { interprocess_marker[f] = 1; });
+ }
+
  for (int i = 0; i < num_integrals_type[interior_facet]; ++i)
  {
  const int id = ids[i];
@@ -631,6 +646,13 @@ Form<T, U> create_form_factory(
  default_facets_int.insert(default_facets_int.end(), pairs.begin(),
  pairs.end());
  }
+ else if (interprocess_marker[f])
+ {
+ throw std::runtime_error(
+ "Cannot compute interior facet integral over interprocess "
+ "facet. Please use ghost mode shared facet when creating the "
+ "mesh");
+ }
  }
  itg.first->second.emplace_back(id, k, default_facets_int,
  active_coeffs);

cpp/dolfinx/mesh/Topology.h

@@ -191,13 +191,19 @@ class Topology
  /// @brief Compute entity permutations and reflections.
  void create_entity_permutations();
 
- /// @brief List of inter-process facets, if facet topology has been
- /// computed.
+ /// @brief List of inter-process facets.
+ ///
+ /// "Inter-process" facets are facets that are connected (1) to a cell
+ /// that is owned by the calling process (rank) and (2) to a cell that
+ /// is owned by another process.
+ ///
+ /// @pre Inter-process facets are available only if facet topology has
+ /// been computed.
  const std::vector<std::int32_t>& interprocess_facets() const;
 
  /// @brief List of inter-process facets, if facet topology has been
- /// computed, for the facet type in `Topology::entity_types` identified by
- /// index
+ /// computed, for the facet type in `Topology::entity_types`
+ /// identified by index.
  /// @param index Index of facet type
  const std::vector<std::int32_t>& interprocess_facets(std::int8_t index) const;
 

cpp/dolfinx/mesh/topologycomputation.cpp

@@ -95,7 +95,7 @@ int get_ownership(const U& processes, const V& vertices)
 /// @param[in] num_entities_per_cell Number of entities per cell
 /// @param[in] entity_index Initial numbering for each row in
 /// entity_list
-/// @returns Local indices and index map
+/// @returns Local indices, the index map and shared entities
 std::tuple<std::vector<int>, common::IndexMap, std::vector<std::int32_t>>
 get_local_indexing(MPI_Comm comm, const common::IndexMap& vertex_map,
  std::span<const std::int32_t> entity_list,
@@ -804,13 +804,13 @@ mesh::compute_entities(MPI_Comm comm, const Topology& topology, int dim,
  cell_lists[i] = {cell_types[i], cells, cell_map};
  }
 
- auto [d0, d1, im, interprocess_facets] = compute_entities_by_key_matching(
+ auto [d0, d1, im, interprocess_entities] = compute_entities_by_key_matching(
  comm, cell_lists, *vertex_map, entity_type, dim);
 
  return {d0,
  std::make_shared<graph::AdjacencyList<std::int32_t>>(std::move(d1)),
  std::make_shared<common::IndexMap>(std::move(im)),
- std::move(interprocess_facets)};
+ std::move(interprocess_entities)};
 }
 //-----------------------------------------------------------------------------
 std::array<std::shared_ptr<graph::AdjacencyList<std::int32_t>>, 2>

cpp/dolfinx/mesh/topologycomputation.h

@@ -41,9 +41,9 @@ class Topology;
 /// `Topology::entity_types(dim)`.
 /// @return Tuple of (cell-entity connectivity, entity-vertex
 /// connectivity, index map, list of interprocess entities).
-/// Interprocess entities lie on the "true" boundary between owned cells of each
-/// process. If the entities already exist, then {nullptr, nullptr, nullptr,
-/// std::vector()} is returned.
+/// Interprocess entities lie on the "true" boundary between owned cells
+/// of each process. If the entities already exists, then {nullptr,
+/// nullptr, nullptr, std::vector()} is returned.
 std::tuple<std::vector<std::shared_ptr<graph::AdjacencyList<std::int32_t>>>,
  std::shared_ptr<graph::AdjacencyList<std::int32_t>>,
  std::shared_ptr<common::IndexMap>, std::vector<std::int32_t>>

python/test/unit/fem/test_assemble_domains.py

@@ -49,12 +49,7 @@ def create_cell_meshtags_from_entities(mesh: Mesh, dim: int, cells: np.ndarray,
  reason="Unghosted interior facets fail in parallel",
  ),
  ),
- pytest.param(
- GhostMode.shared_facet,
- marks=pytest.mark.skipif(
- condition=MPI.COMM_WORLD.size == 1, reason="Shared ghost modes fail in serial"
- ),
- ),
+ GhostMode.shared_facet,
  ],
 )