Reference documentation for deal.II version GIT fac6d58a2d 2022-05-23 17:55:02+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 
40 
47 template <typename VectorType>
48 class MGSmoother : public MGSmootherBase<VectorType>
49 {
50 public:
54  MGSmoother(const unsigned int steps = 1,
55  const bool variable = false,
56  const bool symmetric = false,
57  const bool transpose = false);
58 
62  void
63  set_steps(const unsigned int);
64 
68  void
69  set_variable(const bool);
70 
74  void
75  set_symmetric(const bool);
76 
81  void
82  set_transpose(const bool);
83 
88  void
89  set_debug(const unsigned int level);
90 
91 protected:
99 
104  unsigned int steps;
105 
110  bool variable;
111 
116  bool symmetric;
117 
122  bool transpose;
123 
127  unsigned int debug;
128 };
129 
130 
137 template <typename VectorType>
138 class MGSmootherIdentity : public MGSmootherBase<VectorType>
139 {
140 public:
146  virtual void
147  smooth(const unsigned int level, VectorType &u, const VectorType &rhs) const;
148 
149  virtual void
150  clear();
151 };
152 
153 
154 namespace mg
155 {
182  template <class RelaxationType, typename VectorType>
183  class SmootherRelaxation : public MGLevelObject<RelaxationType>,
184  public MGSmoother<VectorType>
185  {
186  public:
190  SmootherRelaxation(const unsigned int steps = 1,
191  const bool variable = false,
192  const bool symmetric = false,
193  const bool transpose = false);
194 
203  template <typename MatrixType2>
204  void
206  const typename RelaxationType::AdditionalData &additional_data =
207  typename RelaxationType::AdditionalData());
208 
216  template <typename MatrixType2, class DATA>
217  void
219  const MGLevelObject<DATA> & additional_data);
220 
224  void
225  clear() override;
226 
230  virtual void
231  smooth(const unsigned int level,
232  VectorType & u,
233  const VectorType & rhs) const override;
234 
251  virtual void
252  apply(const unsigned int level,
253  VectorType & u,
254  const VectorType & rhs) const override;
255 
259  std::size_t
261  };
262 } // namespace mg
263 
293 template <typename MatrixType, class RelaxationType, typename VectorType>
294 class MGSmootherRelaxation : public MGSmoother<VectorType>
295 {
296 public:
300  MGSmootherRelaxation(const unsigned int steps = 1,
301  const bool variable = false,
302  const bool symmetric = false,
303  const bool transpose = false);
304 
313  template <typename MatrixType2>
314  void
316  const typename RelaxationType::AdditionalData &additional_data =
317  typename RelaxationType::AdditionalData());
318 
327  template <typename MatrixType2, class DATA>
328  void
330  const MGLevelObject<DATA> & additional_data);
331 
341  template <typename MatrixType2, class DATA>
342  void
344  const DATA & additional_data,
345  const unsigned int block_row,
346  const unsigned int block_col);
347 
357  template <typename MatrixType2, class DATA>
358  void
360  const MGLevelObject<DATA> & additional_data,
361  const unsigned int block_row,
362  const unsigned int block_col);
363 
367  void
368  clear();
369 
373  virtual void
374  smooth(const unsigned int level, VectorType &u, const VectorType &rhs) const;
375 
391  virtual void
392  apply(const unsigned int level, VectorType &u, const VectorType &rhs) const;
393 
398 
402  std::size_t
404 
405 
406 private:
411 };
412 
413 
414 
443 template <typename MatrixType, typename PreconditionerType, typename VectorType>
444 class MGSmootherPrecondition : public MGSmoother<VectorType>
445 {
446 public:
450  MGSmootherPrecondition(const unsigned int steps = 1,
451  const bool variable = false,
452  const bool symmetric = false,
453  const bool transpose = false);
454 
464  template <typename MatrixType2>
465  void
467  const typename PreconditionerType::AdditionalData &
468  additional_data = typename PreconditionerType::AdditionalData());
469 
477  template <typename MatrixType2>
478  void
480 
490  template <typename MatrixType2, class DATA>
491  void
493  const MGLevelObject<DATA> & additional_data);
494 
505  template <typename MatrixType2, class DATA>
506  void
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
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 
567 
571  std::size_t
573 
574 
575 private:
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
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)
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
635 {
636  variable = flag;
637 }
638 
639 
640 template <typename VectorType>
641 inline void
643 {
644  symmetric = flag;
645 }
646 
647 
648 template <typename VectorType>
649 inline void
651 {
652  transpose = flag;
653 }
654 
655 //----------------------------------------------------------------------//
656 
657 namespace mg
658 {
659  template <class RelaxationType, typename VectorType>
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  {
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(
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
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
782  {
783  return sizeof(*this) - sizeof(MGLevelObject<RelaxationType>) +
785  this->vector_memory.memory_consumption();
786  }
787 } // namespace mg
788 
789 
790 //----------------------------------------------------------------------//
791 
792 template <typename MatrixType, class RelaxationType, typename 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
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
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
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
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
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
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
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
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>
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
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
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>(),
1049  smoothers[i].initialize(Utilities::get_underlying_value(m[i]), data);
1050  }
1051 }
1052 
1053 
1054 
1055 template <typename MatrixType, typename PreconditionerType, typename VectorType>
1056 template <typename MatrixType2>
1057 inline void
1060 {
1061  const unsigned int min = m.min_level();
1062  const unsigned int max = m.max_level();
1063 
1064  matrices.resize(min, max);
1065  smoothers.resize(min, max);
1066 
1067  for (unsigned int i = min; i <= max; ++i)
1068  {
1069  // Workaround: Unfortunately, not every "m[i]" object has a rich
1070  // enough interface to populate reinit_(domain|range)_vector. Thus,
1071  // apply an empty LinearOperator exemplar.
1072  matrices[i] =
1073  linear_operator<VectorType>(LinearOperator<VectorType>(),
1075  }
1076 }
1077 
1078 
1079 
1080 template <typename MatrixType, typename PreconditionerType, typename VectorType>
1081 template <typename MatrixType2, class DATA>
1082 inline void
1084  const MGLevelObject<MatrixType2> &m,
1085  const MGLevelObject<DATA> & data)
1086 {
1087  const unsigned int min = m.min_level();
1088  const unsigned int max = m.max_level();
1089 
1090  Assert(data.min_level() == min, ExcDimensionMismatch(data.min_level(), min));
1091  Assert(data.max_level() == max, ExcDimensionMismatch(data.max_level(), max));
1092 
1093  matrices.resize(min, max);
1094  smoothers.resize(min, max);
1095 
1096  for (unsigned int i = min; i <= max; ++i)
1097  {
1098  // Workaround: Unfortunately, not every "m[i]" object has a rich
1099  // enough interface to populate reinit_(domain|range)_vector. Thus,
1100  // apply an empty LinearOperator exemplar.
1101  matrices[i] =
1102  linear_operator<VectorType>(LinearOperator<VectorType>(),
1104  smoothers[i].initialize(Utilities::get_underlying_value(m[i]), data[i]);
1105  }
1106 }
1107 
1108 
1109 
1110 template <typename MatrixType, typename PreconditionerType, typename VectorType>
1111 template <typename MatrixType2, class DATA>
1112 inline void
1114  const MGLevelObject<MatrixType2> &m,
1115  const DATA & data,
1116  const unsigned int row,
1117  const unsigned int col)
1118 {
1119  const unsigned int min = m.min_level();
1120  const unsigned int max = m.max_level();
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);
1129  }
1130 }
1131 
1132 
1133 
1134 template <typename MatrixType, typename PreconditionerType, typename VectorType>
1135 template <typename MatrixType2, class DATA>
1136 inline void
1138  const MGLevelObject<MatrixType2> &m,
1139  const MGLevelObject<DATA> & data,
1140  const unsigned int row,
1141  const unsigned int col)
1142 {
1143  const unsigned int min = m.min_level();
1144  const unsigned int max = m.max_level();
1145 
1146  Assert(data.min_level() == min, ExcDimensionMismatch(data.min_level(), min));
1147  Assert(data.max_level() == max, ExcDimensionMismatch(data.max_level(), max));
1148 
1149  matrices.resize(min, max);
1150  smoothers.resize(min, max);
1151 
1152  for (unsigned int i = min; i <= max; ++i)
1153  {
1154  matrices[i] = &(m[i].block(row, col));
1155  smoothers[i].initialize(m[i].block(row, col), data[i]);
1156  }
1157 }
1158 
1159 
1160 
1161 template <typename MatrixType, typename PreconditionerType, typename VectorType>
1162 inline void
1164  const unsigned int level,
1165  VectorType & u,
1166  const VectorType & rhs) const
1167 {
1168  unsigned int maxlevel = matrices.max_level();
1169  unsigned int steps2 = this->steps;
1170 
1171  if (this->variable)
1172  steps2 *= (1 << (maxlevel - level));
1173 
1174  typename VectorMemory<VectorType>::Pointer r(this->vector_memory);
1175  typename VectorMemory<VectorType>::Pointer d(this->vector_memory);
1176 
1177  r->reinit(u, true);
1178  d->reinit(u, true);
1179 
1180  bool T = this->transpose;
1181  if (this->symmetric && (steps2 % 2 == 0))
1182  T = false;
1183  if (this->debug > 0)
1184  deallog << 'S' << level << ' ';
1185 
1186  for (unsigned int i = 0; i < steps2; ++i)
1187  {
1188  if (T)
1189  {
1190  if (this->debug > 0)
1191  deallog << 'T';
1192  matrices[level].Tvmult(*r, u);
1193  r->sadd(-1., 1., rhs);
1194  if (this->debug > 2)
1195  deallog << ' ' << r->l2_norm() << ' ';
1196  smoothers[level].Tvmult(*d, *r);
1197  if (this->debug > 1)
1198  deallog << ' ' << d->l2_norm() << ' ';
1199  }
1200  else
1201  {
1202  if (this->debug > 0)
1203  deallog << 'N';
1204  matrices[level].vmult(*r, u);
1205  r->sadd(-1., rhs);
1206  if (this->debug > 2)
1207  deallog << ' ' << r->l2_norm() << ' ';
1208  smoothers[level].vmult(*d, *r);
1209  if (this->debug > 1)
1210  deallog << ' ' << d->l2_norm() << ' ';
1211  }
1212  u += *d;
1213  if (this->symmetric)
1214  T = !T;
1215  }
1216  if (this->debug > 0)
1217  deallog << std::endl;
1218 }
1219 
1220 
1221 
1222 template <typename MatrixType, typename PreconditionerType, typename VectorType>
1223 inline void
1225  const unsigned int level,
1226  VectorType & u,
1227  const VectorType & rhs) const
1228 {
1229  unsigned int maxlevel = matrices.max_level();
1230  unsigned int steps2 = this->steps;
1231 
1232  if (this->variable)
1233  steps2 *= (1 << (maxlevel - level));
1234 
1235  bool T = this->transpose;
1236  if (this->symmetric && (steps2 % 2 == 0))
1237  T = false;
1238  if (this->debug > 0)
1239  deallog << 'S' << level << ' ';
1240 
1241  // first step where we overwrite the result
1242  if (this->debug > 2)
1243  deallog << ' ' << rhs.l2_norm() << ' ';
1244  if (this->debug > 0)
1245  deallog << (T ? 'T' : 'N');
1246  if (T)
1247  smoothers[level].Tvmult(u, rhs);
1248  else
1249  smoothers[level].vmult(u, rhs);
1250  if (this->debug > 1)
1251  deallog << ' ' << u.l2_norm() << ' ';
1252  if (this->symmetric)
1253  T = !T;
1254 
1255  typename VectorMemory<VectorType>::Pointer r(this->vector_memory);
1256  typename VectorMemory<VectorType>::Pointer d(this->vector_memory);
1257 
1258  if (steps2 > 1)
1259  {
1260  r->reinit(u, true);
1261  d->reinit(u, true);
1262  }
1263 
1264  for (unsigned int i = 1; i < steps2; ++i)
1265  {
1266  if (T)
1267  {
1268  if (this->debug > 0)
1269  deallog << 'T';
1270  matrices[level].Tvmult(*r, u);
1271  r->sadd(-1., 1., rhs);
1272  if (this->debug > 2)
1273  deallog << ' ' << r->l2_norm() << ' ';
1274  smoothers[level].Tvmult(*d, *r);
1275  if (this->debug > 1)
1276  deallog << ' ' << d->l2_norm() << ' ';
1277  }
1278  else
1279  {
1280  if (this->debug > 0)
1281  deallog << 'N';
1282  matrices[level].vmult(*r, u);
1283  r->sadd(-1., rhs);
1284  if (this->debug > 2)
1285  deallog << ' ' << r->l2_norm() << ' ';
1286  smoothers[level].vmult(*d, *r);
1287  if (this->debug > 1)
1288  deallog << ' ' << d->l2_norm() << ' ';
1289  }
1290  u += *d;
1291  if (this->symmetric)
1292  T = !T;
1293  }
1294  if (this->debug > 0)
1295  deallog << std::endl;
1296 }
1297 
1298 
1299 
1300 template <typename MatrixType, typename PreconditionerType, typename VectorType>
1301 inline std::size_t
1303  memory_consumption() const
1304 {
1305  return sizeof(*this) + matrices.memory_consumption() +
1306  smoothers.memory_consumption() +
1307  this->vector_memory.memory_consumption();
1308 }
1309 
1310 
1311 #endif // DOXYGEN
1312 
1314 
1315 #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:566
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:579
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:410
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:397
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:122
GrowingVectorMemory< VectorType > vector_memory
Definition: mg_smoother.h:98
MGSmoother(const unsigned int steps=1, const bool variable=false, const bool symmetric=false, const bool transpose=false)
bool variable
Definition: mg_smoother.h:110
void set_debug(const unsigned int level)
void set_steps(const unsigned int)
unsigned int debug
Definition: mg_smoother.h:127
bool symmetric
Definition: mg_smoother.h:116
void set_symmetric(const bool)
void set_transpose(const bool)
unsigned int steps
Definition: mg_smoother.h:104
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:416
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:417
DerivativeForm< 1, spacedim, dim, Number > transpose(const DerivativeForm< 1, dim, spacedim, Number > &DF)
unsigned int level
Definition: grid_out.cc:4606
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define Assert(cond, exc)
Definition: exceptions.h:1473
LogStream deallog
Definition: logstream.cc:37
@ symmetric
Matrix is symmetric.
static const char T
std::enable_if< std::is_fundamental< T >::value, std::size_t >::type 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:1427
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