Reference documentation for deal.II version Git 423cd11810 2020-05-28 01:07:25 +0200
\(\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 - 2020 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 #ifdef DEBUG
89  if (this->mg_constrained_dofs != nullptr)
90  Assert(this->mg_constrained_dofs->get_user_constraint_matrix(to_level - 1)
91  .get_local_lines()
92  .size() == 0,
94 #endif
95 
96  prolongation_matrices[to_level - 1]->vmult(dst, src);
97 }
98 
99 
100 
101 template <typename VectorType>
102 void
104  VectorType & dst,
105  const VectorType & src) const
106 {
107  Assert((from_level >= 1) && (from_level <= prolongation_matrices.size()),
108  ExcIndexRange(from_level, 1, prolongation_matrices.size() + 1));
109  (void)from_level;
110 
111  prolongation_matrices[from_level - 1]->Tvmult_add(dst, src);
112 }
113 
114 
115 namespace
116 {
121  void
122  replace(const MGConstrainedDoFs * mg_constrained_dofs,
123  const unsigned int level,
124  std::vector<types::global_dof_index> &dof_indices)
125  {
126  if (mg_constrained_dofs != nullptr &&
127  mg_constrained_dofs->get_level_constraints(level).n_constraints() > 0)
128  for (auto &ind : dof_indices)
129  if (mg_constrained_dofs->get_level_constraints(level)
131  {
132  Assert(mg_constrained_dofs->get_level_constraints(level)
134  ->size() == 1,
135  ExcInternalError());
136  ind = mg_constrained_dofs->get_level_constraints(level)
138  ->front()
139  .first;
140  }
141  }
142 } // namespace
143 
144 template <typename VectorType>
145 template <int dim, int spacedim>
146 void
148  const DoFHandler<dim, spacedim> &dof_handler)
149 {
150  const unsigned int n_levels =
151  dof_handler.get_triangulation().n_global_levels();
152  const unsigned int dofs_per_cell = dof_handler.get_fe().dofs_per_cell;
153 
154  this->sizes.resize(n_levels);
155  for (unsigned int l = 0; l < n_levels; ++l)
156  this->sizes[l] = dof_handler.n_dofs(l);
157 
158  // reset the size of the array of
159  // matrices. call resize(0) first,
160  // in order to delete all elements
161  // and clear their memory. then
162  // repopulate these arrays
163  //
164  // note that on resize(0), the
165  // shared_ptr class takes care of
166  // deleting the object it points to
167  // by itself
168  prolongation_matrices.resize(0);
169  prolongation_sparsities.resize(0);
170  prolongation_matrices.reserve(n_levels - 1);
171  prolongation_sparsities.reserve(n_levels - 1);
172 
173  for (unsigned int i = 0; i < n_levels - 1; ++i)
174  {
175  prolongation_sparsities.emplace_back(
177  prolongation_matrices.emplace_back(
179  }
180 
181  // two fields which will store the
182  // indices of the multigrid dofs
183  // for a cell and one of its children
184  std::vector<types::global_dof_index> dof_indices_parent(dofs_per_cell);
185  std::vector<types::global_dof_index> dof_indices_child(dofs_per_cell);
186  std::vector<types::global_dof_index> entries(dofs_per_cell);
187 
188  // for each level: first build the sparsity
189  // pattern of the matrices and then build the
190  // matrices themselves. note that we only
191  // need to take care of cells on the coarser
192  // level which have children
193  for (unsigned int level = 0; level < n_levels - 1; ++level)
194  {
195  // reset the dimension of the structure. note that for the number of
196  // entries per row, the number of parent dofs coupling to a child dof is
197  // necessary. this, of course, is the number of degrees of freedom per
198  // cell
199  //
200  // increment dofs_per_cell since a useless diagonal element will be
201  // stored
202  IndexSet level_p1_relevant_dofs;
204  level + 1,
205  level_p1_relevant_dofs);
206  DynamicSparsityPattern dsp(this->sizes[level + 1],
207  this->sizes[level],
208  level_p1_relevant_dofs);
209  typename DoFHandler<dim>::cell_iterator cell,
210  endc = dof_handler.end(level);
211  for (cell = dof_handler.begin(level); cell != endc; ++cell)
212  if (cell->has_children() &&
213  (dof_handler.get_triangulation().locally_owned_subdomain() ==
215  cell->level_subdomain_id() ==
216  dof_handler.get_triangulation().locally_owned_subdomain()))
217  {
218  cell->get_mg_dof_indices(dof_indices_parent);
219 
220  replace(this->mg_constrained_dofs, level, dof_indices_parent);
221 
222  Assert(cell->n_children() ==
225  for (unsigned int child = 0; child < cell->n_children(); ++child)
226  {
227  // set an alias to the prolongation matrix for this child
228  const FullMatrix<double> &prolongation =
229  dof_handler.get_fe().get_prolongation_matrix(
230  child, cell->refinement_case());
231 
232  Assert(prolongation.n() != 0, ExcNoProlongation());
233 
234  cell->child(child)->get_mg_dof_indices(dof_indices_child);
235 
236  replace(this->mg_constrained_dofs,
237  level + 1,
238  dof_indices_child);
239 
240  // now tag the entries in the
241  // matrix which will be used
242  // for this pair of parent/child
243  for (unsigned int i = 0; i < dofs_per_cell; ++i)
244  {
245  entries.resize(0);
246  for (unsigned int j = 0; j < dofs_per_cell; ++j)
247  if (prolongation(i, j) != 0)
248  entries.push_back(dof_indices_parent[j]);
249  dsp.add_entries(dof_indices_child[i],
250  entries.begin(),
251  entries.end());
252  }
253  }
254  }
255 
256 #ifdef DEAL_II_WITH_MPI
258  VectorType>::requires_distributed_sparsity_pattern)
259  {
260  // Since PETSc matrices do not offer the functionality to fill up in-
261  // complete sparsity patterns on their own, the sparsity pattern must
262  // be manually distributed.
263 
264  // Retrieve communicator from triangulation if it is parallel
266  dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
267  &(dof_handler.get_triangulation()));
268 
269  MPI_Comm communicator = dist_tria != nullptr ?
270  dist_tria->get_communicator() :
271  MPI_COMM_SELF;
272 
273  // Distribute sparsity pattern
275  dsp,
276  dof_handler.locally_owned_mg_dofs(level + 1),
277  communicator,
278  dsp.row_index_set());
279  }
280 #endif
281 
283  *prolongation_matrices[level],
284  *prolongation_sparsities[level],
285  level,
286  dsp,
287  dof_handler);
288  dsp.reinit(0, 0);
289 
290  // In the end, the entries in this object will only be real valued.
291  // Nevertheless, we have to take the underlying scalar type of the
292  // vector we want to use this class with. The global matrix the entries
293  // of this matrix are copied into has to be able to perform a
294  // matrix-vector multiplication and this is in general only implemented if
295  // the scalar type for matrix and vector is the same. Furthermore,
296  // copying entries between this local object and the global matrix is only
297  // implemented if the objects have the same scalar type.
299 
300  // now actually build the matrices
301  for (cell = dof_handler.begin(level); cell != endc; ++cell)
302  if (cell->has_children() &&
303  (dof_handler.get_triangulation().locally_owned_subdomain() ==
305  cell->level_subdomain_id() ==
306  dof_handler.get_triangulation().locally_owned_subdomain()))
307  {
308  cell->get_mg_dof_indices(dof_indices_parent);
309 
310  replace(this->mg_constrained_dofs, level, dof_indices_parent);
311 
312  Assert(cell->n_children() ==
315  for (unsigned int child = 0; child < cell->n_children(); ++child)
316  {
317  // set an alias to the prolongation matrix for this child
318  prolongation = dof_handler.get_fe().get_prolongation_matrix(
319  child, cell->refinement_case());
320 
321  if (this->mg_constrained_dofs != nullptr &&
322  this->mg_constrained_dofs->have_boundary_indices())
323  for (unsigned int j = 0; j < dofs_per_cell; ++j)
324  if (this->mg_constrained_dofs->is_boundary_index(
325  level, dof_indices_parent[j]))
326  for (unsigned int i = 0; i < dofs_per_cell; ++i)
327  prolongation(i, j) = 0.;
328 
329  cell->child(child)->get_mg_dof_indices(dof_indices_child);
330 
331  replace(this->mg_constrained_dofs,
332  level + 1,
333  dof_indices_child);
334 
335  // now set the entries in the matrix
336  for (unsigned int i = 0; i < dofs_per_cell; ++i)
337  prolongation_matrices[level]->set(dof_indices_child[i],
338  dofs_per_cell,
339  dof_indices_parent.data(),
340  &prolongation(i, 0),
341  true);
342  }
343  }
344  prolongation_matrices[level]->compress(VectorOperation::insert);
345  }
346 
347  this->fill_and_communicate_copy_indices(dof_handler);
348 }
349 
350 
351 
352 template <typename VectorType>
353 template <int dim, int spacedim>
354 void
356  const DoFHandler<dim, spacedim> &dof_handler)
357 {
358  build(dof_handler);
359 }
360 
361 
362 
363 template <typename VectorType>
364 void
366 {
367  for (unsigned int level = 0; level < prolongation_matrices.size(); ++level)
368  {
369  os << "Level " << level << std::endl;
370  prolongation_matrices[level]->print(os);
371  os << std::endl;
372  }
373 }
374 
375 
376 
377 template <typename VectorType>
378 std::size_t
380 {
382  for (unsigned int i = 0; i < prolongation_matrices.size(); ++i)
383  result += prolongation_matrices[i]->memory_consumption() +
384  prolongation_sparsities[i]->memory_consumption();
385 
386  return result;
387 }
388 
389 
390 // explicit instantiation
391 #include "mg_transfer_prebuilt.inst"
392 
393 
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:285
std::size_t memory_consumption() const
cell_iterator begin(const unsigned int level=0) const
Definition: dof_handler.cc:918
static void reinit(Matrix &matrix, Sparsity &sparsity, int level, const SparsityPatternType &sp, const DoFHandlerType &)
Definition: mg_transfer.h:58
cell_iterator end() const
Definition: dof_handler.cc:947
const AffineConstraints< double > & get_level_constraints(const unsigned int level) const
static ::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
const FiniteElement< dim, spacedim > & get_fe(const unsigned int index=0) const
void extract_locally_relevant_level_dofs(const DoFHandlerType &dof_handler, const unsigned int level, IndexSet &dof_set)
Definition: dof_tools.cc:1215
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
std::size_t memory_consumption() const
size_type n() const
#define Assert(cond, exc)
Definition: exceptions.h:1419
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:362
unsigned int level
Definition: grid_out.cc:4355
types::global_dof_index n_dofs() const
void build_matrices(const DoFHandler< dim, spacedim > &dof_handler)
const IndexSet & locally_owned_mg_dofs(const unsigned int level) const
bool is_identity_constrained(const size_type line_n) const
virtual const MPI_Comm & get_communicator() const
Definition: tria_base.cc:138
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:361
static ::ExceptionBase & ExcNotImplemented()
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()