Reference documentation for deal.II version 9.4.1
\(\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\}}\)
Loading...
Searching...
No Matches
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
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
32DEAL_II_NAMESPACE_OPEN
33
34/*
35 * MGSmootherBase is defined in mg_base.h
36 */
37
40
47template <typename VectorType>
48class MGSmoother : public MGSmootherBase<VectorType>
49{
50public:
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
91protected:
98 mutable GrowingVectorMemory<VectorType> vector_memory;
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
137template <typename VectorType>
138class MGSmootherIdentity : public MGSmootherBase<VectorType>
139{
140public:
146 virtual void
147 smooth(const unsigned int level, VectorType &u, const VectorType &rhs) const;
148
149 virtual void
150 clear();
151};
152
153
154namespace 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
205 initialize(const MGLevelObject<MatrixType2> & matrices,
206 const typename RelaxationType::AdditionalData &additional_data =
207 typename RelaxationType::AdditionalData());
208
216 template <typename MatrixType2, class DATA>
217 void
218 initialize(const MGLevelObject<MatrixType2> &matrices,
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
260 memory_consumption() const;
261 };
262} // namespace mg
263
293template <typename MatrixType, class RelaxationType, typename VectorType>
294class MGSmootherRelaxation : public MGSmoother<VectorType>
295{
296public:
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
315 initialize(const MGLevelObject<MatrixType2> & matrices,
316 const typename RelaxationType::AdditionalData &additional_data =
317 typename RelaxationType::AdditionalData());
318
327 template <typename MatrixType2, class DATA>
328 void
329 initialize(const MGLevelObject<MatrixType2> &matrices,
330 const MGLevelObject<DATA> & additional_data);
331
341 template <typename MatrixType2, class DATA>
342 void
343 initialize(const MGLevelObject<MatrixType2> &matrices,
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
359 initialize(const MGLevelObject<MatrixType2> &matrices,
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
397 MGLevelObject<RelaxationType> smoothers;
398
402 std::size_t
403 memory_consumption() const;
404
405
406private:
410 MGLevelObject<LinearOperator<VectorType>> matrices;
411};
412
413
414
443template <typename MatrixType, typename PreconditionerType, typename VectorType>
444class MGSmootherPrecondition : public MGSmoother<VectorType>
445{
446public:
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
466 initialize(const MGLevelObject<MatrixType2> &matrices,
467 const typename PreconditionerType::AdditionalData &
468 additional_data = typename PreconditionerType::AdditionalData());
469
477 template <typename MatrixType2>
478 void
479 initialize_matrices(const MGLevelObject<MatrixType2> &matrices);
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
575private:
579 MGLevelObject<LinearOperator<VectorType>> matrices;
580};
581
584/* ------------------------------- Inline functions --------------------------
585 */
586
587#ifndef DOXYGEN
588
589template <typename VectorType>
590inline void
591MGSmootherIdentity<VectorType>::smooth(const unsigned int,
592 VectorType &,
593 const VectorType &) const
594{}
595
596template <typename VectorType>
597inline void
598MGSmootherIdentity<VectorType>::clear()
599{}
600
601//---------------------------------------------------------------------------
602
603template <typename VectorType>
604inline 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
616template <typename VectorType>
617inline void
618MGSmoother<VectorType>::set_steps(const unsigned int s)
619{
620 steps = s;
621}
622
623
624template <typename VectorType>
625inline void
626MGSmoother<VectorType>::set_debug(const unsigned int s)
627{
628 debug = s;
629}
630
631
632template <typename VectorType>
633inline void
634MGSmoother<VectorType>::set_variable(const bool flag)
635{
636 variable = flag;
637}
638
639
640template <typename VectorType>
641inline void
642MGSmoother<VectorType>::set_symmetric(const bool flag)
643{
644 symmetric = flag;
645}
646
647
648template <typename VectorType>
649inline void
650MGSmoother<VectorType>::set_transpose(const bool flag)
651{
652 transpose = flag;
653}
654
655//----------------------------------------------------------------------//
656
657namespace 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
792template <typename MatrixType, class RelaxationType, typename VectorType>
793inline 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
803template <typename MatrixType, class RelaxationType, typename VectorType>
804inline void
805MGSmootherRelaxation<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
815template <typename MatrixType, class RelaxationType, typename VectorType>
816template <typename MatrixType2>
817inline void
818MGSmootherRelaxation<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
839template <typename MatrixType, class RelaxationType, typename VectorType>
840template <typename MatrixType2, class DATA>
841inline void
842MGSmootherRelaxation<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
866template <typename MatrixType, class RelaxationType, typename VectorType>
867template <typename MatrixType2, class DATA>
868inline void
869MGSmootherRelaxation<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
892template <typename MatrixType, class RelaxationType, typename VectorType>
893template <typename MatrixType2, class DATA>
894inline void
895MGSmootherRelaxation<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
922template <typename MatrixType, class RelaxationType, typename VectorType>
923inline void
924MGSmootherRelaxation<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
953template <typename MatrixType, class RelaxationType, typename VectorType>
954inline void
955MGSmootherRelaxation<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
991template <typename MatrixType, class RelaxationType, typename VectorType>
992inline std::size_t
993MGSmootherRelaxation<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
1004template <typename MatrixType, typename PreconditionerType, typename VectorType>
1005inline 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
1015template <typename MatrixType, typename PreconditionerType, typename VectorType>
1016inline void
1017MGSmootherPrecondition<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
1028template <typename MatrixType, typename PreconditionerType, typename VectorType>
1029template <typename MatrixType2>
1030inline void
1031MGSmootherPrecondition<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>(),
1048 Utilities::get_underlying_value(m[i]));
1049 smoothers[i].initialize(Utilities::get_underlying_value(m[i]), data);
1050 }
1051}
1052
1053
1054
1055template <typename MatrixType, typename PreconditionerType, typename VectorType>
1056template <typename MatrixType2>
1057inline void
1058MGSmootherPrecondition<MatrixType, PreconditionerType, VectorType>::
1059 initialize_matrices(const MGLevelObject<MatrixType2> &m)
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>(),
1074 Utilities::get_underlying_value(m[i]));
1075 }
1076}
1077
1078
1079
1080template <typename MatrixType, typename PreconditionerType, typename VectorType>
1081template <typename MatrixType2, class DATA>
1082inline void
1083MGSmootherPrecondition<MatrixType, PreconditionerType, VectorType>::initialize(
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>(),
1103 Utilities::get_underlying_value(m[i]));
1104 smoothers[i].initialize(Utilities::get_underlying_value(m[i]), data[i]);
1105 }
1106}
1107
1108
1109
1110template <typename MatrixType, typename PreconditionerType, typename VectorType>
1111template <typename MatrixType2, class DATA>
1112inline void
1113MGSmootherPrecondition<MatrixType, PreconditionerType, VectorType>::initialize(
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
1134template <typename MatrixType, typename PreconditionerType, typename VectorType>
1135template <typename MatrixType2, class DATA>
1136inline void
1137MGSmootherPrecondition<MatrixType, PreconditionerType, VectorType>::initialize(
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
1161template <typename MatrixType, typename PreconditionerType, typename VectorType>
1162inline void
1163MGSmootherPrecondition<MatrixType, PreconditionerType, VectorType>::smooth(
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
1222template <typename MatrixType, typename PreconditionerType, typename VectorType>
1223inline void
1224MGSmootherPrecondition<MatrixType, PreconditionerType, VectorType>::apply(
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
1300template <typename MatrixType, typename PreconditionerType, typename VectorType>
1301inline std::size_t
1302MGSmootherPrecondition<MatrixType, PreconditionerType, VectorType>::
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
1313DEAL_II_NAMESPACE_CLOSE
1314
1315#endif
MGLevelObject::memory_consumption
std::size_t memory_consumption() const
Definition: mg_level_object.h:325
MGLevelObject::max_level
unsigned int max_level() const
Definition: mg_level_object.h:298
MGLevelObject::min_level
unsigned int min_level() const
Definition: mg_level_object.h:290
MGSmootherIdentity::clear
virtual void clear()
MGSmootherIdentity::smooth
virtual void smooth(const unsigned int level, VectorType &u, const VectorType &rhs) const
MGSmootherPrecondition::memory_consumption
std::size_t memory_consumption() const
MGSmootherPrecondition::smoothers
MGLevelObject< PreconditionerType > smoothers
Definition: mg_smoother.h:566
MGSmootherPrecondition::initialize_matrices
void initialize_matrices(const MGLevelObject< MatrixType2 > &matrices)
MGSmootherPrecondition::smooth
virtual void smooth(const unsigned int level, VectorType &u, const VectorType &rhs) const override
MGSmootherPrecondition::initialize
void initialize(const MGLevelObject< MatrixType2 > &matrices, const typename PreconditionerType::AdditionalData &additional_data=typename PreconditionerType::AdditionalData())
MGSmootherPrecondition::apply
virtual void apply(const unsigned int level, VectorType &u, const VectorType &rhs) const override
MGSmootherPrecondition::clear
void clear() override
MGSmootherPrecondition::initialize
void initialize(const MGLevelObject< MatrixType2 > &matrices, const MGLevelObject< DATA > &additional_data, const unsigned int block_row, const unsigned int block_col)
MGSmootherPrecondition::initialize
void initialize(const MGLevelObject< MatrixType2 > &matrices, const DATA &additional_data, const unsigned int block_row, const unsigned int block_col)
MGSmootherPrecondition::matrices
MGLevelObject< LinearOperator< VectorType > > matrices
Definition: mg_smoother.h:579
MGSmootherPrecondition::MGSmootherPrecondition
MGSmootherPrecondition(const unsigned int steps=1, const bool variable=false, const bool symmetric=false, const bool transpose=false)
MGSmootherPrecondition::initialize
void initialize(const MGLevelObject< MatrixType2 > &matrices, const MGLevelObject< DATA > &additional_data)
MGSmootherRelaxation::smooth
virtual void smooth(const unsigned int level, VectorType &u, const VectorType &rhs) const
MGSmootherRelaxation::matrices
MGLevelObject< LinearOperator< VectorType > > matrices
Definition: mg_smoother.h:410
MGSmootherRelaxation::initialize
void initialize(const MGLevelObject< MatrixType2 > &matrices, const DATA &additional_data, const unsigned int block_row, const unsigned int block_col)
MGSmootherRelaxation::initialize
void initialize(const MGLevelObject< MatrixType2 > &matrices, const MGLevelObject< DATA > &additional_data, const unsigned int block_row, const unsigned int block_col)
MGSmootherRelaxation::memory_consumption
std::size_t memory_consumption() const
MGSmootherRelaxation::apply
virtual void apply(const unsigned int level, VectorType &u, const VectorType &rhs) const
MGSmootherRelaxation::smoothers
MGLevelObject< RelaxationType > smoothers
Definition: mg_smoother.h:397
MGSmootherRelaxation::initialize
void initialize(const MGLevelObject< MatrixType2 > &matrices, const MGLevelObject< DATA > &additional_data)
MGSmootherRelaxation::MGSmootherRelaxation
MGSmootherRelaxation(const unsigned int steps=1, const bool variable=false, const bool symmetric=false, const bool transpose=false)
MGSmootherRelaxation::initialize
void initialize(const MGLevelObject< MatrixType2 > &matrices, const typename RelaxationType::AdditionalData &additional_data=typename RelaxationType::AdditionalData())
MGSmoother::transpose
bool transpose
Definition: mg_smoother.h:122
MGSmoother::vector_memory
GrowingVectorMemory< VectorType > vector_memory
Definition: mg_smoother.h:98
MGSmoother::MGSmoother
MGSmoother(const unsigned int steps=1, const bool variable=false, const bool symmetric=false, const bool transpose=false)
MGSmoother::variable
bool variable
Definition: mg_smoother.h:110
MGSmoother::set_debug
void set_debug(const unsigned int level)
MGSmoother::set_steps
void set_steps(const unsigned int)
MGSmoother::debug
unsigned int debug
Definition: mg_smoother.h:127
MGSmoother::symmetric
bool symmetric
Definition: mg_smoother.h:116
MGSmoother::set_symmetric
void set_symmetric(const bool)
MGSmoother::set_transpose
void set_transpose(const bool)
MGSmoother::steps
unsigned int steps
Definition: mg_smoother.h:104
MGSmoother::set_variable
void set_variable(const bool)
mg::SmootherRelaxation::memory_consumption
std::size_t memory_consumption() const
mg::SmootherRelaxation::SmootherRelaxation
SmootherRelaxation(const unsigned int steps=1, const bool variable=false, const bool symmetric=false, const bool transpose=false)
mg::SmootherRelaxation::smooth
virtual void smooth(const unsigned int level, VectorType &u, const VectorType &rhs) const override
mg::SmootherRelaxation::apply
virtual void apply(const unsigned int level, VectorType &u, const VectorType &rhs) const override
mg::SmootherRelaxation::initialize
void initialize(const MGLevelObject< MatrixType2 > &matrices, const MGLevelObject< DATA > &additional_data)
mg::SmootherRelaxation::clear
void clear() override
mg::SmootherRelaxation::initialize
void initialize(const MGLevelObject< MatrixType2 > &matrices, const typename RelaxationType::AdditionalData &additional_data=typename RelaxationType::AdditionalData())
DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:442
DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
transpose
DerivativeForm< 1, spacedim, dim, Number > transpose(const DerivativeForm< 1, dim, spacedim, Number > &DF)
Definition: derivative_form.h:538
level
unsigned int level
Definition: grid_out.cc:4606
Assert
#define Assert(cond, exc)
Definition: exceptions.h:1473
StandardExceptions::ExcDimensionMismatch
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
vector_memory.h
linear_operator.h
deallog
LogStream deallog
Definition: logstream.cc:37
mg_base.h
mg_level_object.h
LAPACKSupport::symmetric
@ symmetric
Matrix is symmetric.
Definition: lapack_support.h:113
Physics::Elasticity::Kinematics::d
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
Utilities::get_underlying_value
T & get_underlying_value(T &p)
Definition: utilities.h:1753
mg
Definition: mg.h:82
std::min
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
Definition: vectorization.h:5605
std::max
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
Definition: vectorization.h:5588
Generated by doxygen 1.9.5