]> https://gitweb.dealii.org/ - dealii.git/commitdiff
tests for 1D matrices
authorMichał Wichrowski <mtwichrowski@gmail.com>
Thu, 1 May 2025 18:25:58 +0000 (20:25 +0200)
committerMichał Wichrowski <mtwichrowski@gmail.com>
Thu, 1 May 2025 18:25:58 +0000 (20:25 +0200)
include/deal.II/numerics/tensor_product_matrix_creator.h
tests/numerics/tensor_product_creator_01.cc
tests/numerics/tensor_product_creator_01.output

index c1caa78c6a043fc89bc3d907095df739e674be01..60b672eb591064ce8f85d469395e7f9020b76724 100644 (file)
@@ -113,10 +113,11 @@ namespace TensorProductMatrixCreator
    * @return A FullMatrix<double> representing the assembled mass matrix.
    *
    */
-  FullMatrix<double>
+  template <typename Number = double>
+  FullMatrix<Number>
   create_1d_cell_mass_matrix(
     const FiniteElement<1>     &fe,
-    const double               &h,
+    const Number               &h,
     const std::pair<bool, bool> include_endpoints = {true, true},
     std::vector<unsigned int>   numbering = std::vector<unsigned int>());
 
@@ -140,10 +141,11 @@ namespace TensorProductMatrixCreator
    * numbering of the degrees of freedom. If empty, the
    * numbering of the finite element is used.
    */
-  FullMatrix<double>
+  template <typename Number = double>
+  FullMatrix<Number>
   create_1d_cell_laplace_matrix(
     const FiniteElement<1>     &fe,
-    const double               &h,
+    const Number               &h,
     const std::pair<bool, bool> include_endpoints = {true, true},
     std::vector<unsigned int>   numbering = std::vector<unsigned int>());
 
@@ -179,9 +181,10 @@ namespace TensorProductMatrixCreator
    * if the number of cells is above 10 an exception will be thrown.
    *
    */
-  FullMatrix<double>
+  template <typename Number = double>
+  FullMatrix<Number>
   create_1D_discretization_matrix(
-    FullMatrix<double>         &cell_matrix,
+    FullMatrix<Number>         &cell_matrix,
     const unsigned int         &n_cells,
     const unsigned int         &overlap,
     const std::pair<bool, bool> include_endpoints = {true, true});
@@ -203,11 +206,12 @@ namespace TensorProductMatrixCreator
    *
    * @return A full matrix representing the 1D ghost penalty matrix.
    */
-  FullMatrix<double>
+  template <typename Number = double>
+  FullMatrix<Number>
   create_1d_ghost_penalty_matrix(
     const FiniteElement<1> &fe,
-    const double            h,
-    std::vector<double>     coefficients = std::vector<double>());
+    const Number            h,
+    std::vector<Number>     coefficients = std::vector<Number>());
 
 
 
@@ -228,7 +232,8 @@ namespace TensorProductMatrixCreator
    *
    * @return A full matrix representing the 1D ghost penalty.
    */
-  FullMatrix<double>
+  template <typename Number = double>
+  FullMatrix<Number>
   create_1d_ghost_penalty_matrix(
     const std::vector<Polynomials::Polynomial<double>>
                       &polynomial_basis_derivative,
@@ -526,9 +531,57 @@ namespace TensorProductMatrixCreator
       fe, quadrature, boundary_ids, cell_extent, n_overlap);
   }
 
-  FullMatrix<double>
+  template <typename Number>
+  FullMatrix<Number>
+  create_1d_cell_mass_matrix(const FiniteElement<1>     &fe,
+                             const Number               &h,
+                             const std::pair<bool, bool> include_endpoints,
+                             std::vector<unsigned int>   numbering)
+  {
+    if (dynamic_cast<const FE_DGQ<1> *>(&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<Number> 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 <typename Number>
+  FullMatrix<Number>
   create_1d_cell_laplace_matrix(const FiniteElement<1>     &fe,
-                                const double               &h,
+                                const Number               &h,
                                 const std::pair<bool, bool> include_endpoints,
                                 std::vector<unsigned int>   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<double> cell_matrix(n_dofs_per_cell, n_dofs_per_cell);
+    FullMatrix<Number> 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<double>
-  create_1D_discretization_matrix(FullMatrix<double>         &cell_matrix,
+  template <typename Number>
+  FullMatrix<Number>
+  create_1D_discretization_matrix(FullMatrix<Number>         &cell_matrix,
                                   const unsigned int         &n_cells,
                                   const unsigned int         &overlap,
                                   const std::pair<bool, bool> include_endpoints)
@@ -600,7 +654,7 @@ namespace TensorProductMatrixCreator
     if (!include_endpoints.second)
       n_total_dofs -= 1;
 
-    FullMatrix<double> result_matrix(n_total_dofs, n_total_dofs);
+    FullMatrix<Number> 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<double>
+  template <typename Number>
+  FullMatrix<Number>
   create_1d_ghost_penalty_matrix(const FiniteElement<1> &fe,
-                                 const double            h,
-                                 std::vector<double>     coefficients)
+                                 const Number            h,
+                                 std::vector<Number>     coefficients)
   {
     Assert(dynamic_cast<const FE_Q<1> *>(&fe) != nullptr, ExcNotImplemented());
     Assert(h > 0, ExcMessage("Provided element size h is negative"));
@@ -683,13 +738,13 @@ namespace TensorProductMatrixCreator
       }
 
 
-    FullMatrix<double> penalty_matrix =
+    FullMatrix<Number> penalty_matrix =
       create_1d_ghost_penalty_matrix(polynomial_basis[1]);
     penalty_matrix *= coefficients[0];
 
     for (unsigned int k = 2; k < degree + 1; ++k)
       {
-        FullMatrix<double> kth_matrix =
+        FullMatrix<Number> 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<double>
+  template <typename Number>
+  FullMatrix<Number>
   create_1d_ghost_penalty_matrix(
     const std::vector<Polynomials::Polynomial<double>>
                       &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<double> penalty_matrix(n_total_dofs, n_total_dofs);
+    FullMatrix<Number> penalty_matrix(n_total_dofs, n_total_dofs);
 
     std::vector<double> values_left(n_dofs_per_cell);
     std::vector<double> values_right(n_dofs_per_cell);
index 1a9ca2942028e276197f3eb168a6387ac668dce3..bf30bc1227b9209bf9b575c3b69136ba1bfddde9 100644 (file)
 
 
 void
-print_matrix(FullMatrix<double> &matrix)
+print_matrix(const FullMatrix<double> &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 <int fe_degree>
 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
index f7c81d770729a33faef2d137060b7a0d832150d2..fa3194c2405591df11e0adfeeddf8f5b05c9d406 100644 (file)
@@ -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 
+

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.