Add TensorProductMatrixCreator namespace

include/deal.II/lac/tensor_product_matrix.h

@@ -472,6 +472,61 @@ class TensorProductMatrixSymmetricSumCollection
 
 
 
+/**
+ * A namespace with functions that create input for certain standard matrices
+ * for the classes TensorProductMatrixSymmetricSum and
+ * TensorProductMatrixSymmetricSumCache.
+ */
+namespace TensorProductMatrixCreator
+{
+  /**
+   * Boundary type that can be used in create_laplace_tensor_product_matrix();
+   */
+  enum LaplaceBoundaryType
+  {
+    dirichlet,
+    neumann,
+    internal_boundary,
+  };
+
+  /**
+   * Create 1D mass matrix and 1D derivative matrix for a scalar
+   * consant-coefficient
+   * Laplacian cell for a @p dim dimensional Cartesian cell. Its boundary types
+   * can be specified with @p bids. The cell extend (including the cell extend
+   * of each neighbor) can be specified via @p cell_extend. With @p n_overlap, an
+   * overlap with neighboring cells can be specified. The default value is one,
+   * which correspond to all matrix enties restricted to the cell-local DoFs.
+   */
+  template <int dim, typename Number>
+  std::pair<std::array<FullMatrix<Number>, dim>,
+            std::array<FullMatrix<Number>, dim>>
+  create_laplace_tensor_product_matrix(
+    const FiniteElement<1> &                            fe,
+    const Quadrature<1> &                               quadrature,
+    const dealii::ndarray<LaplaceBoundaryType, dim, 2> &bids,
+    const dealii::ndarray<double, dim, 3> &             cell_extend,
+    const unsigned int                                  n_overlap = 1);
+
+  /**
+   * Same as above but the boundary IDs are extracted from the given @p cell
+   * and are mapped to the boundary type via the sets @p dbcs and @p nbcs.
+   */
+  template <int dim, typename Number>
+  std::pair<std::array<FullMatrix<Number>, dim>,
+            std::array<FullMatrix<Number>, dim>>
+  create_laplace_tensor_product_matrix(
+    const typename Triangulation<dim>::cell_iterator &cell,
+    const std::set<types::boundary_id> &              dbcs,
+    const std::set<types::boundary_id> &              nbcs,
+    const FiniteElement<1> &                          fe,
+    const Quadrature<1> &                             quadrature,
+    const dealii::ndarray<double, dim, 3> &           cell_extend,
+    const unsigned int                                n_overlap = 1);
+
+} // namespace TensorProductMatrixCreator
+
+
 /*----------------------- Inline functions ----------------------------------*/
 
 #ifndef DOXYGEN
@@ -1186,6 +1241,273 @@ TensorProductMatrixSymmetricSumCollection<dim, Number, n_rows_1d>::
 
 
 
