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

@@ -7,39 +7,44 @@
 #pragma once
 
 #include <algorithm>
-#include <cstddef>
-#include <mpi.h>
-
 #include <concepts>
+#include <cstddef>
 #include <cstdint>
 #include <optional>
 #include <stdexcept>
 #include <vector>
 
-#include "dolfinx/mesh/Mesh.h"
-#include "dolfinx/mesh/cell_types.h"
-#include "dolfinx/mesh/utils.h"
-#include "dolfinx/refinement/plaza.h"
+#include <mpi.h>
 
-namespace dolfinx::refinement
-{
+#include "dolfinx/mesh/Mesh.h"
+#include "dolfinx/refinement/option.h"
+#include "dolfinx/refinement/utils.h"
 
-namespace impl
+namespace dolfinx::refinement::interval
 {
 
 /// 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] 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);
@@ -127,8 +132,12 @@ 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)
   {
@@ -150,75 +159,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() == refined_cell_count * 2);
-  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++; });
 
   graph::AdjacencyList cell_adj(std::move(cell_topology), std::move(offsets));
 
   return {std::move(cell_adj), std::move(new_vertex_coords), xshape,
-          std::move(parent_cell)};
-}
-
-} // namespace impl
-
-/// Refines 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)};
-  }
+          std::move(parent_cell), std::nullopt};
 }
 
-} // namespace dolfinx::refinement
+} // namespace dolfinx::refinement::interval

cpp/dolfinx/refinement/option.h

@@ -0,0 +1,50 @@
+// 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,
+    parent_facet = 0b01,
+    parent_cell = 0b10,
+    parent_cell_and_facet = 0b11
+};
+
+/// @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;
+}
+}