/**
* @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.
/**
* @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
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,
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;
--- /dev/null
+#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