Reference documentation for deal.II version 9.5.0
\(\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\}}\)
Loading...
Searching...
No Matches
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
19
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
48template <typename VectorType>
50 const MGConstrainedDoFs &mg_c)
51{
52 this->mg_constrained_dofs = &mg_c;
53}
54
55
56
57template <typename VectorType>
58void
60 const MGConstrainedDoFs &mg_c)
61{
62 this->mg_constrained_dofs = &mg_c;
63}
64
65
66
67template <typename VectorType>
68void
70{
72 prolongation_matrices.resize(0);
73 prolongation_sparsities.resize(0);
74 interface_dofs.resize(0);
75}
76
77
78
79template <typename VectorType>
80void
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
106template <typename VectorType>
107void
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
120namespace
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)
135 .is_identity_constrained(ind))
136 {
137 Assert(mg_constrained_dofs->get_level_constraints(level)
139 ->size() == 1,
141 ind = mg_constrained_dofs->get_level_constraints(level)
143 ->front()
144 .first;
145 }
146 }
147} // namespace
148
149template <typename VectorType>
150template <int dim, int spacedim>
151void
153 const DoFHandler<dim, spacedim> &dof_handler)
154{
155 Assert(dof_handler.has_level_dofs(),
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 for (const auto &cell : dof_handler.cell_iterators_on_level(level))
222 if (cell->has_children() && cell->is_locally_owned_on_level())
223 {
224 cell->get_mg_dof_indices(dof_indices_parent);
225
226 replace(this->mg_constrained_dofs, level, dof_indices_parent);
227
228 Assert(cell->n_children() ==
231 for (unsigned int child = 0; child < cell->n_children(); ++child)
232 {
233 // set an alias to the prolongation matrix for this child
234 const FullMatrix<double> &prolongation =
235 dof_handler.get_fe().get_prolongation_matrix(
236 child, cell->refinement_case());
237
238 Assert(prolongation.n() != 0, ExcNoProlongation());
239
240 cell->child(child)->get_mg_dof_indices(dof_indices_child);
241
242 replace(this->mg_constrained_dofs,
243 level + 1,
244 dof_indices_child);
245
246 // now tag the entries in the
247 // matrix which will be used
248 // for this pair of parent/child
249 for (unsigned int i = 0; i < dofs_per_cell; ++i)
250 {
251 entries.resize(0);
252 for (unsigned int j = 0; j < dofs_per_cell; ++j)
253 if (prolongation(i, j) != 0)
254 entries.push_back(dof_indices_parent[j]);
255 dsp.add_entries(dof_indices_child[i],
256 entries.begin(),
257 entries.end());
258 }
259 }
260 }
261
262#ifdef DEAL_II_WITH_MPI
264 VectorType>::requires_distributed_sparsity_pattern)
265 {
266 // Since PETSc matrices do not offer the functionality to fill up in-
267 // complete sparsity patterns on their own, the sparsity pattern must
268 // be manually distributed.
269
270 // Distribute sparsity pattern
272 dsp,
273 dof_handler.locally_owned_mg_dofs(level + 1),
274 dof_handler.get_communicator(),
275 dsp.row_index_set());
276 }
277#endif
278
280 *prolongation_matrices[level],
281 *prolongation_sparsities[level],
282 level,
283 dsp,
284 dof_handler);
285 dsp.reinit(0, 0);
286
287 // In the end, the entries in this object will only be real valued.
288 // Nevertheless, we have to take the underlying scalar type of the
289 // vector we want to use this class with. The global matrix the entries
290 // of this matrix are copied into has to be able to perform a
291 // matrix-vector multiplication and this is in general only implemented if
292 // the scalar type for matrix and vector is the same. Furthermore,
293 // copying entries between this local object and the global matrix is only
294 // implemented if the objects have the same scalar type.
296
297 // now actually build the matrices
298 for (const auto &cell : dof_handler.cell_iterators_on_level(level))
299 if (cell->has_children() && cell->is_locally_owned_on_level())
300 {
301 cell->get_mg_dof_indices(dof_indices_parent);
302
303 replace(this->mg_constrained_dofs, level, dof_indices_parent);
304
305 Assert(cell->n_children() ==
308 for (unsigned int child = 0; child < cell->n_children(); ++child)
309 {
310 // set an alias to the prolongation matrix for this child
311 prolongation = dof_handler.get_fe().get_prolongation_matrix(
312 child, cell->refinement_case());
313
314 if (this->mg_constrained_dofs != nullptr &&
315 this->mg_constrained_dofs->have_boundary_indices())
316 for (unsigned int j = 0; j < dofs_per_cell; ++j)
317 if (this->mg_constrained_dofs->is_boundary_index(
318 level, dof_indices_parent[j]))
319 for (unsigned int i = 0; i < dofs_per_cell; ++i)
320 prolongation(i, j) = 0.;
321
322 cell->child(child)->get_mg_dof_indices(dof_indices_child);
323
324 replace(this->mg_constrained_dofs,
325 level + 1,
326 dof_indices_child);
327
328 // now set the entries in the matrix
329 for (unsigned int i = 0; i < dofs_per_cell; ++i)
330 prolongation_matrices[level]->set(dof_indices_child[i],
331 dofs_per_cell,
332 dof_indices_parent.data(),
333 &prolongation(i, 0),
334 true);
335 }
336 }
337 prolongation_matrices[level]->compress(VectorOperation::insert);
338 }
339
340 this->fill_and_communicate_copy_indices(dof_handler);
341}
342
343
344
345template <typename VectorType>
346void
348{
349 for (unsigned int level = 0; level < prolongation_matrices.size(); ++level)
350 {
351 os << "Level " << level << std::endl;
352 prolongation_matrices[level]->print(os);
353 os << std::endl;
354 }
355}
356
357
358
359template <typename VectorType>
360std::size_t
362{
364 for (unsigned int i = 0; i < prolongation_matrices.size(); ++i)
365 result += prolongation_matrices[i]->memory_consumption() +
366 prolongation_sparsities[i]->memory_consumption();
367
368 return result;
369}
370
371
372// explicit instantiation
373#include "mg_transfer_prebuilt.inst"
374
375
const std::vector< std::pair< size_type, number > > * get_constraint_entries(const size_type line_n) const
size_type n_constraints() const
const FiniteElement< dim, spacedim > & get_fe(const types::fe_index index=0) const
const IndexSet & locally_owned_mg_dofs(const unsigned int level) const
const Triangulation< dim, spacedim > & get_triangulation() const
types::global_dof_index n_dofs() const
MPI_Comm get_communicator() const
bool has_level_dofs() const
void add_entries(const size_type row, ForwardIterator begin, ForwardIterator end, const bool indices_are_unique_and_sorted=false)
const IndexSet & row_index_set() const
void reinit(const size_type m, const size_type n, const IndexSet &rowset=IndexSet())
unsigned int n_dofs_per_cell() const
virtual const FullMatrix< double > & get_prolongation_matrix(const unsigned int child, const RefinementCase< dim > &refinement_case=RefinementCase< dim >::isotropic_refinement) const
size_type n() const
bool is_boundary_index(const unsigned int level, const types::global_dof_index index) const
bool have_boundary_indices() const
const AffineConstraints< double > & get_level_constraints(const unsigned int level) const
std::size_t memory_consumption() const
virtual void restrict_and_add(const unsigned int from_level, VectorType &dst, const VectorType &src) const override
void build(const DoFHandler< dim, spacedim > &dof_handler)
std::size_t memory_consumption() const
MGTransferPrebuilt()=default
virtual void prolongate(const unsigned int to_level, VectorType &dst, const VectorType &src) const override
void print_matrices(std::ostream &os) const
void initialize_constraints(const MGConstrainedDoFs &mg_constrained_dofs)
virtual unsigned int n_global_levels() const
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
unsigned int level
Definition grid_out.cc:4618
IteratorRange< cell_iterator > cell_iterators_on_level(const unsigned int level) const
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcIndexRange(std::size_t arg1, std::size_t arg2, std::size_t arg3)
static ::ExceptionBase & ExcMessage(std::string arg1)
IndexSet extract_locally_relevant_level_dofs(const DoFHandler< dim, spacedim > &dof_handler, const unsigned int level)
void distribute_sparsity_pattern(DynamicSparsityPattern &dsp, const IndexSet &locally_owned_rows, const MPI_Comm mpi_comm, const IndexSet &locally_relevant_rows)
static void reinit(Matrix &matrix, Sparsity &sparsity, int level, const SparsityPatternType &sp, const DoFHandler< dim, spacedim > &)
Definition mg_transfer.h:58