Reference documentation for deal.II version 9.5.0
\(\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
49template <typename VectorType>
50class MGSmoother : public MGSmootherBase<VectorType>
51{
52public:
56 MGSmoother(const unsigned int steps = 1,
57 const bool variable = false,
58 const bool symmetric = false,
59 const bool transpose = false);
60
64 void
65 set_steps(const unsigned int);
66
70 void
71 set_variable(const bool);
72
76 void
77 set_symmetric(const bool);
78
83 void
84 set_transpose(const bool);
85
90 void
91 set_debug(const unsigned int level);
92
93protected:
100 mutable GrowingVectorMemory<VectorType> vector_memory;
101
106 unsigned int steps;
107
112 bool variable;
113
118 bool symmetric;
119
124 bool transpose;
125
129 unsigned int debug;
130};
131
132
139template <typename VectorType>
140class MGSmootherIdentity : public MGSmootherBase<VectorType>
141{
142public:
148 virtual void
149 smooth(const unsigned int level, VectorType &u, const VectorType &rhs) const;
150
151 virtual void
152 clear();
153};
154
155
156namespace mg
157{
184 template <class RelaxationType, typename VectorType>
185 class SmootherRelaxation : public MGLevelObject<RelaxationType>,
186 public MGSmoother<VectorType>
187 {
188 public:
192 SmootherRelaxation(const unsigned int steps = 1,
193 const bool variable = false,
194 const bool symmetric = false,
195 const bool transpose = false);
196
205 template <typename MatrixType2>
206 void
207 initialize(const MGLevelObject<MatrixType2> & matrices,
208 const typename RelaxationType::AdditionalData &additional_data =
209 typename RelaxationType::AdditionalData());
210
218 template <typename MatrixType2, class DATA>
219 void
220 initialize(const MGLevelObject<MatrixType2> &matrices,
221 const MGLevelObject<DATA> & additional_data);
222
226 void
227 clear() override;
228
232 virtual void
233 smooth(const unsigned int level,
234 VectorType & u,
235 const VectorType & rhs) const override;
236
253 virtual void
254 apply(const unsigned int level,
255 VectorType & u,
256 const VectorType & rhs) const override;
257
261 std::size_t
262 memory_consumption() const;
263 };
264} // namespace mg
265
295template <typename MatrixType, class RelaxationType, typename VectorType>
296class MGSmootherRelaxation : public MGSmoother<VectorType>
297{
298public:
302 MGSmootherRelaxation(const unsigned int steps = 1,
303 const bool variable = false,
304 const bool symmetric = false,
305 const bool transpose = false);
306
315 template <typename MatrixType2>
316 void
317 initialize(const MGLevelObject<MatrixType2> & matrices,
318 const typename RelaxationType::AdditionalData &additional_data =
319 typename RelaxationType::AdditionalData());
320
329 template <typename MatrixType2, class DATA>
330 void
331 initialize(const MGLevelObject<MatrixType2> &matrices,
332 const MGLevelObject<DATA> & additional_data);
333
343 template <typename MatrixType2, class DATA>
344 void
345 initialize(const MGLevelObject<MatrixType2> &matrices,
346 const DATA & additional_data,
347 const unsigned int block_row,
348 const unsigned int block_col);
349
359 template <typename MatrixType2, class DATA>
360 void
361 initialize(const MGLevelObject<MatrixType2> &matrices,
362 const MGLevelObject<DATA> & additional_data,
363 const unsigned int block_row,
364 const unsigned int block_col);
365
369 void
370 clear();
371
375 virtual void
376 smooth(const unsigned int level, VectorType &u, const VectorType &rhs) const;
377
393 virtual void
394 apply(const unsigned int level, VectorType &u, const VectorType &rhs) const;
395
399 MGLevelObject<RelaxationType> smoothers;
400
404 std::size_t
405 memory_consumption() const;
406
407
408private:
412 MGLevelObject<LinearOperator<VectorType>> matrices;
413};
414
415
416
445template <typename MatrixType, typename PreconditionerType, typename VectorType>
446class MGSmootherPrecondition : public MGSmoother<VectorType>
447{
448public:
452 MGSmootherPrecondition(const unsigned int steps = 1,
453 const bool variable = false,
454 const bool symmetric = false,
455 const bool transpose = false);
456
466 template <typename MatrixType2>
467 void
468 initialize(const MGLevelObject<MatrixType2> &matrices,
469 const typename PreconditionerType::AdditionalData &
470 additional_data = typename PreconditionerType::AdditionalData());
471
479 template <typename MatrixType2>
480 void
481 initialize_matrices(const MGLevelObject<MatrixType2> &matrices);
482
492 template <typename MatrixType2, class DATA>
493 void
494 initialize(const MGLevelObject<MatrixType2> &matrices,
495 const MGLevelObject<DATA> & additional_data);
496
507 template <typename MatrixType2, class DATA>
508 void
509 initialize(const MGLevelObject<MatrixType2> &matrices,
510 const DATA & additional_data,
511 const unsigned int block_row,
512 const unsigned int block_col);
513
524 template <typename MatrixType2, class DATA>
525 void
526 initialize(const MGLevelObject<MatrixType2> &matrices,
527 const MGLevelObject<DATA> & additional_data,
528 const unsigned int block_row,
529 const unsigned int block_col);
530
534 void
535 clear() override;
536
540 virtual void
541 smooth(const unsigned int level,
542 VectorType & u,
543 const VectorType & rhs) const override;
544
560 virtual void
561 apply(const unsigned int level,
562 VectorType & u,
563 const VectorType & rhs) const override;
564
568 MGLevelObject<PreconditionerType> smoothers;
569
573 std::size_t
574 memory_consumption() const;
575
576
577private:
581 MGLevelObject<LinearOperator<VectorType>> matrices;
582};
583
586/* ------------------------------- Inline functions --------------------------
587 */
588
589#ifndef DOXYGEN
590
591template <typename VectorType>
592inline void
593MGSmootherIdentity<VectorType>::smooth(const unsigned int,
594 VectorType &,
595 const VectorType &) const
596{}
597
598template <typename VectorType>
599inline void
600MGSmootherIdentity<VectorType>::clear()
601{}
602
603//---------------------------------------------------------------------------
604
605template <typename VectorType>
606inline MGSmoother<VectorType>::MGSmoother(const unsigned int steps,
607 const bool variable,
608 const bool symmetric,
609 const bool transpose)
610 : steps(steps)
611 , variable(variable)
612 , symmetric(symmetric)
613 , transpose(transpose)
614 , debug(0)
615{}
616
617
618template <typename VectorType>
619inline void
620MGSmoother<VectorType>::set_steps(const unsigned int s)
621{
622 steps = s;
623}
624
625
626template <typename VectorType>
627inline void
628MGSmoother<VectorType>::set_debug(const unsigned int s)
629{
630 debug = s;
631}
632
633
634template <typename VectorType>
635inline void
636MGSmoother<VectorType>::set_variable(const bool flag)
637{
638 variable = flag;
639}
640
641
642template <typename VectorType>
643inline void
644MGSmoother<VectorType>::set_symmetric(const bool flag)
645{
646 symmetric = flag;
647}
648
649
650template <typename VectorType>
651inline void
652MGSmoother<VectorType>::set_transpose(const bool flag)
653{
654 transpose = flag;
655}
656
657//----------------------------------------------------------------------//
658
659namespace mg
660{
661 template <class RelaxationType, typename VectorType>
662 inline SmootherRelaxation<RelaxationType, VectorType>::SmootherRelaxation(
663 const unsigned int steps,
664 const bool variable,
665 const bool symmetric,
666 const bool transpose)
667 : MGSmoother<VectorType>(steps, variable, symmetric, transpose)
668 {}
669
670
671 template <class RelaxationType, typename VectorType>
672 inline void
673 SmootherRelaxation<RelaxationType, VectorType>::clear()
674 {
675 MGLevelObject<RelaxationType>::clear_elements();
676 }
677
678
679 template <class RelaxationType, typename VectorType>
680 template <typename MatrixType2>
681 inline void
682 SmootherRelaxation<RelaxationType, VectorType>::initialize(
683 const MGLevelObject<MatrixType2> & m,
684 const typename RelaxationType::AdditionalData &data)
685 {
686 const unsigned int min = m.min_level();
687 const unsigned int max = m.max_level();
688
689 this->resize(min, max);
690
691 for (unsigned int i = min; i <= max; ++i)
692 (*this)[i].initialize(m[i], data);
693 }
694
695
696 template <class RelaxationType, typename VectorType>
697 template <typename MatrixType2, class DATA>
698 inline void
699 SmootherRelaxation<RelaxationType, VectorType>::initialize(
700 const MGLevelObject<MatrixType2> &m,
701 const MGLevelObject<DATA> & data)
702 {
703 const unsigned int min = std::max(m.min_level(), data.min_level());
704 const unsigned int max = std::min(m.max_level(), data.max_level());
705
706 this->resize(min, max);
707
708 for (unsigned int i = min; i <= max; ++i)
709 (*this)[i].initialize(m[i], data[i]);
710 }
711
712
713 template <class RelaxationType, typename VectorType>
714 inline void
715 SmootherRelaxation<RelaxationType, VectorType>::smooth(
716 const unsigned int level,
717 VectorType & u,
718 const VectorType & rhs) const
719 {
720 unsigned int maxlevel = this->max_level();
721 unsigned int steps2 = this->steps;
722
723 if (this->variable)
724 steps2 *= (1 << (maxlevel - level));
725
726 bool T = this->transpose;
727 if (this->symmetric && (steps2 % 2 == 0))
728 T = false;
729 if (this->debug > 0)
730 deallog << 'S' << level << ' ';
731
732 for (unsigned int i = 0; i < steps2; ++i)
733 {
734 if (T)
735 (*this)[level].Tstep(u, rhs);
736 else
737 (*this)[level].step(u, rhs);
738 if (this->symmetric)
739 T = !T;
740 }
741 }
742
743
744 template <class RelaxationType, typename VectorType>
745 inline void
746 SmootherRelaxation<RelaxationType, VectorType>::apply(
747 const unsigned int level,
748 VectorType & u,
749 const VectorType & rhs) const
750 {
751 unsigned int maxlevel = this->max_level();
752 unsigned int steps2 = this->steps;
753
754 if (this->variable)
755 steps2 *= (1 << (maxlevel - level));
756
757 bool T = this->transpose;
758 if (this->symmetric && (steps2 % 2 == 0))
759 T = false;
760 if (this->debug > 0)
761 deallog << 'S' << level << ' ';
762
763 if (T)
764 (*this)[level].Tvmult(u, rhs);
765 else
766 (*this)[level].vmult(u, rhs);
767 if (this->symmetric)
768 T = !T;
769 for (unsigned int i = 1; i < steps2; ++i)
770 {
771 if (T)
772 (*this)[level].Tstep(u, rhs);
773 else
774 (*this)[level].step(u, rhs);
775 if (this->symmetric)
776 T = !T;
777 }
778 }
779
780
781 template <class RelaxationType, typename VectorType>
782 inline std::size_t
783 SmootherRelaxation<RelaxationType, VectorType>::memory_consumption() const
784 {
785 return sizeof(*this) - sizeof(MGLevelObject<RelaxationType>) +
786 MGLevelObject<RelaxationType>::memory_consumption() +
787 this->vector_memory.memory_consumption();
788 }
789} // namespace mg
790
791
792//----------------------------------------------------------------------//
793
794template <typename MatrixType, class RelaxationType, typename VectorType>
795inline MGSmootherRelaxation<MatrixType, RelaxationType, VectorType>::
796 MGSmootherRelaxation(const unsigned int steps,
797 const bool variable,
798 const bool symmetric,
799 const bool transpose)
800 : MGSmoother<VectorType>(steps, variable, symmetric, transpose)
801{}
802
803
804
805template <typename MatrixType, class RelaxationType, typename VectorType>
806inline void
807MGSmootherRelaxation<MatrixType, RelaxationType, VectorType>::clear()
808{
809 smoothers.clear_elements();
810
811 unsigned int i = matrices.min_level(), max_level = matrices.max_level();
812 for (; i <= max_level; ++i)
813 matrices[i] = LinearOperator<VectorType>();
814}
815
816
817template <typename MatrixType, class RelaxationType, typename VectorType>
818template <typename MatrixType2>
819inline void
820MGSmootherRelaxation<MatrixType, RelaxationType, VectorType>::initialize(
821 const MGLevelObject<MatrixType2> & m,
822 const typename RelaxationType::AdditionalData &data)
823{
824 const unsigned int min = m.min_level();
825 const unsigned int max = m.max_level();
826
827 matrices.resize(min, max);
828 smoothers.resize(min, max);
829
830 for (unsigned int i = min; i <= max; ++i)
831 {
832 // Workaround: Unfortunately, not every "m[i]" object has a rich
833 // enough interface to populate reinit_(domain|range)_vector. Thus,
834 // apply an empty LinearOperator exemplar.
835 matrices[i] =
836 linear_operator<VectorType>(LinearOperator<VectorType>(), m[i]);
837 smoothers[i].initialize(m[i], data);
838 }
839}
840
841template <typename MatrixType, class RelaxationType, typename VectorType>
842template <typename MatrixType2, class DATA>
843inline void
844MGSmootherRelaxation<MatrixType, RelaxationType, VectorType>::initialize(
845 const MGLevelObject<MatrixType2> &m,
846 const MGLevelObject<DATA> & data)
847{
848 const unsigned int min = m.min_level();
849 const unsigned int max = m.max_level();
850
851 Assert(data.min_level() == min, ExcDimensionMismatch(data.min_level(), min));
852 Assert(data.max_level() == max, ExcDimensionMismatch(data.max_level(), max));
853
854 matrices.resize(min, max);
855 smoothers.resize(min, max);
856
857 for (unsigned int i = min; i <= max; ++i)
858 {
859 // Workaround: Unfortunately, not every "m[i]" object has a rich
860 // enough interface to populate reinit_(domain|range)_vector. Thus,
861 // apply an empty LinearOperator exemplar.
862 matrices[i] =
863 linear_operator<VectorType>(LinearOperator<VectorType>(), m[i]);
864 smoothers[i].initialize(m[i], data[i]);
865 }
866}
867
868template <typename MatrixType, class RelaxationType, typename VectorType>
869template <typename MatrixType2, class DATA>
870inline void
871MGSmootherRelaxation<MatrixType, RelaxationType, VectorType>::initialize(
872 const MGLevelObject<MatrixType2> &m,
873 const DATA & data,
874 const unsigned int row,
875 const unsigned int col)
876{
877 const unsigned int min = m.min_level();
878 const unsigned int max = m.max_level();
879
880 matrices.resize(min, max);
881 smoothers.resize(min, max);
882
883 for (unsigned int i = min; i <= max; ++i)
884 {
885 // Workaround: Unfortunately, not every "m[i]" object has a rich
886 // enough interface to populate reinit_(domain|range)_vector. Thus,
887 // apply an empty LinearOperator exemplar.
888 matrices[i] = linear_operator<VectorType>(LinearOperator<VectorType>(),
889 m[i].block(row, col));
890 smoothers[i].initialize(m[i].block(row, col), data);
891 }
892}
893
894template <typename MatrixType, class RelaxationType, typename VectorType>
895template <typename MatrixType2, class DATA>
896inline void
897MGSmootherRelaxation<MatrixType, RelaxationType, VectorType>::initialize(
898 const MGLevelObject<MatrixType2> &m,
899 const MGLevelObject<DATA> & data,
900 const unsigned int row,
901 const unsigned int col)
902{
903 const unsigned int min = m.min_level();
904 const unsigned int max = m.max_level();
905
906 Assert(data.min_level() == min, ExcDimensionMismatch(data.min_level(), min));
907 Assert(data.max_level() == max, ExcDimensionMismatch(data.max_level(), max));
908
909 matrices.resize(min, max);
910 smoothers.resize(min, max);
911
912 for (unsigned int i = min; i <= max; ++i)
913 {
914 // Workaround: Unfortunately, not every "m[i]" object has a rich
915 // enough interface to populate reinit_(domain|range)_vector. Thus,
916 // apply an empty LinearOperator exemplar.
917 matrices[i] = linear_operator<VectorType>(LinearOperator<VectorType>(),
918 m[i].block(row, col));
919 smoothers[i].initialize(m[i].block(row, col), data[i]);
920 }
921}
922
923
924template <typename MatrixType, class 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, class 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, class 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, class DATA>
1084inline void
1085MGSmootherPrecondition<MatrixType, PreconditionerType, VectorType>::initialize(
1086 const MGLevelObject<MatrixType2> &m,
1087 const MGLevelObject<DATA> & data)
1088{
1089 const unsigned int min = m.min_level();
1090 const unsigned int max = m.max_level();
1091
1092 Assert(data.min_level() == min, ExcDimensionMismatch(data.min_level(), min));
1093 Assert(data.max_level() == max, ExcDimensionMismatch(data.max_level(), max));
1094
1095 matrices.resize(min, max);
1096 smoothers.resize(min, max);
1097
1098 for (unsigned int i = min; i <= max; ++i)
1099 {
1100 // Workaround: Unfortunately, not every "m[i]" object has a rich
1101 // enough interface to populate reinit_(domain|range)_vector. Thus,
1102 // apply an empty LinearOperator exemplar.
1103 matrices[i] =
1104 linear_operator<VectorType>(LinearOperator<VectorType>(),
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, class DATA>
1114inline void
1115MGSmootherPrecondition<MatrixType, PreconditionerType, VectorType>::initialize(
1116 const MGLevelObject<MatrixType2> &m,
1117 const DATA & data,
1118 const unsigned int row,
1119 const unsigned int col)
1120{
1121 const unsigned int min = m.min_level();
1122 const unsigned int max = m.max_level();
1123
1124 matrices.resize(min, max);
1125 smoothers.resize(min, max);
1126
1127 for (unsigned int i = min; i <= max; ++i)
1128 {
1129 matrices[i] = &(m[i].block(row, col));
1130 smoothers[i].initialize(m[i].block(row, col), data);
1131 }
1132}
1133
1134
1135
1136template <typename MatrixType, typename PreconditionerType, typename VectorType>
1137template <typename MatrixType2, class DATA>
1138inline void
1139MGSmootherPrecondition<MatrixType, PreconditionerType, VectorType>::initialize(
1140 const MGLevelObject<MatrixType2> &m,
1141 const MGLevelObject<DATA> & data,
1142 const unsigned int row,
1143 const unsigned int col)
1144{
1145 const unsigned int min = m.min_level();
1146 const unsigned int max = m.max_level();
1147
1148 Assert(data.min_level() == min, ExcDimensionMismatch(data.min_level(), min));
1149 Assert(data.max_level() == max, ExcDimensionMismatch(data.max_level(), max));
1150
1151 matrices.resize(min, max);
1152 smoothers.resize(min, max);
1153
1154 for (unsigned int i = min; i <= max; ++i)
1155 {
1156 matrices[i] = &(m[i].block(row, col));
1157 smoothers[i].initialize(m[i].block(row, col), data[i]);
1158 }
1159}
1160
1161
1162
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
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:568
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:581
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:412
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:399
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::vector_memory
GrowingVectorMemory< VectorType > vector_memory
Definition mg_smoother.h:100
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:129
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:106
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:472
DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
transpose
DerivativeForm< 1, spacedim, dim, Number > transpose(const DerivativeForm< 1, dim, spacedim, Number > &DF)
Definition derivative_form.h:586
level
unsigned int level
Definition grid_out.cc:4618
Assert
#define Assert(cond, exc)
Definition exceptions.h:1614
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:1612
mg
Definition mg.h:82
std::min
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
Definition vectorization.h:6041
std::max
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
Definition vectorization.h:6024
Generated by doxygen 1.9.7