From: MichaƂ Wichrowski Date: Thu, 1 May 2025 18:25:58 +0000 (+0200) Subject: tests for 1D matrices X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=3d5dd1ff3cd9933cd744e3d212390c064e867002;p=dealii.git tests for 1D matrices --- diff --git a/include/deal.II/numerics/tensor_product_matrix_creator.h b/include/deal.II/numerics/tensor_product_matrix_creator.h index c1caa78c6a..60b672eb59 100644 --- a/include/deal.II/numerics/tensor_product_matrix_creator.h +++ b/include/deal.II/numerics/tensor_product_matrix_creator.h @@ -113,10 +113,11 @@ namespace TensorProductMatrixCreator * @return A FullMatrix representing the assembled mass matrix. * */ - FullMatrix + template + FullMatrix create_1d_cell_mass_matrix( const FiniteElement<1> &fe, - const double &h, + const Number &h, const std::pair include_endpoints = {true, true}, std::vector numbering = std::vector()); @@ -140,10 +141,11 @@ namespace TensorProductMatrixCreator * numbering of the degrees of freedom. If empty, the * numbering of the finite element is used. */ - FullMatrix + template + FullMatrix create_1d_cell_laplace_matrix( const FiniteElement<1> &fe, - const double &h, + const Number &h, const std::pair include_endpoints = {true, true}, std::vector numbering = std::vector()); @@ -179,9 +181,10 @@ namespace TensorProductMatrixCreator * if the number of cells is above 10 an exception will be thrown. * */ - FullMatrix + template + FullMatrix create_1D_discretization_matrix( - FullMatrix &cell_matrix, + FullMatrix &cell_matrix, const unsigned int &n_cells, const unsigned int &overlap, const std::pair include_endpoints = {true, true}); @@ -203,11 +206,12 @@ namespace TensorProductMatrixCreator * * @return A full matrix representing the 1D ghost penalty matrix. */ - FullMatrix + template + FullMatrix create_1d_ghost_penalty_matrix( const FiniteElement<1> &fe, - const double h, - std::vector coefficients = std::vector()); + const Number h, + std::vector coefficients = std::vector()); @@ -228,7 +232,8 @@ namespace TensorProductMatrixCreator * * @return A full matrix representing the 1D ghost penalty. */ - FullMatrix + template + FullMatrix create_1d_ghost_penalty_matrix( const std::vector> &polynomial_basis_derivative, @@ -526,9 +531,57 @@ namespace TensorProductMatrixCreator fe, quadrature, boundary_ids, cell_extent, n_overlap); } - FullMatrix + template + FullMatrix + create_1d_cell_mass_matrix(const FiniteElement<1> &fe, + const Number &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 Number &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_value(numbering[i], quadrature.point(q)) * + fe.shape_value(numbering[j], quadrature.point(q))) * + JxW * quadrature.weight(q); + + return cell_matrix; + } + + template + FullMatrix create_1d_cell_laplace_matrix(const FiniteElement<1> &fe, - const double &h, + const Number &h, const std::pair include_endpoints, std::vector numbering) { @@ -550,10 +603,10 @@ namespace TensorProductMatrixCreator const unsigned int degree = fe.degree; const unsigned int n_dofs_per_cell = fe.dofs_per_cell; - const double &JxW = h; + const Number &JxW = h; QGauss<1> quadrature(degree + 1); - FullMatrix cell_matrix(n_dofs_per_cell, n_dofs_per_cell); + FullMatrix cell_matrix(n_dofs_per_cell, n_dofs_per_cell); cell_matrix = 0; unsigned int start_dof = include_endpoints.first ? 0 : 1; @@ -572,8 +625,9 @@ namespace TensorProductMatrixCreator return cell_matrix; } - FullMatrix - create_1D_discretization_matrix(FullMatrix &cell_matrix, + template + FullMatrix + create_1D_discretization_matrix(FullMatrix &cell_matrix, const unsigned int &n_cells, const unsigned int &overlap, const std::pair include_endpoints) @@ -600,7 +654,7 @@ namespace TensorProductMatrixCreator if (!include_endpoints.second) n_total_dofs -= 1; - FullMatrix result_matrix(n_total_dofs, n_total_dofs); + FullMatrix result_matrix(n_total_dofs, n_total_dofs); result_matrix = 0; const unsigned int left_shift = include_endpoints.first ? 0 : 1; @@ -629,10 +683,11 @@ namespace TensorProductMatrixCreator - FullMatrix + template + FullMatrix create_1d_ghost_penalty_matrix(const FiniteElement<1> &fe, - const double h, - std::vector coefficients) + const Number h, + std::vector coefficients) { Assert(dynamic_cast *>(&fe) != nullptr, ExcNotImplemented()); Assert(h > 0, ExcMessage("Provided element size h is negative")); @@ -683,13 +738,13 @@ namespace TensorProductMatrixCreator } - FullMatrix penalty_matrix = + FullMatrix penalty_matrix = create_1d_ghost_penalty_matrix(polynomial_basis[1]); penalty_matrix *= coefficients[0]; for (unsigned int k = 2; k < degree + 1; ++k) { - FullMatrix kth_matrix = + FullMatrix kth_matrix = create_1d_ghost_penalty_matrix(polynomial_basis[k]); penalty_matrix.add(coefficients[k - 1], kth_matrix); } @@ -699,8 +754,8 @@ namespace TensorProductMatrixCreator } - - FullMatrix + template + FullMatrix create_1d_ghost_penalty_matrix( const std::vector> &polynomial_basis_derivative, @@ -710,7 +765,7 @@ namespace TensorProductMatrixCreator const unsigned int n_total_dofs = 2 * n_dofs_per_cell - overlap; const unsigned int shift = n_dofs_per_cell - overlap; - FullMatrix penalty_matrix(n_total_dofs, n_total_dofs); + FullMatrix penalty_matrix(n_total_dofs, n_total_dofs); std::vector values_left(n_dofs_per_cell); std::vector values_right(n_dofs_per_cell); diff --git a/tests/numerics/tensor_product_creator_01.cc b/tests/numerics/tensor_product_creator_01.cc index 1a9ca29420..bf30bc1227 100644 --- a/tests/numerics/tensor_product_creator_01.cc +++ b/tests/numerics/tensor_product_creator_01.cc @@ -23,32 +23,27 @@ void -print_matrix(FullMatrix &matrix) +print_matrix(const FullMatrix &matrix, const std::string &label) { - for (unsigned int i = 0; i < matrix.m(); ++i) - { - for (unsigned int j = 0; j < matrix.n(); ++j) - deallog << std::setw(10) << std::fixed << std::setprecision(4) - << matrix(i, j) << " "; - - deallog << std::endl; - } - deallog << std::endl; + std::stringstream sstring; + sstring << label << "\n"; + matrix.print_formatted(sstring, 4, true, 10, "0."); + deallog << sstring.str() << std::endl; } template void test() { - FE_Q<1> fe_q(fe_degree); - FE_Q<1> fe_dgq(fe_degree); - auto mass_matrix = + FE_Q<1> fe_q(fe_degree); + FE_DGQ<1> fe_dgq(fe_degree); + auto mass_matrix = TensorProductMatrixCreator::create_1d_cell_mass_matrix(fe_q, - 1, + 1., {true, true}); - mass_matrix.print_formatted(std::cout, 4, true, 10, "0."); + print_matrix(mass_matrix, + "Mass matrix for FE_Q<1>(" + std::to_string(fe_degree) + "):"); - std::cout << std::endl; FEEvaluation<1, fe_degree, fe_degree + 1, 1, double> fe_eval( @@ -58,43 +53,47 @@ test() auto mass_matrix_renumbered = TensorProductMatrixCreator::create_1d_cell_mass_matrix( fe_q, 1., {true, true}, fe_eval.get_internal_dof_numbering()); - mass_matrix_renumbered.print_formatted(std::cout, 3, true, 10, "0."); - std::cout << std::endl; + print_matrix(mass_matrix_renumbered, + "Mass matrix for FE_Q<1>(" + std::to_string(fe_degree) + + ") renumbered with FEEvaluation:"); auto mass_matrix_dgq = TensorProductMatrixCreator::create_1d_cell_mass_matrix(fe_dgq, 1., {true, true}); + print_matrix(mass_matrix_dgq, + "Mass matrix for FE_DGQ<1>(" + std::to_string(fe_degree) + "):"); + for (unsigned int i = 0; i < mass_matrix_dgq.m(); ++i) for (unsigned int j = 0; j < mass_matrix_dgq.n(); ++j) if (std::abs(mass_matrix_dgq(i, j) - mass_matrix_renumbered(i, j)) > 1e-12) - std::cout << "Different at " << i << " " << j << std::endl; + deallog << "Different at " << i << " " << j << std::endl; - std::cout << std::endl; auto mass_matrix_patch = TensorProductMatrixCreator::create_1D_discretization_matrix( - mass_matrix_renumbered, 4, 1, {true, true}); + mass_matrix_renumbered, 2, 1, {true, true}); - mass_matrix_patch.print_formatted(std::cout, 4, true, 10, "0."); - std::cout << std::endl; + print_matrix(mass_matrix_patch, "Mass matrix patch:"); auto laplace_matrix = TensorProductMatrixCreator::create_1d_cell_laplace_matrix(fe_q, 1., {true, true}); - laplace_matrix.print_formatted(std::cout, 4, true, 10, "0."); - std::cout << std::endl; + print_matrix(laplace_matrix, + "Laplace matrix for FE_Q<1>(" + std::to_string(fe_degree) + + "):"); auto laplace_matrix_renumbered = TensorProductMatrixCreator::create_1d_cell_laplace_matrix( fe_q, 1., {true, true}, fe_eval.get_internal_dof_numbering()); - laplace_matrix_renumbered.print_formatted(std::cout, 4, true, 10, "0."); - std::cout << std::endl; + print_matrix(laplace_matrix_renumbered, + "Laplace matrix for FE_Q<1>(" + std::to_string(fe_degree) + + ") renumbered with FEEvaluation:"); } @@ -106,5 +105,6 @@ main() deallog << std::fixed; test<1>(); + test<2>(); test<3>(); } \ No newline at end of file diff --git a/tests/numerics/tensor_product_creator_01.output b/tests/numerics/tensor_product_creator_01.output index f7c81d7707..fa3194c240 100644 --- a/tests/numerics/tensor_product_creator_01.output +++ b/tests/numerics/tensor_product_creator_01.output @@ -1,13 +1,97 @@ -DEAL::1.0000 -2.0000 1.0000 -DEAL::-2.0000 4.0000 -2.0000 -DEAL::1.0000 -2.0000 1.0000 -DEAL:: -DEAL::51.0000 -130.3444 160.3444 -62.0000 -32.3607 12.3607 1.0000 -DEAL::-130.3444 338.1966 -425.0000 148.8854 125.0000 -69.0983 12.3607 -DEAL::160.3444 -425.0000 561.8034 -208.8854 -180.9017 125.0000 -32.3607 -DEAL::-62.0000 148.8854 -208.8854 244.0000 -208.8854 148.8854 -62.0000 -DEAL::-32.3607 125.0000 -180.9017 -208.8854 561.8034 -425.0000 160.3444 -DEAL::12.3607 -69.0983 125.0000 148.8854 -425.0000 338.1966 -130.3444 -DEAL::1.0000 12.3607 -32.3607 -62.0000 160.3444 -130.3444 51.0000 -DEAL:: +DEAL::Mass matrix for FE_Q<1>(1): +3.3333e-01 1.6667e-01 +1.6667e-01 3.3333e-01 + +DEAL::Mass matrix for FE_Q<1>(1) renumbered with FEEvaluation: +3.3333e-01 1.6667e-01 +1.6667e-01 3.3333e-01 + +DEAL::Mass matrix for FE_DGQ<1>(1): +3.3333e-01 1.6667e-01 +1.6667e-01 3.3333e-01 + +DEAL::Mass matrix patch: +3.3333e-01 1.6667e-01 0. +1.6667e-01 6.6667e-01 1.6667e-01 + 0. 1.6667e-01 3.3333e-01 + +DEAL::Laplace matrix for FE_Q<1>(1): +1.0000e+00 -1.0000e+00 +-1.0000e+00 1.0000e+00 + +DEAL::Laplace matrix for FE_Q<1>(1) renumbered with FEEvaluation: +1.0000e+00 -1.0000e+00 +-1.0000e+00 1.0000e+00 + +DEAL::Mass matrix for FE_Q<1>(2): +1.3333e-01 -3.3333e-02 6.6667e-02 +-3.3333e-02 1.3333e-01 6.6667e-02 +6.6667e-02 6.6667e-02 5.3333e-01 + +DEAL::Mass matrix for FE_Q<1>(2) renumbered with FEEvaluation: +1.3333e-01 6.6667e-02 -3.3333e-02 +6.6667e-02 5.3333e-01 6.6667e-02 +-3.3333e-02 6.6667e-02 1.3333e-01 + +DEAL::Mass matrix for FE_DGQ<1>(2): +1.3333e-01 6.6667e-02 -3.3333e-02 +6.6667e-02 5.3333e-01 6.6667e-02 +-3.3333e-02 6.6667e-02 1.3333e-01 + +DEAL::Mass matrix patch: +1.3333e-01 6.6667e-02 -3.3333e-02 0. 0. +6.6667e-02 5.3333e-01 6.6667e-02 0. 0. +-3.3333e-02 6.6667e-02 2.6667e-01 6.6667e-02 -3.3333e-02 + 0. 0. 6.6667e-02 5.3333e-01 6.6667e-02 + 0. 0. -3.3333e-02 6.6667e-02 1.3333e-01 + +DEAL::Laplace matrix for FE_Q<1>(2): +2.3333e+00 3.3333e-01 -2.6667e+00 +3.3333e-01 2.3333e+00 -2.6667e+00 +-2.6667e+00 -2.6667e+00 5.3333e+00 + +DEAL::Laplace matrix for FE_Q<1>(2) renumbered with FEEvaluation: +2.3333e+00 -2.6667e+00 3.3333e-01 +-2.6667e+00 5.3333e+00 -2.6667e+00 +3.3333e-01 -2.6667e+00 2.3333e+00 + +DEAL::Mass matrix for FE_Q<1>(3): +7.1429e-02 1.1905e-02 2.6620e-02 -2.6620e-02 +1.1905e-02 7.1429e-02 -2.6620e-02 2.6620e-02 +2.6620e-02 -2.6620e-02 3.5714e-01 5.9524e-02 +-2.6620e-02 2.6620e-02 5.9524e-02 3.5714e-01 + +DEAL::Mass matrix for FE_Q<1>(3) renumbered with FEEvaluation: +7.1429e-02 2.6620e-02 -2.6620e-02 1.1905e-02 +2.6620e-02 3.5714e-01 5.9524e-02 -2.6620e-02 +-2.6620e-02 5.9524e-02 3.5714e-01 2.6620e-02 +1.1905e-02 -2.6620e-02 2.6620e-02 7.1429e-02 + +DEAL::Mass matrix for FE_DGQ<1>(3): +7.1429e-02 2.6620e-02 -2.6620e-02 1.1905e-02 +2.6620e-02 3.5714e-01 5.9524e-02 -2.6620e-02 +-2.6620e-02 5.9524e-02 3.5714e-01 2.6620e-02 +1.1905e-02 -2.6620e-02 2.6620e-02 7.1429e-02 + +DEAL::Mass matrix patch: +7.1429e-02 2.6620e-02 -2.6620e-02 1.1905e-02 0. 0. 0. +2.6620e-02 3.5714e-01 5.9524e-02 -2.6620e-02 0. 0. 0. +-2.6620e-02 5.9524e-02 3.5714e-01 2.6620e-02 0. 0. 0. +1.1905e-02 -2.6620e-02 2.6620e-02 1.4286e-01 2.6620e-02 -2.6620e-02 1.1905e-02 + 0. 0. 0. 2.6620e-02 3.5714e-01 5.9524e-02 -2.6620e-02 + 0. 0. 0. -2.6620e-02 5.9524e-02 3.5714e-01 2.6620e-02 + 0. 0. 0. 1.1905e-02 -2.6620e-02 2.6620e-02 7.1429e-02 + +DEAL::Laplace matrix for FE_Q<1>(3): +4.3333e+00 -1.6667e-01 -4.8784e+00 7.1175e-01 +-1.6667e-01 4.3333e+00 7.1175e-01 -4.8784e+00 +-4.8784e+00 7.1175e-01 8.3333e+00 -4.1667e+00 +7.1175e-01 -4.8784e+00 -4.1667e+00 8.3333e+00 + +DEAL::Laplace matrix for FE_Q<1>(3) renumbered with FEEvaluation: +4.3333e+00 -4.8784e+00 7.1175e-01 -1.6667e-01 +-4.8784e+00 8.3333e+00 -4.1667e+00 7.1175e-01 +7.1175e-01 -4.1667e+00 8.3333e+00 -4.8784e+00 +-1.6667e-01 7.1175e-01 -4.8784e+00 4.3333e+00 +