Reference documentation for deal.II version Git 3f1f337db3 2021-10-23 13:19:02 -0600
\(\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 - 2021 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 
128  virtual void
130  const unsigned int to_level,
132  const LinearAlgebra::distributed::Vector<Number> &src) const override;
133 
152  virtual void
154  const unsigned int from_level,
156  const LinearAlgebra::distributed::Vector<Number> &src) const override;
157 
172  template <typename Number2, int spacedim>
173  void
175  const DoFHandler<dim, spacedim> & dof_handler,
178 
183 
187  std::size_t
188  memory_consumption() const;
189 
190 private:
196  unsigned int fe_degree;
197 
203 
208  unsigned int n_components;
209 
215  unsigned int n_child_cell_dofs;
216 
227  std::vector<std::vector<unsigned int>> level_dof_indices;
228 
233  std::vector<std::vector<std::pair<unsigned int, unsigned int>>>
235 
240  std::vector<unsigned int> n_owned_level_cells;
241 
247 
252 
264  std::vector<AlignedVector<VectorizedArray<Number>>> weights_on_refined;
265 
271  std::vector<std::vector<std::vector<unsigned short>>> dirichlet_indices;
272 
281 
285  template <int degree>
286  void
288  const unsigned int to_level,
291 
295  template <int degree>
296  void
297  do_restrict_add(const unsigned int from_level,
300 };
301 
302 
316 template <int dim, typename Number>
318  : public MGTransferBase<LinearAlgebra::distributed::BlockVector<Number>>
319 {
320 public:
325  MGTransferBlockMatrixFree() = default;
326 
332 
337  const std::vector<MGConstrainedDoFs> &mg_constrained_dofs);
338 
342  virtual ~MGTransferBlockMatrixFree() override = default;
343 
347  void
348  initialize_constraints(const MGConstrainedDoFs &mg_constrained_dofs);
349 
353  void
355  const std::vector<MGConstrainedDoFs> &mg_constrained_dofs);
356 
360  void
361  clear();
362 
366  void
367  build(const DoFHandler<dim, dim> &dof_handler);
368 
372  void
373  build(const std::vector<const DoFHandler<dim, dim> *> &dof_handler);
374 
389  virtual void
390  prolongate(
391  const unsigned int to_level,
393  const LinearAlgebra::distributed::BlockVector<Number> &src) const override;
394 
395  virtual void
397  const unsigned int to_level,
399  const LinearAlgebra::distributed::BlockVector<Number> &src) const override;
400 
419  virtual void
421  const unsigned int from_level,
423  const LinearAlgebra::distributed::BlockVector<Number> &src) const override;
424 
435  template <typename Number2, int spacedim>
436  void
437  copy_to_mg(
438  const DoFHandler<dim, spacedim> & dof_handler,
441 
445  template <typename Number2, int spacedim>
446  void
447  copy_to_mg(
448  const std::vector<const DoFHandler<dim, spacedim> *> & dof_handler,
451 
455  template <typename Number2, int spacedim>
456  void
457  copy_from_mg(
458  const DoFHandler<dim, spacedim> & dof_handler,
461  const;
462 
466  template <typename Number2, int spacedim>
467  void
468  copy_from_mg(
469  const std::vector<const DoFHandler<dim, spacedim> *> &dof_handler,
472  const;
473 
477  std::size_t
478  memory_consumption() const;
479 
484  static const bool supports_dof_handler_vector = true;
485 
486 private:
490  std::vector<MGTransferMatrixFree<dim, Number>> matrix_free_transfer_vector;
491 
496  const bool same_for_all;
497 };
498 
499 
503 //------------------------ templated functions -------------------------
504 #ifndef DOXYGEN
505 
506 
507 template <int dim, typename Number>
508 template <typename Number2, int spacedim>
509 void
511  const DoFHandler<dim, spacedim> & dof_handler,
514 {
515  const unsigned int min_level = dst.min_level();
516  const unsigned int max_level = dst.max_level();
517 
518  Assert(max_level == dof_handler.get_triangulation().n_global_levels() - 1,
520  max_level, dof_handler.get_triangulation().n_global_levels() - 1));
521 
522  const FiniteElement<dim, spacedim> &fe = dof_handler.get_fe();
523 
524  for (unsigned int level = min_level; level <= max_level; ++level)
525  if (dst[level].size() != dof_handler.n_dofs(level) ||
526  dst[level].locally_owned_size() !=
527  dof_handler.locally_owned_mg_dofs(level).n_elements())
528  dst[level].reinit(this->vector_partitioners[level]);
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 
537  // do the transfer from level to level-1:
538  dst[max_level].update_ghost_values();
539 
540  for (unsigned int level = max_level; level > min_level; --level)
541  {
542  // auxiliary vector which always has ghost elements
543  const LinearAlgebra::distributed::Vector<Number> *input = nullptr;
545  if (dst[level].get_partitioner().get() ==
546  this->vector_partitioners[level].get())
547  input = &dst[level];
548  else
549  {
550  ghosted_fine.reinit(this->vector_partitioners[level]);
551  ghosted_fine.copy_locally_owned_data_from(dst[level]);
552  ghosted_fine.update_ghost_values();
553  input = &ghosted_fine;
554  }
555 
556  std::vector<Number> dof_values_coarse(fe.n_dofs_per_cell());
557  Vector<Number> dof_values_fine(fe.n_dofs_per_cell());
558  Vector<Number> tmp(fe.n_dofs_per_cell());
559  std::vector<types::global_dof_index> dof_indices(fe.n_dofs_per_cell());
560  for (const auto &cell : dof_handler.cell_iterators_on_level(level - 1))
561  if (cell->is_locally_owned_on_level())
562  {
563  // if we get to a cell without children (== active), we can
564  // skip it as there values should be already set by the
565  // equivalent of copy_to_mg()
566  if (cell->is_active())
567  continue;
568 
569  std::fill(dof_values_coarse.begin(), dof_values_coarse.end(), 0.);
570  for (unsigned int child = 0; child < cell->n_children(); ++child)
571  {
572  cell->child(child)->get_mg_dof_indices(dof_indices);
573 
575  this->mg_constrained_dofs, level, dof_indices);
576 
577  for (unsigned int i = 0; i < fe.n_dofs_per_cell(); ++i)
578  dof_values_fine(i) = (*input)(dof_indices[i]);
579  fe.get_restriction_matrix(child, cell->refinement_case())
580  .vmult(tmp, dof_values_fine);
581  for (unsigned int i = 0; i < fe.n_dofs_per_cell(); ++i)
582  if (fe.restriction_is_additive(i))
583  dof_values_coarse[i] += tmp[i];
584  else if (tmp(i) != 0.)
585  dof_values_coarse[i] = tmp[i];
586  }
587  cell->get_mg_dof_indices(dof_indices);
588  for (unsigned int i = 0; i < fe.n_dofs_per_cell(); ++i)
589  if (dof_handler.locally_owned_mg_dofs(level - 1).is_element(
590  dof_indices[i]))
591  dst[level - 1](dof_indices[i]) = dof_values_coarse[i];
592  }
593 
594  dst[level - 1].update_ghost_values();
595  }
596 }
597 
598 
599 
600 template <int dim, typename Number>
601 template <typename Number2, int spacedim>
602 void
604  const DoFHandler<dim, spacedim> & dof_handler,
607 {
608  AssertDimension(matrix_free_transfer_vector.size(), 1);
609  Assert(same_for_all,
610  ExcMessage(
611  "This object was initialized with support for usage with one "
612  "DoFHandler for each block, but this method assumes that "
613  "the same DoFHandler is used for all the blocks!"));
614  const std::vector<const DoFHandler<dim, spacedim> *> mg_dofs(src.n_blocks(),
615  &dof_handler);
616 
617  copy_to_mg(mg_dofs, dst, src);
618 }
619 
620 
621 
622 template <int dim, typename Number>
623 template <typename Number2, int spacedim>
624 void
626  const std::vector<const DoFHandler<dim, spacedim> *> & dof_handler,
629 {
630  const unsigned int n_blocks = src.n_blocks();
631  AssertDimension(dof_handler.size(), n_blocks);
632 
633  if (n_blocks == 0)
634  return;
635 
636  const unsigned int min_level = dst.min_level();
637  const unsigned int max_level = dst.max_level();
638 
639  // this function is normally called within the Multigrid class with
640  // dst == defect level block vector. At first run this vector is not
641  // initialized. Do this below:
642  {
644  (dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
645  &(dof_handler[0]->get_triangulation())));
646  for (unsigned int i = 1; i < n_blocks; ++i)
647  AssertThrow(
648  (dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
649  &(dof_handler[0]->get_triangulation())) == tria),
650  ExcMessage("The DoFHandler use different Triangulations!"));
651 
652  MGLevelObject<bool> do_reinit;
653  do_reinit.resize(min_level, max_level);
654  for (unsigned int level = min_level; level <= max_level; ++level)
655  {
656  do_reinit[level] = false;
657  if (dst[level].n_blocks() != n_blocks)
658  {
659  do_reinit[level] = true;
660  continue; // level
661  }
662  for (unsigned int b = 0; b < n_blocks; ++b)
663  {
665  if (v.size() !=
666  dof_handler[b]->locally_owned_mg_dofs(level).size() ||
667  v.locally_owned_size() !=
668  dof_handler[b]->locally_owned_mg_dofs(level).n_elements())
669  {
670  do_reinit[level] = true;
671  break; // b
672  }
673  }
674  }
675 
676  for (unsigned int level = min_level; level <= max_level; ++level)
677  {
678  if (do_reinit[level])
679  {
680  dst[level].reinit(n_blocks);
681  for (unsigned int b = 0; b < n_blocks; ++b)
682  {
684  dst[level].block(b);
685  v.reinit(dof_handler[b]->locally_owned_mg_dofs(level),
686  dof_handler[b]->get_communicator());
687  }
688  dst[level].collect_sizes();
689  }
690  else
691  dst[level] = 0;
692  }
693  }
694 
695  // FIXME: this a quite ugly as we need a temporary object:
697  min_level, max_level);
698 
699  for (unsigned int b = 0; b < n_blocks; ++b)
700  {
701  for (unsigned int l = min_level; l <= max_level; ++l)
702  dst_non_block[l].reinit(dst[l].block(b));
703  const unsigned int data_block = same_for_all ? 0 : b;
704  matrix_free_transfer_vector[data_block].copy_to_mg(*dof_handler[b],
705  dst_non_block,
706  src.block(b));
707 
708  for (unsigned int l = min_level; l <= max_level; ++l)
709  dst[l].block(b) = dst_non_block[l];
710  }
711 }
712 
713 template <int dim, typename Number>
714 template <typename Number2, int spacedim>
715 void
717  const DoFHandler<dim, spacedim> & dof_handler,
720  const
721 {
722  AssertDimension(matrix_free_transfer_vector.size(), 1);
723  const std::vector<const DoFHandler<dim, spacedim> *> mg_dofs(dst.n_blocks(),
724  &dof_handler);
725 
726  copy_from_mg(mg_dofs, dst, src);
727 }
728 
729 template <int dim, typename Number>
730 template <typename Number2, int spacedim>
731 void
733  const std::vector<const DoFHandler<dim, spacedim> *> &dof_handler,
736  const
737 {
738  const unsigned int n_blocks = dst.n_blocks();
739  AssertDimension(dof_handler.size(), n_blocks);
740 
741  if (n_blocks == 0)
742  return;
743 
744  const unsigned int min_level = src.min_level();
745  const unsigned int max_level = src.max_level();
746 
747  for (unsigned int l = min_level; l <= max_level; ++l)
748  AssertDimension(src[l].n_blocks(), dst.n_blocks());
749 
750  // FIXME: this a quite ugly as we need a temporary object:
752  min_level, max_level);
753 
754  for (unsigned int b = 0; b < n_blocks; ++b)
755  {
756  for (unsigned int l = min_level; l <= max_level; ++l)
757  {
758  src_non_block[l].reinit(src[l].block(b));
759  src_non_block[l] = src[l].block(b);
760  }
761  const unsigned int data_block = same_for_all ? 0 : b;
762  matrix_free_transfer_vector[data_block].copy_from_mg(*dof_handler[b],
763  dst.block(b),
764  src_non_block);
765  }
766 }
767 
768 
769 
770 #endif // DOXYGEN
771 
772 
774 
775 #endif
const Triangulation< dim, spacedim > & get_triangulation() const
virtual void prolongate_and_add(const unsigned int to_level, LinearAlgebra::distributed::Vector< Number > &dst, const LinearAlgebra::distributed::Vector< Number > &src) const override
void copy_locally_owned_data_from(const Vector< Number2, MemorySpace > &src)
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:3271
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:1655
std::vector< unsigned int > n_owned_level_cells
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
const ::Triangulation< dim, spacedim > & tria
std::vector< std::vector< std::vector< unsigned short > > > dirichlet_indices
void reinit(const size_type size, const bool omit_zeroing_entries=false)
void resize(const unsigned int new_minlevel, const unsigned int new_maxlevel, Args &&...args)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1571
const FiniteElement< dim, spacedim > & get_fe(const unsigned int index=0) const
size_type locally_owned_size() const
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:1639
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:1461
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:464
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:401
unsigned int level
Definition: grid_out.cc:4589
std::vector< AlignedVector< VectorizedArray< Number > > > weights_on_refined
types::global_dof_index n_dofs() const
void resolve_identity_constraints(const MGConstrainedDoFs *mg_constrained_dofs, const unsigned int level, std::vector< types::global_dof_index > &dof_indices)
std::enable_if< IsBlockVector< VectorType >::value, unsigned int >::type n_blocks(const VectorType &vector)
Definition: operators.h:50
const IndexSet & locally_owned_mg_dofs(const unsigned int level) const
std::vector< std::vector< unsigned int > > level_dof_indices
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:586
std::vector< std::vector< std::pair< unsigned int, unsigned int > > > parent_child_connect
unsigned int n_dofs_per_cell() const
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:400
std::vector< MGTransferMatrixFree< dim, Number > > matrix_free_transfer_vector
virtual ~MGTransferMatrixFree() override=default
IteratorRange< cell_iterator > cell_iterators_on_level(const unsigned int level) const
bool is_element(const size_type index) const
Definition: index_set.h:1770
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:1837
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)