ryujin 2.1.1 revision c38afcfdec63b34961924f05408760247496a1f0
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 elminiating 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 * Next we can (re)initialize all local matrices:
410 */
411
412 mass_matrix_.reinit(sparsity_pattern_simd_);
413 if (discretization_->have_discontinuous_ansatz())
414 mass_matrix_inverse_.reinit(sparsity_pattern_simd_);
415
416 lumped_mass_matrix_.reinit(scalar_partitioner_);
417 lumped_mass_matrix_inverse_.reinit(scalar_partitioner_);
418
419 cij_matrix_.reinit(sparsity_pattern_simd_);
420 if (discretization_->have_discontinuous_ansatz())
421 incidence_matrix_.reinit(sparsity_pattern_simd_);
422 }
423
424
425 template <int dim, typename Number>
426 void OfflineData<dim, Number>::assemble()
427 {
428#ifdef DEBUG_OUTPUT
429 std::cout << "OfflineData<dim, Number>::assemble()" << std::endl;
430#endif
431
432 auto &dof_handler = *dof_handler_;
433
434 measure_of_omega_ = 0.;
435
436#ifdef DEAL_II_WITH_TRILINOS
437 /* Variant using TrilinosWrappers::SparseMatrix with global numbering */
438
439 AffineConstraints<double> affine_constraints_assembly;
440 /* This small dance is necessary to translate from Number to double: */
441 affine_constraints_assembly.reinit(affine_constraints_.get_local_lines());
442 for (auto line : affine_constraints_.get_lines()) {
443 affine_constraints_assembly.add_line(line.index);
444 for (auto entry : line.entries)
445 affine_constraints_assembly.add_entry(
446 line.index, entry.first, entry.second);
447 affine_constraints_assembly.set_inhomogeneity(line.index,
448 line.inhomogeneity);
449 }
450 affine_constraints_assembly.close();
451
452 const IndexSet &locally_owned = dof_handler.locally_owned_dofs();
453 TrilinosWrappers::SparsityPattern trilinos_sparsity_pattern;
454 trilinos_sparsity_pattern.reinit(locally_owned,
455 sparsity_pattern_,
456 mpi_ensemble_.ensemble_communicator());
457
458 TrilinosWrappers::SparseMatrix mass_matrix_tmp;
459 TrilinosWrappers::SparseMatrix mass_matrix_inverse_tmp;
460 if (discretization_->have_discontinuous_ansatz())
461 mass_matrix_inverse_tmp.reinit(trilinos_sparsity_pattern);
462
463 std::array<TrilinosWrappers::SparseMatrix, dim> cij_matrix_tmp;
464
465 mass_matrix_tmp.reinit(trilinos_sparsity_pattern);
466 for (auto &matrix : cij_matrix_tmp)
467 matrix.reinit(trilinos_sparsity_pattern);
468
469#else
470 /* Variant using deal.II SparseMatrix with local numbering */
471
472 AffineConstraints<Number> affine_constraints_assembly;
473 affine_constraints_assembly.copy_from(affine_constraints_);
474 transform_to_local_range(*scalar_partitioner_, affine_constraints_assembly);
475
476 SparsityPattern sparsity_pattern_assembly;
477 {
478 DynamicSparsityPattern dsp(n_locally_relevant_, n_locally_relevant_);
479 for (const auto &entry : sparsity_pattern_) {
480 const auto i = scalar_partitioner_->global_to_local(entry.row());
481 const auto j = scalar_partitioner_->global_to_local(entry.column());
482 dsp.add(i, j);
483 }
484 sparsity_pattern_assembly.copy_from(dsp);
485 }
486
487 dealii::SparseMatrix<Number> mass_matrix_tmp;
488 dealii::SparseMatrix<Number> mass_matrix_inverse_tmp;
489 if (discretization_->have_discontinuous_ansatz())
490 mass_matrix_inverse_tmp.reinit(sparsity_pattern_assembly);
491
492 std::array<dealii::SparseMatrix<Number>, dim> cij_matrix_tmp;
493
494 mass_matrix_tmp.reinit(sparsity_pattern_assembly);
495 for (auto &matrix : cij_matrix_tmp)
496 matrix.reinit(sparsity_pattern_assembly);
497#endif
498
499 const unsigned int dofs_per_cell =
500 discretization_->finite_element().dofs_per_cell;
501
502 /*
503 * Now, assemble all matrices:
504 */
505
506 /* The local, per-cell assembly routine: */
507 const auto local_assemble_system = [&](const auto &cell,
508 auto &scratch,
509 auto &copy) {
510 /* iterate over locally owned cells and the ghost layer */
511
512 auto &is_locally_owned = copy.is_locally_owned_;
513 auto &local_dof_indices = copy.local_dof_indices_;
514 auto &neighbor_local_dof_indices = copy.neighbor_local_dof_indices_;
515
516 auto &cell_mass_matrix = copy.cell_mass_matrix_;
517 auto &cell_mass_matrix_inverse = copy.cell_mass_matrix_inverse_;
518 auto &cell_cij_matrix = copy.cell_cij_matrix_;
519 auto &interface_cij_matrix = copy.interface_cij_matrix_;
520 auto &cell_measure = copy.cell_measure_;
521
522 auto &fe_values = scratch.fe_values_;
523 auto &fe_face_values = scratch.fe_face_values_;
524 auto &fe_neighbor_face_values = scratch.fe_neighbor_face_values_;
525
526#ifdef DEAL_II_WITH_TRILINOS
527 is_locally_owned = cell->is_locally_owned();
528#else
529 /*
530 * When using a local dealii::SparseMatrix<Number> we don not
531 * have a compress(VectorOperation::add) available. In this case
532 * we assemble contributions over all locally relevant (non
533 * artificial) cells.
534 */
535 is_locally_owned = !cell->is_artificial();
536#endif
537 if (!is_locally_owned)
538 return;
539
540 cell_mass_matrix.reinit(dofs_per_cell, dofs_per_cell);
541 for (auto &matrix : cell_cij_matrix)
542 matrix.reinit(dofs_per_cell, dofs_per_cell);
543 if (discretization_->have_discontinuous_ansatz()) {
544 cell_mass_matrix_inverse.reinit(dofs_per_cell, dofs_per_cell);
545 for (auto &it : interface_cij_matrix)
546 for (auto &matrix : it)
547 matrix.reinit(dofs_per_cell, dofs_per_cell);
548 }
549
550 fe_values.reinit(cell);
551
552 local_dof_indices.resize(dofs_per_cell);
553 cell->get_dof_indices(local_dof_indices);
554
555 /* clear out copy data: */
556 cell_mass_matrix = 0.;
557 for (auto &matrix : cell_cij_matrix)
558 matrix = 0.;
559 if (discretization_->have_discontinuous_ansatz()) {
560 cell_mass_matrix_inverse = 0.;
561 for (auto &it : interface_cij_matrix)
562 for (auto &matrix : it)
563 matrix = 0.;
564 }
565 cell_measure = 0.;
566
567 for (unsigned int q : fe_values.quadrature_point_indices()) {
568 const auto JxW = fe_values.JxW(q);
569
570 if (cell->is_locally_owned())
571 cell_measure += Number(JxW);
572
573 for (unsigned int j : fe_values.dof_indices()) {
574 const auto value_JxW = fe_values.shape_value(j, q) * JxW;
575 const auto grad_JxW = fe_values.shape_grad(j, q) * JxW;
576
577 for (unsigned int i : fe_values.dof_indices()) {
578 const auto value = fe_values.shape_value(i, q);
579
580 cell_mass_matrix(i, j) += Number(value * value_JxW);
581 for (unsigned int d = 0; d < dim; ++d)
582 cell_cij_matrix[d](i, j) += Number((value * grad_JxW)[d]);
583 } /* for i */
584 } /* for j */
585 } /* for q */
586
587 /*
588 * For a discontinuous finite element ansatz we need to assemble
589 * additional face contributions:
590 */
591
592 if (!discretization_->have_discontinuous_ansatz())
593 return;
594
595 for (const auto f_index : cell->face_indices()) {
596 const auto &face = cell->face(f_index);
597
598 /* Skip faces without neighbors... */
599 const bool has_neighbor =
600 !face->at_boundary() || cell->has_periodic_neighbor(f_index);
601 if (!has_neighbor) {
602 // set the vector of local dof indices to 0 to indicate that
603 // there is nothing to do for this face:
604 neighbor_local_dof_indices[f_index].resize(0);
605 continue;
606 }
607
608 /* Avoid artificial cells: */
609 const auto neighbor_cell = cell->neighbor_or_periodic_neighbor(f_index);
610 if (neighbor_cell->is_artificial()) {
611 // set the vector of local dof indices to 0 to indicate that
612 // there is nothing to do for this face:
613 neighbor_local_dof_indices[f_index].resize(0);
614 continue;
615 }
616
617 fe_face_values.reinit(cell, f_index);
618
619 /* Face contribution: */
620
621 for (unsigned int q : fe_face_values.quadrature_point_indices()) {
622 const auto JxW = fe_face_values.JxW(q);
623 const auto &normal = fe_face_values.get_normal_vectors()[q];
624
625 for (unsigned int j : fe_face_values.dof_indices()) {
626 const auto value_JxW = fe_face_values.shape_value(j, q) * JxW;
627
628 for (unsigned int i : fe_face_values.dof_indices()) {
629 const auto value = fe_face_values.shape_value(i, q);
630
631 for (unsigned int d = 0; d < dim; ++d)
632 cell_cij_matrix[d](i, j) -=
633 Number(0.5 * normal[d] * value * value_JxW);
634 } /* for i */
635 } /* for j */
636 } /* for q */
637
638 /* Coupling part: */
639
640 const unsigned int f_index_neighbor =
641 cell->has_periodic_neighbor(f_index)
642 ? cell->periodic_neighbor_of_periodic_neighbor(f_index)
643 : cell->neighbor_of_neighbor(f_index);
644
645 neighbor_local_dof_indices[f_index].resize(dofs_per_cell);
646 neighbor_cell->get_dof_indices(neighbor_local_dof_indices[f_index]);
647
648 fe_neighbor_face_values.reinit(neighbor_cell, f_index_neighbor);
649
650 for (unsigned int q : fe_face_values.quadrature_point_indices()) {
651 const auto JxW = fe_face_values.JxW(q);
652 const auto &normal = fe_face_values.get_normal_vectors()[q];
653
654 /* index j for neighbor, index i for current cell: */
655 for (unsigned int j : fe_face_values.dof_indices()) {
656 const auto value_JxW =
657 fe_neighbor_face_values.shape_value(j, q) * JxW;
658
659 for (unsigned int i : fe_face_values.dof_indices()) {
660 const auto value = fe_face_values.shape_value(i, q);
661
662 for (unsigned int d = 0; d < dim; ++d)
663 interface_cij_matrix[f_index][d](i, j) +=
664 Number(0.5 * normal[d] * value * value_JxW);
665 } /* for i */
666 } /* for j */
667 } /* for q */
668 }
669
670 /*
671 * Compute block inverse of mass matrix:
672 */
673
674 if (discretization_->have_discontinuous_ansatz()) {
675 // FIXME: rewrite with CellwiseInverseMassMatrix
676 cell_mass_matrix_inverse.invert(cell_mass_matrix);
677 }
678 };
679
680 const auto copy_local_to_global = [&](const auto &copy) {
681 const auto &is_locally_owned = copy.is_locally_owned_;
682#ifdef DEAL_II_WITH_TRILINOS
683 const auto &local_dof_indices = copy.local_dof_indices_;
684 const auto &neighbor_local_dof_indices = copy.neighbor_local_dof_indices_;
685#else
686 /*
687 * We have to transform indices to the local index range
688 * [0, n_locally_relevant_) when using the dealii::SparseMatrix.
689 * Thus, copy all index vectors:
690 */
691 auto local_dof_indices = copy.local_dof_indices_;
692 auto neighbor_local_dof_indices = copy.neighbor_local_dof_indices_;
693#endif
694 const auto &cell_mass_matrix = copy.cell_mass_matrix_;
695 const auto &cell_mass_matrix_inverse = copy.cell_mass_matrix_inverse_;
696 const auto &cell_cij_matrix = copy.cell_cij_matrix_;
697 const auto &interface_cij_matrix = copy.interface_cij_matrix_;
698 const auto &cell_measure = copy.cell_measure_;
699
700 if (!is_locally_owned)
701 return;
702
703#ifndef DEAL_II_WITH_TRILINOS
704 transform_to_local_range(*scalar_partitioner_, local_dof_indices);
705 for (auto &indices : neighbor_local_dof_indices)
706 transform_to_local_range(*scalar_partitioner_, indices);
707#endif
708
709 affine_constraints_assembly.distribute_local_to_global(
710 cell_mass_matrix, local_dof_indices, mass_matrix_tmp);
711
712 for (int k = 0; k < dim; ++k) {
713 affine_constraints_assembly.distribute_local_to_global(
714 cell_cij_matrix[k], local_dof_indices, cij_matrix_tmp[k]);
715
716 for (unsigned int f_index = 0; f_index < copy.n_faces; ++f_index) {
717 if (neighbor_local_dof_indices[f_index].size() != 0) {
718 affine_constraints_assembly.distribute_local_to_global(
719 interface_cij_matrix[f_index][k],
720 local_dof_indices,
721 neighbor_local_dof_indices[f_index],
722 cij_matrix_tmp[k]);
723 }
724 }
725 }
726
727 if (discretization_->have_discontinuous_ansatz())
728 affine_constraints_assembly.distribute_local_to_global(
729 cell_mass_matrix_inverse,
730 local_dof_indices,
731 mass_matrix_inverse_tmp);
732
733 measure_of_omega_ += cell_measure;
734 };
735
736 WorkStream::run(dof_handler.begin_active(),
737 dof_handler.end(),
738 local_assemble_system,
739 copy_local_to_global,
740 AssemblyScratchData<dim>(*discretization_),
741#ifdef DEAL_II_WITH_TRILINOS
742 AssemblyCopyData<dim, double>());
743#else
744 AssemblyCopyData<dim, Number>());
745#endif
746
747#ifdef DEAL_II_WITH_TRILINOS
748 mass_matrix_tmp.compress(VectorOperation::add);
749 for (auto &it : cij_matrix_tmp)
750 it.compress(VectorOperation::add);
751
752 mass_matrix_.read_in(mass_matrix_tmp, /*locally_indexed*/ false);
753 if (discretization_->have_discontinuous_ansatz())
754 mass_matrix_inverse_.read_in(mass_matrix_inverse_tmp, /*loc_ind*/ false);
755 cij_matrix_.read_in(cij_matrix_tmp, /*locally_indexed*/ false);
756#else
757 mass_matrix_.read_in(mass_matrix_tmp, /*locally_indexed*/ true);
758 if (discretization_->have_discontinuous_ansatz())
759 mass_matrix_inverse_.read_in(mass_matrix_inverse_tmp, /*loc_ind*/ true);
760 cij_matrix_.read_in(cij_matrix_tmp, /*locally_indexed*/ true);
761#endif
762
763 mass_matrix_.update_ghost_rows();
764 if (discretization_->have_discontinuous_ansatz())
765 mass_matrix_inverse_.update_ghost_rows();
766 cij_matrix_.update_ghost_rows();
767
768 measure_of_omega_ = Utilities::MPI::sum(
769 measure_of_omega_, mpi_ensemble_.ensemble_communicator());
770
771 /*
772 * Create lumped mass matrix:
773 */
774
775 {
776#ifdef DEAL_II_WITH_TRILINOS
777 ScalarVector one(scalar_partitioner_);
778 one = 1.;
779 affine_constraints_assembly.set_zero(one);
780
781 ScalarVector local_lumped_mass_matrix(scalar_partitioner_);
782 mass_matrix_tmp.vmult(local_lumped_mass_matrix, one);
783 local_lumped_mass_matrix.compress(VectorOperation::add);
784
785 for (unsigned int i = 0; i < scalar_partitioner_->locally_owned_size();
786 ++i) {
787 lumped_mass_matrix_.local_element(i) =
788 local_lumped_mass_matrix.local_element(i);
789 lumped_mass_matrix_inverse_.local_element(i) =
790 1. / lumped_mass_matrix_.local_element(i);
791 }
792 lumped_mass_matrix_.update_ghost_values();
793 lumped_mass_matrix_inverse_.update_ghost_values();
794
795#else
796
797 dealii::Vector<Number> one(mass_matrix_tmp.m());
798 one = 1.;
799 affine_constraints_assembly.set_zero(one);
800
801 dealii::Vector<Number> local_lumped_mass_matrix(mass_matrix_tmp.m());
802 mass_matrix_tmp.vmult(local_lumped_mass_matrix, one);
803
804 for (unsigned int i = 0; i < scalar_partitioner_->locally_owned_size();
805 ++i) {
806 lumped_mass_matrix_.local_element(i) = local_lumped_mass_matrix(i);
807 lumped_mass_matrix_inverse_.local_element(i) =
808 1. / lumped_mass_matrix_.local_element(i);
809 }
810 lumped_mass_matrix_.update_ghost_values();
811 lumped_mass_matrix_inverse_.update_ghost_values();
812#endif
813 }
814
815 /*
816 * Assemble incidence matrix:
817 */
818
819 if (discretization_->have_discontinuous_ansatz()) {
820#ifdef DEAL_II_WITH_TRILINOS
821 TrilinosWrappers::SparseMatrix incidence_matrix_tmp;
822 incidence_matrix_tmp.reinit(trilinos_sparsity_pattern);
823#else
824 dealii::SparseMatrix<Number> incidence_matrix_tmp;
825 incidence_matrix_tmp.reinit(sparsity_pattern_assembly);
826#endif
827
828 /* The local, per-cell assembly routine: */
829 const auto local_assemble_system = [&](const auto &cell,
830 auto &scratch,
831 auto &copy) {
832 /* iterate over locally owned cells and the ghost layer */
833
834 auto &is_locally_owned = copy.is_locally_owned_;
835 auto &local_dof_indices = copy.local_dof_indices_;
836 auto &neighbor_local_dof_indices = copy.neighbor_local_dof_indices_;
837 auto &interface_incidence_matrix = copy.interface_incidence_matrix_;
838 auto &fe_face_values_nodal = scratch.fe_face_values_nodal_;
839 auto &fe_neighbor_face_values_nodal =
840 scratch.fe_neighbor_face_values_nodal_;
841
842#ifdef DEAL_II_WITH_TRILINOS
843 is_locally_owned = cell->is_locally_owned();
844#else
845 is_locally_owned = !cell->is_artificial();
846#endif
847 if (!is_locally_owned)
848 return;
849
850 for (auto &matrix : interface_incidence_matrix)
851 matrix.reinit(dofs_per_cell, dofs_per_cell);
852
853 local_dof_indices.resize(dofs_per_cell);
854 cell->get_dof_indices(local_dof_indices);
855
856 /* clear out copy data: */
857 for (auto &matrix : interface_incidence_matrix)
858 matrix = 0.;
859
860 for (const auto f_index : cell->face_indices()) {
861 const auto &face = cell->face(f_index);
862
863 /* Skip faces without neighbors... */
864 const bool has_neighbor =
865 !face->at_boundary() || cell->has_periodic_neighbor(f_index);
866 if (!has_neighbor) {
867 // set the vector of local dof indices to 0 to indicate that
868 // there is nothing to do for this face:
869 neighbor_local_dof_indices[f_index].resize(0);
870 continue;
871 }
872
873 /* Avoid artificial cells: */
874 const auto neighbor_cell =
875 cell->neighbor_or_periodic_neighbor(f_index);
876 if (neighbor_cell->is_artificial()) {
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 const unsigned int f_index_neighbor =
884 cell->has_periodic_neighbor(f_index)
885 ? cell->periodic_neighbor_of_periodic_neighbor(f_index)
886 : cell->neighbor_of_neighbor(f_index);
887
888 neighbor_local_dof_indices[f_index].resize(dofs_per_cell);
889 neighbor_cell->get_dof_indices(neighbor_local_dof_indices[f_index]);
890
891 fe_face_values_nodal.reinit(cell, f_index);
892 fe_neighbor_face_values_nodal.reinit(neighbor_cell, f_index_neighbor);
893
894 /* Lumped incidence matrix: */
895
896 for (unsigned int q :
897 fe_face_values_nodal.quadrature_point_indices()) {
898 /* index j for neighbor, index i for current cell: */
899 for (unsigned int j : fe_neighbor_face_values_nodal.dof_indices()) {
900 const auto v_j = fe_neighbor_face_values_nodal.shape_value(j, q);
901 for (unsigned int i : fe_face_values_nodal.dof_indices()) {
902 const auto v_i = fe_face_values_nodal.shape_value(i, q);
903 constexpr auto eps = std::numeric_limits<Number>::epsilon();
904 if (std::abs(v_i * v_j) > 100. * eps) {
905 const auto &ansatz = discretization_->ansatz();
906
907 const auto glob_i = local_dof_indices[i];
908 const auto glob_j = neighbor_local_dof_indices[f_index][j];
909 const auto m_i = lumped_mass_matrix_[glob_i];
910 const auto m_j = lumped_mass_matrix_[glob_j];
911 const auto hd_ij =
912 Number(0.5) * (m_i + m_j) / measure_of_omega_;
913
914 Number r_ij = 1.0;
915
916 if (ansatz == Ansatz::dg_q0 || ansatz == Ansatz::dg_q2) {
917 /*
918 * For even polynomial degree we normalize the incidence
919 * matrix to (0.5 (m_i + m_j) / |Omega|) ^ (1.5 / d).
920 * Note, that we will visit every coupling
921 * pair of degrees of freedom (i, j) precisely once.
922 */
923 r_ij = std::pow(hd_ij, incidence_relaxation_even_ / dim);
924 } else {
925 /*
926 * For odd polynomial degree we normalize the incidence
927 * matrix to 1. Note, that we will visit every coupling
928 * pair of degrees of freedom (i, j) precisely once.
929 */
930 r_ij = std::pow(hd_ij, incidence_relaxation_odd_ / dim);
931 }
932
933 interface_incidence_matrix[f_index](i, j) += r_ij;
934 }
935 } /* for i */
936 } /* for j */
937 } /* for q */
938 }
939 };
940
941 const auto copy_local_to_global = [&](const auto &copy) {
942 const auto &is_locally_owned = copy.is_locally_owned_;
943#ifdef DEAL_II_WITH_TRILINOS
944 const auto &local_dof_indices = copy.local_dof_indices_;
945 const auto &neighbor_local_dof_indices =
946 copy.neighbor_local_dof_indices_;
947#else
948 /*
949 * We have to transform indices to the local index range
950 * [0, n_locally_relevant_) when using the dealii::SparseMatrix.
951 * Thus, copy all index vectors:
952 */
953 auto local_dof_indices = copy.local_dof_indices_;
954 auto neighbor_local_dof_indices = copy.neighbor_local_dof_indices_;
955#endif
956 const auto &interface_incidence_matrix =
957 copy.interface_incidence_matrix_;
958
959 if (!is_locally_owned)
960 return;
961
962#ifndef DEAL_II_WITH_TRILINOS
963 transform_to_local_range(*scalar_partitioner_, local_dof_indices);
964 for (auto &indices : neighbor_local_dof_indices)
965 transform_to_local_range(*scalar_partitioner_, indices);
966#endif
967
968 for (unsigned int f_index = 0; f_index < copy.n_faces; ++f_index) {
969 if (neighbor_local_dof_indices[f_index].size() != 0) {
970 affine_constraints_assembly.distribute_local_to_global(
971 interface_incidence_matrix[f_index],
972 local_dof_indices,
973 neighbor_local_dof_indices[f_index],
974 incidence_matrix_tmp);
975 }
976 }
977 };
978
979 WorkStream::run(dof_handler.begin_active(),
980 dof_handler.end(),
981 local_assemble_system,
982 copy_local_to_global,
983 AssemblyScratchData<dim>(*discretization_),
984#ifdef DEAL_II_WITH_TRILINOS
985 AssemblyCopyData<dim, double>());
986#else
987 AssemblyCopyData<dim, Number>());
988#endif
989
990#ifdef DEAL_II_WITH_TRILINOS
991 incidence_matrix_.read_in(incidence_matrix_tmp, /*locally_indexe*/ false);
992#else
993 incidence_matrix_.read_in(incidence_matrix_tmp, /*locally_indexed*/ true);
994#endif
995 incidence_matrix_.update_ghost_rows();
996 }
997
998 /*
999 * Populate boundary map and collect coupling boundary pairs:
1000 */
1001 {
1002 boundary_map_ = construct_boundary_map(
1003 dof_handler.begin_active(), dof_handler.end(), *scalar_partitioner_);
1004
1005 coupling_boundary_pairs_ = collect_coupling_boundary_pairs(
1006 dof_handler.begin_active(), dof_handler.end(), *scalar_partitioner_);
1007 }
1008
1009#ifdef DEBUG
1010 /*
1011 * Verify that we have consistent mass:
1012 */
1013
1014 double total_mass = 0.;
1015 for (unsigned int i = 0; i < n_locally_owned_; ++i)
1016 total_mass += lumped_mass_matrix_.local_element(i);
1017 total_mass =
1018 Utilities::MPI::sum(total_mass, mpi_ensemble_.ensemble_communicator());
1019
1020 Assert(std::abs(measure_of_omega_ - total_mass) <
1021 1.e-12 * measure_of_omega_,
1022 dealii::ExcMessage(
1023 "Total mass differs from the measure of the domain."));
1024
1025 /*
1026 * Verify that the mij_matrix_ object is consistent:
1027 */
1028
1029 for (unsigned int i = 0; i < n_locally_owned_; ++i) {
1030 /* Skip constrained degrees of freedom: */
1031 const unsigned int row_length = sparsity_pattern_simd_.row_length(i);
1032 if (row_length == 1)
1033 continue;
1034
1035 auto sum =
1036 mass_matrix_.get_entry(i, 0) - lumped_mass_matrix_.local_element(i);
1037
1038 /* skip diagonal */
1039 constexpr auto simd_length = VectorizedArray<Number>::size();
1040 const unsigned int *js = sparsity_pattern_simd_.columns(i);
1041 for (unsigned int col_idx = 1; col_idx < row_length; ++col_idx) {
1042 const auto j = *(i < n_locally_internal_ ? js + col_idx * simd_length
1043 : js + col_idx);
1044 Assert(j < n_locally_relevant_, dealii::ExcInternalError());
1045
1046 const auto m_ij = mass_matrix_.get_entry(i, col_idx);
1047 if (discretization_->have_discontinuous_ansatz()) {
1048 // Interfacial coupling terms are present in the stencil but zero
1049 // in the mass matrix
1050 Assert(std::abs(m_ij) > -1.e-12, dealii::ExcInternalError());
1051 } else {
1052 Assert(std::abs(m_ij) > 1.e-12, dealii::ExcInternalError());
1053 }
1054 sum += m_ij;
1055
1056 const auto m_ji = mass_matrix_.get_transposed_entry(i, col_idx);
1057 if (std::abs(m_ij - m_ji) >= 1.e-12) {
1058 // The m_ij matrix is not symmetric
1059 std::stringstream ss;
1060 ss << "m_ij matrix is not symmetric: " << m_ij << " <-> " << m_ji;
1061 Assert(false, dealii::ExcMessage(ss.str()));
1062 }
1063 }
1064
1065 Assert(std::abs(sum) < 1.e-12, dealii::ExcInternalError());
1066 }
1067
1068 /*
1069 * Verify that the cij_matrix_ object is consistent:
1070 */
1071
1072 for (unsigned int i = 0; i < n_locally_owned_; ++i) {
1073 /* Skip constrained degrees of freedom: */
1074 const unsigned int row_length = sparsity_pattern_simd_.row_length(i);
1075 if (row_length == 1)
1076 continue;
1077
1078 auto sum = cij_matrix_.get_tensor(i, 0);
1079
1080 /* skip diagonal */
1081 constexpr auto simd_length = VectorizedArray<Number>::size();
1082 const unsigned int *js = sparsity_pattern_simd_.columns(i);
1083 for (unsigned int col_idx = 1; col_idx < row_length; ++col_idx) {
1084 const auto j = *(i < n_locally_internal_ ? js + col_idx * simd_length
1085 : js + col_idx);
1086 Assert(j < n_locally_relevant_, dealii::ExcInternalError());
1087
1088 const auto c_ij = cij_matrix_.get_tensor(i, col_idx);
1089 Assert(c_ij.norm() > 1.e-12, dealii::ExcInternalError());
1090 sum += c_ij;
1091
1092 const auto c_ji = cij_matrix_.get_transposed_tensor(i, col_idx);
1093 if ((c_ij + c_ji).norm() >= 1.e-12) {
1094 // The c_ij matrix is not symmetric, this can only happen if i
1095 // and j are both located on the boundary.
1096
1097 CouplingDescription coupling{i, col_idx, j};
1098 const auto it = std::find(coupling_boundary_pairs_.begin(),
1099 coupling_boundary_pairs_.end(),
1100 coupling);
1101 if (it == coupling_boundary_pairs_.end()) {
1102 std::stringstream ss;
1103 ss << "c_ij matrix is not anti-symmetric: " << c_ij << " <-> "
1104 << c_ji;
1105 Assert(false, dealii::ExcMessage(ss.str()));
1106 }
1107 }
1108 }
1109
1110 Assert(sum.norm() < 1.e-12, dealii::ExcInternalError());
1111 }
1112#endif
1113 }
1114
1115
1116 template <int dim, typename Number>
1117 void OfflineData<dim, Number>::create_multigrid_data()
1118 {
1119#ifdef DEBUG_OUTPUT
1120 std::cout << "OfflineData<dim, Number>::create_multigrid_data()"
1121 << std::endl;
1122#endif
1123
1124 auto &dof_handler = *dof_handler_;
1125
1126 dof_handler.distribute_mg_dofs();
1127
1128 const auto n_levels = dof_handler.get_triangulation().n_global_levels();
1129
1130 AffineConstraints<float> level_constraints;
1131 // TODO not yet thread-parallel and without periodicity
1132
1133 level_boundary_map_.resize(n_levels);
1134 level_lumped_mass_matrix_.resize(n_levels);
1135
1136 for (unsigned int level = 0; level < n_levels; ++level) {
1137 /* Assemble lumped mass matrix vector: */
1138
1139 IndexSet relevant_dofs;
1140 dealii::DoFTools::extract_locally_relevant_level_dofs(
1141 dof_handler, level, relevant_dofs);
1142 const auto partitioner = std::make_shared<Utilities::MPI::Partitioner>(
1143 dof_handler.locally_owned_mg_dofs(level),
1144 relevant_dofs,
1145 mpi_ensemble_.ensemble_communicator());
1146 level_lumped_mass_matrix_[level].reinit(partitioner);
1147 std::vector<types::global_dof_index> dof_indices(
1148 dof_handler.get_fe().dofs_per_cell);
1149 dealii::Vector<Number> mass_values(dof_handler.get_fe().dofs_per_cell);
1150 FEValues<dim> fe_values(discretization_->mapping(),
1151 discretization_->finite_element(),
1152 discretization_->quadrature(),
1153 update_values | update_JxW_values);
1154 for (const auto &cell : dof_handler.cell_iterators_on_level(level))
1155 // TODO for assembly with dealii::SparseMatrix and local
1156 // numbering this probably has to read !cell->is_artificial()
1157 if (cell->is_locally_owned_on_level()) {
1158 fe_values.reinit(cell);
1159 for (unsigned int i = 0; i < mass_values.size(); ++i) {
1160 double sum = 0;
1161 for (unsigned int q = 0; q < fe_values.n_quadrature_points; ++q)
1162 sum += fe_values.shape_value(i, q) * fe_values.JxW(q);
1163 mass_values(i) = sum;
1164 }
1165 cell->get_mg_dof_indices(dof_indices);
1166 level_constraints.distribute_local_to_global(
1167 mass_values, dof_indices, level_lumped_mass_matrix_[level]);
1168 }
1169 level_lumped_mass_matrix_[level].compress(VectorOperation::add);
1170
1171 /* Populate boundary map: */
1172
1173 level_boundary_map_[level] = construct_boundary_map(
1174 dof_handler.begin_mg(level), dof_handler.end_mg(level), *partitioner);
1175 }
1176 }
1177
1178
1179 template <int dim, typename Number>
1180 template <typename ITERATOR1, typename ITERATOR2>
1182 const ITERATOR1 &begin,
1183 const ITERATOR2 &end,
1184 const Utilities::MPI::Partitioner &partitioner) const -> BoundaryMap
1185 {
1186#ifdef DEBUG_OUTPUT
1187 std::cout << "OfflineData<dim, Number>::construct_boundary_map()"
1188 << std::endl;
1189#endif
1190
1191 /*
1192 * Create a temporary multimap with the (local) dof index as key:
1193 */
1194
1195 using BoundaryData = std::tuple<dealii::Tensor<1, dim, Number> /*normal*/,
1196 Number /*normal mass*/,
1197 Number /*boundary mass*/,
1198 dealii::types::boundary_id /*id*/,
1199 dealii::Point<dim>> /*position*/;
1200 std::multimap<unsigned int, BoundaryData> preliminary_map;
1201
1202 std::vector<dealii::types::global_dof_index> local_dof_indices;
1203
1204 const dealii::QGauss<dim - 1> face_quadrature(3);
1205 dealii::FEFaceValues<dim> fe_face_values(discretization_->mapping(),
1206 discretization_->finite_element(),
1207 face_quadrature,
1208 dealii::update_normal_vectors |
1209 dealii::update_values |
1210 dealii::update_JxW_values);
1211
1212 const unsigned int dofs_per_cell =
1213 discretization_->finite_element().dofs_per_cell;
1214
1215 const auto support_points =
1216 discretization_->finite_element().get_unit_support_points();
1217
1218 for (auto cell = begin; cell != end; ++cell) {
1219
1220 /*
1221 * Workaround: Make sure to iterate over the entire locally relevant
1222 * set so that we compute normals between cells with differing owners
1223 * correctly. This strategy works for 2D but will fail in 3D with
1224 * locally refined meshes and hanging nodes situated at the boundary.
1225 */
1226 if ((cell->is_active() && cell->is_artificial()) ||
1227 (!cell->is_active() && cell->is_artificial_on_level()))
1228 continue;
1229
1230 local_dof_indices.resize(dofs_per_cell);
1231 cell->get_active_or_mg_dof_indices(local_dof_indices);
1232
1233 for (auto f : cell->face_indices()) {
1234 const auto face = cell->face(f);
1235 const auto id = face->boundary_id();
1236
1237 if (!face->at_boundary())
1238 continue;
1239
1240 /*
1241 * Skip periodic boundary faces. For our algorithm these are
1242 * interior degrees of freedom (if not simultaneously located at
1243 * another boundary as well).
1244 */
1245 if (id == Boundary::periodic)
1246 continue;
1247
1248 fe_face_values.reinit(cell, f);
1249
1250 for (unsigned int j : fe_face_values.dof_indices()) {
1251 if (!discretization_->finite_element().has_support_on_face(j, f))
1252 continue;
1253
1254 Number boundary_mass = 0.;
1255 dealii::Tensor<1, dim, Number> normal;
1256
1257 for (unsigned int q : fe_face_values.quadrature_point_indices()) {
1258 const auto JxW = fe_face_values.JxW(q);
1259 const auto phi_i = fe_face_values.shape_value(j, q);
1260
1261 boundary_mass += phi_i * JxW;
1262 normal += phi_i * fe_face_values.normal_vector(q) * JxW;
1263 }
1264
1265 const auto global_index = local_dof_indices[j];
1266 const auto index = partitioner.global_to_local(global_index);
1267
1268 /* Skip nonlocal degrees of freedom: */
1269 if (index >= n_locally_owned_)
1270 continue;
1271
1272 /* Skip constrained degrees of freedom: */
1273 const unsigned int row_length =
1274 sparsity_pattern_simd_.row_length(index);
1275 if (row_length == 1)
1276 continue;
1277
1278 Point<dim> position =
1279 discretization_->mapping().transform_unit_to_real_cell(
1280 cell, support_points[j]);
1281
1282 /*
1283 * Temporarily insert a (wrong) boundary mass value for the
1284 * normal mass. We'll fix this later.
1285 */
1286 preliminary_map.insert(
1287 {index, {normal, boundary_mass, boundary_mass, id, position}});
1288 } /* j */
1289 } /* f */
1290 } /* cell */
1291
1292 /*
1293 * Filter boundary map:
1294 *
1295 * At this point we have collected multiple cell contributions for each
1296 * boundary degree of freedom. We now merge all entries that have the
1297 * same boundary id and whose normals describe an acute angle of about
1298 * 60 degrees or less.
1299 *
1300 * FIXME: is this robust in 3D?
1301 */
1302
1303 std::multimap<unsigned int, BoundaryData> filtered_map;
1304 std::set<dealii::types::global_dof_index> boundary_dofs;
1305 for (auto entry : preliminary_map) {
1306 bool inserted = false;
1307 const auto range = filtered_map.equal_range(entry.first);
1308 for (auto it = range.first; it != range.second; ++it) {
1309 auto &[new_normal,
1310 new_normal_mass,
1311 new_boundary_mass,
1312 new_id,
1313 new_point] = entry.second;
1314 auto &[normal, normal_mass, boundary_mass, id, point] = it->second;
1315
1316 if (id != new_id)
1317 continue;
1318
1319 Assert(point.distance(new_point) < 1.0e-14, dealii::ExcInternalError());
1320
1321 if (normal * new_normal / normal.norm() / new_normal.norm() > 0.50) {
1322 /*
1323 * Both normals describe an acute angle of 85 degrees or less.
1324 * Merge the entries and continue.
1325 */
1326 normal += new_normal;
1327 boundary_mass += new_boundary_mass;
1328 inserted = true;
1329 continue;
1330
1331 } else if constexpr (dim == 2) {
1332 /*
1333 * Workaround for 2D: When enforcing slip boundary conditions
1334 * with two noncollinear vectors the resulting momentum must be
1335 * 0. But the normals don't necessarily describe an orthonormal
1336 * basis and we cannot use orthogonal projection. Therefore,
1337 * simply set the boundary type to no slip:
1338 */
1339 if (new_id == Boundary::slip) {
1340 Assert(id == Boundary::slip, dealii::ExcInternalError());
1341 new_id = Boundary::no_slip;
1342 id = Boundary::no_slip;
1343 }
1344 }
1345 }
1346 if (!inserted)
1347 filtered_map.insert(entry);
1348 }
1349
1350 /*
1351 * Normalize all normal vectors and create final boundary_map:
1352 */
1353
1354 BoundaryMap boundary_map;
1355 std::transform(
1356 std::begin(filtered_map),
1357 std::end(filtered_map),
1358 std::back_inserter(boundary_map),
1359 [&](const auto &it) -> BoundaryDescription { //
1360 auto index = it.first;
1361 const auto &[normal, normal_mass, boundary_mass, id, point] =
1362 it.second;
1363
1364 const auto new_normal_mass =
1365 normal.norm() + std::numeric_limits<Number>::epsilon();
1366 const auto new_normal = normal / new_normal_mass;
1367
1368 return {index, new_normal, new_normal_mass, boundary_mass, id, point};
1369 });
1370
1371 return boundary_map;
1372 }
1373
1374
1375 template <int dim, typename Number>
1376 template <typename ITERATOR1, typename ITERATOR2>
1378 const ITERATOR1 &begin,
1379 const ITERATOR2 &end,
1380 const Utilities::MPI::Partitioner &partitioner) const
1381 -> CouplingBoundaryPairs
1382 {
1383#ifdef DEBUG_OUTPUT
1384 std::cout << "OfflineData<dim, Number>::collect_coupling_boundary_pairs()"
1385 << std::endl;
1386#endif
1387
1388 /*
1389 * First, collect *all* locally relevant degrees of freedom that are
1390 * located on a (non periodic) boundary. We also collect constrained
1391 * degrees of freedom for the time being (and filter afterwards).
1392 */
1393
1394 std::set<unsigned int> locally_relevant_boundary_indices;
1395
1396 std::vector<dealii::types::global_dof_index> local_dof_indices;
1397
1398 const auto &finite_element = discretization_->finite_element();
1399 const unsigned int dofs_per_cell = finite_element.dofs_per_cell;
1400 local_dof_indices.resize(dofs_per_cell);
1401
1402 for (auto cell = begin; cell != end; ++cell) {
1403
1404 /* Make sure to iterate over the entire locally relevant set: */
1405 if (cell->is_artificial())
1406 continue;
1407
1408 cell->get_active_or_mg_dof_indices(local_dof_indices);
1409
1410 for (auto f : cell->face_indices()) {
1411 const auto face = cell->face(f);
1412 const auto id = face->boundary_id();
1413
1414 if (!face->at_boundary())
1415 continue;
1416
1417 /* Skip periodic boundary faces; see above. */
1418 if (id == Boundary::periodic)
1419 continue;
1420
1421 for (unsigned int j = 0; j < dofs_per_cell; ++j) {
1422
1423 if (!finite_element.has_support_on_face(j, f))
1424 continue;
1425
1426 const auto global_index = local_dof_indices[j];
1427 const auto index = partitioner.global_to_local(global_index);
1428
1429 /* Skip irrelevant degrees of freedom: */
1430 if (index >= n_locally_relevant_)
1431 continue;
1432
1433 locally_relevant_boundary_indices.insert(index);
1434 } /* j */
1435 } /* f */
1436 } /* cell */
1437
1438 /*
1439 * Now, collect all coupling boundary pairs:
1440 */
1441
1442 CouplingBoundaryPairs result;
1443
1444 for (const auto i : locally_relevant_boundary_indices) {
1445
1446 /* Only record pairs with a left index that is locally owned: */
1447 if (i >= n_locally_owned_)
1448 continue;
1449
1450 const unsigned int row_length = sparsity_pattern_simd_.row_length(i);
1451
1452 /* Skip all constrained degrees of freedom: */
1453 if (row_length == 1)
1454 continue;
1455
1456 const unsigned int *js = sparsity_pattern_simd_.columns(i);
1457 constexpr auto simd_length = VectorizedArray<Number>::size();
1458 /* skip diagonal: */
1459 for (unsigned int col_idx = 1; col_idx < row_length; ++col_idx) {
1460 const auto j = *(i < n_locally_internal_ ? js + col_idx * simd_length
1461 : js + col_idx);
1462
1463 if (locally_relevant_boundary_indices.count(j) != 0) {
1464 result.push_back({i, col_idx, j});
1465 }
1466 }
1467 }
1468
1469 return result;
1470 }
1471
1472} /* 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)