From 55fbe402dd3ebd23b65a6b342229a068617e4217 Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Mon, 19 Sep 2022 19:45:53 +0200 Subject: [PATCH] Add TensorProductMatrixCreator namespace --- include/deal.II/lac/tensor_product_matrix.h | 1 - .../numerics/tensor_product_matrix_creator.h | 375 ++++++++++++++++++ 2 files changed, 375 insertions(+), 1 deletion(-) create mode 100644 include/deal.II/numerics/tensor_product_matrix_creator.h diff --git a/include/deal.II/lac/tensor_product_matrix.h b/include/deal.II/lac/tensor_product_matrix.h index 0c167ea0c9..86a32dcb26 100644 --- a/include/deal.II/lac/tensor_product_matrix.h +++ b/include/deal.II/lac/tensor_product_matrix.h @@ -472,7 +472,6 @@ private: }; - /*----------------------- Inline functions ----------------------------------*/ #ifndef DOXYGEN diff --git a/include/deal.II/numerics/tensor_product_matrix_creator.h b/include/deal.II/numerics/tensor_product_matrix_creator.h new file mode 100644 index 0000000000..1242e64b80 --- /dev/null +++ b/include/deal.II/numerics/tensor_product_matrix_creator.h @@ -0,0 +1,375 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2022 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + +#ifndef dealii_tensor_product_matrix_creator_h +#define dealii_tensor_product_matrix_creator_h + + +#include + +#include + +#include + +#include +#include +#include +#include + +#include +#include + +#include + +DEAL_II_NAMESPACE_OPEN + + +/** + * A namespace with functions that create input for certain standard matrices + * for the classes TensorProductMatrixSymmetricSum and + * TensorProductMatrixSymmetricSumCache. + */ +namespace TensorProductMatrixCreator +{ + /** + * Boundary type that can be used in create_laplace_tensor_product_matrix(); + */ + enum LaplaceBoundaryType + { + dirichlet, + neumann, + internal_boundary, + }; + + /** + * Create 1D mass matrix and 1D derivative matrix for a scalar + * constant-coefficient + * Laplacian for a @p dim dimensional Cartesian cell. Its boundary types + * can be specified with @p boundary_ids. The cell extent (including the cell extent + * of each neighbor) can be specified via @p cell_extent. With @p n_overlap, an + * overlap with neighboring cells can be specified. The default value is one, + * which correspond to all matrix entries restricted to the cell-local DoFs. + */ + template + std::pair, dim>, + std::array, dim>> + create_laplace_tensor_product_matrix( + const FiniteElement<1> & fe, + const Quadrature<1> & quadrature, + const dealii::ndarray &boundary_ids, + const dealii::ndarray & cell_extent, + const unsigned int n_overlap = 1); + + /** + * Same as above but the boundary IDs are extracted from the given @p cell + * and are mapped to the boundary type via the sets @p dirichlet_boundaries and @p neumann_boundaries. + */ + template + std::pair, dim>, + std::array, dim>> + create_laplace_tensor_product_matrix( + const typename Triangulation::cell_iterator &cell, + const std::set & dirichlet_boundaries, + const std::set & neumann_boundaries, + const FiniteElement<1> & fe, + const Quadrature<1> & quadrature, + const dealii::ndarray & cell_extent, + const unsigned int n_overlap = 1); + +} // namespace TensorProductMatrixCreator + + + +/*----------------------- Inline functions ----------------------------------*/ + + +namespace TensorProductMatrixCreator +{ + namespace internal + { + template + std::tuple, FullMatrix, bool> + create_reference_mass_and_stiffness_matrices( + const FiniteElement<1> &fe, + const Quadrature<1> & quadrature) + { + Triangulation<1> tria; + GridGenerator::hyper_cube(tria); + + DoFHandler<1> dof_handler(tria); + dof_handler.distribute_dofs(fe); + + MappingQ1<1> mapping; + + const unsigned int n_dofs_1D = fe.n_dofs_per_cell(); + + FullMatrix mass_matrix_reference(n_dofs_1D, n_dofs_1D); + FullMatrix derivative_matrix_reference(n_dofs_1D, n_dofs_1D); + + FEValues<1> fe_values(mapping, + fe, + quadrature, + update_values | update_gradients | + update_JxW_values); + + fe_values.reinit(tria.begin()); + + const auto lexicographic_to_hierarchic_numbering = + Utilities::invert_permutation( + FETools::hierarchic_to_lexicographic_numbering<1>( + fe.tensor_degree())); + + for (const unsigned int q_index : fe_values.quadrature_point_indices()) + for (const unsigned int i : fe_values.dof_indices()) + for (const unsigned int j : fe_values.dof_indices()) + { + mass_matrix_reference(i, j) += + (fe_values.shape_value(lexicographic_to_hierarchic_numbering[i], + q_index) * + fe_values.shape_value(lexicographic_to_hierarchic_numbering[j], + q_index) * + fe_values.JxW(q_index)); + + derivative_matrix_reference(i, j) += + (fe_values.shape_grad(lexicographic_to_hierarchic_numbering[i], + q_index) * + fe_values.shape_grad(lexicographic_to_hierarchic_numbering[j], + q_index) * + fe_values.JxW(q_index)); + } + + return {mass_matrix_reference, derivative_matrix_reference, false}; + } + } // namespace internal + + + + template + std::pair, dim>, + std::array, dim>> + create_laplace_tensor_product_matrix( + const FiniteElement<1> & fe, + const Quadrature<1> & quadrature, + const dealii::ndarray &boundary_ids, + const dealii::ndarray & cell_extent, + const unsigned int n_overlap) + { + // 1) create element mass and siffness matrix (without overlap) + const auto create_reference_mass_and_stiffness_matrices = + internal::create_reference_mass_and_stiffness_matrices( + fe, quadrature); + + const auto &M_ref = + std::get<0>(create_reference_mass_and_stiffness_matrices); + const auto &K_ref = + std::get<1>(create_reference_mass_and_stiffness_matrices); + const auto &is_dg = + std::get<2>(create_reference_mass_and_stiffness_matrices); + + AssertIndexRange(n_overlap, M_ref.n()); + AssertIndexRange(0, n_overlap); + AssertThrow(is_dg == false, ExcNotImplemented()); + + // 2) loop over all dimensions and create 1D mass and stiffness + // matrices so that boundary conditions and overlap are considered + + const unsigned int n_dofs_1D = M_ref.n(); + const unsigned int n_dofs_1D_with_overlap = M_ref.n() - 2 + 2 * n_overlap; + + std::array, dim> Ms; + std::array, dim> Ks; + + const auto clear_row_and_column = [&](const unsigned int n, auto &matrix) { + for (unsigned int i = 0; i < n_dofs_1D_with_overlap; ++i) + { + matrix[i][n] = 0.0; + matrix[n][i] = 0.0; + } + }; + + for (unsigned int d = 0; d < dim; ++d) + { + Ms[d].reinit(n_dofs_1D_with_overlap, n_dofs_1D_with_overlap); + Ks[d].reinit(n_dofs_1D_with_overlap, n_dofs_1D_with_overlap); + + // inner cell + for (unsigned int i = 0; i < n_dofs_1D; ++i) + for (unsigned int j = 0; j < n_dofs_1D; ++j) + { + const unsigned int i0 = i + n_overlap - 1; + const unsigned int j0 = j + n_overlap - 1; + Ms[d][i0][j0] = M_ref[i][j] * cell_extent[d][1]; + Ks[d][i0][j0] = K_ref[i][j] / cell_extent[d][1]; + } + + // left neighbor or left boundary + if (boundary_ids[d][0] == LaplaceBoundaryType::internal_boundary) + { + // left neighbor + Assert(cell_extent[d][0] > 0.0, ExcInternalError()); + + for (unsigned int i = 0; i < n_overlap; ++i) + for (unsigned int j = 0; j < n_overlap; ++j) + { + const unsigned int i0 = n_dofs_1D - n_overlap + i; + const unsigned int j0 = n_dofs_1D - n_overlap + j; + Ms[d][i][j] += M_ref[i0][j0] * cell_extent[d][0]; + Ks[d][i][j] += K_ref[i0][j0] / cell_extent[d][0]; + } + } + else + { + if (boundary_ids[d][0] == LaplaceBoundaryType::dirichlet) + { + // left DBC + const unsigned i0 = n_overlap - 1; + clear_row_and_column(i0, Ms[d]); + clear_row_and_column(i0, Ks[d]); + } + else if (boundary_ids[d][0] == LaplaceBoundaryType::neumann) + { + // left NBC -> nothing to do + } + else + { + AssertThrow(false, ExcNotImplemented()); + } + } + + // right neighbor or right boundary + if (boundary_ids[d][1] == LaplaceBoundaryType::internal_boundary) + { + Assert(cell_extent[d][2] > 0.0, ExcInternalError()); + + for (unsigned int i = 0; i < n_overlap; ++i) + for (unsigned int j = 0; j < n_overlap; ++j) + { + const unsigned int i0 = n_overlap + n_dofs_1D + i - 2; + const unsigned int j0 = n_overlap + n_dofs_1D + j - 2; + Ms[d][i0][j0] += M_ref[i][j] * cell_extent[d][2]; + Ks[d][i0][j0] += K_ref[i][j] / cell_extent[d][2]; + } + } + else + { + if (boundary_ids[d][1] == LaplaceBoundaryType::dirichlet) + { + // right DBC + const unsigned i0 = n_overlap + n_dofs_1D - 2; + clear_row_and_column(i0, Ms[d]); + clear_row_and_column(i0, Ks[d]); + } + else if (boundary_ids[d][1] == LaplaceBoundaryType::neumann) + { + // right NBC -> nothing to do + } + else + { + AssertThrow(false, ExcNotImplemented()); + } + } + } + + return {Ms, Ks}; + } + + + template + std::pair, dim>, + std::array, dim>> + create_laplace_tensor_product_matrix( + const typename Triangulation::cell_iterator &cell, + const std::set & dirichlet_boundaries, + const std::set & neumann_boundaries, + const FiniteElement<1> & fe, + const Quadrature<1> & quadrature, + const dealii::ndarray & cell_extent, + const unsigned int n_overlap) + { + dealii::ndarray boundary_ids; + + for (unsigned int d = 0; d < dim; ++d) + { + // left neighbor or left boundary + if ((cell->at_boundary(2 * d) == false) || + cell->has_periodic_neighbor(2 * d)) + { + // left neighbor + Assert(cell_extent[d][0] > 0.0, ExcInternalError()); + + boundary_ids[d][0] = LaplaceBoundaryType::internal_boundary; + } + else + { + const auto bid = cell->face(2 * d)->boundary_id(); + if (dirichlet_boundaries.find(bid) != + dirichlet_boundaries.end() /*DBC*/) + { + // left DBC + boundary_ids[d][0] = LaplaceBoundaryType::dirichlet; + } + else if (neumann_boundaries.find(bid) != + neumann_boundaries.end() /*NBC*/) + { + // left NBC + boundary_ids[d][0] = LaplaceBoundaryType::neumann; + } + else + { + AssertThrow(false, ExcNotImplemented()); + } + } + + // right neighbor or right boundary + if ((cell->at_boundary(2 * d + 1) == false) || + cell->has_periodic_neighbor(2 * d + 1)) + { + Assert(cell_extent[d][2] > 0.0, ExcInternalError()); + + boundary_ids[d][1] = LaplaceBoundaryType::internal_boundary; + } + else + { + const auto bid = cell->face(2 * d + 1)->boundary_id(); + if (dirichlet_boundaries.find(bid) != + dirichlet_boundaries.end() /*DBC*/) + { + // right DBC + boundary_ids[d][1] = LaplaceBoundaryType::dirichlet; + } + else if (neumann_boundaries.find(bid) != + neumann_boundaries.end() /*NBC*/) + { + // right NBC + boundary_ids[d][1] = LaplaceBoundaryType::neumann; + } + else + { + AssertThrow(false, ExcNotImplemented()); + } + } + } + + return create_laplace_tensor_product_matrix( + fe, quadrature, boundary_ids, cell_extent, n_overlap); + } + +} // namespace TensorProductMatrixCreator + + + +DEAL_II_NAMESPACE_CLOSE + +#endif -- 2.39.5