+namespace TensorProductMatrixCreator
+{
+  namespace internal
+  {
+    template <typename Number>
+    std::tuple<FullMatrix<Number>, FullMatrix<Number>, bool>
+    create_referece_mass_and_striffness_matrices(
+      const FiniteElement<1> &fe,
+      const Quadrature<1> &   quadrature)
+    {
+      Triangulation<1> tria;
+      GridGenerator::hyper_cube(tria);
+
+      DoFHandler<1> dof_handler(tria);
+      dof_handler.distribute_dofs(fe);
+
+      MappingQ1<1> mapping;
+
+      const unsigned int n_dofs_1D = fe.n_dofs_per_cell();
+
+      FullMatrix<Number> mass_matrix_reference(n_dofs_1D, n_dofs_1D);
+      FullMatrix<Number> derivative_matrix_reference(n_dofs_1D, n_dofs_1D);
+
+      FEValues<1> fe_values(mapping,
+                            fe,
+                            quadrature,
+                            update_values | update_gradients |
+                              update_JxW_values);
+
+      fe_values.reinit(tria.begin());
+
+      const auto lexicographic_to_hierarchic_numbering =
+        Utilities::invert_permutation(
+          FETools::hierarchic_to_lexicographic_numbering<1>(
+            fe.tensor_degree()));
+
+      for (const unsigned int q_index : fe_values.quadrature_point_indices())
+        for (const unsigned int i : fe_values.dof_indices())
+          for (const unsigned int j : fe_values.dof_indices())
+            {
+              mass_matrix_reference(i, j) +=
+                (fe_values.shape_value(lexicographic_to_hierarchic_numbering[i],
+                                       q_index) *
+                 fe_values.shape_value(lexicographic_to_hierarchic_numbering[j],
+                                       q_index) *
+                 fe_values.JxW(q_index));
+
+              derivative_matrix_reference(i, j) +=
+                (fe_values.shape_grad(lexicographic_to_hierarchic_numbering[i],
+                                      q_index) *
+                 fe_values.shape_grad(lexicographic_to_hierarchic_numbering[j],
+                                      q_index) *
+                 fe_values.JxW(q_index));
+            }
+
+      return {mass_matrix_reference, derivative_matrix_reference, false};
+    }
+  } // namespace internal
+
+
+
+  template <int dim, typename Number>
+  std::pair<std::array<FullMatrix<Number>, dim>,
+            std::array<FullMatrix<Number>, dim>>
+  create_laplace_tensor_product_matrix(
+    const FiniteElement<1> &                            fe,
+    const Quadrature<1> &                               quadrature,
+    const dealii::ndarray<LaplaceBoundaryType, dim, 2> &bids,
+    const dealii::ndarray<double, dim, 3> &             cell_extend,
+    const unsigned int                                  n_overlap)
+  {
+    // 1) create element mass and siffness matrix (without overlap)
+    const auto [M_ref, K_ref, is_dg] =
+      internal::create_referece_mass_and_striffness_matrices<Number>(
+        fe, quadrature);
+
+    AssertIndexRange(n_overlap, M_ref.n());
+    AssertIndexRange(0, n_overlap);
+    AssertThrow(is_dg == false, ExcNotImplemented());
+
+    // 2) loop over all dimensions and create 1D mass and stiffness
+    // matrices so that boundary conditions and overlap are considered
+
+    const unsigned int n_dofs_1D              = M_ref.n();
+    const unsigned int n_dofs_1D_with_overlap = M_ref.n() - 2 + 2 * n_overlap;
+
+    std::array<FullMatrix<Number>, dim> Ms;
+    std::array<FullMatrix<Number>, dim> Ks;
+
+    const auto clear_row_and_column = [&](const unsigned int n, auto &matrix) {
+      for (unsigned int i = 0; i < n_dofs_1D_with_overlap; ++i)
+        {
+          matrix[i][n] = 0.0;
+          matrix[n][i] = 0.0;
+        }
+    };
+
+    for (unsigned int d = 0; d < dim; ++d)
+      {
+        FullMatrix<Number> M(n_dofs_1D_with_overlap, n_dofs_1D_with_overlap);
+        FullMatrix<Number> K(n_dofs_1D_with_overlap, n_dofs_1D_with_overlap);
+
+        // inner cell
+        for (unsigned int i = 0; i < n_dofs_1D; ++i)
+          for (unsigned int j = 0; j < n_dofs_1D; ++j)
+            {
+              const unsigned int i0 = i + n_overlap - 1;
+              const unsigned int j0 = j + n_overlap - 1;
+              M[i0][j0]             = M_ref[i][j] * cell_extend[d][1];
+              K[i0][j0]             = K_ref[i][j] / cell_extend[d][1];
+            }
+
+        // left neighbor or left boundary
+        if (bids[d][0] == LaplaceBoundaryType::internal_boundary)
+          {
+            // left neighbor
+            Assert(cell_extend[d][0] > 0.0, ExcInternalError());
+
+            for (unsigned int i = 0; i < n_overlap; ++i)
+              for (unsigned int j = 0; j < n_overlap; ++j)
+                {
+                  const unsigned int i0 = n_dofs_1D - n_overlap + i;
+                  const unsigned int j0 = n_dofs_1D - n_overlap + j;
+                  M[i][j] += M_ref[i0][j0] * cell_extend[d][0];
+                  K[i][j] += K_ref[i0][j0] / cell_extend[d][0];
+                }
+          }
+        else
+          {
+            if (bids[d][0] == LaplaceBoundaryType::dirichlet)
+              {
+                // left DBC
+                const unsigned i0 = n_overlap - 1;
+                clear_row_and_column(i0, M);
+                clear_row_and_column(i0, K);
+              }
+            else if (bids[d][0] == LaplaceBoundaryType::neumann)
+              {
+                // left NBC -> nothing to do
+              }
+            else
+              {
+                AssertThrow(false, ExcNotImplemented());
+              }
+          }
+
+        // reight neighbor or right boundary
+        if (bids[d][1] == LaplaceBoundaryType::internal_boundary)
+          {
+            Assert(cell_extend[d][2] > 0.0, ExcInternalError());
+
+            for (unsigned int i = 0; i < n_overlap; ++i)
+              for (unsigned int j = 0; j < n_overlap; ++j)
+                {
+                  const unsigned int i0 = n_overlap + n_dofs_1D + i - 2;
+                  const unsigned int j0 = n_overlap + n_dofs_1D + j - 2;
+                  M[i0][j0] += M_ref[i][j] * cell_extend[d][2];
+                  K[i0][j0] += K_ref[i][j] / cell_extend[d][2];
+                }
+          }
+        else
+          {
+            if (bids[d][1] == LaplaceBoundaryType::dirichlet)
+              {
+                // right DBC
+                const unsigned i0 = n_overlap + n_dofs_1D - 2;
+                clear_row_and_column(i0, M);
+                clear_row_and_column(i0, K);
+              }
+            else if (bids[d][1] == LaplaceBoundaryType::neumann)
+              {
+                // right NBC -> nothing to do
+              }
+            else
+              {
+                AssertThrow(false, ExcNotImplemented());
+              }
+          }
+
+        Ms[d] = M;
+        Ks[d] = K;
+      }
+
+    return {Ms, Ks};
+  }
+
+
+  template <int dim, typename Number>
+  std::pair<std::array<FullMatrix<Number>, dim>,
+            std::array<FullMatrix<Number>, dim>>
+  create_laplace_tensor_product_matrix(
+    const typename Triangulation<dim>::cell_iterator &cell,
+    const std::set<types::boundary_id> &              dbcs,
+    const std::set<types::boundary_id> &              nbcs,
+    const FiniteElement<1> &                          fe,
+    const Quadrature<1> &                             quadrature,
+    const dealii::ndarray<double, dim, 3> &           cell_extend,
+    const unsigned int                                n_overlap)
+  {
+    dealii::ndarray<LaplaceBoundaryType, dim, 2> bids;
+
+    for (unsigned int d = 0; d < dim; ++d)
+      {
+        // left neighbor or left boundary
+        if ((cell->at_boundary(2 * d) == false) ||
+            cell->has_periodic_neighbor(2 * d))
+          {
+            // left neighbor
+            Assert(cell_extend[d][0] > 0.0, ExcInternalError());
+
+            bids[d][0] = LaplaceBoundaryType::internal_boundary;
+          }
+        else
+          {
+            const auto bid = cell->face(2 * d)->boundary_id();
+            if (dbcs.contains(bid) /*DBC*/)
+              {
+                // left DBC
+                bids[d][0] = LaplaceBoundaryType::dirichlet;
+              }
+            else if (nbcs.contains(bid) /*NBC*/)
+              {
+                // left NBC
+                bids[d][0] = LaplaceBoundaryType::neumann;
+              }
+            else
+              {
+                AssertThrow(false, ExcNotImplemented());
+              }
+          }
+
+        // reight neighbor or right boundary
+        if ((cell->at_boundary(2 * d + 1) == false) ||
+            cell->has_periodic_neighbor(2 * d + 1))
+          {
+            Assert(cell_extend[d][2] > 0.0, ExcInternalError());
+
+            bids[d][1] = LaplaceBoundaryType::internal_boundary;
+          }
+        else
+          {
+            const auto bid = cell->face(2 * d + 1)->boundary_id();
+            if (dbcs.contains(bid) /*DBC*/)
+              {
+                // right DBC
+                bids[d][1] = LaplaceBoundaryType::dirichlet;
+              }
+            else if (nbcs.contains(bid) /*NBC*/)
+              {
+                // right NBC
+                bids[d][1] = LaplaceBoundaryType::neumann;
+              }
+            else
+              {
+                AssertThrow(false, ExcNotImplemented());
+              }
+          }
+      }
+
+    return create_laplace_tensor_product_matrix<dim, Number>(
+      fe, quadrature, bids, cell_extend, n_overlap);
+  }
+
+} // namespace TensorProductMatrixCreator
+
+
+
 #endif
 
 DEAL_II_NAMESPACE_CLOSE