Reference documentation for deal.II version 9.6.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_block.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
17
20
21#include <deal.II/fe/fe.h>
22
23#include <deal.II/grid/tria.h>
25
29#include <deal.II/lac/vector.h>
30
33#include <deal.II/multigrid/mg_transfer_block.templates.h>
34
35#include <algorithm>
36#include <iostream>
37#include <numeric>
38#include <utility>
39
41
42namespace
43{
51 template <int dim, typename number, int spacedim>
52 void
53 reinit_vector_by_blocks(
54 const DoFHandler<dim, spacedim> &dof_handler,
56 const std::vector<bool> &sel,
57 std::vector<std::vector<types::global_dof_index>> &ndofs)
58 {
59 std::vector<bool> selected = sel;
60 // Compute the number of blocks needed
61 const unsigned int n_selected =
62 std::accumulate(selected.begin(), selected.end(), 0u);
63
64 if (ndofs.empty())
65 {
66 std::vector<std::vector<types::global_dof_index>> new_dofs(
67 dof_handler.get_triangulation().n_levels(),
68 std::vector<types::global_dof_index>(selected.size()));
69 std::swap(ndofs, new_dofs);
70 MGTools::count_dofs_per_block(dof_handler, ndofs);
71 }
72
73 for (unsigned int level = v.min_level(); level <= v.max_level(); ++level)
74 {
75 v[level].reinit(n_selected, 0);
76 unsigned int k = 0;
77 for (unsigned int i = 0;
78 i < selected.size() && (k < v[level].n_blocks());
79 ++i)
80 {
81 if (selected[i])
82 {
83 v[level].block(k++).reinit(ndofs[level][i]);
84 }
85 v[level].collect_sizes();
86 }
87 }
88 }
89
90
97 template <int dim, typename number, int spacedim>
98 void
99 reinit_vector_by_blocks(
100 const DoFHandler<dim, spacedim> &dof_handler,
102 const unsigned int selected_block,
103 std::vector<std::vector<types::global_dof_index>> &ndofs)
104 {
105 const unsigned int n_blocks = dof_handler.get_fe().n_blocks();
106 AssertIndexRange(selected_block, n_blocks);
107
108 std::vector<bool> selected(n_blocks, false);
109 selected[selected_block] = true;
110
111 if (ndofs.empty())
112 {
113 std::vector<std::vector<types::global_dof_index>> new_dofs(
114 dof_handler.get_triangulation().n_levels(),
115 std::vector<types::global_dof_index>(selected.size()));
116 std::swap(ndofs, new_dofs);
117 MGTools::count_dofs_per_block(dof_handler, ndofs);
118 }
119
120 for (unsigned int level = v.min_level(); level <= v.max_level(); ++level)
121 {
122 v[level].reinit(ndofs[level][selected_block]);
123 }
124 }
125} // namespace
126
127
128template <typename number>
129template <int dim, typename number2, int spacedim>
130void
132 const DoFHandler<dim, spacedim> &dof_handler,
134 const BlockVector<number2> &src) const
135{
136 reinit_vector_by_blocks(dof_handler, dst, selected_block, sizes);
137 // For MGTransferBlockSelect, the
138 // multilevel block is always the
139 // first, since only one block is
140 // selected.
141 for (unsigned int level = dof_handler.get_triangulation().n_levels();
142 level != 0;)
143 {
144 --level;
145 for (IT i = copy_indices[selected_block][level].begin();
146 i != copy_indices[selected_block][level].end();
147 ++i)
148 dst[level](i->second) = src.block(selected_block)(i->first);
149 }
150}
151
152
153
154template <typename number>
155template <int dim, typename number2, int spacedim>
156void
158 const DoFHandler<dim, spacedim> &dof_handler,
160 const Vector<number2> &src) const
161{
162 reinit_vector_by_blocks(dof_handler, dst, selected_block, sizes);
163 // For MGTransferBlockSelect, the
164 // multilevel block is always the
165 // first, since only one block is selected.
166 for (unsigned int level = dof_handler.get_triangulation().n_levels();
167 level != 0;)
168 {
169 --level;
170 for (IT i = copy_indices[selected_block][level].begin();
171 i != copy_indices[selected_block][level].end();
172 ++i)
173 dst[level](i->second) = src(i->first);
174 }
175}
176
177
178
179template <typename number>
180template <int dim, typename number2, int spacedim>
181void
183 const DoFHandler<dim, spacedim> &dof_handler,
185 const BlockVector<number2> &src) const
186{
187 reinit_vector_by_blocks(dof_handler, dst, selected, sizes);
188 for (unsigned int level = dof_handler.get_triangulation().n_levels();
189 level != 0;)
190 {
191 --level;
192 for (unsigned int block = 0; block < selected.size(); ++block)
193 if (selected[block])
194 for (IT i = copy_indices[block][level].begin();
195 i != copy_indices[block][level].end();
196 ++i)
197 dst[level].block(mg_block[block])(i->second) =
198 src.block(block)(i->first);
199 }
200}
201
202
203
204template <int dim, int spacedim>
205void
207{
208 const FiniteElement<dim> &fe = dof_handler.get_fe();
209 const unsigned int n_blocks = fe.n_blocks();
210 const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
211 const unsigned int n_levels = dof_handler.get_triangulation().n_levels();
212
213 Assert(selected.size() == n_blocks,
214 ExcDimensionMismatch(selected.size(), n_blocks));
215
216 // Compute the mapping between real
217 // blocks and blocks used for
218 // multigrid computations.
219 mg_block.resize(n_blocks);
220 n_mg_blocks = 0;
221 for (unsigned int i = 0; i < n_blocks; ++i)
222 if (selected[i])
223 mg_block[i] = n_mg_blocks++;
224 else
226
227 // Compute the lengths of all blocks
228 sizes.clear();
229 sizes.resize(n_levels, std::vector<types::global_dof_index>(fe.n_blocks()));
231
232 // Fill some index vectors
233 // for later use.
235 // Compute start indices from sizes
236 for (auto &level_blocks : mg_block_start)
237 {
239 for (types::global_dof_index &level_block_start : level_blocks)
240 {
241 const types::global_dof_index t = level_block_start;
242 level_block_start = k;
243 k += t;
244 }
245 }
246
248 static_cast<const DoFHandler<dim, spacedim> &>(dof_handler));
249
251 for (types::global_dof_index &first_index : block_start)
252 {
253 const types::global_dof_index t = first_index;
254 first_index = k;
255 k += t;
256 }
257 // Build index vectors for
258 // copy_to_mg and
259 // copy_from_mg. These vectors must
260 // be prebuilt, since the
261 // get_dof_indices functions are
262 // too slow
263 copy_indices.resize(n_blocks);
264 for (unsigned int block = 0; block < n_blocks; ++block)
265 if (selected[block])
266 copy_indices[block].resize(n_levels);
267
268 // Building the prolongation matrices starts here!
269
270 // reset the size of the array of
271 // matrices. call resize(0) first,
272 // in order to delete all elements
273 // and clear their memory. then
274 // repopulate these arrays
275 //
276 // note that on resize(0), the
277 // shared_ptr class takes care of
278 // deleting the object it points to
279 // by itself
280 prolongation_matrices.resize(0);
281 prolongation_sparsities.resize(0);
282 prolongation_matrices.reserve(n_levels - 1);
283 prolongation_sparsities.reserve(n_levels - 1);
284
285 for (unsigned int i = 0; i < n_levels - 1; ++i)
286 {
289 }
290
291 // two fields which will store the
292 // indices of the multigrid dofs
293 // for a cell and one of its children
294 std::vector<types::global_dof_index> dof_indices_parent(dofs_per_cell);
295 std::vector<types::global_dof_index> dof_indices_child(dofs_per_cell);
296
297 // for each level: first build the
298 // sparsity pattern of the matrices
299 // and then build the matrices
300 // themselves. note that we only
301 // need to take care of cells on
302 // the coarser level which have
303 // children
304
305 for (unsigned int level = 0; level < n_levels - 1; ++level)
306 {
307 // reset the dimension of the
308 // structure. note that for
309 // the number of entries per
310 // row, the number of parent
311 // dofs coupling to a child dof
312 // is necessary. this, is the
313 // number of degrees of freedom
314 // per cell
315 prolongation_sparsities[level]->reinit(n_blocks, n_blocks);
316 for (unsigned int i = 0; i < n_blocks; ++i)
317 for (unsigned int j = 0; j < n_blocks; ++j)
318 if (i == j)
319 prolongation_sparsities[level]->block(i, j).reinit(
320 sizes[level + 1][i], sizes[level][j], dofs_per_cell + 1);
321 else
322 prolongation_sparsities[level]->block(i, j).reinit(
323 sizes[level + 1][i], sizes[level][j], 0);
324
325 prolongation_sparsities[level]->collect_sizes();
326
328 dof_handler.begin(level);
329 cell != dof_handler.end(level);
330 ++cell)
331 if (cell->has_children())
332 {
333 cell->get_mg_dof_indices(dof_indices_parent);
334
335 Assert(cell->n_children() ==
338 for (unsigned int child = 0; child < cell->n_children(); ++child)
339 {
340 // set an alias to the
341 // prolongation matrix for
342 // this child
343 const FullMatrix<double> &prolongation =
344 dof_handler.get_fe().get_prolongation_matrix(
345 child, cell->refinement_case());
346
347 cell->child(child)->get_mg_dof_indices(dof_indices_child);
348
349 // now tag the entries in the
350 // matrix which will be used
351 // for this pair of parent/child
352 for (unsigned int i = 0; i < dofs_per_cell; ++i)
353 for (unsigned int j = 0; j < dofs_per_cell; ++j)
354 if (prolongation(i, j) != 0)
355 {
356 const unsigned int icomp =
357 fe.system_to_block_index(i).first;
358 const unsigned int jcomp =
359 fe.system_to_block_index(j).first;
360 if ((icomp == jcomp) && selected[icomp])
362 dof_indices_child[i], dof_indices_parent[j]);
363 }
364 }
365 }
366 prolongation_sparsities[level]->compress();
367
369 // now actually build the matrices
371 dof_handler.begin(level);
372 cell != dof_handler.end(level);
373 ++cell)
374 if (cell->has_children())
375 {
376 cell->get_mg_dof_indices(dof_indices_parent);
377
378 Assert(cell->n_children() ==
381 for (unsigned int child = 0; child < cell->n_children(); ++child)
382 {
383 // set an alias to the
384 // prolongation matrix for
385 // this child
386 const FullMatrix<double> &prolongation =
387 dof_handler.get_fe().get_prolongation_matrix(
388 child, cell->refinement_case());
389
390 cell->child(child)->get_mg_dof_indices(dof_indices_child);
391
392 // now set the entries in the
393 // matrix
394 for (unsigned int i = 0; i < dofs_per_cell; ++i)
395 for (unsigned int j = 0; j < dofs_per_cell; ++j)
396 if (prolongation(i, j) != 0)
397 {
398 const unsigned int icomp =
399 fe.system_to_block_index(i).first;
400 const unsigned int jcomp =
401 fe.system_to_block_index(j).first;
402 if ((icomp == jcomp) && selected[icomp])
404 dof_indices_child[i],
405 dof_indices_parent[j],
406 prolongation(i, j));
407 }
408 }
409 }
410 }
411 // impose boundary conditions
412 // but only in the column of
413 // the prolongation matrix
414 if (mg_constrained_dofs != nullptr &&
416 {
417 std::vector<types::global_dof_index> constrain_indices;
418 std::vector<std::vector<bool>> constraints_per_block(n_blocks);
419 for (int level = n_levels - 2; level >= 0; --level)
420 {
422 0)
423 continue;
424
425 // need to delete all the columns in the
426 // matrix that are on the boundary. to achieve
427 // this, create an array as long as there are
428 // matrix columns, and find which columns we
429 // need to filter away.
430 constrain_indices.resize(0);
431 constrain_indices.resize(prolongation_matrices[level]->n(), 0);
435 for (; dof != endd; ++dof)
436 constrain_indices[*dof] = 1;
437
438 unsigned int index = 0;
439 for (unsigned int block = 0; block < n_blocks; ++block)
440 {
441 const types::global_dof_index n_dofs =
442 prolongation_matrices[level]->block(block, block).m();
443 constraints_per_block[block].resize(0);
444 constraints_per_block[block].resize(n_dofs, false);
445 for (types::global_dof_index i = 0; i < n_dofs; ++i, ++index)
446 constraints_per_block[block][i] =
447 (constrain_indices[index] == 1);
448
449 for (types::global_dof_index i = 0; i < n_dofs; ++i)
450 {
452 start_row = prolongation_matrices[level]
453 ->block(block, block)
454 .begin(i),
455 end_row =
456 prolongation_matrices[level]->block(block, block).end(i);
457 for (; start_row != end_row; ++start_row)
458 {
459 if (constraints_per_block[block][start_row->column()])
460 start_row->value() = 0.;
461 }
462 }
463 }
464 }
465 }
466}
467
468template <typename number>
469template <int dim, int spacedim>
470void
472 const DoFHandler<dim, spacedim> &dof_handler,
473 unsigned int select)
474{
475 const FiniteElement<dim> &fe = dof_handler.get_fe();
476 unsigned int n_blocks = dof_handler.get_fe().n_blocks();
477
478 selected_block = select;
479 selected.resize(n_blocks, false);
480 selected[select] = true;
481
482 MGTransferBlockBase::build(dof_handler);
483
484 std::vector<types::global_dof_index> temp_copy_indices;
485 std::vector<types::global_dof_index> global_dof_indices(fe.n_dofs_per_cell());
486 std::vector<types::global_dof_index> level_dof_indices(fe.n_dofs_per_cell());
487
488 for (int level = dof_handler.get_triangulation().n_levels() - 1; level >= 0;
489 --level)
490 {
492 dof_handler.begin_active(level);
493 const typename DoFHandler<dim, spacedim>::active_cell_iterator level_end =
494 dof_handler.end_active(level);
495
496 temp_copy_indices.resize(0);
497 temp_copy_indices.resize(sizes[level][selected_block],
499
500 // Compute coarse level right hand side
501 // by restricting from fine level.
502 for (; level_cell != level_end; ++level_cell)
503 {
504 // get the dof numbers of
505 // this cell for the global
506 // and the level-wise
507 // numbering
508 level_cell->get_dof_indices(global_dof_indices);
509 level_cell->get_mg_dof_indices(level_dof_indices);
510
511 for (unsigned int i = 0; i < fe.n_dofs_per_cell(); ++i)
512 {
513 const unsigned int block = fe.system_to_block_index(i).first;
514 if (selected[block])
515 {
516 if (mg_constrained_dofs != nullptr)
517 {
518 if (!mg_constrained_dofs->at_refinement_edge(
519 level, level_dof_indices[i]))
520 temp_copy_indices[level_dof_indices[i] -
521 mg_block_start[level][block]] =
522 global_dof_indices[i] - block_start[block];
523 }
524 else
525 temp_copy_indices[level_dof_indices[i] -
526 mg_block_start[level][block]] =
527 global_dof_indices[i] - block_start[block];
528 }
529 }
530 }
531
532 // now all the active dofs got a valid entry,
533 // the other ones have an invalid entry. Count
534 // the invalid entries and then resize the
535 // copy_indices object. Then, insert the pairs
536 // of global index and level index into
537 // copy_indices.
538 const types::global_dof_index n_active_dofs =
539 std::count_if(temp_copy_indices.begin(),
540 temp_copy_indices.end(),
541 [](const types::global_dof_index index) {
542 return index != numbers::invalid_dof_index;
543 });
544 copy_indices[selected_block][level].resize(n_active_dofs);
545 types::global_dof_index counter = 0;
546 for (types::global_dof_index i = 0; i < temp_copy_indices.size(); ++i)
547 if (temp_copy_indices[i] != numbers::invalid_dof_index)
548 copy_indices[selected_block][level][counter++] =
549 std::pair<types::global_dof_index, unsigned int>(
550 temp_copy_indices[i], i);
551 Assert(counter == n_active_dofs, ExcInternalError());
552 }
553}
554
555
556template <typename number>
557template <int dim, int spacedim>
558void
560 const std::vector<bool> &sel)
561{
562 const FiniteElement<dim> &fe = dof_handler.get_fe();
563 unsigned int n_blocks = dof_handler.get_fe().n_blocks();
564
565 if (sel.size() != 0)
566 {
567 Assert(sel.size() == n_blocks,
568 ExcDimensionMismatch(sel.size(), n_blocks));
569 selected = sel;
570 }
571 if (selected.empty())
572 selected = std::vector<bool>(n_blocks, true);
573
574 MGTransferBlockBase::build(dof_handler);
575
576 std::vector<std::vector<types::global_dof_index>> temp_copy_indices(n_blocks);
577 std::vector<types::global_dof_index> global_dof_indices(fe.n_dofs_per_cell());
578 std::vector<types::global_dof_index> level_dof_indices(fe.n_dofs_per_cell());
579 for (int level = dof_handler.get_triangulation().n_levels() - 1; level >= 0;
580 --level)
581 {
583 dof_handler.begin_active(level);
584 const typename DoFHandler<dim, spacedim>::active_cell_iterator level_end =
585 dof_handler.end_active(level);
586
587 for (unsigned int block = 0; block < n_blocks; ++block)
588 if (selected[block])
589 {
590 temp_copy_indices[block].resize(0);
591 temp_copy_indices[block].resize(sizes[level][block],
593 }
594
595 // Compute coarse level right hand side
596 // by restricting from fine level.
597 for (; level_cell != level_end; ++level_cell)
598 {
599 // get the dof numbers of
600 // this cell for the global
601 // and the level-wise
602 // numbering
603 level_cell->get_dof_indices(global_dof_indices);
604 level_cell->get_mg_dof_indices(level_dof_indices);
605
606 for (unsigned int i = 0; i < fe.n_dofs_per_cell(); ++i)
607 {
608 const unsigned int block = fe.system_to_block_index(i).first;
609 if (selected[block])
610 temp_copy_indices[block][level_dof_indices[i] -
611 mg_block_start[level][block]] =
612 global_dof_indices[i] - block_start[block];
613 }
614 }
615
616 for (unsigned int block = 0; block < n_blocks; ++block)
617 if (selected[block])
618 {
619 const types::global_dof_index n_active_dofs =
620 std::count_if(temp_copy_indices[block].begin(),
621 temp_copy_indices[block].end(),
622 [](const types::global_dof_index index) {
623 return index != numbers::invalid_dof_index;
624 });
625 copy_indices[block][level].resize(n_active_dofs);
626 types::global_dof_index counter = 0;
627 for (types::global_dof_index i = 0;
628 i < temp_copy_indices[block].size();
629 ++i)
630 if (temp_copy_indices[block][i] != numbers::invalid_dof_index)
631 copy_indices[block][level][counter++] =
632 std::pair<types::global_dof_index, unsigned int>(
633 temp_copy_indices[block][i], i);
634 Assert(counter == n_active_dofs, ExcInternalError());
635 }
636 }
637}
638
639
640
641// explicit instantiations
642#include "mg_transfer_block.inst"
643
644
BlockType & block(const unsigned int i)
cell_iterator end() const
const FiniteElement< dim, spacedim > & get_fe(const types::fe_index index=0) const
const Triangulation< dim, spacedim > & get_triangulation() const
active_cell_iterator begin_active(const unsigned int level=0) const
active_cell_iterator end_active(const unsigned int level) const
cell_iterator begin(const unsigned int level=0) const
unsigned int n_dofs_per_cell() const
unsigned int n_blocks() const
std::pair< unsigned int, types::global_dof_index > system_to_block_index(const unsigned int component) 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_elements() const
Definition index_set.h:1934
ElementIterator begin() const
Definition index_set.h:1710
ElementIterator end() const
Definition index_set.h:1722
bool have_boundary_indices() const
const IndexSet & get_boundary_indices(const unsigned int level) const
std::vector< unsigned int > mg_block
std::vector< std::vector< types::global_dof_index > > sizes
void build(const DoFHandler< dim, spacedim > &dof_handler)
std::vector< bool > selected
std::vector< types::global_dof_index > block_start
std::vector< std::shared_ptr< BlockSparseMatrix< double > > > prolongation_matrices
std::vector< std::shared_ptr< BlockSparsityPattern > > prolongation_sparsities
SmartPointer< const MGConstrainedDoFs, MGTransferBlockBase > mg_constrained_dofs
std::vector< std::vector< types::global_dof_index > > mg_block_start
std::vector< std::vector< std::vector< std::pair< unsigned int, unsigned int > > > > copy_indices
void copy_to_mg(const DoFHandler< dim, spacedim > &dof_handler, MGLevelObject< Vector< number > > &dst, const Vector< number2 > &src) const
void build(const DoFHandler< dim, spacedim > &dof_handler, unsigned int selected)
void build(const DoFHandler< dim, spacedim > &dof_handler, const std::vector< bool > &selected)
void copy_to_mg(const DoFHandler< dim, spacedim > &dof_handler, MGLevelObject< BlockVector< number > > &dst, const BlockVector< number2 > &src) const
unsigned int n_levels() const
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:503
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:504
unsigned int level
Definition grid_out.cc:4626
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
typename ActiveSelector::cell_iterator cell_iterator
typename ActiveSelector::active_cell_iterator active_cell_iterator
std::vector< types::global_dof_index > count_dofs_per_fe_block(const DoFHandler< dim, spacedim > &dof, const std::vector< unsigned int > &target_block=std::vector< unsigned int >())
void count_dofs_per_block(const DoFHandler< dim, spacedim > &dof_handler, std::vector< std::vector< types::global_dof_index > > &dofs_per_block, std::vector< unsigned int > target_block={})
Definition mg_tools.cc:1156
std::enable_if_t< IsBlockVector< VectorType >::value, unsigned int > n_blocks(const VectorType &vector)
Definition operators.h:49
static const unsigned int invalid_unsigned_int
Definition types.h:220
const types::global_dof_index invalid_dof_index
Definition types.h:252