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>
28#ifdef FORCE_DEAL_II_SPARSE_MATRIX
29#undef DEAL_II_WITH_TRILINOS
34 using namespace dealii;
37 template <
int dim,
typename Number>
39 const MPI_Comm &mpi_communicator,
41 const std::string &subsection )
42 : ParameterAcceptor(subsection)
43 , discretization_(&discretization)
44 , mpi_communicator_(mpi_communicator)
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.");
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.");
62 template <
int dim,
typename Number>
71 auto &dof_handler = *dof_handler_;
72 const IndexSet &locally_owned = dof_handler.locally_owned_dofs();
74 IndexSet locally_relevant;
75 DoFTools::extract_locally_relevant_dofs(dof_handler, locally_relevant);
77 affine_constraints_.reinit(locally_relevant);
78 DoFTools::make_hanging_node_constraints(dof_handler, affine_constraints_);
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"));
91 const auto &periodic_faces =
92 discretization_->triangulation().get_periodic_face_map();
94 for (
const auto &[left, value] : periodic_faces) {
95 const auto &[right, orientation] = value;
97 typename DoFHandler<dim>::cell_iterator dof_cell_left(
98 &left.first->get_triangulation(),
103 typename DoFHandler<dim>::cell_iterator dof_cell_right(
104 &right.first->get_triangulation(),
105 right.first->level(),
106 right.first->index(),
109 if constexpr (std::is_same<Number, double>::value) {
110 DoFTools::make_periodicity_constraints(
111 dof_cell_left->face(left.second),
112 dof_cell_right->face(right.second),
115#
if DEAL_II_VERSION_GTE(9, 6, 0)
123 AssertThrow(
false, dealii::ExcNotImplemented());
128 affine_constraints_.close();
133 const std::vector<IndexSet> &locally_owned_dofs =
134 Utilities::MPI::all_gather(mpi_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(locally_owned_dofs,
146 sparsity_pattern_.reinit(
147 dof_handler.n_dofs(), dof_handler.n_dofs(), locally_relevant);
149 if (discretization_->have_discontinuous_ansatz()) {
154 dof_handler, sparsity_pattern_, affine_constraints_,
false);
159#ifdef DEAL_II_WITH_TRILINOS
160 DoFTools::make_sparsity_pattern(
161 dof_handler, sparsity_pattern_, affine_constraints_,
false);
170 dof_handler, sparsity_pattern_, affine_constraints_,
false);
180 SparsityTools::distribute_sparsity_pattern(
181 sparsity_pattern_, locally_owned, mpi_communicator_, locally_relevant);
185 template <
int dim,
typename Number>
186 void OfflineData<dim, Number>::setup(
const unsigned int problem_dimension,
187 const unsigned int n_precomputed_values)
190 std::cout <<
"OfflineData<dim, Number>::setup()" << std::endl;
197 const auto &triangulation = discretization_->triangulation();
199 dof_handler_ = std::make_unique<dealii::DoFHandler<dim>>(triangulation);
200 auto &dof_handler = *dof_handler_;
202 dof_handler.distribute_dofs(discretization_->finite_element());
204 n_locally_owned_ = dof_handler.locally_owned_dofs().n_elements();
210 DoFRenumbering::Cuthill_McKee(dof_handler);
221 dof_handler, mpi_communicator_, n_locally_owned_, 1);
230 create_constraints_and_sparsity_pattern();
232 dof_handler, sparsity_pattern_, VectorizedArray<Number>::size());
249 VectorizedArray<Number>::size());
255 const auto consistent_stride_range [[maybe_unused]] = [&]() {
256 constexpr auto group_size = VectorizedArray<Number>::size();
257 const IndexSet &locally_owned = dof_handler.locally_owned_dofs();
258 const auto offset = n_locally_owned_ != 0 ? *locally_owned.begin() : 0;
260 unsigned int group_row_length = 0;
262 for (; i < n_locally_internal_; ++i) {
263 if (i % group_size == 0) {
264 group_row_length = sparsity_pattern_.row_length(offset + i);
266 if (group_row_length != sparsity_pattern_.row_length(offset + i)) {
271 return i / group_size * group_size;
277 const auto mpi_allreduce_logical_or = [&](
const bool local_value) {
278 std::function<bool(
const bool &,
const bool &)> comparator =
279 [](
const bool &left,
const bool &right) ->
bool {
280 return left || right;
282 return Utilities::MPI::all_reduce(
283 local_value, mpi_communicator_, comparator);
290 create_constraints_and_sparsity_pattern();
300#if DEAL_II_VERSION_GTE(9, 5, 0)
301 if (mpi_allreduce_logical_or(affine_constraints_.n_constraints() > 0)) {
302 if (mpi_allreduce_logical_or(
303 consistent_stride_range() != n_locally_internal_)) {
313 VectorizedArray<Number>::size());
314 create_constraints_and_sparsity_pattern();
315 n_locally_internal_ = consistent_stride_range();
325 Assert(consistent_stride_range() == n_locally_internal_,
326 dealii::ExcInternalError());
332 const IndexSet &locally_owned = dof_handler.locally_owned_dofs();
333 Assert(n_locally_owned_ == locally_owned.n_elements(),
334 dealii::ExcInternalError());
336 IndexSet locally_relevant;
337 DoFTools::extract_locally_relevant_dofs(dof_handler, locally_relevant);
340 IndexSet additional_dofs(dof_handler.n_dofs());
341 for (
auto &entry : sparsity_pattern_)
342 if (!locally_relevant.is_element(entry.column())) {
343 Assert(locally_owned.is_element(entry.row()), ExcInternalError());
344 additional_dofs.add_index(entry.column());
346 additional_dofs.compress();
347 locally_relevant.add_indices(additional_dofs);
348 locally_relevant.compress();
351 n_locally_relevant_ = locally_relevant.n_elements();
353 scalar_partitioner_ = std::make_shared<dealii::Utilities::MPI::Partitioner>(
354 locally_owned, locally_relevant, mpi_communicator_);
357 scalar_partitioner_, problem_dimension);
360 scalar_partitioner_, n_precomputed_values);
368 if (mpi_allreduce_logical_or(affine_constraints_.n_constraints() > 0)) {
372 n_export_indices_ = 0;
373 for (
const auto &it : scalar_partitioner_->import_indices())
374 if (it.second <= n_locally_internal_)
375 n_export_indices_ = std::max(n_export_indices_, it.second);
377 constexpr auto simd_length = VectorizedArray<Number>::size();
379 (n_export_indices_ + simd_length - 1) / simd_length * simd_length;
384 unsigned int control = 0;
385 for (
const auto &it : scalar_partitioner_->import_indices())
386 if (it.second <= n_locally_internal_)
387 control = std::max(control, it.second);
389 Assert(control <= n_export_indices_, ExcInternalError());
390 Assert(n_export_indices_ <= n_locally_internal_, ExcInternalError());
399 sparsity_pattern_simd_.reinit(
400 n_locally_internal_, sparsity_pattern_, scalar_partitioner_);
406 mass_matrix_.reinit(sparsity_pattern_simd_);
407 if (discretization_->have_discontinuous_ansatz())
408 mass_matrix_inverse_.reinit(sparsity_pattern_simd_);
410 lumped_mass_matrix_.reinit(scalar_partitioner_);
411 lumped_mass_matrix_inverse_.reinit(scalar_partitioner_);
413 cij_matrix_.reinit(sparsity_pattern_simd_);
414 if (discretization_->have_discontinuous_ansatz())
415 incidence_matrix_.reinit(sparsity_pattern_simd_);
419 template <
int dim,
typename Number>
420 void OfflineData<dim, Number>::assemble()
423 std::cout <<
"OfflineData<dim, Number>::assemble()" << std::endl;
426 auto &dof_handler = *dof_handler_;
428 measure_of_omega_ = 0.;
430#ifdef DEAL_II_WITH_TRILINOS
433 AffineConstraints<double> affine_constraints_assembly;
435 affine_constraints_assembly.reinit(affine_constraints_.get_local_lines());
436 for (
auto line : affine_constraints_.get_lines()) {
437 affine_constraints_assembly.add_line(line.index);
438 for (
auto entry : line.entries)
439 affine_constraints_assembly.add_entry(
440 line.index, entry.first, entry.second);
441 affine_constraints_assembly.set_inhomogeneity(line.index,
444 affine_constraints_assembly.close();
446 const IndexSet &locally_owned = dof_handler.locally_owned_dofs();
447 TrilinosWrappers::SparsityPattern trilinos_sparsity_pattern;
448 trilinos_sparsity_pattern.reinit(
449 locally_owned, sparsity_pattern_, mpi_communicator_);
451 TrilinosWrappers::SparseMatrix mass_matrix_tmp;
452 TrilinosWrappers::SparseMatrix mass_matrix_inverse_tmp;
453 if (discretization_->have_discontinuous_ansatz())
454 mass_matrix_inverse_tmp.reinit(trilinos_sparsity_pattern);
456 std::array<TrilinosWrappers::SparseMatrix, dim> cij_matrix_tmp;
458 mass_matrix_tmp.reinit(trilinos_sparsity_pattern);
459 for (
auto &matrix : cij_matrix_tmp)
460 matrix.reinit(trilinos_sparsity_pattern);
465 AffineConstraints<Number> affine_constraints_assembly;
466 affine_constraints_assembly.copy_from(affine_constraints_);
469 SparsityPattern sparsity_pattern_assembly;
471 DynamicSparsityPattern dsp(n_locally_relevant_, n_locally_relevant_);
472 for (
const auto &entry : sparsity_pattern_) {
473 const auto i = scalar_partitioner_->global_to_local(entry.row());
474 const auto j = scalar_partitioner_->global_to_local(entry.column());
477 sparsity_pattern_assembly.copy_from(dsp);
480 dealii::SparseMatrix<Number> mass_matrix_tmp;
481 dealii::SparseMatrix<Number> mass_matrix_inverse_tmp;
482 if (discretization_->have_discontinuous_ansatz())
483 mass_matrix_inverse_tmp.reinit(sparsity_pattern_assembly);
485 std::array<dealii::SparseMatrix<Number>, dim> cij_matrix_tmp;
487 mass_matrix_tmp.reinit(sparsity_pattern_assembly);
488 for (
auto &matrix : cij_matrix_tmp)
489 matrix.reinit(sparsity_pattern_assembly);
492 const unsigned int dofs_per_cell =
493 discretization_->finite_element().dofs_per_cell;
495 const unsigned int n_q_points = discretization_->quadrature().size();
496 const unsigned int n_face_q_points =
497 discretization_->face_quadrature().size();
498 const unsigned int n_face_q_points_nodal =
499 discretization_->face_nodal_quadrature().size();
506 const auto local_assemble_system = [&](
const auto &cell,
511 auto &is_locally_owned = copy.is_locally_owned_;
512 auto &local_dof_indices = copy.local_dof_indices_;
513 auto &neighbor_local_dof_indices = copy.neighbor_local_dof_indices_;
515 auto &cell_mass_matrix = copy.cell_mass_matrix_;
516 auto &cell_mass_matrix_inverse = copy.cell_mass_matrix_inverse_;
517 auto &cell_cij_matrix = copy.cell_cij_matrix_;
518 auto &interface_cij_matrix = copy.interface_cij_matrix_;
519 auto &cell_measure = copy.cell_measure_;
521 auto &fe_values = scratch.fe_values_;
522 auto &fe_face_values = scratch.fe_face_values_;
523 auto &fe_neighbor_face_values = scratch.fe_neighbor_face_values_;
525#ifdef DEAL_II_WITH_TRILINOS
526 is_locally_owned = cell->is_locally_owned();
534 is_locally_owned = !cell->is_artificial();
536 if (!is_locally_owned)
539 cell_mass_matrix.reinit(dofs_per_cell, dofs_per_cell);
540 for (
auto &matrix : cell_cij_matrix)
541 matrix.reinit(dofs_per_cell, dofs_per_cell);
542 if (discretization_->have_discontinuous_ansatz()) {
543 cell_mass_matrix_inverse.reinit(dofs_per_cell, dofs_per_cell);
544 for (
auto &it : interface_cij_matrix)
545 for (
auto &matrix : it)
546 matrix.reinit(dofs_per_cell, dofs_per_cell);
549 fe_values.reinit(cell);
551 local_dof_indices.resize(dofs_per_cell);
552 cell->get_dof_indices(local_dof_indices);
555 cell_mass_matrix = 0.;
556 for (
auto &matrix : cell_cij_matrix)
558 if (discretization_->have_discontinuous_ansatz()) {
559 cell_mass_matrix_inverse = 0.;
560 for (
auto &it : interface_cij_matrix)
561 for (
auto &matrix : it)
566 for (
unsigned int q_point = 0; q_point < n_q_points; ++q_point) {
567 const auto JxW = fe_values.JxW(q_point);
569 if (cell->is_locally_owned())
570 cell_measure += Number(JxW);
572 for (
unsigned int j = 0; j < dofs_per_cell; ++j) {
573 const auto value_JxW = fe_values.shape_value(j, q_point) * JxW;
574 const auto grad_JxW = fe_values.shape_grad(j, q_point) * JxW;
576 for (
unsigned int i = 0; i < dofs_per_cell; ++i) {
578 const auto value = fe_values.shape_value(i, q_point);
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]);
593 if (!discretization_->have_discontinuous_ansatz())
596 for (
const auto f_index : cell->face_indices()) {
597 const auto &face = cell->face(f_index);
600 const bool has_neighbor =
601 !face->at_boundary() || cell->has_periodic_neighbor(f_index);
605 neighbor_local_dof_indices[f_index].resize(0);
610 const auto neighbor_cell = cell->neighbor_or_periodic_neighbor(f_index);
611 if (neighbor_cell->is_artificial()) {
614 neighbor_local_dof_indices[f_index].resize(0);
618 fe_face_values.reinit(cell, f_index);
622 for (
unsigned int q = 0; q < n_face_q_points; ++q) {
623 const auto JxW = fe_face_values.JxW(q);
624 const auto &normal = fe_face_values.get_normal_vectors()[q];
626 for (
unsigned int j = 0; j < dofs_per_cell; ++j) {
627 const auto value_JxW = fe_face_values.shape_value(j, q) * JxW;
629 for (
unsigned int i = 0; i < dofs_per_cell; ++i) {
630 const auto value = fe_face_values.shape_value(i, q);
632 for (
unsigned int d = 0; d < dim; ++d)
633 cell_cij_matrix[d](i, j) -=
634 Number(0.5 * normal[d] * value * value_JxW);
641 const unsigned int f_index_neighbor =
642 cell->has_periodic_neighbor(f_index)
643 ? cell->periodic_neighbor_of_periodic_neighbor(f_index)
644 : cell->neighbor_of_neighbor(f_index);
646 neighbor_local_dof_indices[f_index].resize(dofs_per_cell);
647 neighbor_cell->get_dof_indices(neighbor_local_dof_indices[f_index]);
649 fe_neighbor_face_values.reinit(neighbor_cell, f_index_neighbor);
651 for (
unsigned int q = 0; q < n_face_q_points; ++q) {
652 const auto JxW = fe_face_values.JxW(q);
653 const auto &normal = fe_face_values.get_normal_vectors()[q];
656 for (
unsigned int j = 0; j < dofs_per_cell; ++j) {
657 const auto value_JxW =
658 fe_neighbor_face_values.shape_value(j, q) * JxW;
660 for (
unsigned int i = 0; i < dofs_per_cell; ++i) {
661 const auto value = fe_face_values.shape_value(i, q);
663 for (
unsigned int d = 0; d < dim; ++d)
664 interface_cij_matrix[f_index][d](i, j) +=
665 Number(0.5 * normal[d] * value * value_JxW);
675 if (discretization_->have_discontinuous_ansatz()) {
677 cell_mass_matrix_inverse.invert(cell_mass_matrix);
681 const auto copy_local_to_global = [&](
const auto ©) {
682 const auto &is_locally_owned = copy.is_locally_owned_;
683#ifdef DEAL_II_WITH_TRILINOS
684 const auto &local_dof_indices = copy.local_dof_indices_;
685 const auto &neighbor_local_dof_indices = copy.neighbor_local_dof_indices_;
692 auto local_dof_indices = copy.local_dof_indices_;
693 auto neighbor_local_dof_indices = copy.neighbor_local_dof_indices_;
695 const auto &cell_mass_matrix = copy.cell_mass_matrix_;
696 const auto &cell_mass_matrix_inverse = copy.cell_mass_matrix_inverse_;
697 const auto &cell_cij_matrix = copy.cell_cij_matrix_;
698 const auto &interface_cij_matrix = copy.interface_cij_matrix_;
699 const auto &cell_measure = copy.cell_measure_;
701 if (!is_locally_owned)
704#ifndef DEAL_II_WITH_TRILINOS
706 for (
auto &indices : neighbor_local_dof_indices)
710 affine_constraints_assembly.distribute_local_to_global(
711 cell_mass_matrix, local_dof_indices, mass_matrix_tmp);
713 for (
int k = 0; k < dim; ++k) {
714 affine_constraints_assembly.distribute_local_to_global(
715 cell_cij_matrix[k], local_dof_indices, cij_matrix_tmp[k]);
717 for (
unsigned int f_index = 0; f_index < copy.n_faces; ++f_index) {
718 if (neighbor_local_dof_indices[f_index].size() != 0) {
719 affine_constraints_assembly.distribute_local_to_global(
720 interface_cij_matrix[f_index][k],
722 neighbor_local_dof_indices[f_index],
728 if (discretization_->have_discontinuous_ansatz())
729 affine_constraints_assembly.distribute_local_to_global(
730 cell_mass_matrix_inverse,
732 mass_matrix_inverse_tmp);
734 measure_of_omega_ += cell_measure;
737 WorkStream::run(dof_handler.begin_active(),
739 local_assemble_system,
740 copy_local_to_global,
741 AssemblyScratchData<dim>(*discretization_),
742#ifdef DEAL_II_WITH_TRILINOS
743 AssemblyCopyData<dim, double>());
745 AssemblyCopyData<dim, Number>());
748#ifdef DEAL_II_WITH_TRILINOS
750 for (
auto &it : cij_matrix_tmp)
753 mass_matrix_.read_in(mass_matrix_tmp,
false);
754 if (discretization_->have_discontinuous_ansatz())
755 mass_matrix_inverse_.read_in(mass_matrix_inverse_tmp,
false);
756 cij_matrix_.read_in(cij_matrix_tmp,
false);
758 mass_matrix_.read_in(mass_matrix_tmp,
true);
759 if (discretization_->have_discontinuous_ansatz())
760 mass_matrix_inverse_.read_in(mass_matrix_inverse_tmp,
true);
761 cij_matrix_.read_in(cij_matrix_tmp,
true);
764 mass_matrix_.update_ghost_rows();
765 if (discretization_->have_discontinuous_ansatz())
766 mass_matrix_inverse_.update_ghost_rows();
767 cij_matrix_.update_ghost_rows();
770 Utilities::MPI::sum(measure_of_omega_, mpi_communicator_);
777#ifdef DEAL_II_WITH_TRILINOS
780 affine_constraints_assembly.set_zero(one);
782 ScalarVector local_lumped_mass_matrix(scalar_partitioner_);
783 mass_matrix_tmp.vmult(local_lumped_mass_matrix, one);
786 for (
unsigned int i = 0; i < scalar_partitioner_->locally_owned_size();
788 lumped_mass_matrix_.local_element(i) =
789 local_lumped_mass_matrix.local_element(i);
790 lumped_mass_matrix_inverse_.local_element(i) =
791 1. / lumped_mass_matrix_.local_element(i);
793 lumped_mass_matrix_.update_ghost_values();
794 lumped_mass_matrix_inverse_.update_ghost_values();
798 dealii::Vector<Number> one(mass_matrix_tmp.m());
800 affine_constraints_assembly.set_zero(one);
802 dealii::Vector<Number> local_lumped_mass_matrix(mass_matrix_tmp.m());
803 mass_matrix_tmp.vmult(local_lumped_mass_matrix, one);
805 for (
unsigned int i = 0; i < scalar_partitioner_->locally_owned_size();
807 lumped_mass_matrix_.local_element(i) = local_lumped_mass_matrix(i);
808 lumped_mass_matrix_inverse_.local_element(i) =
809 1. / lumped_mass_matrix_.local_element(i);
811 lumped_mass_matrix_.update_ghost_values();
812 lumped_mass_matrix_inverse_.update_ghost_values();
820 if (discretization_->have_discontinuous_ansatz()) {
821#ifdef DEAL_II_WITH_TRILINOS
822 TrilinosWrappers::SparseMatrix incidence_matrix_tmp;
823 incidence_matrix_tmp.reinit(trilinos_sparsity_pattern);
825 dealii::SparseMatrix<Number> incidence_matrix_tmp;
826 incidence_matrix_tmp.reinit(sparsity_pattern_assembly);
830 const auto local_assemble_system = [&](
const auto &cell,
835 auto &is_locally_owned = copy.is_locally_owned_;
836 auto &local_dof_indices = copy.local_dof_indices_;
837 auto &neighbor_local_dof_indices = copy.neighbor_local_dof_indices_;
838 auto &interface_incidence_matrix = copy.interface_incidence_matrix_;
839 auto &fe_face_values_nodal = scratch.fe_face_values_nodal_;
840 auto &fe_neighbor_face_values_nodal =
841 scratch.fe_neighbor_face_values_nodal_;
843#ifdef DEAL_II_WITH_TRILINOS
844 is_locally_owned = cell->is_locally_owned();
846 is_locally_owned = !cell->is_artificial();
848 if (!is_locally_owned)
851 for (
auto &matrix : interface_incidence_matrix)
852 matrix.reinit(dofs_per_cell, dofs_per_cell);
854 local_dof_indices.resize(dofs_per_cell);
855 cell->get_dof_indices(local_dof_indices);
858 for (
auto &matrix : interface_incidence_matrix)
861 for (
const auto f_index : cell->face_indices()) {
862 const auto &face = cell->face(f_index);
865 const bool has_neighbor =
866 !face->at_boundary() || cell->has_periodic_neighbor(f_index);
870 neighbor_local_dof_indices[f_index].resize(0);
875 const auto neighbor_cell =
876 cell->neighbor_or_periodic_neighbor(f_index);
877 if (neighbor_cell->is_artificial()) {
880 neighbor_local_dof_indices[f_index].resize(0);
884 const unsigned int f_index_neighbor =
885 cell->has_periodic_neighbor(f_index)
886 ? cell->periodic_neighbor_of_periodic_neighbor(f_index)
887 : cell->neighbor_of_neighbor(f_index);
889 neighbor_local_dof_indices[f_index].resize(dofs_per_cell);
890 neighbor_cell->get_dof_indices(neighbor_local_dof_indices[f_index]);
892 fe_face_values_nodal.reinit(cell, f_index);
893 fe_neighbor_face_values_nodal.reinit(neighbor_cell, f_index_neighbor);
897 for (
unsigned int q = 0; q < n_face_q_points_nodal; ++q) {
899 for (
unsigned int j = 0; j < dofs_per_cell; ++j) {
900 const auto v_j = fe_neighbor_face_values_nodal.shape_value(j, q);
901 for (
unsigned int i = 0; i < dofs_per_cell; ++i) {
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();
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];
912 Number(0.5) * (m_i + m_j) / measure_of_omega_;
923 r_ij =
std::pow(hd_ij, incidence_relaxation_even_ / dim);
930 r_ij =
std::pow(hd_ij, incidence_relaxation_odd_ / dim);
933 interface_incidence_matrix[f_index](i, j) += r_ij;
941 const auto copy_local_to_global = [&](
const auto ©) {
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_;
953 auto local_dof_indices = copy.local_dof_indices_;
954 auto neighbor_local_dof_indices = copy.neighbor_local_dof_indices_;
956 const auto &interface_incidence_matrix =
957 copy.interface_incidence_matrix_;
959 if (!is_locally_owned)
962#ifndef DEAL_II_WITH_TRILINOS
964 for (
auto &indices : neighbor_local_dof_indices)
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],
973 neighbor_local_dof_indices[f_index],
974 incidence_matrix_tmp);
979 WorkStream::run(dof_handler.begin_active(),
981 local_assemble_system,
982 copy_local_to_global,
983 AssemblyScratchData<dim>(*discretization_),
984#ifdef DEAL_II_WITH_TRILINOS
985 AssemblyCopyData<dim, double>());
987 AssemblyCopyData<dim, Number>());
990#ifdef DEAL_II_WITH_TRILINOS
991 incidence_matrix_.read_in(incidence_matrix_tmp,
false);
993 incidence_matrix_.read_in(incidence_matrix_tmp,
true);
995 incidence_matrix_.update_ghost_rows();
1002 boundary_map_ = construct_boundary_map(
1003 dof_handler.begin_active(), dof_handler.end(), *scalar_partitioner_);
1005 coupling_boundary_pairs_ = collect_coupling_boundary_pairs(
1006 dof_handler.begin_active(), dof_handler.end(), *scalar_partitioner_);
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 = Utilities::MPI::sum(total_mass, mpi_communicator_);
1019 Assert(std::abs(measure_of_omega_ - total_mass) < 1.e-12,
1021 "Total mass differs from the measure of the domain."));
1027 for (
unsigned int i = 0; i < n_locally_owned_; ++i) {
1029 const unsigned int row_length = sparsity_pattern_simd_.row_length(i);
1030 if (row_length == 1)
1034 mass_matrix_.get_entry(i, 0) - lumped_mass_matrix_.local_element(i);
1037 constexpr auto simd_length = VectorizedArray<Number>::size();
1038 const unsigned int *js = sparsity_pattern_simd_.columns(i);
1039 for (
unsigned int col_idx = 1; col_idx < row_length; ++col_idx) {
1040 const auto j = *(i < n_locally_internal_ ? js + col_idx * simd_length
1042 Assert(j < n_locally_relevant_, dealii::ExcInternalError());
1044 const auto m_ij = mass_matrix_.get_entry(i, col_idx);
1045 Assert(std::abs(m_ij) > 1.e-12, dealii::ExcInternalError());
1048 const auto m_ji = mass_matrix_.get_transposed_entry(i, col_idx);
1049 if (std::abs(m_ij - m_ji) >= 1.e-12) {
1051 std::stringstream ss;
1052 ss <<
"m_ij matrix is not symmetric: " << m_ij <<
" <-> " << m_ji;
1053 Assert(
false, dealii::ExcMessage(ss.str()));
1057 Assert(std::abs(sum) < 1.e-12, dealii::ExcInternalError());
1064 for (
unsigned int i = 0; i < n_locally_owned_; ++i) {
1066 const unsigned int row_length = sparsity_pattern_simd_.row_length(i);
1067 if (row_length == 1)
1070 auto sum = cij_matrix_.get_tensor(i, 0);
1073 constexpr auto simd_length = VectorizedArray<Number>::size();
1074 const unsigned int *js = sparsity_pattern_simd_.columns(i);
1075 for (
unsigned int col_idx = 1; col_idx < row_length; ++col_idx) {
1076 const auto j = *(i < n_locally_internal_ ? js + col_idx * simd_length
1078 Assert(j < n_locally_relevant_, dealii::ExcInternalError());
1080 const auto c_ij = cij_matrix_.get_tensor(i, col_idx);
1081 Assert(c_ij.norm() > 1.e-12, dealii::ExcInternalError());
1084 const auto c_ji = cij_matrix_.get_transposed_tensor(i, col_idx);
1085 if ((c_ij + c_ji).norm() >= 1.e-12) {
1089 CouplingDescription coupling{i, col_idx, j};
1090 const auto it = std::find(coupling_boundary_pairs_.begin(),
1091 coupling_boundary_pairs_.end(),
1093 if (it == coupling_boundary_pairs_.end()) {
1094 std::stringstream ss;
1095 ss <<
"c_ij matrix is not anti-symmetric: " << c_ij <<
" <-> "
1097 Assert(
false, dealii::ExcMessage(ss.str()));
1102 Assert(sum.norm() < 1.e-12, dealii::ExcInternalError());
1108 template <
int dim,
typename Number>
1109 void OfflineData<dim, Number>::create_multigrid_data()
1112 std::cout <<
"OfflineData<dim, Number>::create_multigrid_data()"
1116 auto &dof_handler = *dof_handler_;
1118 dof_handler.distribute_mg_dofs();
1120 const auto n_levels = dof_handler.get_triangulation().n_global_levels();
1122 AffineConstraints<float> level_constraints;
1125 level_boundary_map_.resize(n_levels);
1126 level_lumped_mass_matrix_.resize(n_levels);
1128 for (
unsigned int level = 0; level < n_levels; ++level) {
1131 IndexSet relevant_dofs;
1132 dealii::DoFTools::extract_locally_relevant_level_dofs(
1133 dof_handler, level, relevant_dofs);
1134 const auto partitioner = std::make_shared<Utilities::MPI::Partitioner>(
1135 dof_handler.locally_owned_mg_dofs(level),
1137 lumped_mass_matrix_.get_mpi_communicator());
1138 level_lumped_mass_matrix_[level].reinit(partitioner);
1139 std::vector<types::global_dof_index> dof_indices(
1140 dof_handler.get_fe().dofs_per_cell);
1141 dealii::Vector<Number> mass_values(dof_handler.get_fe().dofs_per_cell);
1142 FEValues<dim> fe_values(discretization_->mapping(),
1143 discretization_->finite_element(),
1144 discretization_->quadrature(),
1145 update_values | update_JxW_values);
1146 for (
const auto &cell : dof_handler.cell_iterators_on_level(level))
1149 if (cell->is_locally_owned_on_level()) {
1150 fe_values.reinit(cell);
1151 for (
unsigned int i = 0; i < mass_values.size(); ++i) {
1153 for (
unsigned int q = 0; q < fe_values.n_quadrature_points; ++q)
1154 sum += fe_values.shape_value(i, q) * fe_values.JxW(q);
1155 mass_values(i) = sum;
1157 cell->get_mg_dof_indices(dof_indices);
1158 level_constraints.distribute_local_to_global(
1159 mass_values, dof_indices, level_lumped_mass_matrix_[level]);
1165 level_boundary_map_[level] = construct_boundary_map(
1166 dof_handler.begin_mg(level), dof_handler.end_mg(level), *partitioner);
1171 template <
int dim,
typename Number>
1172 template <
typename ITERATOR1,
typename ITERATOR2>
1174 const ITERATOR1 &begin,
1175 const ITERATOR2 &end,
1176 const Utilities::MPI::Partitioner &partitioner)
const -> BoundaryMap
1179 std::cout <<
"OfflineData<dim, Number>::construct_boundary_map()"
1187 using BoundaryData = std::tuple<dealii::Tensor<1, dim, Number> ,
1190 dealii::types::boundary_id ,
1191 dealii::Point<dim>> ;
1192 std::multimap<unsigned int, BoundaryData> preliminary_map;
1194 std::vector<dealii::types::global_dof_index> local_dof_indices;
1196 const dealii::QGauss<dim - 1> face_quadrature(3);
1197 dealii::FEFaceValues<dim> fe_face_values(discretization_->mapping(),
1198 discretization_->finite_element(),
1200 dealii::update_normal_vectors |
1201 dealii::update_values |
1202 dealii::update_JxW_values);
1204 const unsigned int dofs_per_cell =
1205 discretization_->finite_element().dofs_per_cell;
1207 const auto support_points =
1208 discretization_->finite_element().get_unit_support_points();
1210 for (
auto cell = begin; cell != end; ++cell) {
1218 if ((cell->is_active() && cell->is_artificial()) ||
1219 (!cell->is_active() && cell->is_artificial_on_level()))
1222 local_dof_indices.resize(dofs_per_cell);
1223 cell->get_active_or_mg_dof_indices(local_dof_indices);
1225 for (
auto f : GeometryInfo<dim>::face_indices()) {
1226 const auto face = cell->face(f);
1227 const auto id = face->boundary_id();
1229 if (!face->at_boundary())
1240 fe_face_values.reinit(cell, f);
1241 const unsigned int n_face_q_points = face_quadrature.size();
1243 for (
unsigned int j = 0; j < dofs_per_cell; ++j) {
1245 if (!discretization_->finite_element().has_support_on_face(j, f))
1248 Number boundary_mass = 0.;
1249 dealii::Tensor<1, dim, Number> normal;
1251 for (
unsigned int q = 0; q < n_face_q_points; ++q) {
1252 const auto JxW = fe_face_values.JxW(q);
1253 const auto phi_i = fe_face_values.shape_value(j, q);
1255 boundary_mass += phi_i * JxW;
1256 normal += phi_i * fe_face_values.normal_vector(q) * JxW;
1259 const auto global_index = local_dof_indices[j];
1260 const auto index = partitioner.global_to_local(global_index);
1263 if (index >= n_locally_owned_)
1267 const unsigned int row_length =
1268 sparsity_pattern_simd_.row_length(index);
1269 if (row_length == 1)
1272 Point<dim> position =
1273 discretization_->mapping().transform_unit_to_real_cell(
1274 cell, support_points[j]);
1280 preliminary_map.insert(
1281 {index, {normal, boundary_mass, boundary_mass, id, position}});
1297 std::multimap<unsigned int, BoundaryData> filtered_map;
1298 std::set<dealii::types::global_dof_index> boundary_dofs;
1299 for (
auto entry : preliminary_map) {
1300 bool inserted =
false;
1301 const auto range = filtered_map.equal_range(entry.first);
1302 for (
auto it = range.first; it != range.second; ++it) {
1307 new_point] = entry.second;
1308 auto &[normal, normal_mass, boundary_mass, id, point] = it->second;
1313 Assert(point.distance(new_point) < 1.0e-14, dealii::ExcInternalError());
1315 if (normal * new_normal / normal.norm() / new_normal.norm() > 0.50) {
1320 normal += new_normal;
1321 boundary_mass += new_boundary_mass;
1325 }
else if constexpr (dim == 2) {
1341 filtered_map.insert(entry);
1348 BoundaryMap boundary_map;
1350 std::begin(filtered_map),
1351 std::end(filtered_map),
1352 std::back_inserter(boundary_map),
1354 auto index = it.first;
1355 const auto &[normal, normal_mass, boundary_mass, id, point] =
1358 const auto new_normal_mass =
1359 normal.norm() + std::numeric_limits<Number>::epsilon();
1360 const auto new_normal = normal / new_normal_mass;
1362 return {index, new_normal, new_normal_mass, boundary_mass, id, point};
1365 return boundary_map;
1369 template <
int dim,
typename Number>
1370 template <
typename ITERATOR1,
typename ITERATOR2>
1372 const ITERATOR1 &begin,
1373 const ITERATOR2 &end,
1374 const Utilities::MPI::Partitioner &partitioner)
const
1375 -> CouplingBoundaryPairs
1378 std::cout <<
"OfflineData<dim, Number>::collect_coupling_boundary_pairs()"
1388 std::set<unsigned int> locally_relevant_boundary_indices;
1390 std::vector<dealii::types::global_dof_index> local_dof_indices;
1392 const unsigned int dofs_per_cell =
1393 discretization_->finite_element().dofs_per_cell;
1395 for (
auto cell = begin; cell != end; ++cell) {
1398 if (cell->is_artificial())
1401 local_dof_indices.resize(dofs_per_cell);
1402 cell->get_active_or_mg_dof_indices(local_dof_indices);
1404 for (
auto f : GeometryInfo<dim>::face_indices()) {
1405 const auto face = cell->face(f);
1406 const auto id = face->boundary_id();
1408 if (!face->at_boundary())
1415 for (
unsigned int j = 0; j < dofs_per_cell; ++j) {
1417 if (!discretization_->finite_element().has_support_on_face(j, f))
1420 const auto global_index = local_dof_indices[j];
1421 const auto index = partitioner.global_to_local(global_index);
1424 if (index >= n_locally_relevant_)
1427 locally_relevant_boundary_indices.insert(index);
1436 CouplingBoundaryPairs result;
1438 for (
const auto i : locally_relevant_boundary_indices) {
1441 if (i >= n_locally_owned_)
1444 const unsigned int row_length = sparsity_pattern_simd_.row_length(i);
1447 if (row_length == 1)
1450 const unsigned int *js = sparsity_pattern_simd_.columns(i);
1451 constexpr auto simd_length = VectorizedArray<Number>::size();
1453 for (
unsigned int col_idx = 1; col_idx < row_length; ++col_idx) {
1454 const auto j = *(i < n_locally_internal_ ? js + col_idx * simd_length
1457 if (locally_relevant_boundary_indices.count(j) != 0) {
1458 result.push_back({i, col_idx, j});
std::tuple< unsigned int, dealii::Tensor< 1, dim, Number >, Number, Number, dealii::types::boundary_id, dealii::Point< dim > > BoundaryDescription
OfflineData(const MPI_Comm &mpi_communicator, 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
DEAL_II_ALWAYS_INLINE FT add(const FT &flux_left_ij, const FT &flux_right_ij)