Reference documentation for deal.II version 9.6.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// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2022 - 2023 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15#ifndef dealii_tensor_product_matrix_creator_h
16#define dealii_tensor_product_matrix_creator_h
17
18
19#include <deal.II/base/config.h>
20
22
24
25#include <deal.II/fe/fe.h>
26#include <deal.II/fe/fe_tools.h>
29
31#include <deal.II/grid/tria.h>
32
33#include <set>
34
36
37
44{
54
64 template <int dim, typename Number>
65 std::pair<std::array<FullMatrix<Number>, dim>,
66 std::array<FullMatrix<Number>, dim>>
68 const FiniteElement<1> &fe,
69 const Quadrature<1> &quadrature,
70 const ::ndarray<LaplaceBoundaryType, dim, 2> &boundary_ids,
71 const ::ndarray<double, dim, 3> &cell_extent,
72 const unsigned int n_overlap = 1);
73
78 template <int dim, typename Number>
79 std::pair<std::array<FullMatrix<Number>, dim>,
80 std::array<FullMatrix<Number>, dim>>
82 const typename Triangulation<dim>::cell_iterator &cell,
83 const std::set<types::boundary_id> &dirichlet_boundaries,
84 const std::set<types::boundary_id> &neumann_boundaries,
85 const FiniteElement<1> &fe,
86 const Quadrature<1> &quadrature,
87 const ::ndarray<double, dim, 3> &cell_extent,
88 const unsigned int n_overlap = 1);
89
90} // namespace TensorProductMatrixCreator
91
92
93
94/*----------------------- Inline functions ----------------------------------*/
95
96
98{
99 namespace internal
100 {
101 template <typename Number>
102 void
103 clear_row_and_column(const unsigned int n_dofs_1D_with_overlap,
104 const unsigned int n,
105 FullMatrix<Number> &matrix)
106 {
107 for (unsigned int i = 0; i < n_dofs_1D_with_overlap; ++i)
108 {
109 matrix[i][n] = 0.0;
110 matrix[n][i] = 0.0;
111 }
112 }
113
114 template <typename Number>
115 std::tuple<FullMatrix<Number>, FullMatrix<Number>, bool>
117 const FiniteElement<1> &fe,
118 const Quadrature<1> &quadrature)
119 {
120 Triangulation<1> tria;
122
123 DoFHandler<1> dof_handler(tria);
124 dof_handler.distribute_dofs(fe);
125
126 MappingQ1<1> mapping;
127
128 const unsigned int n_dofs_1D = fe.n_dofs_per_cell();
129
130 FullMatrix<Number> mass_matrix_reference(n_dofs_1D, n_dofs_1D);
131 FullMatrix<Number> derivative_matrix_reference(n_dofs_1D, n_dofs_1D);
132
133 FEValues<1> fe_values(mapping,
134 fe,
135 quadrature,
138
139 fe_values.reinit(tria.begin());
140
141 const auto lexicographic_to_hierarchic_numbering =
144 fe.tensor_degree()));
145
146 for (const unsigned int q_index : fe_values.quadrature_point_indices())
147 for (const unsigned int i : fe_values.dof_indices())
148 for (const unsigned int j : fe_values.dof_indices())
149 {
150 mass_matrix_reference(i, j) +=
151 (fe_values.shape_value(lexicographic_to_hierarchic_numbering[i],
152 q_index) *
153 fe_values.shape_value(lexicographic_to_hierarchic_numbering[j],
154 q_index) *
155 fe_values.JxW(q_index));
156
157 derivative_matrix_reference(i, j) +=
158 (fe_values.shape_grad(lexicographic_to_hierarchic_numbering[i],
159 q_index) *
160 fe_values.shape_grad(lexicographic_to_hierarchic_numbering[j],
161 q_index) *
162 fe_values.JxW(q_index));
163 }
164
165 return std::tuple<FullMatrix<Number>, FullMatrix<Number>, bool>{
166 mass_matrix_reference, derivative_matrix_reference, false};
167 }
168 } // namespace internal
169
170
171
172 template <int dim, typename Number>
173 std::pair<std::array<FullMatrix<Number>, dim>,
174 std::array<FullMatrix<Number>, dim>>
176 const FiniteElement<1> &fe,
177 const Quadrature<1> &quadrature,
178 const ::ndarray<LaplaceBoundaryType, dim, 2> &boundary_ids,
179 const ::ndarray<double, dim, 3> &cell_extent,
180 const unsigned int n_overlap)
181 {
182 // 1) create element mass and siffness matrix (without overlap)
183 const auto create_reference_mass_and_stiffness_matrices =
185 fe, quadrature);
186
187 const auto &M_ref =
188 std::get<0>(create_reference_mass_and_stiffness_matrices);
189 const auto &K_ref =
190 std::get<1>(create_reference_mass_and_stiffness_matrices);
191 const auto &is_dg =
192 std::get<2>(create_reference_mass_and_stiffness_matrices);
193
194 AssertIndexRange(n_overlap, M_ref.n());
195 AssertIndexRange(0, n_overlap);
196 AssertThrow(is_dg == false, ExcNotImplemented());
197
198 // 2) loop over all dimensions and create 1d mass and stiffness
199 // matrices so that boundary conditions and overlap are considered
200
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;
203
204 std::array<FullMatrix<Number>, dim> Ms;
205 std::array<FullMatrix<Number>, dim> Ks;
206
207 for (unsigned int d = 0; d < dim; ++d)
208 {
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);
211
212 // inner cell
213 for (unsigned int i = 0; i < n_dofs_1D; ++i)
214 for (unsigned int j = 0; j < n_dofs_1D; ++j)
215 {
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];
220 }
221
222 // left neighbor or left boundary
223 if (boundary_ids[d][0] == LaplaceBoundaryType::internal_boundary)
224 {
225 // left neighbor
226 Assert(cell_extent[d][0] > 0.0, ExcInternalError());
227
228 for (unsigned int i = 0; i < n_overlap; ++i)
229 for (unsigned int j = 0; j < n_overlap; ++j)
230 {
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];
235 }
236 }
237 else
238 {
239 if (boundary_ids[d][0] == LaplaceBoundaryType::dirichlet)
240 {
241 // left DBC
242 const unsigned i0 = n_overlap - 1;
243 internal::clear_row_and_column(n_dofs_1D_with_overlap,
244 i0,
245 Ms[d]);
246 internal::clear_row_and_column(n_dofs_1D_with_overlap,
247 i0,
248 Ks[d]);
249 }
250 else if (boundary_ids[d][0] == LaplaceBoundaryType::neumann)
251 {
252 // left NBC -> nothing to do
253 }
254 else
255 {
257 }
258 }
259
260 // right neighbor or right boundary
261 if (boundary_ids[d][1] == LaplaceBoundaryType::internal_boundary)
262 {
263 Assert(cell_extent[d][2] > 0.0, ExcInternalError());
264
265 for (unsigned int i = 0; i < n_overlap; ++i)
266 for (unsigned int j = 0; j < n_overlap; ++j)
267 {
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];
272 }
273 }
274 else
275 {
276 if (boundary_ids[d][1] == LaplaceBoundaryType::dirichlet)
277 {
278 // right DBC
279 const unsigned i0 = n_overlap + n_dofs_1D - 2;
280 internal::clear_row_and_column(n_dofs_1D_with_overlap,
281 i0,
282 Ms[d]);
283 internal::clear_row_and_column(n_dofs_1D_with_overlap,
284 i0,
285 Ks[d]);
286 }
287 else if (boundary_ids[d][1] == LaplaceBoundaryType::neumann)
288 {
289 // right NBC -> nothing to do
290 }
291 else
292 {
294 }
295 }
296 }
297
298 return {Ms, Ks};
299 }
300
301
302 template <int dim, typename Number>
303 std::pair<std::array<FullMatrix<Number>, dim>,
304 std::array<FullMatrix<Number>, dim>>
306 const typename Triangulation<dim>::cell_iterator &cell,
307 const std::set<types::boundary_id> &dirichlet_boundaries,
308 const std::set<types::boundary_id> &neumann_boundaries,
309 const FiniteElement<1> &fe,
310 const Quadrature<1> &quadrature,
311 const ::ndarray<double, dim, 3> &cell_extent,
312 const unsigned int n_overlap)
313 {
315
316 for (unsigned int d = 0; d < dim; ++d)
317 {
318 // left neighbor or left boundary
319 if ((cell->at_boundary(2 * d) == false) ||
320 cell->has_periodic_neighbor(2 * d))
321 {
322 // left neighbor
323 Assert(cell_extent[d][0] > 0.0, ExcInternalError());
324
325 boundary_ids[d][0] = LaplaceBoundaryType::internal_boundary;
326 }
327 else
328 {
329 const auto bid = cell->face(2 * d)->boundary_id();
330 if (dirichlet_boundaries.find(bid) !=
331 dirichlet_boundaries.end() /*DBC*/)
332 {
333 // left DBC
334 boundary_ids[d][0] = LaplaceBoundaryType::dirichlet;
335 }
336 else if (neumann_boundaries.find(bid) !=
337 neumann_boundaries.end() /*NBC*/)
338 {
339 // left NBC
340 boundary_ids[d][0] = LaplaceBoundaryType::neumann;
341 }
342 else
343 {
345 }
346 }
347
348 // right neighbor or right boundary
349 if ((cell->at_boundary(2 * d + 1) == false) ||
350 cell->has_periodic_neighbor(2 * d + 1))
351 {
352 Assert(cell_extent[d][2] > 0.0, ExcInternalError());
353
354 boundary_ids[d][1] = LaplaceBoundaryType::internal_boundary;
355 }
356 else
357 {
358 const auto bid = cell->face(2 * d + 1)->boundary_id();
359 if (dirichlet_boundaries.find(bid) !=
360 dirichlet_boundaries.end() /*DBC*/)
361 {
362 // right DBC
363 boundary_ids[d][1] = LaplaceBoundaryType::dirichlet;
364 }
365 else if (neumann_boundaries.find(bid) !=
366 neumann_boundaries.end() /*NBC*/)
367 {
368 // right NBC
369 boundary_ids[d][1] = LaplaceBoundaryType::neumann;
370 }
371 else
372 {
374 }
375 }
376 }
377
379 fe, quadrature, boundary_ids, cell_extent, n_overlap);
380 }
381
382} // namespace TensorProductMatrixCreator
383
384
385
387
388#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:503
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:504
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.
std::vector< unsigned int > hierarchic_to_lexicographic_numbering(unsigned int degree)
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:1699
typename internal::ndarray::HelperArray< T, Ns... >::type ndarray
Definition ndarray.h:107