37auto compute_parent_facets(std::span<const std::int32_t> simplex_set)
39 static_assert(tdim == 2 or tdim == 3);
40 assert(simplex_set.size() % (tdim + 1) == 0);
42 = std::conditional_t<tdim == 2, std::array<std::int8_t, 12>,
43 std::array<std::int8_t, 32>>;
51 constexpr std::array<std::array<int, 3>, 3> facet_table_2d{
52 {{1, 2, 3}, {0, 2, 4}, {0, 1, 5}}};
54 constexpr std::array<std::array<int, 6>, 4> facet_table_3d{
60 const int ncells = simplex_set.size() / (tdim + 1);
61 for (
int fpi = 0; fpi < (tdim + 1); ++fpi)
64 for (
int cc = 0; cc < ncells; ++cc)
66 for (
int fci = 0; fci < (tdim + 1); ++fci)
69 std::array<int, tdim> cf, set_output;
71 int num_common_vertices;
72 if constexpr (tdim == 2)
74 for (
int j = 0; j < tdim; ++j)
75 cf[j] = simplex_set[cc * 3 + facet_table_2d[fci][j]];
77 std::ranges::sort(cf);
78 auto [last1, last2, it_last] = std::ranges::set_intersection(
79 facet_table_2d[fpi], cf, set_output.begin());
80 num_common_vertices = std::distance(set_output.begin(), it_last);
84 for (
int j = 0; j < tdim; ++j)
85 cf[j] = simplex_set[cc * 4 + facet_table_3d[fci][j]];
87 std::ranges::sort(cf);
88 auto [last1, last2, it_last] = std::ranges::set_intersection(
89 facet_table_3d[fpi], cf, set_output.begin());
90 num_common_vertices = std::distance(set_output.begin(), it_last);
93 if (num_common_vertices == tdim)
124std::pair<std::array<std::int32_t, 32>, std::size_t>
125get_simplices(std::span<const std::int64_t> indices,
126 std::span<const std::int32_t> longest_edge,
int tdim,
132 std::span<std::int8_t> marked_edges,
134 std::span<const std::int32_t> long_edge);
145template <std::
floating_po
int T>
146std::pair<std::vector<std::int32_t>, std::vector<std::int8_t>>
149 const int tdim =
mesh.topology()->dim();
151 mesh.topology_mutable()->create_entities(1);
152 mesh.topology_mutable()->create_entities(2);
153 mesh.topology_mutable()->create_connectivity(2, 1);
154 mesh.topology_mutable()->create_connectivity(1, tdim);
155 mesh.topology_mutable()->create_connectivity(tdim, 2);
157 std::int64_t num_faces =
mesh.topology()->index_map(2)->size_local()
158 +
mesh.topology()->index_map(2)->num_ghosts();
161 std::vector<std::int32_t> long_edge(num_faces);
162 std::vector<std::int8_t> edge_ratio_ok;
166 const T min_ratio = std::sqrt(2.0) / 2.0;
168 edge_ratio_ok.resize(num_faces);
170 auto x_dofmap =
mesh.geometry().dofmap();
172 auto c_to_v =
mesh.topology()->connectivity(tdim, 0);
174 auto e_to_c =
mesh.topology()->connectivity(1, tdim);
176 auto e_to_v =
mesh.topology()->connectivity(1, 0);
180 auto map_e =
mesh.topology()->index_map(1);
182 std::vector<T> edge_length(map_e->size_local() + map_e->num_ghosts());
183 for (std::size_t e = 0; e < edge_length.size(); ++e)
186 auto cells = e_to_c->links(e);
187 assert(!cells.empty());
188 auto cell_vertices = c_to_v->links(cells.front());
189 auto edge_vertices = e_to_v->links(e);
192 auto it0 = std::find(cell_vertices.begin(), cell_vertices.end(),
194 assert(it0 != cell_vertices.end());
195 const std::size_t local0 = std::distance(cell_vertices.begin(), it0);
196 auto it1 = std::find(cell_vertices.begin(), cell_vertices.end(),
198 assert(it1 != cell_vertices.end());
199 const std::size_t local1 = std::distance(cell_vertices.begin(), it1);
201 auto x_dofs = md::submdspan(x_dofmap, cells.front(), md::full_extent);
202 std::span<const T, 3> x0(
mesh.geometry().x().data() + 3 * x_dofs[local0],
204 std::span<const T, 3> x1(
mesh.geometry().x().data() + 3 * x_dofs[local1],
208 edge_length[e] = std::sqrt(std::transform_reduce(
209 x0.begin(), x0.end(), x1.begin(), 0.0, std::plus<>(),
210 [](
auto x0,
auto x1) { return (x0 - x1) * (x0 - x1); }));
214 auto f_to_v =
mesh.topology()->connectivity(2, 0);
216 auto f_to_e =
mesh.topology()->connectivity(2, 1);
218 const std::vector global_indices
219 =
mesh.topology()->index_map(0)->global_indices();
220 for (
int f = 0; f < f_to_v->num_nodes(); ++f)
222 auto face_edges = f_to_e->links(f);
224 std::int32_t imax = 0;
226 T min_len = std::numeric_limits<T>::max();
228 for (
int i = 0; i < 3; ++i)
230 const T e_len = edge_length[face_edges[i]];
231 min_len = std::min(e_len, min_len);
237 else if (tdim == 3 and e_len == max_len)
242 auto vertices = f_to_v->links(f);
243 const int vmax = vertices[imax];
244 const int vi = vertices[i];
245 if (global_indices[vi] > global_indices[vmax])
252 edge_ratio_ok[f] = (min_len / max_len >= min_ratio);
254 long_edge[f] = face_edges[imax];
257 return std::pair(std::move(long_edge), std::move(edge_ratio_ok));
284template <std::
floating_po
int T>
285std::tuple<graph::AdjacencyList<std::int64_t>, std::vector<T>,
286 std::array<std::size_t, 2>, std::optional<std::vector<std::int32_t>>,
287 std::optional<std::vector<std::int8_t>>>
288compute_refinement(MPI_Comm neighbor_comm,
289 std::span<const std::int8_t> marked_edges,
292 std::span<const std::int32_t> long_edge,
293 std::span<const std::int8_t> edge_ratio_ok,
Option option)
295 int tdim =
mesh.topology()->dim();
296 int num_cell_edges = tdim * 3 - 3;
297 int num_cell_vertices = tdim + 1;
303 const auto [new_vertex_map, new_vertex_coords, xshape]
306 std::optional<std::vector<std::int32_t>>
parent_cell(std::nullopt);
307 if (compute_parent_cell)
310 std::optional<std::vector<std::int8_t>>
parent_facet(std::nullopt);
314 std::vector<std::int64_t> indices(num_cell_vertices + num_cell_edges);
315 std::vector<std::int32_t> simplex_set;
317 auto map_c =
mesh.topology()->index_map(tdim);
319 auto c_to_v =
mesh.topology()->connectivity(tdim, 0);
321 auto c_to_e =
mesh.topology()->connectivity(tdim, 1);
323 auto c_to_f =
mesh.topology()->connectivity(tdim, 2);
326 std::int32_t num_new_vertices_local = std::count(
327 marked_edges.begin(),
328 marked_edges.begin() +
mesh.topology()->index_map(1)->size_local(),
true);
330 std::vector<std::int64_t> global_indices
333 const std::int32_t num_cells = map_c->size_local();
336 std::vector<std::int64_t> cell_topology;
337 for (
int c = 0; c < num_cells; ++c)
343 auto vertices = c_to_v->links(c);
344 for (std::size_t v = 0; v < vertices.size(); ++v)
345 indices[v] = global_indices[vertices[v]];
348 auto edges = c_to_e->links(c);
349 bool no_edge_marked =
true;
350 for (std::size_t ei = 0; ei < edges.size(); ++ei)
352 if (marked_edges[edges[ei]])
354 no_edge_marked =
false;
355 auto it = new_vertex_map.find(edges[ei]);
356 assert(it != new_vertex_map.end());
357 indices[num_cell_vertices + ei] = it->second;
360 indices[num_cell_vertices + ei] = -1;
366 for (
auto v : vertices)
367 cell_topology.push_back(global_indices[v]);
369 if (compute_parent_cell)
384 std::vector<std::int32_t> longest_edge;
385 for (
auto f : c_to_f->links(c))
386 longest_edge.push_back(long_edge[f]);
389 for (std::int32_t& p : longest_edge)
391 for (std::size_t ej = 0; ej < edges.size(); ++ej)
401 const bool uniform = (tdim == 2) ? edge_ratio_ok[c] :
false;
402 const auto [simplex_set_b, simplex_set_size]
403 = get_simplices(indices, longest_edge, tdim, uniform);
404 std::span<const std::int32_t> simplex_set(simplex_set_b.data(),
408 const std::int32_t ncells = simplex_set.size() / num_cell_vertices;
409 if (compute_parent_cell)
411 for (std::int32_t i = 0; i < ncells; ++i)
419 auto npf = compute_parent_facets<3>(simplex_set);
421 std::next(npf.begin(), simplex_set.size()));
425 auto npf = compute_parent_facets<2>(simplex_set);
427 std::next(npf.begin(), simplex_set.size()));
432 for (std::int32_t v : simplex_set)
433 cell_topology.push_back(indices[v]);
437 assert(cell_topology.size() % num_cell_vertices == 0);
438 std::vector<std::int32_t> offsets(
439 cell_topology.size() / num_cell_vertices + 1, 0);
440 for (std::size_t i = 0; i < offsets.size() - 1; ++i)
441 offsets[i + 1] = offsets[i] + num_cell_vertices;
444 return {std::move(cell_adj), std::move(new_vertex_coords), xshape,
458template <std::
floating_po
int T>
459std::tuple<graph::AdjacencyList<std::int64_t>, std::vector<T>,
460 std::array<std::size_t, 2>, std::optional<std::vector<std::int32_t>>,
461 std::optional<std::vector<std::int8_t>>>
463 std::optional<std::span<const std::int32_t>> edges,
467 auto topology =
mesh.topology();
470 if (topology->cell_type() != mesh::CellType::triangle
471 and topology->cell_type() != mesh::CellType::tetrahedron)
473 throw std::runtime_error(
"Cell type not supported");
476 auto map_e = topology->index_map(1);
478 throw std::runtime_error(
"Edges must be initialised");
485 std::vector<int> ranks(edge_ranks.
array().begin(), edge_ranks.
array().end());
486 std::ranges::sort(ranks);
487 auto [unique_end, range_end] = std::ranges::unique(ranks);
488 ranks.erase(unique_end, range_end);
491 std::ranges::transform(edge_ranks.
array(), edge_ranks.
array().begin(),
494 auto it = std::ranges::lower_bound(ranks, r);
495 assert(it != ranks.end() and *it == r);
496 return std::distance(ranks.begin(), it);
500 std::vector<std::int8_t> marked_edges(
501 map_e->size_local() + map_e->num_ghosts(), !edges.has_value());
502 std::vector<std::vector<std::int32_t>> marked_for_update(ranks.size());
506 for (
auto edge : edges.value())
508 if (!marked_edges[edge])
510 marked_edges[edge] =
true;
513 for (
int rank : edge_ranks.
links(edge))
514 marked_for_update[rank].push_back(edge);
520 MPI_Dist_graph_create_adjacent(
mesh.comm(), ranks.size(), ranks.data(),
521 MPI_UNWEIGHTED, ranks.size(), ranks.data(),
522 MPI_UNWEIGHTED, MPI_INFO_NULL,
false, &comm);
529 const auto [long_edge, edge_ratio_ok] = impl::face_long_edge(
mesh);
530 impl::enforce_rules(comm, edge_ranks, marked_edges, *topology, long_edge);
533 = impl::compute_refinement(comm, marked_edges, edge_ranks,
mesh,
534 long_edge, edge_ratio_ok, option);
535 MPI_Comm_free(&comm);
537 return {std::move(cell_adj), std::move(new_vertex_coords), xshape,