Unify refine interface

cpp/dolfinx/refinement/CMakeLists.txt

@@ -4,6 +4,7 @@ set(HEADERS_refinement
     ${CMAKE_CURRENT_SOURCE_DIR}/plaza.h
     ${CMAKE_CURRENT_SOURCE_DIR}/refine.h
     ${CMAKE_CURRENT_SOURCE_DIR}/utils.h
+    ${CMAKE_CURRENT_SOURCE_DIR}/option.h
     PARENT_SCOPE
 )
 

cpp/dolfinx/refinement/interval.h

@@ -19,23 +19,33 @@
 #include <stdexcept>
 #include <vector>
 
-namespace dolfinx::refinement
-{
-namespace impl
+#include "dolfinx/refinement/option.h"
+#include "dolfinx/refinement/utils.h"
+
+namespace dolfinx::refinement::interval
 {
 /// @brief Refine with markers returning new mesh data.
 ///
-/// @param[in] mesh Input mesh to be refined.
-/// @param[in] cells Indices of the cells that are marked for refinement.
-///
+/// @param[in] mesh Input mesh to be refined
+/// @param[in] cells Indices of the cells that are marked for refinement
+/// @param[in] option Refinement option indicating if parent cells
+/// and/or facets are to be computed.
 /// @return New mesh data: cell topology, vertex coordinates and parent
 /// cell indices.
 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>>
-compute_interval_refinement(const mesh::Mesh<T>& mesh,
-                            std::optional<std::span<const std::int32_t>> cells)
+           std::array<std::size_t, 2>, std::optional<std::vector<std::int32_t>>,
+           std::optional<std::vector<std::int8_t>>>
+compute_refinement_data(const mesh::Mesh<T>& mesh,
+                        std::optional<std::span<const std::int32_t>> cells,
+                        Option option)
 {
+  bool compute_parent_facet = option_parent_facet(option);
+  bool compute_parent_cell = option_parent_cell(option);
+
+  if (compute_parent_facet)
+    throw std::runtime_error("Parent facet computation not yet supported!");
+
   auto topology = mesh.topology();
   assert(topology);
   assert(topology->dim() == 1);
@@ -124,8 +134,14 @@ compute_interval_refinement(const mesh::Mesh<T>& mesh,
 
   std::vector<std::int64_t> cell_topology;
   cell_topology.reserve(refined_cell_count * 2);
-  std::vector<std::int32_t> parent_cell;
-  parent_cell.reserve(refined_cell_count);
+
+  std::optional<std::vector<std::int32_t>> parent_cell(std::nullopt);
+  if (compute_parent_cell)
+  {
+    parent_cell.emplace();
+    parent_cell->reserve(refined_cell_count);
+  }
+
   for (std::int32_t cell = 0; cell < map_c->size_local(); ++cell)
   {
     const auto& vertices = c_to_v->links(cell);
@@ -145,71 +161,30 @@ compute_interval_refinement(const mesh::Mesh<T>& mesh,
 
       // Add new cells/edges to refined topology
       cell_topology.insert(cell_topology.end(), {a, c, c, b});
-      parent_cell.insert(parent_cell.end(), {cell, cell});
+
+      if (compute_parent_cell)
+        parent_cell->insert(parent_cell->end(), {cell, cell});
     }
     else
     {
       // Copy the previous cell
       cell_topology.insert(cell_topology.end(), {a, b});
-      parent_cell.push_back(cell);
+
+      if (compute_parent_cell)
+        parent_cell->push_back(cell);
     }
   }
 
   assert(cell_topology.size() == 2 * refined_cell_count);
-  assert(parent_cell.size() == refined_cell_count);
+  assert(parent_cell->size() == (compute_parent_cell ? refined_cell_count : 0));
 
   std::vector<std::int32_t> offsets(refined_cell_count + 1);
   std::ranges::generate(offsets, [i = 0]() mutable { return 2 * i++; });
 
-  return {graph::AdjacencyList(std::move(cell_topology), std::move(offsets)),
-          std::move(new_vertex_coords), xshape, std::move(parent_cell)};
-}
-
-} // namespace impl
+  graph::AdjacencyList cell_adj(std::move(cell_topology), std::move(offsets));
 
-/// @brief Refine a (topologically) one-dimensional mesh by splitting
-/// cells, i.e. edges.
-///
-/// @param[in] mesh Mesh to be refined.
-/// @param[in] cells Optional indices of the cells that should be
-/// refined by this refinement. If not provided, all cells are
-/// considered marked for refinement, i.e. a uniform refinement is
-/// performed.
-/// @param[in] redistribute Option to enable redistribution of the
-/// refined mesh across processes.
-/// @param[in] ghost_mode Ghost mode of the refined mesh, default is
-/// ghost mode none.
-///
-/// @return Refined mesh, and list of parent cells and an array mapping
-/// the child cell index of the refined mesh to its parent cell index in
-/// the unrefined mesh.
-template <std::floating_point T>
-std::tuple<mesh::Mesh<T>, std::vector<std::int32_t>>
-refine_interval(const mesh::Mesh<T>& mesh,
-                std::optional<std::span<const std::int32_t>> cells,
-                bool redistribute,
-                mesh::GhostMode ghost_mode = mesh::GhostMode::shared_facet)
-{
-  if (mesh.topology()->cell_type() != mesh::CellType::interval)
-    throw std::runtime_error("Cell type not supported");
-  assert(mesh.topology()->dim() == 1);
-  assert(mesh.topology()->index_map(1));
-
-  auto [cell_adj, new_coords, xshape, parent_cell]
-      = impl::compute_interval_refinement(mesh, cells);
-  if (dolfinx::MPI::size(mesh.comm()) == 1)
-  {
-    return {mesh::create_mesh(mesh.comm(), cell_adj.array(),
-                              mesh.geometry().cmap(), new_coords, xshape,
-                              mesh::GhostMode::none),
-            std::move(parent_cell)};
-  }
-  else
-  {
-    return {partition<T>(mesh, cell_adj, std::span(new_coords), xshape,
-                         redistribute, ghost_mode),
-            std::move(parent_cell)};
-  }
+  return {std::move(cell_adj), std::move(new_vertex_coords), xshape,
+          std::move(parent_cell), std::nullopt};
 }
 
-} // namespace dolfinx::refinement
+} // namespace dolfinx::refinement::interval

