Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
mg_transfer_matrix_free.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2016 - 2018 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 
21 #include <deal.II/base/mg_level_object.h>
22 #include <deal.II/base/vectorization.h>
23 
24 #include <deal.II/dofs/dof_handler.h>
25 
26 #include <deal.II/lac/la_parallel_block_vector.h>
27 #include <deal.II/lac/la_parallel_vector.h>
28 
29 #include <deal.II/multigrid/mg_base.h>
30 #include <deal.II/multigrid/mg_constrained_dofs.h>
31 #include <deal.II/multigrid/mg_transfer.h>
32 #include <deal.II/multigrid/mg_transfer_internal.h>
33 
34 
35 DEAL_II_NAMESPACE_OPEN
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
83 
87  void
88  clear();
89 
93  void
94  build(const DoFHandler<dim, dim> &mg_dof);
95 
110  virtual void
111  prolongate(
112  const unsigned int to_level,
114  const LinearAlgebra::distributed::Vector<Number> &src) const override;
115 
134  virtual void
136  const unsigned int from_level,
138  const LinearAlgebra::distributed::Vector<Number> &src) const override;
139 
147  template <typename Number2, int spacedim>
148  void
150  const DoFHandler<dim, spacedim> & mg_dof,
153 
158 
162  std::size_t
163  memory_consumption() const;
164 
165 private:
171  unsigned int fe_degree;
172 
178 
183  unsigned int n_components;
184 
190  unsigned int n_child_cell_dofs;
191 
202  std::vector<std::vector<unsigned int>> level_dof_indices;
203 
208  std::vector<std::vector<std::pair<unsigned int, unsigned int>>>
210 
215  std::vector<unsigned int> n_owned_level_cells;
216 
222 
227 
239  std::vector<AlignedVector<VectorizedArray<Number>>> weights_on_refined;
240 
246  std::vector<std::vector<std::vector<unsigned short>>> dirichlet_indices;
247 
251  template <int degree>
252  void
254  const unsigned int to_level,
257 
261  template <int degree>
262  void
263  do_restrict_add(const unsigned int from_level,
266 };
267 
268 
285 template <int dim, typename Number>
287  : public MGTransferBase<LinearAlgebra::distributed::BlockVector<Number>>
288 {
289 public:
294  MGTransferBlockMatrixFree() = default;
295 
300  MGTransferBlockMatrixFree(const MGConstrainedDoFs &mg_constrained_dofs);
301 
306  const std::vector<MGConstrainedDoFs> &mg_constrained_dofs);
307 
311  virtual ~MGTransferBlockMatrixFree() override = default;
312 
316  void
317  initialize_constraints(const MGConstrainedDoFs &mg_constrained_dofs);
318 
322  void
324  const std::vector<MGConstrainedDoFs> &mg_constrained_dofs);
325 
329  void
330  clear();
331 
335  void
336  build(const DoFHandler<dim, dim> &mg_dof);
337 
341  void
342  build(const std::vector<const DoFHandler<dim, dim> *> &mg_dof);
343 
358  virtual void
359  prolongate(
360  const unsigned int to_level,
362  const LinearAlgebra::distributed::BlockVector<Number> &src) const override;
363 
382  virtual void
384  const unsigned int from_level,
386  const LinearAlgebra::distributed::BlockVector<Number> &src) const override;
387 
398  template <typename Number2, int spacedim>
399  void
400  copy_to_mg(
401  const DoFHandler<dim, spacedim> & mg_dof,
404 
408  template <typename Number2, int spacedim>
409  void
410  copy_to_mg(
411  const std::vector<const DoFHandler<dim, spacedim> *> & mg_dof,
414 
418  template <typename Number2, int spacedim>
419  void
420  copy_from_mg(
421  const DoFHandler<dim, spacedim> & mg_dof,
424  const;
425 
429  template <typename Number2, int spacedim>
430  void
431  copy_from_mg(
432  const std::vector<const DoFHandler<dim, spacedim> *> &mg_dof,
435  const;
436 
440  std::size_t
441  memory_consumption() const;
442 
447  static const bool supports_dof_handler_vector = true;
448 
449 private:
453  std::vector<MGTransferMatrixFree<dim, Number>> matrix_free_transfer_vector;
454 
459  const bool same_for_all;
460 };
461 
462 
466 //------------------------ templated functions -------------------------
467 #ifndef DOXYGEN
468 
469 
470 template <int dim, typename Number>
471 template <typename Number2, int spacedim>
472 void
474  const DoFHandler<dim, spacedim> & dof_handler,
477 {
478  const unsigned int min_level = dst.min_level();
479  const unsigned int max_level = dst.max_level();
480 
481  Assert(max_level == dof_handler.get_triangulation().n_global_levels() - 1,
483  max_level, dof_handler.get_triangulation().n_global_levels() - 1));
484 
486  (dynamic_cast<const parallel::Triangulation<dim, spacedim> *>(
487  &dof_handler.get_triangulation()));
488  MPI_Comm mpi_communicator =
489  p_tria != nullptr ? p_tria->get_communicator() : MPI_COMM_SELF;
490 
491  // resize the dst vector if it's empty or has incorrect size
492  MGLevelObject<IndexSet> relevant_dofs(min_level, max_level);
493  for (unsigned int level = min_level; level <= max_level; ++level)
494  {
496  level,
497  relevant_dofs[level]);
498  if (dst[level].size() !=
499  dof_handler.locally_owned_mg_dofs(level).size() ||
500  dst[level].local_size() !=
501  dof_handler.locally_owned_mg_dofs(level).n_elements())
502  dst[level].reinit(dof_handler.locally_owned_mg_dofs(level),
503  relevant_dofs[level],
504  mpi_communicator);
505  }
506 
507  const FiniteElement<dim, spacedim> &fe = dof_handler.get_fe();
508 
509  // copy fine level vector to active cells in MG hierarchy
510  this->copy_to_mg(dof_handler, dst, src, true);
511 
512  // FIXME: maybe need to store hanging nodes constraints per level?
513  // MGConstrainedDoFs does NOT keep this info right now, only periodicity
514  // constraints...
515  dst[max_level].update_ghost_values();
516  // do the transfer from level to level-1:
517  for (unsigned int level = max_level; level > min_level; --level)
518  {
519  // auxiliary vector which always has ghost elements
521  dof_handler.locally_owned_mg_dofs(level),
522  relevant_dofs[level],
523  mpi_communicator);
524  ghosted_vector = dst[level];
525  ghosted_vector.update_ghost_values();
526 
527  std::vector<Number> dof_values_coarse(fe.dofs_per_cell);
528  Vector<Number> dof_values_fine(fe.dofs_per_cell);
529  Vector<Number> tmp(fe.dofs_per_cell);
530  std::vector<types::global_dof_index> dof_indices(fe.dofs_per_cell);
531  typename DoFHandler<dim>::cell_iterator cell =
532  dof_handler.begin(level - 1);
533  typename DoFHandler<dim>::cell_iterator endc = dof_handler.end(level - 1);
534  for (; cell != endc; ++cell)
535  if (cell->is_locally_owned_on_level())
536  {
537  // if we get to a cell without children (== active), we can
538  // skip it as there values should be already set by the
539  // equivalent of copy_to_mg()
540  if (!cell->has_children())
541  continue;
542 
543  std::fill(dof_values_coarse.begin(), dof_values_coarse.end(), 0.);
544  for (unsigned int child = 0; child < cell->n_children(); ++child)
545  {
546  cell->child(child)->get_mg_dof_indices(dof_indices);
547  for (unsigned int i = 0; i < fe.dofs_per_cell; ++i)
548  dof_values_fine(i) = ghosted_vector(dof_indices[i]);
549  fe.get_restriction_matrix(child, cell->refinement_case())
550  .vmult(tmp, dof_values_fine);
551  for (unsigned int i = 0; i < fe.dofs_per_cell; ++i)
552  if (fe.restriction_is_additive(i))
553  dof_values_coarse[i] += tmp[i];
554  else if (tmp(i) != 0.)
555  dof_values_coarse[i] = tmp[i];
556  }
557  cell->get_mg_dof_indices(dof_indices);
558  for (unsigned int i = 0; i < fe.dofs_per_cell; ++i)
559  dst[level - 1](dof_indices[i]) = dof_values_coarse[i];
560  }
561 
562  dst[level - 1].compress(VectorOperation::insert);
563  dst[level - 1].update_ghost_values();
564  }
565 }
566 
567 
568 
569 template <int dim, typename Number>
570 template <typename Number2, int spacedim>
571 void
573  const DoFHandler<dim, spacedim> & mg_dof,
576 {
577  AssertDimension(matrix_free_transfer_vector.size(), 1);
578  Assert(same_for_all,
579  ExcMessage(
580  "This object was initialized with support for usage with one "
581  "DoFHandler for each block, but this method assumes that "
582  "the same DoFHandler is used for all the blocks!"));
583  const std::vector<const DoFHandler<dim, spacedim> *> mg_dofs(src.n_blocks(),
584  &mg_dof);
585 
586  copy_to_mg(mg_dofs, dst, src);
587 }
588 
589 template <int dim, typename Number>
590 template <typename Number2, int spacedim>
591 void
593  const std::vector<const DoFHandler<dim, spacedim> *> & mg_dof,
596 {
597  const unsigned int n_blocks = src.n_blocks();
598  AssertDimension(mg_dof.size(), n_blocks);
599 
600  if (n_blocks == 0)
601  return;
602 
603  const unsigned int min_level = dst.min_level();
604  const unsigned int max_level = dst.max_level();
605 
606  // this function is normally called within the Multigrid class with
607  // dst == defect level block vector. At first run this vector is not
608  // initialized. Do this below:
609  {
611  (dynamic_cast<const parallel::Triangulation<dim, spacedim> *>(
612  &(mg_dof[0]->get_triangulation())));
613  for (unsigned int i = 1; i < n_blocks; ++i)
615  &(mg_dof[0]->get_triangulation())) == tria),
616  ExcMessage("The DoFHandler use different Triangulations!"));
617 
618  MGLevelObject<bool> do_reinit;
619  do_reinit.resize(min_level, max_level);
620  for (unsigned int level = min_level; level <= max_level; ++level)
621  {
622  do_reinit[level] = false;
623  if (dst[level].n_blocks() != n_blocks)
624  {
625  do_reinit[level] = true;
626  continue; // level
627  }
628  for (unsigned int b = 0; b < n_blocks; ++b)
629  {
630  LinearAlgebra::distributed::Vector<Number> &v = dst[level].block(b);
631  if (v.size() != mg_dof[b]->locally_owned_mg_dofs(level).size() ||
632  v.local_size() !=
633  mg_dof[b]->locally_owned_mg_dofs(level).n_elements())
634  {
635  do_reinit[level] = true;
636  break; // b
637  }
638  }
639  }
640 
641  for (unsigned int level = min_level; level <= max_level; ++level)
642  {
643  if (do_reinit[level])
644  {
645  dst[level].reinit(n_blocks);
646  for (unsigned int b = 0; b < n_blocks; ++b)
647  {
649  dst[level].block(b);
650  v.reinit(mg_dof[b]->locally_owned_mg_dofs(level),
651  tria != nullptr ? tria->get_communicator() :
652  MPI_COMM_SELF);
653  }
654  dst[level].collect_sizes();
655  }
656  else
657  dst[level] = 0;
658  }
659  }
660 
661  // FIXME: this a quite ugly as we need a temporary object:
663  min_level, max_level);
664 
665  for (unsigned int b = 0; b < n_blocks; ++b)
666  {
667  for (unsigned int l = min_level; l <= max_level; ++l)
668  dst_non_block[l].reinit(dst[l].block(b));
669  const unsigned int data_block = same_for_all ? 0 : b;
670  matrix_free_transfer_vector[data_block].copy_to_mg(*mg_dof[b],
671  dst_non_block,
672  src.block(b));
673 
674  for (unsigned int l = min_level; l <= max_level; ++l)
675  dst[l].block(b) = dst_non_block[l];
676  }
677 }
678 
679 template <int dim, typename Number>
680 template <typename Number2, int spacedim>
681 void
683  const DoFHandler<dim, spacedim> & mg_dof,
686  const
687 {
688  AssertDimension(matrix_free_transfer_vector.size(), 1);
689  const std::vector<const DoFHandler<dim, spacedim> *> mg_dofs(dst.n_blocks(),
690  &mg_dof);
691 
692  copy_from_mg(mg_dofs, dst, src);
693 }
694 
695 template <int dim, typename Number>
696 template <typename Number2, int spacedim>
697 void
699  const std::vector<const DoFHandler<dim, spacedim> *> & mg_dof,
702  const
703 {
704  const unsigned int n_blocks = dst.n_blocks();
705  AssertDimension(mg_dof.size(), n_blocks);
706 
707  if (n_blocks == 0)
708  return;
709 
710  const unsigned int min_level = src.min_level();
711  const unsigned int max_level = src.max_level();
712 
713  for (unsigned int l = min_level; l <= max_level; ++l)
714  AssertDimension(src[l].n_blocks(), dst.n_blocks());
715 
716  // FIXME: this a quite ugly as we need a temporary object:
718  min_level, max_level);
719 
720  for (unsigned int b = 0; b < n_blocks; ++b)
721  {
722  for (unsigned int l = min_level; l <= max_level; ++l)
723  {
724  src_non_block[l].reinit(src[l].block(b));
725  src_non_block[l] = src[l].block(b);
726  }
727  const unsigned int data_block = same_for_all ? 0 : b;
728  matrix_free_transfer_vector[data_block].copy_from_mg(*mg_dof[b],
729  dst.block(b),
730  src_non_block);
731  }
732 }
733 
734 
735 
736 #endif // DOXYGEN
737 
738 
739 DEAL_II_NAMESPACE_CLOSE
740 
741 #endif
const Triangulation< dim, spacedim > & get_triangulation() const
bool restriction_is_additive(const unsigned int index) const
Definition: fe.h:3276
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:1567
std::vector< unsigned int > n_owned_level_cells
cell_iterator begin(const unsigned int level=0) const
Definition: dof_handler.cc:930
AlignedVector< VectorizedArray< Number > > evaluation_data
static ::ExceptionBase & ExcNoProlongation()
std::size_t memory_consumption() const
std::vector< std::vector< std::vector< unsigned short > > > dirichlet_indices
void build(const DoFHandler< dim, dim > &mg_dof)
void initialize_constraints(const MGConstrainedDoFs &mg_constrained_dofs)
cell_iterator end() const
Definition: dof_handler.cc:959
void reinit(const size_type size, const bool omit_zeroing_entries=false)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1519
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:1179
std::size_t memory_consumption() const
size_type size() const
Definition: index_set.h:1600
void interpolate_to_mg(const DoFHandler< dim, spacedim > &mg_dof, MGLevelObject< LinearAlgebra::distributed::Vector< Number >> &dst, const LinearAlgebra::distributed::Vector< Number2 > &src) const
static ::ExceptionBase & ExcMessage(std::string arg1)
AlignedVector< VectorizedArray< Number > > prolongation_matrix_1d
#define Assert(cond, exc)
Definition: exceptions.h:1407
virtual const FullMatrix< double > & get_restriction_matrix(const unsigned int child, const RefinementCase< dim > &refinement_case=RefinementCase< dim >::isotropic_refinement) const
Definition: fe.cc:306
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 restrict_and_add(const unsigned int from_level, LinearAlgebra::distributed::BlockVector< Number > &dst, const LinearAlgebra::distributed::BlockVector< 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
virtual MPI_Comm get_communicator() const
Definition: tria_base.cc:155
void copy_to_mg(const DoFHandler< dim, spacedim > &mg_dof, MGLevelObject< LinearAlgebra::distributed::BlockVector< Number >> &dst, const LinearAlgebra::distributed::BlockVector< Number2 > &src) const
#define DeclException0(Exception0)
Definition: exceptions.h:473
std::vector< AlignedVector< VectorizedArray< Number > > > weights_on_refined
const IndexSet & locally_owned_mg_dofs(const unsigned int level) const
std::vector< std::vector< unsigned int > > level_dof_indices
virtual ~MGTransferBlockMatrixFree() override=default
const unsigned int dofs_per_cell
Definition: fe_base.h:282
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
static const bool supports_dof_handler_vector
void copy_from_mg(const DoFHandler< dim, spacedim > &mg_dof, 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:616
std::vector< std::vector< std::pair< unsigned int, unsigned int > > > parent_child_connect
virtual void prolongate(const unsigned int to_level, LinearAlgebra::distributed::BlockVector< Number > &dst, const LinearAlgebra::distributed::BlockVector< Number > &src) const override
std::vector< MGTransferMatrixFree< dim, Number > > matrix_free_transfer_vector
MGTransferBlockMatrixFree()=default
virtual ~MGTransferMatrixFree() override=default
void resize(const unsigned int new_minlevel, const unsigned int new_maxlevel)
void initialize_constraints(const MGConstrainedDoFs &mg_constrained_dofs)
BlockType & block(const unsigned int i)
size_type n_elements() const
Definition: index_set.h:1732
void build(const DoFHandler< dim, dim > &mg_dof)
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)