8#include <compile_time_options.h>
10#include <deal.II/base/partitioner.h>
11#include <deal.II/dofs/dof_handler.h>
12#include <deal.II/dofs/dof_renumbering.h>
13#include <deal.II/dofs/dof_tools.h>
14#include <deal.II/lac/affine_constraints.h>
15#include <deal.II/lac/dynamic_sparsity_pattern.h>
16#include <deal.II/lac/sparsity_tools.h>
28 template <
typename Number>
30 const dealii::Utilities::MPI::Partitioner &partitioner,
31 dealii::AffineConstraints<Number> &affine_constraints)
33 affine_constraints.close();
35 dealii::AffineConstraints<Number> temporary;
37 for (
auto line : affine_constraints.get_lines()) {
39 line.index = partitioner.global_to_local(line.index);
40 std::transform(line.entries.begin(),
44 return std::make_pair(
45 partitioner.global_to_local(entry.first),
49 temporary.add_line(line.index);
50 temporary.add_entries(line.index, line.entries);
51 temporary.set_inhomogeneity(line.index, line.inhomogeneity);
56 affine_constraints = std::move(temporary);
68 template <
typename VECTOR>
70 const dealii::Utilities::MPI::Partitioner &partitioner, VECTOR &vector)
73 vector.begin(), vector.end(), vector.begin(), [&](
auto index) {
74 return partitioner.global_to_local(index);
85 namespace DoFRenumbering
91 using dealii::DoFRenumbering::Cuthill_McKee;
104 const MPI_Comm &mpi_communicator,
105 const unsigned int n_locally_internal,
106 const std::size_t group_size)
108 using namespace dealii;
110 const IndexSet &locally_owned = dof_handler.locally_owned_dofs();
111 const auto n_locally_owned = locally_owned.n_elements();
114 Assert(locally_owned.is_contiguous() ==
true,
116 "Need a contiguous set of locally owned indices."));
119 const auto offset = n_locally_owned != 0 ? *locally_owned.begin() : 0;
121#if DEAL_II_VERSION_GTE(9, 6, 0)
122 const auto locally_relevant =
123 DoFTools::extract_locally_relevant_dofs(dof_handler);
125 IndexSet locally_relevant;
126 DoFTools::extract_locally_relevant_dofs(dof_handler, locally_relevant);
131 Utilities::MPI::Partitioner partitioner(
132 locally_owned, locally_relevant, mpi_communicator);
134 IndexSet export_indices(n_locally_owned);
135 for (
const auto &it : partitioner.import_indices()) {
136 export_indices.add_range(it.first, it.second);
139 std::vector<dealii::types::global_dof_index> new_order(n_locally_owned);
146 unsigned int n_export_indices = 0;
148 Assert(n_locally_internal <= n_locally_owned, dealii::ExcInternalError());
150 for (
unsigned int i = 0; i < n_locally_internal; i += group_size) {
151 bool export_index_present =
false;
152 for (
unsigned int j = 0; j < group_size; ++j) {
153 if (export_indices.is_element(i + j)) {
154 export_index_present =
true;
159 if (export_index_present) {
160 Assert(n_export_indices % group_size == 0,
161 dealii::ExcInternalError());
162 for (
unsigned int j = 0; j < group_size; ++j) {
163 new_order[i + j] = offset + n_export_indices++;
166 for (
unsigned int j = 0; j < group_size; ++j)
167 new_order[i + j] = dealii::numbers::invalid_dof_index;
172 unsigned int n_other = 0;
173 for (
unsigned int i = n_locally_internal; i < n_locally_owned; ++i)
174 if (export_indices.is_element(i))
177 Assert(n_other + n_export_indices >= export_indices.n_elements(),
178 dealii::ExcInternalError());
181 unsigned int running_index = n_export_indices;
187 for (
unsigned int i = 0; i < n_locally_internal; i += group_size) {
188 if (new_order[i] == dealii::numbers::invalid_dof_index) {
189 for (
unsigned int j = 0; j < group_size; ++j) {
190 Assert(new_order[i + j] == dealii::numbers::invalid_dof_index,
191 dealii::ExcInternalError());
192 new_order[i + j] = offset + running_index++;
197 Assert(running_index == n_locally_internal, dealii::ExcInternalError());
199 for (
unsigned int i = n_locally_internal; i < n_locally_owned; i++) {
200 new_order[i] = offset + running_index++;
203 Assert(running_index == n_locally_owned, dealii::ExcInternalError());
205 dof_handler.renumber_dofs(new_order);
207 Assert(n_export_indices % group_size == 0, dealii::ExcInternalError());
208 Assert(n_export_indices <= n_locally_internal,
209 dealii::ExcInternalError());
210 return n_export_indices;
223 const dealii::DynamicSparsityPattern &sparsity,
224 const unsigned int n_locally_internal,
225 const std::size_t group_size)
227 using namespace dealii;
229 const IndexSet &locally_owned = dof_handler.locally_owned_dofs();
230 const auto n_locally_owned = locally_owned.n_elements();
233 Assert(locally_owned.is_contiguous() ==
true,
235 "Need a contiguous set of locally owned indices."));
238 const auto offset = n_locally_owned != 0 ? *locally_owned.begin() : 0;
240 std::vector<dealii::types::global_dof_index> new_order(n_locally_owned);
248 unsigned int n_consistent_range = 0;
250 Assert(n_locally_internal <= n_locally_owned, dealii::ExcInternalError());
252 for (
unsigned int i = 0; i < n_locally_internal; i += group_size) {
254 bool stride_is_consistent =
true;
255 const auto group_row_length = sparsity.row_length(offset + i);
256 for (
unsigned int j = 0; j < group_size; ++j) {
257 if (group_row_length != sparsity.row_length(offset + i + j)) {
258 stride_is_consistent =
false;
263 if (stride_is_consistent) {
264 for (
unsigned int j = 0; j < group_size; ++j) {
265 new_order[i + j] = offset + n_consistent_range++;
268 for (
unsigned int j = 0; j < group_size; ++j)
269 new_order[i + j] = dealii::numbers::invalid_dof_index;
277 unsigned int running_index = n_consistent_range;
279 for (
unsigned int i = 0; i < n_locally_internal; i += group_size) {
280 if (new_order[i] == dealii::numbers::invalid_dof_index) {
281 for (
unsigned int j = 0; j < group_size; ++j) {
282 Assert(new_order[i + j] == dealii::numbers::invalid_dof_index,
283 dealii::ExcInternalError());
284 new_order[i + j] = offset + running_index++;
289 Assert(running_index == n_locally_internal, dealii::ExcInternalError());
291 for (
unsigned int i = n_locally_internal; i < n_locally_owned; i++) {
292 new_order[i] = offset + running_index++;
295 Assert(running_index == n_locally_owned, dealii::ExcInternalError());
297 dof_handler.renumber_dofs(new_order);
299 Assert(n_consistent_range % group_size == 0, dealii::ExcInternalError());
300 Assert(n_consistent_range <= n_locally_internal,
301 dealii::ExcInternalError());
302 return n_consistent_range;
322 const dealii::DynamicSparsityPattern &sparsity,
323 const std::size_t group_size)
325 using namespace dealii;
327 const auto &locally_owned = dof_handler.locally_owned_dofs();
328 const auto n_locally_owned = locally_owned.n_elements();
332 Assert(locally_owned.is_contiguous() ==
true,
334 "Need a contiguous set of locally owned indices."));
337 const auto offset = n_locally_owned != 0 ? *locally_owned.begin() : 0;
339 using dof_type = dealii::types::global_dof_index;
340 std::vector<dof_type> new_order(n_locally_owned);
341 dof_type current_index = offset;
349 std::map<unsigned int, std::set<dof_type>> bins;
351 for (
unsigned int i = 0; i < n_locally_owned; ++i) {
352 const dof_type index = i;
353 const unsigned int row_length = sparsity.row_length(offset + index);
354 bins[row_length].insert(index);
356 if (bins[row_length].size() == group_size) {
357 for (
const auto &index : bins[row_length])
358 new_order[index] = current_index++;
359 bins.erase(row_length);
363 unsigned int n_locally_internal = current_index - offset;
367 for (
const auto &entries : bins) {
368 Assert(entries.second.size() > 0, ExcInternalError());
369 for (
const auto &index : entries.second)
370 new_order[index] = current_index++;
372 Assert(current_index == offset + n_locally_owned, ExcInternalError());
374 dof_handler.renumber_dofs(new_order);
376 Assert(n_locally_internal % group_size == 0, ExcInternalError());
377 return n_locally_internal;
391 using dealii::DoFTools::extract_locally_relevant_dofs;
394 using dealii::DoFTools::make_hanging_node_constraints;
397 using dealii::DoFTools::make_periodicity_constraints;
400 using dealii::DoFTools::make_sparsity_pattern;
410 template <
int dim,
typename Number,
typename SPARSITY>
412 const dealii::DoFHandler<dim> &dof_handler,
414 const dealii::AffineConstraints<Number> &affine_constraints,
415 bool keep_constrained)
417 std::vector<dealii::types::global_dof_index> dof_indices;
419 for (
auto cell : dof_handler.active_cell_iterators()) {
421 if (cell->is_artificial())
424 const unsigned int dofs_per_cell = cell->get_fe().n_dofs_per_cell();
425 dof_indices.resize(dofs_per_cell);
426 cell->get_dof_indices(dof_indices);
428 affine_constraints.add_entries_local_to_global(
429 dof_indices, dsp, keep_constrained);
442 template <
int dim,
typename Number,
typename SPARSITY>
444 const dealii::DoFHandler<dim> &dof_handler,
446 const dealii::AffineConstraints<Number> &affine_constraints,
447 bool keep_constrained)
449 Assert(affine_constraints.n_constraints() == 0,
450 dealii::ExcMessage(
"I don't think constraints make sense for dG"));
452 std::vector<dealii::types::global_dof_index> dof_indices;
453 std::vector<dealii::types::global_dof_index> neighbor_dof_indices;
459 std::vector<dealii::types::global_dof_index> coupling_indices;
460 std::vector<dealii::types::global_dof_index> neighbor_coupling_indices;
463 for (
auto cell : dof_handler.active_cell_iterators()) {
464 if (cell->is_artificial())
467 const unsigned int dofs_per_cell = cell->get_fe().n_dofs_per_cell();
468 dof_indices.resize(dofs_per_cell);
469 cell->get_dof_indices(dof_indices);
471 affine_constraints.add_entries_local_to_global(
472 dof_indices, dsp, keep_constrained);
474 for (
const auto f_index : cell->face_indices()) {
475 const auto &face = cell->face(f_index);
478 const bool has_neighbor =
479 !face->at_boundary() || cell->has_periodic_neighbor(f_index);
484 const auto neighbor_cell =
485 cell->neighbor_or_periodic_neighbor(f_index);
486 if (neighbor_cell->is_artificial())
489 const unsigned int neighbor_dofs_per_cell =
490 neighbor_cell->get_fe().n_dofs_per_cell();
491 neighbor_dof_indices.resize(neighbor_dofs_per_cell);
492 neighbor_cell->get_dof_indices(neighbor_dof_indices);
494 const unsigned int f_index_neighbor =
495 cell->has_periodic_neighbor(f_index)
496 ? cell->periodic_neighbor_of_periodic_neighbor(f_index)
497 : cell->neighbor_of_neighbor(f_index);
504 coupling_indices.resize(0);
505 for (
unsigned int i = 0; i < dofs_per_cell; ++i)
506 if (cell->get_fe().has_support_on_face(i, f_index))
507 coupling_indices.push_back(dof_indices[i]);
509 neighbor_coupling_indices.resize(0);
510 for (
unsigned int j = 0; j < neighbor_dofs_per_cell; ++j)
511 if (neighbor_cell->get_fe().has_support_on_face(j,
513 neighbor_coupling_indices.push_back(neighbor_dof_indices[j]);
515 affine_constraints.add_entries_local_to_global(
517 neighbor_coupling_indices,
void make_extended_sparsity_pattern(const dealii::DoFHandler< dim > &dof_handler, SPARSITY &dsp, const dealii::AffineConstraints< Number > &affine_constraints, bool keep_constrained)
void transform_to_local_range(const dealii::Utilities::MPI::Partitioner &partitioner, dealii::AffineConstraints< Number > &affine_constraints)
unsigned int internal_range(dealii::DoFHandler< dim > &dof_handler, const dealii::DynamicSparsityPattern &sparsity, const std::size_t group_size)
unsigned int export_indices_first(dealii::DoFHandler< dim > &dof_handler, const MPI_Comm &mpi_communicator, const unsigned int n_locally_internal, const std::size_t group_size)
unsigned int inconsistent_strides_last(dealii::DoFHandler< dim > &dof_handler, const dealii::DynamicSparsityPattern &sparsity, const unsigned int n_locally_internal, const std::size_t group_size)
void make_extended_sparsity_pattern_dg(const dealii::DoFHandler< dim > &dof_handler, SPARSITY &dsp, const dealii::AffineConstraints< Number > &affine_constraints, bool keep_constrained)