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