Fix dofmap builder performance fixes (#4081)

* Dofmap builder performance fixes

* Fix

* Simplify

* Simplify

* Update cpp/dolfinx/mesh/cell_types.cpp

Co-authored-by: Jørgen Schartum Dokken <dokken92@gmail.com>

* Apply suggestions from code review

Co-authored-by: Paul T. Kühner <56360279+schnellerhase@users.noreply.github.com>

---------

Co-authored-by: Jørgen Schartum Dokken <dokken92@gmail.com>
Co-authored-by: Paul T. Kühner <56360279+schnellerhase@users.noreply.github.com>

Author: 3 people

cpp/dolfinx/fem/dofmapbuilder.cpp

@@ -155,7 +155,7 @@ build_basic_dofmaps(
     const std::vector<fem::ElementDofLayout>& element_dof_layouts)
 {
   // Start timer for dofmap initialization
-  common::Timer t0("Init dofmap from element dofmap");
+  common::Timer t0("Dofmap builder: init dofmap from element dofmap");
 
   // Topological dimension
   const std::size_t D = topology.dim();
@@ -208,7 +208,6 @@ build_basic_dofmaps(
             local_entity_offsets.push_back(
                 local_entity_offsets.back()
                 + num_entity_dofs * (im->size_local() + im->num_ghosts()));
-
             if (d < D
                 and !topology.connectivity({int(D), int(i)},
                                            {int(d), int(et_index)}))
@@ -263,6 +262,15 @@ build_basic_dofmaps(
     dofs[i].array.resize(num_cells * dofmap_width);
     spdlog::info("Cell type: {} dofmap: {}x{}", i, num_cells, dofmap_width);
 
+    std::vector<std::vector<mesh::CellType>> cell_entity_types(D + 1);
+    for (std::size_t d = 0; d < D + 1; ++d)
+    {
+      int entities_d = mesh::cell_num_entities(cell_type, d);
+      cell_entity_types[d].reserve(entities_d);
+      for (int e = 0; e < entities_d; ++e)
+        cell_entity_types[d].push_back(mesh::cell_entity_type(cell_type, d, e));
+    }
+
     std::int32_t dofmap_offset = 0;
     for (std::int32_t c = 0; c < num_cells; ++c)
     {
@@ -279,12 +287,12 @@ build_basic_dofmaps(
         std::size_t et = required_dim_et[k].second;
         mesh::CellType e_type = topology.entity_types(d)[et];
 
-        const std::vector<std::vector<int>>& e_dofs_d = entity_dofs[d];
-
         // Skip over undefined topology, e.g. quad facets of tetrahedra
         if (d < D
             and !topology.connectivity({int(D), int(i)}, {int(d), int(et)}))
+        {
           continue;
+        }
 
         // Iterate over each entity of current dimension d and type et
         std::span<const std::int32_t> c_to_e
@@ -293,13 +301,15 @@ build_basic_dofmaps(
                         ->links(c)
                   : std::span<const std::int32_t>(&c, 1);
 
+        const std::vector<std::vector<int>>& e_dofs_d = entity_dofs[d];
+        const std::vector<mesh::CellType>& e_types = cell_entity_types[d];
         int w = 0;
         for (std::size_t e = 0; e < e_dofs_d.size(); ++e)
         {
           // Skip entities of wrong type (e.g. for facets of prism)
           // Use separate connectivity index 'w' which only advances for
           // correct entities
-          if (mesh::cell_entity_type(cell_type, d, e) == e_type)
+          if (e_type == e_types[e])
           {
             const std::vector<int>& e_dofs_d_e = e_dofs_d[e];
             std::size_t num_entity_dofs = e_dofs_d_e.size();
@@ -347,7 +357,6 @@ build_basic_dofmaps(
     auto& map = topo_index_maps[k];
     assert(map);
     std::vector<std::int64_t> global_indices = map->global_indices();
-
     for (std::size_t e_index = 0; e_index < global_indices.size(); ++e_index)
     {
       auto e_index_global = global_indices[e_index];
@@ -369,16 +378,19 @@ build_basic_dofmaps(
 }
 //-----------------------------------------------------------------------------
 
-/// Compute re-ordering map from old local index to new local index. The
-/// M dofs owned by this process are reordered for locality and fill the
-/// positions [0, ..., M). Dof owned by another process are placed at
-/// the end, i.e. in the positions [M, ..., N), where N is the total
+/// @brief Compute re-ordering map from old local index to new local
+/// index.
+///
+/// The M dofs owned by this process are reordered for locality and fill
+/// the positions [0, ..., M). Dof owned by another process are placed
+/// at the end, i.e. in the positions [M, ..., N), where N is the total
 /// number of dofs on this process.
 ///
-/// @param [in] dofmaps The basic dofmap data in multiple dofmaps sharing the
-/// same range
-/// @param [in] dof_entity Map from dof index to (index_map, entity_index),
-/// where entity_index is the local mesh entity index in the given index_map
+/// @param [in] dofmaps The basic dofmap data in multiple dofmaps
+/// sharing the same range
+/// @param [in] dof_entity Map from dof index to (index_map,
+/// entity_index), where entity_index is the local mesh entity index in
+/// the given index_map
 /// @param [in] index_maps The set of IndexMaps, one for each
 /// topological entity type used in the dofmap. The location in this
 /// array is referred to by the first item in each entry of
@@ -394,7 +406,7 @@ std::pair<std::vector<std::int32_t>, std::int32_t> compute_reordering_map(
     const std::function<std::vector<int>(
         const graph::AdjacencyList<std::int32_t>&)>& reorder_fn)
 {
-  common::Timer t0("Compute dof reordering map");
+  common::Timer t0("Dofmap builder: compute dof reordering map");
 
   // Get mesh entity ownership offset for each IndexMap
   std::vector<std::int32_t> offset(index_maps.size(), -1);
@@ -427,10 +439,10 @@ std::pair<std::vector<std::int32_t>, std::int32_t> compute_reordering_map(
     }
   }
 
-  // Check for any -1's remaining in `original_to_contiguous` due to vertices
-  // on the process that don't belong to a cell. Determine if the dof is owned
-  // or a ghost and map to the ends of the owned and ghost "parts" of the
-  // contiguous array respectively.
+  // Check for any -1's remaining in `original_to_contiguous` due to
+  // vertices on the process that don't belong to a cell. Determine if
+  // the dof is owned or a ghost and map to the ends of the owned and
+  // ghost "parts" of the contiguous array respectively.
   for (std::size_t dof = 0; dof < original_to_contiguous.size(); ++dof)
   {
     if (original_to_contiguous[dof] == -1)
@@ -459,18 +471,18 @@ std::pair<std::vector<std::int32_t>, std::int32_t> compute_reordering_map(
 }
 //-----------------------------------------------------------------------------
 
-/// Get global indices for unowned dofs
-/// @param [in] index_maps Set of index maps corresponding to dofs in @p
-/// dof_entity, below.
-/// @param [in] num_owned The number of nodes owned by this process
-/// @param [in] process_offset The node offset for this process, i.e.
-/// the global index of owned node i is i + process_offset
-/// @param [in] global_indices_old The old global index of the old local
+/// @brief Get global indices for unowned dofs
+///
+/// @param[in] index_maps Set of index maps corresponding to dofs in
+/// `dof_entity`, below.
+/// @param[in] num_owned The number of nodes owned by this process
+/// @param[in] process_offset The node offset for this process, i.e. the
+/// global index of owned node i is i + process_offset
+/// @param[in] global_indices_old The old global index of the old local
 /// node i
-/// @param [in] old_to_new The old local index to new local index map
-/// @param [in] dof_entity The ith entry gives (index_map, local
-/// index) of the mesh entity to which node i (old local index) is
-/// associated.
+/// @param[in] old_to_new The old local index to new local index map
+/// @param[in] dof_entity The ith entry gives (index_map, local index)
+/// of the mesh entity to which node i (old local index) is associated.
 /// @returns The (0) global indices for unowned dofs, (1) owner rank of
 /// each unowned dof
 std::pair<std::vector<std::int64_t>, std::vector<int>> get_global_indices(
@@ -480,6 +492,8 @@ std::pair<std::vector<std::int64_t>, std::vector<int>> get_global_indices(
     const std::vector<std::int32_t>& old_to_new,
     const std::vector<std::pair<std::int8_t, std::int32_t>>& dof_entity)
 {
+  common::Timer t0("Dofmap builder: get dofmap global indices");
+
   assert(dof_entity.size() == global_indices_old.size());
 
   // Build list of flags for owned mesh entities that are shared, i.e.
@@ -630,7 +644,7 @@ fem::build_dofmap_data(
     const std::function<std::vector<int>(
         const graph::AdjacencyList<std::int32_t>&)>& reorder_fn)
 {
-  common::Timer t0("Build dofmap data");
+  common::Timer t0("Dofmap builder: build dofmap data");
 
   // Build a simple dofmap based on mesh entity numbering, returning (i)
   // a local dofmap, (ii) local-to-global map for dof indices, and (iii)

cpp/dolfinx/mesh/cell_types.cpp

@@ -1,4 +1,4 @@
-// Copyright (C) 2006-2019 Anders Logg and Garth N. Wells
+// Copyright (C) 2006-2026 Anders Logg and Garth N. Wells
 //
 // This file is part of DOLFINx (https://www.fenicsproject.org)
 //
@@ -61,53 +61,9 @@ mesh::CellType mesh::to_type(const std::string& cell)
     throw std::runtime_error("Unknown cell type (" + cell + ")");
 }
 //-----------------------------------------------------------------------------
-mesh::CellType mesh::cell_entity_type(CellType type, int d, int index)
-{
-  const int dim = cell_dim(type);
-  if (d == dim)
-    return type;
-  else if (d == 1)
-    return CellType::interval;
-  else if (d == (dim - 1))
-    return cell_facet_type(type, index);
-  else
-    return CellType::point;
-}
-//-----------------------------------------------------------------------------
-mesh::CellType mesh::cell_facet_type(CellType type, int index)
-{
-  switch (type)
-  {
-  case CellType::point:
-    return CellType::point;
-  case CellType::interval:
-    return CellType::point;
-  case CellType::triangle:
-    return CellType::interval;
-  case CellType::tetrahedron:
-    return CellType::triangle;
-  case CellType::quadrilateral:
-    return CellType::interval;
-  case CellType::pyramid:
-    if (index == 0)
-      return CellType::quadrilateral;
-    else
-      return CellType::triangle;
-  case CellType::prism:
-    if (index == 0 or index == 4)
-      return CellType::triangle;
-    else
-      return CellType::quadrilateral;
-  case CellType::hexahedron:
-    return CellType::quadrilateral;
-  default:
-    throw std::runtime_error("Unknown cell type.");
-  }
-}
-//-----------------------------------------------------------------------------
 graph::AdjacencyList<int> mesh::get_entity_vertices(CellType type, int dim)
 {
-  const std::vector<std::vector<int>> topology
+  std::vector<std::vector<int>> topology
       = basix::cell::topology(cell_type_to_basix_type(type))[dim];
   return graph::AdjacencyList<int>(topology);
 }
@@ -121,7 +77,7 @@ graph::AdjacencyList<int> mesh::get_sub_entities(CellType type, int dim0,
   else if (type == CellType::point)
     return graph::AdjacencyList<int>(0);
 
-  const std::vector<std::vector<std::vector<int>>> connectivity
+  std::vector<std::vector<std::vector<int>>> connectivity
       = basix::cell::sub_entity_connectivity(
           cell_type_to_basix_type(type))[dim0];
   std::vector<std::vector<int>> subset;
@@ -131,11 +87,6 @@ graph::AdjacencyList<int> mesh::get_sub_entities(CellType type, int dim0,
   return graph::AdjacencyList<int>(subset);
 }
 //-----------------------------------------------------------------------------
-int mesh::cell_dim(CellType type)
-{
-  return basix::cell::topological_dimension(cell_type_to_basix_type(type));
-}
-//-----------------------------------------------------------------------------
 int mesh::cell_num_entities(CellType type, int dim)
 {
   assert(dim <= 3);

cpp/dolfinx/mesh/cell_types.h

@@ -1,4 +1,4 @@
-// Copyright (C) 2019-2020 Garth N. Wells
+// Copyright (C) 2019-2026 Garth N. Wells
 //
 // This file is part of DOLFINx (https://www.fenicsproject.org)
 //
@@ -41,42 +41,110 @@ std::string to_string(CellType type);
 /// @return The cell type
 CellType to_type(const std::string& cell);
 
-/// Return type of cell for entity of dimension d at given entity index.
-CellType cell_entity_type(CellType type, int d, int index);
+/// Return topological dimension of cell type
+inline int cell_dim(CellType type)
+{
+  switch (type)
+  {
+  case CellType::point:
+    return 0;
+  case CellType::interval:
+    return 1;
+  case CellType::triangle:
+    return 2;
+  case CellType::quadrilateral:
+    return 2;
+  case CellType::tetrahedron:
+    return 3;
+  case CellType::hexahedron:
+    return 3;
+  case CellType::prism:
+    return 3;
+  case CellType::pyramid:
+    return 3;
+  default:
+    throw std::runtime_error("Unsupported cell type");
+  }
+}
+
+/// @brief Return facet type of cell.
+///
+/// For simplex and hypercube cell types, this is independent of the
+/// facet index, but for prism and pyramid, it can be triangle or
+/// quadrilateral.
+///
+/// @param[in] type Cell type.
+/// @param[in] index Facet index (relative to the cell).
+/// @return Type of facet for this cell at this index.
+inline CellType cell_facet_type(CellType type, int index)
+{
+  switch (type)
+  {
+  case CellType::point:
+    return CellType::point;
+  case CellType::interval:
+    return CellType::point;
+  case CellType::triangle:
+    return CellType::interval;
+  case CellType::tetrahedron:
+    return CellType::triangle;
+  case CellType::quadrilateral:
+    return CellType::interval;
+  case CellType::pyramid:
+    if (index == 0)
+      return CellType::quadrilateral;
+    else
+      return CellType::triangle;
+  case CellType::prism:
+    if (index == 0 or index == 4)
+      return CellType::triangle;
+    else
+      return CellType::quadrilateral;
+  case CellType::hexahedron:
+    return CellType::quadrilateral;
+  default:
+    throw std::runtime_error("Unknown cell type.");
+  }
+}
 
-/// Return facet type of cell
-/// For simplex and hypercube cell types, this is independent of the facet
-/// index, but for prism and pyramid, it can be triangle or quadrilateral.
-/// @param[in] type The cell type
-/// @param[in] index The facet index
-/// @return The type of facet for this cell at this index
-CellType cell_facet_type(CellType type, int index);
+/// Return type of cell for entity of dimension d at given entity index.
+inline CellType cell_entity_type(CellType type, int d, int index)
+{
+  if (int dim = mesh::cell_dim(type); d == dim)
+    return type;
+  else if (d == 1)
+    return CellType::interval;
+  else if (d == (dim - 1))
+    return mesh::cell_facet_type(type, index);
+  else
+    return CellType::point;
+}
 
 /// Return list of entities, where entities(e, k) is the local vertex
 /// index for the kth vertex of entity e of dimension dim
 graph::AdjacencyList<int> get_entity_vertices(CellType type, int dim);
 
 /// Get entities of dimension dim1 and that make up entities of dimension
-/// dim0
+/// dim0.
 graph::AdjacencyList<int> get_sub_entities(CellType type, int dim0, int dim1);
 
-/// Return topological dimension of cell type
-int cell_dim(CellType type);
-
-/// Number of entities of dimension dim
-/// @param[in] dim Entity dimension
-/// @param[in] type Cell type
-/// @return Number of entities in cell
+/// @brief Number of entities of dimension.
+///
+/// @param[in] dim Entity dimension.
+/// @param[in] type Cell type.
+/// @return Number of entities in cell.
 int cell_num_entities(CellType type, int dim);
 
-/// Check if cell is a simplex
-/// @param[in] type Cell type
-/// @return True is the cell type is a simplex
+/// @brief Check if cell is a simplex.
+///
+/// @param[in] type Cell type.
+/// @return True is the cell type is a simplex.
 bool is_simplex(CellType type);
 
-/// Number vertices for a cell type
+/// @brief Number vertices for a cell type.
+///
 /// @param[in] type Cell type
-/// @return The number of cell vertices
+/// @return Number of cell vertices
 int num_cell_vertices(CellType type);
 
 // [dim, entity] -> closure{sub_dim, (sub_entities)}