Improve mesh refinement documentation (#3010)

* Document refinement functions

* Explain parent_facet_index

* Fix Garth's comments

* More comments

* Revert spaces

* Apply suggestions from code review

* Fix docs

* Apply suggestions from code review

Co-authored-by: Chris Richardson <chris@bpi.cam.ac.uk>

* Tweak formatting.

* Typo.

* Typos

---------

Co-authored-by: Chris Richardson <chris@bpi.cam.ac.uk>
Co-authored-by: Jack S. Hale <mail@jackhale.co.uk>

cpp/dolfinx/refinement/plaza.cpp

@@ -24,13 +24,23 @@ using namespace dolfinx::refinement;
 namespace
 {
 //-----------------------------------------------------------------------------
-// 2D version of subdivision allowing for uniform subdivision (flag)
+/// 2D version of subdivision allowing for uniform subdivision (flag)
+/// @param[in] indices Vector containing the
+/// global indices for the original vertices and potential new vertices at each
+/// edge. If an edge is not refined its corresponding entry is -1. Size
+/// `num_vertices + num_edges`
+/// @param[in] longest_edge Local index of the longest edge in the triangle.
+/// @param[in] uniform If true, the triangle is subdivided into four similar
+/// sub-triangles.
+/// @returns Local indices for each sub-divived triangle
 std::vector<std::int32_t> get_triangles(std::span<const std::int64_t> indices,
                                         const std::int32_t longest_edge,
                                         bool uniform)
 {
-  // v0 and v1 are at ends of longest_edge (e2) opposite vertex has same
-  // index as longest_edge
+  // NOTE: The assumption below is based on the UFC ordering of a triangle, i.e.
+  // that the N-th edge of a triangle is the edge where the N-th vertex is not
+  // part of the set. v0 and v1 are at ends of longest_edge (e2) opposite vertex
+  // has same index as longest_edge
   const std::int32_t v0 = (longest_edge + 1) % 3;
   const std::int32_t v1 = (longest_edge + 2) % 3;
   const std::int32_t v2 = longest_edge;

cpp/dolfinx/refinement/plaza.h

@@ -116,13 +116,16 @@ compute_parent_facets(std::span<const std::int32_t> simplex_set)
 /// (cell local indexing). A flag indicates if a uniform subdivision is
 /// preferable in 2D.
 ///
-/// @param[in] indices Vector indicating which edges are to be
-///   split (value >=0)
+/// @param[in] indices Vector containing the global indices for the original
+/// vertices and potential new vertices at each edge. Size (num_vertices +
+/// num_edges). If an edge is not refined its corresponding entry is -1
 /// @param[in] longest_edge Vector indicating the longest edge for each
-///   triangle. For tdim=2, one entry, for tdim=3, four entries.
+/// triangle in the cell. For triangular cells (2D) there is only one value,
+/// and for tetrahedra (3D) there are four values, one for each facet. The
+/// values give the local edge indices of the cell.
 /// @param[in] tdim Topological dimension (2 or 3)
-/// @param[in] uniform Make a "uniform" subdivision with all triangles
-///   being similar shape
+/// @param[in] uniform Make a "uniform" subdivision with all triangles being
+/// similar shape
 /// @return
 std::vector<std::int32_t>
 get_simplices(std::span<const std::int64_t> indices,
@@ -136,7 +139,15 @@ void enforce_rules(MPI_Comm comm, const graph::AdjacencyList<int>& shared_edges,
                    const mesh::Topology& topology,
                    std::span<const std::int32_t> long_edge);
 
-/// Get the longest edge of each face (using local mesh index)
+/// @brief Get the longest edge of each face (using local mesh index)
+///
+/// @note Edge ratio ok returns an array in 2D, where it checks if the ratio
+/// between the shortest and longest edge of a cell is bigger than sqrt(2)/2. In
+/// 3D an empty array is returned
+///
+/// @param[in] mesh The mesh
+/// @return A tuple (longest edge, edge ratio ok) where longest edge gives the
+/// local index of the longest edge for each face.
 template <std::floating_point T>
 std::pair<std::vector<std::int32_t>, std::vector<std::int8_t>>
 face_long_edge(const mesh::Mesh<T>& mesh)
@@ -253,7 +264,30 @@ face_long_edge(const mesh::Mesh<T>& mesh)
   return std::pair(std::move(long_edge), std::move(edge_ratio_ok));
 }
 
-/// Convenient interface for both uniform and marker refinement
+/// @brief Convenient interface for both uniform and marker refinement
+///
+/// @note The parent facet map gives you the map from a cell given by parent
+/// cell map to the local index (relative to the cell), e.g. the i-th entry of
+/// parent facets relates to the local facet index of the i-th entry parent
+/// cell.
+///
+/// @param[in] neighbor_comm Neighbourhood communciator scattering owned edges
+/// to processes with ghosts
+/// @param[in] marked_edges A marker for all edges on the process (local +
+/// ghosts) indicating if an edge should be refined
+/// @param[in] shared_edges For each local edge on a process map to ghost
+/// processes
+/// @param[in] mesh The mesh
+/// @param[in] long_edge A map from each face to its longest edge. Index is
+/// local to the process.
+/// @param[in] edge_ratio_ok For each face in a 2D mesh this error contains a
+/// marker indicating if the ratio between smallest and largest edge is bigger
+/// than sqrt(2)/2
+/// @param[in] option Option to compute additional information relating refined
+/// and original mesh entities
+/// @return (0) The new mesh topology, (1) the new flattened mesh geometry, (3)
+/// Shape of the new geometry_shape, (4) Map from new cells to parent cells
+/// and (5) map from refined facets to parent facets.
 template <std::floating_point T>
 std::tuple<graph::AdjacencyList<std::int64_t>, std::vector<T>,
            std::array<std::size_t, 2>, std::vector<std::int32_t>,

cpp/dolfinx/refinement/utils.cpp

@@ -103,9 +103,6 @@ refinement::adjust_indices(const common::IndexMap& map, std::int32_t n)
 {
   // NOTE: Is this effectively concatenating index maps?
 
-  // Add in an extra "n" indices at the end of the current local_range
-  // of "index_map", and adjust existing indices to match.
-
   // Get offset for 'n' for this process
   const std::int64_t num_local = n;
   std::int64_t global_offset = 0;

cpp/dolfinx/refinement/utils.h

@@ -46,9 +46,12 @@ std::int64_t local_to_global(std::int32_t local_index,
 /// Create geometric points of new Mesh, from current Mesh and a
 /// edge_to_vertex map listing the new local points (midpoints of those
 /// edges)
-/// @param mesh
-/// @param local_edge_to_new_vertex
-/// @return array of points
+///
+/// @param mesh Current mesh
+/// @param local_edge_to_new_vertex A map from a local edge to the new
+/// global vertex index that will be inserted at its midpoint
+/// @return (1) Array of new (flattened) mesh geometry and (2) its
+/// multi-dimensional shape
 template <std::floating_point T>
 std::pair<std::vector<T>, std::array<std::size_t, 2>> create_new_geometry(
     const mesh::Mesh<T>& mesh,
@@ -135,12 +138,14 @@ void update_logical_edgefunction(
 /// Communicate new vertices with MPI to all affected processes.
 ///
 /// @param[in] comm MPI Communicator for neighborhood
-/// @param[in] shared_edges
+/// @param[in] shared_edges For each local edge on a process map to ghost
+/// processes
 /// @param[in] mesh Existing mesh
-/// @param[in] marked_edges
-/// @return (0) map from local edge index to new vertex global index,
-/// (1) the coordinates of the new vertices (row-major storage) and (2)
-/// the shape of the new coordinates.
+/// @param[in] marked_edges A marker for all edges on the process (local +
+/// ghosts) indicating if an edge should be refined
+/// @return (0) map from local edge index to new vertex global index, (1) the
+/// coordinates of the new vertices (row-major storage) and (2) the shape of the
+/// new coordinates.
 template <std::floating_point T>
 std::tuple<std::map<std::int32_t, std::int64_t>, std::vector<T>,
            std::array<std::size_t, 2>>
@@ -153,7 +158,8 @@ create_new_vertices(MPI_Comm comm,
   std::shared_ptr<const common::IndexMap> edge_index_map
       = mesh.topology()->index_map(1);
 
-  // Add new edge midpoints to list of vertices
+  // Add new edge midpoints to list of vertices. The new vertex will be owned by
+  // the process owning the edge.
   int n = 0;
   std::map<std::int32_t, std::int64_t> local_edge_to_new_vertex;
   for (int local_i = 0; local_i < edge_index_map->size_local(); ++local_i)
@@ -303,18 +309,14 @@ mesh::Mesh<T> partition(const mesh::Mesh<T>& old_mesh,
   }
 }
 
-/// @todo Fix docstring. It is unclear.
-///
-/// @brief Add indices to account for extra n values on this process.
+/// @brief Given an index map, add "n" extra indices at the end of local range
 ///
-/// This is a utility to help add new topological vertices on each
-/// process into the space of the index map.
-///
-/// @param[in] map Index map for the current mesh vertices
-/// @param[in] n Number of new entries to be accommodated on this
-/// process
-/// @return Global indices as if "n" extra values are appended on each
-/// process
+/// @note The returned global indices (local and ghosts) are adjust for the
+/// addition of new local indices.
+/// @param[in] map Index map
+/// @param[in] n Number of new local indices
+/// @return New global indices for both owned and ghosted indices in input
+/// index map.
 std::vector<std::int64_t> adjust_indices(const common::IndexMap& map,
                                          std::int32_t n);