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_component.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
19
22
23#include <deal.II/fe/fe.h>
24
25#include <deal.II/grid/tria.h>
27
32#include <deal.II/lac/vector.h>
33
36#include <deal.II/multigrid/mg_transfer_component.templates.h>
37
38#include <algorithm>
39#include <iostream>
40#include <numeric>
41
43
44
45namespace
46{
79 template <int dim, typename number, int spacedim>
80 void
81 reinit_vector_by_components(
82 const ::DoFHandler<dim, spacedim> & mg_dof,
84 const std::vector<bool> & sel,
85 const std::vector<unsigned int> & target_comp,
86 std::vector<std::vector<types::global_dof_index>> &ndofs)
87 {
88 std::vector<bool> selected = sel;
89 std::vector<unsigned int> target_component = target_comp;
90 const unsigned int ncomp = mg_dof.get_fe(0).n_components();
91
92 // If the selected and
93 // target_component have size 0,
94 // they must be replaced by default
95 // values.
96 //
97 // Since we already made copies
98 // directly after this function was
99 // called, we use the arguments
100 // directly.
101 if (target_component.size() == 0)
102 {
103 target_component.resize(ncomp);
104 for (unsigned int i = 0; i < ncomp; ++i)
105 target_component[i] = i;
106 }
107
108 // If selected is an empty vector,
109 // all components are selected.
110 if (selected.size() == 0)
111 {
112 selected.resize(target_component.size());
113 std::fill_n(selected.begin(), ncomp, false);
114 for (const unsigned int component : target_component)
115 selected[component] = true;
116 }
117
118 Assert(selected.size() == target_component.size(),
119 ExcDimensionMismatch(selected.size(), target_component.size()));
120
121 // Compute the number of blocks needed
122 const unsigned int n_selected =
123 std::accumulate(selected.begin(), selected.end(), 0u);
124
125 if (ndofs.size() == 0)
126 {
127 std::vector<std::vector<types::global_dof_index>> new_dofs(
128 mg_dof.get_triangulation().n_levels(),
129 std::vector<types::global_dof_index>(target_component.size()));
130 std::swap(ndofs, new_dofs);
131 MGTools::count_dofs_per_block(mg_dof, ndofs, target_component);
132 }
133
134 for (unsigned int level = v.min_level(); level <= v.max_level(); ++level)
135 {
136 v[level].reinit(n_selected, 0);
137 unsigned int k = 0;
138 for (unsigned int i = 0;
139 i < selected.size() && (k < v[level].n_blocks());
140 ++i)
141 {
142 if (selected[i])
143 {
144 v[level].block(k++).reinit(ndofs[level][i]);
145 }
146 v[level].collect_sizes();
147 }
148 }
149 }
150
151
178 template <int dim, typename number, int spacedim>
179 void
180 reinit_vector_by_components(
181 const ::DoFHandler<dim, spacedim> & mg_dof,
183 const ComponentMask & component_mask,
184 const std::vector<unsigned int> & target_component,
185 std::vector<std::vector<types::global_dof_index>> &ndofs)
186 {
187 Assert(component_mask.represents_n_components(target_component.size()),
188 ExcMessage("The component mask does not have the correct size."));
189
190 unsigned int selected_block = 0;
191 for (unsigned int i = 0; i < target_component.size(); ++i)
192 if (component_mask[i])
193 selected_block = target_component[i];
194
195 if (ndofs.size() == 0)
196 {
197 std::vector<std::vector<types::global_dof_index>> new_dofs(
198 mg_dof.get_triangulation().n_levels(),
199 std::vector<types::global_dof_index>(target_component.size()));
200 std::swap(ndofs, new_dofs);
201 MGTools::count_dofs_per_block(mg_dof, ndofs, target_component);
202 }
203
204 for (unsigned int level = v.min_level(); level <= v.max_level(); ++level)
205 {
206 v[level].reinit(ndofs[level][selected_block]);
207 }
208 }
209} // namespace
210
211
212template <typename number>
213template <int dim, class InVector, int spacedim>
214void
216 const DoFHandler<dim, spacedim> &mg_dof_handler,
218 const InVector & src) const
219{
220 dst = 0;
221
222 Assert(sizes.size() == mg_dof_handler.get_triangulation().n_levels(),
223 ExcMatricesNotBuilt());
224
225 reinit_vector_by_components(
226 mg_dof_handler, dst, mg_component_mask, mg_target_component, sizes);
227
228 // traverse the grid top-down
229 // (i.e. starting with the most
230 // refined grid). this way, we can
231 // always get that part of one
232 // level of the output vector which
233 // corresponds to a region which is
234 // more refined, by restriction of
235 // the respective vector on the
236 // next finer level, which we then
237 // already have built.
238
239 for (unsigned int level = mg_dof_handler.get_triangulation().n_levels();
240 level != 0;)
241 {
242 --level;
243
244 using IT = std::vector<
245 std::pair<types::global_dof_index, unsigned int>>::const_iterator;
246 for (IT i = copy_to_and_from_indices[level].begin();
247 i != copy_to_and_from_indices[level].end();
248 ++i)
249 dst[level](i->second) = src(i->first);
250 }
251}
252
253
254
255template <int dim, int spacedim>
256void
258 const DoFHandler<dim, spacedim> &mg_dof)
259{
260 build(mg_dof);
261}
262
263
264
265template <int dim, int spacedim>
266void
268{
269 // Fill target component with
270 // standard values (identity) if it
271 // is empty
272 if (target_component.size() == 0)
273 {
274 target_component.resize(mg_dof.get_fe(0).n_components());
275 for (unsigned int i = 0; i < target_component.size(); ++i)
276 target_component[i] = i;
277 }
278 else
279 {
280 // otherwise, check it for consistency
281 Assert(target_component.size() == mg_dof.get_fe(0).n_components(),
283 mg_dof.get_fe(0).n_components()));
284
285 for (unsigned int i = 0; i < target_component.size(); ++i)
286 {
288 }
289 }
290 // Do the same for the multilevel
291 // components. These may be
292 // different.
293 if (mg_target_component.size() == 0)
294 {
295 mg_target_component.resize(mg_dof.get_fe(0).n_components());
296 for (unsigned int i = 0; i < mg_target_component.size(); ++i)
298 }
299 else
300 {
301 Assert(mg_target_component.size() == mg_dof.get_fe(0).n_components(),
303 mg_dof.get_fe(0).n_components()));
304
305 for (unsigned int i = 0; i < mg_target_component.size(); ++i)
306 {
308 }
309 }
310
311 const FiniteElement<dim> &fe = mg_dof.get_fe();
312
313 // Effective number of components
314 // is the maximum entry in
315 // mg_target_component. This
316 // assumes that the values in that
317 // vector don't have holes.
318 const unsigned int n_components =
319 *std::max_element(mg_target_component.begin(), mg_target_component.end()) +
320 1;
321 const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
322 const unsigned int n_levels = mg_dof.get_triangulation().n_levels();
323
325 ExcMessage("Component mask has wrong size."));
326
327 // Compute the lengths of all blocks
328 sizes.resize(n_levels);
329 for (unsigned int l = 0; l < n_levels; ++l)
330 sizes[l].resize(n_components);
331
333
334 // Fill some index vectors
335 // for later use.
337 // Compute start indices from sizes
338 for (auto &level_components : mg_component_start)
339 {
341 for (types::global_dof_index &level_component_start : level_components)
342 {
343 const types::global_dof_index t = level_component_start;
344 level_component_start = k;
345 k += t;
346 }
347 }
348
350
352 for (types::global_dof_index &first_index : component_start)
353 {
354 const types::global_dof_index t = first_index;
355 first_index = k;
356 k += t;
357 }
358
359 // Build index vectors for
360 // copy_to_mg and
361 // copy_from_mg. These vectors must
362 // be prebuilt, since the
363 // get_dof_indices functions are
364 // too slow
365
366 copy_to_and_from_indices.resize(n_levels);
367
368 // Building the prolongation matrices starts here!
369
370 // reset the size of the array of
371 // matrices. call resize(0) first,
372 // in order to delete all elements
373 // and clear their memory. then
374 // repopulate these arrays
375 //
376 // note that on resize(0), the
377 // shared_ptr class takes care of
378 // deleting the object it points to
379 // by itself
380 prolongation_matrices.resize(0);
381 prolongation_sparsities.resize(0);
382 prolongation_matrices.reserve(n_levels - 1);
383 prolongation_sparsities.reserve(n_levels - 1);
384
385 for (unsigned int i = 0; i < n_levels - 1; ++i)
386 {
389 }
390
391 // two fields which will store the
392 // indices of the multigrid dofs
393 // for a cell and one of its children
394 std::vector<types::global_dof_index> dof_indices_parent(dofs_per_cell);
395 std::vector<types::global_dof_index> dof_indices_child(dofs_per_cell);
396
397 // for each level: first build the
398 // sparsity pattern of the matrices
399 // and then build the matrices
400 // themselves. note that we only
401 // need to take care of cells on
402 // the coarser level which have
403 // children
404 for (unsigned int level = 0; level < n_levels - 1; ++level)
405 {
406 // reset the dimension of the
407 // structure. note that for
408 // the number of entries per
409 // row, the number of parent
410 // dofs coupling to a child dof
411 // is necessary. this, is the
412 // number of degrees of freedom
413 // per cell
414 prolongation_sparsities[level]->reinit(n_components, n_components);
415 for (unsigned int i = 0; i < n_components; ++i)
416 for (unsigned int j = 0; j < n_components; ++j)
417 if (i == j)
418 prolongation_sparsities[level]->block(i, j).reinit(
419 sizes[level + 1][i], sizes[level][j], dofs_per_cell + 1);
420 else
421 prolongation_sparsities[level]->block(i, j).reinit(
422 sizes[level + 1][i], sizes[level][j], 0);
423
424 prolongation_sparsities[level]->collect_sizes();
425
427 mg_dof.begin(level);
428 cell != mg_dof.end(level);
429 ++cell)
430 if (cell->has_children())
431 {
432 cell->get_mg_dof_indices(dof_indices_parent);
433
434 for (unsigned int child = 0; child < cell->n_children(); ++child)
435 {
436 // set an alias to the
437 // prolongation matrix for
438 // this child
439 const FullMatrix<double> &prolongation =
441 child, cell->refinement_case());
442
443 cell->child(child)->get_mg_dof_indices(dof_indices_child);
444
445 // now tag the entries in the
446 // matrix which will be used
447 // for this pair of parent/child
448 for (unsigned int i = 0; i < dofs_per_cell; ++i)
449 for (unsigned int j = 0; j < dofs_per_cell; ++j)
450 if (prolongation(i, j) != 0)
451 {
452 const unsigned int icomp =
453 fe.system_to_component_index(i).first;
454 const unsigned int jcomp =
455 fe.system_to_component_index(j).first;
456 if ((icomp == jcomp) && mg_component_mask[icomp])
458 dof_indices_child[i], dof_indices_parent[j]);
459 }
460 }
461 }
462 prolongation_sparsities[level]->compress();
463
465 // now actually build the matrices
467 mg_dof.begin(level);
468 cell != mg_dof.end(level);
469 ++cell)
470 if (cell->has_children())
471 {
472 cell->get_mg_dof_indices(dof_indices_parent);
473
474 for (unsigned int child = 0; child < cell->n_children(); ++child)
475 {
476 // set an alias to the
477 // prolongation matrix for
478 // this child
479 const FullMatrix<double> &prolongation =
481 child, cell->refinement_case());
482
483 cell->child(child)->get_mg_dof_indices(dof_indices_child);
484
485 // now set the entries in the
486 // matrix
487 for (unsigned int i = 0; i < dofs_per_cell; ++i)
488 for (unsigned int j = 0; j < dofs_per_cell; ++j)
489 if (prolongation(i, j) != 0)
490 {
491 const unsigned int icomp =
492 fe.system_to_component_index(i).first;
493 const unsigned int jcomp =
494 fe.system_to_component_index(j).first;
495 if ((icomp == jcomp) && mg_component_mask[icomp])
497 dof_indices_child[i],
498 dof_indices_parent[j],
499 prolongation(i, j));
500 }
501 }
502 }
503 }
504 // impose boundary conditions
505 // but only in the column of
506 // the prolongation matrix
507 // TODO: this way is not very efficient
508
509 if (boundary_indices.size() != 0)
510 {
511 std::vector<std::vector<types::global_dof_index>> dofs_per_component(
512 mg_dof.get_triangulation().n_levels(),
513 std::vector<types::global_dof_index>(n_components));
514
516 dofs_per_component,
518 for (unsigned int level = 0; level < n_levels - 1; ++level)
519 {
520 if (boundary_indices[level].size() == 0)
521 continue;
522
523 for (unsigned int iblock = 0; iblock < n_components; ++iblock)
524 for (unsigned int jblock = 0; jblock < n_components; ++jblock)
525 if (iblock == jblock)
526 {
527 const types::global_dof_index n_dofs =
528 prolongation_matrices[level]->block(iblock, jblock).m();
529 for (types::global_dof_index i = 0; i < n_dofs; ++i)
530 {
533 ->block(iblock, jblock)
534 .begin(i),
536 ->block(iblock, jblock)
537 .end(i);
538 for (; begin != end; ++begin)
539 {
540 const types::global_dof_index column_number =
541 begin->column();
542
543 // convert global indices into local ones
544 const BlockIndices block_indices_coarse(
545 dofs_per_component[level]);
546 const types::global_dof_index global_j =
547 block_indices_coarse.local_to_global(iblock,
548 column_number);
549
550 std::set<types::global_dof_index>::const_iterator
551 found_dof = boundary_indices[level].find(global_j);
552
553 const bool is_boundary_index =
554 (found_dof != boundary_indices[level].end());
555
556 if (is_boundary_index)
557 {
559 ->block(iblock, jblock)
560 .set(i, column_number, 0);
561 }
562 }
563 }
564 }
565 }
566 }
567}
568
569
570
571template <typename number>
572template <int dim, int spacedim>
573void
575 const DoFHandler<dim, spacedim> & /*dof*/,
576 const DoFHandler<dim, spacedim> & mg_dof,
577 unsigned int select,
578 unsigned int mg_select,
579 const std::vector<unsigned int> & t_component,
580 const std::vector<unsigned int> & mg_t_component,
581 const std::vector<std::set<types::global_dof_index>> &bdry_indices)
582{
583 build(mg_dof, select, mg_select, t_component, mg_t_component, bdry_indices);
584}
585
586
587
588template <typename number>
589template <int dim, int spacedim>
590void
592 const DoFHandler<dim, spacedim> & mg_dof,
593 unsigned int select,
594 unsigned int mg_select,
595 const std::vector<unsigned int> & t_component,
596 const std::vector<unsigned int> & mg_t_component,
597 const std::vector<std::set<types::global_dof_index>> &bdry_indices)
598{
599 const FiniteElement<dim> &fe = mg_dof.get_fe();
600 unsigned int ncomp = mg_dof.get_fe(0).n_components();
601
602 target_component = t_component;
603 mg_target_component = mg_t_component;
604 boundary_indices = bdry_indices;
605
606 selected_component = select;
607 mg_selected_component = mg_select;
608
609 {
610 std::vector<bool> tmp(ncomp, false);
611 for (unsigned int c = 0; c < ncomp; ++c)
612 if (t_component[c] == selected_component)
613 tmp[c] = true;
614 component_mask = ComponentMask(tmp);
615 }
616
617 {
618 std::vector<bool> tmp(ncomp, false);
619 for (unsigned int c = 0; c < ncomp; ++c)
620 if (mg_t_component[c] == mg_selected_component)
621 tmp[c] = true;
622 mg_component_mask = ComponentMask(tmp);
623 }
624
625 // If components are renumbered,
626 // find the first original
627 // component corresponding to the
628 // target component.
629 for (unsigned int i = 0; i < target_component.size(); ++i)
630 {
631 if (target_component[i] == select)
632 {
633 selected_component = i;
634 break;
635 }
636 }
637
638 for (unsigned int i = 0; i < mg_target_component.size(); ++i)
639 {
640 if (mg_target_component[i] == mg_select)
641 {
642 mg_selected_component = i;
643 break;
644 }
645 }
646
648
649 interface_dofs.resize(mg_dof.get_triangulation().n_levels());
650 for (unsigned int l = 0; l < mg_dof.get_triangulation().n_levels(); ++l)
651 {
652 interface_dofs[l].clear();
653 interface_dofs[l].set_size(mg_dof.n_dofs(l));
654 }
655 MGTools::extract_inner_interface_dofs(mg_dof, interface_dofs);
656
657 // use a temporary vector to create the
658 // relation between global and level dofs
659 std::vector<types::global_dof_index> temp_copy_indices;
660 std::vector<types::global_dof_index> global_dof_indices(fe.n_dofs_per_cell());
661 std::vector<types::global_dof_index> level_dof_indices(fe.n_dofs_per_cell());
662 for (int level = mg_dof.get_triangulation().n_levels() - 1; level >= 0;
663 --level)
664 {
665 copy_to_and_from_indices[level].clear();
667 mg_dof.begin_active(level);
668 const typename DoFHandler<dim, spacedim>::active_cell_iterator level_end =
669 mg_dof.end_active(level);
670
671 temp_copy_indices.resize(0);
672 temp_copy_indices.resize(mg_dof.n_dofs(level),
674
675 // Compute coarse level right hand side
676 // by restricting from fine level.
677 for (; level_cell != level_end; ++level_cell)
678 {
679 // get the dof numbers of
680 // this cell for the global
681 // and the level-wise
682 // numbering
683 level_cell->get_dof_indices(global_dof_indices);
684 level_cell->get_mg_dof_indices(level_dof_indices);
685
686 for (unsigned int i = 0; i < fe.n_dofs_per_cell(); ++i)
687 {
688 const unsigned int component =
689 fe.system_to_component_index(i).first;
690 if (component_mask[component] &&
691 !interface_dofs[level].is_element(level_dof_indices[i]))
692 {
693 const types::global_dof_index level_start =
694 mg_component_start[level][mg_target_component[component]];
695 const types::global_dof_index global_start =
696 component_start[target_component[component]];
697 temp_copy_indices[level_dof_indices[i] - level_start] =
698 global_dof_indices[i] - global_start;
699 }
700 }
701 }
702
703 // write indices from vector into the map from
704 // global to level dofs
705 const types::global_dof_index n_active_dofs =
706 std::count_if(temp_copy_indices.begin(),
707 temp_copy_indices.end(),
708 [](const types::global_dof_index index) {
709 return index != numbers::invalid_dof_index;
710 });
711 copy_to_and_from_indices[level].resize(n_active_dofs);
712 types::global_dof_index counter = 0;
713 for (types::global_dof_index i = 0; i < temp_copy_indices.size(); ++i)
714 if (temp_copy_indices[i] != numbers::invalid_dof_index)
715 copy_to_and_from_indices[level][counter++] =
716 std::pair<types::global_dof_index, unsigned int>(
717 temp_copy_indices[i], i);
718 Assert(counter == n_active_dofs, ExcInternalError());
719 }
720}
721
722
723
724// explicit instantiations
725#include "mg_transfer_component.inst"
726
727
size_type local_to_global(const unsigned int block, const size_type index) const
bool represents_n_components(const unsigned int n) const
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
types::global_dof_index n_dofs() const
cell_iterator begin(const unsigned int level=0) const
unsigned int n_dofs_per_cell() const
unsigned int n_components() const
std::pair< unsigned int, unsigned int > system_to_component_index(const unsigned int index) const
virtual const FullMatrix< double > & get_prolongation_matrix(const unsigned int child, const RefinementCase< dim > &refinement_case=RefinementCase< dim >::isotropic_refinement) const
std::vector< unsigned int > mg_target_component
std::vector< unsigned int > target_component
void build_matrices(const DoFHandler< dim, spacedim > &dof, const DoFHandler< dim, spacedim > &mg_dof)
std::vector< types::global_dof_index > component_start
std::vector< std::shared_ptr< BlockSparseMatrix< double > > > prolongation_matrices
std::vector< std::vector< std::pair< types::global_dof_index, unsigned int > > > copy_to_and_from_indices
std::vector< std::vector< types::global_dof_index > > mg_component_start
std::vector< std::shared_ptr< BlockSparsityPattern > > prolongation_sparsities
void build(const DoFHandler< dim, spacedim > &dof_handler)
std::vector< std::vector< types::global_dof_index > > sizes
std::vector< std::set< types::global_dof_index > > boundary_indices
void do_copy_to_mg(const DoFHandler< dim, spacedim > &mg_dof, MGLevelObject< Vector< number > > &dst, const InVector &src) const
void build_matrices(const DoFHandler< dim, spacedim > &dof, const DoFHandler< dim, spacedim > &mg_dof, unsigned int selected, unsigned int mg_selected, const std::vector< unsigned int > &target_component=std::vector< unsigned int >(), const std::vector< unsigned int > &mg_target_component=std::vector< unsigned int >(), const std::vector< std::set< types::global_dof_index > > &boundary_indices=std::vector< std::set< types::global_dof_index > >())
void build(const DoFHandler< dim, spacedim > &dof, unsigned int selected, unsigned int mg_selected, const std::vector< unsigned int > &target_component=std::vector< unsigned int >(), const std::vector< unsigned int > &mg_target_component=std::vector< unsigned int >(), const std::vector< std::set< types::global_dof_index > > &boundary_indices=std::vector< std::set< types::global_dof_index > >())
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
#define Assert(cond, exc)
Definition: exceptions.h:1465
#define AssertIndexRange(index, range)
Definition: exceptions.h:1690
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcMessage(std::string arg1)
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 extract_inner_interface_dofs(const DoFHandler< dim, spacedim > &mg_dof_handler, std::vector< IndexSet > &interface_dofs)
Definition: mg_tools.cc:1480
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 > &)
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
VectorType::value_type * end(VectorType &V)
VectorType::value_type * begin(VectorType &V)
const types::global_dof_index invalid_dof_index
Definition: types.h:211