]> https://gitweb.dealii.org/ - dealii.git/commitdiff
updated documentation moved implementation to .cc file
authorMichał Wichrowski <mtwichrowski@gmail.com>
Thu, 1 May 2025 16:45:48 +0000 (18:45 +0200)
committerMichał Wichrowski <mtwichrowski@gmail.com>
Thu, 1 May 2025 16:45:48 +0000 (18:45 +0200)
include/deal.II/numerics/tensor_product_matrix_creator.h
scratch/mwichro/dealii/src/numerics/tensor_product_matrix_creator_test.cc [new file with mode: 0644]

index be3f9bc15121c276984f666a94baea4a072360d8..5ea67a6ab672af8fc4eaf4225ec3281ce6171db5 100644 (file)
@@ -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<double>
-  create_1d_cell_mass_matrix(const FiniteElement<1>     &fe,
-                             const double               &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 double      &JxW             = h;
-    QGauss<1>          quadrature(degree + 1);
-
-    FullMatrix<double> 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<double>
   create_1d_cell_derivative_matrix(
     const FiniteElement<1>     &fe,
@@ -648,10 +602,11 @@ namespace TensorProductMatrixCreator
     FullMatrix<double> 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 (file)
index 0000000..fd1a39a
--- /dev/null
@@ -0,0 +1,70 @@
+#include <deal.II/base/exceptions.h>
+#include <deal.II/base/quadrature.h>
+
+#include <deal.II/fe/fe.h>
+
+#include <deal.II/lac/full_matrix.h>
+
+#include <deal.II/numerics/tensor_product_matrix_creator.h>
+
+#include <numeric>
+
+
+DEAL_II_NAMESPACE_OPEN
+
+
+namespace TensorProductMatrixCreator
+{
+
+
+  FullMatrix<double>
+  create_1d_cell_mass_matrix(const FiniteElement<1>     &fe,
+                             const double               &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 double      &JxW             = h;
+    QGauss<1>          quadrature(degree + 1);
+
+    FullMatrix<double> 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

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.