From 2a88f9fc67daef48a266a1e7a22d025fd1173257 Mon Sep 17 00:00:00 2001 From: =?utf8?q?Micha=C5=82=20Wichrowski?= Date: Thu, 17 Apr 2025 20:59:53 +0200 Subject: [PATCH] move functionality to the correct place --- include/deal.II/numerics/matrix_creator.h | 98 ------- .../numerics/tensor_product_matrix_creator.h | 246 ++++++++++++++++++ source/numerics/matrix_creator.cc | 155 ----------- 3 files changed, 246 insertions(+), 253 deletions(-) diff --git a/include/deal.II/numerics/matrix_creator.h b/include/deal.II/numerics/matrix_creator.h index 864d2818bd..9af6e15155 100644 --- a/include/deal.II/numerics/matrix_creator.h +++ b/include/deal.II/numerics/matrix_creator.h @@ -602,104 +602,6 @@ namespace MatrixCreator - /** - * Computes a 1D cell mass matrix for a given finite element.This function is - * intended to provide a mass matrix for tensor product operators. - * - * @param fe The finite element object defining the shape functions. - * - * @param h The cell size used to scale the integration (typically the - * element size). - * - * @param include_endpoints A pair of booleans indicating whether to include - * the first and last dof in the matrix. By excluding those DoFs, one can set - * a Dirichlet BC on either side. Setting either one to false only makes sense - * if the indexing of DoFs in the matrix is lexicographic, i.e., for FE_Q - * numbering is required, while FE_DGQ can be processed as it is. - * - * @param numbering A vector of unsigned integers representing the - * numbering of the degrees of freedom. If empty, the - * numbering of the finite element is used. - * - * @return A FullMatrix representing the assembled mass matrix. - * - */ - FullMatrix - create_1d_cell_mass_matrix( - const FiniteElement<1> &fe, - const double &h, - const std::pair include_endpoints = {true, true}, - std::vector numbering = std::vector()); - - - - /** - * Computes a 1D cell derivative matrix for a given finite element. - * - * @param fe The finite element object defining the shape functions. - * - * @param h The cell size used to scale the integration (typically the - * element size). - * - * @param include_endpoints A pair of booleans indicating whether to include - * the first and last dof in the matrix. By excluding those DoFs, one can set - * a Dirichlet BC on either side. Setting either one to false only makes sense - * if the indexing of DoFs in the matrix is lexicographic, i.e., for FE_Q - * numbering is required, while FE_DGQ can be processed as it is. - * - * @param numbering A vector of unsigned integers representing the - * numbering of the degrees of freedom. If empty, the - * numbering of the finite element is used. - */ - FullMatrix - create_1d_cell_derivative_matrix( - const FiniteElement<1> &fe, - const double &h, - const std::pair include_endpoints = {true, true}, - std::vector numbering = std::vector()); - - - /** - * This function assembles a global matrix from a given cell matrix, assuming - * a 1D discretization with a specified number of cells and overlap between - * them. - * - * @param cell_matrix The local cell matrix to be assembled into the - * global matrix. It is assumed that the cell matrix has a size corresponding - * to the number of degrees of freedom per cell. The cell matrix should be - * ordered lexicographically, meaning that FE_Q should be renumbered - * - * @param n_cells The number of cells in the 1D discretization. - * - * @param overlap The number of degrees of freedom that overlap between - * adjacent cells. For discretization with FE_Q elements, the overlap should - * be set to 1, while for FE_DGQ elements, it should be set to 0. - * - * @param include_endpoints A pair of booleans indicating whether to - * include the left and right most dofs of the domain in the constructed - * matrix. The default value is {true, true}, which includes both endpoints. - * - * @return The assembled global matrix. The size of the matrix is determined by - * the number of cells, the overlap, and whether the endpoints are included. - * - * @note The size of the cell matrix must be consistent with the overlap - * parameter. Specifically, if the cell matrix is $n \times n$, then the - * overlap should be less than $n$. - * - * @warning This function returns a full matrix as it is intended to - * be used inside local tensor product operators. - * If the number of cells is above 10 an exception will be thrown. - * - */ - FullMatrix - create_1D_discretization_matrix( - FullMatrix &cell_matrix, - const unsigned int &n_cells, - const unsigned int &overlap, - const std::pair include_endpoints = {true, true}); - - - /** * Exception */ diff --git a/include/deal.II/numerics/tensor_product_matrix_creator.h b/include/deal.II/numerics/tensor_product_matrix_creator.h index e871054719..a4cc4b9353 100644 --- a/include/deal.II/numerics/tensor_product_matrix_creator.h +++ b/include/deal.II/numerics/tensor_product_matrix_creator.h @@ -87,6 +87,103 @@ namespace TensorProductMatrixCreator const dealii::ndarray &cell_extent, const unsigned int n_overlap = 1); + + + /** + * Computes a 1D cell mass matrix for a given finite element.This function is + * intended to provide a mass matrix for tensor product operators. + * + * @param fe The finite element object defining the shape functions. + * + * @param h The cell size used to scale the integration (typically the + * element size). + * + * @param include_endpoints A pair of booleans indicating whether to include + * the first and last dof in the matrix. By excluding those DoFs, one can set + * a Dirichlet BC on either side. Setting either one to false only makes sense + * if the indexing of DoFs in the matrix is lexicographic, i.e., for FE_Q + * numbering is required, while FE_DGQ can be processed as it is. + * + * @param numbering A vector of unsigned integers representing the + * numbering of the degrees of freedom. If empty, the + * numbering of the finite element is used. + * + * @return A FullMatrix representing the assembled mass matrix. + * + */ + FullMatrix + create_1d_cell_mass_matrix( + const FiniteElement<1> &fe, + const double &h, + const std::pair include_endpoints = {true, true}, + std::vector numbering = std::vector()); + + + + /** + * Computes a 1D cell derivative matrix for a given finite element. + * + * @param fe The finite element object defining the shape functions. + * + * @param h The cell size used to scale the integration (typically the + * element size). + * + * @param include_endpoints A pair of booleans indicating whether to include + * the first and last dof in the matrix. By excluding those DoFs, one can set + * a Dirichlet BC on either side. Setting either one to false only makes sense + * if the indexing of DoFs in the matrix is lexicographic, i.e., for FE_Q + * numbering is required, while FE_DGQ can be processed as it is. + * + * @param numbering A vector of unsigned integers representing the + * numbering of the degrees of freedom. If empty, the + * numbering of the finite element is used. + */ + FullMatrix + create_1d_cell_derivative_matrix( + const FiniteElement<1> &fe, + const double &h, + const std::pair include_endpoints = {true, true}, + std::vector numbering = std::vector()); + + + /** + * This function assembles a global matrix from a given cell matrix, assuming + * a 1D discretization with a specified number of cells and overlap between + * them. + * + * @param cell_matrix The local cell matrix to be assembled into the + * global matrix. It is assumed that the cell matrix has a size corresponding + * to the number of degrees of freedom per cell. The cell matrix should be + * ordered lexicographically, meaning that FE_Q should be renumbered + * + * @param n_cells The number of cells in the 1D discretization. + * + * @param overlap The number of degrees of freedom that overlap between + * adjacent cells. For discretization with FE_Q elements, the overlap should + * be set to 1, while for FE_DGQ elements, it should be set to 0. + * + * @param include_endpoints A pair of booleans indicating whether to + * include the left and right most dofs of the domain in the constructed + * matrix. The default value is {true, true}, which includes both endpoints. + * + * @return The assembled global matrix. The size of the matrix is determined by + * the number of cells, the overlap, and whether the endpoints are included. + * + * @note The size of the cell matrix must be consistent with the overlap + * parameter. Specifically, if the cell matrix is $n \times n$, then the + * overlap should be less than $n$. + * + * @warning This function returns a full matrix, + * if the number of cells is above 10 an exception will be thrown. + * + */ + FullMatrix + create_1D_discretization_matrix( + FullMatrix &cell_matrix, + const unsigned int &n_cells, + const unsigned int &overlap, + const std::pair include_endpoints = {true, true}); + } // namespace TensorProductMatrixCreator @@ -379,6 +476,155 @@ namespace TensorProductMatrixCreator fe, quadrature, boundary_ids, cell_extent, n_overlap); } + + FullMatrix + create_1d_cell_mass_matrix(const FiniteElement<1> &fe, + const double &h, + const std::pair include_endpoints, + std::vector numbering) + { + if (dynamic_cast *>(&fe) == nullptr && + numbering.size() == 0) + { + Assert( + include_endpoints.first == true && include_endpoints.second == true, + ExcMessage( + "You tried to generate a 1D mass matrix with excluding boundary " + "dofs for a non-DGQ element without providing a numbering.")); + } + + + if (numbering.size() == 0) + { + numbering.resize(fe.dofs_per_cell); + std::iota(numbering.begin(), numbering.end(), 0); + } + const unsigned int degree = fe.degree; + const unsigned int n_dofs_per_cell = fe.dofs_per_cell; + const double &JxW = h; + QGauss<1> quadrature(degree + 1); + + FullMatrix cell_mass_matrix(n_dofs_per_cell, n_dofs_per_cell); + cell_mass_matrix = 0; + + unsigned int start_dof = include_endpoints.first ? 0 : 1; + unsigned int end_dof = + include_endpoints.second ? n_dofs_per_cell : n_dofs_per_cell - 1; + const unsigned int shift = include_endpoints.first ? 0 : 1; + + for (unsigned int i = start_dof; i < end_dof; ++i) + for (unsigned int j = start_dof; j < end_dof; ++j) + for (unsigned int q = 0; q < quadrature.size(); ++q) + cell_mass_matrix(i - shift, j - shift) += + (fe.shape_value(numbering[i], quadrature.point(q)) * + fe.shape_value(numbering[j], quadrature.point(q))) * + JxW * quadrature.weight(q); + + return cell_mass_matrix; + } + + + FullMatrix + create_1d_cell_derivative_matrix( + const FiniteElement<1> &fe, + const double &h, + const std::pair include_endpoints, + std::vector numbering) + { + if (dynamic_cast *>(&fe) == nullptr && + numbering.size() == 0) + { + Assert( + include_endpoints.first == true && include_endpoints.second == true, + ExcMessage( + "You tried to generate a 1D derivative matrix with excluding boundary " + "dofs for a non-DGQ element without providing a numbering.")); + } + + if (numbering.size() == 0) + { + numbering.resize(fe.dofs_per_cell); + std::iota(numbering.begin(), numbering.end(), 0); + } + + const unsigned int degree = fe.degree; + const unsigned int n_dofs_per_cell = fe.dofs_per_cell; + const double &JxW = h; + QGauss<1> quadrature(degree + 1); + + FullMatrix cell_matrix(n_dofs_per_cell, n_dofs_per_cell); + cell_matrix = 0; + + unsigned int start_dof = include_endpoints.first ? 0 : 1; + unsigned int end_dof = + include_endpoints.second ? n_dofs_per_cell : n_dofs_per_cell - 1; + const unsigned int shift = include_endpoints.first ? 0 : 1; + + for (unsigned int i = start_dof; i < end_dof; ++i) + for (unsigned int j = start_dof; j < end_dof; ++j) + for (unsigned int q = 0; q < quadrature.size(); ++q) + cell_matrix(i - shift, j - shift) += + (fe.shape_grad(numbering[i], quadrature.point(q)) / h * + fe.shape_grad(numbering[j], quadrature.point(q))) / + h * JxW * quadrature.weight(q); + + return cell_matrix; + } + + FullMatrix + create_1D_discretization_matrix(FullMatrix &cell_matrix, + const unsigned int &n_cells, + const unsigned int &overlap, + const std::pair include_endpoints) + { + const unsigned int n_dofs_per_cell = cell_matrix.n(); + + Assert(cell_matrix.m() == n_dofs_per_cell, + ExcMessage( + "The provided cell mass matrix must be a square matrix.")); + AssertThrow( + n_cells <= 10, + ExcMessage( + "create_1D_discretization_matrix() returns a full matrix and is not meant to be used with a larger number of cells. ")); + Assert(n_cells > 0, + ExcMessage("You are trying to get a mass matrix of zero cells.")); + Assert(overlap < n_dofs_per_cell, + ExcMessage("The overlap must be smaller than the number of dofs.")); + + unsigned int n_total_dofs = + n_cells * n_dofs_per_cell - overlap * (n_cells - 1); + + if (!include_endpoints.first) + n_total_dofs -= 1; + if (!include_endpoints.second) + n_total_dofs -= 1; + + FullMatrix result_matrix(n_total_dofs, n_total_dofs); + result_matrix = 0; + + for (unsigned int cell = 0; cell < n_cells; ++cell) + { + const unsigned int dof_shift = + cell * overlap + !include_endpoints.first; + + const unsigned int start_dof = + (cell == 0 && !include_endpoints.first) ? 1 : 0; + + const unsigned int end_dof = + (cell == n_cells - 1 && !include_endpoints.second) ? + n_dofs_per_cell - 1 : + n_dofs_per_cell; + for (unsigned int i = start_dof; i < end_dof; ++i) + for (unsigned int j = start_dof; j < end_dof; ++j) + { + result_matrix(i + cell * n_dofs_per_cell - dof_shift, + j + cell * n_dofs_per_cell - dof_shift) += + cell_matrix(i, j); + } + } + return result_matrix; + } + } // namespace TensorProductMatrixCreator diff --git a/source/numerics/matrix_creator.cc b/source/numerics/matrix_creator.cc index fbb9972864..6e5187f529 100644 --- a/source/numerics/matrix_creator.cc +++ b/source/numerics/matrix_creator.cc @@ -24,159 +24,4 @@ DEAL_II_NAMESPACE_OPEN #endif #include "numerics/matrix_creator.inst" -#if defined(SPLIT_INSTANTIATIONS_INDEX) && SPLIT_INSTANTIATIONS_INDEX == 0 -namespace MatrixCreator -{ - - FullMatrix - create_1d_cell_mass_matrix(const FiniteElement<1> &fe, - const double &h, - const std::pair include_endpoints, - std::vector numbering) - { - if (dynamic_cast *>(&fe) == nullptr && - numbering.size() == 0) - { - Assert( - include_endpoints.first == true && include_endpoints.second == true, - ExcMessage( - "You tried to generate a 1D mass matrix with excluding boundary " - "dofs for a non-DGQ element without providing a numbering.")); - } - - - if (numbering.size() == 0) - { - numbering.resize(fe.dofs_per_cell); - std::iota(numbering.begin(), numbering.end(), 0); - } - const unsigned int degree = fe.degree; - const unsigned int n_dofs_per_cell = fe.dofs_per_cell; - const double &JxW = h; - QGauss<1> quadrature(degree + 1); - - FullMatrix cell_mass_matrix(n_dofs_per_cell, n_dofs_per_cell); - cell_mass_matrix = 0; - - unsigned int start_dof = include_endpoints.first ? 0 : 1; - unsigned int end_dof = - include_endpoints.second ? n_dofs_per_cell : n_dofs_per_cell - 1; - const unsigned int shift = include_endpoints.first ? 0 : 1; - - for (unsigned int i = start_dof; i < end_dof; ++i) - for (unsigned int j = start_dof; j < end_dof; ++j) - for (unsigned int q = 0; q < quadrature.size(); ++q) - cell_mass_matrix(i - shift, j - shift) += - (fe.shape_value(numbering[i], quadrature.point(q)) * - fe.shape_value(numbering[j], quadrature.point(q))) * - JxW * quadrature.weight(q); - - return cell_mass_matrix; - } - - - FullMatrix - create_1d_cell_derivative_matrix( - const FiniteElement<1> &fe, - const double &h, - const std::pair include_endpoints, - std::vector numbering) - { - if (dynamic_cast *>(&fe) == nullptr && - numbering.size() == 0) - { - Assert( - include_endpoints.first == true && include_endpoints.second == true, - ExcMessage( - "You tried to generate a 1D derivative matrix with excluding boundary " - "dofs for a non-DGQ element without providing a numbering.")); - } - - if (numbering.size() == 0) - { - numbering.resize(fe.dofs_per_cell); - std::iota(numbering.begin(), numbering.end(), 0); - } - - const unsigned int degree = fe.degree; - const unsigned int n_dofs_per_cell = fe.dofs_per_cell; - const double &JxW = h; - QGauss<1> quadrature(degree + 1); - - FullMatrix cell_matrix(n_dofs_per_cell, n_dofs_per_cell); - cell_matrix = 0; - - unsigned int start_dof = include_endpoints.first ? 0 : 1; - unsigned int end_dof = - include_endpoints.second ? n_dofs_per_cell : n_dofs_per_cell - 1; - const unsigned int shift = include_endpoints.first ? 0 : 1; - - for (unsigned int i = start_dof; i < end_dof; ++i) - for (unsigned int j = start_dof; j < end_dof; ++j) - for (unsigned int q = 0; q < quadrature.size(); ++q) - cell_matrix(i - shift, j - shift) += - (fe.shape_grad(numbering[i], quadrature.point(q)) / h * - fe.shape_grad(numbering[j], quadrature.point(q))) / - h * JxW * quadrature.weight(q); - - return cell_matrix; - } - - FullMatrix - create_1D_discretization_matrix(FullMatrix &cell_matrix, - const unsigned int &n_cells, - const unsigned int &overlap, - const std::pair include_endpoints) - { - const unsigned int n_dofs_per_cell = cell_matrix.n(); - - Assert(cell_matrix.m() == n_dofs_per_cell, - ExcMessage( - "The provided cell mass matrix must be a square matrix.")); - AssertThrow( - n_cells <= 10, - ExcMessage( - "create_1D_discretization_matrix() returns a full matrix and is not meant to be used with a larger number of cells. ")); - Assert(n_cells > 0, - ExcMessage("You are trying to get a mass matrix of zero cells.")); - Assert(overlap < n_dofs_per_cell, - ExcMessage("The overlap must be smaller than the number of dofs.")); - - unsigned int n_total_dofs = - n_cells * n_dofs_per_cell - overlap * (n_cells - 1); - - if (!include_endpoints.first) - n_total_dofs -= 1; - if (!include_endpoints.second) - n_total_dofs -= 1; - - FullMatrix result_matrix(n_total_dofs, n_total_dofs); - result_matrix = 0; - - for (unsigned int cell = 0; cell < n_cells; ++cell) - { - const unsigned int dof_shift = - cell * overlap + !include_endpoints.first; - - const unsigned int start_dof = - (cell == 0 && !include_endpoints.first) ? 1 : 0; - - const unsigned int end_dof = - (cell == n_cells - 1 && !include_endpoints.second) ? - n_dofs_per_cell - 1 : - n_dofs_per_cell; - for (unsigned int i = start_dof; i < end_dof; ++i) - for (unsigned int j = start_dof; j < end_dof; ++j) - { - result_matrix(i + cell * n_dofs_per_cell - dof_shift, - j + cell * n_dofs_per_cell - dof_shift) += - cell_matrix(i, j); - } - } - return result_matrix; - } -} // namespace MatrixCreator - -#endif // SPLIT_INSTANTIATIONS_INDEX == 0 - DEAL_II_NAMESPACE_CLOSE -- 2.39.5