48std::vector<std::int8_t>
49compute_parent_facets(std::span<const std::int32_t> simplex_set)
52 assert(simplex_set.size() % (tdim + 1) == 0);
53 std::vector<std::int8_t>
parent_facet(simplex_set.size(), -1);
58 constexpr std::array<std::array<int, 3>, 3> facet_table_2d{
59 {{1, 2, 3}, {0, 2, 4}, {0, 1, 5}}};
61 constexpr std::array<std::array<int, 6>, 4> facet_table_3d{
67 const int ncells = simplex_set.size() / (tdim + 1);
68 for (
int fpi = 0; fpi < (tdim + 1); ++fpi)
71 for (
int cc = 0; cc < ncells; ++cc)
73 for (
int fci = 0; fci < (tdim + 1); ++fci)
76 std::array<int, tdim> cf, set_output;
78 int num_common_vertices;
79 if constexpr (tdim == 2)
81 for (
int j = 0; j < tdim; ++j)
82 cf[j] = simplex_set[cc * 3 + facet_table_2d[fci][j]];
83 std::sort(cf.begin(), cf.end());
84 auto it = std::set_intersection(
85 facet_table_2d[fpi].begin(), facet_table_2d[fpi].end(),
86 cf.begin(), cf.end(), set_output.begin());
87 num_common_vertices = std::distance(set_output.begin(), it);
91 for (
int j = 0; j < tdim; ++j)
92 cf[j] = simplex_set[cc * 4 + facet_table_3d[fci][j]];
93 std::sort(cf.begin(), cf.end());
94 auto it = std::set_intersection(
95 facet_table_3d[fpi].begin(), facet_table_3d[fpi].end(),
96 cf.begin(), cf.end(), set_output.begin());
97 num_common_vertices = std::distance(set_output.begin(), it);
100 if (num_common_vertices == tdim)
127std::vector<std::int32_t>
128get_simplices(std::span<const std::int64_t> indices,
129 std::span<const std::int32_t> longest_edge,
int tdim,
134void enforce_rules(MPI_Comm comm,
const graph::AdjacencyList<int>& shared_edges,
135 std::vector<std::int8_t>& marked_edges,
136 const mesh::Topology& topology,
137 std::span<const std::int32_t> long_edge);
140template <std::
floating_po
int T>
141std::pair<std::vector<std::int32_t>, std::vector<std::int8_t>>
142face_long_edge(
const mesh::Mesh<T>& mesh)
144 const int tdim = mesh.topology()->dim();
146 mesh.topology_mutable()->create_entities(1);
147 mesh.topology_mutable()->create_entities(2);
148 mesh.topology_mutable()->create_connectivity(2, 1);
149 mesh.topology_mutable()->create_connectivity(1, tdim);
150 mesh.topology_mutable()->create_connectivity(tdim, 2);
152 std::int64_t num_faces = mesh.topology()->index_map(2)->size_local()
153 + mesh.topology()->index_map(2)->num_ghosts();
156 std::vector<std::int32_t> long_edge(num_faces);
157 std::vector<std::int8_t> edge_ratio_ok;
161 const T min_ratio = sqrt(2.0) / 2.0;
163 edge_ratio_ok.resize(num_faces);
165 auto x_dofmap = mesh.geometry().dofmap();
167 auto c_to_v = mesh.topology()->connectivity(tdim, 0);
169 auto e_to_c = mesh.topology()->connectivity(1, tdim);
171 auto e_to_v = mesh.topology()->connectivity(1, 0);
175 auto map_e = mesh.topology()->index_map(1);
177 std::vector<T> edge_length(map_e->size_local() + map_e->num_ghosts());
178 for (std::size_t e = 0; e < edge_length.size(); ++e)
181 auto cells = e_to_c->links(e);
182 assert(!
cells.empty());
183 auto cell_vertices = c_to_v->links(
cells.front());
184 auto edge_vertices = e_to_v->links(e);
187 auto it0 = std::find(cell_vertices.begin(), cell_vertices.end(),
189 assert(it0 != cell_vertices.end());
190 const std::size_t local0 = std::distance(cell_vertices.begin(), it0);
191 auto it1 = std::find(cell_vertices.begin(), cell_vertices.end(),
193 assert(it1 != cell_vertices.end());
194 const std::size_t local1 = std::distance(cell_vertices.begin(), it1);
196 auto x_dofs = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan(
197 x_dofmap,
cells.front(), MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent);
198 std::span<const T, 3> x0(mesh.geometry().x().data() + 3 * x_dofs[local0],
200 std::span<const T, 3> x1(mesh.geometry().x().data() + 3 * x_dofs[local1],
204 edge_length[e] = std::sqrt(std::transform_reduce(
205 x0.begin(), x0.end(), x1.begin(), 0.0, std::plus<>(),
206 [](
auto x0,
auto x1) { return (x0 - x1) * (x0 - x1); }));
210 auto f_to_v = mesh.topology()->connectivity(2, 0);
212 auto f_to_e = mesh.topology()->connectivity(2, 1);
214 const std::vector global_indices
215 = mesh.topology()->index_map(0)->global_indices();
216 for (
int f = 0; f < f_to_v->num_nodes(); ++f)
218 auto face_edges = f_to_e->links(f);
220 std::int32_t imax = 0;
222 T min_len = std::numeric_limits<T>::max();
224 for (
int i = 0; i < 3; ++i)
226 const T e_len = edge_length[face_edges[i]];
227 min_len = std::min(e_len, min_len);
233 else if (tdim == 3 and e_len == max_len)
238 auto vertices = f_to_v->links(f);
239 const int vmax = vertices[imax];
240 const int vi = vertices[i];
241 if (global_indices[vi] > global_indices[vmax])
248 edge_ratio_ok[f] = (min_len / max_len >= min_ratio);
250 long_edge[f] = face_edges[imax];
253 return std::pair(std::move(long_edge), std::move(edge_ratio_ok));
257template <std::
floating_po
int T>
258std::tuple<graph::AdjacencyList<std::int64_t>, std::vector<T>,
259 std::array<std::size_t, 2>, std::vector<std::int32_t>,
260 std::vector<std::int8_t>>
261compute_refinement(MPI_Comm neighbor_comm,
262 const std::vector<std::int8_t>& marked_edges,
263 const graph::AdjacencyList<int>& shared_edges,
264 const mesh::Mesh<T>& mesh,
265 const std::vector<std::int32_t>& long_edge,
266 const std::vector<std::int8_t>& edge_ratio_ok,
269 int tdim = mesh.topology()->dim();
270 int num_cell_edges = tdim * 3 - 3;
274 bool compute_parent_cell
279 const auto [new_vertex_map, new_vertex_coords, xshape]
284 std::vector<std::int64_t> indices(num_cell_vertices + num_cell_edges);
285 std::vector<std::int32_t> simplex_set;
287 auto map_c = mesh.topology()->index_map(tdim);
289 auto c_to_v = mesh.topology()->connectivity(tdim, 0);
291 auto c_to_e = mesh.topology()->connectivity(tdim, 1);
293 auto c_to_f = mesh.topology()->connectivity(tdim, 2);
296 std::int32_t num_new_vertices_local = std::count(
297 marked_edges.begin(),
298 marked_edges.begin() + mesh.topology()->index_map(1)->size_local(),
true);
300 std::vector<std::int64_t> global_indices
301 =
adjust_indices(*mesh.topology()->index_map(0), num_new_vertices_local);
303 const int num_cells = map_c->size_local();
306 std::vector<std::int64_t> cell_topology;
307 for (
int c = 0; c < num_cells; ++c)
313 auto vertices = c_to_v->links(c);
314 for (std::size_t v = 0; v < vertices.size(); ++v)
315 indices[v] = global_indices[vertices[v]];
318 auto edges = c_to_e->links(c);
319 bool no_edge_marked =
true;
320 for (std::size_t ei = 0; ei < edges.size(); ++ei)
322 if (marked_edges[edges[ei]])
324 no_edge_marked =
false;
325 auto it = new_vertex_map.find(edges[ei]);
326 assert(it != new_vertex_map.end());
336 for (
auto v : vertices)
337 cell_topology.push_back(global_indices[v]);
339 if (compute_parent_cell)
354 std::vector<std::int32_t> longest_edge;
355 for (
auto f : c_to_f->links(c))
356 longest_edge.push_back(long_edge[f]);
359 for (std::int32_t& p : longest_edge)
361 for (std::size_t ej = 0; ej < edges.size(); ++ej)
371 const bool uniform = (tdim == 2) ? edge_ratio_ok[c] : false;
374 simplex_set = get_simplices(indices, longest_edge, tdim, uniform);
379 if (compute_parent_cell)
381 for (std::int32_t i = 0; i < ncells; ++i)
387 std::vector<std::int8_t> npf;
389 npf = compute_parent_facets<3>(simplex_set);
391 npf = compute_parent_facets<2>(simplex_set);
396 for (std::int32_t v : simplex_set)
397 cell_topology.push_back(indices[v]);
401 assert(cell_topology.size() % num_cell_vertices == 0);
402 std::vector<std::int32_t> offsets(
403 cell_topology.size() / num_cell_vertices + 1, 0);
404 for (std::size_t i = 0; i < offsets.size() - 1; ++i)
405 offsets[i + 1] = offsets[i] + num_cell_vertices;
406 graph::AdjacencyList cell_adj(std::move(cell_topology), std::move(offsets));
408 return {std::move(cell_adj), std::move(new_vertex_coords), xshape,
423template <std::
floating_po
int T>
424std::tuple<mesh::Mesh<T>, std::vector<std::int32_t>, std::vector<std::int8_t>>
433 mesh.geometry().cmap(), new_coords, xshape,
434 mesh::GhostMode::none),
439 std::shared_ptr<const common::IndexMap> map_c
440 = mesh.topology()->index_map(mesh.topology()->dim());
441 const int num_ghost_cells = map_c->num_ghosts();
444 int max_ghost_cells = 0;
445 MPI_Allreduce(&num_ghost_cells, &max_ghost_cells, 1, MPI_INT, MPI_MAX,
450 ? mesh::GhostMode::none
451 : mesh::GhostMode::shared_facet;
452 return {partition<T>(mesh, cell_adj, std::span(new_coords), xshape,
453 redistribute, ghost_mode),
469template <std::
floating_po
int T>
470std::tuple<mesh::Mesh<T>, std::vector<std::int32_t>, std::vector<std::int8_t>>
472 bool redistribute,
Option option)
480 mesh.geometry().cmap(), new_vertex_coords, xshape,
481 mesh::GhostMode::none),
486 std::shared_ptr<const common::IndexMap> map_c
487 = mesh.topology()->index_map(mesh.topology()->dim());
488 const int num_ghost_cells = map_c->num_ghosts();
491 int max_ghost_cells = 0;
492 MPI_Allreduce(&num_ghost_cells, &max_ghost_cells, 1, MPI_INT, MPI_MAX,
497 ? mesh::GhostMode::none
498 : mesh::GhostMode::shared_facet;
500 return {partition<T>(mesh, cell_adj, new_vertex_coords, xshape,
501 redistribute, ghost_mode),
514template <std::
floating_po
int T>
515std::tuple<graph::AdjacencyList<std::int64_t>, std::vector<T>,
516 std::array<std::size_t, 2>, std::vector<std::int32_t>,
517 std::vector<std::int8_t>>
521 auto topology = mesh.topology();
524 if (topology->cell_type() != mesh::CellType::triangle
525 and topology->cell_type() != mesh::CellType::tetrahedron)
527 throw std::runtime_error(
"Cell type not supported");
530 auto map_e = topology->index_map(1);
532 throw std::runtime_error(
"Edges must be initialised");
539 std::vector<int> ranks(edge_ranks.
array().begin(), edge_ranks.
array().end());
540 std::sort(ranks.begin(), ranks.end());
541 ranks.erase(std::unique(ranks.begin(), ranks.end()), ranks.end());
544 std::transform(edge_ranks.
array().begin(), edge_ranks.
array().end(),
545 edge_ranks.
array().begin(),
548 auto it = std::lower_bound(ranks.begin(), ranks.end(), r);
549 assert(it != ranks.end() and *it == r);
550 return std::distance(ranks.begin(), it);
554 MPI_Dist_graph_create_adjacent(mesh.comm(), ranks.size(), ranks.data(),
555 MPI_UNWEIGHTED, ranks.size(), ranks.data(),
556 MPI_UNWEIGHTED, MPI_INFO_NULL,
false, &comm);
558 const auto [long_edge, edge_ratio_ok] = impl::face_long_edge(mesh);
560 = impl::compute_refinement(
562 std::vector<std::int8_t>(map_e->size_local() + map_e->num_ghosts(),
564 edge_ranks, mesh, long_edge, edge_ratio_ok, option);
565 MPI_Comm_free(&comm);
567 return {std::move(cell_adj), std::move(new_vertex_coords), xshape,
580template <std::
floating_po
int T>
581std::tuple<graph::AdjacencyList<std::int64_t>, std::vector<T>,
582 std::array<std::size_t, 2>, std::vector<std::int32_t>,
583 std::vector<std::int8_t>>
585 std::span<const std::int32_t> edges,
Option option)
588 auto topology = mesh.topology();
591 if (topology->cell_type() != mesh::CellType::triangle
592 and topology->cell_type() != mesh::CellType::tetrahedron)
594 throw std::runtime_error(
"Cell type not supported");
597 auto map_e = topology->index_map(1);
599 throw std::runtime_error(
"Edges must be initialised");
606 std::vector<int> ranks(edge_ranks.
array().begin(), edge_ranks.
array().end());
607 std::sort(ranks.begin(), ranks.end());
608 ranks.erase(std::unique(ranks.begin(), ranks.end()), ranks.end());
611 std::transform(edge_ranks.
array().begin(), edge_ranks.
array().end(),
612 edge_ranks.
array().begin(),
615 auto it = std::lower_bound(ranks.begin(), ranks.end(), r);
616 assert(it != ranks.end() and *it == r);
617 return std::distance(ranks.begin(), it);
621 std::vector<std::int8_t> marked_edges(
622 map_e->size_local() + map_e->num_ghosts(),
false);
623 std::vector<std::vector<std::int32_t>> marked_for_update(ranks.size());
624 for (
auto edge : edges)
626 if (!marked_edges[edge])
628 marked_edges[edge] =
true;
631 for (
int rank : edge_ranks.
links(edge))
632 marked_for_update[rank].push_back(edge);
637 MPI_Dist_graph_create_adjacent(mesh.comm(), ranks.size(), ranks.data(),
638 MPI_UNWEIGHTED, ranks.size(), ranks.data(),
639 MPI_UNWEIGHTED, MPI_INFO_NULL,
false, &comm);
646 const auto [long_edge, edge_ratio_ok] = impl::face_long_edge(mesh);
647 impl::enforce_rules(comm, edge_ranks, marked_edges, *topology, long_edge);
650 = impl::compute_refinement(comm, marked_edges, edge_ranks, mesh,
651 long_edge, edge_ratio_ok, option);
652 MPI_Comm_free(&comm);
654 return {std::move(cell_adj), std::move(new_vertex_coords), xshape,