DOLFINx 0.9.0
DOLFINx C++ interface
Loading...
Searching...
No Matches
interval.h
1// Copyright (C) 2024 Paul T. Kühner
2//
3// This file is part of DOLFINX (https://www.fenicsproject.org)
4//
5// SPDX-License-Identifier: LGPL-3.0-or-later
6
7#pragma once
8
9#include "dolfinx/mesh/Mesh.h"
10#include "dolfinx/mesh/cell_types.h"
11#include "dolfinx/mesh/utils.h"
12#include "dolfinx/refinement/plaza.h"
13#include <algorithm>
14#include <concepts>
15#include <cstddef>
16#include <cstdint>
17#include <mpi.h>
18#include <optional>
19#include <stdexcept>
20#include <vector>
21
22#include "dolfinx/refinement/option.h"
23#include "dolfinx/refinement/utils.h"
24
25namespace dolfinx::refinement::interval
26{
35template <std::floating_point T>
36std::tuple<graph::AdjacencyList<std::int64_t>, std::vector<T>,
37 std::array<std::size_t, 2>, std::optional<std::vector<std::int32_t>>,
38 std::optional<std::vector<std::int8_t>>>
39compute_refinement_data(const mesh::Mesh<T>& mesh,
40 std::optional<std::span<const std::int32_t>> cells,
41 Option option)
42{
43 bool compute_parent_facet = option_parent_facet(option);
44 bool compute_parent_cell = option_parent_cell(option);
45
46 if (compute_parent_facet)
47 throw std::runtime_error("Parent facet computation not yet supported!");
48
49 auto topology = mesh.topology();
50 assert(topology);
51 assert(topology->dim() == 1);
52 auto map_c = topology->index_map(1);
53 assert(map_c);
54
55 // TODO: creation of sharing ranks in external function? Also same
56 // code in use for plaza
57 // Get sharing ranks for each cell
58 graph::AdjacencyList<int> cell_ranks = map_c->index_to_dest_ranks();
59
60 // Create unique list of ranks that share cells (owners of ghosts plus
61 // ranks that ghost owned indices)
62 std::vector<int> ranks = cell_ranks.array();
63 std::ranges::sort(ranks);
64 auto to_remove = std::ranges::unique(ranks);
65 ranks.erase(to_remove.begin(), to_remove.end());
66
67 // Convert cell_ranks from global rank to to neighbourhood ranks
68 std::ranges::transform(cell_ranks.array(), cell_ranks.array().begin(),
69 [&ranks](auto r)
70 {
71 auto it = std::lower_bound(ranks.begin(),
72 ranks.end(), r);
73 assert(it != ranks.end() and *it == r);
74 return std::distance(ranks.begin(), it);
75 });
76
77 // Create refinement flag for cells
78 std::vector<std::int8_t> refinement_marker(
79 map_c->size_local() + map_c->num_ghosts(), !cells.has_value());
80
81 // Mark cells for refinement
82 std::vector<std::vector<std::int32_t>> marked_for_update(ranks.size());
83 if (cells.has_value())
84 {
85 std::ranges::for_each(cells.value(),
86 [&](auto cell)
87 {
88 if (!refinement_marker[cell])
89 {
90 refinement_marker[cell] = true;
91 for (int rank : cell_ranks.links(cell))
92 marked_for_update[rank].push_back(cell);
93 }
94 });
95 }
96
97 // Create neighborhood communicator for vertex creation
98 MPI_Comm neighbor_comm;
99 MPI_Dist_graph_create_adjacent(
100 mesh.comm(), ranks.size(), ranks.data(), MPI_UNWEIGHTED, ranks.size(),
101 ranks.data(), MPI_UNWEIGHTED, MPI_INFO_NULL, false, &neighbor_comm);
102
103 // Communicate ghost cells that might have been marked. This is not
104 // necessary for a uniform refinement.
105 if (cells.has_value())
106 {
107 update_logical_edgefunction(neighbor_comm, marked_for_update,
108 refinement_marker, *map_c);
109 }
110
111 // Construct the new vertices
112 const auto [new_vertex_map, new_vertex_coords, xshape]
113 = create_new_vertices(neighbor_comm, cell_ranks, mesh, refinement_marker);
114 MPI_Comm_free(&neighbor_comm);
115
116 auto c_to_v = mesh.topology()->connectivity(1, 0);
117 assert(c_to_v);
118
119 // Get the count of cells to refine, note: we only consider non-ghost
120 // cells
121 const std::int32_t number_of_refined_cells
122 = std::count(refinement_marker.begin(),
123 std::next(refinement_marker.begin(),
124 mesh.topology()->index_map(1)->size_local()),
125 true);
126
127 // Produce local global indices, by padding out the previous index map
128 std::vector<std::int64_t> global_indices
129 = adjust_indices(*mesh.topology()->index_map(0), number_of_refined_cells);
130
131 // Build the topology on the new vertices
132 const std::int32_t refined_cell_count
133 = mesh.topology()->index_map(1)->size_local() + number_of_refined_cells;
134
135 std::vector<std::int64_t> cell_topology;
136 cell_topology.reserve(refined_cell_count * 2);
137
138 std::optional<std::vector<std::int32_t>> parent_cell(std::nullopt);
139 if (compute_parent_cell)
140 {
141 parent_cell.emplace();
142 parent_cell->reserve(refined_cell_count);
143 }
144
145 for (std::int32_t cell = 0; cell < map_c->size_local(); ++cell)
146 {
147 const auto& vertices = c_to_v->links(cell);
148 assert(vertices.size() == 2);
149
150 // We consider a cell (defined by global vertices)
151 // a ----------- b
152 const std::int64_t a = global_indices[vertices[0]];
153 const std::int64_t b = global_indices[vertices[1]];
154 if (refinement_marker[cell])
155 {
156 // Find (global) index of new midpoint vertex:
157 // a --- c --- b
158 auto it = new_vertex_map.find(cell);
159 assert(it != new_vertex_map.end());
160 const std::int64_t c = it->second;
161
162 // Add new cells/edges to refined topology
163 cell_topology.insert(cell_topology.end(), {a, c, c, b});
164
165 if (compute_parent_cell)
166 parent_cell->insert(parent_cell->end(), {cell, cell});
167 }
168 else
169 {
170 // Copy the previous cell
171 cell_topology.insert(cell_topology.end(), {a, b});
172
173 if (compute_parent_cell)
174 parent_cell->push_back(cell);
175 }
176 }
177
178 assert(cell_topology.size() == 2 * refined_cell_count);
179 assert(parent_cell->size() == (compute_parent_cell ? refined_cell_count : 0));
180
181 std::vector<std::int32_t> offsets(refined_cell_count + 1);
182 std::ranges::generate(offsets, [i = 0]() mutable { return 2 * i++; });
183
184 graph::AdjacencyList cell_adj(std::move(cell_topology), std::move(offsets));
185
186 return {std::move(cell_adj), std::move(new_vertex_coords), xshape,
187 std::move(parent_cell), std::nullopt};
188}
189
190} // namespace dolfinx::refinement::interval
Functions supporting mesh operations.
void cells(la::SparsityPattern &pattern, std::array< std::span< const std::int32_t >, 2 > cells, std::array< std::reference_wrapper< const DofMap >, 2 > dofmaps)
Iterate over cells and insert entries into sparsity pattern.
Definition sparsitybuild.cpp:16
constexpr bool option_parent_cell(Option a)
Check if parent_cell flag is set.
Definition option.h:45
Option
Options for data to compute during mesh refinement.
Definition option.h:16
constexpr bool option_parent_facet(Option a)
Check if parent_facet flag is set.
Definition option.h:36
std::tuple< std::map< std::int32_t, std::int64_t >, std::vector< T >, std::array< std::size_t, 2 > > create_new_vertices(MPI_Comm comm, const graph::AdjacencyList< int > &shared_edges, const mesh::Mesh< T > &mesh, std::span< const std::int8_t > marked_edges)
Add new vertex for each marked edge, and create new_vertex_coordinates and global_edge->new_vertex ma...
Definition utils.h:150
std::vector< std::int64_t > adjust_indices(const common::IndexMap &map, std::int32_t n)
Given an index map, add "n" extra indices at the end of local range.
Definition utils.cpp:103
void update_logical_edgefunction(MPI_Comm comm, const std::vector< std::vector< std::int32_t > > &marked_for_update, std::span< std::int8_t > marked_edges, const common::IndexMap &map)
Communicate edge markers between processes that share edges.
Definition utils.cpp:47