Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
mg_transfer_prebuilt.cc
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2003 - 2019 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 
20 #include <deal.II/dofs/dof_accessor.h>
21 #include <deal.II/dofs/dof_tools.h>
22 
23 #include <deal.II/fe/fe.h>
24 
25 #include <deal.II/grid/tria.h>
26 #include <deal.II/grid/tria_iterator.h>
27 
28 #include <deal.II/lac/block_sparse_matrix.h>
29 #include <deal.II/lac/block_vector.h>
30 #include <deal.II/lac/dynamic_sparsity_pattern.h>
31 #include <deal.II/lac/la_parallel_block_vector.h>
32 #include <deal.II/lac/la_parallel_vector.h>
33 #include <deal.II/lac/sparse_matrix.h>
34 #include <deal.II/lac/sparsity_tools.h>
35 #include <deal.II/lac/trilinos_epetra_vector.h>
36 #include <deal.II/lac/trilinos_parallel_block_vector.h>
37 #include <deal.II/lac/trilinos_vector.h>
38 #include <deal.II/lac/vector.h>
39 
40 #include <deal.II/multigrid/mg_tools.h>
41 #include <deal.II/multigrid/mg_transfer.h>
42 
43 #include <algorithm>
44 
45 DEAL_II_NAMESPACE_OPEN
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>
59  const AffineConstraints<double> & /*c*/,
60  const MGConstrainedDoFs &mg_c)
61 {
62  this->mg_constrained_dofs = &mg_c;
63 }
64 
65 
66 
67 template <typename VectorType>
68 void
70  const MGConstrainedDoFs &mg_c)
71 {
72  this->mg_constrained_dofs = &mg_c;
73 }
74 
75 
76 
77 template <typename VectorType>
78 void
80  const AffineConstraints<double> & /*c*/,
81  const MGConstrainedDoFs &mg_c)
82 {
83  initialize_constraints(mg_c);
84 }
85 
86 
87 
88 template <typename VectorType>
89 void
91 {
93  prolongation_matrices.resize(0);
94  prolongation_sparsities.resize(0);
95  interface_dofs.resize(0);
96 }
97 
98 
99 
100 template <typename VectorType>
101 void
102 MGTransferPrebuilt<VectorType>::prolongate(const unsigned int to_level,
103  VectorType & dst,
104  const VectorType & src) const
105 {
106  Assert((to_level >= 1) && (to_level <= prolongation_matrices.size()),
107  ExcIndexRange(to_level, 1, prolongation_matrices.size() + 1));
108 
109  prolongation_matrices[to_level - 1]->vmult(dst, src);
110 }
111 
112 
113 
114 template <typename VectorType>
115 void
117  VectorType & dst,
118  const VectorType & src) const
119 {
120  Assert((from_level >= 1) && (from_level <= prolongation_matrices.size()),
121  ExcIndexRange(from_level, 1, prolongation_matrices.size() + 1));
122  (void)from_level;
123 
124  prolongation_matrices[from_level - 1]->Tvmult_add(dst, src);
125 }
126 
127 
128 namespace
129 {
134  void
135  replace(const MGConstrainedDoFs * mg_constrained_dofs,
136  const unsigned int level,
137  std::vector<types::global_dof_index> &dof_indices)
138  {
139  if (mg_constrained_dofs != nullptr &&
140  mg_constrained_dofs->get_level_constraints(level).n_constraints() > 0)
141  for (auto &ind : dof_indices)
142  if (mg_constrained_dofs->get_level_constraints(level)
144  {
145  Assert(mg_constrained_dofs->get_level_constraints(level)
147  ->size() == 1,
148  ExcInternalError());
149  ind = mg_constrained_dofs->get_level_constraints(level)
151  ->front()
152  .first;
153  }
154  }
155 } // namespace
156 
157 template <typename VectorType>
158 template <int dim, int spacedim>
159 void
161  const DoFHandler<dim, spacedim> &mg_dof)
162 {
163  const unsigned int n_levels = mg_dof.get_triangulation().n_global_levels();
164  const unsigned int dofs_per_cell = mg_dof.get_fe().dofs_per_cell;
165 
166  this->sizes.resize(n_levels);
167  for (unsigned int l = 0; l < n_levels; ++l)
168  this->sizes[l] = mg_dof.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, endc = mg_dof.end(level);
222  for (cell = mg_dof.begin(level); cell != endc; ++cell)
223  if (cell->has_children() &&
224  (mg_dof.get_triangulation().locally_owned_subdomain() ==
226  cell->level_subdomain_id() ==
227  mg_dof.get_triangulation().locally_owned_subdomain()))
228  {
229  cell->get_mg_dof_indices(dof_indices_parent);
230 
231  replace(this->mg_constrained_dofs, level, dof_indices_parent);
232 
233  Assert(cell->n_children() ==
236  for (unsigned int child = 0; child < cell->n_children(); ++child)
237  {
238  // set an alias to the prolongation matrix for this child
239  const FullMatrix<double> &prolongation =
240  mg_dof.get_fe().get_prolongation_matrix(
241  child, cell->refinement_case());
242 
243  Assert(prolongation.n() != 0, ExcNoProlongation());
244 
245  cell->child(child)->get_mg_dof_indices(dof_indices_child);
246 
247  replace(this->mg_constrained_dofs,
248  level + 1,
249  dof_indices_child);
250 
251  // now tag the entries in the
252  // matrix which will be used
253  // for this pair of parent/child
254  for (unsigned int i = 0; i < dofs_per_cell; ++i)
255  {
256  entries.resize(0);
257  for (unsigned int j = 0; j < dofs_per_cell; ++j)
258  if (prolongation(i, j) != 0)
259  entries.push_back(dof_indices_parent[j]);
260  dsp.add_entries(dof_indices_child[i],
261  entries.begin(),
262  entries.end());
263  }
264  }
265  }
266 
267 #ifdef DEAL_II_WITH_MPI
268  if (internal::MatrixSelector<
269  VectorType>::requires_distributed_sparsity_pattern)
270  {
271  // Since PETSc matrices do not offer the functionality to fill up in-
272  // complete sparsity patterns on their own, the sparsity pattern must
273  // be manually distributed.
274 
275  // Retrieve communicator from triangulation if it is parallel
276  const parallel::Triangulation<dim, spacedim> *dist_tria =
277  dynamic_cast<const parallel::Triangulation<dim, spacedim> *>(
278  &(mg_dof.get_triangulation()));
279 
280  MPI_Comm communicator = dist_tria != nullptr ?
281  dist_tria->get_communicator() :
282  MPI_COMM_SELF;
283 
284  // Compute # of locally owned MG dofs / processor for distribution
285  const std::vector<::IndexSet>
286  &locally_owned_mg_dofs_per_processor =
287  mg_dof.locally_owned_mg_dofs_per_processor(level + 1);
288  std::vector<::types::global_dof_index>
289  n_locally_owned_mg_dofs_per_processor(
290  locally_owned_mg_dofs_per_processor.size(), 0);
291 
292  for (std::size_t index = 0;
293  index < n_locally_owned_mg_dofs_per_processor.size();
294  ++index)
295  {
296  n_locally_owned_mg_dofs_per_processor[index] =
297  locally_owned_mg_dofs_per_processor[index].n_elements();
298  }
299 
300  // Distribute sparsity pattern
301  ::SparsityTools::distribute_sparsity_pattern(
302  dsp,
303  n_locally_owned_mg_dofs_per_processor,
304  communicator,
305  dsp.row_index_set());
306  }
307 #endif
308 
309  internal::MatrixSelector<VectorType>::reinit(
310  *prolongation_matrices[level],
311  *prolongation_sparsities[level],
312  level,
313  dsp,
314  mg_dof);
315  dsp.reinit(0, 0);
316 
317  // In the end, the entries in this object will only be real valued.
318  // Nevertheless, we have to take the underlying scalar type of the
319  // vector we want to use this class with. The global matrix the entries
320  // of this matrix are copied into has to be able to perform a
321  // matrix-vector multiplication and this is in general only implemented if
322  // the scalar type for matrix and vector is the same. Furthermore,
323  // copying entries between this local object and the global matrix is only
324  // implemented if the objects have the same scalar type.
326 
327  // now actually build the matrices
328  for (cell = mg_dof.begin(level); cell != endc; ++cell)
329  if (cell->has_children() &&
330  (mg_dof.get_triangulation().locally_owned_subdomain() ==
332  cell->level_subdomain_id() ==
333  mg_dof.get_triangulation().locally_owned_subdomain()))
334  {
335  cell->get_mg_dof_indices(dof_indices_parent);
336 
337  replace(this->mg_constrained_dofs, level, dof_indices_parent);
338 
339  Assert(cell->n_children() ==
342  for (unsigned int child = 0; child < cell->n_children(); ++child)
343  {
344  // set an alias to the prolongation matrix for this child
345  prolongation = mg_dof.get_fe().get_prolongation_matrix(
346  child, cell->refinement_case());
347 
348  if (this->mg_constrained_dofs != nullptr &&
349  this->mg_constrained_dofs->have_boundary_indices())
350  for (unsigned int j = 0; j < dofs_per_cell; ++j)
351  if (this->mg_constrained_dofs->is_boundary_index(
352  level, dof_indices_parent[j]))
353  for (unsigned int i = 0; i < dofs_per_cell; ++i)
354  prolongation(i, j) = 0.;
355 
356  cell->child(child)->get_mg_dof_indices(dof_indices_child);
357 
358  replace(this->mg_constrained_dofs,
359  level + 1,
360  dof_indices_child);
361 
362  // now set the entries in the matrix
363  for (unsigned int i = 0; i < dofs_per_cell; ++i)
364  prolongation_matrices[level]->set(dof_indices_child[i],
365  dofs_per_cell,
366  dof_indices_parent.data(),
367  &prolongation(i, 0),
368  true);
369  }
370  }
371  prolongation_matrices[level]->compress(VectorOperation::insert);
372  }
373 
374  this->fill_and_communicate_copy_indices(mg_dof);
375 }
376 
377 
378 
379 template <typename VectorType>
380 void
382 {
383  for (unsigned int level = 0; level < prolongation_matrices.size(); ++level)
384  {
385  os << "Level " << level << std::endl;
386  prolongation_matrices[level]->print(os);
387  os << std::endl;
388  }
389 }
390 
391 
392 
393 template <typename VectorType>
394 std::size_t
396 {
398  for (unsigned int i = 0; i < prolongation_matrices.size(); ++i)
399  result += prolongation_matrices[i]->memory_consumption() +
400  prolongation_sparsities[i]->memory_consumption();
401 
402  return result;
403 }
404 
405 
406 // explicit instantiation
407 #include "mg_transfer_prebuilt.inst"
408 
409 
410 DEAL_II_NAMESPACE_CLOSE
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:258
std::size_t memory_consumption() const
cell_iterator begin(const unsigned int level=0) const
Definition: dof_handler.cc:930
void build_matrices(const DoFHandler< dim, spacedim > &mg_dof)
bool is_boundary_index(const unsigned int level, const types::global_dof_index index) const
cell_iterator end() const
Definition: dof_handler.cc:959
const std::vector< IndexSet > & locally_owned_mg_dofs_per_processor(const unsigned int level) const
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:1179
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:1407
void reinit(const size_type m, const size_type n, const IndexSet &rowset=IndexSet())
virtual MPI_Comm get_communicator() const
Definition: tria_base.cc:155
const IndexSet & row_index_set() const
types::global_dof_index n_dofs() 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
static ::ExceptionBase & ExcNotImplemented()
bool have_boundary_indices() const
static ::ExceptionBase & ExcInternalError()