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