Reference documentation for deal.II version 9.2.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\}}\)
mg_transfer_matrix_free.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2016 - 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 #ifndef dealii_mg_transfer_matrix_free_h
17 #define dealii_mg_transfer_matrix_free_h
18 
19 #include <deal.II/base/config.h>
20 
23 
25 
28 
33 
34 
36 
37 
40 
56 template <int dim, typename Number>
58  : public MGLevelGlobalTransfer<LinearAlgebra::distributed::Vector<Number>>
59 {
60 public:
66 
72 
76  virtual ~MGTransferMatrixFree() override = default;
77 
81  void
82  initialize_constraints(const MGConstrainedDoFs &mg_constrained_dofs);
83 
87  void
88  clear();
89 
105  void
106  build(const DoFHandler<dim, dim> &dof_handler,
107  const std::vector<std::shared_ptr<const Utilities::MPI::Partitioner>>
108  &external_partitioners =
109  std::vector<std::shared_ptr<const Utilities::MPI::Partitioner>>());
110 
125  virtual void
126  prolongate(
127  const unsigned int to_level,
129  const LinearAlgebra::distributed::Vector<Number> &src) const override;
130 
149  virtual void
151  const unsigned int from_level,
153  const LinearAlgebra::distributed::Vector<Number> &src) const override;
154 
165  template <typename Number2, int spacedim>
166  void
168  const DoFHandler<dim, spacedim> & dof_handler,
171 
176 
180  std::size_t
181  memory_consumption() const;
182 
183 private:
189  unsigned int fe_degree;
190 
196 
201  unsigned int n_components;
202 
208  unsigned int n_child_cell_dofs;
209 
220  std::vector<std::vector<unsigned int>> level_dof_indices;
221 
226  std::vector<std::vector<std::pair<unsigned int, unsigned int>>>
228 
233  std::vector<unsigned int> n_owned_level_cells;
234 
240 
245 
257  std::vector<AlignedVector<VectorizedArray<Number>>> weights_on_refined;
258 
264  std::vector<std::vector<std::vector<unsigned short>>> dirichlet_indices;
265 
274 
278  template <int degree>
279  void
281  const unsigned int to_level,
284 
288  template <int degree>
289  void
290  do_restrict_add(const unsigned int from_level,
293 };
294 
295 
312 template <int dim, typename Number>
314  : public MGTransferBase<LinearAlgebra::distributed::BlockVector<Number>>
315 {
316 public:
321  MGTransferBlockMatrixFree() = default;
322 
328 
333  const std::vector<MGConstrainedDoFs> &mg_constrained_dofs);
334 
338  virtual ~MGTransferBlockMatrixFree() override = default;
339 
343  void
344  initialize_constraints(const MGConstrainedDoFs &mg_constrained_dofs);
345 
349  void
351  const std::vector<MGConstrainedDoFs> &mg_constrained_dofs);
352 
356  void
357  clear();
358 
362  void
363  build(const DoFHandler<dim, dim> &dof_handler);
364 
368  void
369  build(const std::vector<const DoFHandler<dim, dim> *> &dof_handler);
370 
385  virtual void
386  prolongate(
387  const unsigned int to_level,
389  const LinearAlgebra::distributed::BlockVector<Number> &src) const override;
390 
409  virtual void
411  const unsigned int from_level,
413  const LinearAlgebra::distributed::BlockVector<Number> &src) const override;
414 
425  template <typename Number2, int spacedim>
426  void
427  copy_to_mg(
428  const DoFHandler<dim, spacedim> & dof_handler,
431 
435  template <typename Number2, int spacedim>
436  void
437  copy_to_mg(
438  const std::vector<const DoFHandler<dim, spacedim> *> & dof_handler,
441 
445  template <typename Number2, int spacedim>
446  void
447  copy_from_mg(
448  const DoFHandler<dim, spacedim> & dof_handler,
451  const;
452 
456  template <typename Number2, int spacedim>
457  void
458  copy_from_mg(
459  const std::vector<const DoFHandler<dim, spacedim> *> &dof_handler,
462  const;
463 
467  std::size_t
468  memory_consumption() const;
469 
474  static const bool supports_dof_handler_vector = true;
475 
476 private:
480  std::vector<MGTransferMatrixFree<dim, Number>> matrix_free_transfer_vector;
481 
486  const bool same_for_all;
487 };
488 
489 
493 //------------------------ templated functions -------------------------
494 #ifndef DOXYGEN
495 
496 
497 template <int dim, typename Number>
498 template <typename Number2, int spacedim>
499 void
501  const DoFHandler<dim, spacedim> & dof_handler,
504 {
505  const unsigned int min_level = dst.min_level();
506  const unsigned int max_level = dst.max_level();
507 
508  Assert(max_level == dof_handler.get_triangulation().n_global_levels() - 1,
510  max_level, dof_handler.get_triangulation().n_global_levels() - 1));
511 
513  (dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
514  &dof_handler.get_triangulation()));
515  MPI_Comm mpi_communicator =
516  p_tria != nullptr ? p_tria->get_communicator() : MPI_COMM_SELF;
517 
518  // resize the dst vector if it's empty or has incorrect size
519  MGLevelObject<IndexSet> relevant_dofs(min_level, max_level);
520  for (unsigned int level = min_level; level <= max_level; ++level)
521  {
523  level,
524  relevant_dofs[level]);
525  if (dst[level].size() !=
526  dof_handler.locally_owned_mg_dofs(level).size() ||
527  dst[level].local_size() !=
528  dof_handler.locally_owned_mg_dofs(level).n_elements())
529  dst[level].reinit(dof_handler.locally_owned_mg_dofs(level),
530  relevant_dofs[level],
531  mpi_communicator);
532  }
533 
534  const FiniteElement<dim, spacedim> &fe = dof_handler.get_fe();
535 
536  // copy fine level vector to active cells in MG hierarchy
537  this->copy_to_mg(dof_handler, dst, src, true);
538 
539  // FIXME: maybe need to store hanging nodes constraints per level?
540  // MGConstrainedDoFs does NOT keep this info right now, only periodicity
541  // constraints...
542  dst[max_level].update_ghost_values();
543  // do the transfer from level to level-1:
544  for (unsigned int level = max_level; level > min_level; --level)
545  {
546  // auxiliary vector which always has ghost elements
548  dof_handler.locally_owned_mg_dofs(level),
549  relevant_dofs[level],
550  mpi_communicator);
551  ghosted_vector = dst[level];
552  ghosted_vector.update_ghost_values();
553 
554  std::vector<Number> dof_values_coarse(fe.dofs_per_cell);
555  Vector<Number> dof_values_fine(fe.dofs_per_cell);
556  Vector<Number> tmp(fe.dofs_per_cell);
557  std::vector<types::global_dof_index> dof_indices(fe.dofs_per_cell);
558  typename DoFHandler<dim>::cell_iterator cell =
559  dof_handler.begin(level - 1);
560  typename DoFHandler<dim>::cell_iterator endc = dof_handler.end(level - 1);
561  for (; cell != endc; ++cell)
562  if (cell->is_locally_owned_on_level())
563  {
564  // if we get to a cell without children (== active), we can
565  // skip it as there values should be already set by the
566  // equivalent of copy_to_mg()
567  if (cell->is_active())
568  continue;
569 
570  std::fill(dof_values_coarse.begin(), dof_values_coarse.end(), 0.);
571  for (unsigned int child = 0; child < cell->n_children(); ++child)
572  {
573  cell->child(child)->get_mg_dof_indices(dof_indices);
574  for (unsigned int i = 0; i < fe.dofs_per_cell; ++i)
575  dof_values_fine(i) = ghosted_vector(dof_indices[i]);
576  fe.get_restriction_matrix(child, cell->refinement_case())
577  .vmult(tmp, dof_values_fine);
578  for (unsigned int i = 0; i < fe.dofs_per_cell; ++i)
579  if (fe.restriction_is_additive(i))
580  dof_values_coarse[i] += tmp[i];
581  else if (tmp(i) != 0.)
582  dof_values_coarse[i] = tmp[i];
583  }
584  cell->get_mg_dof_indices(dof_indices);
585  for (unsigned int i = 0; i < fe.dofs_per_cell; ++i)
586  dst[level - 1](dof_indices[i]) = dof_values_coarse[i];
587  }
588 
589  dst[level - 1].compress(VectorOperation::insert);
590  dst[level - 1].update_ghost_values();
591  }
592 }
593 
594 
595 
596 template <int dim, typename Number>
597 template <typename Number2, int spacedim>
598 void
600  const DoFHandler<dim, spacedim> & dof_handler,
603 {
604  AssertDimension(matrix_free_transfer_vector.size(), 1);
605  Assert(same_for_all,
606  ExcMessage(
607  "This object was initialized with support for usage with one "
608  "DoFHandler for each block, but this method assumes that "
609  "the same DoFHandler is used for all the blocks!"));
610  const std::vector<const DoFHandler<dim, spacedim> *> mg_dofs(src.n_blocks(),
611  &dof_handler);
612 
613  copy_to_mg(mg_dofs, dst, src);
614 }
615 
616 template <int dim, typename Number>
617 template <typename Number2, int spacedim>
618 void
620  const std::vector<const DoFHandler<dim, spacedim> *> & dof_handler,
623 {
624  const unsigned int n_blocks = src.n_blocks();
625  AssertDimension(dof_handler.size(), n_blocks);
626 
627  if (n_blocks == 0)
628  return;
629 
630  const unsigned int min_level = dst.min_level();
631  const unsigned int max_level = dst.max_level();
632 
633  // this function is normally called within the Multigrid class with
634  // dst == defect level block vector. At first run this vector is not
635  // initialized. Do this below:
636  {
638  (dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
639  &(dof_handler[0]->get_triangulation())));
640  for (unsigned int i = 1; i < n_blocks; ++i)
641  AssertThrow(
642  (dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
643  &(dof_handler[0]->get_triangulation())) == tria),
644  ExcMessage("The DoFHandler use different Triangulations!"));
645 
646  MGLevelObject<bool> do_reinit;
647  do_reinit.resize(min_level, max_level);
648  for (unsigned int level = min_level; level <= max_level; ++level)
649  {
650  do_reinit[level] = false;
651  if (dst[level].n_blocks() != n_blocks)
652  {
653  do_reinit[level] = true;
654  continue; // level
655  }
656  for (unsigned int b = 0; b < n_blocks; ++b)
657  {
659  if (v.size() !=
660  dof_handler[b]->locally_owned_mg_dofs(level).size() ||
661  v.local_size() !=
662  dof_handler[b]->locally_owned_mg_dofs(level).n_elements())
663  {
664  do_reinit[level] = true;
665  break; // b
666  }
667  }
668  }
669 
670  for (unsigned int level = min_level; level <= max_level; ++level)
671  {
672  if (do_reinit[level])
673  {
674  dst[level].reinit(n_blocks);
675  for (unsigned int b = 0; b < n_blocks; ++b)
676  {
678  dst[level].block(b);
679  v.reinit(dof_handler[b]->locally_owned_mg_dofs(level),
680  tria != nullptr ? tria->get_communicator() :
681  MPI_COMM_SELF);
682  }
683  dst[level].collect_sizes();
684  }
685  else
686  dst[level] = 0;
687  }
688  }
689 
690  // FIXME: this a quite ugly as we need a temporary object:
692  min_level, max_level);
693 
694  for (unsigned int b = 0; b < n_blocks; ++b)
695  {
696  for (unsigned int l = min_level; l <= max_level; ++l)
697  dst_non_block[l].reinit(dst[l].block(b));
698  const unsigned int data_block = same_for_all ? 0 : b;
699  matrix_free_transfer_vector[data_block].copy_to_mg(*dof_handler[b],
700  dst_non_block,
701  src.block(b));
702 
703  for (unsigned int l = min_level; l <= max_level; ++l)
704  dst[l].block(b) = dst_non_block[l];
705  }
706 }
707 
708 template <int dim, typename Number>
709 template <typename Number2, int spacedim>
710 void
712  const DoFHandler<dim, spacedim> & dof_handler,
715  const
716 {
717  AssertDimension(matrix_free_transfer_vector.size(), 1);
718  const std::vector<const DoFHandler<dim, spacedim> *> mg_dofs(dst.n_blocks(),
719  &dof_handler);
720 
721  copy_from_mg(mg_dofs, dst, src);
722 }
723 
724 template <int dim, typename Number>
725 template <typename Number2, int spacedim>
726 void
728  const std::vector<const DoFHandler<dim, spacedim> *> &dof_handler,
731  const
732 {
733  const unsigned int n_blocks = dst.n_blocks();
734  AssertDimension(dof_handler.size(), n_blocks);
735 
736  if (n_blocks == 0)
737  return;
738 
739  const unsigned int min_level = src.min_level();
740  const unsigned int max_level = src.max_level();
741 
742  for (unsigned int l = min_level; l <= max_level; ++l)
743  AssertDimension(src[l].n_blocks(), dst.n_blocks());
744 
745  // FIXME: this a quite ugly as we need a temporary object:
747  min_level, max_level);
748 
749  for (unsigned int b = 0; b < n_blocks; ++b)
750  {
751  for (unsigned int l = min_level; l <= max_level; ++l)
752  {
753  src_non_block[l].reinit(src[l].block(b));
754  src_non_block[l] = src[l].block(b);
755  }
756  const unsigned int data_block = same_for_all ? 0 : b;
757  matrix_free_transfer_vector[data_block].copy_from_mg(*dof_handler[b],
758  dst.block(b),
759  src_non_block);
760  }
761 }
762 
763 
764 
765 #endif // DOXYGEN
766 
767 
769 
770 #endif
const Triangulation< dim, spacedim > & get_triangulation() const
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
Definition: matrix_block.h:621
bool restriction_is_additive(const unsigned int index) const
Definition: fe.h:3242
void do_restrict_add(const unsigned int from_level, LinearAlgebra::distributed::Vector< Number > &dst, const LinearAlgebra::distributed::Vector< Number > &src) const
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1579
std::vector< unsigned int > n_owned_level_cells
cell_iterator begin(const unsigned int level=0) const
Definition: dof_handler.cc:922
void copy_from_mg(const DoFHandler< dim, spacedim > &dof_handler, LinearAlgebra::distributed::Vector< Number2 > &dst, const MGLevelObject< LinearAlgebra::distributed::Vector< Number >> &src) const
void build(const DoFHandler< dim, dim > &dof_handler, const std::vector< std::shared_ptr< const Utilities::MPI::Partitioner >> &external_partitioners=std::vector< std::shared_ptr< const Utilities::MPI::Partitioner >>())
AlignedVector< VectorizedArray< Number > > evaluation_data
static ::ExceptionBase & ExcNoProlongation()
std::size_t memory_consumption() const
std::vector< std::vector< std::vector< unsigned short > > > dirichlet_indices
cell_iterator end() const
Definition: dof_handler.cc:951
void reinit(const size_type size, const bool omit_zeroing_entries=false)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1531
const FiniteElement< dim, spacedim > & get_fe(const unsigned int index=0) const
void extract_locally_relevant_level_dofs(const DoFHandlerType &dof_handler, const unsigned int level, IndexSet &dof_set)
Definition: dof_tools.cc:1215
void interpolate_to_mg(const DoFHandler< dim, spacedim > &dof_handler, MGLevelObject< LinearAlgebra::distributed::Vector< Number >> &dst, const LinearAlgebra::distributed::Vector< Number2 > &src) const
size_type size() const
Definition: index_set.h:1635
static ::ExceptionBase & ExcMessage(std::string arg1)
AlignedVector< VectorizedArray< Number > > prolongation_matrix_1d
void copy_to_mg(const DoFHandler< dim, spacedim > &dof_handler, MGLevelObject< LinearAlgebra::distributed::BlockVector< Number >> &dst, const LinearAlgebra::distributed::BlockVector< Number2 > &src) const
#define Assert(cond, exc)
Definition: exceptions.h:1419
virtual const FullMatrix< double > & get_restriction_matrix(const unsigned int child, const RefinementCase< dim > &refinement_case=RefinementCase< dim >::isotropic_refinement) const
Definition: fe.cc:305
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
virtual void restrict_and_add(const unsigned int from_level, LinearAlgebra::distributed::Vector< Number > &dst, const LinearAlgebra::distributed::Vector< Number > &src) const override
virtual void prolongate(const unsigned int to_level, LinearAlgebra::distributed::Vector< Number > &dst, const LinearAlgebra::distributed::Vector< Number > &src) const override
#define DeclException0(Exception0)
Definition: exceptions.h:473
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:359
unsigned int level
Definition: grid_out.cc:4355
std::vector< AlignedVector< VectorizedArray< Number > > > weights_on_refined
std::enable_if< IsBlockVector< VectorType >::value, unsigned int >::type n_blocks(const VectorType &vector)
Definition: operators.h:49
const IndexSet & locally_owned_mg_dofs(const unsigned int level) const
std::vector< std::vector< unsigned int > > level_dof_indices
virtual const MPI_Comm & get_communicator() const
Definition: tria_base.cc:138
const unsigned int dofs_per_cell
Definition: fe_base.h:282
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
void copy_from_mg(const DoFHandler< dim, spacedim > &dof_handler, LinearAlgebra::distributed::BlockVector< Number2 > &dst, const MGLevelObject< LinearAlgebra::distributed::BlockVector< Number >> &src) const
virtual size_type size() const override
void do_prolongate_add(const unsigned int to_level, LinearAlgebra::distributed::Vector< Number > &dst, const LinearAlgebra::distributed::Vector< Number > &src) const
SmartPointer< const MGConstrainedDoFs, MGLevelGlobalTransfer< LinearAlgebra::distributed::Vector< Number > > > mg_constrained_dofs
Definition: mg_transfer.h:618
std::vector< std::vector< std::pair< unsigned int, unsigned int > > > parent_child_connect
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:358
std::vector< MGTransferMatrixFree< dim, Number > > matrix_free_transfer_vector
virtual ~MGTransferMatrixFree() override=default
void resize(const unsigned int new_minlevel, const unsigned int new_maxlevel)
MGLevelObject< std::shared_ptr< const Utilities::MPI::Partitioner > > vector_partitioners
void initialize_constraints(const MGConstrainedDoFs &mg_constrained_dofs)
void copy_to_mg(const DoFHandler< dim, spacedim > &dof_handler, MGLevelObject< LinearAlgebra::distributed::Vector< Number >> &dst, const LinearAlgebra::distributed::Vector< Number2 > &src) const
BlockType & block(const unsigned int i)
size_type n_elements() const
Definition: index_set.h:1833
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)