DOLFINx 0.10.0.0
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 "dolfinx/graph/AdjacencyList.h"
8#include "dolfinx/mesh/Mesh.h"
9#include "dolfinx/mesh/Topology.h"
10#include "dolfinx/mesh/utils.h"
11#include "option.h"
12#include "utils.h"
13#include <algorithm>
14#include <cmath>
15#include <cstdint>
16#include <optional>
17#include <span>
18#include <tuple>
19#include <utility>
20#include <vector>
21
22#pragma once
23
30{
31namespace impl
32{
36template <int tdim>
37auto compute_parent_facets(std::span<const std::int32_t> simplex_set)
38{
39 static_assert(tdim == 2 or tdim == 3);
40 assert(simplex_set.size() % (tdim + 1) == 0);
41 using parent_facet_t
42 = std::conditional_t<tdim == 2, std::array<std::int8_t, 12>,
43 std::array<std::int8_t, 32>>;
44 parent_facet_t parent_facet;
45 parent_facet.fill(-1);
46 assert(simplex_set.size() <= parent_facet.size());
47
48 // Index lookups in 'indices' for the child vertices that occur on
49 // each parent facet in 2D and 3D. In 2D each edge has 3 child
50 // vertices, and in 3D each triangular facet has six child vertices.
51 constexpr std::array<std::array<int, 3>, 3> facet_table_2d{
52 {{1, 2, 3}, {0, 2, 4}, {0, 1, 5}}};
53
54 constexpr std::array<std::array<int, 6>, 4> facet_table_3d{
55 {{1, 2, 3, 4, 5, 6},
56 {0, 2, 3, 4, 7, 8},
57 {0, 1, 3, 5, 7, 9},
58 {0, 1, 2, 6, 8, 9}}};
59
60 const int ncells = simplex_set.size() / (tdim + 1);
61 for (int fpi = 0; fpi < (tdim + 1); ++fpi)
62 {
63 // For each child cell, consider all facets
64 for (int cc = 0; cc < ncells; ++cc)
65 {
66 for (int fci = 0; fci < (tdim + 1); ++fci)
67 {
68 // Indices of all vertices on child facet, sorted
69 std::array<int, tdim> cf, set_output;
70
71 int num_common_vertices;
72 if constexpr (tdim == 2)
73 {
74 for (int j = 0; j < tdim; ++j)
75 cf[j] = simplex_set[cc * 3 + facet_table_2d[fci][j]];
76
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);
81 }
82 else
83 {
84 for (int j = 0; j < tdim; ++j)
85 cf[j] = simplex_set[cc * 4 + facet_table_3d[fci][j]];
86
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);
91 }
92
93 if (num_common_vertices == tdim)
94 {
95 assert(parent_facet[cc * (tdim + 1) + fci] == -1);
96 // Child facet "fci" of cell cc, lies on parent facet "fpi"
97 parent_facet[cc * (tdim + 1) + fci] = fpi;
98 }
99 }
100 }
101 }
102
103 return parent_facet;
104}
105
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,
127 bool uniform);
128
131void enforce_rules(MPI_Comm comm, const graph::AdjacencyList<int>& shared_edges,
132 std::span<std::int8_t> marked_edges,
133 const mesh::Topology& topology,
134 std::span<const std::int32_t> long_edge);
135
145template <std::floating_point T>
146std::pair<std::vector<std::int32_t>, std::vector<std::int8_t>>
147face_long_edge(const mesh::Mesh<T>& mesh)
148{
149 const int tdim = mesh.topology()->dim();
150 // FIXME: cleanup these calls? Some of the happen internally again.
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);
156
157 std::int64_t num_faces = mesh.topology()->index_map(2)->size_local()
158 + mesh.topology()->index_map(2)->num_ghosts();
159
160 // Storage for face-local index of longest edge
161 std::vector<std::int32_t> long_edge(num_faces);
162 std::vector<std::int8_t> edge_ratio_ok;
163
164 // Check mesh face quality (may be used in 2D to switch to "uniform"
165 // refinement)
166 const T min_ratio = std::sqrt(2.0) / 2.0;
167 if (tdim == 2)
168 edge_ratio_ok.resize(num_faces);
169
170 auto x_dofmap = mesh.geometry().dofmap();
171
172 auto c_to_v = mesh.topology()->connectivity(tdim, 0);
173 assert(c_to_v);
174 auto e_to_c = mesh.topology()->connectivity(1, tdim);
175 assert(e_to_c);
176 auto e_to_v = mesh.topology()->connectivity(1, 0);
177 assert(e_to_v);
178
179 // Store all edge lengths in Mesh to save recalculating for each Face
180 auto map_e = mesh.topology()->index_map(1);
181 assert(map_e);
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)
184 {
185 // Get first attached cell
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);
190
191 // Find local index of edge vertices in the cell geometry map
192 auto it0 = std::find(cell_vertices.begin(), cell_vertices.end(),
193 edge_vertices[0]);
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(),
197 edge_vertices[1]);
198 assert(it1 != cell_vertices.end());
199 const std::size_t local1 = std::distance(cell_vertices.begin(), it1);
200
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],
203 3);
204 std::span<const T, 3> x1(mesh.geometry().x().data() + 3 * x_dofs[local1],
205 3);
206
207 // Compute length of edge between vertex x0 and x1
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); }));
211 }
212
213 // Get longest edge of each face
214 auto f_to_v = mesh.topology()->connectivity(2, 0);
215 assert(f_to_v);
216 auto f_to_e = mesh.topology()->connectivity(2, 1);
217 assert(f_to_e);
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)
221 {
222 auto face_edges = f_to_e->links(f);
223
224 std::int32_t imax = 0;
225 T max_len = 0.0;
226 T min_len = std::numeric_limits<T>::max();
227
228 for (int i = 0; i < 3; ++i)
229 {
230 const T e_len = edge_length[face_edges[i]];
231 min_len = std::min(e_len, min_len);
232 if (e_len > max_len)
233 {
234 max_len = e_len;
235 imax = i;
236 }
237 else if (tdim == 3 and e_len == max_len)
238 {
239 // If edges are the same length, compare global index of
240 // opposite vertex. Only important so that tetrahedral faces
241 // have a matching refinement pattern across processes.
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])
246 imax = i;
247 }
248 }
249
250 // Only save edge ratio in 2D
251 if (tdim == 2)
252 edge_ratio_ok[f] = (min_len / max_len >= min_ratio);
253
254 long_edge[f] = face_edges[imax];
255 }
256
257 return std::pair(std::move(long_edge), std::move(edge_ratio_ok));
258}
259
284template <std::floating_point 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,
290 const graph::AdjacencyList<int>& shared_edges,
291 const mesh::Mesh<T>& mesh,
292 std::span<const std::int32_t> long_edge,
293 std::span<const std::int8_t> edge_ratio_ok, Option option)
294{
295 int tdim = mesh.topology()->dim();
296 int num_cell_edges = tdim * 3 - 3;
297 int num_cell_vertices = tdim + 1;
298
299 bool compute_facets = option_parent_facet(option);
300 bool compute_parent_cell = option_parent_cell(option);
301
302 // Make new vertices in parallel
303 const auto [new_vertex_map, new_vertex_coords, xshape]
304 = create_new_vertices(neighbor_comm, shared_edges, mesh, marked_edges);
305
306 std::optional<std::vector<std::int32_t>> parent_cell(std::nullopt);
307 if (compute_parent_cell)
308 parent_cell.emplace();
309
310 std::optional<std::vector<std::int8_t>> parent_facet(std::nullopt);
311 if (compute_facets)
312 parent_facet.emplace();
313
314 std::vector<std::int64_t> indices(num_cell_vertices + num_cell_edges);
315 std::vector<std::int32_t> simplex_set;
316
317 auto map_c = mesh.topology()->index_map(tdim);
318 assert(map_c);
319 auto c_to_v = mesh.topology()->connectivity(tdim, 0);
320 assert(c_to_v);
321 auto c_to_e = mesh.topology()->connectivity(tdim, 1);
322 assert(c_to_e);
323 auto c_to_f = mesh.topology()->connectivity(tdim, 2);
324 assert(c_to_f);
325
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);
329
330 std::vector<std::int64_t> global_indices
331 = adjust_indices(*mesh.topology()->index_map(0), num_new_vertices_local);
332
333 const std::int32_t num_cells = map_c->size_local();
334
335 // Iterate over all cells, and refine if cell has a marked edge
336 std::vector<std::int64_t> cell_topology;
337 for (int c = 0; c < num_cells; ++c)
338 {
339 // Create vector of indices in the order [vertices][edges], 3+3 in
340 // 2D, 4+6 in 3D
341
342 // Copy vertices
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]];
346
347 // Get cell-local indices of marked edges
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)
351 {
352 if (marked_edges[edges[ei]])
353 {
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;
358 }
359 else
360 indices[num_cell_vertices + ei] = -1;
361 }
362
363 if (no_edge_marked)
364 {
365 // Copy over existing cell to new topology
366 for (auto v : vertices)
367 cell_topology.push_back(global_indices[v]);
368
369 if (compute_parent_cell)
370 parent_cell->push_back(c);
371
372 if (compute_facets)
373 {
374 if (tdim == 3)
375 parent_facet->insert(parent_facet->end(), {0, 1, 2, 3});
376 else
377 parent_facet->insert(parent_facet->end(), {0, 1, 2});
378 }
379 }
380 else
381 {
382 // Need longest edges of each face in cell local indexing. NB in
383 // 2D the face is the cell itself, and there is just one entry.
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]);
387
388 // Convert to cell local index
389 for (std::int32_t& p : longest_edge)
390 {
391 for (std::size_t ej = 0; ej < edges.size(); ++ej)
392 {
393 if (p == edges[ej])
394 {
395 p = ej;
396 break;
397 }
398 }
399 }
400
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(),
405 simplex_set_size);
406
407 // Save parent index
408 const std::int32_t ncells = simplex_set.size() / num_cell_vertices;
409 if (compute_parent_cell)
410 {
411 for (std::int32_t i = 0; i < ncells; ++i)
412 parent_cell->push_back(c);
413 }
414
415 if (compute_facets)
416 {
417 if (tdim == 3)
418 {
419 auto npf = compute_parent_facets<3>(simplex_set);
420 parent_facet->insert(parent_facet->end(), npf.begin(),
421 std::next(npf.begin(), simplex_set.size()));
422 }
423 else
424 {
425 auto npf = compute_parent_facets<2>(simplex_set);
426 parent_facet->insert(parent_facet->end(), npf.begin(),
427 std::next(npf.begin(), simplex_set.size()));
428 }
429 }
430
431 // Convert from cell local index to mesh index and add to cells
432 for (std::int32_t v : simplex_set)
433 cell_topology.push_back(indices[v]);
434 }
435 }
436
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;
442 graph::AdjacencyList cell_adj(std::move(cell_topology), std::move(offsets));
443
444 return {std::move(cell_adj), std::move(new_vertex_coords), xshape,
445 std::move(parent_cell), std::move(parent_facet)};
446}
447} // namespace impl
448
458template <std::floating_point 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,
464 Option option)
465{
466 common::Timer t0("PLAZA: refine");
467 auto topology = mesh.topology();
468 assert(topology);
469
470 if (topology->cell_type() != mesh::CellType::triangle
471 and topology->cell_type() != mesh::CellType::tetrahedron)
472 {
473 throw std::runtime_error("Cell type not supported");
474 }
475
476 auto map_e = topology->index_map(1);
477 if (!map_e)
478 throw std::runtime_error("Edges must be initialised");
479
480 // Get sharing ranks for each edge
481 graph::AdjacencyList<int> edge_ranks = map_e->index_to_dest_ranks();
482
483 // Create unique list of ranks that share edges (owners of ghosts plus
484 // ranks that ghost owned indices)
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);
489
490 // Convert edge_ranks from global rank to to neighbourhood ranks
491 std::ranges::transform(edge_ranks.array(), edge_ranks.array().begin(),
492 [&ranks](auto r)
493 {
494 auto it = std::ranges::lower_bound(ranks, r);
495 assert(it != ranks.end() and *it == r);
496 return std::distance(ranks.begin(), it);
497 });
498
499 // Get number of neighbors
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());
503
504 if (edges)
505 {
506 for (auto edge : edges.value())
507 {
508 if (!marked_edges[edge])
509 {
510 marked_edges[edge] = true;
511
512 // If it is a shared edge, add all sharing neighbors to update set
513 for (int rank : edge_ranks.links(edge))
514 marked_for_update[rank].push_back(edge);
515 }
516 }
517 }
518
519 MPI_Comm comm;
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);
523
524 // Communicate any shared edges
525 update_logical_edgefunction(comm, marked_for_update, marked_edges, *map_e);
526
527 // Enforce rules about refinement (i.e. if any edge is marked in a
528 // triangle, then the longest edge must also be marked).
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);
531
532 auto [cell_adj, new_vertex_coords, xshape, parent_cell, parent_facet]
533 = impl::compute_refinement(comm, marked_edges, edge_ranks, mesh,
534 long_edge, edge_ratio_ok, option);
535 MPI_Comm_free(&comm);
536
537 return {std::move(cell_adj), std::move(new_vertex_coords), xshape,
538 std::move(parent_cell), std::move(parent_facet)};
539}
540
541} // namespace dolfinx::refinement::plaza
Timer for measuring and logging elapsed time durations.
Definition Timer.h:41
Definition AdjacencyList.h:27
std::span< T > links(std::size_t node)
Definition AdjacencyList.h:111
const std::vector< T > & array() const
Return contiguous array of links for all nodes (const version)
Definition AdjacencyList.h:128
A Mesh consists of a set of connected and numbered mesh topological entities, and geometry data.
Definition Mesh.h:23
Topology stores the topology of a mesh, consisting of mesh entities and connectivity (incidence relat...
Definition Topology.h:46
Functions supporting mesh operations.
Mesh data structures and algorithms on meshes.
Definition DofMap.h:32
Plaza mesh refinement.
Definition plaza.h:30
std::tuple< graph::AdjacencyList< std::int64_t >, std::vector< T >, 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 > > edges, Option option)
Definition plaza.h:462
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
@ parent_cell
Definition option.h:20
@ parent_facet
Definition option.h:18
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:149
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