Reference documentation for deal.II version 9.5.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
tensor_product_matrix_creator.h
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2022 - 2023 by the deal.II authors
4//
5// This file is part of the deal.II library.
6//
7// The deal.II library is free software; you can use it, redistribute
8// it, and/or modify it under the terms of the GNU Lesser General
9// Public License as published by the Free Software Foundation; either
10// version 2.1 of the License, or (at your option) any later version.
11// The full text of the license can be found in the file LICENSE.md at
12// the top level directory of deal.II.
13//
14// ---------------------------------------------------------------------
15
16#ifndef dealii_tensor_product_matrix_creator_h
17#define dealii_tensor_product_matrix_creator_h
18
19
20#include <deal.II/base/config.h>
21
23
25
26#include <deal.II/fe/fe.h>
27#include <deal.II/fe/fe_tools.h>
30
32#include <deal.II/grid/tria.h>
33
34#include <set>
35
37
38
45{
50 {
54 };
55
65 template <int dim, typename Number>
66 std::pair<std::array<FullMatrix<Number>, dim>,
67 std::array<FullMatrix<Number>, dim>>
69 const FiniteElement<1> & fe,
70 const Quadrature<1> & quadrature,
71 const ::ndarray<LaplaceBoundaryType, dim, 2> &boundary_ids,
72 const ::ndarray<double, dim, 3> & cell_extent,
73 const unsigned int n_overlap = 1);
74
79 template <int dim, typename Number>
80 std::pair<std::array<FullMatrix<Number>, dim>,
81 std::array<FullMatrix<Number>, dim>>
83 const typename Triangulation<dim>::cell_iterator &cell,
84 const std::set<types::boundary_id> & dirichlet_boundaries,
85 const std::set<types::boundary_id> & neumann_boundaries,
86 const FiniteElement<1> & fe,
87 const Quadrature<1> & quadrature,
88 const ::ndarray<double, dim, 3> & cell_extent,
89 const unsigned int n_overlap = 1);
90
91} // namespace TensorProductMatrixCreator
92
93
94
95/*----------------------- Inline functions ----------------------------------*/
96
97
99{
100 namespace internal
101 {
102 template <typename Number>
103 void
104 clear_row_and_column(const unsigned int n_dofs_1D_with_overlap,
105 const unsigned int n,
106 FullMatrix<Number> &matrix)
107 {
108 for (unsigned int i = 0; i < n_dofs_1D_with_overlap; ++i)
109 {
110 matrix[i][n] = 0.0;
111 matrix[n][i] = 0.0;
112 }
113 }
114
115 template <typename Number>
116 std::tuple<FullMatrix<Number>, FullMatrix<Number>, bool>
118 const FiniteElement<1> &fe,
119 const Quadrature<1> & quadrature)
120 {
123
124 DoFHandler<1> dof_handler(tria);
125 dof_handler.distribute_dofs(fe);
126
127 MappingQ1<1> mapping;
128
129 const unsigned int n_dofs_1D = fe.n_dofs_per_cell();
130
131 FullMatrix<Number> mass_matrix_reference(n_dofs_1D, n_dofs_1D);
132 FullMatrix<Number> derivative_matrix_reference(n_dofs_1D, n_dofs_1D);
133
134 FEValues<1> fe_values(mapping,
135 fe,
136 quadrature,
139
140 fe_values.reinit(tria.begin());
141
142 const auto lexicographic_to_hierarchic_numbering =
144 FETools::hierarchic_to_lexicographic_numbering<1>(
145 fe.tensor_degree()));
146
147 for (const unsigned int q_index : fe_values.quadrature_point_indices())
148 for (const unsigned int i : fe_values.dof_indices())
149 for (const unsigned int j : fe_values.dof_indices())
150 {
151 mass_matrix_reference(i, j) +=
152 (fe_values.shape_value(lexicographic_to_hierarchic_numbering[i],
153 q_index) *
154 fe_values.shape_value(lexicographic_to_hierarchic_numbering[j],
155 q_index) *
156 fe_values.JxW(q_index));
157
158 derivative_matrix_reference(i, j) +=
159 (fe_values.shape_grad(lexicographic_to_hierarchic_numbering[i],
160 q_index) *
161 fe_values.shape_grad(lexicographic_to_hierarchic_numbering[j],
162 q_index) *
163 fe_values.JxW(q_index));
164 }
165
166 return std::tuple<FullMatrix<Number>, FullMatrix<Number>, bool>{
167 mass_matrix_reference, derivative_matrix_reference, false};
168 }
169 } // namespace internal
170
171
172
173 template <int dim, typename Number>
174 std::pair<std::array<FullMatrix<Number>, dim>,
175 std::array<FullMatrix<Number>, dim>>
177 const FiniteElement<1> & fe,
178 const Quadrature<1> & quadrature,
179 const ::ndarray<LaplaceBoundaryType, dim, 2> &boundary_ids,
180 const ::ndarray<double, dim, 3> & cell_extent,
181 const unsigned int n_overlap)
182 {
183 // 1) create element mass and siffness matrix (without overlap)
184 const auto create_reference_mass_and_stiffness_matrices =
185 internal::create_reference_mass_and_stiffness_matrices<Number>(
186 fe, quadrature);
187
188 const auto &M_ref =
189 std::get<0>(create_reference_mass_and_stiffness_matrices);
190 const auto &K_ref =
191 std::get<1>(create_reference_mass_and_stiffness_matrices);
192 const auto &is_dg =
193 std::get<2>(create_reference_mass_and_stiffness_matrices);
194
195 AssertIndexRange(n_overlap, M_ref.n());
196 AssertIndexRange(0, n_overlap);
197 AssertThrow(is_dg == false, ExcNotImplemented());
198
199 // 2) loop over all dimensions and create 1d mass and stiffness
200 // matrices so that boundary conditions and overlap are considered
201
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;
204
205 std::array<FullMatrix<Number>, dim> Ms;
206 std::array<FullMatrix<Number>, dim> Ks;
207
208 for (unsigned int d = 0; d < dim; ++d)
209 {
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);
212
213 // inner cell
214 for (unsigned int i = 0; i < n_dofs_1D; ++i)
215 for (unsigned int j = 0; j < n_dofs_1D; ++j)
216 {
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];
221 }
222
223 // left neighbor or left boundary
224 if (boundary_ids[d][0] == LaplaceBoundaryType::internal_boundary)
225 {
226 // left neighbor
227 Assert(cell_extent[d][0] > 0.0, ExcInternalError());
228
229 for (unsigned int i = 0; i < n_overlap; ++i)
230 for (unsigned int j = 0; j < n_overlap; ++j)
231 {
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];
236 }
237 }
238 else
239 {
240 if (boundary_ids[d][0] == LaplaceBoundaryType::dirichlet)
241 {
242 // left DBC
243 const unsigned i0 = n_overlap - 1;
244 internal::clear_row_and_column(n_dofs_1D_with_overlap,
245 i0,
246 Ms[d]);
247 internal::clear_row_and_column(n_dofs_1D_with_overlap,
248 i0,
249 Ks[d]);
250 }
251 else if (boundary_ids[d][0] == LaplaceBoundaryType::neumann)
252 {
253 // left NBC -> nothing to do
254 }
255 else
256 {
258 }
259 }
260
261 // right neighbor or right boundary
262 if (boundary_ids[d][1] == LaplaceBoundaryType::internal_boundary)
263 {
264 Assert(cell_extent[d][2] > 0.0, ExcInternalError());
265
266 for (unsigned int i = 0; i < n_overlap; ++i)
267 for (unsigned int j = 0; j < n_overlap; ++j)
268 {
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];
273 }
274 }
275 else
276 {
277 if (boundary_ids[d][1] == LaplaceBoundaryType::dirichlet)
278 {
279 // right DBC
280 const unsigned i0 = n_overlap + n_dofs_1D - 2;
281 internal::clear_row_and_column(n_dofs_1D_with_overlap,
282 i0,
283 Ms[d]);
284 internal::clear_row_and_column(n_dofs_1D_with_overlap,
285 i0,
286 Ks[d]);
287 }
288 else if (boundary_ids[d][1] == LaplaceBoundaryType::neumann)
289 {
290 // right NBC -> nothing to do
291 }
292 else
293 {
295 }
296 }
297 }
298
299 return {Ms, Ks};
300 }
301
302
303 template <int dim, typename Number>
304 std::pair<std::array<FullMatrix<Number>, dim>,
305 std::array<FullMatrix<Number>, dim>>
307 const typename Triangulation<dim>::cell_iterator &cell,
308 const std::set<types::boundary_id> & dirichlet_boundaries,
309 const std::set<types::boundary_id> & neumann_boundaries,
310 const FiniteElement<1> & fe,
311 const Quadrature<1> & quadrature,
312 const ::ndarray<double, dim, 3> & cell_extent,
313 const unsigned int n_overlap)
314 {
316
317 for (unsigned int d = 0; d < dim; ++d)
318 {
319 // left neighbor or left boundary
320 if ((cell->at_boundary(2 * d) == false) ||
321 cell->has_periodic_neighbor(2 * d))
322 {
323 // left neighbor
324 Assert(cell_extent[d][0] > 0.0, ExcInternalError());
325
326 boundary_ids[d][0] = LaplaceBoundaryType::internal_boundary;
327 }
328 else
329 {
330 const auto bid = cell->face(2 * d)->boundary_id();
331 if (dirichlet_boundaries.find(bid) !=
332 dirichlet_boundaries.end() /*DBC*/)
333 {
334 // left DBC
335 boundary_ids[d][0] = LaplaceBoundaryType::dirichlet;
336 }
337 else if (neumann_boundaries.find(bid) !=
338 neumann_boundaries.end() /*NBC*/)
339 {
340 // left NBC
341 boundary_ids[d][0] = LaplaceBoundaryType::neumann;
342 }
343 else
344 {
346 }
347 }
348
349 // right neighbor or right boundary
350 if ((cell->at_boundary(2 * d + 1) == false) ||
351 cell->has_periodic_neighbor(2 * d + 1))
352 {
353 Assert(cell_extent[d][2] > 0.0, ExcInternalError());
354
355 boundary_ids[d][1] = LaplaceBoundaryType::internal_boundary;
356 }
357 else
358 {
359 const auto bid = cell->face(2 * d + 1)->boundary_id();
360 if (dirichlet_boundaries.find(bid) !=
361 dirichlet_boundaries.end() /*DBC*/)
362 {
363 // right DBC
364 boundary_ids[d][1] = LaplaceBoundaryType::dirichlet;
365 }
366 else if (neumann_boundaries.find(bid) !=
367 neumann_boundaries.end() /*NBC*/)
368 {
369 // right NBC
370 boundary_ids[d][1] = LaplaceBoundaryType::neumann;
371 }
372 else
373 {
375 }
376 }
377 }
378
379 return create_laplace_tensor_product_matrix<dim, Number>(
380 fe, quadrature, boundary_ids, cell_extent, n_overlap);
381 }
382
383} // namespace TensorProductMatrixCreator
384
385
386
388
389#endif
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
Definition config.h:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
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)
Definition utilities.h:1672
typename internal::ndarray::HelperArray< T, Ns... >::type ndarray
Definition ndarray.h:108
const ::Triangulation< dim, spacedim > & tria