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