15#ifndef dealii_tensor_product_matrix_creator_h
16#define dealii_tensor_product_matrix_creator_h
64 template <
int dim,
typename Number>
65 std::pair<std::array<FullMatrix<Number>, dim>,
66 std::array<FullMatrix<Number>, dim>>
70 const ::ndarray<LaplaceBoundaryType, dim, 2> &boundary_ids,
71 const ::ndarray<double, dim, 3> &cell_extent,
72 const unsigned int n_overlap = 1);
78 template <
int dim,
typename Number>
79 std::pair<std::array<FullMatrix<Number>, dim>,
80 std::array<FullMatrix<Number>, dim>>
83 const std::set<types::boundary_id> &dirichlet_boundaries,
84 const std::set<types::boundary_id> &neumann_boundaries,
87 const ::ndarray<double, dim, 3> &cell_extent,
88 const unsigned int n_overlap = 1);
101 template <
typename Number>
104 const unsigned int n,
107 for (
unsigned int i = 0; i < n_dofs_1D_with_overlap; ++i)
114 template <
typename Number>
141 const auto lexicographic_to_hierarchic_numbering =
143 FETools::hierarchic_to_lexicographic_numbering<1>(
147 for (
const unsigned int i : fe_values.
dof_indices())
148 for (
const unsigned int j : fe_values.
dof_indices())
150 mass_matrix_reference(i, j) +=
151 (fe_values.
shape_value(lexicographic_to_hierarchic_numbering[i],
153 fe_values.
shape_value(lexicographic_to_hierarchic_numbering[j],
155 fe_values.
JxW(q_index));
157 derivative_matrix_reference(i, j) +=
158 (fe_values.
shape_grad(lexicographic_to_hierarchic_numbering[i],
160 fe_values.
shape_grad(lexicographic_to_hierarchic_numbering[j],
162 fe_values.
JxW(q_index));
166 mass_matrix_reference, derivative_matrix_reference,
false};
172 template <
int dim,
typename Number>
173 std::pair<std::array<FullMatrix<Number>, dim>,
174 std::array<FullMatrix<Number>, dim>>
178 const ::ndarray<LaplaceBoundaryType, dim, 2> &boundary_ids,
179 const ::ndarray<double, dim, 3> &cell_extent,
180 const unsigned int n_overlap)
183 const auto create_reference_mass_and_stiffness_matrices =
184 internal::create_reference_mass_and_stiffness_matrices<Number>(
188 std::get<0>(create_reference_mass_and_stiffness_matrices);
190 std::get<1>(create_reference_mass_and_stiffness_matrices);
192 std::get<2>(create_reference_mass_and_stiffness_matrices);
201 const unsigned int n_dofs_1D = M_ref.n();
202 const unsigned int n_dofs_1D_with_overlap = M_ref.n() - 2 + 2 * n_overlap;
204 std::array<FullMatrix<Number>, dim> Ms;
205 std::array<FullMatrix<Number>, dim> Ks;
207 for (
unsigned int d = 0; d < dim; ++d)
209 Ms[d].reinit(n_dofs_1D_with_overlap, n_dofs_1D_with_overlap);
210 Ks[d].reinit(n_dofs_1D_with_overlap, n_dofs_1D_with_overlap);
213 for (
unsigned int i = 0; i < n_dofs_1D; ++i)
214 for (
unsigned int j = 0; j < n_dofs_1D; ++j)
216 const unsigned int i0 = i + n_overlap - 1;
217 const unsigned int j0 = j + n_overlap - 1;
218 Ms[d][i0][j0] = M_ref[i][j] * cell_extent[d][1];
219 Ks[d][i0][j0] = K_ref[i][j] / cell_extent[d][1];
228 for (
unsigned int i = 0; i < n_overlap; ++i)
229 for (
unsigned int j = 0; j < n_overlap; ++j)
231 const unsigned int i0 = n_dofs_1D - n_overlap + i;
232 const unsigned int j0 = n_dofs_1D - n_overlap + j;
233 Ms[d][i][j] += M_ref[i0][j0] * cell_extent[d][0];
234 Ks[d][i][j] += K_ref[i0][j0] / cell_extent[d][0];
242 const unsigned i0 = n_overlap - 1;
265 for (
unsigned int i = 0; i < n_overlap; ++i)
266 for (
unsigned int j = 0; j < n_overlap; ++j)
268 const unsigned int i0 = n_overlap + n_dofs_1D + i - 2;
269 const unsigned int j0 = n_overlap + n_dofs_1D + j - 2;
270 Ms[d][i0][j0] += M_ref[i][j] * cell_extent[d][2];
271 Ks[d][i0][j0] += K_ref[i][j] / cell_extent[d][2];
279 const unsigned i0 = n_overlap + n_dofs_1D - 2;
302 template <
int dim,
typename Number>
303 std::pair<std::array<FullMatrix<Number>, dim>,
304 std::array<FullMatrix<Number>, dim>>
307 const std::set<types::boundary_id> &dirichlet_boundaries,
308 const std::set<types::boundary_id> &neumann_boundaries,
311 const ::ndarray<double, dim, 3> &cell_extent,
312 const unsigned int n_overlap)
316 for (
unsigned int d = 0; d < dim; ++d)
319 if ((cell->at_boundary(2 * d) ==
false) ||
320 cell->has_periodic_neighbor(2 * d))
329 const auto bid = cell->face(2 * d)->boundary_id();
330 if (dirichlet_boundaries.find(bid) !=
331 dirichlet_boundaries.end() )
336 else if (neumann_boundaries.find(bid) !=
337 neumann_boundaries.end() )
349 if ((cell->at_boundary(2 * d + 1) ==
false) ||
350 cell->has_periodic_neighbor(2 * d + 1))
358 const auto bid = cell->face(2 * d + 1)->boundary_id();
359 if (dirichlet_boundaries.find(bid) !=
360 dirichlet_boundaries.end() )
365 else if (neumann_boundaries.find(bid) !=
366 neumann_boundaries.end() )
378 return create_laplace_tensor_product_matrix<dim, Number>(
379 fe, quadrature, boundary_ids, cell_extent, n_overlap);
void distribute_dofs(const FiniteElement< dim, spacedim > &fe)
std_cxx20::ranges::iota_view< unsigned int, unsigned int > dof_indices() const
std_cxx20::ranges::iota_view< unsigned int, unsigned int > quadrature_point_indices() const
const Tensor< 1, spacedim > & shape_grad(const unsigned int i, const unsigned int q_point) const
double JxW(const unsigned int q_point) const
const double & shape_value(const unsigned int i, const unsigned int q_point) const
void reinit(const TriaIterator< DoFCellAccessor< dim, spacedim, level_dof_access > > &cell)
unsigned int n_dofs_per_cell() const
unsigned int tensor_degree() const
cell_iterator begin(const unsigned int level=0) const
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_CLOSE
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
#define AssertThrow(cond, exc)
@ update_values
Shape function values.
@ update_JxW_values
Transformed quadrature weights.
@ update_gradients
Shape function gradients.
void hyper_cube(Triangulation< dim, spacedim > &tria, const double left=0., const double right=1., const bool colorize=false)
void clear_row_and_column(const unsigned int n_dofs_1D_with_overlap, const unsigned int n, FullMatrix< Number > &matrix)
std::tuple< FullMatrix< Number >, FullMatrix< Number >, bool > create_reference_mass_and_stiffness_matrices(const FiniteElement< 1 > &fe, const Quadrature< 1 > &quadrature)
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 ::ndarray< LaplaceBoundaryType, dim, 2 > &boundary_ids, const ::ndarray< double, dim, 3 > &cell_extent, const unsigned int n_overlap=1)
std::vector< Integer > invert_permutation(const std::vector< Integer > &permutation)
typename internal::ndarray::HelperArray< T, Ns... >::type ndarray