Reference documentation for deal.II version GIT ca9e7e8105 2023-03-24 14:15:03+00:00
\(\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 - 2022 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 
24 
27 
29 
30 #include <vector>
31 
33 
34 /*
35  * MGSmootherBase is defined in mg_base.h
36  */
37 
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:
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 
139 template <typename VectorType>
140 class MGSmootherIdentity : public MGSmootherBase<VectorType>
141 {
142 public:
148  virtual void
149  smooth(const unsigned int level, VectorType &u, const VectorType &rhs) const;
150 
151  virtual void
152  clear();
153 };
154 
155 
156 namespace mg
157 {
184  template <class RelaxationType, typename VectorType>
185  class SmootherRelaxation : public MGLevelObject<RelaxationType>,
186  public MGSmoother<VectorType>
187  {
188  public:
192  SmootherRelaxation(const unsigned int steps = 1,
193  const bool variable = false,
194  const bool symmetric = false,
195  const bool transpose = false);
196 
205  template <typename MatrixType2>
206  void
208  const typename RelaxationType::AdditionalData &additional_data =
209  typename RelaxationType::AdditionalData());
210 
218  template <typename MatrixType2, class DATA>
219  void
221  const MGLevelObject<DATA> & additional_data);
222 
226  void
227  clear() override;
228 
232  virtual void
233  smooth(const unsigned int level,
234  VectorType & u,
235  const VectorType & rhs) const override;
236 
253  virtual void
254  apply(const unsigned int level,
255  VectorType & u,
256  const VectorType & rhs) const override;
257 
261  std::size_t
263  };
264 } // namespace mg
265 
295 template <typename MatrixType, class RelaxationType, typename VectorType>
296 class MGSmootherRelaxation : public MGSmoother<VectorType>
297 {
298 public:
302  MGSmootherRelaxation(const unsigned int steps = 1,
303  const bool variable = false,
304  const bool symmetric = false,
305  const bool transpose = false);
306 
315  template <typename MatrixType2>
316  void
318  const typename RelaxationType::AdditionalData &additional_data =
319  typename RelaxationType::AdditionalData());
320 
329  template <typename MatrixType2, class DATA>
330  void
332  const MGLevelObject<DATA> & additional_data);
333 
343  template <typename MatrixType2, class DATA>
344  void
346  const DATA & additional_data,
347  const unsigned int block_row,
348  const unsigned int block_col);
349 
359  template <typename MatrixType2, class DATA>
360  void
362  const MGLevelObject<DATA> & additional_data,
363  const unsigned int block_row,
364  const unsigned int block_col);
365 
369  void
370  clear();
371 
375  virtual void
376  smooth(const unsigned int level, VectorType &u, const VectorType &rhs) const;
377 
393  virtual void
394  apply(const unsigned int level, VectorType &u, const VectorType &rhs) const;
395 
400 
404  std::size_t
406 
407 
408 private:
413 };
414 
415 
416 
445 template <typename MatrixType, typename PreconditionerType, typename VectorType>
446 class MGSmootherPrecondition : public MGSmoother<VectorType>
447 {
448 public:
452  MGSmootherPrecondition(const unsigned int steps = 1,
453  const bool variable = false,
454  const bool symmetric = false,
455  const bool transpose = false);
456 
466  template <typename MatrixType2>
467  void
469  const typename PreconditionerType::AdditionalData &
470  additional_data = typename PreconditionerType::AdditionalData());
471 
479  template <typename MatrixType2>
480  void
482 
492  template <typename MatrixType2, class DATA>
493  void
495  const MGLevelObject<DATA> & additional_data);
496 
507  template <typename MatrixType2, class DATA>
508  void
510  const DATA & additional_data,
511  const unsigned int block_row,
512  const unsigned int block_col);
513 
524  template <typename MatrixType2, class DATA>
525  void
527  const MGLevelObject<DATA> & additional_data,
528  const unsigned int block_row,
529  const unsigned int block_col);
530 
534  void
535  clear() override;
536 
540  virtual void
541  smooth(const unsigned int level,
542  VectorType & u,
543  const VectorType & rhs) const override;
544 
560  virtual void
561  apply(const unsigned int level,
562  VectorType & u,
563  const VectorType & rhs) const override;
564 
569 
573  std::size_t
575 
576 
577 private:
582 };
583 
586 /* ------------------------------- Inline functions --------------------------
587  */
588 
589 #ifndef DOXYGEN
590 
591 template <typename VectorType>
592 inline void
593 MGSmootherIdentity<VectorType>::smooth(const unsigned int,
594  VectorType &,
595  const VectorType &) const
596 {}
597 
598 template <typename VectorType>
599 inline void
601 {}
602 
603 //---------------------------------------------------------------------------
604 
605 template <typename VectorType>
606 inline MGSmoother<VectorType>::MGSmoother(const unsigned int steps,
607  const bool variable,
608  const bool symmetric,
609  const bool transpose)
610  : steps(steps)
611  , variable(variable)
614  , debug(0)
615 {}
616 
617 
618 template <typename VectorType>
619 inline void
620 MGSmoother<VectorType>::set_steps(const unsigned int s)
621 {
622  steps = s;
623 }
624 
625 
626 template <typename VectorType>
627 inline void
628 MGSmoother<VectorType>::set_debug(const unsigned int s)
629 {
630  debug = s;
631 }
632 
633 
634 template <typename VectorType>
635 inline void
637 {
638  variable = flag;
639 }
640 
641 
642 template <typename VectorType>
643 inline void
645 {
646  symmetric = flag;
647 }
648 
649 
650 template <typename VectorType>
651 inline void
653 {
654  transpose = flag;
655 }
656 
657 //----------------------------------------------------------------------//
658 
659 namespace mg
660 {
661  template <class RelaxationType, typename VectorType>
663  const unsigned int steps,
664  const bool variable,
665  const bool symmetric,
666  const bool transpose)
667  : MGSmoother<VectorType>(steps, variable, symmetric, transpose)
668  {}
669 
670 
671  template <class RelaxationType, typename VectorType>
672  inline void
673  SmootherRelaxation<RelaxationType, VectorType>::clear()
674  {
676  }
677 
678 
679  template <class RelaxationType, typename VectorType>
680  template <typename MatrixType2>
681  inline void
682  SmootherRelaxation<RelaxationType, VectorType>::initialize(
683  const MGLevelObject<MatrixType2> & m,
684  const typename RelaxationType::AdditionalData &data)
685  {
686  const unsigned int min = m.min_level();
687  const unsigned int max = m.max_level();
688 
689  this->resize(min, max);
690 
691  for (unsigned int i = min; i <= max; ++i)
692  (*this)[i].initialize(m[i], data);
693  }
694 
695 
696  template <class RelaxationType, typename VectorType>
697  template <typename MatrixType2, class DATA>
698  inline void
699  SmootherRelaxation<RelaxationType, VectorType>::initialize(
701  const MGLevelObject<DATA> & data)
702  {
703  const unsigned int min = std::max(m.min_level(), data.min_level());
704  const unsigned int max = std::min(m.max_level(), data.max_level());
705 
706  this->resize(min, max);
707 
708  for (unsigned int i = min; i <= max; ++i)
709  (*this)[i].initialize(m[i], data[i]);
710  }
711 
712 
713  template <class RelaxationType, typename VectorType>
714  inline void
715  SmootherRelaxation<RelaxationType, VectorType>::smooth(
716  const unsigned int level,
717  VectorType & u,
718  const VectorType & rhs) const
719  {
720  unsigned int maxlevel = this->max_level();
721  unsigned int steps2 = this->steps;
722 
723  if (this->variable)
724  steps2 *= (1 << (maxlevel - level));
725 
726  bool T = this->transpose;
727  if (this->symmetric && (steps2 % 2 == 0))
728  T = false;
729  if (this->debug > 0)
730  deallog << 'S' << level << ' ';
731 
732  for (unsigned int i = 0; i < steps2; ++i)
733  {
734  if (T)
735  (*this)[level].Tstep(u, rhs);
736  else
737  (*this)[level].step(u, rhs);
738  if (this->symmetric)
739  T = !T;
740  }
741  }
742 
743 
744  template <class RelaxationType, typename VectorType>
745  inline void
747  const unsigned int level,
748  VectorType & u,
749  const VectorType & rhs) const
750  {
751  unsigned int maxlevel = this->max_level();
752  unsigned int steps2 = this->steps;
753 
754  if (this->variable)
755  steps2 *= (1 << (maxlevel - level));
756 
757  bool T = this->transpose;
758  if (this->symmetric && (steps2 % 2 == 0))
759  T = false;
760  if (this->debug > 0)
761  deallog << 'S' << level << ' ';
762 
763  if (T)
764  (*this)[level].Tvmult(u, rhs);
765  else
766  (*this)[level].vmult(u, rhs);
767  if (this->symmetric)
768  T = !T;
769  for (unsigned int i = 1; i < steps2; ++i)
770  {
771  if (T)
772  (*this)[level].Tstep(u, rhs);
773  else
774  (*this)[level].step(u, rhs);
775  if (this->symmetric)
776  T = !T;
777  }
778  }
779 
780 
781  template <class RelaxationType, typename VectorType>
782  inline std::size_t
784  {
785  return sizeof(*this) - sizeof(MGLevelObject<RelaxationType>) +
787  this->vector_memory.memory_consumption();
788  }
789 } // namespace mg
790 
791 
792 //----------------------------------------------------------------------//
793 
794 template <typename MatrixType, class RelaxationType, typename VectorType>
796  MGSmootherRelaxation(const unsigned int steps,
797  const bool variable,
798  const bool symmetric,
799  const bool transpose)
800  : MGSmoother<VectorType>(steps, variable, symmetric, transpose)
801 {}
802 
803 
804 
805 template <typename MatrixType, class RelaxationType, typename VectorType>
806 inline void
808 {
809  smoothers.clear_elements();
810 
811  unsigned int i = matrices.min_level(), max_level = matrices.max_level();
812  for (; i <= max_level; ++i)
813  matrices[i] = LinearOperator<VectorType>();
814 }
815 
816 
817 template <typename MatrixType, class RelaxationType, typename VectorType>
818 template <typename MatrixType2>
819 inline void
821  const MGLevelObject<MatrixType2> & m,
822  const typename RelaxationType::AdditionalData &data)
823 {
824  const unsigned int min = m.min_level();
825  const unsigned int max = m.max_level();
826 
827  matrices.resize(min, max);
828  smoothers.resize(min, max);
829 
830  for (unsigned int i = min; i <= max; ++i)
831  {
832  // Workaround: Unfortunately, not every "m[i]" object has a rich
833  // enough interface to populate reinit_(domain|range)_vector. Thus,
834  // apply an empty LinearOperator exemplar.
835  matrices[i] =
836  linear_operator<VectorType>(LinearOperator<VectorType>(), m[i]);
837  smoothers[i].initialize(m[i], data);
838  }
839 }
840 
841 template <typename MatrixType, class RelaxationType, typename VectorType>
842 template <typename MatrixType2, class DATA>
843 inline void
846  const MGLevelObject<DATA> & data)
847 {
848  const unsigned int min = m.min_level();
849  const unsigned int max = m.max_level();
850 
851  Assert(data.min_level() == min, ExcDimensionMismatch(data.min_level(), min));
852  Assert(data.max_level() == max, ExcDimensionMismatch(data.max_level(), max));
853 
854  matrices.resize(min, max);
855  smoothers.resize(min, max);
856 
857  for (unsigned int i = min; i <= max; ++i)
858  {
859  // Workaround: Unfortunately, not every "m[i]" object has a rich
860  // enough interface to populate reinit_(domain|range)_vector. Thus,
861  // apply an empty LinearOperator exemplar.
862  matrices[i] =
863  linear_operator<VectorType>(LinearOperator<VectorType>(), m[i]);
864  smoothers[i].initialize(m[i], data[i]);
865  }
866 }
867 
868 template <typename MatrixType, class RelaxationType, typename VectorType>
869 template <typename MatrixType2, class DATA>
870 inline void
873  const DATA & data,
874  const unsigned int row,
875  const unsigned int col)
876 {
877  const unsigned int min = m.min_level();
878  const unsigned int max = m.max_level();
879 
880  matrices.resize(min, max);
881  smoothers.resize(min, max);
882 
883  for (unsigned int i = min; i <= max; ++i)
884  {
885  // Workaround: Unfortunately, not every "m[i]" object has a rich
886  // enough interface to populate reinit_(domain|range)_vector. Thus,
887  // apply an empty LinearOperator exemplar.
888  matrices[i] = linear_operator<VectorType>(LinearOperator<VectorType>(),
889  m[i].block(row, col));
890  smoothers[i].initialize(m[i].block(row, col), data);
891  }
892 }
893 
894 template <typename MatrixType, class RelaxationType, typename VectorType>
895 template <typename MatrixType2, class DATA>
896 inline void
899  const MGLevelObject<DATA> & data,
900  const unsigned int row,
901  const unsigned int col)
902 {
903  const unsigned int min = m.min_level();
904  const unsigned int max = m.max_level();
905 
906  Assert(data.min_level() == min, ExcDimensionMismatch(data.min_level(), min));
907  Assert(data.max_level() == max, ExcDimensionMismatch(data.max_level(), max));
908 
909  matrices.resize(min, max);
910  smoothers.resize(min, max);
911 
912  for (unsigned int i = min; i <= max; ++i)
913  {
914  // Workaround: Unfortunately, not every "m[i]" object has a rich
915  // enough interface to populate reinit_(domain|range)_vector. Thus,
916  // apply an empty LinearOperator exemplar.
917  matrices[i] = linear_operator<VectorType>(LinearOperator<VectorType>(),
918  m[i].block(row, col));
919  smoothers[i].initialize(m[i].block(row, col), data[i]);
920  }
921 }
922 
923 
924 template <typename MatrixType, class RelaxationType, typename VectorType>
925 inline void
927  const unsigned int level,
928  VectorType & u,
929  const VectorType & rhs) const
930 {
931  unsigned int maxlevel = smoothers.max_level();
932  unsigned int steps2 = this->steps;
933 
934  if (this->variable)
935  steps2 *= (1 << (maxlevel - level));
936 
937  bool T = this->transpose;
938  if (this->symmetric && (steps2 % 2 == 0))
939  T = false;
940  if (this->debug > 0)
941  deallog << 'S' << level << ' ';
942 
943  for (unsigned int i = 0; i < steps2; ++i)
944  {
945  if (T)
946  smoothers[level].Tstep(u, rhs);
947  else
948  smoothers[level].step(u, rhs);
949  if (this->symmetric)
950  T = !T;
951  }
952 }
953 
954 
955 template <typename MatrixType, class RelaxationType, typename VectorType>
956 inline void
958  const unsigned int level,
959  VectorType & u,
960  const VectorType & rhs) const
961 {
962  unsigned int maxlevel = smoothers.max_level();
963  unsigned int steps2 = this->steps;
964 
965  if (this->variable)
966  steps2 *= (1 << (maxlevel - level));
967 
968  bool T = this->transpose;
969  if (this->symmetric && (steps2 % 2 == 0))
970  T = false;
971  if (this->debug > 0)
972  deallog << 'S' << level << ' ';
973 
974  if (T)
975  smoothers[level].Tvmult(u, rhs);
976  else
977  smoothers[level].vmult(u, rhs);
978  if (this->symmetric)
979  T = !T;
980  for (unsigned int i = 1; i < steps2; ++i)
981  {
982  if (T)
983  smoothers[level].Tstep(u, rhs);
984  else
985  smoothers[level].step(u, rhs);
986  if (this->symmetric)
987  T = !T;
988  }
989 }
990 
991 
992 
993 template <typename MatrixType, class RelaxationType, typename VectorType>
994 inline std::size_t
996  memory_consumption() const
997 {
998  return sizeof(*this) + matrices.memory_consumption() +
999  smoothers.memory_consumption() +
1000  this->vector_memory.memory_consumption();
1001 }
1002 
1003 
1004 //----------------------------------------------------------------------//
1005 
1006 template <typename MatrixType, typename PreconditionerType, typename VectorType>
1008  MGSmootherPrecondition(const unsigned int steps,
1009  const bool variable,
1010  const bool symmetric,
1011  const bool transpose)
1012  : MGSmoother<VectorType>(steps, variable, symmetric, transpose)
1013 {}
1014 
1015 
1016 
1017 template <typename MatrixType, typename PreconditionerType, typename VectorType>
1018 inline void
1020 {
1021  smoothers.clear_elements();
1022 
1023  unsigned int i = matrices.min_level(), max_level = matrices.max_level();
1024  for (; i <= max_level; ++i)
1025  matrices[i] = LinearOperator<VectorType>();
1026 }
1027 
1028 
1029 
1030 template <typename MatrixType, typename PreconditionerType, typename VectorType>
1031 template <typename MatrixType2>
1032 inline void
1034  const MGLevelObject<MatrixType2> & m,
1035  const typename PreconditionerType::AdditionalData &data)
1036 {
1037  const unsigned int min = m.min_level();
1038  const unsigned int max = m.max_level();
1039 
1040  matrices.resize(min, max);
1041  smoothers.resize(min, max);
1042 
1043  for (unsigned int i = min; i <= max; ++i)
1044  {
1045  // Workaround: Unfortunately, not every "m[i]" object has a rich
1046  // enough interface to populate reinit_(domain|range)_vector. Thus,
1047  // apply an empty LinearOperator exemplar.
1048  matrices[i] =
1049  linear_operator<VectorType>(LinearOperator<VectorType>(),
1051  smoothers[i].initialize(Utilities::get_underlying_value(m[i]), data);
1052  }
1053 }
1054 
1055 
1056 
1057 template <typename MatrixType, typename PreconditionerType, typename VectorType>
1058 template <typename MatrixType2>
1059 inline void
1062 {
1063  const unsigned int min = m.min_level();
1064  const unsigned int max = m.max_level();
1065 
1066  matrices.resize(min, max);
1067  smoothers.resize(min, max);
1068 
1069  for (unsigned int i = min; i <= max; ++i)
1070  {
1071  // Workaround: Unfortunately, not every "m[i]" object has a rich
1072  // enough interface to populate reinit_(domain|range)_vector. Thus,
1073  // apply an empty LinearOperator exemplar.
1074  matrices[i] =
1075  linear_operator<VectorType>(LinearOperator<VectorType>(),
1077  }
1078 }
1079 
1080 
1081 
1082 template <typename MatrixType, typename PreconditionerType, typename VectorType>
1083 template <typename MatrixType2, class DATA>
1084 inline void
1086  const MGLevelObject<MatrixType2> &m,
1087  const MGLevelObject<DATA> & data)
1088 {
1089  const unsigned int min = m.min_level();
1090  const unsigned int max = m.max_level();
1091 
1092  Assert(data.min_level() == min, ExcDimensionMismatch(data.min_level(), min));
1093  Assert(data.max_level() == max, ExcDimensionMismatch(data.max_level(), max));
1094 
1095  matrices.resize(min, max);
1096  smoothers.resize(min, max);
1097 
1098  for (unsigned int i = min; i <= max; ++i)
1099  {
1100  // Workaround: Unfortunately, not every "m[i]" object has a rich
1101  // enough interface to populate reinit_(domain|range)_vector. Thus,
1102  // apply an empty LinearOperator exemplar.
1103  matrices[i] =
1104  linear_operator<VectorType>(LinearOperator<VectorType>(),
1106  smoothers[i].initialize(Utilities::get_underlying_value(m[i]), data[i]);
1107  }
1108 }
1109 
1110 
1111 
1112 template <typename MatrixType, typename PreconditionerType, typename VectorType>
1113 template <typename MatrixType2, class DATA>
1114 inline void
1116  const MGLevelObject<MatrixType2> &m,
1117  const DATA & data,
1118  const unsigned int row,
1119  const unsigned int col)
1120 {
1121  const unsigned int min = m.min_level();
1122  const unsigned int max = m.max_level();
1123 
1124  matrices.resize(min, max);
1125  smoothers.resize(min, max);
1126 
1127  for (unsigned int i = min; i <= max; ++i)
1128  {
1129  matrices[i] = &(m[i].block(row, col));
1130  smoothers[i].initialize(m[i].block(row, col), data);
1131  }
1132 }
1133 
1134 
1135 
1136 template <typename MatrixType, typename PreconditionerType, typename VectorType>
1137 template <typename MatrixType2, class DATA>
1138 inline void
1140  const MGLevelObject<MatrixType2> &m,
1141  const MGLevelObject<DATA> & data,
1142  const unsigned int row,
1143  const unsigned int col)
1144 {
1145  const unsigned int min = m.min_level();
1146  const unsigned int max = m.max_level();
1147 
1148  Assert(data.min_level() == min, ExcDimensionMismatch(data.min_level(), min));
1149  Assert(data.max_level() == max, ExcDimensionMismatch(data.max_level(), max));
1150 
1151  matrices.resize(min, max);
1152  smoothers.resize(min, max);
1153 
1154  for (unsigned int i = min; i <= max; ++i)
1155  {
1156  matrices[i] = &(m[i].block(row, col));
1157  smoothers[i].initialize(m[i].block(row, col), data[i]);
1158  }
1159 }
1160 
1161 
1162 
1163 template <typename MatrixType, typename PreconditionerType, typename VectorType>
1164 inline void
1166  const unsigned int level,
1167  VectorType & u,
1168  const VectorType & rhs) const
1169 {
1170  unsigned int maxlevel = matrices.max_level();
1171  unsigned int steps2 = this->steps;
1172 
1173  if (this->variable)
1174  steps2 *= (1 << (maxlevel - level));
1175 
1176  typename VectorMemory<VectorType>::Pointer r(this->vector_memory);
1177  typename VectorMemory<VectorType>::Pointer d(this->vector_memory);
1178 
1179  r->reinit(u, true);
1180  d->reinit(u, true);
1181 
1182  bool T = this->transpose;
1183  if (this->symmetric && (steps2 % 2 == 0))
1184  T = false;
1185  if (this->debug > 0)
1186  deallog << 'S' << level << ' ';
1187 
1188  for (unsigned int i = 0; i < steps2; ++i)
1189  {
1190  if (T)
1191  {
1192  if (this->debug > 0)
1193  deallog << 'T';
1194  matrices[level].Tvmult(*r, u);
1195  r->sadd(-1., 1., rhs);
1196  if (this->debug > 2)
1197  deallog << ' ' << r->l2_norm() << ' ';
1198  smoothers[level].Tvmult(*d, *r);
1199  if (this->debug > 1)
1200  deallog << ' ' << d->l2_norm() << ' ';
1201  }
1202  else
1203  {
1204  if (this->debug > 0)
1205  deallog << 'N';
1206  matrices[level].vmult(*r, u);
1207  r->sadd(-1., rhs);
1208  if (this->debug > 2)
1209  deallog << ' ' << r->l2_norm() << ' ';
1210  smoothers[level].vmult(*d, *r);
1211  if (this->debug > 1)
1212  deallog << ' ' << d->l2_norm() << ' ';
1213  }
1214  u += *d;
1215  if (this->symmetric)
1216  T = !T;
1217  }
1218  if (this->debug > 0)
1219  deallog << std::endl;
1220 }
1221 
1222 
1223 
1224 template <typename MatrixType, typename PreconditionerType, typename VectorType>
1225 inline void
1227  const unsigned int level,
1228  VectorType & u,
1229  const VectorType & rhs) const
1230 {
1231  unsigned int maxlevel = matrices.max_level();
1232  unsigned int steps2 = this->steps;
1233 
1234  if (this->variable)
1235  steps2 *= (1 << (maxlevel - level));
1236 
1237  bool T = this->transpose;
1238  if (this->symmetric && (steps2 % 2 == 0))
1239  T = false;
1240  if (this->debug > 0)
1241  deallog << 'S' << level << ' ';
1242 
1243  // first step where we overwrite the result
1244  if (this->debug > 2)
1245  deallog << ' ' << rhs.l2_norm() << ' ';
1246  if (this->debug > 0)
1247  deallog << (T ? 'T' : 'N');
1248  if (T)
1249  smoothers[level].Tvmult(u, rhs);
1250  else
1251  smoothers[level].vmult(u, rhs);
1252  if (this->debug > 1)
1253  deallog << ' ' << u.l2_norm() << ' ';
1254  if (this->symmetric)
1255  T = !T;
1256 
1257  typename VectorMemory<VectorType>::Pointer r(this->vector_memory);
1258  typename VectorMemory<VectorType>::Pointer d(this->vector_memory);
1259 
1260  if (steps2 > 1)
1261  {
1262  r->reinit(u, true);
1263  d->reinit(u, true);
1264  }
1265 
1266  for (unsigned int i = 1; i < steps2; ++i)
1267  {
1268  if (T)
1269  {
1270  if (this->debug > 0)
1271  deallog << 'T';
1272  matrices[level].Tvmult(*r, u);
1273  r->sadd(-1., 1., rhs);
1274  if (this->debug > 2)
1275  deallog << ' ' << r->l2_norm() << ' ';
1276  smoothers[level].Tvmult(*d, *r);
1277  if (this->debug > 1)
1278  deallog << ' ' << d->l2_norm() << ' ';
1279  }
1280  else
1281  {
1282  if (this->debug > 0)
1283  deallog << 'N';
1284  matrices[level].vmult(*r, u);
1285  r->sadd(-1., rhs);
1286  if (this->debug > 2)
1287  deallog << ' ' << r->l2_norm() << ' ';
1288  smoothers[level].vmult(*d, *r);
1289  if (this->debug > 1)
1290  deallog << ' ' << d->l2_norm() << ' ';
1291  }
1292  u += *d;
1293  if (this->symmetric)
1294  T = !T;
1295  }
1296  if (this->debug > 0)
1297  deallog << std::endl;
1298 }
1299 
1300 
1301 
1302 template <typename MatrixType, typename PreconditionerType, typename VectorType>
1303 inline std::size_t
1305  memory_consumption() const
1306 {
1307  return sizeof(*this) + matrices.memory_consumption() +
1308  smoothers.memory_consumption() +
1309  this->vector_memory.memory_consumption();
1310 }
1311 
1312 
1313 #endif // DOXYGEN
1314 
1316 
1317 #endif
DerivativeForm< 1, spacedim, dim, Number > transpose(const DerivativeForm< 1, dim, spacedim, Number > &DF)
std::size_t memory_consumption() const
unsigned int max_level() const
unsigned int min_level() const
virtual void clear()
virtual void smooth(const unsigned int level, VectorType &u, const VectorType &rhs) const
std::size_t memory_consumption() const
MGLevelObject< PreconditionerType > smoothers
Definition: mg_smoother.h:568
void initialize_matrices(const MGLevelObject< MatrixType2 > &matrices)
virtual void smooth(const unsigned int level, VectorType &u, const VectorType &rhs) const override
void initialize(const MGLevelObject< MatrixType2 > &matrices, const typename PreconditionerType::AdditionalData &additional_data=typename PreconditionerType::AdditionalData())
virtual void apply(const unsigned int level, VectorType &u, const VectorType &rhs) const override
void clear() override
void initialize(const MGLevelObject< MatrixType2 > &matrices, const MGLevelObject< DATA > &additional_data, const unsigned int block_row, const unsigned int block_col)
void initialize(const MGLevelObject< MatrixType2 > &matrices, const DATA &additional_data, const unsigned int block_row, const unsigned int block_col)
MGLevelObject< LinearOperator< VectorType > > matrices
Definition: mg_smoother.h:581
MGSmootherPrecondition(const unsigned int steps=1, const bool variable=false, const bool symmetric=false, const bool transpose=false)
void initialize(const MGLevelObject< MatrixType2 > &matrices, const MGLevelObject< DATA > &additional_data)
virtual void smooth(const unsigned int level, VectorType &u, const VectorType &rhs) const
MGLevelObject< LinearOperator< VectorType > > matrices
Definition: mg_smoother.h:412
void initialize(const MGLevelObject< MatrixType2 > &matrices, const DATA &additional_data, const unsigned int block_row, const unsigned int block_col)
void initialize(const MGLevelObject< MatrixType2 > &matrices, const MGLevelObject< DATA > &additional_data, const unsigned int block_row, const unsigned int block_col)
std::size_t memory_consumption() const
virtual void apply(const unsigned int level, VectorType &u, const VectorType &rhs) const
MGLevelObject< RelaxationType > smoothers
Definition: mg_smoother.h:399
void initialize(const MGLevelObject< MatrixType2 > &matrices, const MGLevelObject< DATA > &additional_data)
MGSmootherRelaxation(const unsigned int steps=1, const bool variable=false, const bool symmetric=false, const bool transpose=false)
void initialize(const MGLevelObject< MatrixType2 > &matrices, const typename RelaxationType::AdditionalData &additional_data=typename RelaxationType::AdditionalData())
bool transpose
Definition: mg_smoother.h:124
GrowingVectorMemory< VectorType > vector_memory
Definition: mg_smoother.h:100
MGSmoother(const unsigned int steps=1, const bool variable=false, const bool symmetric=false, const bool transpose=false)
bool variable
Definition: mg_smoother.h:112
void set_debug(const unsigned int level)
void set_steps(const unsigned int)
unsigned int debug
Definition: mg_smoother.h:129
bool symmetric
Definition: mg_smoother.h:118
void set_symmetric(const bool)
void set_transpose(const bool)
unsigned int steps
Definition: mg_smoother.h:106
void set_variable(const bool)
VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &x, const ::VectorizedArray< Number, width > &y)
VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &x, const ::VectorizedArray< Number, width > &y)
std::size_t memory_consumption() const
SmootherRelaxation(const unsigned int steps=1, const bool variable=false, const bool symmetric=false, const bool transpose=false)
virtual void smooth(const unsigned int level, VectorType &u, const VectorType &rhs) const override
virtual void apply(const unsigned int level, VectorType &u, const VectorType &rhs) const override
void initialize(const MGLevelObject< MatrixType2 > &matrices, const MGLevelObject< DATA > &additional_data)
void clear() override
void initialize(const MGLevelObject< MatrixType2 > &matrices, const typename RelaxationType::AdditionalData &additional_data=typename RelaxationType::AdditionalData())
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:474
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:475
DerivativeForm< 1, spacedim, dim, Number > transpose(const DerivativeForm< 1, dim, spacedim, Number > &DF)
unsigned int level
Definition: grid_out.cc:4618
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define Assert(cond, exc)
Definition: exceptions.h:1586
LogStream deallog
Definition: logstream.cc:37
@ symmetric
Matrix is symmetric.
static const char T
std::enable_if_t< std::is_fundamental< T >::value, std::size_t > memory_consumption(const T &t)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
T & get_underlying_value(T &p)
Definition: utilities.h:1595
Definition: mg.h:82
auto apply(F &&fn, Tuple &&t) -> decltype(apply_impl(std::forward< F >(fn), std::forward< Tuple >(t), std::make_index_sequence< std::tuple_size< typename std::remove_reference< Tuple >::type >::value >()))
Definition: tuple.h:36