Reference documentation for deal.II version Git e038590ef8 2020-10-28 14:38:53 +0100
\(\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 
53 template <int dim, typename Number>
55  : public MGLevelGlobalTransfer<LinearAlgebra::distributed::Vector<Number>>
56 {
57 public:
63 
69 
73  virtual ~MGTransferMatrixFree() override = default;
74 
78  void
79  initialize_constraints(const MGConstrainedDoFs &mg_constrained_dofs);
80 
84  void
85  clear();
86 
102  void
103  build(const DoFHandler<dim, dim> &dof_handler,
104  const std::vector<std::shared_ptr<const Utilities::MPI::Partitioner>>
105  &external_partitioners =
106  std::vector<std::shared_ptr<const Utilities::MPI::Partitioner>>());
107 
122  virtual void
123  prolongate(
124  const unsigned int to_level,
126  const LinearAlgebra::distributed::Vector<Number> &src) const override;
127 
146  virtual void
148  const unsigned int from_level,
150  const LinearAlgebra::distributed::Vector<Number> &src) const override;
151 
162  template <typename Number2, int spacedim>
163  void
165  const DoFHandler<dim, spacedim> & dof_handler,
168 
173 
177  std::size_t
178  memory_consumption() const;
179 
180 private:
186  unsigned int fe_degree;
187 
193 
198  unsigned int n_components;
199 
205  unsigned int n_child_cell_dofs;
206 
217  std::vector<std::vector<unsigned int>> level_dof_indices;
218 
223  std::vector<std::vector<std::pair<unsigned int, unsigned int>>>
225 
230  std::vector<unsigned int> n_owned_level_cells;
231 
237 
242 
254  std::vector<AlignedVector<VectorizedArray<Number>>> weights_on_refined;
255 
261  std::vector<std::vector<std::vector<unsigned short>>> dirichlet_indices;
262 
271 
275  template <int degree>
276  void
278  const unsigned int to_level,
281 
285  template <int degree>
286  void
287  do_restrict_add(const unsigned int from_level,
290 };
291 
292 
306 template <int dim, typename Number>
308  : public MGTransferBase<LinearAlgebra::distributed::BlockVector<Number>>
309 {
310 public:
315  MGTransferBlockMatrixFree() = default;
316 
322 
327  const std::vector<MGConstrainedDoFs> &mg_constrained_dofs);
328 
332  virtual ~MGTransferBlockMatrixFree() override = default;
333 
337  void
338  initialize_constraints(const MGConstrainedDoFs &mg_constrained_dofs);
339 
343  void
345  const std::vector<MGConstrainedDoFs> &mg_constrained_dofs);
346 
350  void
351  clear();
352 
356  void
357  build(const DoFHandler<dim, dim> &dof_handler);
358 
362  void
363  build(const std::vector<const DoFHandler<dim, dim> *> &dof_handler);
364 
379  virtual void
380  prolongate(
381  const unsigned int to_level,
383  const LinearAlgebra::distributed::BlockVector<Number> &src) const override;
384 
403  virtual void
405  const unsigned int from_level,
407  const LinearAlgebra::distributed::BlockVector<Number> &src) const override;
408 
419  template <typename Number2, int spacedim>
420  void
421  copy_to_mg(
422  const DoFHandler<dim, spacedim> & dof_handler,
425 
429  template <typename Number2, int spacedim>
430  void
431  copy_to_mg(
432  const std::vector<const DoFHandler<dim, spacedim> *> & dof_handler,
435 
439  template <typename Number2, int spacedim>
440  void
441  copy_from_mg(
442  const DoFHandler<dim, spacedim> & dof_handler,
445  const;
446 
450  template <typename Number2, int spacedim>
451  void
452  copy_from_mg(
453  const std::vector<const DoFHandler<dim, spacedim> *> &dof_handler,
456  const;
457 
461  std::size_t
462  memory_consumption() const;
463 
468  static const bool supports_dof_handler_vector = true;
469 
470 private:
474  std::vector<MGTransferMatrixFree<dim, Number>> matrix_free_transfer_vector;
475 
480  const bool same_for_all;
481 };
482 
483 
487 //------------------------ templated functions -------------------------
488 #ifndef DOXYGEN
489 
490 
491 template <int dim, typename Number>
492 template <typename Number2, int spacedim>
493 void
495  const DoFHandler<dim, spacedim> & dof_handler,
498 {
499  const unsigned int min_level = dst.min_level();
500  const unsigned int max_level = dst.max_level();
501 
502  Assert(max_level == dof_handler.get_triangulation().n_global_levels() - 1,
504  max_level, dof_handler.get_triangulation().n_global_levels() - 1));
505 
507  (dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
508  &dof_handler.get_triangulation()));
509  MPI_Comm mpi_communicator =
510  p_tria != nullptr ? p_tria->get_communicator() : MPI_COMM_SELF;
511 
512  // resize the dst vector if it's empty or has incorrect size
513  MGLevelObject<IndexSet> relevant_dofs(min_level, max_level);
514  for (unsigned int level = min_level; level <= max_level; ++level)
515  {
517  level,
518  relevant_dofs[level]);
519  if (dst[level].size() !=
520  dof_handler.locally_owned_mg_dofs(level).size() ||
521  dst[level].local_size() !=
522  dof_handler.locally_owned_mg_dofs(level).n_elements())
523  dst[level].reinit(dof_handler.locally_owned_mg_dofs(level),
524  relevant_dofs[level],
525  mpi_communicator);
526  }
527 
528  const FiniteElement<dim, spacedim> &fe = dof_handler.get_fe();
529 
530  // copy fine level vector to active cells in MG hierarchy
531  this->copy_to_mg(dof_handler, dst, src, true);
532 
533  // FIXME: maybe need to store hanging nodes constraints per level?
534  // MGConstrainedDoFs does NOT keep this info right now, only periodicity
535  // constraints...
536  dst[max_level].update_ghost_values();
537  // do the transfer from level to level-1:
538  for (unsigned int level = max_level; level > min_level; --level)
539  {
540  // auxiliary vector which always has ghost elements
542  dof_handler.locally_owned_mg_dofs(level),
543  relevant_dofs[level],
544  mpi_communicator);
545  ghosted_vector = dst[level];
546  ghosted_vector.update_ghost_values();
547 
548  std::vector<Number> dof_values_coarse(fe.n_dofs_per_cell());
549  Vector<Number> dof_values_fine(fe.n_dofs_per_cell());
550  Vector<Number> tmp(fe.n_dofs_per_cell());
551  std::vector<types::global_dof_index> dof_indices(fe.n_dofs_per_cell());
552  typename DoFHandler<dim>::cell_iterator cell =
553  dof_handler.begin(level - 1);
554  typename DoFHandler<dim>::cell_iterator endc = dof_handler.end(level - 1);
555  for (; cell != endc; ++cell)
556  if (cell->is_locally_owned_on_level())
557  {
558  // if we get to a cell without children (== active), we can
559  // skip it as there values should be already set by the
560  // equivalent of copy_to_mg()
561  if (cell->is_active())
562  continue;
563 
564  std::fill(dof_values_coarse.begin(), dof_values_coarse.end(), 0.);
565  for (unsigned int child = 0; child < cell->n_children(); ++child)
566  {
567  cell->child(child)->get_mg_dof_indices(dof_indices);
568  for (unsigned int i = 0; i < fe.n_dofs_per_cell(); ++i)
569  dof_values_fine(i) = ghosted_vector(dof_indices[i]);
570  fe.get_restriction_matrix(child, cell->refinement_case())
571  .vmult(tmp, dof_values_fine);
572  for (unsigned int i = 0; i < fe.n_dofs_per_cell(); ++i)
573  if (fe.restriction_is_additive(i))
574  dof_values_coarse[i] += tmp[i];
575  else if (tmp(i) != 0.)
576  dof_values_coarse[i] = tmp[i];
577  }
578  cell->get_mg_dof_indices(dof_indices);
579  for (unsigned int i = 0; i < fe.n_dofs_per_cell(); ++i)
580  dst[level - 1](dof_indices[i]) = dof_values_coarse[i];
581  }
582 
583  dst[level - 1].compress(VectorOperation::insert);
584  dst[level - 1].update_ghost_values();
585  }
586 }
587 
588 
589 
590 template <int dim, typename Number>
591 template <typename Number2, int spacedim>
592 void
594  const DoFHandler<dim, spacedim> & dof_handler,
597 {
598  AssertDimension(matrix_free_transfer_vector.size(), 1);
599  Assert(same_for_all,
600  ExcMessage(
601  "This object was initialized with support for usage with one "
602  "DoFHandler for each block, but this method assumes that "
603  "the same DoFHandler is used for all the blocks!"));
604  const std::vector<const DoFHandler<dim, spacedim> *> mg_dofs(src.n_blocks(),
605  &dof_handler);
606 
607  copy_to_mg(mg_dofs, dst, src);
608 }
609 
610 template <int dim, typename Number>
611 template <typename Number2, int spacedim>
612 void
614  const std::vector<const DoFHandler<dim, spacedim> *> & dof_handler,
617 {
618  const unsigned int n_blocks = src.n_blocks();
619  AssertDimension(dof_handler.size(), n_blocks);
620 
621  if (n_blocks == 0)
622  return;
623 
624  const unsigned int min_level = dst.min_level();
625  const unsigned int max_level = dst.max_level();
626 
627  // this function is normally called within the Multigrid class with
628  // dst == defect level block vector. At first run this vector is not
629  // initialized. Do this below:
630  {
632  (dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
633  &(dof_handler[0]->get_triangulation())));
634  for (unsigned int i = 1; i < n_blocks; ++i)
635  AssertThrow(
636  (dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
637  &(dof_handler[0]->get_triangulation())) == tria),
638  ExcMessage("The DoFHandler use different Triangulations!"));
639 
640  MGLevelObject<bool> do_reinit;
641  do_reinit.resize(min_level, max_level);
642  for (unsigned int level = min_level; level <= max_level; ++level)
643  {
644  do_reinit[level] = false;
645  if (dst[level].n_blocks() != n_blocks)
646  {
647  do_reinit[level] = true;
648  continue; // level
649  }
650  for (unsigned int b = 0; b < n_blocks; ++b)
651  {
653  if (v.size() !=
654  dof_handler[b]->locally_owned_mg_dofs(level).size() ||
655  v.local_size() !=
656  dof_handler[b]->locally_owned_mg_dofs(level).n_elements())
657  {
658  do_reinit[level] = true;
659  break; // b
660  }
661  }
662  }
663 
664  for (unsigned int level = min_level; level <= max_level; ++level)
665  {
666  if (do_reinit[level])
667  {
668  dst[level].reinit(n_blocks);
669  for (unsigned int b = 0; b < n_blocks; ++b)
670  {
672  dst[level].block(b);
673  v.reinit(dof_handler[b]->locally_owned_mg_dofs(level),
674  tria != nullptr ? tria->get_communicator() :
675  MPI_COMM_SELF);
676  }
677  dst[level].collect_sizes();
678  }
679  else
680  dst[level] = 0;
681  }
682  }
683 
684  // FIXME: this a quite ugly as we need a temporary object:
686  min_level, max_level);
687 
688  for (unsigned int b = 0; b < n_blocks; ++b)
689  {
690  for (unsigned int l = min_level; l <= max_level; ++l)
691  dst_non_block[l].reinit(dst[l].block(b));
692  const unsigned int data_block = same_for_all ? 0 : b;
693  matrix_free_transfer_vector[data_block].copy_to_mg(*dof_handler[b],
694  dst_non_block,
695  src.block(b));
696 
697  for (unsigned int l = min_level; l <= max_level; ++l)
698  dst[l].block(b) = dst_non_block[l];
699  }
700 }
701 
702 template <int dim, typename Number>
703 template <typename Number2, int spacedim>
704 void
706  const DoFHandler<dim, spacedim> & dof_handler,
709  const
710 {
711  AssertDimension(matrix_free_transfer_vector.size(), 1);
712  const std::vector<const DoFHandler<dim, spacedim> *> mg_dofs(dst.n_blocks(),
713  &dof_handler);
714 
715  copy_from_mg(mg_dofs, dst, src);
716 }
717 
718 template <int dim, typename Number>
719 template <typename Number2, int spacedim>
720 void
722  const std::vector<const DoFHandler<dim, spacedim> *> &dof_handler,
725  const
726 {
727  const unsigned int n_blocks = dst.n_blocks();
728  AssertDimension(dof_handler.size(), n_blocks);
729 
730  if (n_blocks == 0)
731  return;
732 
733  const unsigned int min_level = src.min_level();
734  const unsigned int max_level = src.max_level();
735 
736  for (unsigned int l = min_level; l <= max_level; ++l)
737  AssertDimension(src[l].n_blocks(), dst.n_blocks());
738 
739  // FIXME: this a quite ugly as we need a temporary object:
741  min_level, max_level);
742 
743  for (unsigned int b = 0; b < n_blocks; ++b)
744  {
745  for (unsigned int l = min_level; l <= max_level; ++l)
746  {
747  src_non_block[l].reinit(src[l].block(b));
748  src_non_block[l] = src[l].block(b);
749  }
750  const unsigned int data_block = same_for_all ? 0 : b;
751  matrix_free_transfer_vector[data_block].copy_from_mg(*dof_handler[b],
752  dst.block(b),
753  src_non_block);
754  }
755 }
756 
757 
758 
759 #endif // DOXYGEN
760 
761 
763 
764 #endif
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
Definition: matrix_block.h:618
bool restriction_is_additive(const unsigned int index) const
Definition: fe.h:3241
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:1580
std::vector< unsigned int > n_owned_level_cells
cell_iterator begin(const unsigned int level=0) const
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
void reinit(const size_type size, const bool omit_zeroing_entries=false)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1533
void resize(const unsigned int new_minlevel, const unsigned int new_maxlevel, Args &&... args)
void interpolate_to_mg(const DoFHandler< dim, spacedim > &dof_handler, MGLevelObject< LinearAlgebra::distributed::Vector< Number >> &dst, const LinearAlgebra::distributed::Vector< Number2 > &src) const
const FiniteElement< dim, spacedim > & get_fe(const unsigned int index=0) const
void extract_locally_relevant_level_dofs(const DoFHandler< dim, spacedim > &dof_handler, const unsigned int level, IndexSet &dof_set)
Definition: dof_tools.cc:1180
size_type size() const
Definition: index_set.h:1632
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:1423
virtual const FullMatrix< double > & get_restriction_matrix(const unsigned int child, const RefinementCase< dim > &refinement_case=RefinementCase< dim >::isotropic_refinement) const
Definition: fe.cc:335
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:470
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:369
unsigned int level
Definition: grid_out.cc:4337
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
std::vector< std::vector< unsigned int > > level_dof_indices
virtual const MPI_Comm & get_communicator() const
Definition: tria_base.cc:138
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:607
std::vector< std::vector< std::pair< unsigned int, unsigned int > > > parent_child_connect
unsigned int n_dofs_per_cell() const
const Triangulation< dim, spacedim > & get_triangulation() const
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:368
std::vector< MGTransferMatrixFree< dim, Number > > matrix_free_transfer_vector
virtual ~MGTransferMatrixFree() override=default
const IndexSet & locally_owned_mg_dofs(const unsigned int level) const
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:1830
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)