deal.II version GIT relicensing-2225-g361923b506 2024-12-06 21:30:00+00:00
\(\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// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2003 - 2023 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15
18
21
22#include <deal.II/fe/fe.h>
23
24#include <deal.II/grid/tria.h>
26
37#include <deal.II/lac/vector.h>
38
41
42#include <algorithm>
43
45
46
47template <typename VectorType>
49 const MGConstrainedDoFs &mg_c)
50{
51 this->mg_constrained_dofs = &mg_c;
52}
53
54
55
56template <typename VectorType>
57void
59 const MGConstrainedDoFs &mg_c)
60{
61 this->mg_constrained_dofs = &mg_c;
62}
63
64
65
66template <typename VectorType>
67void
69{
71 prolongation_matrices.resize(0);
72 prolongation_sparsities.resize(0);
73 interface_dofs.resize(0);
74}
75
76
77
78template <typename VectorType>
79void
81 VectorType &dst,
82 const VectorType &src) const
83{
84 Assert((to_level >= 1) && (to_level <= prolongation_matrices.size()),
85 ExcIndexRange(to_level, 1, prolongation_matrices.size() + 1));
86
87 if (this->mg_constrained_dofs != nullptr &&
88 this->mg_constrained_dofs->get_user_constraint_matrix(to_level - 1)
89 .get_local_lines()
90 .size() > 0)
91 {
92 VectorType copy_src(src);
93 this->mg_constrained_dofs->get_user_constraint_matrix(to_level - 1)
94 .distribute(copy_src);
95 prolongation_matrices[to_level - 1]->vmult(dst, copy_src);
96 }
97 else
98 {
99 prolongation_matrices[to_level - 1]->vmult(dst, src);
100 }
101}
102
103
104
105template <typename VectorType>
106void
108 VectorType &dst,
109 const VectorType &src) const
110{
111 Assert((from_level >= 1) && (from_level <= prolongation_matrices.size()),
112 ExcIndexRange(from_level, 1, prolongation_matrices.size() + 1));
113 (void)from_level;
114
115 prolongation_matrices[from_level - 1]->Tvmult_add(dst, src);
116}
117
118
119namespace
120{
125 void
126 replace(const MGConstrainedDoFs *mg_constrained_dofs,
127 const unsigned int level,
128 std::vector<types::global_dof_index> &dof_indices)
129 {
130 if (mg_constrained_dofs != nullptr &&
131 mg_constrained_dofs->get_level_constraints(level).n_constraints() > 0)
132 for (auto &ind : dof_indices)
133 if (mg_constrained_dofs->get_level_constraints(level)
134 .is_identity_constrained(ind))
135 {
136 Assert(mg_constrained_dofs->get_level_constraints(level)
138 ->size() == 1,
140 ind = mg_constrained_dofs->get_level_constraints(level)
142 ->front()
143 .first;
144 }
145 }
146} // namespace
147
148template <typename VectorType>
149template <int dim, int spacedim>
150void
152 const DoFHandler<dim, spacedim> &dof_handler)
153{
154 Assert(dof_handler.has_level_dofs(),
156 "The underlying DoFHandler object has not had its "
157 "distribute_mg_dofs() function called, but this is a prerequisite "
158 "for multigrid transfers. You will need to call this function, "
159 "probably close to where you already call distribute_dofs()."));
160
161 const unsigned int n_levels =
162 dof_handler.get_triangulation().n_global_levels();
163 const unsigned int dofs_per_cell = dof_handler.get_fe().n_dofs_per_cell();
164
165 this->sizes.resize(n_levels);
166 for (unsigned int l = 0; l < n_levels; ++l)
167 this->sizes[l] = dof_handler.n_dofs(l);
168
169 // reset the size of the array of
170 // matrices. call resize(0) first,
171 // in order to delete all elements
172 // and clear their memory. then
173 // repopulate these arrays
174 //
175 // note that on resize(0), the
176 // shared_ptr class takes care of
177 // deleting the object it points to
178 // by itself
179 prolongation_matrices.resize(0);
180 prolongation_sparsities.resize(0);
181 prolongation_matrices.reserve(n_levels - 1);
182 prolongation_sparsities.reserve(n_levels - 1);
183
184 for (unsigned int i = 0; i < n_levels - 1; ++i)
185 {
186 prolongation_sparsities.emplace_back(
188 prolongation_matrices.emplace_back(
190 }
191
192 // two fields which will store the
193 // indices of the multigrid dofs
194 // for a cell and one of its children
195 std::vector<types::global_dof_index> dof_indices_parent(dofs_per_cell);
196 std::vector<types::global_dof_index> dof_indices_child(dofs_per_cell);
197 std::vector<types::global_dof_index> entries(dofs_per_cell);
198
199 // for each level: first build the sparsity
200 // pattern of the matrices and then build the
201 // matrices themselves. note that we only
202 // need to take care of cells on the coarser
203 // level which have children
204 for (unsigned int level = 0; level < n_levels - 1; ++level)
205 {
206 // reset the dimension of the structure. note that for the number of
207 // entries per row, the number of parent dofs coupling to a child dof is
208 // necessary. this, of course, is the number of degrees of freedom per
209 // cell
210 //
211 // increment dofs_per_cell since a useless diagonal element will be
212 // stored
213 const IndexSet level_p1_relevant_dofs =
215 DynamicSparsityPattern dsp(this->sizes[level + 1],
216 this->sizes[level],
217 level_p1_relevant_dofs);
218 for (const auto &cell : dof_handler.cell_iterators_on_level(level))
219 if (cell->has_children() && cell->is_locally_owned_on_level())
220 {
221 cell->get_mg_dof_indices(dof_indices_parent);
222
223 replace(this->mg_constrained_dofs, level, dof_indices_parent);
224
225 Assert(cell->n_children() ==
228 for (unsigned int child = 0; child < cell->n_children(); ++child)
229 {
230 // set an alias to the prolongation matrix for this child
231 const FullMatrix<double> &prolongation =
232 dof_handler.get_fe().get_prolongation_matrix(
233 child, cell->refinement_case());
234
235 Assert(prolongation.n() != 0, ExcNoProlongation());
236
237 cell->child(child)->get_mg_dof_indices(dof_indices_child);
238
239 replace(this->mg_constrained_dofs,
240 level + 1,
241 dof_indices_child);
242
243 // now tag the entries in the
244 // matrix which will be used
245 // for this pair of parent/child
246 for (unsigned int i = 0; i < dofs_per_cell; ++i)
247 {
248 entries.resize(0);
249 for (unsigned int j = 0; j < dofs_per_cell; ++j)
250 if (prolongation(i, j) != 0)
251 entries.push_back(dof_indices_parent[j]);
252 dsp.add_entries(dof_indices_child[i],
253 entries.begin(),
254 entries.end());
255 }
256 }
257 }
258
259#ifdef DEAL_II_WITH_MPI
261 VectorType>::requires_distributed_sparsity_pattern)
262 {
263 // Since PETSc matrices do not offer the functionality to fill up in-
264 // complete sparsity patterns on their own, the sparsity pattern must
265 // be manually distributed.
266
267 // Distribute sparsity pattern
269 dsp,
270 dof_handler.locally_owned_mg_dofs(level + 1),
271 dof_handler.get_mpi_communicator(),
272 dsp.row_index_set());
273 }
274#endif
275
277 *prolongation_matrices[level],
278 *prolongation_sparsities[level],
279 level,
280 dsp,
281 dof_handler);
282 dsp.reinit(0, 0);
283
284 // In the end, the entries in this object will only be real valued.
285 // Nevertheless, we have to take the underlying scalar type of the
286 // vector we want to use this class with. The global matrix the entries
287 // of this matrix are copied into has to be able to perform a
288 // matrix-vector multiplication and this is in general only implemented if
289 // the scalar type for matrix and vector is the same. Furthermore,
290 // copying entries between this local object and the global matrix is only
291 // implemented if the objects have the same scalar type.
293
294 // now actually build the matrices
295 for (const auto &cell : dof_handler.cell_iterators_on_level(level))
296 if (cell->has_children() && cell->is_locally_owned_on_level())
297 {
298 cell->get_mg_dof_indices(dof_indices_parent);
299
300 replace(this->mg_constrained_dofs, level, dof_indices_parent);
301
302 Assert(cell->n_children() ==
305 for (unsigned int child = 0; child < cell->n_children(); ++child)
306 {
307 // set an alias to the prolongation matrix for this child
308 prolongation = dof_handler.get_fe().get_prolongation_matrix(
309 child, cell->refinement_case());
310
311 if (this->mg_constrained_dofs != nullptr &&
312 this->mg_constrained_dofs->have_boundary_indices())
313 for (unsigned int j = 0; j < dofs_per_cell; ++j)
314 if (this->mg_constrained_dofs->is_boundary_index(
315 level, dof_indices_parent[j]))
316 for (unsigned int i = 0; i < dofs_per_cell; ++i)
317 prolongation(i, j) = 0.;
318
319 cell->child(child)->get_mg_dof_indices(dof_indices_child);
320
321 replace(this->mg_constrained_dofs,
322 level + 1,
323 dof_indices_child);
324
325 // now set the entries in the matrix
326 for (unsigned int i = 0; i < dofs_per_cell; ++i)
327 prolongation_matrices[level]->set(dof_indices_child[i],
328 dofs_per_cell,
329 dof_indices_parent.data(),
330 &prolongation(i, 0),
331 true);
332 }
333 }
334 prolongation_matrices[level]->compress(VectorOperation::insert);
335 }
336
337 this->fill_and_communicate_copy_indices(dof_handler);
338}
339
340
341
342template <typename VectorType>
343void
345{
346 for (unsigned int level = 0; level < prolongation_matrices.size(); ++level)
347 {
348 os << "Level " << level << std::endl;
349 prolongation_matrices[level]->print(os);
350 os << std::endl;
351 }
352}
353
354
355
356template <typename VectorType>
357std::size_t
359{
361 for (unsigned int i = 0; i < prolongation_matrices.size(); ++i)
362 result += prolongation_matrices[i]->memory_consumption() +
363 prolongation_sparsities[i]->memory_consumption();
364
365 return result;
366}
367
368
369// explicit instantiation
370#include "mg_transfer_prebuilt.inst"
371
372
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_mpi_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:498
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:499
unsigned int level
Definition grid_out.cc:4632
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:57