Reference documentation for deal.II version GIT 2f5445400b 2023-02-05 22:30:02+00:00
\(\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\}}\)
tensor_product_matrix_creator.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2022 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>
28 #include <deal.II/fe/fe_values.h>
29 #include <deal.II/fe/mapping_q1.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,
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 
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) +=
153  q_index) *
155  q_index) *
156  fe_values.JxW(q_index));
157 
158  derivative_matrix_reference(i, j) +=
160  q_index) *
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)
185  internal::create_reference_mass_and_stiffness_matrices<Number>(
186  fe, quadrature);
187 
188  const auto &M_ref =
190  const auto &K_ref =
192  const auto &is_dg =
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  {
257  AssertThrow(false, ExcNotImplemented());
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  {
294  AssertThrow(false, ExcNotImplemented());
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  {
345  AssertThrow(false, ExcNotImplemented());
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  {
374  AssertThrow(false, ExcNotImplemented());
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)
const double & shape_value(const unsigned int function_no, const unsigned int point_no) const
const Tensor< 1, spacedim > & shape_grad(const unsigned int function_no, const unsigned int quadrature_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 quadrature_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:461
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:462
static ::ExceptionBase & ExcInternalError()
#define Assert(cond, exc)
Definition: exceptions.h:1583
static ::ExceptionBase & ExcNotImplemented()
#define AssertIndexRange(index, range)
Definition: exceptions.h:1821
#define AssertThrow(cond, exc)
Definition: exceptions.h:1672
@ update_values
Shape function values.
@ update_JxW_values
Transformed quadrature weights.
@ update_gradients
Shape function gradients.
std::vector< unsigned int > lexicographic_to_hierarchic_numbering(unsigned int degree)
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)
Definition: utilities.h:1655
typename internal::ndarray::HelperArray< T, Ns... >::type ndarray
Definition: ndarray.h:108
const ::Triangulation< dim, spacedim > & tria