cpp/dolfinx/refinement/option.h

@@ -0,0 +1,52 @@
+// Copyright (C) 2024 Paul Kühner
+//
+// This file is part of DOLFINX (https://www.fenicsproject.org)
+//
+// SPDX-License-Identifier:    LGPL-3.0-or-later
+
+#pragma once
+
+#include <cstdint>
+#include <type_traits>
+
+namespace dolfinx::refinement
+{
+/// @brief Options for data to compute during mesh refinement.
+enum class Option : std::uint8_t
+{
+  none = 0b00,         /*!< No extra data */
+  parent_facet = 0b01, /*!< Compute list of the cell-local facet indices in the
+          parent cell of each facet in each new cell (or -1 if no match)  */
+  parent_cell
+  = 0b10, /*!< Compute list with the parent cell index for each new cell */
+  parent_cell_and_facet = 0b11 /*< Both cell and facet parent data */
+};
+
+/// @brief Combine two refinement options into one, both flags will be
+/// set for the resulting option.
+inline constexpr Option operator|(Option a, Option b)
+{
+  using bitmask_t = std::underlying_type_t<Option>;
+  bitmask_t a_native = static_cast<bitmask_t>(a);
+  bitmask_t b_native = static_cast<bitmask_t>(b);
+  return static_cast<Option>(a_native | b_native);
+}
+
+/// @brief Check if parent_facet flag is set
+inline constexpr bool option_parent_facet(Option a)
+{
+  using bitmask_t = std::underlying_type_t<Option>;
+  bitmask_t a_native = static_cast<bitmask_t>(a);
+  bitmask_t facet_native = static_cast<bitmask_t>(Option::parent_facet);
+  return (a_native & facet_native) == facet_native;
+}
+
+/// @brief Check if parent_cell flag is set
+inline constexpr bool option_parent_cell(Option a)
+{
+  using bitmask_t = std::underlying_type_t<Option>;
+  bitmask_t a_native = static_cast<bitmask_t>(a);
+  bitmask_t facet_native = static_cast<bitmask_t>(Option::parent_cell);
+  return (a_native & facet_native) == facet_native;
+}
+} // namespace dolfinx::refinement