16 #ifndef dealii_tensor_product_matrix_creator_h
17 #define dealii_tensor_product_matrix_creator_h
65 template <
int dim,
typename Number>
66 std::pair<std::array<FullMatrix<Number>, dim>,
67 std::array<FullMatrix<Number>, dim>>
71 const ::ndarray<LaplaceBoundaryType, dim, 2> &boundary_ids,
72 const ::ndarray<double, dim, 3> &cell_extent,
73 const unsigned int n_overlap = 1);
79 template <
int dim,
typename Number>
80 std::pair<std::array<FullMatrix<Number>, dim>,
81 std::array<FullMatrix<Number>, dim>>
84 const std::set<types::boundary_id> &dirichlet_boundaries,
85 const std::set<types::boundary_id> &neumann_boundaries,
88 const ::ndarray<double, dim, 3> &cell_extent,
89 const unsigned int n_overlap = 1);
102 template <
typename Number>
105 const unsigned int n,
108 for (
unsigned int i = 0; i < n_dofs_1D_with_overlap; ++i)
115 template <
typename Number>
144 FETools::hierarchic_to_lexicographic_numbering<1>(
148 for (
const unsigned int i : fe_values.
dof_indices())
149 for (
const unsigned int j : fe_values.
dof_indices())
151 mass_matrix_reference(i, j) +=
156 fe_values.
JxW(q_index));
158 derivative_matrix_reference(i, j) +=
163 fe_values.
JxW(q_index));
167 mass_matrix_reference, derivative_matrix_reference,
false};
173 template <
int dim,
typename Number>
174 std::pair<std::array<FullMatrix<Number>, dim>,
175 std::array<FullMatrix<Number>, dim>>
179 const ::ndarray<LaplaceBoundaryType, dim, 2> &boundary_ids,
180 const ::ndarray<double, dim, 3> &cell_extent,
181 const unsigned int n_overlap)
185 internal::create_reference_mass_and_stiffness_matrices<Number>(
202 const unsigned int n_dofs_1D = M_ref.n();
203 const unsigned int n_dofs_1D_with_overlap = M_ref.n() - 2 + 2 * n_overlap;
205 std::array<FullMatrix<Number>, dim> Ms;
206 std::array<FullMatrix<Number>, dim> Ks;
208 for (
unsigned int d = 0;
d < dim; ++
d)
210 Ms[
d].reinit(n_dofs_1D_with_overlap, n_dofs_1D_with_overlap);
211 Ks[
d].reinit(n_dofs_1D_with_overlap, n_dofs_1D_with_overlap);
214 for (
unsigned int i = 0; i < n_dofs_1D; ++i)
215 for (
unsigned int j = 0; j < n_dofs_1D; ++j)
217 const unsigned int i0 = i + n_overlap - 1;
218 const unsigned int j0 = j + n_overlap - 1;
219 Ms[
d][i0][j0] = M_ref[i][j] * cell_extent[
d][1];
220 Ks[
d][i0][j0] = K_ref[i][j] / cell_extent[
d][1];
229 for (
unsigned int i = 0; i < n_overlap; ++i)
230 for (
unsigned int j = 0; j < n_overlap; ++j)
232 const unsigned int i0 = n_dofs_1D - n_overlap + i;
233 const unsigned int j0 = n_dofs_1D - n_overlap + j;
234 Ms[
d][i][j] += M_ref[i0][j0] * cell_extent[
d][0];
235 Ks[
d][i][j] += K_ref[i0][j0] / cell_extent[
d][0];
243 const unsigned i0 = n_overlap - 1;
266 for (
unsigned int i = 0; i < n_overlap; ++i)
267 for (
unsigned int j = 0; j < n_overlap; ++j)
269 const unsigned int i0 = n_overlap + n_dofs_1D + i - 2;
270 const unsigned int j0 = n_overlap + n_dofs_1D + j - 2;
271 Ms[
d][i0][j0] += M_ref[i][j] * cell_extent[
d][2];
272 Ks[
d][i0][j0] += K_ref[i][j] / cell_extent[
d][2];
280 const unsigned i0 = n_overlap + n_dofs_1D - 2;
303 template <
int dim,
typename Number>
304 std::pair<std::array<FullMatrix<Number>, dim>,
305 std::array<FullMatrix<Number>, dim>>
308 const std::set<types::boundary_id> &dirichlet_boundaries,
309 const std::set<types::boundary_id> &neumann_boundaries,
312 const ::ndarray<double, dim, 3> &cell_extent,
313 const unsigned int n_overlap)
317 for (
unsigned int d = 0;
d < dim; ++
d)
320 if ((cell->at_boundary(2 *
d) ==
false) ||
321 cell->has_periodic_neighbor(2 *
d))
330 const auto bid = cell->face(2 *
d)->boundary_id();
331 if (dirichlet_boundaries.find(bid) !=
332 dirichlet_boundaries.end() )
337 else if (neumann_boundaries.find(bid) !=
338 neumann_boundaries.end() )
350 if ((cell->at_boundary(2 *
d + 1) ==
false) ||
351 cell->has_periodic_neighbor(2 *
d + 1))
359 const auto bid = cell->face(2 *
d + 1)->boundary_id();
360 if (dirichlet_boundaries.find(bid) !=
361 dirichlet_boundaries.end() )
366 else if (neumann_boundaries.find(bid) !=
367 neumann_boundaries.end() )
379 return create_laplace_tensor_product_matrix<dim, Number>(
380 fe, quadrature, boundary_ids, cell_extent, n_overlap);
void distribute_dofs(const FiniteElement< dim, spacedim > &fe)
const double & shape_value(const unsigned int i, const unsigned int q_point) const
const Tensor< 1, spacedim > & shape_grad(const unsigned int i, const unsigned int q_point) const
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
double JxW(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 & ExcInternalError()
#define Assert(cond, exc)
static ::ExceptionBase & ExcNotImplemented()
#define AssertIndexRange(index, range)
#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)
@ matrix
Contents is actually a matrix.
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
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
const ::Triangulation< dim, spacedim > & tria