Reference documentation for deal.II version 8.4.1
mg_transfer_prebuilt.cc
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2003 - 2016 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 at
12 // the top level of the deal.II distribution.
13 //
14 // ---------------------------------------------------------------------
15 
16 
17 #include <deal.II/base/logstream.h>
18 #include <deal.II/base/function.h>
19 
20 #include <deal.II/lac/vector.h>
21 #include <deal.II/lac/block_vector.h>
22 #include <deal.II/lac/parallel_vector.h>
23 #include <deal.II/lac/parallel_block_vector.h>
24 #include <deal.II/lac/petsc_vector.h>
25 #include <deal.II/lac/petsc_block_vector.h>
26 #include <deal.II/lac/trilinos_vector.h>
27 #include <deal.II/lac/trilinos_block_vector.h>
28 #include <deal.II/lac/sparse_matrix.h>
29 #include <deal.II/lac/block_sparse_matrix.h>
30 #include <deal.II/lac/dynamic_sparsity_pattern.h>
31 #include <deal.II/grid/tria.h>
32 #include <deal.II/grid/tria_iterator.h>
33 #include <deal.II/dofs/dof_tools.h>
34 #include <deal.II/fe/fe.h>
35 #include <deal.II/dofs/dof_accessor.h>
36 #include <deal.II/multigrid/mg_tools.h>
37 #include <deal.II/multigrid/mg_transfer.h>
38 
39 #include <algorithm>
40 
41 DEAL_II_NAMESPACE_OPEN
42 
43 
44 template<typename VectorType>
46 {}
47 
48 
49 
50 template<typename VectorType>
52  :
53  constraints(&c)
54 {
55  this->mg_constrained_dofs = &mg_c;
56 }
57 
58 
59 
60 template <typename VectorType>
62 {}
63 
64 
65 
66 template <typename VectorType>
68 (const ConstraintMatrix &c, const MGConstrainedDoFs &mg_c)
69 {
70  constraints = &c;
71  this->mg_constrained_dofs = &mg_c;
72 }
73 
74 
75 
76 template <typename VectorType>
78 {
80  prolongation_matrices.resize(0);
81  prolongation_sparsities.resize(0);
82  interface_dofs.resize(0);
83  constraints = 0;
84 }
85 
86 
87 
88 template <typename VectorType>
89 void MGTransferPrebuilt<VectorType>::prolongate (const unsigned int to_level,
90  VectorType &dst,
91  const VectorType &src) const
92 {
93  Assert ((to_level >= 1) && (to_level<=prolongation_matrices.size()),
94  ExcIndexRange (to_level, 1, prolongation_matrices.size()+1));
95 
96  prolongation_matrices[to_level-1]->vmult (dst, src);
97 }
98 
99 
100 
101 template <typename VectorType>
102 void MGTransferPrebuilt<VectorType>::restrict_and_add (const unsigned int from_level,
103  VectorType &dst,
104  const VectorType &src) const
105 {
106  Assert ((from_level >= 1) && (from_level<=prolongation_matrices.size()),
107  ExcIndexRange (from_level, 1, prolongation_matrices.size()+1));
108  (void)from_level;
109 
110  prolongation_matrices[from_level-1]->Tvmult_add (dst, src);
111 }
112 
113 
114 
115 template <typename VectorType>
116 template <int dim, int spacedim>
119 {
120  const unsigned int n_levels = mg_dof.get_triangulation().n_global_levels();
121  const unsigned int dofs_per_cell = mg_dof.get_fe().dofs_per_cell;
122 
123  this->sizes.resize(n_levels);
124  for (unsigned int l=0; l<n_levels; ++l)
125  this->sizes[l] = mg_dof.n_dofs(l);
126 
127  // reset the size of the array of
128  // matrices. call resize(0) first,
129  // in order to delete all elements
130  // and clear their memory. then
131  // repopulate these arrays
132  //
133  // note that on resize(0), the
134  // shared_ptr class takes care of
135  // deleting the object it points to
136  // by itself
137  prolongation_matrices.resize (0);
138  prolongation_sparsities.resize (0);
139 
140  for (unsigned int i=0; i<n_levels-1; ++i)
141  {
142  prolongation_sparsities.push_back
144  prolongation_matrices.push_back
146  }
147 
148  // two fields which will store the
149  // indices of the multigrid dofs
150  // for a cell and one of its children
151  std::vector<types::global_dof_index> dof_indices_parent (dofs_per_cell);
152  std::vector<types::global_dof_index> dof_indices_child (dofs_per_cell);
153  std::vector<types::global_dof_index> entries (dofs_per_cell);
154 
155  // for each level: first build the sparsity
156  // pattern of the matrices and then build the
157  // matrices themselves. note that we only
158  // need to take care of cells on the coarser
159  // level which have children
160  for (unsigned int level=0; level<n_levels-1; ++level)
161  {
162  // reset the dimension of the structure. note that for the number of
163  // entries per row, the number of parent dofs coupling to a child dof is
164  // necessary. this, of course, is the number of degrees of freedom per
165  // cell
166  //
167  // increment dofs_per_cell since a useless diagonal element will be
168  // stored
169  IndexSet level_p1_relevant_dofs;
171  level_p1_relevant_dofs);
172  DynamicSparsityPattern dsp (this->sizes[level+1],
173  this->sizes[level],
174  level_p1_relevant_dofs);
175  for (typename DoFHandler<dim,spacedim>::cell_iterator cell=mg_dof.begin(level);
176  cell != mg_dof.end(level); ++cell)
177  if (cell->has_children() &&
178  ( mg_dof.get_triangulation().locally_owned_subdomain()==numbers::invalid_subdomain_id
179  || cell->level_subdomain_id()==mg_dof.get_triangulation().locally_owned_subdomain()
180  ))
181  {
182  cell->get_mg_dof_indices (dof_indices_parent);
183 
185  ExcNotImplemented());
186  for (unsigned int child=0; child<cell->n_children(); ++child)
187  {
188  // set an alias to the prolongation matrix for this child
189  const FullMatrix<double> &prolongation
190  = mg_dof.get_fe().get_prolongation_matrix (child,
191  cell->refinement_case());
192 
193  Assert (prolongation.n() != 0, ExcNoProlongation());
194 
195  cell->child(child)->get_mg_dof_indices (dof_indices_child);
196 
197  // now tag the entries in the
198  // matrix which will be used
199  // for this pair of parent/child
200  for (unsigned int i=0; i<dofs_per_cell; ++i)
201  {
202  entries.resize(0);
203  for (unsigned int j=0; j<dofs_per_cell; ++j)
204  if (prolongation(i,j) != 0)
205  entries.push_back (dof_indices_parent[j]);
206  dsp.add_entries (dof_indices_child[i],
207  entries.begin(), entries.end());
208  }
209  }
210  }
211 
212  internal::MatrixSelector<VectorType>::reinit(*prolongation_matrices[level],
213  *prolongation_sparsities[level],
214  level,
215  dsp,
216  mg_dof);
217  dsp.reinit(0,0);
218 
219  FullMatrix<double> prolongation;
220 
221  // now actually build the matrices
222  for (typename DoFHandler<dim,spacedim>::cell_iterator cell=mg_dof.begin(level);
223  cell != mg_dof.end(level); ++cell)
224  if (cell->has_children() &&
225  (mg_dof.get_triangulation().locally_owned_subdomain()==numbers::invalid_subdomain_id
226  || cell->level_subdomain_id()==mg_dof.get_triangulation().locally_owned_subdomain())
227  )
228  {
229  cell->get_mg_dof_indices (dof_indices_parent);
230 
232  ExcNotImplemented());
233  for (unsigned int child=0; child<cell->n_children(); ++child)
234  {
235  // set an alias to the prolongation matrix for this child
236  prolongation
237  = mg_dof.get_fe().get_prolongation_matrix (child,
238  cell->refinement_case());
239 
240  if (this->mg_constrained_dofs != 0 &&
241  this->mg_constrained_dofs->have_boundary_indices())
242  for (unsigned int j=0; j<dofs_per_cell; ++j)
243  if (this->mg_constrained_dofs->is_boundary_index(level, dof_indices_parent[j]))
244  for (unsigned int i=0; i<dofs_per_cell; ++i)
245  prolongation(i,j) = 0.;
246 
247  cell->child(child)->get_mg_dof_indices (dof_indices_child);
248 
249  // now set the entries in the matrix
250  for (unsigned int i=0; i<dofs_per_cell; ++i)
251  prolongation_matrices[level]->set (dof_indices_child[i],
252  dofs_per_cell,
253  &dof_indices_parent[0],
254  &prolongation(i,0),
255  true);
256  }
257  }
258  prolongation_matrices[level]->compress(VectorOperation::insert);
259  }
260 
261  this->fill_and_communicate_copy_indices(mg_dof);
262 }
263 
264 
265 
266 template <typename VectorType>
267 void
269 {
270  for (unsigned int level = 0; level<prolongation_matrices.size(); ++level)
271  {
272  os << "Level " << level << std::endl;
273  prolongation_matrices[level]->print(os);
274  os << std::endl;
275  }
276 }
277 
278 
279 
280 template <typename VectorType>
281 std::size_t
283 {
285  for (unsigned int i=0; i<prolongation_matrices.size(); ++i)
286  result += prolongation_matrices[i]->memory_consumption()
287  + prolongation_sparsities[i]->memory_consumption();
288 
289  return result;
290 }
291 
292 
293 // explicit instantiation
294 #include "mg_transfer_prebuilt.inst"
295 
296 
297 DEAL_II_NAMESPACE_CLOSE
std::size_t memory_consumption() const
void add_entries(const size_type row, ForwardIterator begin, ForwardIterator end, const bool indices_are_unique_and_sorted=false)
virtual void restrict_and_add(const unsigned int from_level, VectorType &dst, const VectorType &src) const
const types::subdomain_id invalid_subdomain_id
Definition: types.h:239
void build_matrices(const DoFHandler< dim, spacedim > &mg_dof)
const Triangulation< dim, spacedim > & get_triangulation() const
void extract_locally_relevant_level_dofs(const DoFHandlerType &dof_handler, const unsigned int level, IndexSet &dof_set)
Definition: dof_tools.cc:983
void initialize_constraints(const ConstraintMatrix &constraints, const MGConstrainedDoFs &mg_constrained_dofs)
size_type n() const
void print_matrices(std::ostream &os) const
ActiveSelector::cell_iterator cell_iterator
Definition: dof_handler.h:249
#define Assert(cond, exc)
Definition: exceptions.h:294
types::global_dof_index n_dofs() const
void reinit(const size_type m, const size_type n, const IndexSet &rowset=IndexSet())
SmartPointer< const MGConstrainedDoFs, MGLevelGlobalTransfer< VectorType > > mg_constrained_dofs
Definition: mg_transfer.h:253
cell_iterator end() const
Definition: dof_handler.cc:861
virtual void prolongate(const unsigned int to_level, VectorType &dst, const VectorType &src) const
cell_iterator begin(const unsigned int level=0) const
Definition: dof_handler.cc:833
std::size_t memory_consumption() const
const FiniteElement< dim, spacedim > & get_fe() const