From e2bb9cec7bda430d63d56223509700a77c803a8f Mon Sep 17 00:00:00 2001 From: =?utf8?q?Micha=C5=82=20Wichrowski?= Date: Thu, 1 May 2025 18:45:48 +0200 Subject: [PATCH] updated documentation moved implementation to .cc file --- .../numerics/tensor_product_matrix_creator.h | 59 ++-------------- .../tensor_product_matrix_creator_test.cc | 70 +++++++++++++++++++ 2 files changed, 77 insertions(+), 52 deletions(-) create mode 100644 scratch/mwichro/dealii/src/numerics/tensor_product_matrix_creator_test.cc diff --git a/include/deal.II/numerics/tensor_product_matrix_creator.h b/include/deal.II/numerics/tensor_product_matrix_creator.h index be3f9bc151..5ea67a6ab6 100644 --- a/include/deal.II/numerics/tensor_product_matrix_creator.h +++ b/include/deal.II/numerics/tensor_product_matrix_creator.h @@ -187,6 +187,7 @@ namespace TensorProductMatrixCreator /** * @brief Create a 1D ghost penalty matrix for a given finite element. + * Ghost penalty is used for stabilization in cutFEM, see step-85. * Implemented only for FE_Q. * * @param fe The finite element space used for discretization. @@ -211,8 +212,9 @@ namespace TensorProductMatrixCreator /** * @brief Create a 1D ghost penalty matrix. * - * This function creates a ghost penalty matrix in 1D. The matrix is created - * using the derivatives of the polynomial basis functions. + * This function creates a single component of the sum over all derivatives in + * the 1D ghost penalty matrix. The matrix is created using the + * derivatives of the polynomial basis functions. * * @param polynomial_basis_derivative A vector of polynomial basis functions. * These should be the derivatives of the basis functions used to represent @@ -522,54 +524,6 @@ 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, @@ -648,10 +602,11 @@ namespace TensorProductMatrixCreator FullMatrix result_matrix(n_total_dofs, n_total_dofs); result_matrix = 0; + const unsigned int left_shift = include_endpoints.first ? 0 : 1; + for (unsigned int cell = 0; cell < n_cells; ++cell) { - const unsigned int dof_shift = - cell * overlap + !include_endpoints.first; + const unsigned int dof_shift = cell * overlap + left_shift; const unsigned int start_dof = (cell == 0 && !include_endpoints.first) ? 1 : 0; diff --git a/scratch/mwichro/dealii/src/numerics/tensor_product_matrix_creator_test.cc b/scratch/mwichro/dealii/src/numerics/tensor_product_matrix_creator_test.cc new file mode 100644 index 0000000000..fd1a39ad23 --- /dev/null +++ b/scratch/mwichro/dealii/src/numerics/tensor_product_matrix_creator_test.cc @@ -0,0 +1,70 @@ +#include +#include + +#include + +#include + +#include + +#include + + +DEAL_II_NAMESPACE_OPEN + + +namespace TensorProductMatrixCreator +{ + + + 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; + } + + +} // namespace TensorProductMatrixCreator + + +DEAL_II_NAMESPACE_CLOSE \ No newline at end of file -- 2.39.5