Reference documentation for deal.II version 9.3.3
\(\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 - 2021 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
479 template <typename MatrixType2, class DATA>
480 void
481 initialize(const MGLevelObject<MatrixType2> &matrices,
482 const MGLevelObject<DATA> & additional_data);
483
494 template <typename MatrixType2, class DATA>
495 void
496 initialize(const MGLevelObject<MatrixType2> &matrices,
497 const DATA & additional_data,
498 const unsigned int block_row,
499 const unsigned int block_col);
500
511 template <typename MatrixType2, class DATA>
512 void
513 initialize(const MGLevelObject<MatrixType2> &matrices,
514 const MGLevelObject<DATA> & additional_data,
515 const unsigned int block_row,
516 const unsigned int block_col);
517
521 void
522 clear() override;
523
527 virtual void
528 smooth(const unsigned int level,
529 VectorType & u,
530 const VectorType & rhs) const override;
531
547 virtual void
548 apply(const unsigned int level,
549 VectorType & u,
550 const VectorType & rhs) const override;
551
555 MGLevelObject<PreconditionerType> smoothers;
556
560 std::size_t
561 memory_consumption() const;
562
563
564private:
568 MGLevelObject<LinearOperator<VectorType>> matrices;
569};
570
573/* ------------------------------- Inline functions --------------------------
574 */
575
576#ifndef DOXYGEN
577
578template <typename VectorType>
579inline void
580MGSmootherIdentity<VectorType>::smooth(const unsigned int,
581 VectorType &,
582 const VectorType &) const
583{}
584
585template <typename VectorType>
586inline void
587MGSmootherIdentity<VectorType>::clear()
588{}
589
590//---------------------------------------------------------------------------
591
592template <typename VectorType>
593inline MGSmoother<VectorType>::MGSmoother(const unsigned int steps,
594 const bool variable,
595 const bool symmetric,
596 const bool transpose)
597 : steps(steps)
598 , variable(variable)
599 , symmetric(symmetric)
600 , transpose(transpose)
601 , debug(0)
602{}
603
604
605template <typename VectorType>
606inline void
607MGSmoother<VectorType>::set_steps(const unsigned int s)
608{
609 steps = s;
610}
611
612
613template <typename VectorType>
614inline void
615MGSmoother<VectorType>::set_debug(const unsigned int s)
616{
617 debug = s;
618}
619
620
621template <typename VectorType>
622inline void
623MGSmoother<VectorType>::set_variable(const bool flag)
624{
625 variable = flag;
626}
627
628
629template <typename VectorType>
630inline void
631MGSmoother<VectorType>::set_symmetric(const bool flag)
632{
633 symmetric = flag;
634}
635
636
637template <typename VectorType>
638inline void
639MGSmoother<VectorType>::set_transpose(const bool flag)
640{
641 transpose = flag;
642}
643
644//----------------------------------------------------------------------//
645
646namespace mg
647{
648 template <class RelaxationType, typename VectorType>
649 inline SmootherRelaxation<RelaxationType, VectorType>::SmootherRelaxation(
650 const unsigned int steps,
651 const bool variable,
652 const bool symmetric,
653 const bool transpose)
654 : MGSmoother<VectorType>(steps, variable, symmetric, transpose)
655 {}
656
657
658 template <class RelaxationType, typename VectorType>
659 inline void
660 SmootherRelaxation<RelaxationType, VectorType>::clear()
661 {
662 MGLevelObject<RelaxationType>::clear_elements();
663 }
664
665
666 template <class RelaxationType, typename VectorType>
667 template <typename MatrixType2>
668 inline void
669 SmootherRelaxation<RelaxationType, VectorType>::initialize(
670 const MGLevelObject<MatrixType2> & m,
671 const typename RelaxationType::AdditionalData &data)
672 {
673 const unsigned int min = m.min_level();
674 const unsigned int max = m.max_level();
675
676 this->resize(min, max);
677
678 for (unsigned int i = min; i <= max; ++i)
679 (*this)[i].initialize(m[i], data);
680 }
681
682
683 template <class RelaxationType, typename VectorType>
684 template <typename MatrixType2, class DATA>
685 inline void
686 SmootherRelaxation<RelaxationType, VectorType>::initialize(
687 const MGLevelObject<MatrixType2> &m,
688 const MGLevelObject<DATA> & data)
689 {
690 const unsigned int min = std::max(m.min_level(), data.min_level());
691 const unsigned int max = std::min(m.max_level(), data.max_level());
692
693 this->resize(min, max);
694
695 for (unsigned int i = min; i <= max; ++i)
696 (*this)[i].initialize(m[i], data[i]);
697 }
698
699
700 template <class RelaxationType, typename VectorType>
701 inline void
702 SmootherRelaxation<RelaxationType, VectorType>::smooth(
703 const unsigned int level,
704 VectorType & u,
705 const VectorType & rhs) const
706 {
707 unsigned int maxlevel = this->max_level();
708 unsigned int steps2 = this->steps;
709
710 if (this->variable)
711 steps2 *= (1 << (maxlevel - level));
712
713 bool T = this->transpose;
714 if (this->symmetric && (steps2 % 2 == 0))
715 T = false;
716 if (this->debug > 0)
717 deallog << 'S' << level << ' ';
718
719 for (unsigned int i = 0; i < steps2; ++i)
720 {
721 if (T)
722 (*this)[level].Tstep(u, rhs);
723 else
724 (*this)[level].step(u, rhs);
725 if (this->symmetric)
726 T = !T;
727 }
728 }
729
730
731 template <class RelaxationType, typename VectorType>
732 inline void
733 SmootherRelaxation<RelaxationType, VectorType>::apply(
734 const unsigned int level,
735 VectorType & u,
736 const VectorType & rhs) const
737 {
738 unsigned int maxlevel = this->max_level();
739 unsigned int steps2 = this->steps;
740
741 if (this->variable)
742 steps2 *= (1 << (maxlevel - level));
743
744 bool T = this->transpose;
745 if (this->symmetric && (steps2 % 2 == 0))
746 T = false;
747 if (this->debug > 0)
748 deallog << 'S' << level << ' ';
749
750 if (T)
751 (*this)[level].Tvmult(u, rhs);
752 else
753 (*this)[level].vmult(u, rhs);
754 if (this->symmetric)
755 T = !T;
756 for (unsigned int i = 1; i < steps2; ++i)
757 {
758 if (T)
759 (*this)[level].Tstep(u, rhs);
760 else
761 (*this)[level].step(u, rhs);
762 if (this->symmetric)
763 T = !T;
764 }
765 }
766
767
768 template <class RelaxationType, typename VectorType>
769 inline std::size_t
770 SmootherRelaxation<RelaxationType, VectorType>::memory_consumption() const
771 {
772 return sizeof(*this) - sizeof(MGLevelObject<RelaxationType>) +
773 MGLevelObject<RelaxationType>::memory_consumption() +
774 this->vector_memory.memory_consumption();
775 }
776} // namespace mg
777
778
779//----------------------------------------------------------------------//
780
781template <typename MatrixType, class RelaxationType, typename VectorType>
782inline MGSmootherRelaxation<MatrixType, RelaxationType, VectorType>::
783 MGSmootherRelaxation(const unsigned int steps,
784 const bool variable,
785 const bool symmetric,
786 const bool transpose)
787 : MGSmoother<VectorType>(steps, variable, symmetric, transpose)
788{}
789
790
791
792template <typename MatrixType, class RelaxationType, typename VectorType>
793inline void
794MGSmootherRelaxation<MatrixType, RelaxationType, VectorType>::clear()
795{
796 smoothers.clear_elements();
797
798 unsigned int i = matrices.min_level(), max_level = matrices.max_level();
799 for (; i <= max_level; ++i)
800 matrices[i] = LinearOperator<VectorType>();
801}
802
803
804template <typename MatrixType, class RelaxationType, typename VectorType>
805template <typename MatrixType2>
806inline void
807MGSmootherRelaxation<MatrixType, RelaxationType, VectorType>::initialize(
808 const MGLevelObject<MatrixType2> & m,
809 const typename RelaxationType::AdditionalData &data)
810{
811 const unsigned int min = m.min_level();
812 const unsigned int max = m.max_level();
813
814 matrices.resize(min, max);
815 smoothers.resize(min, max);
816
817 for (unsigned int i = min; i <= max; ++i)
818 {
819 // Workaround: Unfortunately, not every "m[i]" object has a rich
820 // enough interface to populate reinit_(domain|range)_vector. Thus,
821 // apply an empty LinearOperator exemplar.
822 matrices[i] =
823 linear_operator<VectorType>(LinearOperator<VectorType>(), m[i]);
824 smoothers[i].initialize(m[i], data);
825 }
826}
827
828template <typename MatrixType, class RelaxationType, typename VectorType>
829template <typename MatrixType2, class DATA>
830inline void
831MGSmootherRelaxation<MatrixType, RelaxationType, VectorType>::initialize(
832 const MGLevelObject<MatrixType2> &m,
833 const MGLevelObject<DATA> & data)
834{
835 const unsigned int min = m.min_level();
836 const unsigned int max = m.max_level();
837
838 Assert(data.min_level() == min, ExcDimensionMismatch(data.min_level(), min));
839 Assert(data.max_level() == max, ExcDimensionMismatch(data.max_level(), max));
840
841 matrices.resize(min, max);
842 smoothers.resize(min, max);
843
844 for (unsigned int i = min; i <= max; ++i)
845 {
846 // Workaround: Unfortunately, not every "m[i]" object has a rich
847 // enough interface to populate reinit_(domain|range)_vector. Thus,
848 // apply an empty LinearOperator exemplar.
849 matrices[i] =
850 linear_operator<VectorType>(LinearOperator<VectorType>(), m[i]);
851 smoothers[i].initialize(m[i], data[i]);
852 }
853}
854
855template <typename MatrixType, class RelaxationType, typename VectorType>
856template <typename MatrixType2, class DATA>
857inline void
858MGSmootherRelaxation<MatrixType, RelaxationType, VectorType>::initialize(
859 const MGLevelObject<MatrixType2> &m,
860 const DATA & data,
861 const unsigned int row,
862 const unsigned int col)
863{
864 const unsigned int min = m.min_level();
865 const unsigned int max = m.max_level();
866
867 matrices.resize(min, max);
868 smoothers.resize(min, max);
869
870 for (unsigned int i = min; i <= max; ++i)
871 {
872 // Workaround: Unfortunately, not every "m[i]" object has a rich
873 // enough interface to populate reinit_(domain|range)_vector. Thus,
874 // apply an empty LinearOperator exemplar.
875 matrices[i] = linear_operator<VectorType>(LinearOperator<VectorType>(),
876 m[i].block(row, col));
877 smoothers[i].initialize(m[i].block(row, col), data);
878 }
879}
880
881template <typename MatrixType, class RelaxationType, typename VectorType>
882template <typename MatrixType2, class DATA>
883inline void
884MGSmootherRelaxation<MatrixType, RelaxationType, VectorType>::initialize(
885 const MGLevelObject<MatrixType2> &m,
886 const MGLevelObject<DATA> & data,
887 const unsigned int row,
888 const unsigned int col)
889{
890 const unsigned int min = m.min_level();
891 const unsigned int max = m.max_level();
892
893 Assert(data.min_level() == min, ExcDimensionMismatch(data.min_level(), min));
894 Assert(data.max_level() == max, ExcDimensionMismatch(data.max_level(), max));
895
896 matrices.resize(min, max);
897 smoothers.resize(min, max);
898
899 for (unsigned int i = min; i <= max; ++i)
900 {
901 // Workaround: Unfortunately, not every "m[i]" object has a rich
902 // enough interface to populate reinit_(domain|range)_vector. Thus,
903 // apply an empty LinearOperator exemplar.
904 matrices[i] = linear_operator<VectorType>(LinearOperator<VectorType>(),
905 m[i].block(row, col));
906 smoothers[i].initialize(m[i].block(row, col), data[i]);
907 }
908}
909
910
911template <typename MatrixType, class RelaxationType, typename VectorType>
912inline void
913MGSmootherRelaxation<MatrixType, RelaxationType, VectorType>::smooth(
914 const unsigned int level,
915 VectorType & u,
916 const VectorType & rhs) const
917{
918 unsigned int maxlevel = smoothers.max_level();
919 unsigned int steps2 = this->steps;
920
921 if (this->variable)
922 steps2 *= (1 << (maxlevel - level));
923
924 bool T = this->transpose;
925 if (this->symmetric && (steps2 % 2 == 0))
926 T = false;
927 if (this->debug > 0)
928 deallog << 'S' << level << ' ';
929
930 for (unsigned int i = 0; i < steps2; ++i)
931 {
932 if (T)
933 smoothers[level].Tstep(u, rhs);
934 else
935 smoothers[level].step(u, rhs);
936 if (this->symmetric)
937 T = !T;
938 }
939}
940
941
942template <typename MatrixType, class RelaxationType, typename VectorType>
943inline void
944MGSmootherRelaxation<MatrixType, RelaxationType, VectorType>::apply(
945 const unsigned int level,
946 VectorType & u,
947 const VectorType & rhs) const
948{
949 unsigned int maxlevel = smoothers.max_level();
950 unsigned int steps2 = this->steps;
951
952 if (this->variable)
953 steps2 *= (1 << (maxlevel - level));
954
955 bool T = this->transpose;
956 if (this->symmetric && (steps2 % 2 == 0))
957 T = false;
958 if (this->debug > 0)
959 deallog << 'S' << level << ' ';
960
961 if (T)
962 smoothers[level].Tvmult(u, rhs);
963 else
964 smoothers[level].vmult(u, rhs);
965 if (this->symmetric)
966 T = !T;
967 for (unsigned int i = 1; i < steps2; ++i)
968 {
969 if (T)
970 smoothers[level].Tstep(u, rhs);
971 else
972 smoothers[level].step(u, rhs);
973 if (this->symmetric)
974 T = !T;
975 }
976}
977
978
979
980template <typename MatrixType, class RelaxationType, typename VectorType>
981inline std::size_t
982MGSmootherRelaxation<MatrixType, RelaxationType, VectorType>::
983 memory_consumption() const
984{
985 return sizeof(*this) + matrices.memory_consumption() +
986 smoothers.memory_consumption() +
987 this->vector_memory.memory_consumption();
988}
989
990
991//----------------------------------------------------------------------//
992
993template <typename MatrixType, typename PreconditionerType, typename VectorType>
994inline MGSmootherPrecondition<MatrixType, PreconditionerType, VectorType>::
995 MGSmootherPrecondition(const unsigned int steps,
996 const bool variable,
997 const bool symmetric,
998 const bool transpose)
999 : MGSmoother<VectorType>(steps, variable, symmetric, transpose)
1000{}
1001
1002
1003
1004template <typename MatrixType, typename PreconditionerType, typename VectorType>
1005inline void
1006MGSmootherPrecondition<MatrixType, PreconditionerType, VectorType>::clear()
1007{
1008 smoothers.clear_elements();
1009
1010 unsigned int i = matrices.min_level(), max_level = matrices.max_level();
1011 for (; i <= max_level; ++i)
1012 matrices[i] = LinearOperator<VectorType>();
1013}
1014
1015
1016
1017template <typename MatrixType, typename PreconditionerType, typename VectorType>
1018template <typename MatrixType2>
1019inline void
1020MGSmootherPrecondition<MatrixType, PreconditionerType, VectorType>::initialize(
1021 const MGLevelObject<MatrixType2> & m,
1022 const typename PreconditionerType::AdditionalData &data)
1023{
1024 const unsigned int min = m.min_level();
1025 const unsigned int max = m.max_level();
1026
1027 matrices.resize(min, max);
1028 smoothers.resize(min, max);
1029
1030 for (unsigned int i = min; i <= max; ++i)
1031 {
1032 // Workaround: Unfortunately, not every "m[i]" object has a rich
1033 // enough interface to populate reinit_(domain|range)_vector. Thus,
1034 // apply an empty LinearOperator exemplar.
1035 matrices[i] =
1036 linear_operator<VectorType>(LinearOperator<VectorType>(), m[i]);
1037 smoothers[i].initialize(m[i], data);
1038 }
1039}
1040
1041
1042
1043template <typename MatrixType, typename PreconditionerType, typename VectorType>
1044template <typename MatrixType2, class DATA>
1045inline void
1046MGSmootherPrecondition<MatrixType, PreconditionerType, VectorType>::initialize(
1047 const MGLevelObject<MatrixType2> &m,
1048 const MGLevelObject<DATA> & data)
1049{
1050 const unsigned int min = m.min_level();
1051 const unsigned int max = m.max_level();
1052
1053 Assert(data.min_level() == min, ExcDimensionMismatch(data.min_level(), min));
1054 Assert(data.max_level() == max, ExcDimensionMismatch(data.max_level(), max));
1055
1056 matrices.resize(min, max);
1057 smoothers.resize(min, max);
1058
1059 for (unsigned int i = min; i <= max; ++i)
1060 {
1061 // Workaround: Unfortunately, not every "m[i]" object has a rich
1062 // enough interface to populate reinit_(domain|range)_vector. Thus,
1063 // apply an empty LinearOperator exemplar.
1064 matrices[i] =
1065 linear_operator<VectorType>(LinearOperator<VectorType>(),
1066 Utilities::get_underlying_value(m[i]));
1067 smoothers[i].initialize(Utilities::get_underlying_value(m[i]), data[i]);
1068 }
1069}
1070
1071
1072
1073template <typename MatrixType, typename PreconditionerType, typename VectorType>
1074template <typename MatrixType2, class DATA>
1075inline void
1076MGSmootherPrecondition<MatrixType, PreconditionerType, VectorType>::initialize(
1077 const MGLevelObject<MatrixType2> &m,
1078 const DATA & data,
1079 const unsigned int row,
1080 const unsigned int col)
1081{
1082 const unsigned int min = m.min_level();
1083 const unsigned int max = m.max_level();
1084
1085 matrices.resize(min, max);
1086 smoothers.resize(min, max);
1087
1088 for (unsigned int i = min; i <= max; ++i)
1089 {
1090 matrices[i] = &(m[i].block(row, col));
1091 smoothers[i].initialize(m[i].block(row, col), data);
1092 }
1093}
1094
1095
1096
1097template <typename MatrixType, typename PreconditionerType, typename VectorType>
1098template <typename MatrixType2, class DATA>
1099inline void
1100MGSmootherPrecondition<MatrixType, PreconditionerType, VectorType>::initialize(
1101 const MGLevelObject<MatrixType2> &m,
1102 const MGLevelObject<DATA> & data,
1103 const unsigned int row,
1104 const unsigned int col)
1105{
1106 const unsigned int min = m.min_level();
1107 const unsigned int max = m.max_level();
1108
1109 Assert(data.min_level() == min, ExcDimensionMismatch(data.min_level(), min));
1110 Assert(data.max_level() == max, ExcDimensionMismatch(data.max_level(), max));
1111
1112 matrices.resize(min, max);
1113 smoothers.resize(min, max);
1114
1115 for (unsigned int i = min; i <= max; ++i)
1116 {
1117 matrices[i] = &(m[i].block(row, col));
1118 smoothers[i].initialize(m[i].block(row, col), data[i]);
1119 }
1120}
1121
1122
1123
1124template <typename MatrixType, typename PreconditionerType, typename VectorType>
1125inline void
1126MGSmootherPrecondition<MatrixType, PreconditionerType, VectorType>::smooth(
1127 const unsigned int level,
1128 VectorType & u,
1129 const VectorType & rhs) const
1130{
1131 unsigned int maxlevel = matrices.max_level();
1132 unsigned int steps2 = this->steps;
1133
1134 if (this->variable)
1135 steps2 *= (1 << (maxlevel - level));
1136
1137 typename VectorMemory<VectorType>::Pointer r(this->vector_memory);
1138 typename VectorMemory<VectorType>::Pointer d(this->vector_memory);
1139
1140 r->reinit(u, true);
1141 d->reinit(u, true);
1142
1143 bool T = this->transpose;
1144 if (this->symmetric && (steps2 % 2 == 0))
1145 T = false;
1146 if (this->debug > 0)
1147 deallog << 'S' << level << ' ';
1148
1149 for (unsigned int i = 0; i < steps2; ++i)
1150 {
1151 if (T)
1152 {
1153 if (this->debug > 0)
1154 deallog << 'T';
1155 matrices[level].Tvmult(*r, u);
1156 r->sadd(-1., 1., rhs);
1157 if (this->debug > 2)
1158 deallog << ' ' << r->l2_norm() << ' ';
1159 smoothers[level].Tvmult(*d, *r);
1160 if (this->debug > 1)
1161 deallog << ' ' << d->l2_norm() << ' ';
1162 }
1163 else
1164 {
1165 if (this->debug > 0)
1166 deallog << 'N';
1167 matrices[level].vmult(*r, u);
1168 r->sadd(-1., rhs);
1169 if (this->debug > 2)
1170 deallog << ' ' << r->l2_norm() << ' ';
1171 smoothers[level].vmult(*d, *r);
1172 if (this->debug > 1)
1173 deallog << ' ' << d->l2_norm() << ' ';
1174 }
1175 u += *d;
1176 if (this->symmetric)
1177 T = !T;
1178 }
1179 if (this->debug > 0)
1180 deallog << std::endl;
1181}
1182
1183
1184
1185template <typename MatrixType, typename PreconditionerType, typename VectorType>
1186inline void
1187MGSmootherPrecondition<MatrixType, PreconditionerType, VectorType>::apply(
1188 const unsigned int level,
1189 VectorType & u,
1190 const VectorType & rhs) const
1191{
1192 unsigned int maxlevel = matrices.max_level();
1193 unsigned int steps2 = this->steps;
1194
1195 if (this->variable)
1196 steps2 *= (1 << (maxlevel - level));
1197
1198 bool T = this->transpose;
1199 if (this->symmetric && (steps2 % 2 == 0))
1200 T = false;
1201 if (this->debug > 0)
1202 deallog << 'S' << level << ' ';
1203
1204 // first step where we overwrite the result
1205 if (this->debug > 2)
1206 deallog << ' ' << rhs.l2_norm() << ' ';
1207 if (this->debug > 0)
1208 deallog << (T ? 'T' : 'N');
1209 if (T)
1210 smoothers[level].Tvmult(u, rhs);
1211 else
1212 smoothers[level].vmult(u, rhs);
1213 if (this->debug > 1)
1214 deallog << ' ' << u.l2_norm() << ' ';
1215 if (this->symmetric)
1216 T = !T;
1217
1218 typename VectorMemory<VectorType>::Pointer r(this->vector_memory);
1219 typename VectorMemory<VectorType>::Pointer d(this->vector_memory);
1220
1221 if (steps2 > 1)
1222 {
1223 r->reinit(u, true);
1224 d->reinit(u, true);
1225 }
1226
1227 for (unsigned int i = 1; i < steps2; ++i)
1228 {
1229 if (T)
1230 {
1231 if (this->debug > 0)
1232 deallog << 'T';
1233 matrices[level].Tvmult(*r, u);
1234 r->sadd(-1., 1., rhs);
1235 if (this->debug > 2)
1236 deallog << ' ' << r->l2_norm() << ' ';
1237 smoothers[level].Tvmult(*d, *r);
1238 if (this->debug > 1)
1239 deallog << ' ' << d->l2_norm() << ' ';
1240 }
1241 else
1242 {
1243 if (this->debug > 0)
1244 deallog << 'N';
1245 matrices[level].vmult(*r, u);
1246 r->sadd(-1., rhs);
1247 if (this->debug > 2)
1248 deallog << ' ' << r->l2_norm() << ' ';
1249 smoothers[level].vmult(*d, *r);
1250 if (this->debug > 1)
1251 deallog << ' ' << d->l2_norm() << ' ';
1252 }
1253 u += *d;
1254 if (this->symmetric)
1255 T = !T;
1256 }
1257 if (this->debug > 0)
1258 deallog << std::endl;
1259}
1260
1261
1262
1263template <typename MatrixType, typename PreconditionerType, typename VectorType>
1264inline std::size_t
1265MGSmootherPrecondition<MatrixType, PreconditionerType, VectorType>::
1266 memory_consumption() const
1267{
1268 return sizeof(*this) + matrices.memory_consumption() +
1269 smoothers.memory_consumption() +
1270 this->vector_memory.memory_consumption();
1271}
1272
1273
1274#endif // DOXYGEN
1275
1276DEAL_II_NAMESPACE_CLOSE
1277
1278#endif
DerivativeForm::transpose
DerivativeForm< 1, spacedim, dim, Number > transpose(const DerivativeForm< 1, dim, spacedim, Number > &DF)
Definition: derivative_form.h:536
MGLevelObject::memory_consumption
std::size_t memory_consumption() const
Definition: mg_level_object.h:293
MGLevelObject::max_level
unsigned int max_level() const
Definition: mg_level_object.h:274
MGLevelObject::min_level
unsigned int min_level() const
Definition: mg_level_object.h:266
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:555
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:568
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)
VectorizedArray::max
VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &x, const ::VectorizedArray< Number, width > &y)
Definition: vectorization.h:5484
VectorizedArray::min
VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &x, const ::VectorizedArray< Number, width > &y)
Definition: vectorization.h:5501
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:402
DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:403
transpose
DerivativeForm< 1, spacedim, dim, Number > transpose(const DerivativeForm< 1, dim, spacedim, Number > &DF)
Definition: derivative_form.h:536
level
unsigned int level
Definition: grid_out.cc:4590
Assert
#define Assert(cond, exc)
Definition: exceptions.h:1465
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
LAPACKSupport::T
static const char T
Definition: lapack_support.h:165
MemoryConsumption::memory_consumption
std::enable_if< std::is_fundamental< T >::value, std::size_t >::type memory_consumption(const T &t)
Definition: memory_consumption.h:266
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:1422
mg
Definition: mg.h:82
std::min
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
Definition: vectorization.h:5501
std::max
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
Definition: vectorization.h:5484
Generated by doxygen 1.9.4