Reference documentation for deal.II version Git 3f1f337db3 2021-10-23 13:19:02 -0600
\(\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\}}\)
mg_transfer_prebuilt.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2003 - 2021 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 
17 #include <deal.II/base/function.h>
18 #include <deal.II/base/logstream.h>
19 
21 #include <deal.II/dofs/dof_tools.h>
22 
23 #include <deal.II/fe/fe.h>
24 
25 #include <deal.II/grid/tria.h>
27 
38 #include <deal.II/lac/vector.h>
39 
42 
43 #include <algorithm>
44 
46 
47 
48 template <typename VectorType>
50  const MGConstrainedDoFs &mg_c)
51 {
52  this->mg_constrained_dofs = &mg_c;
53 }
54 
55 
56 
57 template <typename VectorType>
58 void
60  const MGConstrainedDoFs &mg_c)
61 {
62  this->mg_constrained_dofs = &mg_c;
63 }
64 
65 
66 
67 template <typename VectorType>
68 void
70 {
72  prolongation_matrices.resize(0);
73  prolongation_sparsities.resize(0);
74  interface_dofs.resize(0);
75 }
76 
77 
78 
79 template <typename VectorType>
80 void
81 MGTransferPrebuilt<VectorType>::prolongate(const unsigned int to_level,
82  VectorType & dst,
83  const VectorType & src) const
84 {
85  Assert((to_level >= 1) && (to_level <= prolongation_matrices.size()),
86  ExcIndexRange(to_level, 1, prolongation_matrices.size() + 1));
87 
88  if (this->mg_constrained_dofs != nullptr &&
89  this->mg_constrained_dofs->get_user_constraint_matrix(to_level - 1)
90  .get_local_lines()
91  .size() > 0)
92  {
93  VectorType copy_src(src);
94  this->mg_constrained_dofs->get_user_constraint_matrix(to_level - 1)
95  .distribute(copy_src);
96  prolongation_matrices[to_level - 1]->vmult(dst, copy_src);
97  }
98  else
99  {
100  prolongation_matrices[to_level - 1]->vmult(dst, src);
101  }
102 }
103 
104 
105 
106 template <typename VectorType>
107 void
109  VectorType & dst,
110  const VectorType & src) const
111 {
112  Assert((from_level >= 1) && (from_level <= prolongation_matrices.size()),
113  ExcIndexRange(from_level, 1, prolongation_matrices.size() + 1));
114  (void)from_level;
115 
116  prolongation_matrices[from_level - 1]->Tvmult_add(dst, src);
117 }
118 
119 
120 namespace
121 {
126  void
127  replace(const MGConstrainedDoFs * mg_constrained_dofs,
128  const unsigned int level,
129  std::vector<types::global_dof_index> &dof_indices)
130  {
131  if (mg_constrained_dofs != nullptr &&
132  mg_constrained_dofs->get_level_constraints(level).n_constraints() > 0)
133  for (auto &ind : dof_indices)
134  if (mg_constrained_dofs->get_level_constraints(level)
136  {
137  Assert(mg_constrained_dofs->get_level_constraints(level)
139  ->size() == 1,
140  ExcInternalError());
141  ind = mg_constrained_dofs->get_level_constraints(level)
143  ->front()
144  .first;
145  }
146  }
147 } // namespace
148 
149 template <typename VectorType>
150 template <int dim, int spacedim>
151 void
153  const DoFHandler<dim, spacedim> &dof_handler)
154 {
155  Assert(dof_handler.has_level_dofs(),
156  ExcMessage(
157  "The underlying DoFHandler object has not had its "
158  "distribute_mg_dofs() function called, but this is a prerequisite "
159  "for multigrid transfers. You will need to call this function, "
160  "probably close to where you already call distribute_dofs()."));
161 
162  const unsigned int n_levels =
163  dof_handler.get_triangulation().n_global_levels();
164  const unsigned int dofs_per_cell = dof_handler.get_fe().n_dofs_per_cell();
165 
166  this->sizes.resize(n_levels);
167  for (unsigned int l = 0; l < n_levels; ++l)
168  this->sizes[l] = dof_handler.n_dofs(l);
169 
170  // reset the size of the array of
171  // matrices. call resize(0) first,
172  // in order to delete all elements
173  // and clear their memory. then
174  // repopulate these arrays
175  //
176  // note that on resize(0), the
177  // shared_ptr class takes care of
178  // deleting the object it points to
179  // by itself
180  prolongation_matrices.resize(0);
181  prolongation_sparsities.resize(0);
182  prolongation_matrices.reserve(n_levels - 1);
183  prolongation_sparsities.reserve(n_levels - 1);
184 
185  for (unsigned int i = 0; i < n_levels - 1; ++i)
186  {
187  prolongation_sparsities.emplace_back(
189  prolongation_matrices.emplace_back(
191  }
192 
193  // two fields which will store the
194  // indices of the multigrid dofs
195  // for a cell and one of its children
196  std::vector<types::global_dof_index> dof_indices_parent(dofs_per_cell);
197  std::vector<types::global_dof_index> dof_indices_child(dofs_per_cell);
198  std::vector<types::global_dof_index> entries(dofs_per_cell);
199 
200  // for each level: first build the sparsity
201  // pattern of the matrices and then build the
202  // matrices themselves. note that we only
203  // need to take care of cells on the coarser
204  // level which have children
205  for (unsigned int level = 0; level < n_levels - 1; ++level)
206  {
207  // reset the dimension of the structure. note that for the number of
208  // entries per row, the number of parent dofs coupling to a child dof is
209  // necessary. this, of course, is the number of degrees of freedom per
210  // cell
211  //
212  // increment dofs_per_cell since a useless diagonal element will be
213  // stored
214  IndexSet level_p1_relevant_dofs;
216  level + 1,
217  level_p1_relevant_dofs);
218  DynamicSparsityPattern dsp(this->sizes[level + 1],
219  this->sizes[level],
220  level_p1_relevant_dofs);
221  typename DoFHandler<dim>::cell_iterator cell,
222  endc = dof_handler.end(level);
223  for (cell = dof_handler.begin(level); cell != endc; ++cell)
224  if (cell->has_children() &&
225  (dof_handler.get_triangulation().locally_owned_subdomain() ==
227  cell->level_subdomain_id() ==
228  dof_handler.get_triangulation().locally_owned_subdomain()))
229  {
230  cell->get_mg_dof_indices(dof_indices_parent);
231 
232  replace(this->mg_constrained_dofs, level, dof_indices_parent);
233 
234  Assert(cell->n_children() ==
237  for (unsigned int child = 0; child < cell->n_children(); ++child)
238  {
239  // set an alias to the prolongation matrix for this child
240  const FullMatrix<double> &prolongation =
241  dof_handler.get_fe().get_prolongation_matrix(
242  child, cell->refinement_case());
243 
244  Assert(prolongation.n() != 0, ExcNoProlongation());
245 
246  cell->child(child)->get_mg_dof_indices(dof_indices_child);
247 
248  replace(this->mg_constrained_dofs,
249  level + 1,
250  dof_indices_child);
251 
252  // now tag the entries in the
253  // matrix which will be used
254  // for this pair of parent/child
255  for (unsigned int i = 0; i < dofs_per_cell; ++i)
256  {
257  entries.resize(0);
258  for (unsigned int j = 0; j < dofs_per_cell; ++j)
259  if (prolongation(i, j) != 0)
260  entries.push_back(dof_indices_parent[j]);
261  dsp.add_entries(dof_indices_child[i],
262  entries.begin(),
263  entries.end());
264  }
265  }
266  }
267 
268 #ifdef DEAL_II_WITH_MPI
270  VectorType>::requires_distributed_sparsity_pattern)
271  {
272  // Since PETSc matrices do not offer the functionality to fill up in-
273  // complete sparsity patterns on their own, the sparsity pattern must
274  // be manually distributed.
275 
276  // Distribute sparsity pattern
278  dsp,
279  dof_handler.locally_owned_mg_dofs(level + 1),
280  dof_handler.get_communicator(),
281  dsp.row_index_set());
282  }
283 #endif
284 
286  *prolongation_matrices[level],
287  *prolongation_sparsities[level],
288  level,
289  dsp,
290  dof_handler);
291  dsp.reinit(0, 0);
292 
293  // In the end, the entries in this object will only be real valued.
294  // Nevertheless, we have to take the underlying scalar type of the
295  // vector we want to use this class with. The global matrix the entries
296  // of this matrix are copied into has to be able to perform a
297  // matrix-vector multiplication and this is in general only implemented if
298  // the scalar type for matrix and vector is the same. Furthermore,
299  // copying entries between this local object and the global matrix is only
300  // implemented if the objects have the same scalar type.
302 
303  // now actually build the matrices
304  for (cell = dof_handler.begin(level); cell != endc; ++cell)
305  if (cell->has_children() &&
306  (dof_handler.get_triangulation().locally_owned_subdomain() ==
308  cell->level_subdomain_id() ==
309  dof_handler.get_triangulation().locally_owned_subdomain()))
310  {
311  cell->get_mg_dof_indices(dof_indices_parent);
312 
313  replace(this->mg_constrained_dofs, level, dof_indices_parent);
314 
315  Assert(cell->n_children() ==
318  for (unsigned int child = 0; child < cell->n_children(); ++child)
319  {
320  // set an alias to the prolongation matrix for this child
321  prolongation = dof_handler.get_fe().get_prolongation_matrix(
322  child, cell->refinement_case());
323 
324  if (this->mg_constrained_dofs != nullptr &&
325  this->mg_constrained_dofs->have_boundary_indices())
326  for (unsigned int j = 0; j < dofs_per_cell; ++j)
327  if (this->mg_constrained_dofs->is_boundary_index(
328  level, dof_indices_parent[j]))
329  for (unsigned int i = 0; i < dofs_per_cell; ++i)
330  prolongation(i, j) = 0.;
331 
332  cell->child(child)->get_mg_dof_indices(dof_indices_child);
333 
334  replace(this->mg_constrained_dofs,
335  level + 1,
336  dof_indices_child);
337 
338  // now set the entries in the matrix
339  for (unsigned int i = 0; i < dofs_per_cell; ++i)
340  prolongation_matrices[level]->set(dof_indices_child[i],
341  dofs_per_cell,
342  dof_indices_parent.data(),
343  &prolongation(i, 0),
344  true);
345  }
346  }
347  prolongation_matrices[level]->compress(VectorOperation::insert);
348  }
349 
350  this->fill_and_communicate_copy_indices(dof_handler);
351 }
352 
353 
354 
355 template <typename VectorType>
356 void
358 {
359  for (unsigned int level = 0; level < prolongation_matrices.size(); ++level)
360  {
361  os << "Level " << level << std::endl;
362  prolongation_matrices[level]->print(os);
363  os << std::endl;
364  }
365 }
366 
367 
368 
369 template <typename VectorType>
370 std::size_t
372 {
374  for (unsigned int i = 0; i < prolongation_matrices.size(); ++i)
375  result += prolongation_matrices[i]->memory_consumption() +
376  prolongation_sparsities[i]->memory_consumption();
377 
378  return result;
379 }
380 
381 
382 // explicit instantiation
383 #include "mg_transfer_prebuilt.inst"
384 
385 
void initialize_constraints(const MGConstrainedDoFs &mg_constrained_dofs)
MGTransferPrebuilt()=default
void add_entries(const size_type row, ForwardIterator begin, ForwardIterator end, const bool indices_are_unique_and_sorted=false)
const Triangulation< dim, spacedim > & get_triangulation() const
const types::subdomain_id invalid_subdomain_id
Definition: types.h:276
std::size_t memory_consumption() const
cell_iterator begin(const unsigned int level=0) const
cell_iterator end() const
MPI_Comm get_communicator() const
const AffineConstraints< double > & get_level_constraints(const unsigned int level) const
static void reinit(Matrix &matrix, Sparsity &sparsity, int level, const SparsityPatternType &sp, const DoFHandler< dim, spacedim > &)
Definition: mg_transfer.h:58
const FiniteElement< dim, spacedim > & get_fe(const unsigned int index=0) const
void extract_locally_relevant_level_dofs(const DoFHandler< dim, spacedim > &dof_handler, const unsigned int level, IndexSet &dof_set)
Definition: dof_tools.cc:1175
size_type n_constraints() const
virtual void restrict_and_add(const unsigned int from_level, VectorType &dst, const VectorType &src) const override
void print_matrices(std::ostream &os) const
static ::ExceptionBase & ExcMessage(std::string arg1)
std::size_t memory_consumption() const
size_type n() const
#define Assert(cond, exc)
Definition: exceptions.h:1461
void reinit(const size_type m, const size_type n, const IndexSet &rowset=IndexSet())
const IndexSet & row_index_set() const
void build(const DoFHandler< dim, spacedim > &dof_handler)
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:401
unsigned int level
Definition: grid_out.cc:4589
types::global_dof_index n_dofs() const
const IndexSet & locally_owned_mg_dofs(const unsigned int level) const
bool is_identity_constrained(const size_type line_n) const
virtual void prolongate(const unsigned int to_level, VectorType &dst, const VectorType &src) const override
const std::vector< std::pair< size_type, number > > * get_constraint_entries(const size_type line_n) const
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:400
static ::ExceptionBase & ExcIndexRange(std::size_t arg1, std::size_t arg2, std::size_t arg3)
static ::ExceptionBase & ExcNotImplemented()
bool has_level_dofs() const
void distribute_sparsity_pattern(DynamicSparsityPattern &dsp, const IndexSet &locally_owned_rows, const MPI_Comm &mpi_comm, const IndexSet &locally_relevant_rows)
std::enable_if< std::is_fundamental< T >::value, std::size_t >::type memory_consumption(const T &t)
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
static ::ExceptionBase & ExcInternalError()