deal.II version GIT relicensing-6842-g793a97d2aa 2026-10-02 14:00:01+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\}}\)
Loading...
Searching...
No Matches
mg_smoother.h
Go to the documentation of this file.
1// -----------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
4// Copyright (C) 1999 - 2025 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Detailed license information governing the source code and contributions
9// can be found in LICENSE.md and CONTRIBUTING.md at the top level directory.
10//
11// -----------------------------------------------------------------------------
12
13#ifndef dealii_mg_smoother_h
14#define dealii_mg_smoother_h
15
16
17#include <deal.II/base/config.h>
18
19#include <deal.II/base/logstream.h>
20#include <deal.II/base/mg_level_object.h>
21#include <deal.II/base/observer_pointer.h>
22
23#include <deal.II/lac/linear_operator.h>
24#include <deal.II/lac/vector_memory.h>
25
26#include <deal.II/multigrid/mg_base.h>
27
28#include <vector>
29
30DEAL_II_NAMESPACE_OPEN
31
32/*
33 * MGSmootherBase is defined in mg_base.h
34 */
35
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 <typename 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, typename DataType>
217 void
218 initialize(const MGLevelObject<MatrixType2> &matrices,
219 const MGLevelObject<DataType> &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, typename 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, typename DataType>
328 void
329 initialize(const MGLevelObject<MatrixType2> &matrices,
330 const MGLevelObject<DataType> &additional_data);
331
341 template <typename MatrixType2, typename DataType>
342 void
343 initialize(const MGLevelObject<MatrixType2> &matrices,
344 const DataType &additional_data,
345 const unsigned int block_row,
346 const unsigned int block_col);
347
357 template <typename MatrixType2, typename DataType>
358 void
359 initialize(const MGLevelObject<MatrixType2> &matrices,
360 const MGLevelObject<DataType> &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, typename DataType>
491 void
492 initialize(const MGLevelObject<MatrixType2> &matrices,
493 const MGLevelObject<DataType> &additional_data);
494
505 template <typename MatrixType2, typename DataType>
506 void
507 initialize(const MGLevelObject<MatrixType2> &matrices,
508 const DataType &additional_data,
509 const unsigned int block_row,
510 const unsigned int block_col);
511
522 template <typename MatrixType2, typename DataType>
523 void
524 initialize(const MGLevelObject<MatrixType2> &matrices,
525 const MGLevelObject<DataType> &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 <typename 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 <typename RelaxationType, typename VectorType>
670 inline void
671 SmootherRelaxation<RelaxationType, VectorType>::clear()
672 {
673 MGLevelObject<RelaxationType>::clear_elements();
674 }
675
676
677 template <typename 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(Utilities::get_underlying_value(m[i]), data);
691 }
692
693
694 template <typename RelaxationType, typename VectorType>
695 template <typename MatrixType2, typename DataType>
696 inline void
697 SmootherRelaxation<RelaxationType, VectorType>::initialize(
698 const MGLevelObject<MatrixType2> &m,
699 const MGLevelObject<DataType> &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(Utilities::get_underlying_value(m[i]), data[i]);
708 }
709
710
711 template <typename 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 <typename 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 <typename 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, typename 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, typename 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, typename 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>(),
835 Utilities::get_underlying_value(m[i]));
836 smoothers[i].initialize(Utilities::get_underlying_value(m[i]), data);
837 }
838}
839
840template <typename MatrixType, typename RelaxationType, typename VectorType>
841template <typename MatrixType2, typename DataType>
842inline void
843MGSmootherRelaxation<MatrixType, RelaxationType, VectorType>::initialize(
844 const MGLevelObject<MatrixType2> &m,
845 const MGLevelObject<DataType> &data)
846{
847 const unsigned int min = m.min_level();
848 const unsigned int max = m.max_level();
849
850 Assert(data.min_level() == min, ExcDimensionMismatch(data.min_level(), min));
851 Assert(data.max_level() == max, ExcDimensionMismatch(data.max_level(), max));
852
853 matrices.resize(min, max);
854 smoothers.resize(min, max);
855
856 for (unsigned int i = min; i <= max; ++i)
857 {
858 // Workaround: Unfortunately, not every "m[i]" object has a rich
859 // enough interface to populate reinit_(domain|range)_vector. Thus,
860 // apply an empty LinearOperator exemplar.
861 matrices[i] =
862 linear_operator<VectorType>(LinearOperator<VectorType>(),
863 Utilities::get_underlying_value(m[i]));
864 smoothers[i].initialize(Utilities::get_underlying_value(m[i]), data[i]);
865 }
866}
867
868template <typename MatrixType, typename RelaxationType, typename VectorType>
869template <typename MatrixType2, typename DataType>
870inline void
871MGSmootherRelaxation<MatrixType, RelaxationType, VectorType>::initialize(
872 const MGLevelObject<MatrixType2> &m,
873 const DataType &data,
874 const unsigned int row,
875 const unsigned int col)
876{
877 const unsigned int min = m.min_level();
878 const unsigned int max = m.max_level();
879
880 matrices.resize(min, max);
881 smoothers.resize(min, max);
882
883 for (unsigned int i = min; i <= max; ++i)
884 {
885 // Workaround: Unfortunately, not every "m[i]" object has a rich
886 // enough interface to populate reinit_(domain|range)_vector. Thus,
887 // apply an empty LinearOperator exemplar.
888 matrices[i] = linear_operator<VectorType>(LinearOperator<VectorType>(),
889 m[i].block(row, col));
890 smoothers[i].initialize(m[i].block(row, col), data);
891 }
892}
893
894template <typename MatrixType, typename RelaxationType, typename VectorType>
895template <typename MatrixType2, typename DataType>
896inline void
897MGSmootherRelaxation<MatrixType, RelaxationType, VectorType>::initialize(
898 const MGLevelObject<MatrixType2> &m,
899 const MGLevelObject<DataType> &data,
900 const unsigned int row,
901 const unsigned int col)
902{
903 const unsigned int min = m.min_level();
904 const unsigned int max = m.max_level();
905
906 Assert(data.min_level() == min, ExcDimensionMismatch(data.min_level(), min));
907 Assert(data.max_level() == max, ExcDimensionMismatch(data.max_level(), max));
908
909 matrices.resize(min, max);
910 smoothers.resize(min, max);
911
912 for (unsigned int i = min; i <= max; ++i)
913 {
914 // Workaround: Unfortunately, not every "m[i]" object has a rich
915 // enough interface to populate reinit_(domain|range)_vector. Thus,
916 // apply an empty LinearOperator exemplar.
917 matrices[i] = linear_operator<VectorType>(LinearOperator<VectorType>(),
918 m[i].block(row, col));
919 smoothers[i].initialize(m[i].block(row, col), data[i]);
920 }
921}
922
923
924template <typename MatrixType, typename RelaxationType, typename VectorType>
925inline void
926MGSmootherRelaxation<MatrixType, RelaxationType, VectorType>::smooth(
927 const unsigned int level,
928 VectorType &u,
929 const VectorType &rhs) const
930{
931 unsigned int maxlevel = smoothers.max_level();
932 unsigned int steps2 = this->steps;
933
934 if (this->variable)
935 steps2 *= (1 << (maxlevel - level));
936
937 bool T = this->transpose;
938 if (this->symmetric && (steps2 % 2 == 0))
939 T = false;
940 if (this->debug > 0)
941 deallog << 'S' << level << ' ';
942
943 for (unsigned int i = 0; i < steps2; ++i)
944 {
945 if (T)
946 smoothers[level].Tstep(u, rhs);
947 else
948 smoothers[level].step(u, rhs);
949 if (this->symmetric)
950 T = !T;
951 }
952}
953
954
955template <typename MatrixType, typename RelaxationType, typename VectorType>
956inline void
957MGSmootherRelaxation<MatrixType, RelaxationType, VectorType>::apply(
958 const unsigned int level,
959 VectorType &u,
960 const VectorType &rhs) const
961{
962 unsigned int maxlevel = smoothers.max_level();
963 unsigned int steps2 = this->steps;
964
965 if (this->variable)
966 steps2 *= (1 << (maxlevel - level));
967
968 bool T = this->transpose;
969 if (this->symmetric && (steps2 % 2 == 0))
970 T = false;
971 if (this->debug > 0)
972 deallog << 'S' << level << ' ';
973
974 if (T)
975 smoothers[level].Tvmult(u, rhs);
976 else
977 smoothers[level].vmult(u, rhs);
978 if (this->symmetric)
979 T = !T;
980 for (unsigned int i = 1; i < steps2; ++i)
981 {
982 if (T)
983 smoothers[level].Tstep(u, rhs);
984 else
985 smoothers[level].step(u, rhs);
986 if (this->symmetric)
987 T = !T;
988 }
989}
990
991
992
993template <typename MatrixType, typename RelaxationType, typename VectorType>
994inline std::size_t
995MGSmootherRelaxation<MatrixType, RelaxationType, VectorType>::
996 memory_consumption() const
997{
998 return sizeof(*this) + matrices.memory_consumption() +
999 smoothers.memory_consumption() +
1000 this->vector_memory.memory_consumption();
1001}
1002
1003
1004//----------------------------------------------------------------------//
1005
1006template <typename MatrixType, typename PreconditionerType, typename VectorType>
1007inline MGSmootherPrecondition<MatrixType, PreconditionerType, VectorType>::
1008 MGSmootherPrecondition(const unsigned int steps,
1009 const bool variable,
1010 const bool symmetric,
1011 const bool transpose)
1012 : MGSmoother<VectorType>(steps, variable, symmetric, transpose)
1013{}
1014
1015
1016
1017template <typename MatrixType, typename PreconditionerType, typename VectorType>
1018inline void
1019MGSmootherPrecondition<MatrixType, PreconditionerType, VectorType>::clear()
1020{
1021 smoothers.clear_elements();
1022
1023 unsigned int i = matrices.min_level(), max_level = matrices.max_level();
1024 for (; i <= max_level; ++i)
1025 matrices[i] = LinearOperator<VectorType>();
1026}
1027
1028
1029
1030template <typename MatrixType, typename PreconditionerType, typename VectorType>
1031template <typename MatrixType2>
1032inline void
1033MGSmootherPrecondition<MatrixType, PreconditionerType, VectorType>::initialize(
1034 const MGLevelObject<MatrixType2> &m,
1035 const typename PreconditionerType::AdditionalData &data)
1036{
1037 const unsigned int min = m.min_level();
1038 const unsigned int max = m.max_level();
1039
1040 matrices.resize(min, max);
1041 smoothers.resize(min, max);
1042
1043 for (unsigned int i = min; i <= max; ++i)
1044 {
1045 // Workaround: Unfortunately, not every "m[i]" object has a rich
1046 // enough interface to populate reinit_(domain|range)_vector. Thus,
1047 // apply an empty LinearOperator exemplar.
1048 matrices[i] =
1049 linear_operator<VectorType>(LinearOperator<VectorType>(),
1050 Utilities::get_underlying_value(m[i]));
1051 smoothers[i].initialize(Utilities::get_underlying_value(m[i]), data);
1052 }
1053}
1054
1055
1056
1057template <typename MatrixType, typename PreconditionerType, typename VectorType>
1058template <typename MatrixType2>
1059inline void
1060MGSmootherPrecondition<MatrixType, PreconditionerType, VectorType>::
1061 initialize_matrices(const MGLevelObject<MatrixType2> &m)
1062{
1063 const unsigned int min = m.min_level();
1064 const unsigned int max = m.max_level();
1065
1066 matrices.resize(min, max);
1067 smoothers.resize(min, max);
1068
1069 for (unsigned int i = min; i <= max; ++i)
1070 {
1071 // Workaround: Unfortunately, not every "m[i]" object has a rich
1072 // enough interface to populate reinit_(domain|range)_vector. Thus,
1073 // apply an empty LinearOperator exemplar.
1074 matrices[i] =
1075 linear_operator<VectorType>(LinearOperator<VectorType>(),
1076 Utilities::get_underlying_value(m[i]));
1077 }
1078}
1079
1080
1081
1082template <typename MatrixType, typename PreconditionerType, typename VectorType>
1083template <typename MatrixType2, typename DataType>
1084inline void
1085MGSmootherPrecondition<MatrixType, PreconditionerType, VectorType>::initialize(
1086 const MGLevelObject<MatrixType2> &m,
1087 const MGLevelObject<DataType> &data)
1088{
1089 const unsigned int min = m.min_level();
1090 const unsigned int max = m.max_level();
1091
1092 Assert(data.min_level() == min, ExcDimensionMismatch(data.min_level(), min));
1093 Assert(data.max_level() == max, ExcDimensionMismatch(data.max_level(), max));
1094
1095 matrices.resize(min, max);
1096 smoothers.resize(min, max);
1097
1098 for (unsigned int i = min; i <= max; ++i)
1099 {
1100 // Workaround: Unfortunately, not every "m[i]" object has a rich
1101 // enough interface to populate reinit_(domain|range)_vector. Thus,
1102 // apply an empty LinearOperator exemplar.
1103 matrices[i] =
1104 linear_operator<VectorType>(LinearOperator<VectorType>(),
1105 Utilities::get_underlying_value(m[i]));
1106 smoothers[i].initialize(Utilities::get_underlying_value(m[i]), data[i]);
1107 }
1108}
1109
1110
1111
1112template <typename MatrixType, typename PreconditionerType, typename VectorType>
1113template <typename MatrixType2, typename DataType>
1114inline void
1115MGSmootherPrecondition<MatrixType, PreconditionerType, VectorType>::initialize(
1116 const MGLevelObject<MatrixType2> &m,
1117 const DataType &data,
1118 const unsigned int row,
1119 const unsigned int col)
1120{
1121 const unsigned int min = m.min_level();
1122 const unsigned int max = m.max_level();
1123
1124 matrices.resize(min, max);
1125 smoothers.resize(min, max);
1126
1127 for (unsigned int i = min; i <= max; ++i)
1128 {
1129 matrices[i] = &(m[i].block(row, col));
1130 smoothers[i].initialize(m[i].block(row, col), data);
1131 }
1132}
1133
1134
1135
1136template <typename MatrixType, typename PreconditionerType, typename VectorType>
1137template <typename MatrixType2, typename DataType>
1138inline void
1139MGSmootherPrecondition<MatrixType, PreconditionerType, VectorType>::initialize(
1140 const MGLevelObject<MatrixType2> &m,
1141 const MGLevelObject<DataType> &data,
1142 const unsigned int row,
1143 const unsigned int col)
1144{
1145 const unsigned int min = m.min_level();
1146 const unsigned int max = m.max_level();
1147
1148 Assert(data.min_level() == min, ExcDimensionMismatch(data.min_level(), min));
1149 Assert(data.max_level() == max, ExcDimensionMismatch(data.max_level(), max));
1150
1151 matrices.resize(min, max);
1152 smoothers.resize(min, max);
1153
1154 for (unsigned int i = min; i <= max; ++i)
1155 {
1156 matrices[i] = &(m[i].block(row, col));
1157 smoothers[i].initialize(m[i].block(row, col), data[i]);
1158 }
1159}
1160
1161
1162
1163template <typename MatrixType, typename PreconditionerType, typename VectorType>
1164inline void
1165MGSmootherPrecondition<MatrixType, PreconditionerType, VectorType>::smooth(
1166 const unsigned int level,
1167 VectorType &u,
1168 const VectorType &rhs) const
1169{
1170 unsigned int maxlevel = matrices.max_level();
1171 unsigned int steps2 = this->steps;
1172
1173 if (this->variable)
1174 steps2 *= (1 << (maxlevel - level));
1175
1176 typename VectorMemory<VectorType>::Pointer r(this->vector_memory);
1177 typename VectorMemory<VectorType>::Pointer d(this->vector_memory);
1178
1179 r->reinit(u, true);
1180 d->reinit(u, true);
1181
1182 bool T = this->transpose;
1183 if (this->symmetric && (steps2 % 2 == 0))
1184 T = false;
1185 if (this->debug > 0)
1186 deallog << 'S' << level << ' ';
1187
1188 for (unsigned int i = 0; i < steps2; ++i)
1189 {
1190 if (T)
1191 {
1192 if (this->debug > 0)
1193 deallog << 'T';
1194 matrices[level].Tvmult(*r, u);
1195 r->sadd(-1., 1., rhs);
1196 if (this->debug > 2)
1197 deallog << ' ' << r->l2_norm() << ' ';
1198 smoothers[level].Tvmult(*d, *r);
1199 if (this->debug > 1)
1200 deallog << ' ' << d->l2_norm() << ' ';
1201 }
1202 else
1203 {
1204 if (this->debug > 0)
1205 deallog << 'N';
1206 matrices[level].vmult(*r, u);
1207 r->sadd(-1., rhs);
1208 if (this->debug > 2)
1209 deallog << ' ' << r->l2_norm() << ' ';
1210 smoothers[level].vmult(*d, *r);
1211 if (this->debug > 1)
1212 deallog << ' ' << d->l2_norm() << ' ';
1213 }
1214 u += *d;
1215 if (this->symmetric)
1216 T = !T;
1217 }
1218 if (this->debug > 0)
1219 deallog << std::endl;
1220}
1221
1222
1223
1224template <typename MatrixType, typename PreconditionerType, typename VectorType>
1225inline void
1226MGSmootherPrecondition<MatrixType, PreconditionerType, VectorType>::apply(
1227 const unsigned int level,
1228 VectorType &u,
1229 const VectorType &rhs) const
1230{
1231 unsigned int maxlevel = matrices.max_level();
1232 unsigned int steps2 = this->steps;
1233
1234 if (this->variable)
1235 steps2 *= (1 << (maxlevel - level));
1236
1237 bool T = this->transpose;
1238 if (this->symmetric && (steps2 % 2 == 0))
1239 T = false;
1240 if (this->debug > 0)
1241 deallog << 'S' << level << ' ';
1242
1243 // first step where we overwrite the result
1244 if (this->debug > 2)
1245 deallog << ' ' << rhs.l2_norm() << ' ';
1246 if (this->debug > 0)
1247 deallog << (T ? 'T' : 'N');
1248 if (T)
1249 smoothers[level].Tvmult(u, rhs);
1250 else
1251 smoothers[level].vmult(u, rhs);
1252 if (this->debug > 1)
1253 deallog << ' ' << u.l2_norm() << ' ';
1254 if (this->symmetric)
1255 T = !T;
1256
1257 typename VectorMemory<VectorType>::Pointer r(this->vector_memory);
1258 typename VectorMemory<VectorType>::Pointer d(this->vector_memory);
1259
1260 if (steps2 > 1)
1261 {
1262 r->reinit(u, true);
1263 d->reinit(u, true);
1264 }
1265
1266 for (unsigned int i = 1; i < steps2; ++i)
1267 {
1268 if (T)
1269 {
1270 if (this->debug > 0)
1271 deallog << 'T';
1272 matrices[level].Tvmult(*r, u);
1273 r->sadd(-1., 1., rhs);
1274 if (this->debug > 2)
1275 deallog << ' ' << r->l2_norm() << ' ';
1276 smoothers[level].Tvmult(*d, *r);
1277 if (this->debug > 1)
1278 deallog << ' ' << d->l2_norm() << ' ';
1279 }
1280 else
1281 {
1282 if (this->debug > 0)
1283 deallog << 'N';
1284 matrices[level].vmult(*r, u);
1285 r->sadd(-1., rhs);
1286 if (this->debug > 2)
1287 deallog << ' ' << r->l2_norm() << ' ';
1288 smoothers[level].vmult(*d, *r);
1289 if (this->debug > 1)
1290 deallog << ' ' << d->l2_norm() << ' ';
1291 }
1292 u += *d;
1293 if (this->symmetric)
1294 T = !T;
1295 }
1296 if (this->debug > 0)
1297 deallog << std::endl;
1298}
1299
1300
1301
1302template <typename MatrixType, typename PreconditionerType, typename VectorType>
1303inline std::size_t
1304MGSmootherPrecondition<MatrixType, PreconditionerType, VectorType>::
1305 memory_consumption() const
1306{
1307 return sizeof(*this) + matrices.memory_consumption() +
1308 smoothers.memory_consumption() +
1309 this->vector_memory.memory_consumption();
1310}
1311
1312
1313#endif // DOXYGEN
1314
1315DEAL_II_NAMESPACE_CLOSE
1316
1317#endif
LogStream::memory_consumption
std::size_t memory_consumption() const
MGLevelObject::memory_consumption
std::size_t memory_consumption() const
Definition mg_level_object.h:338
MGLevelObject::max_level
unsigned int max_level() const
Definition mg_level_object.h:311
MGLevelObject::min_level
unsigned int min_level() const
Definition mg_level_object.h:303
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
void initialize(const MGLevelObject< MatrixType2 > &matrices, const MGLevelObject< DataType > &additional_data)
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 DataType &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< DataType > &additional_data, const unsigned int block_row, const unsigned int block_col)
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 DataType &additional_data, const unsigned int block_row, const unsigned int block_col)
MGSmootherRelaxation::initialize
void initialize(const MGLevelObject< MatrixType2 > &matrices, const MGLevelObject< DataType > &additional_data)
MGSmootherRelaxation::initialize
void initialize(const MGLevelObject< MatrixType2 > &matrices, const MGLevelObject< DataType > &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::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::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::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::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::initialize
void initialize(const MGLevelObject< MatrixType2 > &matrices, const MGLevelObject< DataType > &additional_data)
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::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:38
DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:39
transpose
DerivativeForm< 1, spacedim, dim, Number > transpose(const DerivativeForm< 1, dim, spacedim, Number > &DF)
Definition derivative_form.h:629
level
unsigned int level
Definition grid_out.cc:4642
Assert
#define Assert(cond, exc)
Definition exception_macros.h:543
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:36
mg_base.h
mg_level_object.h
data
std::vector< index_type > data
Definition mpi.cc:734
LAPACKSupport::symmetric
@ symmetric
Matrix is symmetric.
Definition lapack_support.h:111
LAPACKSupport::T
constexpr char T
Definition lapack_support.h:163
LinearAlgebra::TpetraWrappers::TpetraTypes::VectorType
Tpetra::Vector< Number, LO, GO, NodeType< MemorySpace > > VectorType
Definition trilinos_tpetra_types.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:1610
mg
Definition mg.h:79
std::min
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
Definition vectorization.h:6596
std::max
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
Definition vectorization.h:6579
observer_pointer.h
Generated by doxygen 1.9.8