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_smoother.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1999 - 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_smoother_h
17 #define dealii_mg_smoother_h
18 
19 
20 #include <deal.II/base/config.h>
21 
22 #include <deal.II/base/mg_level_object.h>
23 #include <deal.II/base/smartpointer.h>
24 
25 #include <deal.II/lac/linear_operator.h>
26 #include <deal.II/lac/vector_memory.h>
27 
28 #include <deal.II/multigrid/mg_base.h>
29 
30 #include <vector>
31 
32 DEAL_II_NAMESPACE_OPEN
33 
34 /*
35  * MGSmootherBase is defined in mg_base.h
36  */
37 
40 
49 template <typename VectorType>
50 class MGSmoother : public MGSmootherBase<VectorType>
51 {
52 public:
56  MGSmoother(const unsigned int steps = 1,
57  const bool variable = false,
58  const bool symmetric = false,
59  const bool transpose = false);
60 
64  void
65  set_steps(const unsigned int);
66 
70  void
71  set_variable(const bool);
72 
76  void
77  set_symmetric(const bool);
78 
83  void
84  set_transpose(const bool);
85 
90  void
91  set_debug(const unsigned int level);
92 
93 protected:
100  mutable GrowingVectorMemory<VectorType> vector_memory;
101 
106  unsigned int steps;
107 
112  bool variable;
113 
118  bool symmetric;
119 
124  bool transpose;
125 
129  unsigned int debug;
130 };
131 
132 
141 template <typename VectorType>
142 class MGSmootherIdentity : public MGSmootherBase<VectorType>
143 {
144 public:
150  virtual void
151  smooth(const unsigned int level, VectorType &u, const VectorType &rhs) const;
152 
153  virtual void
154  clear();
155 };
156 
157 
158 namespace mg
159 {
189  template <class RelaxationType, typename VectorType>
190  class SmootherRelaxation : public MGLevelObject<RelaxationType>,
191  public MGSmoother<VectorType>
192  {
193  public:
197  SmootherRelaxation(const unsigned int steps = 1,
198  const bool variable = false,
199  const bool symmetric = false,
200  const bool transpose = false);
201 
210  template <typename MatrixType2>
211  void
212  initialize(const MGLevelObject<MatrixType2> & matrices,
213  const typename RelaxationType::AdditionalData &additional_data =
214  typename RelaxationType::AdditionalData());
215 
223  template <typename MatrixType2, class DATA>
224  void
225  initialize(const MGLevelObject<MatrixType2> &matrices,
226  const MGLevelObject<DATA> & additional_data);
227 
231  void
232  clear() override;
233 
237  virtual void
238  smooth(const unsigned int level,
239  VectorType & u,
240  const VectorType & rhs) const override;
241 
258  virtual void
259  apply(const unsigned int level,
260  VectorType & u,
261  const VectorType & rhs) const override;
262 
266  std::size_t
267  memory_consumption() const;
268  };
269 } // namespace mg
270 
302 template <typename MatrixType, class RelaxationType, typename VectorType>
303 class MGSmootherRelaxation : public MGSmoother<VectorType>
304 {
305 public:
309  MGSmootherRelaxation(const unsigned int steps = 1,
310  const bool variable = false,
311  const bool symmetric = false,
312  const bool transpose = false);
313 
322  template <typename MatrixType2>
323  void
324  initialize(const MGLevelObject<MatrixType2> & matrices,
325  const typename RelaxationType::AdditionalData &additional_data =
326  typename RelaxationType::AdditionalData());
327 
336  template <typename MatrixType2, class DATA>
337  void
338  initialize(const MGLevelObject<MatrixType2> &matrices,
339  const MGLevelObject<DATA> & additional_data);
340 
350  template <typename MatrixType2, class DATA>
351  void
352  initialize(const MGLevelObject<MatrixType2> &matrices,
353  const DATA & additional_data,
354  const unsigned int block_row,
355  const unsigned int block_col);
356 
366  template <typename MatrixType2, class DATA>
367  void
368  initialize(const MGLevelObject<MatrixType2> &matrices,
369  const MGLevelObject<DATA> & additional_data,
370  const unsigned int block_row,
371  const unsigned int block_col);
372 
376  void
377  clear();
378 
382  virtual void
383  smooth(const unsigned int level, VectorType &u, const VectorType &rhs) const;
384 
400  virtual void
401  apply(const unsigned int level, VectorType &u, const VectorType &rhs) const;
402 
406  MGLevelObject<RelaxationType> smoothers;
407 
411  std::size_t
412  memory_consumption() const;
413 
414 
415 private:
419  MGLevelObject<LinearOperator<VectorType>> matrices;
420 };
421 
422 
423 
454 template <typename MatrixType, typename PreconditionerType, typename VectorType>
455 class MGSmootherPrecondition : public MGSmoother<VectorType>
456 {
457 public:
461  MGSmootherPrecondition(const unsigned int steps = 1,
462  const bool variable = false,
463  const bool symmetric = false,
464  const bool transpose = false);
465 
475  template <typename MatrixType2>
476  void
477  initialize(const MGLevelObject<MatrixType2> &matrices,
478  const typename PreconditionerType::AdditionalData &
479  additional_data = typename PreconditionerType::AdditionalData());
480 
490  template <typename MatrixType2, class DATA>
491  void
492  initialize(const MGLevelObject<MatrixType2> &matrices,
493  const MGLevelObject<DATA> & additional_data);
494 
505  template <typename MatrixType2, class DATA>
506  void
507  initialize(const MGLevelObject<MatrixType2> &matrices,
508  const DATA & additional_data,
509  const unsigned int block_row,
510  const unsigned int block_col);
511 
522  template <typename MatrixType2, class DATA>
523  void
524  initialize(const MGLevelObject<MatrixType2> &matrices,
525  const MGLevelObject<DATA> & additional_data,
526  const unsigned int block_row,
527  const unsigned int block_col);
528 
532  void
533  clear() override;
534 
538  virtual void
539  smooth(const unsigned int level,
540  VectorType & u,
541  const VectorType & rhs) const override;
542 
558  virtual void
559  apply(const unsigned int level,
560  VectorType & u,
561  const VectorType & rhs) const override;
562 
566  MGLevelObject<PreconditionerType> smoothers;
567 
571  std::size_t
572  memory_consumption() const;
573 
574 
575 private:
579  MGLevelObject<LinearOperator<VectorType>> matrices;
580 };
581 
584 /* ------------------------------- Inline functions --------------------------
585  */
586 
587 #ifndef DOXYGEN
588 
589 template <typename VectorType>
590 inline void
591 MGSmootherIdentity<VectorType>::smooth(const unsigned int,
592  VectorType &,
593  const VectorType &) const
594 {}
595 
596 template <typename VectorType>
597 inline void
598 MGSmootherIdentity<VectorType>::clear()
599 {}
600 
601 //---------------------------------------------------------------------------
602 
603 template <typename VectorType>
604 inline MGSmoother<VectorType>::MGSmoother(const unsigned int steps,
605  const bool variable,
606  const bool symmetric,
607  const bool transpose)
608  : steps(steps)
609  , variable(variable)
610  , symmetric(symmetric)
611  , transpose(transpose)
612  , debug(0)
613 {}
614 
615 
616 template <typename VectorType>
617 inline void
618 MGSmoother<VectorType>::set_steps(const unsigned int s)
619 {
620  steps = s;
621 }
622 
623 
624 template <typename VectorType>
625 inline void
626 MGSmoother<VectorType>::set_debug(const unsigned int s)
627 {
628  debug = s;
629 }
630 
631 
632 template <typename VectorType>
633 inline void
634 MGSmoother<VectorType>::set_variable(const bool flag)
635 {
636  variable = flag;
637 }
638 
639 
640 template <typename VectorType>
641 inline void
642 MGSmoother<VectorType>::set_symmetric(const bool flag)
643 {
644  symmetric = flag;
645 }
646 
647 
648 template <typename VectorType>
649 inline void
650 MGSmoother<VectorType>::set_transpose(const bool flag)
651 {
652  transpose = flag;
653 }
654 
655 //----------------------------------------------------------------------//
656 
657 namespace mg
658 {
659  template <class RelaxationType, typename VectorType>
660  inline SmootherRelaxation<RelaxationType, VectorType>::SmootherRelaxation(
661  const unsigned int steps,
662  const bool variable,
663  const bool symmetric,
664  const bool transpose)
665  : MGSmoother<VectorType>(steps, variable, symmetric, transpose)
666  {}
667 
668 
669  template <class RelaxationType, typename VectorType>
670  inline void
671  SmootherRelaxation<RelaxationType, VectorType>::clear()
672  {
673  MGLevelObject<RelaxationType>::clear_elements();
674  }
675 
676 
677  template <class RelaxationType, typename VectorType>
678  template <typename MatrixType2>
679  inline void
680  SmootherRelaxation<RelaxationType, VectorType>::initialize(
681  const MGLevelObject<MatrixType2> & m,
682  const typename RelaxationType::AdditionalData &data)
683  {
684  const unsigned int min = m.min_level();
685  const unsigned int max = m.max_level();
686 
687  this->resize(min, max);
688 
689  for (unsigned int i = min; i <= max; ++i)
690  (*this)[i].initialize(m[i], data);
691  }
692 
693 
694  template <class RelaxationType, typename VectorType>
695  template <typename MatrixType2, class DATA>
696  inline void
697  SmootherRelaxation<RelaxationType, VectorType>::initialize(
698  const MGLevelObject<MatrixType2> &m,
699  const MGLevelObject<DATA> & data)
700  {
701  const unsigned int min = std::max(m.min_level(), data.min_level());
702  const unsigned int max = std::min(m.max_level(), data.max_level());
703 
704  this->resize(min, max);
705 
706  for (unsigned int i = min; i <= max; ++i)
707  (*this)[i].initialize(m[i], data[i]);
708  }
709 
710 
711  template <class RelaxationType, typename VectorType>
712  inline void
713  SmootherRelaxation<RelaxationType, VectorType>::smooth(
714  const unsigned int level,
715  VectorType & u,
716  const VectorType & rhs) const
717  {
718  unsigned int maxlevel = this->max_level();
719  unsigned int steps2 = this->steps;
720 
721  if (this->variable)
722  steps2 *= (1 << (maxlevel - level));
723 
724  bool T = this->transpose;
725  if (this->symmetric && (steps2 % 2 == 0))
726  T = false;
727  if (this->debug > 0)
728  deallog << 'S' << level << ' ';
729 
730  for (unsigned int i = 0; i < steps2; ++i)
731  {
732  if (T)
733  (*this)[level].Tstep(u, rhs);
734  else
735  (*this)[level].step(u, rhs);
736  if (this->symmetric)
737  T = !T;
738  }
739  }
740 
741 
742  template <class RelaxationType, typename VectorType>
743  inline void
744  SmootherRelaxation<RelaxationType, VectorType>::apply(
745  const unsigned int level,
746  VectorType & u,
747  const VectorType & rhs) const
748  {
749  unsigned int maxlevel = this->max_level();
750  unsigned int steps2 = this->steps;
751 
752  if (this->variable)
753  steps2 *= (1 << (maxlevel - level));
754 
755  bool T = this->transpose;
756  if (this->symmetric && (steps2 % 2 == 0))
757  T = false;
758  if (this->debug > 0)
759  deallog << 'S' << level << ' ';
760 
761  if (T)
762  (*this)[level].Tvmult(u, rhs);
763  else
764  (*this)[level].vmult(u, rhs);
765  if (this->symmetric)
766  T = !T;
767  for (unsigned int i = 1; i < steps2; ++i)
768  {
769  if (T)
770  (*this)[level].Tstep(u, rhs);
771  else
772  (*this)[level].step(u, rhs);
773  if (this->symmetric)
774  T = !T;
775  }
776  }
777 
778 
779  template <class RelaxationType, typename VectorType>
780  inline std::size_t
781  SmootherRelaxation<RelaxationType, VectorType>::memory_consumption() const
782  {
783  return sizeof(*this) - sizeof(MGLevelObject<RelaxationType>) +
784  MGLevelObject<RelaxationType>::memory_consumption() +
785  this->vector_memory.memory_consumption();
786  }
787 } // namespace mg
788 
789 
790 //----------------------------------------------------------------------//
791 
792 template <typename MatrixType, class RelaxationType, typename VectorType>
793 inline MGSmootherRelaxation<MatrixType, RelaxationType, VectorType>::
794  MGSmootherRelaxation(const unsigned int steps,
795  const bool variable,
796  const bool symmetric,
797  const bool transpose)
798  : MGSmoother<VectorType>(steps, variable, symmetric, transpose)
799 {}
800 
801 
802 
803 template <typename MatrixType, class RelaxationType, typename VectorType>
804 inline void
805 MGSmootherRelaxation<MatrixType, RelaxationType, VectorType>::clear()
806 {
807  smoothers.clear_elements();
808 
809  unsigned int i = matrices.min_level(), max_level = matrices.max_level();
810  for (; i <= max_level; ++i)
811  matrices[i] = LinearOperator<VectorType>();
812 }
813 
814 
815 template <typename MatrixType, class RelaxationType, typename VectorType>
816 template <typename MatrixType2>
817 inline void
818 MGSmootherRelaxation<MatrixType, RelaxationType, VectorType>::initialize(
819  const MGLevelObject<MatrixType2> & m,
820  const typename RelaxationType::AdditionalData &data)
821 {
822  const unsigned int min = m.min_level();
823  const unsigned int max = m.max_level();
824 
825  matrices.resize(min, max);
826  smoothers.resize(min, max);
827 
828  for (unsigned int i = min; i <= max; ++i)
829  {
830  // Workaround: Unfortunately, not every "m[i]" object has a rich
831  // enough interface to populate reinit_(domain|range)_vector. Thus,
832  // apply an empty LinearOperator exemplar.
833  matrices[i] =
834  linear_operator<VectorType>(LinearOperator<VectorType>(), m[i]);
835  smoothers[i].initialize(m[i], data);
836  }
837 }
838 
839 template <typename MatrixType, class RelaxationType, typename VectorType>
840 template <typename MatrixType2, class DATA>
841 inline void
842 MGSmootherRelaxation<MatrixType, RelaxationType, VectorType>::initialize(
843  const MGLevelObject<MatrixType2> &m,
844  const MGLevelObject<DATA> & data)
845 {
846  const unsigned int min = m.min_level();
847  const unsigned int max = m.max_level();
848 
849  Assert(data.min_level() == min, ExcDimensionMismatch(data.min_level(), min));
850  Assert(data.max_level() == max, ExcDimensionMismatch(data.max_level(), max));
851 
852  matrices.resize(min, max);
853  smoothers.resize(min, max);
854 
855  for (unsigned int i = min; i <= max; ++i)
856  {
857  // Workaround: Unfortunately, not every "m[i]" object has a rich
858  // enough interface to populate reinit_(domain|range)_vector. Thus,
859  // apply an empty LinearOperator exemplar.
860  matrices[i] =
861  linear_operator<VectorType>(LinearOperator<VectorType>(), m[i]);
862  smoothers[i].initialize(m[i], data[i]);
863  }
864 }
865 
866 template <typename MatrixType, class RelaxationType, typename VectorType>
867 template <typename MatrixType2, class DATA>
868 inline void
869 MGSmootherRelaxation<MatrixType, RelaxationType, VectorType>::initialize(
870  const MGLevelObject<MatrixType2> &m,
871  const DATA & data,
872  const unsigned int row,
873  const unsigned int col)
874 {
875  const unsigned int min = m.min_level();
876  const unsigned int max = m.max_level();
877 
878  matrices.resize(min, max);
879  smoothers.resize(min, max);
880 
881  for (unsigned int i = min; i <= max; ++i)
882  {
883  // Workaround: Unfortunately, not every "m[i]" object has a rich
884  // enough interface to populate reinit_(domain|range)_vector. Thus,
885  // apply an empty LinearOperator exemplar.
886  matrices[i] = linear_operator<VectorType>(LinearOperator<VectorType>(),
887  m[i].block(row, col));
888  smoothers[i].initialize(m[i].block(row, col), data);
889  }
890 }
891 
892 template <typename MatrixType, class RelaxationType, typename VectorType>
893 template <typename MatrixType2, class DATA>
894 inline void
895 MGSmootherRelaxation<MatrixType, RelaxationType, VectorType>::initialize(
896  const MGLevelObject<MatrixType2> &m,
897  const MGLevelObject<DATA> & data,
898  const unsigned int row,
899  const unsigned int col)
900 {
901  const unsigned int min = m.min_level();
902  const unsigned int max = m.max_level();
903 
904  Assert(data.min_level() == min, ExcDimensionMismatch(data.min_level(), min));
905  Assert(data.max_level() == max, ExcDimensionMismatch(data.max_level(), max));
906 
907  matrices.resize(min, max);
908  smoothers.resize(min, max);
909 
910  for (unsigned int i = min; i <= max; ++i)
911  {
912  // Workaround: Unfortunately, not every "m[i]" object has a rich
913  // enough interface to populate reinit_(domain|range)_vector. Thus,
914  // apply an empty LinearOperator exemplar.
915  matrices[i] = linear_operator<VectorType>(LinearOperator<VectorType>(),
916  m[i].block(row, col));
917  smoothers[i].initialize(m[i].block(row, col), data[i]);
918  }
919 }
920 
921 
922 template <typename MatrixType, class RelaxationType, typename VectorType>
923 inline void
924 MGSmootherRelaxation<MatrixType, RelaxationType, VectorType>::smooth(
925  const unsigned int level,
926  VectorType & u,
927  const VectorType & rhs) const
928 {
929  unsigned int maxlevel = smoothers.max_level();
930  unsigned int steps2 = this->steps;
931 
932  if (this->variable)
933  steps2 *= (1 << (maxlevel - level));
934 
935  bool T = this->transpose;
936  if (this->symmetric && (steps2 % 2 == 0))
937  T = false;
938  if (this->debug > 0)
939  deallog << 'S' << level << ' ';
940 
941  for (unsigned int i = 0; i < steps2; ++i)
942  {
943  if (T)
944  smoothers[level].Tstep(u, rhs);
945  else
946  smoothers[level].step(u, rhs);
947  if (this->symmetric)
948  T = !T;
949  }
950 }
951 
952 
953 template <typename MatrixType, class RelaxationType, typename VectorType>
954 inline void
955 MGSmootherRelaxation<MatrixType, RelaxationType, VectorType>::apply(
956  const unsigned int level,
957  VectorType & u,
958  const VectorType & rhs) const
959 {
960  unsigned int maxlevel = smoothers.max_level();
961  unsigned int steps2 = this->steps;
962 
963  if (this->variable)
964  steps2 *= (1 << (maxlevel - level));
965 
966  bool T = this->transpose;
967  if (this->symmetric && (steps2 % 2 == 0))
968  T = false;
969  if (this->debug > 0)
970  deallog << 'S' << level << ' ';
971 
972  if (T)
973  smoothers[level].Tvmult(u, rhs);
974  else
975  smoothers[level].vmult(u, rhs);
976  if (this->symmetric)
977  T = !T;
978  for (unsigned int i = 1; i < steps2; ++i)
979  {
980  if (T)
981  smoothers[level].Tstep(u, rhs);
982  else
983  smoothers[level].step(u, rhs);
984  if (this->symmetric)
985  T = !T;
986  }
987 }
988 
989 
990 
991 template <typename MatrixType, class RelaxationType, typename VectorType>
992 inline std::size_t
993 MGSmootherRelaxation<MatrixType, RelaxationType, VectorType>::
994  memory_consumption() const
995 {
996  return sizeof(*this) + matrices.memory_consumption() +
997  smoothers.memory_consumption() +
998  this->vector_memory.memory_consumption();
999 }
1000 
1001 
1002 //----------------------------------------------------------------------//
1003 
1004 template <typename MatrixType, typename PreconditionerType, typename VectorType>
1005 inline MGSmootherPrecondition<MatrixType, PreconditionerType, VectorType>::
1006  MGSmootherPrecondition(const unsigned int steps,
1007  const bool variable,
1008  const bool symmetric,
1009  const bool transpose)
1010  : MGSmoother<VectorType>(steps, variable, symmetric, transpose)
1011 {}
1012 
1013 
1014 
1015 template <typename MatrixType, typename PreconditionerType, typename VectorType>
1016 inline void
1017 MGSmootherPrecondition<MatrixType, PreconditionerType, VectorType>::clear()
1018 {
1019  smoothers.clear_elements();
1020 
1021  unsigned int i = matrices.min_level(), max_level = matrices.max_level();
1022  for (; i <= max_level; ++i)
1023  matrices[i] = LinearOperator<VectorType>();
1024 }
1025 
1026 
1027 
1028 template <typename MatrixType, typename PreconditionerType, typename VectorType>
1029 template <typename MatrixType2>
1030 inline void
1031 MGSmootherPrecondition<MatrixType, PreconditionerType, VectorType>::initialize(
1032  const MGLevelObject<MatrixType2> & m,
1033  const typename PreconditionerType::AdditionalData &data)
1034 {
1035  const unsigned int min = m.min_level();
1036  const unsigned int max = m.max_level();
1037 
1038  matrices.resize(min, max);
1039  smoothers.resize(min, max);
1040 
1041  for (unsigned int i = min; i <= max; ++i)
1042  {
1043  // Workaround: Unfortunately, not every "m[i]" object has a rich
1044  // enough interface to populate reinit_(domain|range)_vector. Thus,
1045  // apply an empty LinearOperator exemplar.
1046  matrices[i] =
1047  linear_operator<VectorType>(LinearOperator<VectorType>(), m[i]);
1048  smoothers[i].initialize(m[i], data);
1049  }
1050 }
1051 
1052 
1053 
1054 template <typename MatrixType, typename PreconditionerType, typename VectorType>
1055 template <typename MatrixType2, class DATA>
1056 inline void
1057 MGSmootherPrecondition<MatrixType, PreconditionerType, VectorType>::initialize(
1058  const MGLevelObject<MatrixType2> &m,
1059  const MGLevelObject<DATA> & data)
1060 {
1061  const unsigned int min = m.min_level();
1062  const unsigned int max = m.max_level();
1063 
1064  Assert(data.min_level() == min, ExcDimensionMismatch(data.min_level(), min));
1065  Assert(data.max_level() == max, ExcDimensionMismatch(data.max_level(), max));
1066 
1067  matrices.resize(min, max);
1068  smoothers.resize(min, max);
1069 
1070  for (unsigned int i = min; i <= max; ++i)
1071  {
1072  // Workaround: Unfortunately, not every "m[i]" object has a rich
1073  // enough interface to populate reinit_(domain|range)_vector. Thus,
1074  // apply an empty LinearOperator exemplar.
1075  matrices[i] =
1076  linear_operator<VectorType>(LinearOperator<VectorType>(), m[i]);
1077  smoothers[i].initialize(m[i], data[i]);
1078  }
1079 }
1080 
1081 
1082 
1083 template <typename MatrixType, typename PreconditionerType, typename VectorType>
1084 template <typename MatrixType2, class DATA>
1085 inline void
1086 MGSmootherPrecondition<MatrixType, PreconditionerType, VectorType>::initialize(
1087  const MGLevelObject<MatrixType2> &m,
1088  const DATA & data,
1089  const unsigned int row,
1090  const unsigned int col)
1091 {
1092  const unsigned int min = m.min_level();
1093  const unsigned int max = m.max_level();
1094 
1095  matrices.resize(min, max);
1096  smoothers.resize(min, max);
1097 
1098  for (unsigned int i = min; i <= max; ++i)
1099  {
1100  matrices[i] = &(m[i].block(row, col));
1101  smoothers[i].initialize(m[i].block(row, col), data);
1102  }
1103 }
1104 
1105 
1106 
1107 template <typename MatrixType, typename PreconditionerType, typename VectorType>
1108 template <typename MatrixType2, class DATA>
1109 inline void
1110 MGSmootherPrecondition<MatrixType, PreconditionerType, VectorType>::initialize(
1111  const MGLevelObject<MatrixType2> &m,
1112  const MGLevelObject<DATA> & data,
1113  const unsigned int row,
1114  const unsigned int col)
1115 {
1116  const unsigned int min = m.min_level();
1117  const unsigned int max = m.max_level();
1118 
1119  Assert(data.min_level() == min, ExcDimensionMismatch(data.min_level(), min));
1120  Assert(data.max_level() == max, ExcDimensionMismatch(data.max_level(), max));
1121 
1122  matrices.resize(min, max);
1123  smoothers.resize(min, max);
1124 
1125  for (unsigned int i = min; i <= max; ++i)
1126  {
1127  matrices[i] = &(m[i].block(row, col));
1128  smoothers[i].initialize(m[i].block(row, col), data[i]);
1129  }
1130 }
1131 
1132 
1133 
1134 template <typename MatrixType, typename PreconditionerType, typename VectorType>
1135 inline void
1136 MGSmootherPrecondition<MatrixType, PreconditionerType, VectorType>::smooth(
1137  const unsigned int level,
1138  VectorType & u,
1139  const VectorType & rhs) const
1140 {
1141  unsigned int maxlevel = matrices.max_level();
1142  unsigned int steps2 = this->steps;
1143 
1144  if (this->variable)
1145  steps2 *= (1 << (maxlevel - level));
1146 
1147  typename VectorMemory<VectorType>::Pointer r(this->vector_memory);
1148  typename VectorMemory<VectorType>::Pointer d(this->vector_memory);
1149 
1150  r->reinit(u, true);
1151  d->reinit(u, true);
1152 
1153  bool T = this->transpose;
1154  if (this->symmetric && (steps2 % 2 == 0))
1155  T = false;
1156  if (this->debug > 0)
1157  deallog << 'S' << level << ' ';
1158 
1159  for (unsigned int i = 0; i < steps2; ++i)
1160  {
1161  if (T)
1162  {
1163  if (this->debug > 0)
1164  deallog << 'T';
1165  matrices[level].Tvmult(*r, u);
1166  r->sadd(-1., 1., rhs);
1167  if (this->debug > 2)
1168  deallog << ' ' << r->l2_norm() << ' ';
1169  smoothers[level].Tvmult(*d, *r);
1170  if (this->debug > 1)
1171  deallog << ' ' << d->l2_norm() << ' ';
1172  }
1173  else
1174  {
1175  if (this->debug > 0)
1176  deallog << 'N';
1177  matrices[level].vmult(*r, u);
1178  r->sadd(-1., rhs);
1179  if (this->debug > 2)
1180  deallog << ' ' << r->l2_norm() << ' ';
1181  smoothers[level].vmult(*d, *r);
1182  if (this->debug > 1)
1183  deallog << ' ' << d->l2_norm() << ' ';
1184  }
1185  u += *d;
1186  if (this->symmetric)
1187  T = !T;
1188  }
1189  if (this->debug > 0)
1190  deallog << std::endl;
1191 }
1192 
1193 
1194 
1195 template <typename MatrixType, typename PreconditionerType, typename VectorType>
1196 inline void
1197 MGSmootherPrecondition<MatrixType, PreconditionerType, VectorType>::apply(
1198  const unsigned int level,
1199  VectorType & u,
1200  const VectorType & rhs) const
1201 {
1202  unsigned int maxlevel = matrices.max_level();
1203  unsigned int steps2 = this->steps;
1204 
1205  if (this->variable)
1206  steps2 *= (1 << (maxlevel - level));
1207 
1208  bool T = this->transpose;
1209  if (this->symmetric && (steps2 % 2 == 0))
1210  T = false;
1211  if (this->debug > 0)
1212  deallog << 'S' << level << ' ';
1213 
1214  // first step where we overwrite the result
1215  if (this->debug > 2)
1216  deallog << ' ' << rhs.l2_norm() << ' ';
1217  if (this->debug > 0)
1218  deallog << (T ? 'T' : 'N');
1219  if (T)
1220  smoothers[level].Tvmult(u, rhs);
1221  else
1222  smoothers[level].vmult(u, rhs);
1223  if (this->debug > 1)
1224  deallog << ' ' << u.l2_norm() << ' ';
1225  if (this->symmetric)
1226  T = !T;
1227 
1228  typename VectorMemory<VectorType>::Pointer r(this->vector_memory);
1229  typename VectorMemory<VectorType>::Pointer d(this->vector_memory);
1230 
1231  if (steps2 > 1)
1232  {
1233  r->reinit(u, true);
1234  d->reinit(u, true);
1235  }
1236 
1237  for (unsigned int i = 1; i < steps2; ++i)
1238  {
1239  if (T)
1240  {
1241  if (this->debug > 0)
1242  deallog << 'T';
1243  matrices[level].Tvmult(*r, u);
1244  r->sadd(-1., 1., rhs);
1245  if (this->debug > 2)
1246  deallog << ' ' << r->l2_norm() << ' ';
1247  smoothers[level].Tvmult(*d, *r);
1248  if (this->debug > 1)
1249  deallog << ' ' << d->l2_norm() << ' ';
1250  }
1251  else
1252  {
1253  if (this->debug > 0)
1254  deallog << 'N';
1255  matrices[level].vmult(*r, u);
1256  r->sadd(-1., rhs);
1257  if (this->debug > 2)
1258  deallog << ' ' << r->l2_norm() << ' ';
1259  smoothers[level].vmult(*d, *r);
1260  if (this->debug > 1)
1261  deallog << ' ' << d->l2_norm() << ' ';
1262  }
1263  u += *d;
1264  if (this->symmetric)
1265  T = !T;
1266  }
1267  if (this->debug > 0)
1268  deallog << std::endl;
1269 }
1270 
1271 
1272 
1273 template <typename MatrixType, typename PreconditionerType, typename VectorType>
1274 inline std::size_t
1275 MGSmootherPrecondition<MatrixType, PreconditionerType, VectorType>::
1276  memory_consumption() const
1277 {
1278  return sizeof(*this) + matrices.memory_consumption() +
1279  smoothers.memory_consumption() +
1280  this->vector_memory.memory_consumption();
1281 }
1282 
1283 
1284 #endif // DOXYGEN
1285 
1286 DEAL_II_NAMESPACE_CLOSE
1287 
1288 #endif
MGSmootherPrecondition::smooth
virtual void smooth(const unsigned int level, VectorType &u, const VectorType &rhs) const override
mg_base.h
MGSmootherPrecondition::MGSmootherPrecondition
MGSmootherPrecondition(const unsigned int steps=1, const bool variable=false, const bool symmetric=false, const bool transpose=false)
MGSmoother
Definition: mg_smoother.h:50
MGLevelObject::max_level
unsigned int max_level() const
Definition: mg_level_object.h:246
mg::SmootherRelaxation::clear
void clear() override
MGLevelObject::min_level
unsigned int min_level() const
Definition: mg_level_object.h:238
VectorizedArray::max
VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &x, const ::VectorizedArray< Number, width > &y)
Definition: vectorization.h:5466
MGSmootherRelaxation::clear
void clear()
VectorizedArray::min
VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &x, const ::VectorizedArray< Number, width > &y)
Definition: vectorization.h:5483
linear_operator.h
VectorType
MGSmoother::symmetric
bool symmetric
Definition: mg_smoother.h:118
MGSmootherRelaxation::memory_consumption
std::size_t memory_consumption() const
mg::SmootherRelaxation
Definition: mg_smoother.h:190
mg
Definition: mg.h:81
MGSmootherPrecondition::matrices
MGLevelObject< LinearOperator< VectorType > > matrices
Definition: mg_smoother.h:579
MGSmootherRelaxation::apply
virtual void apply(const unsigned int level, VectorType &u, const VectorType &rhs) const
Physics::Elasticity::Kinematics::d
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
deallog
LogStream deallog
Definition: logstream.cc:37
MGSmoother::set_variable
void set_variable(const bool)
MGSmoother::set_debug
void set_debug(const unsigned int level)
MGSmootherRelaxation::matrices
MGLevelObject< LinearOperator< VectorType > > matrices
Definition: mg_smoother.h:419
GrowingVectorMemory
Definition: vector_memory.h:320
MGSmoother::vector_memory
GrowingVectorMemory< VectorType > vector_memory
Definition: mg_smoother.h:100
level
unsigned int level
Definition: grid_out.cc:4355
MGSmoother::set_symmetric
void set_symmetric(const bool)
MGSmootherIdentity::clear
virtual void clear()
MGSmootherPrecondition::initialize
void initialize(const MGLevelObject< MatrixType2 > &matrices, const typename PreconditionerType::AdditionalData &additional_data=typename PreconditionerType::AdditionalData())
LAPACKSupport::T
static const char T
Definition: lapack_support.h:163
mg::SmootherRelaxation::initialize
void initialize(const MGLevelObject< MatrixType2 > &matrices, const typename RelaxationType::AdditionalData &additional_data=typename RelaxationType::AdditionalData())
LAPACKSupport::symmetric
@ symmetric
Matrix is symmetric.
Definition: lapack_support.h:115
MGSmootherIdentity::smooth
virtual void smooth(const unsigned int level, VectorType &u, const VectorType &rhs) const
MGSmootherRelaxation
Definition: mg_smoother.h:303
VectorMemory::Pointer
Definition: vector_memory.h:192
MGSmoother::set_steps
void set_steps(const unsigned int)
DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:358
transpose
DerivativeForm< 1, spacedim, dim, Number > transpose(const DerivativeForm< 1, dim, spacedim, Number > &DF)
Definition: derivative_form.h:470
MGSmootherPrecondition::memory_consumption
std::size_t memory_consumption() const
MemoryConsumption::memory_consumption
std::enable_if< std::is_fundamental< T >::value, std::size_t >::type memory_consumption(const T &t)
Definition: memory_consumption.h:268
MGSmoother::debug
unsigned int debug
Definition: mg_smoother.h:129
MGLevelObject::clear_elements
void clear_elements()
Definition: mg_level_object.h:228
std_cxx17::apply
auto apply(F &&fn, Tuple &&t) -> decltype(apply_impl(std::forward< F >(fn), std::forward< Tuple >(t), std_cxx14::make_index_sequence< std::tuple_size< typename std::remove_reference< Tuple >::type >::value >()))
Definition: tuple.h:40
MGSmootherPrecondition::smoothers
MGLevelObject< PreconditionerType > smoothers
Definition: mg_smoother.h:566
MGSmootherBase
Definition: mg_base.h:245
smartpointer.h
mg_level_object.h
MGSmoother::transpose
bool transpose
Definition: mg_smoother.h:124
MGSmootherRelaxation::initialize
void initialize(const MGLevelObject< MatrixType2 > &matrices, const typename RelaxationType::AdditionalData &additional_data=typename RelaxationType::AdditionalData())
MGSmootherRelaxation::smoothers
MGLevelObject< RelaxationType > smoothers
Definition: mg_smoother.h:406
MGSmootherPrecondition
Definition: mg_smoother.h:455
MGSmootherRelaxation::smooth
virtual void smooth(const unsigned int level, VectorType &u, const VectorType &rhs) const
MGSmoother::variable
bool variable
Definition: mg_smoother.h:112
Assert
#define Assert(cond, exc)
Definition: exceptions.h:1419
MGSmootherPrecondition::clear
void clear() override
MGSmoother::set_transpose
void set_transpose(const bool)
MGLevelObject
Definition: mg_level_object.h:50
mg::SmootherRelaxation::memory_consumption
std::size_t memory_consumption() const
MGLevelObject::memory_consumption
std::size_t memory_consumption() const
Definition: mg_level_object.h:265
mg::SmootherRelaxation::apply
virtual void apply(const unsigned int level, VectorType &u, const VectorType &rhs) const override
Utilities::MPI::min
T min(const T &t, const MPI_Comm &mpi_communicator)
StandardExceptions::ExcDimensionMismatch
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
MGSmootherPrecondition::apply
virtual void apply(const unsigned int level, VectorType &u, const VectorType &rhs) const override
vector_memory.h
MGSmoother::MGSmoother
MGSmoother(const unsigned int steps=1, const bool variable=false, const bool symmetric=false, const bool transpose=false)
config.h
MGSmootherRelaxation::MGSmootherRelaxation
MGSmootherRelaxation(const unsigned int steps=1, const bool variable=false, const bool symmetric=false, const bool transpose=false)
mg::SmootherRelaxation::SmootherRelaxation
SmootherRelaxation(const unsigned int steps=1, const bool variable=false, const bool symmetric=false, const bool transpose=false)
DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:359
mg::SmootherRelaxation::smooth
virtual void smooth(const unsigned int level, VectorType &u, const VectorType &rhs) const override
LinearOperator< VectorType >
DerivativeForm::transpose
DerivativeForm< 1, spacedim, dim, Number > transpose(const DerivativeForm< 1, dim, spacedim, Number > &DF)
Definition: derivative_form.h:470
Utilities::MPI::max
T max(const T &t, const MPI_Comm &mpi_communicator)
MGSmoother::steps
unsigned int steps
Definition: mg_smoother.h:106
MGSmootherIdentity
Definition: mg_smoother.h:142
Generated by   doxygen 1.8.17