ryujin 2.1.1 revision 350e54cc11f3d21282bddcf3e6517944dbc508bf
offline_data.template.h
Go to the documentation of this file.
1//
2// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
3// Copyright (C) 2020 - 2024 by the ryujin authors
4//
5
6#pragma once
7
8#include "discretization.h"
11#include "offline_data.h"
12#include "scratch_data.h"
13#include "sparse_matrix_simd.template.h" /* instantiate read_in */
14
15#include <deal.II/base/graph_coloring.h>
16#include <deal.II/base/parallel.h>
17#include <deal.II/base/work_stream.h>
18#include <deal.II/dofs/dof_renumbering.h>
19#include <deal.II/dofs/dof_tools.h>
20#include <deal.II/fe/fe_values.h>
21#include <deal.II/grid/grid_tools.h>
22#include <deal.II/lac/dynamic_sparsity_pattern.h>
23#include <deal.II/lac/la_parallel_vector.h>
24#ifdef DEAL_II_WITH_TRILINOS
25#include <deal.II/lac/trilinos_sparse_matrix.h>
26#endif
27
28#ifdef FORCE_DEAL_II_SPARSE_MATRIX
29#undef DEAL_II_WITH_TRILINOS
30#endif
31
32namespace ryujin
33{
34 using namespace dealii;
35
36
37 template <int dim, typename Number>
39 const MPIEnsemble &mpi_ensemble,
40 const Discretization<dim> &discretization,
41 const std::string &subsection /*= "OfflineData"*/)
42 : ParameterAcceptor(subsection)
43 , mpi_ensemble_(mpi_ensemble)
44 , discretization_(&discretization)
45 {
46 incidence_relaxation_even_ = 0.5;
47 add_parameter("incidence matrix relaxation even degree",
48 incidence_relaxation_even_,
49 "Scaling exponent for incidence matrix used for "
50 "discontinuous finite elements with even degree. The default "
51 "value 0.5 scales the jump penalization with (h_i+h_j)^0.5.");
52
53 incidence_relaxation_odd_ = 0.0;
54 add_parameter("incidence matrix relaxation odd degree",
55 incidence_relaxation_odd_,
56 "Scaling exponent for incidence matrix used for "
57 "discontinuous finite elements with even degree. The default "
58 "value of 0.0 sets the jump penalization to a constant 1.");
59 }
60
61
62 template <int dim, typename Number>
64 {
65 /*
66 * First, we set up the locally_relevant index set, determine (globally
67 * indexed) affine constraints and create a (globally indexed) sparsity
68 * pattern:
69 */
70
71 auto &dof_handler = *dof_handler_;
72 const IndexSet &locally_owned = dof_handler.locally_owned_dofs();
73
74 IndexSet locally_relevant;
75 DoFTools::extract_locally_relevant_dofs(dof_handler, locally_relevant);
76
77 affine_constraints_.reinit(locally_relevant);
78 DoFTools::make_hanging_node_constraints(dof_handler, affine_constraints_);
79
80#ifndef DEAL_II_WITH_TRILINOS
81 AssertThrow(affine_constraints_.n_constraints() == 0,
82 ExcMessage("ryujin was built without Trilinos support - no "
83 "hanging node support available"));
84#endif
85
86 /*
87 * Enforce periodic boundary conditions. We assume that the mesh is in
88 * "normal configuration."
89 */
90
91 const auto &periodic_faces =
92 discretization_->triangulation().get_periodic_face_map();
93
94 for (const auto &[left, value] : periodic_faces) {
95 const auto &[right, orientation] = value;
96
97 typename DoFHandler<dim>::cell_iterator dof_cell_left(
98 &left.first->get_triangulation(),
99 left.first->level(),
100 left.first->index(),
101 &dof_handler);
102
103 typename DoFHandler<dim>::cell_iterator dof_cell_right(
104 &right.first->get_triangulation(),
105 right.first->level(),
106 right.first->index(),
107 &dof_handler);
108
109 if constexpr (std::is_same_v<Number, double>) {
110 DoFTools::make_periodicity_constraints(
111 dof_cell_left->face(left.second),
112 dof_cell_right->face(right.second),
113 affine_constraints_,
114 ComponentMask(),
115#if DEAL_II_VERSION_GTE(9, 6, 0)
116 orientation);
117#else
118 /* orientation */ orientation[0],
119 /* flip */ orientation[1],
120 /* rotation */ orientation[2]);
121#endif
122 } else {
123 AssertThrow(false, dealii::ExcNotImplemented());
124 __builtin_trap();
125 }
126 }
127
128 affine_constraints_.close();
129
130#ifdef DEBUG
131 {
132 /* Check that constraints are consistent in parallel: */
133 const std::vector<IndexSet> &locally_owned_dofs =
134 Utilities::MPI::all_gather(mpi_ensemble_.ensemble_communicator(),
135 dof_handler.locally_owned_dofs());
136 const IndexSet locally_active =
137 dealii::DoFTools::extract_locally_active_dofs(dof_handler);
138 Assert(affine_constraints_.is_consistent_in_parallel(
139 locally_owned_dofs,
140 locally_active,
141 mpi_ensemble_.ensemble_communicator(),
142 /*verbose*/ true),
143 ExcInternalError());
144 }
145#endif
146
147 sparsity_pattern_.reinit(
148 dof_handler.n_dofs(), dof_handler.n_dofs(), locally_relevant);
149
150 if (discretization_->have_discontinuous_ansatz()) {
151 /*
152 * Create dG sparsity pattern:
153 */
155 dof_handler, sparsity_pattern_, affine_constraints_, false);
156 } else {
157 /*
158 * Create cG sparsity pattern:
159 */
160#ifdef DEAL_II_WITH_TRILINOS
161 DoFTools::make_sparsity_pattern(
162 dof_handler, sparsity_pattern_, affine_constraints_, false);
163#else
164 /*
165 * In case we use dealii::SparseMatrix<Number> for assembly we need a
166 * sparsity pattern that also includes the full locally relevant -
167 * locally relevant coupling block. This gets thrown out again later,
168 * but nevertheless we have to add it.
169 */
171 dof_handler, sparsity_pattern_, affine_constraints_, false);
172#endif
173 }
174
175 /*
176 * We have to complete the local stencil to have consistent size over
177 * all MPI ranks. Otherwise, MPI synchronization in our
178 * SparseMatrixSIMD class will fail.
179 */
180
181 SparsityTools::distribute_sparsity_pattern(
182 sparsity_pattern_,
183 locally_owned,
184 mpi_ensemble_.ensemble_communicator(),
185 locally_relevant);
186 }
187
188
189 template <int dim, typename Number>
190 void OfflineData<dim, Number>::setup(const unsigned int problem_dimension,
191 const unsigned int n_precomputed_values)
192 {
193#ifdef DEBUG_OUTPUT
194 std::cout << "OfflineData<dim, Number>::setup()" << std::endl;
195#endif
196
197 /*
198 * Initialize dof handler:
199 */
200
201 const auto &triangulation = discretization_->triangulation();
202 if (!dof_handler_)
203 dof_handler_ = std::make_unique<dealii::DoFHandler<dim>>(triangulation);
204 auto &dof_handler = *dof_handler_;
205
206 dof_handler.distribute_dofs(discretization_->finite_element());
207
208 n_locally_owned_ = dof_handler.locally_owned_dofs().n_elements();
209
210 /*
211 * Renumbering:
212 */
213
214 DoFRenumbering::Cuthill_McKee(dof_handler);
215
216 /*
217 * Reorder all (individual) export indices at the beginning of the
218 * locally_internal index range to achieve a better packing:
219 *
220 * Note: This function might miss export indices that come from
221 * eliminating hanging node and periodicity constraints (which we do
222 * not know at this point because they depend on the renumbering...).
223 */
225 mpi_ensemble_.ensemble_communicator(),
226 n_locally_owned_,
227 1);
228
229 /*
230 * Group degrees of freedom that have the same stencil size in groups
231 * of multiples of the VectorizedArray<Number>::size().
232 *
233 * In order to determine the stencil size we have to create a first,
234 * temporary sparsity pattern:
235 */
236 create_constraints_and_sparsity_pattern();
237 n_locally_internal_ = DoFRenumbering::internal_range(
238 dof_handler, sparsity_pattern_, VectorizedArray<Number>::size());
239
240 /*
241 * Reorder all (strides of) locally internal indices that contain
242 * export indices to the start of the index range. This reordering
243 * preserves the binning introduced by
244 * DoFRenumbering::internal_range().
245 *
246 * Note: This function might miss export indices that come from
247 * eliminating hanging node and periodicity constraints (which we do
248 * not know at this point because they depend on the renumbering...).
249 * We therefore have to update n_export_indices_ later again.
250 */
251 n_export_indices_ = DoFRenumbering::export_indices_first(
252 dof_handler,
253 mpi_ensemble_.ensemble_communicator(),
254 n_locally_internal_,
255 VectorizedArray<Number>::size());
256
257 /*
258 * A small lambda to check for stride-level consistency of the internal
259 * index range:
260 */
261 const auto consistent_stride_range [[maybe_unused]] = [&]() {
262 constexpr auto group_size = VectorizedArray<Number>::size();
263 const IndexSet &locally_owned = dof_handler.locally_owned_dofs();
264 const auto offset = n_locally_owned_ != 0 ? *locally_owned.begin() : 0;
265
266 unsigned int group_row_length = 0;
267 unsigned int i = 0;
268 for (; i < n_locally_internal_; ++i) {
269 if (i % group_size == 0) {
270 group_row_length = sparsity_pattern_.row_length(offset + i);
271 } else {
272 if (group_row_length != sparsity_pattern_.row_length(offset + i)) {
273 break;
274 }
275 }
276 }
277 return i / group_size * group_size;
278 };
279
280 /*
281 * A small lambda that performs a "logical or" over all MPI ranks:
282 */
283 const auto mpi_allreduce_logical_or = [&](const bool local_value) {
284 std::function<bool(const bool &, const bool &)> comparator =
285 [](const bool &left, const bool &right) -> bool {
286 return left || right;
287 };
288 return Utilities::MPI::all_reduce(
289 local_value, mpi_ensemble_.ensemble_communicator(), comparator);
290 };
291
292 /*
293 * Create final sparsity pattern:
294 */
295
296 create_constraints_and_sparsity_pattern();
297
298 /*
299 * We have to ensure that the locally internal numbering range is still
300 * consistent, meaning that all strides have the same stencil size.
301 * This property might not hold any more after the elimination
302 * procedure of constrained degrees of freedom (periodicity, or hanging
303 * node constraints). Therefore, the following little dance:
304 */
305
306#if DEAL_II_VERSION_GTE(9, 5, 0)
307 if (mpi_allreduce_logical_or(affine_constraints_.n_constraints() > 0)) {
308 if (mpi_allreduce_logical_or( //
309 consistent_stride_range() != n_locally_internal_)) {
310 /*
311 * In this case we try to fix up the numbering by pushing affected
312 * strides to the end and slightly lowering the n_locally_internal_
313 * marker.
314 */
315 n_locally_internal_ = DoFRenumbering::inconsistent_strides_last(
316 dof_handler,
317 sparsity_pattern_,
318 n_locally_internal_,
319 VectorizedArray<Number>::size());
320 create_constraints_and_sparsity_pattern();
321 n_locally_internal_ = consistent_stride_range();
322 }
323 }
324#endif
325
326 /*
327 * Check that after all the dof manipulation and setup we still end up
328 * with indices in [0, locally_internal) that have uniform stencil size
329 * within a stride.
330 */
331 Assert(consistent_stride_range() == n_locally_internal_,
332 dealii::ExcInternalError());
333
334 /*
335 * Set up partitioner:
336 */
337
338 const IndexSet &locally_owned = dof_handler.locally_owned_dofs();
339 Assert(n_locally_owned_ == locally_owned.n_elements(),
340 dealii::ExcInternalError());
341
342 IndexSet locally_relevant;
343 DoFTools::extract_locally_relevant_dofs(dof_handler, locally_relevant);
344 /* Enlarge the locally relevant set to include all additional couplings: */
345 {
346 IndexSet additional_dofs(dof_handler.n_dofs());
347 for (auto &entry : sparsity_pattern_)
348 if (!locally_relevant.is_element(entry.column())) {
349 Assert(locally_owned.is_element(entry.row()), ExcInternalError());
350 additional_dofs.add_index(entry.column());
351 }
352 additional_dofs.compress();
353 locally_relevant.add_indices(additional_dofs);
354 locally_relevant.compress();
355 }
356
357 n_locally_relevant_ = locally_relevant.n_elements();
358
359 scalar_partitioner_ = std::make_shared<dealii::Utilities::MPI::Partitioner>(
360 locally_owned, locally_relevant, mpi_ensemble_.ensemble_communicator());
361
362 hyperbolic_vector_partitioner_ = Vectors::create_vector_partitioner(
363 scalar_partitioner_, problem_dimension);
364
365 precomputed_vector_partitioner_ = Vectors::create_vector_partitioner(
366 scalar_partitioner_, n_precomputed_values);
367
368 /*
369 * After eliminiating periodicity and hanging node constraints we need
370 * to update n_export_indices_ again. This happens because we need to
371 * call export_indices_first() with incomplete information (missing
372 * eliminated degrees of freedom).
373 */
374 if (mpi_allreduce_logical_or(affine_constraints_.n_constraints() > 0)) {
375 /*
376 * Recalculate n_export_indices_:
377 */
378 n_export_indices_ = 0;
379 for (const auto &it : scalar_partitioner_->import_indices())
380 if (it.second <= n_locally_internal_)
381 n_export_indices_ = std::max(n_export_indices_, it.second);
382
383 constexpr auto simd_length = VectorizedArray<Number>::size();
384 n_export_indices_ =
385 (n_export_indices_ + simd_length - 1) / simd_length * simd_length;
386 }
387
388#ifdef DEBUG
389 /* Check that n_export_indices_ is valid: */
390 unsigned int control = 0;
391 for (const auto &it : scalar_partitioner_->import_indices())
392 if (it.second <= n_locally_internal_)
393 control = std::max(control, it.second);
394
395 Assert(control <= n_export_indices_, ExcInternalError());
396 Assert(n_export_indices_ <= n_locally_internal_, ExcInternalError());
397#endif
398
399 /*
400 * Set up SIMD sparsity pattern in local numbering. Nota bene: The
401 * SparsityPatternSIMD::reinit() function will translates the pattern
402 * from global deal.II (typical) dof indexing to local indices.
403 */
404
405 sparsity_pattern_simd_.reinit(
406 n_locally_internal_, sparsity_pattern_, scalar_partitioner_);
407 }
408
409
410 template <int dim, typename Number>
411 void OfflineData<dim, Number>::create_matrices()
412 {
413#ifdef DEBUG_OUTPUT
414 std::cout << "OfflineData<dim, Number>::create_matrices()" << std::endl;
415#endif
416
417 /*
418 * First, (re)initialize all local matrices:
419 */
420
421 mass_matrix_.reinit(sparsity_pattern_simd_);
422 if (discretization_->have_discontinuous_ansatz())
423 mass_matrix_inverse_.reinit(sparsity_pattern_simd_);
424
425 lumped_mass_matrix_.reinit(scalar_partitioner_);
426 lumped_mass_matrix_inverse_.reinit(scalar_partitioner_);
427
428 cij_matrix_.reinit(sparsity_pattern_simd_);
429 if (discretization_->have_discontinuous_ansatz())
430 incidence_matrix_.reinit(sparsity_pattern_simd_);
431
432 /*
433 * Then, assemble:
434 */
435
436 auto &dof_handler = *dof_handler_;
437
438 measure_of_omega_ = 0.;
439
440#ifdef DEAL_II_WITH_TRILINOS
441 /* Variant using TrilinosWrappers::SparseMatrix with global numbering */
442
443 AffineConstraints<double> affine_constraints_assembly;
444 /* This small dance is necessary to translate from Number to double: */
445 affine_constraints_assembly.reinit(affine_constraints_.get_local_lines());
446 for (auto line : affine_constraints_.get_lines()) {
447 affine_constraints_assembly.add_line(line.index);
448 for (auto entry : line.entries)
449 affine_constraints_assembly.add_entry(
450 line.index, entry.first, entry.second);
451 affine_constraints_assembly.set_inhomogeneity(line.index,
452 line.inhomogeneity);
453 }
454 affine_constraints_assembly.close();
455
456 const IndexSet &locally_owned = dof_handler.locally_owned_dofs();
457 TrilinosWrappers::SparsityPattern trilinos_sparsity_pattern;
458 trilinos_sparsity_pattern.reinit(locally_owned,
459 sparsity_pattern_,
460 mpi_ensemble_.ensemble_communicator());
461
462 TrilinosWrappers::SparseMatrix mass_matrix_tmp;
463 TrilinosWrappers::SparseMatrix mass_matrix_inverse_tmp;
464 if (discretization_->have_discontinuous_ansatz())
465 mass_matrix_inverse_tmp.reinit(trilinos_sparsity_pattern);
466
467 std::array<TrilinosWrappers::SparseMatrix, dim> cij_matrix_tmp;
468
469 mass_matrix_tmp.reinit(trilinos_sparsity_pattern);
470 for (auto &matrix : cij_matrix_tmp)
471 matrix.reinit(trilinos_sparsity_pattern);
472
473#else
474 /* Variant using deal.II SparseMatrix with local numbering */
475
476 AffineConstraints<Number> affine_constraints_assembly;
477 affine_constraints_assembly.copy_from(affine_constraints_);
478 transform_to_local_range(*scalar_partitioner_, affine_constraints_assembly);
479
480 SparsityPattern sparsity_pattern_assembly;
481 {
482 DynamicSparsityPattern dsp(n_locally_relevant_, n_locally_relevant_);
483 for (const auto &entry : sparsity_pattern_) {
484 const auto i = scalar_partitioner_->global_to_local(entry.row());
485 const auto j = scalar_partitioner_->global_to_local(entry.column());
486 dsp.add(i, j);
487 }
488 sparsity_pattern_assembly.copy_from(dsp);
489 }
490
491 dealii::SparseMatrix<Number> mass_matrix_tmp;
492 dealii::SparseMatrix<Number> mass_matrix_inverse_tmp;
493 if (discretization_->have_discontinuous_ansatz())
494 mass_matrix_inverse_tmp.reinit(sparsity_pattern_assembly);
495
496 std::array<dealii::SparseMatrix<Number>, dim> cij_matrix_tmp;
497
498 mass_matrix_tmp.reinit(sparsity_pattern_assembly);
499 for (auto &matrix : cij_matrix_tmp)
500 matrix.reinit(sparsity_pattern_assembly);
501#endif
502
503 /*
504 * Now, assemble all matrices:
505 */
506
507 /* The local, per-cell assembly routine: */
508 const auto local_assemble_system = [&](const auto &cell,
509 auto &scratch,
510 auto &copy) {
511 /* iterate over locally owned cells and the ghost layer */
512
513 auto &is_locally_owned = copy.is_locally_owned_;
514 auto &local_dof_indices = copy.local_dof_indices_;
515 auto &neighbor_local_dof_indices = copy.neighbor_local_dof_indices_;
516
517 auto &cell_mass_matrix = copy.cell_mass_matrix_;
518 auto &cell_mass_matrix_inverse = copy.cell_mass_matrix_inverse_;
519 auto &cell_cij_matrix = copy.cell_cij_matrix_;
520 auto &interface_cij_matrix = copy.interface_cij_matrix_;
521 auto &cell_measure = copy.cell_measure_;
522
523 auto &hp_fe_values = scratch.hp_fe_values_;
524 auto &hp_fe_face_values = scratch.hp_fe_face_values_;
525 auto &hp_fe_neighbor_face_values = scratch.hp_fe_neighbor_face_values_;
526
527#ifdef DEAL_II_WITH_TRILINOS
528 is_locally_owned = cell->is_locally_owned();
529#else
530 /*
531 * When using a local dealii::SparseMatrix<Number> we don not
532 * have a compress(VectorOperation::add) available. In this case
533 * we assemble contributions over all locally relevant (non
534 * artificial) cells.
535 */
536 is_locally_owned = !cell->is_artificial();
537#endif
538 if (!is_locally_owned)
539 return;
540
541 const unsigned int dofs_per_cell = cell->get_fe().n_dofs_per_cell();
542
543 cell_mass_matrix.reinit(dofs_per_cell, dofs_per_cell);
544 for (auto &matrix : cell_cij_matrix)
545 matrix.reinit(dofs_per_cell, dofs_per_cell);
546 if (discretization_->have_discontinuous_ansatz()) {
547 cell_mass_matrix_inverse.reinit(dofs_per_cell, dofs_per_cell);
548 for (auto &it : interface_cij_matrix)
549 for (auto &matrix : it)
550 matrix.reinit(dofs_per_cell, dofs_per_cell);
551 }
552
553 hp_fe_values.reinit(cell);
554 const auto &fe_values = hp_fe_values.get_present_fe_values();
555
556 local_dof_indices.resize(dofs_per_cell);
557 cell->get_dof_indices(local_dof_indices);
558
559 /* clear out copy data: */
560
561 cell_mass_matrix = 0.;
562 for (auto &matrix : cell_cij_matrix)
563 matrix = 0.;
564 if (discretization_->have_discontinuous_ansatz()) {
565 cell_mass_matrix_inverse = 0.;
566 for (auto &it : interface_cij_matrix)
567 for (auto &matrix : it)
568 matrix = 0.;
569 }
570 cell_measure = 0.;
571
572 for (unsigned int q : fe_values.quadrature_point_indices()) {
573 const auto JxW = fe_values.JxW(q);
574
575 if (cell->is_locally_owned())
576 cell_measure += Number(JxW);
577
578 for (unsigned int j : fe_values.dof_indices()) {
579 const auto value_JxW = fe_values.shape_value(j, q) * JxW;
580 const auto grad_JxW = fe_values.shape_grad(j, q) * JxW;
581
582 for (unsigned int i : fe_values.dof_indices()) {
583 const auto value = fe_values.shape_value(i, q);
584
585 cell_mass_matrix(i, j) += Number(value * value_JxW);
586 for (unsigned int d = 0; d < dim; ++d)
587 cell_cij_matrix[d](i, j) += Number((value * grad_JxW)[d]);
588 } /* for i */
589 } /* for j */
590 } /* for q */
591
592 /*
593 * For a discontinuous finite element ansatz we need to assemble
594 * additional face contributions:
595 */
596
597 if (!discretization_->have_discontinuous_ansatz())
598 return;
599
600 for (const auto f_index : cell->face_indices()) {
601 const auto &face = cell->face(f_index);
602
603 /* Skip faces without neighbors... */
604 const bool has_neighbor =
605 !face->at_boundary() || cell->has_periodic_neighbor(f_index);
606 if (!has_neighbor) {
607 // set the vector of local dof indices to 0 to indicate that
608 // there is nothing to do for this face:
609 neighbor_local_dof_indices[f_index].resize(0);
610 continue;
611 }
612
613 /* Avoid artificial cells: */
614 const auto neighbor_cell = cell->neighbor_or_periodic_neighbor(f_index);
615 if (neighbor_cell->is_artificial()) {
616 // set the vector of local dof indices to 0 to indicate that
617 // there is nothing to do for this face:
618 neighbor_local_dof_indices[f_index].resize(0);
619 continue;
620 }
621
622 hp_fe_face_values.reinit(cell, f_index);
623 const auto &fe_face_values = hp_fe_face_values.get_present_fe_values();
624
625 /* Face contribution: */
626
627 for (unsigned int q : fe_face_values.quadrature_point_indices()) {
628 const auto JxW = fe_face_values.JxW(q);
629 const auto &normal = fe_face_values.get_normal_vectors()[q];
630
631 for (unsigned int j : fe_face_values.dof_indices()) {
632 const auto value_JxW = fe_face_values.shape_value(j, q) * JxW;
633
634 for (unsigned int i : fe_face_values.dof_indices()) {
635 const auto value = fe_face_values.shape_value(i, q);
636
637 for (unsigned int d = 0; d < dim; ++d)
638 cell_cij_matrix[d](i, j) -=
639 Number(0.5 * normal[d] * value * value_JxW);
640 } /* for i */
641 } /* for j */
642 } /* for q */
643
644 /* Coupling part: */
645
646 const unsigned int f_index_neighbor =
647 cell->has_periodic_neighbor(f_index)
648 ? cell->periodic_neighbor_of_periodic_neighbor(f_index)
649 : cell->neighbor_of_neighbor(f_index);
650
651 neighbor_local_dof_indices[f_index].resize(dofs_per_cell);
652 neighbor_cell->get_dof_indices(neighbor_local_dof_indices[f_index]);
653
654 hp_fe_neighbor_face_values.reinit(neighbor_cell, f_index_neighbor);
655 const auto &fe_neighbor_face_values =
656 hp_fe_neighbor_face_values.get_present_fe_values();
657
658 for (unsigned int q : fe_face_values.quadrature_point_indices()) {
659 const auto JxW = fe_face_values.JxW(q);
660 const auto &normal = fe_face_values.get_normal_vectors()[q];
661
662 /* index j for neighbor, index i for current cell: */
663 for (unsigned int j : fe_face_values.dof_indices()) {
664 const auto value_JxW =
665 fe_neighbor_face_values.shape_value(j, q) * JxW;
666
667 for (unsigned int i : fe_face_values.dof_indices()) {
668 const auto value = fe_face_values.shape_value(i, q);
669
670 for (unsigned int d = 0; d < dim; ++d)
671 interface_cij_matrix[f_index][d](i, j) +=
672 Number(0.5 * normal[d] * value * value_JxW);
673 } /* for i */
674 } /* for j */
675 } /* for q */
676 }
677
678 /*
679 * Compute block inverse of mass matrix:
680 */
681
682 if (discretization_->have_discontinuous_ansatz()) {
683 // FIXME: rewrite with CellwiseInverseMassMatrix
684 cell_mass_matrix_inverse.invert(cell_mass_matrix);
685 }
686 };
687
688 const auto copy_local_to_global = [&](const auto &copy) {
689 const auto &is_locally_owned = copy.is_locally_owned_;
690#ifdef DEAL_II_WITH_TRILINOS
691 const auto &local_dof_indices = copy.local_dof_indices_;
692 const auto &neighbor_local_dof_indices = copy.neighbor_local_dof_indices_;
693#else
694 /*
695 * We have to transform indices to the local index range
696 * [0, n_locally_relevant_) when using the dealii::SparseMatrix.
697 * Thus, copy all index vectors:
698 */
699 auto local_dof_indices = copy.local_dof_indices_;
700 auto neighbor_local_dof_indices = copy.neighbor_local_dof_indices_;
701#endif
702 const auto &cell_mass_matrix = copy.cell_mass_matrix_;
703 const auto &cell_mass_matrix_inverse = copy.cell_mass_matrix_inverse_;
704 const auto &cell_cij_matrix = copy.cell_cij_matrix_;
705 const auto &interface_cij_matrix = copy.interface_cij_matrix_;
706 const auto &cell_measure = copy.cell_measure_;
707
708 if (!is_locally_owned)
709 return;
710
711#ifndef DEAL_II_WITH_TRILINOS
712 transform_to_local_range(*scalar_partitioner_, local_dof_indices);
713 for (auto &indices : neighbor_local_dof_indices)
714 transform_to_local_range(*scalar_partitioner_, indices);
715#endif
716
717 affine_constraints_assembly.distribute_local_to_global(
718 cell_mass_matrix, local_dof_indices, mass_matrix_tmp);
719
720 for (int k = 0; k < dim; ++k) {
721 affine_constraints_assembly.distribute_local_to_global(
722 cell_cij_matrix[k], local_dof_indices, cij_matrix_tmp[k]);
723
724 for (unsigned int f_index = 0; f_index < copy.n_faces; ++f_index) {
725 if (neighbor_local_dof_indices[f_index].size() != 0) {
726 affine_constraints_assembly.distribute_local_to_global(
727 interface_cij_matrix[f_index][k],
728 local_dof_indices,
729 neighbor_local_dof_indices[f_index],
730 cij_matrix_tmp[k]);
731 }
732 }
733 }
734
735 if (discretization_->have_discontinuous_ansatz())
736 affine_constraints_assembly.distribute_local_to_global(
737 cell_mass_matrix_inverse,
738 local_dof_indices,
739 mass_matrix_inverse_tmp);
740
741 measure_of_omega_ += cell_measure;
742 };
743
744 WorkStream::run(dof_handler.begin_active(),
745 dof_handler.end(),
746 local_assemble_system,
747 copy_local_to_global,
748 AssemblyScratchData<dim>(*discretization_),
749#ifdef DEAL_II_WITH_TRILINOS
750 AssemblyCopyData<dim, double>());
751#else
752 AssemblyCopyData<dim, Number>());
753#endif
754
755#ifdef DEAL_II_WITH_TRILINOS
756 mass_matrix_tmp.compress(VectorOperation::add);
757 for (auto &it : cij_matrix_tmp)
758 it.compress(VectorOperation::add);
759
760 mass_matrix_.read_in(mass_matrix_tmp, /*locally_indexed*/ false);
761 if (discretization_->have_discontinuous_ansatz())
762 mass_matrix_inverse_.read_in(mass_matrix_inverse_tmp, /*loc_ind*/ false);
763 cij_matrix_.read_in(cij_matrix_tmp, /*locally_indexed*/ false);
764#else
765 mass_matrix_.read_in(mass_matrix_tmp, /*locally_indexed*/ true);
766 if (discretization_->have_discontinuous_ansatz())
767 mass_matrix_inverse_.read_in(mass_matrix_inverse_tmp, /*loc_ind*/ true);
768 cij_matrix_.read_in(cij_matrix_tmp, /*locally_indexed*/ true);
769#endif
770
771 mass_matrix_.update_ghost_rows();
772 if (discretization_->have_discontinuous_ansatz())
773 mass_matrix_inverse_.update_ghost_rows();
774 cij_matrix_.update_ghost_rows();
775
776 measure_of_omega_ = Utilities::MPI::sum(
777 measure_of_omega_, mpi_ensemble_.ensemble_communicator());
778
779 /*
780 * Create lumped mass matrix:
781 */
782
783 {
784#ifdef DEAL_II_WITH_TRILINOS
785 ScalarVector one(scalar_partitioner_);
786 one = 1.;
787 affine_constraints_assembly.set_zero(one);
788
789 ScalarVector local_lumped_mass_matrix(scalar_partitioner_);
790 mass_matrix_tmp.vmult(local_lumped_mass_matrix, one);
791 local_lumped_mass_matrix.compress(VectorOperation::add);
792
793 for (unsigned int i = 0; i < scalar_partitioner_->locally_owned_size();
794 ++i) {
795 lumped_mass_matrix_.local_element(i) =
796 local_lumped_mass_matrix.local_element(i);
797 lumped_mass_matrix_inverse_.local_element(i) =
798 1. / lumped_mass_matrix_.local_element(i);
799 }
800 lumped_mass_matrix_.update_ghost_values();
801 lumped_mass_matrix_inverse_.update_ghost_values();
802
803#else
804
805 dealii::Vector<Number> one(mass_matrix_tmp.m());
806 one = 1.;
807 affine_constraints_assembly.set_zero(one);
808
809 dealii::Vector<Number> local_lumped_mass_matrix(mass_matrix_tmp.m());
810 mass_matrix_tmp.vmult(local_lumped_mass_matrix, one);
811
812 for (unsigned int i = 0; i < scalar_partitioner_->locally_owned_size();
813 ++i) {
814 lumped_mass_matrix_.local_element(i) = local_lumped_mass_matrix(i);
815 lumped_mass_matrix_inverse_.local_element(i) =
816 1. / lumped_mass_matrix_.local_element(i);
817 }
818 lumped_mass_matrix_.update_ghost_values();
819 lumped_mass_matrix_inverse_.update_ghost_values();
820#endif
821 }
822
823 /*
824 * Assemble incidence matrix:
825 */
826
827 if (discretization_->have_discontinuous_ansatz()) {
828#ifdef DEAL_II_WITH_TRILINOS
829 TrilinosWrappers::SparseMatrix incidence_matrix_tmp;
830 incidence_matrix_tmp.reinit(trilinos_sparsity_pattern);
831#else
832 dealii::SparseMatrix<Number> incidence_matrix_tmp;
833 incidence_matrix_tmp.reinit(sparsity_pattern_assembly);
834#endif
835
836 /* The local, per-cell assembly routine: */
837 const auto local_assemble_system = [&](const auto &cell,
838 auto &scratch,
839 auto &copy) {
840 /* iterate over locally owned cells and the ghost layer */
841
842 auto &is_locally_owned = copy.is_locally_owned_;
843 auto &local_dof_indices = copy.local_dof_indices_;
844 auto &neighbor_local_dof_indices = copy.neighbor_local_dof_indices_;
845 auto &interface_incidence_matrix = copy.interface_incidence_matrix_;
846 auto &hp_fe_face_values_nodal = scratch.hp_fe_face_values_nodal_;
847 auto &hp_fe_neighbor_face_values_nodal =
848 scratch.hp_fe_neighbor_face_values_nodal_;
849
850#ifdef DEAL_II_WITH_TRILINOS
851 is_locally_owned = cell->is_locally_owned();
852#else
853 is_locally_owned = !cell->is_artificial();
854#endif
855 if (!is_locally_owned)
856 return;
857
858 const unsigned int dofs_per_cell = cell->get_fe().n_dofs_per_cell();
859
860 for (auto &matrix : interface_incidence_matrix)
861 matrix.reinit(dofs_per_cell, dofs_per_cell);
862
863 local_dof_indices.resize(dofs_per_cell);
864 cell->get_dof_indices(local_dof_indices);
865
866 /* clear out copy data: */
867 for (auto &matrix : interface_incidence_matrix)
868 matrix = 0.;
869
870 for (const auto f_index : cell->face_indices()) {
871 const auto &face = cell->face(f_index);
872
873 /* Skip faces without neighbors... */
874 const bool has_neighbor =
875 !face->at_boundary() || cell->has_periodic_neighbor(f_index);
876 if (!has_neighbor) {
877 // set the vector of local dof indices to 0 to indicate that
878 // there is nothing to do for this face:
879 neighbor_local_dof_indices[f_index].resize(0);
880 continue;
881 }
882
883 /* Avoid artificial cells: */
884 const auto neighbor_cell =
885 cell->neighbor_or_periodic_neighbor(f_index);
886 if (neighbor_cell->is_artificial()) {
887 // set the vector of local dof indices to 0 to indicate that
888 // there is nothing to do for this face:
889 neighbor_local_dof_indices[f_index].resize(0);
890 continue;
891 }
892
893 const unsigned int f_index_neighbor =
894 cell->has_periodic_neighbor(f_index)
895 ? cell->periodic_neighbor_of_periodic_neighbor(f_index)
896 : cell->neighbor_of_neighbor(f_index);
897
898 neighbor_local_dof_indices[f_index].resize(dofs_per_cell);
899 neighbor_cell->get_dof_indices(neighbor_local_dof_indices[f_index]);
900
901 hp_fe_face_values_nodal.reinit(cell, f_index);
902 const auto &fe_face_values_nodal =
903 hp_fe_face_values_nodal.get_present_fe_values();
904 hp_fe_neighbor_face_values_nodal.reinit(neighbor_cell,
905 f_index_neighbor);
906 const auto &fe_neighbor_face_values_nodal =
907 hp_fe_neighbor_face_values_nodal.get_present_fe_values();
908
909 /* Lumped incidence matrix: */
910
911 for (unsigned int q :
912 fe_face_values_nodal.quadrature_point_indices()) {
913 /* index j for neighbor, index i for current cell: */
914 for (unsigned int j : fe_neighbor_face_values_nodal.dof_indices()) {
915 const auto v_j = fe_neighbor_face_values_nodal.shape_value(j, q);
916 for (unsigned int i : fe_face_values_nodal.dof_indices()) {
917 const auto v_i = fe_face_values_nodal.shape_value(i, q);
918 constexpr auto eps = std::numeric_limits<Number>::epsilon();
919 if (std::abs(v_i * v_j) > 100. * eps) {
920 const auto &ansatz = discretization_->ansatz();
921
922 const auto glob_i = local_dof_indices[i];
923 const auto glob_j = neighbor_local_dof_indices[f_index][j];
924 const auto m_i = lumped_mass_matrix_[glob_i];
925 const auto m_j = lumped_mass_matrix_[glob_j];
926 const auto hd_ij =
927 Number(0.5) * (m_i + m_j) / measure_of_omega_;
928
929 Number r_ij = 1.0;
930
931 if (ansatz == Ansatz::dg_q2) {
932 /*
933 * For even polynomial degree we normalize the incidence
934 * matrix to (0.5 (m_i + m_j) / |Omega|) ^ (1.5 / d).
935 * Note, that we will visit every coupling
936 * pair of degrees of freedom (i, j) precisely once.
937 */
938 r_ij = std::pow(hd_ij, incidence_relaxation_even_ / dim);
939 } else {
940 /*
941 * For odd polynomial degree we normalize the incidence
942 * matrix to 1. Note, that we will visit every coupling
943 * pair of degrees of freedom (i, j) precisely once.
944 */
945 r_ij = std::pow(hd_ij, incidence_relaxation_odd_ / dim);
946 }
947
948 interface_incidence_matrix[f_index](i, j) += r_ij;
949 }
950 } /* for i */
951 } /* for j */
952 } /* for q */
953 }
954 };
955
956 const auto copy_local_to_global = [&](const auto &copy) {
957 const auto &is_locally_owned = copy.is_locally_owned_;
958#ifdef DEAL_II_WITH_TRILINOS
959 const auto &local_dof_indices = copy.local_dof_indices_;
960 const auto &neighbor_local_dof_indices =
961 copy.neighbor_local_dof_indices_;
962#else
963 /*
964 * We have to transform indices to the local index range
965 * [0, n_locally_relevant_) when using the dealii::SparseMatrix.
966 * Thus, copy all index vectors:
967 */
968 auto local_dof_indices = copy.local_dof_indices_;
969 auto neighbor_local_dof_indices = copy.neighbor_local_dof_indices_;
970#endif
971 const auto &interface_incidence_matrix =
972 copy.interface_incidence_matrix_;
973
974 if (!is_locally_owned)
975 return;
976
977#ifndef DEAL_II_WITH_TRILINOS
978 transform_to_local_range(*scalar_partitioner_, local_dof_indices);
979 for (auto &indices : neighbor_local_dof_indices)
980 transform_to_local_range(*scalar_partitioner_, indices);
981#endif
982
983 for (unsigned int f_index = 0; f_index < copy.n_faces; ++f_index) {
984 if (neighbor_local_dof_indices[f_index].size() != 0) {
985 affine_constraints_assembly.distribute_local_to_global(
986 interface_incidence_matrix[f_index],
987 local_dof_indices,
988 neighbor_local_dof_indices[f_index],
989 incidence_matrix_tmp);
990 }
991 }
992 };
993
994 WorkStream::run(dof_handler.begin_active(),
995 dof_handler.end(),
996 local_assemble_system,
997 copy_local_to_global,
998 AssemblyScratchData<dim>(*discretization_),
999#ifdef DEAL_II_WITH_TRILINOS
1000 AssemblyCopyData<dim, double>());
1001#else
1002 AssemblyCopyData<dim, Number>());
1003#endif
1004
1005#ifdef DEAL_II_WITH_TRILINOS
1006 incidence_matrix_.read_in(incidence_matrix_tmp, /*locally_indexe*/ false);
1007#else
1008 incidence_matrix_.read_in(incidence_matrix_tmp, /*locally_indexed*/ true);
1009#endif
1010 incidence_matrix_.update_ghost_rows();
1011 }
1012
1013 /*
1014 * Populate boundary map and collect coupling boundary pairs:
1015 */
1016 {
1017 boundary_map_ = construct_boundary_map(
1018 dof_handler.begin_active(), dof_handler.end(), *scalar_partitioner_);
1019
1020 coupling_boundary_pairs_ = collect_coupling_boundary_pairs(
1021 dof_handler.begin_active(), dof_handler.end(), *scalar_partitioner_);
1022 }
1023
1024#ifdef DEBUG_SYMMETRY_CHECK
1025 /*
1026 * Verify that we have consistent mass:
1027 */
1028
1029 double total_mass = 0.;
1030 for (unsigned int i = 0; i < n_locally_owned_; ++i)
1031 total_mass += lumped_mass_matrix_.local_element(i);
1032 total_mass =
1033 Utilities::MPI::sum(total_mass, mpi_ensemble_.ensemble_communicator());
1034
1035 Assert(std::abs(measure_of_omega_ - total_mass) <
1036 1.e-12 * measure_of_omega_,
1037 dealii::ExcMessage(
1038 "Total mass differs from the measure of the domain."));
1039
1040 /*
1041 * Verify that the mij_matrix_ object is consistent:
1042 */
1043
1044 for (unsigned int i = 0; i < n_locally_owned_; ++i) {
1045 /* Skip constrained degrees of freedom: */
1046 const unsigned int row_length = sparsity_pattern_simd_.row_length(i);
1047 if (row_length == 1)
1048 continue;
1049
1050 auto sum =
1051 mass_matrix_.get_entry(i, 0) - lumped_mass_matrix_.local_element(i);
1052
1053 /* skip diagonal */
1054 constexpr auto simd_length = VectorizedArray<Number>::size();
1055 const unsigned int *js = sparsity_pattern_simd_.columns(i);
1056 for (unsigned int col_idx = 1; col_idx < row_length; ++col_idx) {
1057 const auto j = *(i < n_locally_internal_ ? js + col_idx * simd_length
1058 : js + col_idx);
1059 Assert(j < n_locally_relevant_, dealii::ExcInternalError());
1060
1061 const auto m_ij = mass_matrix_.get_entry(i, col_idx);
1062 if (discretization_->have_discontinuous_ansatz()) {
1063 // Interfacial coupling terms are present in the stencil but zero
1064 // in the mass matrix
1065 Assert(std::abs(m_ij) > -1.e-12, dealii::ExcInternalError());
1066 } else {
1067 Assert(std::abs(m_ij) > 1.e-12, dealii::ExcInternalError());
1068 }
1069 sum += m_ij;
1070
1071 const auto m_ji = mass_matrix_.get_transposed_entry(i, col_idx);
1072 if (std::abs(m_ij - m_ji) >= 1.e-12) {
1073 // The m_ij matrix is not symmetric
1074 std::stringstream ss;
1075 ss << "m_ij matrix is not symmetric: " << m_ij << " <-> " << m_ji;
1076 Assert(false, dealii::ExcMessage(ss.str()));
1077 }
1078 }
1079
1080 Assert(std::abs(sum) < 1.e-12, dealii::ExcInternalError());
1081 }
1082
1083 /*
1084 * Verify that the cij_matrix_ object is consistent:
1085 */
1086
1087 for (unsigned int i = 0; i < n_locally_owned_; ++i) {
1088 /* Skip constrained degrees of freedom: */
1089 const unsigned int row_length = sparsity_pattern_simd_.row_length(i);
1090 if (row_length == 1)
1091 continue;
1092
1093 auto sum = cij_matrix_.get_tensor(i, 0);
1094
1095 /* skip diagonal */
1096 constexpr auto simd_length = VectorizedArray<Number>::size();
1097 const unsigned int *js = sparsity_pattern_simd_.columns(i);
1098 for (unsigned int col_idx = 1; col_idx < row_length; ++col_idx) {
1099 const auto j = *(i < n_locally_internal_ ? js + col_idx * simd_length
1100 : js + col_idx);
1101 Assert(j < n_locally_relevant_, dealii::ExcInternalError());
1102
1103 const auto c_ij = cij_matrix_.get_tensor(i, col_idx);
1104 Assert(c_ij.norm() > 1.e-12, dealii::ExcInternalError());
1105 sum += c_ij;
1106
1107 const auto c_ji = cij_matrix_.get_transposed_tensor(i, col_idx);
1108 if ((c_ij + c_ji).norm() >= 1.e-12) {
1109 // The c_ij matrix is not symmetric, this can only happen if i
1110 // and j are both located on the boundary.
1111
1112 CouplingDescription coupling{i, col_idx, j};
1113 const auto it = std::find(coupling_boundary_pairs_.begin(),
1114 coupling_boundary_pairs_.end(),
1115 coupling);
1116 if (it == coupling_boundary_pairs_.end()) {
1117 std::stringstream ss;
1118 ss << "c_ij matrix is not anti-symmetric: " << c_ij << " <-> "
1119 << c_ji;
1120 Assert(false, dealii::ExcMessage(ss.str()));
1121 }
1122 }
1123 }
1124
1125 Assert(sum.norm() < 1.e-12, dealii::ExcInternalError());
1126 }
1127#endif
1128 }
1129
1130
1131 template <int dim, typename Number>
1132 void OfflineData<dim, Number>::create_multigrid_data()
1133 {
1134#ifdef DEBUG_OUTPUT
1135 std::cout << "OfflineData<dim, Number>::create_multigrid_data()"
1136 << std::endl;
1137#endif
1138
1139 auto &dof_handler = *dof_handler_;
1140
1141 Assert(!dof_handler.has_hp_capabilities(), dealii::ExcInternalError());
1142
1143 dof_handler.distribute_mg_dofs();
1144
1145 const auto n_levels = dof_handler.get_triangulation().n_global_levels();
1146
1147 AffineConstraints<float> level_constraints;
1148 // TODO not yet thread-parallel and without periodicity
1149
1150 level_boundary_map_.resize(n_levels);
1151 level_lumped_mass_matrix_.resize(n_levels);
1152
1153 for (unsigned int level = 0; level < n_levels; ++level) {
1154 /* Assemble lumped mass matrix vector: */
1155
1156 IndexSet relevant_dofs;
1157 dealii::DoFTools::extract_locally_relevant_level_dofs(
1158 dof_handler, level, relevant_dofs);
1159 const auto partitioner = std::make_shared<Utilities::MPI::Partitioner>(
1160 dof_handler.locally_owned_mg_dofs(level),
1161 relevant_dofs,
1162 mpi_ensemble_.ensemble_communicator());
1163 level_lumped_mass_matrix_[level].reinit(partitioner);
1164 std::vector<types::global_dof_index> dof_indices(
1165 dof_handler.get_fe().dofs_per_cell);
1166 dealii::Vector<Number> mass_values(dof_handler.get_fe().dofs_per_cell);
1167 dealii::hp::FEValues<dim> hp_fe_values(discretization_->mapping(),
1168 discretization_->finite_element(),
1169 discretization_->quadrature(),
1170 update_values | update_JxW_values);
1171 for (const auto &cell : dof_handler.cell_iterators_on_level(level))
1172 // TODO for assembly with dealii::SparseMatrix and local
1173 // numbering this probably has to read !cell->is_artificial()
1174 if (cell->is_locally_owned_on_level()) {
1175 hp_fe_values.reinit(cell);
1176 const auto &fe_values = hp_fe_values.get_present_fe_values();
1177 for (unsigned int i = 0; i < mass_values.size(); ++i) {
1178 double sum = 0;
1179 for (unsigned int q = 0; q < fe_values.n_quadrature_points; ++q)
1180 sum += fe_values.shape_value(i, q) * fe_values.JxW(q);
1181 mass_values(i) = sum;
1182 }
1183 cell->get_mg_dof_indices(dof_indices);
1184 level_constraints.distribute_local_to_global(
1185 mass_values, dof_indices, level_lumped_mass_matrix_[level]);
1186 }
1187 level_lumped_mass_matrix_[level].compress(VectorOperation::add);
1188
1189 /* Populate boundary map: */
1190
1191 level_boundary_map_[level] = construct_boundary_map(
1192 dof_handler.begin_mg(level), dof_handler.end_mg(level), *partitioner);
1193 }
1194 }
1195
1196
1197 template <int dim, typename Number>
1198 template <typename ITERATOR1, typename ITERATOR2>
1200 const ITERATOR1 &begin,
1201 const ITERATOR2 &end,
1202 const Utilities::MPI::Partitioner &partitioner) const -> BoundaryMap
1203 {
1204#ifdef DEBUG_OUTPUT
1205 std::cout << "OfflineData<dim, Number>::construct_boundary_map()"
1206 << std::endl;
1207#endif
1208
1209 /*
1210 * Create a temporary multimap with the (local) dof index as key:
1211 */
1212
1213 using BoundaryData = std::tuple<dealii::Tensor<1, dim, Number> /*normal*/,
1214 Number /*normal mass*/,
1215 Number /*boundary mass*/,
1216 dealii::types::boundary_id /*id*/,
1217 dealii::Point<dim>> /*position*/;
1218 std::multimap<unsigned int, BoundaryData> preliminary_map;
1219
1220 std::vector<dealii::types::global_dof_index> local_dof_indices;
1221
1222 dealii::hp::FEFaceValues<dim> hp_fe_face_values(
1223 discretization_->mapping(),
1224 discretization_->finite_element(),
1225 discretization_->face_quadrature(),
1226 dealii::update_normal_vectors | dealii::update_values |
1227 dealii::update_JxW_values);
1228
1229 for (auto cell = begin; cell != end; ++cell) {
1230
1231 /*
1232 * Workaround: Make sure to iterate over the entire locally relevant
1233 * set so that we compute normals between cells with differing owners
1234 * correctly. This strategy works for 2D but will fail in 3D with
1235 * locally refined meshes and hanging nodes situated at the boundary.
1236 */
1237 if ((cell->is_active() && cell->is_artificial()) ||
1238 (!cell->is_active() && cell->is_artificial_on_level()))
1239 continue;
1240
1241 const unsigned int dofs_per_cell = cell->get_fe().n_dofs_per_cell();
1242 const auto &support_points = cell->get_fe().get_unit_support_points();
1243
1244 local_dof_indices.resize(dofs_per_cell);
1245 cell->get_active_or_mg_dof_indices(local_dof_indices);
1246
1247 for (auto f : cell->face_indices()) {
1248 const auto face = cell->face(f);
1249 const auto id = face->boundary_id();
1250
1251 // FIXME: support interior boundary with FE_Nothing
1252
1253 if (!face->at_boundary())
1254 continue;
1255
1256 /*
1257 * Skip periodic boundary faces. For our algorithm these are
1258 * interior degrees of freedom (if not simultaneously located at
1259 * another boundary as well).
1260 */
1261 if (id == Boundary::periodic)
1262 continue;
1263
1264 hp_fe_face_values.reinit(cell, f);
1265 const auto &fe_face_values = hp_fe_face_values.get_present_fe_values();
1266 const auto &mapping =
1267 hp_fe_face_values.get_mapping_collection()[cell->active_fe_index()];
1268
1269 for (unsigned int j : fe_face_values.dof_indices()) {
1270 if (!cell->get_fe().has_support_on_face(j, f))
1271 continue;
1272
1273 Number boundary_mass = 0.;
1274 dealii::Tensor<1, dim, Number> normal;
1275
1276 for (unsigned int q : fe_face_values.quadrature_point_indices()) {
1277 const auto JxW = fe_face_values.JxW(q);
1278 const auto phi_i = fe_face_values.shape_value(j, q);
1279
1280 boundary_mass += phi_i * JxW;
1281 normal += phi_i * fe_face_values.normal_vector(q) * JxW;
1282 }
1283
1284 const auto global_index = local_dof_indices[j];
1285 const auto index = partitioner.global_to_local(global_index);
1286
1287 /* Skip nonlocal degrees of freedom: */
1288 if (index >= n_locally_owned_)
1289 continue;
1290
1291 /* Skip constrained degrees of freedom: */
1292 const unsigned int row_length =
1293 sparsity_pattern_simd_.row_length(index);
1294 if (row_length == 1)
1295 continue;
1296
1297 Point<dim> position =
1298 mapping.transform_unit_to_real_cell(cell, support_points[j]);
1299
1300 /*
1301 * Temporarily insert a (wrong) boundary mass value for the
1302 * normal mass. We'll fix this later.
1303 */
1304 preliminary_map.insert(
1305 {index, {normal, boundary_mass, boundary_mass, id, position}});
1306 } /* j */
1307 } /* f */
1308 } /* cell */
1309
1310 /*
1311 * Filter boundary map:
1312 *
1313 * At this point we have collected multiple cell contributions for each
1314 * boundary degree of freedom. We now merge all entries that have the
1315 * same boundary id and whose normals describe an acute angle of about
1316 * 60 degrees or less.
1317 *
1318 * FIXME: is this robust in 3D?
1319 */
1320
1321 std::multimap<unsigned int, BoundaryData> filtered_map;
1322 std::set<dealii::types::global_dof_index> boundary_dofs;
1323 for (auto entry : preliminary_map) {
1324 bool inserted = false;
1325 const auto range = filtered_map.equal_range(entry.first);
1326 for (auto it = range.first; it != range.second; ++it) {
1327 auto &[new_normal,
1328 new_normal_mass,
1329 new_boundary_mass,
1330 new_id,
1331 new_point] = entry.second;
1332 auto &[normal, normal_mass, boundary_mass, id, point] = it->second;
1333
1334 if (id != new_id)
1335 continue;
1336
1337 Assert(point.distance(new_point) < 1.0e-14, dealii::ExcInternalError());
1338
1339 if (normal * new_normal / normal.norm() / new_normal.norm() > 0.50) {
1340 /*
1341 * Both normals describe an acute angle of 85 degrees or less.
1342 * Merge the entries and continue.
1343 */
1344 normal += new_normal;
1345 boundary_mass += new_boundary_mass;
1346 inserted = true;
1347 continue;
1348
1349 } else if constexpr (dim == 2) {
1350 /*
1351 * Workaround for 2D: When enforcing slip boundary conditions
1352 * with two noncollinear vectors the resulting momentum must be
1353 * 0. But the normals don't necessarily describe an orthonormal
1354 * basis and we cannot use orthogonal projection. Therefore,
1355 * simply set the boundary type to no slip:
1356 */
1357 if (new_id == Boundary::slip) {
1358 Assert(id == Boundary::slip, dealii::ExcInternalError());
1359 new_id = Boundary::no_slip;
1360 id = Boundary::no_slip;
1361 }
1362 }
1363 }
1364 if (!inserted)
1365 filtered_map.insert(entry);
1366 }
1367
1368 /*
1369 * Normalize all normal vectors and create final boundary_map:
1370 */
1371
1372 BoundaryMap boundary_map;
1373 std::transform(
1374 std::begin(filtered_map),
1375 std::end(filtered_map),
1376 std::back_inserter(boundary_map),
1377 [&](const auto &it) -> BoundaryDescription { //
1378 auto index = it.first;
1379 const auto &[normal, normal_mass, boundary_mass, id, point] =
1380 it.second;
1381
1382 const auto new_normal_mass =
1383 normal.norm() + std::numeric_limits<Number>::epsilon();
1384 const auto new_normal = normal / new_normal_mass;
1385
1386 return {index, new_normal, new_normal_mass, boundary_mass, id, point};
1387 });
1388
1389 return boundary_map;
1390 }
1391
1392
1393 template <int dim, typename Number>
1394 template <typename ITERATOR1, typename ITERATOR2>
1396 const ITERATOR1 &begin,
1397 const ITERATOR2 &end,
1398 const Utilities::MPI::Partitioner &partitioner) const
1399 -> CouplingBoundaryPairs
1400 {
1401#ifdef DEBUG_OUTPUT
1402 std::cout << "OfflineData<dim, Number>::collect_coupling_boundary_pairs()"
1403 << std::endl;
1404#endif
1405
1406 /*
1407 * First, collect *all* locally relevant degrees of freedom that are
1408 * located on a (non periodic) boundary. We also collect constrained
1409 * degrees of freedom for the time being (and filter afterwards).
1410 */
1411
1412 std::set<unsigned int> locally_relevant_boundary_indices;
1413
1414 std::vector<dealii::types::global_dof_index> local_dof_indices;
1415
1416 for (auto cell = begin; cell != end; ++cell) {
1417
1418 /* Make sure to iterate over the entire locally relevant set: */
1419 if (cell->is_artificial())
1420 continue;
1421
1422 const auto &finite_element = cell->get_fe();
1423 const unsigned int dofs_per_cell = finite_element.dofs_per_cell;
1424 local_dof_indices.resize(dofs_per_cell);
1425 cell->get_active_or_mg_dof_indices(local_dof_indices);
1426
1427 for (auto f : cell->face_indices()) {
1428 const auto face = cell->face(f);
1429 const auto id = face->boundary_id();
1430
1431 if (!face->at_boundary())
1432 continue;
1433
1434 /* Skip periodic boundary faces; see above. */
1435 if (id == Boundary::periodic)
1436 continue;
1437
1438 for (unsigned int j = 0; j < dofs_per_cell; ++j) {
1439
1440 if (!cell->get_fe().has_support_on_face(j, f))
1441 continue;
1442
1443 const auto global_index = local_dof_indices[j];
1444 const auto index = partitioner.global_to_local(global_index);
1445
1446 /* Skip irrelevant degrees of freedom: */
1447 if (index >= n_locally_relevant_)
1448 continue;
1449
1450 locally_relevant_boundary_indices.insert(index);
1451 } /* j */
1452 } /* f */
1453 } /* cell */
1454
1455 /*
1456 * Now, collect all coupling boundary pairs:
1457 */
1458
1459 CouplingBoundaryPairs result;
1460
1461 for (const auto i : locally_relevant_boundary_indices) {
1462
1463 /* Only record pairs with a left index that is locally owned: */
1464 if (i >= n_locally_owned_)
1465 continue;
1466
1467 const unsigned int row_length = sparsity_pattern_simd_.row_length(i);
1468
1469 /* Skip all constrained degrees of freedom: */
1470 if (row_length == 1)
1471 continue;
1472
1473 const unsigned int *js = sparsity_pattern_simd_.columns(i);
1474 constexpr auto simd_length = VectorizedArray<Number>::size();
1475 /* skip diagonal: */
1476 for (unsigned int col_idx = 1; col_idx < row_length; ++col_idx) {
1477 const auto j = *(i < n_locally_internal_ ? js + col_idx * simd_length
1478 : js + col_idx);
1479
1480 if (locally_relevant_boundary_indices.count(j) != 0) {
1481 result.push_back({i, col_idx, j});
1482 }
1483 }
1484 }
1485
1486 return result;
1487 }
1488
1489} /* namespace ryujin */
std::tuple< unsigned int, dealii::Tensor< 1, dim, Number >, Number, Number, dealii::types::boundary_id, dealii::Point< dim > > BoundaryDescription
Definition: offline_data.h:76
OfflineData(const MPIEnsemble &mpi_ensemble, const Discretization< dim > &discretization, const std::string &subsection="/OfflineData")
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)
std::shared_ptr< const dealii::Utilities::MPI::Partitioner > create_vector_partitioner(const std::shared_ptr< const dealii::Utilities::MPI::Partitioner > &scalar_partitioner, const unsigned int n_components)
T pow(const T x, const T b)
dealii::LinearAlgebra::distributed::Vector< Number > ScalarVector
Definition: state_vector.h:31
DEAL_II_ALWAYS_INLINE FT add(const FT &flux_left_ij, const FT &flux_right_ij)