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_component.cc
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2003 - 2023 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{
259 // Fill target component with
260 // standard values (identity) if it
261 // is empty
262 if (target_component.size() == 0)
263 {
264 target_component.resize(mg_dof.get_fe(0).n_components());
265 for (unsigned int i = 0; i < target_component.size(); ++i)
266 target_component[i] = i;
267 }
268 else
269 {
270 // otherwise, check it for consistency
271 Assert(target_component.size() == mg_dof.get_fe(0).n_components(),
273 mg_dof.get_fe(0).n_components()));
274
275 for (unsigned int i = 0; i < target_component.size(); ++i)
276 {
278 }
279 }
280 // Do the same for the multilevel
281 // components. These may be
282 // different.
283 if (mg_target_component.size() == 0)
284 {
285 mg_target_component.resize(mg_dof.get_fe(0).n_components());
286 for (unsigned int i = 0; i < mg_target_component.size(); ++i)
288 }
289 else
290 {
291 Assert(mg_target_component.size() == mg_dof.get_fe(0).n_components(),
293 mg_dof.get_fe(0).n_components()));
294
295 for (unsigned int i = 0; i < mg_target_component.size(); ++i)
296 {
298 }
299 }
300
301 const FiniteElement<dim> &fe = mg_dof.get_fe();
302
303 // Effective number of components
304 // is the maximum entry in
305 // mg_target_component. This
306 // assumes that the values in that
307 // vector don't have holes.
308 const unsigned int n_components =
309 *std::max_element(mg_target_component.begin(), mg_target_component.end()) +
310 1;
311 const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
312 const unsigned int n_levels = mg_dof.get_triangulation().n_levels();
313
315 ExcMessage("Component mask has wrong size."));
316
317 // Compute the lengths of all blocks
318 sizes.resize(n_levels);
319 for (unsigned int l = 0; l < n_levels; ++l)
320 sizes[l].resize(n_components);
321
323
324 // Fill some index vectors
325 // for later use.
327 // Compute start indices from sizes
328 for (auto &level_components : mg_component_start)
329 {
331 for (types::global_dof_index &level_component_start : level_components)
332 {
333 const types::global_dof_index t = level_component_start;
334 level_component_start = k;
335 k += t;
336 }
337 }
338
340
342 for (types::global_dof_index &first_index : component_start)
343 {
344 const types::global_dof_index t = first_index;
345 first_index = k;
346 k += t;
347 }
348
349 // Build index vectors for
350 // copy_to_mg and
351 // copy_from_mg. These vectors must
352 // be prebuilt, since the
353 // get_dof_indices functions are
354 // too slow
355
356 copy_to_and_from_indices.resize(n_levels);
357
358 // Building the prolongation matrices starts here!
359
360 // reset the size of the array of
361 // matrices. call resize(0) first,
362 // in order to delete all elements
363 // and clear their memory. then
364 // repopulate these arrays
365 //
366 // note that on resize(0), the
367 // shared_ptr class takes care of
368 // deleting the object it points to
369 // by itself
370 prolongation_matrices.resize(0);
371 prolongation_sparsities.resize(0);
372 prolongation_matrices.reserve(n_levels - 1);
373 prolongation_sparsities.reserve(n_levels - 1);
374
375 for (unsigned int i = 0; i < n_levels - 1; ++i)
376 {
379 }
380
381 // two fields which will store the
382 // indices of the multigrid dofs
383 // for a cell and one of its children
384 std::vector<types::global_dof_index> dof_indices_parent(dofs_per_cell);
385 std::vector<types::global_dof_index> dof_indices_child(dofs_per_cell);
386
387 // for each level: first build the
388 // sparsity pattern of the matrices
389 // and then build the matrices
390 // themselves. note that we only
391 // need to take care of cells on
392 // the coarser level which have
393 // children
394 for (unsigned int level = 0; level < n_levels - 1; ++level)
395 {
396 // reset the dimension of the
397 // structure. note that for
398 // the number of entries per
399 // row, the number of parent
400 // dofs coupling to a child dof
401 // is necessary. this, is the
402 // number of degrees of freedom
403 // per cell
404 prolongation_sparsities[level]->reinit(n_components, n_components);
405 for (unsigned int i = 0; i < n_components; ++i)
406 for (unsigned int j = 0; j < n_components; ++j)
407 if (i == j)
408 prolongation_sparsities[level]->block(i, j).reinit(
409 sizes[level + 1][i], sizes[level][j], dofs_per_cell + 1);
410 else
411 prolongation_sparsities[level]->block(i, j).reinit(
412 sizes[level + 1][i], sizes[level][j], 0);
413
414 prolongation_sparsities[level]->collect_sizes();
415
417 mg_dof.begin(level);
418 cell != mg_dof.end(level);
419 ++cell)
420 if (cell->has_children())
421 {
422 cell->get_mg_dof_indices(dof_indices_parent);
423
424 for (unsigned int child = 0; child < cell->n_children(); ++child)
425 {
426 // set an alias to the
427 // prolongation matrix for
428 // this child
429 const FullMatrix<double> &prolongation =
431 child, cell->refinement_case());
432
433 cell->child(child)->get_mg_dof_indices(dof_indices_child);
434
435 // now tag the entries in the
436 // matrix which will be used
437 // for this pair of parent/child
438 for (unsigned int i = 0; i < dofs_per_cell; ++i)
439 for (unsigned int j = 0; j < dofs_per_cell; ++j)
440 if (prolongation(i, j) != 0)
441 {
442 const unsigned int icomp =
443 fe.system_to_component_index(i).first;
444 const unsigned int jcomp =
445 fe.system_to_component_index(j).first;
446 if ((icomp == jcomp) && mg_component_mask[icomp])
448 dof_indices_child[i], dof_indices_parent[j]);
449 }
450 }
451 }
452 prolongation_sparsities[level]->compress();
453
455 // now actually build the matrices
457 mg_dof.begin(level);
458 cell != mg_dof.end(level);
459 ++cell)
460 if (cell->has_children())
461 {
462 cell->get_mg_dof_indices(dof_indices_parent);
463
464 for (unsigned int child = 0; child < cell->n_children(); ++child)
465 {
466 // set an alias to the
467 // prolongation matrix for
468 // this child
469 const FullMatrix<double> &prolongation =
471 child, cell->refinement_case());
472
473 cell->child(child)->get_mg_dof_indices(dof_indices_child);
474
475 // now set the entries in the
476 // matrix
477 for (unsigned int i = 0; i < dofs_per_cell; ++i)
478 for (unsigned int j = 0; j < dofs_per_cell; ++j)
479 if (prolongation(i, j) != 0)
480 {
481 const unsigned int icomp =
482 fe.system_to_component_index(i).first;
483 const unsigned int jcomp =
484 fe.system_to_component_index(j).first;
485 if ((icomp == jcomp) && mg_component_mask[icomp])
487 dof_indices_child[i],
488 dof_indices_parent[j],
489 prolongation(i, j));
490 }
491 }
492 }
493 }
494 // impose boundary conditions
495 // but only in the column of
496 // the prolongation matrix
497 // TODO: this way is not very efficient
498
499 if (boundary_indices.size() != 0)
500 {
501 std::vector<std::vector<types::global_dof_index>> dofs_per_component(
502 mg_dof.get_triangulation().n_levels(),
503 std::vector<types::global_dof_index>(n_components));
504
506 dofs_per_component,
508 for (unsigned int level = 0; level < n_levels - 1; ++level)
509 {
510 if (boundary_indices[level].size() == 0)
511 continue;
512
513 for (unsigned int iblock = 0; iblock < n_components; ++iblock)
514 for (unsigned int jblock = 0; jblock < n_components; ++jblock)
515 if (iblock == jblock)
516 {
517 const types::global_dof_index n_dofs =
518 prolongation_matrices[level]->block(iblock, jblock).m();
519 for (types::global_dof_index i = 0; i < n_dofs; ++i)
520 {
523 ->block(iblock, jblock)
524 .begin(i),
526 ->block(iblock, jblock)
527 .end(i);
528 for (; begin != end; ++begin)
529 {
530 const types::global_dof_index column_number =
531 begin->column();
532
533 // convert global indices into local ones
534 const BlockIndices block_indices_coarse(
535 dofs_per_component[level]);
536 const types::global_dof_index global_j =
537 block_indices_coarse.local_to_global(iblock,
538 column_number);
539
540 std::set<types::global_dof_index>::const_iterator
541 found_dof = boundary_indices[level].find(global_j);
542
543 const bool is_boundary_index =
544 (found_dof != boundary_indices[level].end());
545
546 if (is_boundary_index)
547 {
549 ->block(iblock, jblock)
550 .set(i, column_number, 0);
551 }
552 }
553 }
554 }
555 }
556 }
557}
558
559
560
561template <typename number>
562template <int dim, int spacedim>
563void
565 const DoFHandler<dim, spacedim> & mg_dof,
566 unsigned int select,
567 unsigned int mg_select,
568 const std::vector<unsigned int> & t_component,
569 const std::vector<unsigned int> & mg_t_component,
570 const std::vector<std::set<types::global_dof_index>> &bdry_indices)
571{
572 const FiniteElement<dim> &fe = mg_dof.get_fe();
573 unsigned int ncomp = mg_dof.get_fe(0).n_components();
574
575 target_component = t_component;
576 mg_target_component = mg_t_component;
577 boundary_indices = bdry_indices;
578
579 selected_component = select;
580 mg_selected_component = mg_select;
581
582 {
583 std::vector<bool> tmp(ncomp, false);
584 for (unsigned int c = 0; c < ncomp; ++c)
585 if (t_component[c] == selected_component)
586 tmp[c] = true;
587 component_mask = ComponentMask(tmp);
588 }
589
590 {
591 std::vector<bool> tmp(ncomp, false);
592 for (unsigned int c = 0; c < ncomp; ++c)
593 if (mg_t_component[c] == mg_selected_component)
594 tmp[c] = true;
595 mg_component_mask = ComponentMask(tmp);
596 }
597
598 // If components are renumbered,
599 // find the first original
600 // component corresponding to the
601 // target component.
602 for (unsigned int i = 0; i < target_component.size(); ++i)
603 {
604 if (target_component[i] == select)
605 {
606 selected_component = i;
607 break;
608 }
609 }
610
611 for (unsigned int i = 0; i < mg_target_component.size(); ++i)
612 {
613 if (mg_target_component[i] == mg_select)
614 {
615 mg_selected_component = i;
616 break;
617 }
618 }
619
621
622 interface_dofs.resize(mg_dof.get_triangulation().n_levels());
623 for (unsigned int l = 0; l < mg_dof.get_triangulation().n_levels(); ++l)
624 {
625 interface_dofs[l].clear();
626 interface_dofs[l].set_size(mg_dof.n_dofs(l));
627 }
628 MGTools::extract_inner_interface_dofs(mg_dof, interface_dofs);
629
630 // use a temporary vector to create the
631 // relation between global and level dofs
632 std::vector<types::global_dof_index> temp_copy_indices;
633 std::vector<types::global_dof_index> global_dof_indices(fe.n_dofs_per_cell());
634 std::vector<types::global_dof_index> level_dof_indices(fe.n_dofs_per_cell());
635 for (int level = mg_dof.get_triangulation().n_levels() - 1; level >= 0;
636 --level)
637 {
638 copy_to_and_from_indices[level].clear();
640 mg_dof.begin_active(level);
641 const typename DoFHandler<dim, spacedim>::active_cell_iterator level_end =
642 mg_dof.end_active(level);
643
644 temp_copy_indices.resize(0);
645 temp_copy_indices.resize(mg_dof.n_dofs(level),
647
648 // Compute coarse level right hand side
649 // by restricting from fine level.
650 for (; level_cell != level_end; ++level_cell)
651 {
652 // get the dof numbers of
653 // this cell for the global
654 // and the level-wise
655 // numbering
656 level_cell->get_dof_indices(global_dof_indices);
657 level_cell->get_mg_dof_indices(level_dof_indices);
658
659 for (unsigned int i = 0; i < fe.n_dofs_per_cell(); ++i)
660 {
661 const unsigned int component =
662 fe.system_to_component_index(i).first;
663 if (component_mask[component] &&
664 !interface_dofs[level].is_element(level_dof_indices[i]))
665 {
666 const types::global_dof_index level_start =
667 mg_component_start[level][mg_target_component[component]];
668 const types::global_dof_index global_start =
669 component_start[target_component[component]];
670 temp_copy_indices[level_dof_indices[i] - level_start] =
671 global_dof_indices[i] - global_start;
672 }
673 }
674 }
675
676 // write indices from vector into the map from
677 // global to level dofs
678 const types::global_dof_index n_active_dofs =
679 std::count_if(temp_copy_indices.begin(),
680 temp_copy_indices.end(),
681 [](const types::global_dof_index index) {
682 return index != numbers::invalid_dof_index;
683 });
684 copy_to_and_from_indices[level].resize(n_active_dofs);
685 types::global_dof_index counter = 0;
686 for (types::global_dof_index i = 0; i < temp_copy_indices.size(); ++i)
687 if (temp_copy_indices[i] != numbers::invalid_dof_index)
688 copy_to_and_from_indices[level][counter++] =
689 std::pair<types::global_dof_index, unsigned int>(
690 temp_copy_indices[i], i);
691 Assert(counter == n_active_dofs, ExcInternalError());
692 }
693}
694
695
696
697// explicit instantiations
698#include "mg_transfer_component.inst"
699
700
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 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
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
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(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:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
unsigned int level
Definition grid_out.cc:4618
#define Assert(cond, exc)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcMessage(std::string arg1)
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 extract_inner_interface_dofs(const DoFHandler< dim, spacedim > &mg_dof_handler, std::vector< IndexSet > &interface_dofs)
Definition mg_tools.cc:1446
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:1157
const types::global_dof_index invalid_dof_index
Definition types.h:233