DOLFINx
DOLFINx C++ interface
Loading...
Searching...
No Matches
plaza.h
1// Copyright (C) 2014-2018 Chris Richardson
2//
3// This file is part of DOLFINx (https://www.fenicsproject.org)
4//
5// SPDX-License-Identifier: LGPL-3.0-or-later
6
7#include "utils.h"
8#include <cstdint>
9#include <dolfinx/common/MPI.h>
10#include <dolfinx/graph/AdjacencyList.h>
11#include <dolfinx/mesh/Geometry.h>
12#include <dolfinx/mesh/Mesh.h>
13#include <dolfinx/mesh/Topology.h>
14#include <dolfinx/mesh/utils.h>
15#include <span>
16#include <tuple>
17#include <utility>
18#include <vector>
19
20#pragma once
21
28{
29
31enum class Option : int
32{
33 none = 0,
35 = 1,
37 = 2,
40};
41
42namespace impl
43{
47template <int tdim>
48std::vector<std::int8_t>
49compute_parent_facets(std::span<const std::int32_t> simplex_set)
50{
51 {
52 assert(simplex_set.size() % (tdim + 1) == 0);
53 std::vector<std::int8_t> parent_facet(simplex_set.size(), -1);
54
55 // Index lookups in 'indices' for the child vertices that occur on
56 // each parent facet in 2D and 3D. In 2D each edge has 3 child
57 // vertices, and in 3D each triangular facet has six child vertices.
58 constexpr std::array<std::array<int, 3>, 3> facet_table_2d{
59 {{1, 2, 3}, {0, 2, 4}, {0, 1, 5}}};
60
61 constexpr std::array<std::array<int, 6>, 4> facet_table_3d{
62 {{1, 2, 3, 4, 5, 6},
63 {0, 2, 3, 4, 7, 8},
64 {0, 1, 3, 5, 7, 9},
65 {0, 1, 2, 6, 8, 9}}};
66
67 const int ncells = simplex_set.size() / (tdim + 1);
68 for (int fpi = 0; fpi < (tdim + 1); ++fpi)
69 {
70 // For each child cell, consider all facets
71 for (int cc = 0; cc < ncells; ++cc)
72 {
73 for (int fci = 0; fci < (tdim + 1); ++fci)
74 {
75 // Indices of all vertices on child facet, sorted
76 std::array<int, tdim> cf, set_output;
77
78 int num_common_vertices;
79 if constexpr (tdim == 2)
80 {
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);
88 }
89 else
90 {
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);
98 }
99
100 if (num_common_vertices == tdim)
101 {
102 assert(parent_facet[cc * (tdim + 1) + fci] == -1);
103 // Child facet "fci" of cell cc, lies on parent facet "fpi"
104 parent_facet[cc * (tdim + 1) + fci] = fpi;
105 }
106 }
107 }
108 }
109
110 return parent_facet;
111 }
112}
113
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,
130 bool uniform);
131
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);
138
140template <std::floating_point T>
141std::pair<std::vector<std::int32_t>, std::vector<std::int8_t>>
142face_long_edge(const mesh::Mesh<T>& mesh)
143{
144 const int tdim = mesh.topology()->dim();
145 // FIXME: cleanup these calls? Some of the happen internally again.
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);
151
152 std::int64_t num_faces = mesh.topology()->index_map(2)->size_local()
153 + mesh.topology()->index_map(2)->num_ghosts();
154
155 // Storage for face-local index of longest edge
156 std::vector<std::int32_t> long_edge(num_faces);
157 std::vector<std::int8_t> edge_ratio_ok;
158
159 // Check mesh face quality (may be used in 2D to switch to "uniform"
160 // refinement)
161 const T min_ratio = sqrt(2.0) / 2.0;
162 if (tdim == 2)
163 edge_ratio_ok.resize(num_faces);
164
165 auto x_dofmap = mesh.geometry().dofmap();
166
167 auto c_to_v = mesh.topology()->connectivity(tdim, 0);
168 assert(c_to_v);
169 auto e_to_c = mesh.topology()->connectivity(1, tdim);
170 assert(e_to_c);
171 auto e_to_v = mesh.topology()->connectivity(1, 0);
172 assert(e_to_v);
173
174 // Store all edge lengths in Mesh to save recalculating for each Face
175 auto map_e = mesh.topology()->index_map(1);
176 assert(map_e);
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)
179 {
180 // Get first attached cell
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);
185
186 // Find local index of edge vertices in the cell geometry map
187 auto it0 = std::find(cell_vertices.begin(), cell_vertices.end(),
188 edge_vertices[0]);
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(),
192 edge_vertices[1]);
193 assert(it1 != cell_vertices.end());
194 const std::size_t local1 = std::distance(cell_vertices.begin(), it1);
195
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],
199 3);
200 std::span<const T, 3> x1(mesh.geometry().x().data() + 3 * x_dofs[local1],
201 3);
202
203 // Compute length of edge between vertex x0 and x1
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); }));
207 }
208
209 // Get longest edge of each face
210 auto f_to_v = mesh.topology()->connectivity(2, 0);
211 assert(f_to_v);
212 auto f_to_e = mesh.topology()->connectivity(2, 1);
213 assert(f_to_e);
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)
217 {
218 auto face_edges = f_to_e->links(f);
219
220 std::int32_t imax = 0;
221 T max_len = 0.0;
222 T min_len = std::numeric_limits<T>::max();
223
224 for (int i = 0; i < 3; ++i)
225 {
226 const T e_len = edge_length[face_edges[i]];
227 min_len = std::min(e_len, min_len);
228 if (e_len > max_len)
229 {
230 max_len = e_len;
231 imax = i;
232 }
233 else if (tdim == 3 and e_len == max_len)
234 {
235 // If edges are the same length, compare global index of
236 // opposite vertex. Only important so that tetrahedral faces
237 // have a matching refinement pattern across processes.
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])
242 imax = i;
243 }
244 }
245
246 // Only save edge ratio in 2D
247 if (tdim == 2)
248 edge_ratio_ok[f] = (min_len / max_len >= min_ratio);
249
250 long_edge[f] = face_edges[imax];
251 }
252
253 return std::pair(std::move(long_edge), std::move(edge_ratio_ok));
254}
255
257template <std::floating_point 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,
267 plaza::Option option)
268{
269 int tdim = mesh.topology()->dim();
270 int num_cell_edges = tdim * 3 - 3;
271 int num_cell_vertices = tdim + 1;
272 bool compute_facets = (option == plaza::Option::parent_facet
274 bool compute_parent_cell
275 = (option == plaza::Option::parent_cell
277
278 // Make new vertices in parallel
279 const auto [new_vertex_map, new_vertex_coords, xshape]
280 = create_new_vertices(neighbor_comm, shared_edges, mesh, marked_edges);
281
282 std::vector<std::int32_t> parent_cell;
283 std::vector<std::int8_t> parent_facet;
284 std::vector<std::int64_t> indices(num_cell_vertices + num_cell_edges);
285 std::vector<std::int32_t> simplex_set;
286
287 auto map_c = mesh.topology()->index_map(tdim);
288 assert(map_c);
289 auto c_to_v = mesh.topology()->connectivity(tdim, 0);
290 assert(c_to_v);
291 auto c_to_e = mesh.topology()->connectivity(tdim, 1);
292 assert(c_to_e);
293 auto c_to_f = mesh.topology()->connectivity(tdim, 2);
294 assert(c_to_f);
295
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);
299
300 std::vector<std::int64_t> global_indices
301 = adjust_indices(*mesh.topology()->index_map(0), num_new_vertices_local);
302
303 const int num_cells = map_c->size_local();
304
305 // Iterate over all cells, and refine if cell has a marked edge
306 std::vector<std::int64_t> cell_topology;
307 for (int c = 0; c < num_cells; ++c)
308 {
309 // Create vector of indices in the order [vertices][edges], 3+3 in
310 // 2D, 4+6 in 3D
311
312 // Copy vertices
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]];
316
317 // Get cell-local indices of marked edges
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)
321 {
322 if (marked_edges[edges[ei]])
323 {
324 no_edge_marked = false;
325 auto it = new_vertex_map.find(edges[ei]);
326 assert(it != new_vertex_map.end());
327 indices[num_cell_vertices + ei] = it->second;
328 }
329 else
330 indices[num_cell_vertices + ei] = -1;
331 }
332
333 if (no_edge_marked)
334 {
335 // Copy over existing cell to new topology
336 for (auto v : vertices)
337 cell_topology.push_back(global_indices[v]);
338
339 if (compute_parent_cell)
340 parent_cell.push_back(c);
341
342 if (compute_facets)
343 {
344 if (tdim == 3)
345 parent_facet.insert(parent_facet.end(), {0, 1, 2, 3});
346 else
347 parent_facet.insert(parent_facet.end(), {0, 1, 2});
348 }
349 }
350 else
351 {
352 // Need longest edges of each face in cell local indexing. NB in
353 // 2D the face is the cell itself, and there is just one entry.
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]);
357
358 // Convert to cell local index
359 for (std::int32_t& p : longest_edge)
360 {
361 for (std::size_t ej = 0; ej < edges.size(); ++ej)
362 {
363 if (p == edges[ej])
364 {
365 p = ej;
366 break;
367 }
368 }
369 }
370
371 const bool uniform = (tdim == 2) ? edge_ratio_ok[c] : false;
372
373 // FIXME: this has an expensive dynamic memory allocation
374 simplex_set = get_simplices(indices, longest_edge, tdim, uniform);
375
376 // Save parent index
377 const std::int32_t ncells = simplex_set.size() / num_cell_vertices;
378
379 if (compute_parent_cell)
380 {
381 for (std::int32_t i = 0; i < ncells; ++i)
382 parent_cell.push_back(c);
383 }
384
385 if (compute_facets)
386 {
387 std::vector<std::int8_t> npf;
388 if (tdim == 3)
389 npf = compute_parent_facets<3>(simplex_set);
390 else
391 npf = compute_parent_facets<2>(simplex_set);
392 parent_facet.insert(parent_facet.end(), npf.begin(), npf.end());
393 }
394
395 // Convert from cell local index to mesh index and add to cells
396 for (std::int32_t v : simplex_set)
397 cell_topology.push_back(indices[v]);
398 }
399 }
400
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));
407
408 return {std::move(cell_adj), std::move(new_vertex_coords), xshape,
409 std::move(parent_cell), std::move(parent_facet)};
410}
411} // namespace impl
412
423template <std::floating_point T>
424std::tuple<mesh::Mesh<T>, std::vector<std::int32_t>, std::vector<std::int8_t>>
425refine(const mesh::Mesh<T>& mesh, bool redistribute, Option option)
426{
427 auto [cell_adj, new_coords, xshape, parent_cell, parent_facet]
428 = compute_refinement_data(mesh, option);
429
430 if (dolfinx::MPI::size(mesh.comm()) == 1)
431 {
432 return {mesh::create_mesh(mesh.comm(), cell_adj.array(),
433 mesh.geometry().cmap(), new_coords, xshape,
434 mesh::GhostMode::none),
435 std::move(parent_cell), std::move(parent_facet)};
436 }
437 else
438 {
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();
442 // Check if mesh has ghost cells on any rank
443 // FIXME: this is not a robust test. Should be user option.
444 int max_ghost_cells = 0;
445 MPI_Allreduce(&num_ghost_cells, &max_ghost_cells, 1, MPI_INT, MPI_MAX,
446 mesh.comm());
447
448 // Build mesh
449 const mesh::GhostMode ghost_mode = max_ghost_cells == 0
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),
454 std::move(parent_cell), std::move(parent_facet)};
455 }
456}
457
469template <std::floating_point T>
470std::tuple<mesh::Mesh<T>, std::vector<std::int32_t>, std::vector<std::int8_t>>
471refine(const mesh::Mesh<T>& mesh, std::span<const std::int32_t> edges,
472 bool redistribute, Option option)
473{
474 auto [cell_adj, new_vertex_coords, xshape, parent_cell, parent_facet]
475 = compute_refinement_data(mesh, edges, option);
476
477 if (dolfinx::MPI::size(mesh.comm()) == 1)
478 {
479 return {mesh::create_mesh(mesh.comm(), cell_adj.array(),
480 mesh.geometry().cmap(), new_vertex_coords, xshape,
481 mesh::GhostMode::none),
482 std::move(parent_cell), std::move(parent_facet)};
483 }
484 else
485 {
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();
489 // Check if mesh has ghost cells on any rank
490 // FIXME: this is not a robust test. Should be user option.
491 int max_ghost_cells = 0;
492 MPI_Allreduce(&num_ghost_cells, &max_ghost_cells, 1, MPI_INT, MPI_MAX,
493 mesh.comm());
494
495 // Build mesh
496 const mesh::GhostMode ghost_mode = max_ghost_cells == 0
497 ? mesh::GhostMode::none
498 : mesh::GhostMode::shared_facet;
499
500 return {partition<T>(mesh, cell_adj, new_vertex_coords, xshape,
501 redistribute, ghost_mode),
502 std::move(parent_cell), std::move(parent_facet)};
503 }
504}
505
514template <std::floating_point 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>>
519{
520 common::Timer t0("PLAZA: refine");
521 auto topology = mesh.topology();
522 assert(topology);
523
524 if (topology->cell_type() != mesh::CellType::triangle
525 and topology->cell_type() != mesh::CellType::tetrahedron)
526 {
527 throw std::runtime_error("Cell type not supported");
528 }
529
530 auto map_e = topology->index_map(1);
531 if (!map_e)
532 throw std::runtime_error("Edges must be initialised");
533
534 // Get sharing ranks for each edge
535 graph::AdjacencyList<int> edge_ranks = map_e->index_to_dest_ranks();
536
537 // Create unique list of ranks that share edges (owners of ghosts
538 // plus ranks that ghost owned indices)
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());
542
543 // Convert edge_ranks from global rank to to neighbourhood ranks
544 std::transform(edge_ranks.array().begin(), edge_ranks.array().end(),
545 edge_ranks.array().begin(),
546 [&ranks](auto r)
547 {
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);
551 });
552
553 MPI_Comm comm;
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);
557
558 const auto [long_edge, edge_ratio_ok] = impl::face_long_edge(mesh);
559 auto [cell_adj, new_vertex_coords, xshape, parent_cell, parent_facet]
560 = impl::compute_refinement(
561 comm,
562 std::vector<std::int8_t>(map_e->size_local() + map_e->num_ghosts(),
563 true),
564 edge_ranks, mesh, long_edge, edge_ratio_ok, option);
565 MPI_Comm_free(&comm);
566
567 return {std::move(cell_adj), std::move(new_vertex_coords), xshape,
568 std::move(parent_cell), std::move(parent_facet)};
569}
570
580template <std::floating_point 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)
586{
587 common::Timer t0("PLAZA: refine");
588 auto topology = mesh.topology();
589 assert(topology);
590
591 if (topology->cell_type() != mesh::CellType::triangle
592 and topology->cell_type() != mesh::CellType::tetrahedron)
593 {
594 throw std::runtime_error("Cell type not supported");
595 }
596
597 auto map_e = topology->index_map(1);
598 if (!map_e)
599 throw std::runtime_error("Edges must be initialised");
600
601 // Get sharing ranks for each edge
602 graph::AdjacencyList<int> edge_ranks = map_e->index_to_dest_ranks();
603
604 // Create unique list of ranks that share edges (owners of ghosts plus
605 // ranks that ghost owned indices)
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());
609
610 // Convert edge_ranks from global rank to to neighbourhood ranks
611 std::transform(edge_ranks.array().begin(), edge_ranks.array().end(),
612 edge_ranks.array().begin(),
613 [&ranks](auto r)
614 {
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);
618 });
619
620 // Get number of neighbors
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)
625 {
626 if (!marked_edges[edge])
627 {
628 marked_edges[edge] = true;
629
630 // If it is a shared edge, add all sharing neighbors to update set
631 for (int rank : edge_ranks.links(edge))
632 marked_for_update[rank].push_back(edge);
633 }
634 }
635
636 MPI_Comm comm;
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);
640
641 // Communicate any shared edges
642 update_logical_edgefunction(comm, marked_for_update, marked_edges, *map_e);
643
644 // Enforce rules about refinement (i.e. if any edge is marked in a
645 // triangle, then the longest edge must also be marked).
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);
648
649 auto [cell_adj, new_vertex_coords, xshape, parent_cell, parent_facet]
650 = impl::compute_refinement(comm, marked_edges, edge_ranks, mesh,
651 long_edge, edge_ratio_ok, option);
652 MPI_Comm_free(&comm);
653
654 return {std::move(cell_adj), std::move(new_vertex_coords), xshape,
655 std::move(parent_cell), std::move(parent_facet)};
656}
657
658} // namespace dolfinx::refinement::plaza
A timer can be used for timing tasks. The basic usage is.
Definition Timer.h:31
This class provides a static adjacency list data structure. It is commonly used to store directed gra...
Definition AdjacencyList.h:28
std::span< T > links(std::size_t node)
Get the links (edges) for given node.
Definition AdjacencyList.h:112
const std::vector< T > & array() const
Return contiguous array of links for all nodes (const version)
Definition AdjacencyList.h:129
A Mesh consists of a set of connected and numbered mesh topological entities, and geometry data.
Definition Mesh.h:23
Functions supporting mesh operations.
int size(MPI_Comm comm)
Return size of the group (number of processes) associated with the communicator.
Definition MPI.cpp:72
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:15
GhostMode
Enum for different partitioning ghost modes.
Definition utils.h:35
Mesh< typename std::remove_reference_t< typename U::value_type > > create_mesh(MPI_Comm comm, MPI_Comm commt, std::span< const std::int64_t > cells, const fem::CoordinateElement< typename std::remove_reference_t< typename U::value_type > > &element, MPI_Comm commg, const U &x, std::array< std::size_t, 2 > xshape, const CellPartitionFunction &partitioner)
Create a distributed mesh from mesh data using a provided graph partitioning function for determining...
Definition utils.h:785
int num_cell_vertices(CellType type)
Number vertices for a cell type.
Definition cell_types.cpp:147
Plaza mesh refinement.
Definition plaza.h:28
std::tuple< graph::AdjacencyList< std::int64_t >, std::vector< T >, std::array< std::size_t, 2 >, std::vector< std::int32_t >, std::vector< std::int8_t > > compute_refinement_data(const mesh::Mesh< T > &mesh, Option option)
Refine mesh returning new mesh data.
Definition plaza.h:518
Option
Options for data to compute during mesh refinement.
Definition plaza.h:32
std::tuple< mesh::Mesh< T >, std::vector< std::int32_t >, std::vector< std::int8_t > > refine(const mesh::Mesh< T > &mesh, bool redistribute, Option option)
Uniform refine, optionally redistributing and optionally calculating the parent-child relationships`.
Definition plaza.h:425
void update_logical_edgefunction(MPI_Comm comm, const std::vector< std::vector< std::int32_t > > &marked_for_update, std::vector< std::int8_t > &marked_edges, const common::IndexMap &map)
Communicate edge markers between processes that share edges.
Definition utils.cpp:46
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:147
std::vector< std::int64_t > adjust_indices(const common::IndexMap &map, std::int32_t n)
Add indices to account for extra n values on this process.
Definition utils.cpp:102