deal.II version GIT relicensing-2291-g668cd86249 2024-12-24 11:30:00+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: LGPL-2.1-or-later
4// Copyright (C) 1999 - 2023 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15#ifndef dealii_mg_smoother_h
16#define dealii_mg_smoother_h
17
18
19#include <deal.II/base/config.h>
20
21#include <deal.II/base/mg_level_object.h>
22#include <deal.II/base/observer_pointer.h>
23
24#include <deal.II/lac/linear_operator.h>
25#include <deal.II/lac/vector_memory.h>
26
27#include <deal.II/multigrid/mg_base.h>
28
29#include <vector>
30
31DEAL_II_NAMESPACE_OPEN
32
33/*
34 * MGSmootherBase is defined in mg_base.h
35 */
36
48template <typename VectorType>
49class MGSmoother : public MGSmootherBase<VectorType>
50{
51public:
55 MGSmoother(const unsigned int steps = 1,
56 const bool variable = false,
57 const bool symmetric = false,
58 const bool transpose = false);
59
63 void
64 set_steps(const unsigned int);
65
69 void
70 set_variable(const bool);
71
75 void
76 set_symmetric(const bool);
77
82 void
83 set_transpose(const bool);
84
89 void
90 set_debug(const unsigned int level);
91
92protected:
99 mutable GrowingVectorMemory<VectorType> vector_memory;
100
105 unsigned int steps;
106
111 bool variable;
112
117 bool symmetric;
118
123 bool transpose;
124
128 unsigned int debug;
129};
130
131
138template <typename VectorType>
139class MGSmootherIdentity : public MGSmootherBase<VectorType>
140{
141public:
147 virtual void
148 smooth(const unsigned int level, VectorType &u, const VectorType &rhs) const;
149
150 virtual void
151 clear();
152};
153
154
155namespace mg
156{
183 template <typename RelaxationType, typename VectorType>
184 class SmootherRelaxation : public MGLevelObject<RelaxationType>,
185 public MGSmoother<VectorType>
186 {
187 public:
191 SmootherRelaxation(const unsigned int steps = 1,
192 const bool variable = false,
193 const bool symmetric = false,
194 const bool transpose = false);
195
204 template <typename MatrixType2>
205 void
206 initialize(const MGLevelObject<MatrixType2> &matrices,
207 const typename RelaxationType::AdditionalData &additional_data =
208 typename RelaxationType::AdditionalData());
209
217 template <typename MatrixType2, typename DataType>
218 void
219 initialize(const MGLevelObject<MatrixType2> &matrices,
220 const MGLevelObject<DataType> &additional_data);
221
225 void
226 clear() override;
227
231 virtual void
232 smooth(const unsigned int level,
233 VectorType &u,
234 const VectorType &rhs) const override;
235
252 virtual void
253 apply(const unsigned int level,
254 VectorType &u,
255 const VectorType &rhs) const override;
256
260 std::size_t
261 memory_consumption() const;
262 };
263} // namespace mg
264
294template <typename MatrixType, typename RelaxationType, typename VectorType>
295class MGSmootherRelaxation : public MGSmoother<VectorType>
296{
297public:
301 MGSmootherRelaxation(const unsigned int steps = 1,
302 const bool variable = false,
303 const bool symmetric = false,
304 const bool transpose = false);
305
314 template <typename MatrixType2>
315 void
316 initialize(const MGLevelObject<MatrixType2> &matrices,
317 const typename RelaxationType::AdditionalData &additional_data =
318 typename RelaxationType::AdditionalData());
319
328 template <typename MatrixType2, typename DataType>
329 void
330 initialize(const MGLevelObject<MatrixType2> &matrices,
331 const MGLevelObject<DataType> &additional_data);
332
342 template <typename MatrixType2, typename DataType>
343 void
344 initialize(const MGLevelObject<MatrixType2> &matrices,
345 const DataType &additional_data,
346 const unsigned int block_row,
347 const unsigned int block_col);
348
358 template <typename MatrixType2, typename DataType>
359 void
360 initialize(const MGLevelObject<MatrixType2> &matrices,
361 const MGLevelObject<DataType> &additional_data,
362 const unsigned int block_row,
363 const unsigned int block_col);
364
368 void
369 clear();
370
374 virtual void
375 smooth(const unsigned int level, VectorType &u, const VectorType &rhs) const;
376
392 virtual void
393 apply(const unsigned int level, VectorType &u, const VectorType &rhs) const;
394
398 MGLevelObject<RelaxationType> smoothers;
399
403 std::size_t
404 memory_consumption() const;
405
406
407private:
411 MGLevelObject<LinearOperator<VectorType>> matrices;
412};
413
414
415
444template <typename MatrixType, typename PreconditionerType, typename VectorType>
445class MGSmootherPrecondition : public MGSmoother<VectorType>
446{
447public:
451 MGSmootherPrecondition(const unsigned int steps = 1,
452 const bool variable = false,
453 const bool symmetric = false,
454 const bool transpose = false);
455
465 template <typename MatrixType2>
466 void
467 initialize(const MGLevelObject<MatrixType2> &matrices,
468 const typename PreconditionerType::AdditionalData &
469 additional_data = typename PreconditionerType::AdditionalData());
470
478 template <typename MatrixType2>
479 void
480 initialize_matrices(const MGLevelObject<MatrixType2> &matrices);
481
491 template <typename MatrixType2, typename DataType>
492 void
493 initialize(const MGLevelObject<MatrixType2> &matrices,
494 const MGLevelObject<DataType> &additional_data);
495
506 template <typename MatrixType2, typename DataType>
507 void
508 initialize(const MGLevelObject<MatrixType2> &matrices,
509 const DataType &additional_data,
510 const unsigned int block_row,
511 const unsigned int block_col);
512
523 template <typename MatrixType2, typename DataType>
524 void
525 initialize(const MGLevelObject<MatrixType2> &matrices,
526 const MGLevelObject<DataType> &additional_data,
527 const unsigned int block_row,
528 const unsigned int block_col);
529
533 void
534 clear() override;
535
539 virtual void
540 smooth(const unsigned int level,
541 VectorType &u,
542 const VectorType &rhs) const override;
543
559 virtual void
560 apply(const unsigned int level,
561 VectorType &u,
562 const VectorType &rhs) const override;
563
567 MGLevelObject<PreconditionerType> smoothers;
568
572 std::size_t
573 memory_consumption() const;
574
575
576private:
580 MGLevelObject<LinearOperator<VectorType>> matrices;
581};
582
585/* ------------------------------- Inline functions --------------------------
586 */
587
588#ifndef DOXYGEN
589
590template <typename VectorType>
591inline void
592MGSmootherIdentity<VectorType>::smooth(const unsigned int,
593 VectorType &,
594 const VectorType &) const
595{}
596
597template <typename VectorType>
598inline void
599MGSmootherIdentity<VectorType>::clear()
600{}
601
602//---------------------------------------------------------------------------
603
604template <typename VectorType>
605inline MGSmoother<VectorType>::MGSmoother(const unsigned int steps,
606 const bool variable,
607 const bool symmetric,
608 const bool transpose)
609 : steps(steps)
610 , variable(variable)
611 , symmetric(symmetric)
612 , transpose(transpose)
613 , debug(0)
614{}
615
616
617template <typename VectorType>
618inline void
619MGSmoother<VectorType>::set_steps(const unsigned int s)
620{
621 steps = s;
622}
623
624
625template <typename VectorType>
626inline void
627MGSmoother<VectorType>::set_debug(const unsigned int s)
628{
629 debug = s;
630}
631
632
633template <typename VectorType>
634inline void
635MGSmoother<VectorType>::set_variable(const bool flag)
636{
637 variable = flag;
638}
639
640
641template <typename VectorType>
642inline void
643MGSmoother<VectorType>::set_symmetric(const bool flag)
644{
645 symmetric = flag;
646}
647
648
649template <typename VectorType>
650inline void
651MGSmoother<VectorType>::set_transpose(const bool flag)
652{
653 transpose = flag;
654}
655
656//----------------------------------------------------------------------//
657
658namespace mg
659{
660 template <typename RelaxationType, typename VectorType>
661 inline SmootherRelaxation<RelaxationType, VectorType>::SmootherRelaxation(
662 const unsigned int steps,
663 const bool variable,
664 const bool symmetric,
665 const bool transpose)
666 : MGSmoother<VectorType>(steps, variable, symmetric, transpose)
667 {}
668
669
670 template <typename RelaxationType, typename VectorType>
671 inline void
672 SmootherRelaxation<RelaxationType, VectorType>::clear()
673 {
674 MGLevelObject<RelaxationType>::clear_elements();
675 }
676
677
678 template <typename RelaxationType, typename VectorType>
679 template <typename MatrixType2>
680 inline void
681 SmootherRelaxation<RelaxationType, VectorType>::initialize(
682 const MGLevelObject<MatrixType2> &m,
683 const typename RelaxationType::AdditionalData &data)
684 {
685 const unsigned int min = m.min_level();
686 const unsigned int max = m.max_level();
687
688 this->resize(min, max);
689
690 for (unsigned int i = min; i <= max; ++i)
691 (*this)[i].initialize(Utilities::get_underlying_value(m[i]), data);
692 }
693
694
695 template <typename RelaxationType, typename VectorType>
696 template <typename MatrixType2, typename DataType>
697 inline void
698 SmootherRelaxation<RelaxationType, VectorType>::initialize(
699 const MGLevelObject<MatrixType2> &m,
700 const MGLevelObject<DataType> &data)
701 {
702 const unsigned int min = std::max(m.min_level(), data.min_level());
703 const unsigned int max = std::min(m.max_level(), data.max_level());
704
705 this->resize(min, max);
706
707 for (unsigned int i = min; i <= max; ++i)
708 (*this)[i].initialize(Utilities::get_underlying_value(m[i]), data[i]);
709 }
710
711
712 template <typename RelaxationType, typename VectorType>
713 inline void
714 SmootherRelaxation<RelaxationType, VectorType>::smooth(
715 const unsigned int level,
716 VectorType &u,
717 const VectorType &rhs) const
718 {
719 unsigned int maxlevel = this->max_level();
720 unsigned int steps2 = this->steps;
721
722 if (this->variable)
723 steps2 *= (1 << (maxlevel - level));
724
725 bool T = this->transpose;
726 if (this->symmetric && (steps2 % 2 == 0))
727 T = false;
728 if (this->debug > 0)
729 deallog << 'S' << level << ' ';
730
731 for (unsigned int i = 0; i < steps2; ++i)
732 {
733 if (T)
734 (*this)[level].Tstep(u, rhs);
735 else
736 (*this)[level].step(u, rhs);
737 if (this->symmetric)
738 T = !T;
739 }
740 }
741
742
743 template <typename RelaxationType, typename VectorType>
744 inline void
745 SmootherRelaxation<RelaxationType, VectorType>::apply(
746 const unsigned int level,
747 VectorType &u,
748 const VectorType &rhs) const
749 {
750 unsigned int maxlevel = this->max_level();
751 unsigned int steps2 = this->steps;
752
753 if (this->variable)
754 steps2 *= (1 << (maxlevel - level));
755
756 bool T = this->transpose;
757 if (this->symmetric && (steps2 % 2 == 0))
758 T = false;
759 if (this->debug > 0)
760 deallog << 'S' << level << ' ';
761
762 if (T)
763 (*this)[level].Tvmult(u, rhs);
764 else
765 (*this)[level].vmult(u, rhs);
766 if (this->symmetric)
767 T = !T;
768 for (unsigned int i = 1; i < steps2; ++i)
769 {
770 if (T)
771 (*this)[level].Tstep(u, rhs);
772 else
773 (*this)[level].step(u, rhs);
774 if (this->symmetric)
775 T = !T;
776 }
777 }
778
779
780 template <typename RelaxationType, typename VectorType>
781 inline std::size_t
782 SmootherRelaxation<RelaxationType, VectorType>::memory_consumption() const
783 {
784 return sizeof(*this) - sizeof(MGLevelObject<RelaxationType>) +
785 MGLevelObject<RelaxationType>::memory_consumption() +
786 this->vector_memory.memory_consumption();
787 }
788} // namespace mg
789
790
791//----------------------------------------------------------------------//
792
793template <typename MatrixType, typename RelaxationType, typename VectorType>
794inline MGSmootherRelaxation<MatrixType, RelaxationType, VectorType>::
795 MGSmootherRelaxation(const unsigned int steps,
796 const bool variable,
797 const bool symmetric,
798 const bool transpose)
799 : MGSmoother<VectorType>(steps, variable, symmetric, transpose)
800{}
801
802
803
804template <typename MatrixType, typename RelaxationType, typename VectorType>
805inline void
806MGSmootherRelaxation<MatrixType, RelaxationType, VectorType>::clear()
807{
808 smoothers.clear_elements();
809
810 unsigned int i = matrices.min_level(), max_level = matrices.max_level();
811 for (; i <= max_level; ++i)
812 matrices[i] = LinearOperator<VectorType>();
813}
814
815
816template <typename MatrixType, typename RelaxationType, typename VectorType>
817template <typename MatrixType2>
818inline void
819MGSmootherRelaxation<MatrixType, RelaxationType, VectorType>::initialize(
820 const MGLevelObject<MatrixType2> &m,
821 const typename RelaxationType::AdditionalData &data)
822{
823 const unsigned int min = m.min_level();
824 const unsigned int max = m.max_level();
825
826 matrices.resize(min, max);
827 smoothers.resize(min, max);
828
829 for (unsigned int i = min; i <= max; ++i)
830 {
831 // Workaround: Unfortunately, not every "m[i]" object has a rich
832 // enough interface to populate reinit_(domain|range)_vector. Thus,
833 // apply an empty LinearOperator exemplar.
834 matrices[i] =
835 linear_operator<VectorType>(LinearOperator<VectorType>(),
836 Utilities::get_underlying_value(m[i]));
837 smoothers[i].initialize(Utilities::get_underlying_value(m[i]), data);
838 }
839}
840
841template <typename MatrixType, typename RelaxationType, typename VectorType>
842template <typename MatrixType2, typename DataType>
843inline void
844MGSmootherRelaxation<MatrixType, RelaxationType, VectorType>::initialize(
845 const MGLevelObject<MatrixType2> &m,
846 const MGLevelObject<DataType> &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>(),
864 Utilities::get_underlying_value(m[i]));
865 smoothers[i].initialize(Utilities::get_underlying_value(m[i]), data[i]);
866 }
867}
868
869template <typename MatrixType, typename RelaxationType, typename VectorType>
870template <typename MatrixType2, typename DataType>
871inline void
872MGSmootherRelaxation<MatrixType, RelaxationType, VectorType>::initialize(
873 const MGLevelObject<MatrixType2> &m,
874 const DataType &data,
875 const unsigned int row,
876 const unsigned int col)
877{
878 const unsigned int min = m.min_level();
879 const unsigned int max = m.max_level();
880
881 matrices.resize(min, max);
882 smoothers.resize(min, max);
883
884 for (unsigned int i = min; i <= max; ++i)
885 {
886 // Workaround: Unfortunately, not every "m[i]" object has a rich
887 // enough interface to populate reinit_(domain|range)_vector. Thus,
888 // apply an empty LinearOperator exemplar.
889 matrices[i] = linear_operator<VectorType>(LinearOperator<VectorType>(),
890 m[i].block(row, col));
891 smoothers[i].initialize(m[i].block(row, col), data);
892 }
893}
894
895template <typename MatrixType, typename RelaxationType, typename VectorType>
896template <typename MatrixType2, typename DataType>
897inline void
898MGSmootherRelaxation<MatrixType, RelaxationType, VectorType>::initialize(
899 const MGLevelObject<MatrixType2> &m,
900 const MGLevelObject<DataType> &data,
901 const unsigned int row,
902 const unsigned int col)
903{
904 const unsigned int min = m.min_level();
905 const unsigned int max = m.max_level();
906
907 Assert(data.min_level() == min, ExcDimensionMismatch(data.min_level(), min));
908 Assert(data.max_level() == max, ExcDimensionMismatch(data.max_level(), max));
909
910 matrices.resize(min, max);
911 smoothers.resize(min, max);
912
913 for (unsigned int i = min; i <= max; ++i)
914 {
915 // Workaround: Unfortunately, not every "m[i]" object has a rich
916 // enough interface to populate reinit_(domain|range)_vector. Thus,
917 // apply an empty LinearOperator exemplar.
918 matrices[i] = linear_operator<VectorType>(LinearOperator<VectorType>(),
919 m[i].block(row, col));
920 smoothers[i].initialize(m[i].block(row, col), data[i]);
921 }
922}
923
924
925template <typename MatrixType, typename RelaxationType, typename VectorType>
926inline void
927MGSmootherRelaxation<MatrixType, RelaxationType, VectorType>::smooth(
928 const unsigned int level,
929 VectorType &u,
930 const VectorType &rhs) const
931{
932 unsigned int maxlevel = smoothers.max_level();
933 unsigned int steps2 = this->steps;
934
935 if (this->variable)
936 steps2 *= (1 << (maxlevel - level));
937
938 bool T = this->transpose;
939 if (this->symmetric && (steps2 % 2 == 0))
940 T = false;
941 if (this->debug > 0)
942 deallog << 'S' << level << ' ';
943
944 for (unsigned int i = 0; i < steps2; ++i)
945 {
946 if (T)
947 smoothers[level].Tstep(u, rhs);
948 else
949 smoothers[level].step(u, rhs);
950 if (this->symmetric)
951 T = !T;
952 }
953}
954
955
956template <typename MatrixType, typename RelaxationType, typename VectorType>
957inline void
958MGSmootherRelaxation<MatrixType, RelaxationType, VectorType>::apply(
959 const unsigned int level,
960 VectorType &u,
961 const VectorType &rhs) const
962{
963 unsigned int maxlevel = smoothers.max_level();
964 unsigned int steps2 = this->steps;
965
966 if (this->variable)
967 steps2 *= (1 << (maxlevel - level));
968
969 bool T = this->transpose;
970 if (this->symmetric && (steps2 % 2 == 0))
971 T = false;
972 if (this->debug > 0)
973 deallog << 'S' << level << ' ';
974
975 if (T)
976 smoothers[level].Tvmult(u, rhs);
977 else
978 smoothers[level].vmult(u, rhs);
979 if (this->symmetric)
980 T = !T;
981 for (unsigned int i = 1; i < steps2; ++i)
982 {
983 if (T)
984 smoothers[level].Tstep(u, rhs);
985 else
986 smoothers[level].step(u, rhs);
987 if (this->symmetric)
988 T = !T;
989 }
990}
991
992
993
994template <typename MatrixType, typename RelaxationType, typename VectorType>
995inline std::size_t
996MGSmootherRelaxation<MatrixType, RelaxationType, VectorType>::
997 memory_consumption() const
998{
999 return sizeof(*this) + matrices.memory_consumption() +
1000 smoothers.memory_consumption() +
1001 this->vector_memory.memory_consumption();
1002}
1003
1004
1005//----------------------------------------------------------------------//
1006
1007template <typename MatrixType, typename PreconditionerType, typename VectorType>
1008inline MGSmootherPrecondition<MatrixType, PreconditionerType, VectorType>::
1009 MGSmootherPrecondition(const unsigned int steps,
1010 const bool variable,
1011 const bool symmetric,
1012 const bool transpose)
1013 : MGSmoother<VectorType>(steps, variable, symmetric, transpose)
1014{}
1015
1016
1017
1018template <typename MatrixType, typename PreconditionerType, typename VectorType>
1019inline void
1020MGSmootherPrecondition<MatrixType, PreconditionerType, VectorType>::clear()
1021{
1022 smoothers.clear_elements();
1023
1024 unsigned int i = matrices.min_level(), max_level = matrices.max_level();
1025 for (; i <= max_level; ++i)
1026 matrices[i] = LinearOperator<VectorType>();
1027}
1028
1029
1030
1031template <typename MatrixType, typename PreconditionerType, typename VectorType>
1032template <typename MatrixType2>
1033inline void
1034MGSmootherPrecondition<MatrixType, PreconditionerType, VectorType>::initialize(
1035 const MGLevelObject<MatrixType2> &m,
1036 const typename PreconditionerType::AdditionalData &data)
1037{
1038 const unsigned int min = m.min_level();
1039 const unsigned int max = m.max_level();
1040
1041 matrices.resize(min, max);
1042 smoothers.resize(min, max);
1043
1044 for (unsigned int i = min; i <= max; ++i)
1045 {
1046 // Workaround: Unfortunately, not every "m[i]" object has a rich
1047 // enough interface to populate reinit_(domain|range)_vector. Thus,
1048 // apply an empty LinearOperator exemplar.
1049 matrices[i] =
1050 linear_operator<VectorType>(LinearOperator<VectorType>(),
1051 Utilities::get_underlying_value(m[i]));
1052 smoothers[i].initialize(Utilities::get_underlying_value(m[i]), data);
1053 }
1054}
1055
1056
1057
1058template <typename MatrixType, typename PreconditionerType, typename VectorType>
1059template <typename MatrixType2>
1060inline void
1061MGSmootherPrecondition<MatrixType, PreconditionerType, VectorType>::
1062 initialize_matrices(const MGLevelObject<MatrixType2> &m)
1063{
1064 const unsigned int min = m.min_level();
1065 const unsigned int max = m.max_level();
1066
1067 matrices.resize(min, max);
1068 smoothers.resize(min, max);
1069
1070 for (unsigned int i = min; i <= max; ++i)
1071 {
1072 // Workaround: Unfortunately, not every "m[i]" object has a rich
1073 // enough interface to populate reinit_(domain|range)_vector. Thus,
1074 // apply an empty LinearOperator exemplar.
1075 matrices[i] =
1076 linear_operator<VectorType>(LinearOperator<VectorType>(),
1077 Utilities::get_underlying_value(m[i]));
1078 }
1079}
1080
1081
1082
1083template <typename MatrixType, typename PreconditionerType, typename VectorType>
1084template <typename MatrixType2, typename DataType>
1085inline void
1086MGSmootherPrecondition<MatrixType, PreconditionerType, VectorType>::initialize(
1087 const MGLevelObject<MatrixType2> &m,
1088 const MGLevelObject<DataType> &data)
1089{
1090 const unsigned int min = m.min_level();
1091 const unsigned int max = m.max_level();
1092
1093 Assert(data.min_level() == min, ExcDimensionMismatch(data.min_level(), min));
1094 Assert(data.max_level() == max, ExcDimensionMismatch(data.max_level(), max));
1095
1096 matrices.resize(min, max);
1097 smoothers.resize(min, max);
1098
1099 for (unsigned int i = min; i <= max; ++i)
1100 {
1101 // Workaround: Unfortunately, not every "m[i]" object has a rich
1102 // enough interface to populate reinit_(domain|range)_vector. Thus,
1103 // apply an empty LinearOperator exemplar.
1104 matrices[i] =
1105 linear_operator<VectorType>(LinearOperator<VectorType>(),
1106 Utilities::get_underlying_value(m[i]));
1107 smoothers[i].initialize(Utilities::get_underlying_value(m[i]), data[i]);
1108 }
1109}
1110
1111
1112
1113template <typename MatrixType, typename PreconditionerType, typename VectorType>
1114template <typename MatrixType2, typename DataType>
1115inline void
1116MGSmootherPrecondition<MatrixType, PreconditionerType, VectorType>::initialize(
1117 const MGLevelObject<MatrixType2> &m,
1118 const DataType &data,
1119 const unsigned int row,
1120 const unsigned int col)
1121{
1122 const unsigned int min = m.min_level();
1123 const unsigned int max = m.max_level();
1124
1125 matrices.resize(min, max);
1126 smoothers.resize(min, max);
1127
1128 for (unsigned int i = min; i <= max; ++i)
1129 {
1130 matrices[i] = &(m[i].block(row, col));
1131 smoothers[i].initialize(m[i].block(row, col), data);
1132 }
1133}
1134
1135
1136
1137template <typename MatrixType, typename PreconditionerType, typename VectorType>
1138template <typename MatrixType2, typename DataType>
1139inline void
1140MGSmootherPrecondition<MatrixType, PreconditionerType, VectorType>::initialize(
1141 const MGLevelObject<MatrixType2> &m,
1142 const MGLevelObject<DataType> &data,
1143 const unsigned int row,
1144 const unsigned int col)
1145{
1146 const unsigned int min = m.min_level();
1147 const unsigned int max = m.max_level();
1148
1149 Assert(data.min_level() == min, ExcDimensionMismatch(data.min_level(), min));
1150 Assert(data.max_level() == max, ExcDimensionMismatch(data.max_level(), max));
1151
1152 matrices.resize(min, max);
1153 smoothers.resize(min, max);
1154
1155 for (unsigned int i = min; i <= max; ++i)
1156 {
1157 matrices[i] = &(m[i].block(row, col));
1158 smoothers[i].initialize(m[i].block(row, col), data[i]);
1159 }
1160}
1161
1162
1163
1164template <typename MatrixType, typename PreconditionerType, typename VectorType>
1165inline void
1166MGSmootherPrecondition<MatrixType, PreconditionerType, VectorType>::smooth(
1167 const unsigned int level,
1168 VectorType &u,
1169 const VectorType &rhs) const
1170{
1171 unsigned int maxlevel = matrices.max_level();
1172 unsigned int steps2 = this->steps;
1173
1174 if (this->variable)
1175 steps2 *= (1 << (maxlevel - level));
1176
1177 typename VectorMemory<VectorType>::Pointer r(this->vector_memory);
1178 typename VectorMemory<VectorType>::Pointer d(this->vector_memory);
1179
1180 r->reinit(u, true);
1181 d->reinit(u, true);
1182
1183 bool T = this->transpose;
1184 if (this->symmetric && (steps2 % 2 == 0))
1185 T = false;
1186 if (this->debug > 0)
1187 deallog << 'S' << level << ' ';
1188
1189 for (unsigned int i = 0; i < steps2; ++i)
1190 {
1191 if (T)
1192 {
1193 if (this->debug > 0)
1194 deallog << 'T';
1195 matrices[level].Tvmult(*r, u);
1196 r->sadd(-1., 1., rhs);
1197 if (this->debug > 2)
1198 deallog << ' ' << r->l2_norm() << ' ';
1199 smoothers[level].Tvmult(*d, *r);
1200 if (this->debug > 1)
1201 deallog << ' ' << d->l2_norm() << ' ';
1202 }
1203 else
1204 {
1205 if (this->debug > 0)
1206 deallog << 'N';
1207 matrices[level].vmult(*r, u);
1208 r->sadd(-1., rhs);
1209 if (this->debug > 2)
1210 deallog << ' ' << r->l2_norm() << ' ';
1211 smoothers[level].vmult(*d, *r);
1212 if (this->debug > 1)
1213 deallog << ' ' << d->l2_norm() << ' ';
1214 }
1215 u += *d;
1216 if (this->symmetric)
1217 T = !T;
1218 }
1219 if (this->debug > 0)
1220 deallog << std::endl;
1221}
1222
1223
1224
1225template <typename MatrixType, typename PreconditionerType, typename VectorType>
1226inline void
1227MGSmootherPrecondition<MatrixType, PreconditionerType, VectorType>::apply(
1228 const unsigned int level,
1229 VectorType &u,
1230 const VectorType &rhs) const
1231{
1232 unsigned int maxlevel = matrices.max_level();
1233 unsigned int steps2 = this->steps;
1234
1235 if (this->variable)
1236 steps2 *= (1 << (maxlevel - level));
1237
1238 bool T = this->transpose;
1239 if (this->symmetric && (steps2 % 2 == 0))
1240 T = false;
1241 if (this->debug > 0)
1242 deallog << 'S' << level << ' ';
1243
1244 // first step where we overwrite the result
1245 if (this->debug > 2)
1246 deallog << ' ' << rhs.l2_norm() << ' ';
1247 if (this->debug > 0)
1248 deallog << (T ? 'T' : 'N');
1249 if (T)
1250 smoothers[level].Tvmult(u, rhs);
1251 else
1252 smoothers[level].vmult(u, rhs);
1253 if (this->debug > 1)
1254 deallog << ' ' << u.l2_norm() << ' ';
1255 if (this->symmetric)
1256 T = !T;
1257
1258 typename VectorMemory<VectorType>::Pointer r(this->vector_memory);
1259 typename VectorMemory<VectorType>::Pointer d(this->vector_memory);
1260
1261 if (steps2 > 1)
1262 {
1263 r->reinit(u, true);
1264 d->reinit(u, true);
1265 }
1266
1267 for (unsigned int i = 1; i < steps2; ++i)
1268 {
1269 if (T)
1270 {
1271 if (this->debug > 0)
1272 deallog << 'T';
1273 matrices[level].Tvmult(*r, u);
1274 r->sadd(-1., 1., rhs);
1275 if (this->debug > 2)
1276 deallog << ' ' << r->l2_norm() << ' ';
1277 smoothers[level].Tvmult(*d, *r);
1278 if (this->debug > 1)
1279 deallog << ' ' << d->l2_norm() << ' ';
1280 }
1281 else
1282 {
1283 if (this->debug > 0)
1284 deallog << 'N';
1285 matrices[level].vmult(*r, u);
1286 r->sadd(-1., rhs);
1287 if (this->debug > 2)
1288 deallog << ' ' << r->l2_norm() << ' ';
1289 smoothers[level].vmult(*d, *r);
1290 if (this->debug > 1)
1291 deallog << ' ' << d->l2_norm() << ' ';
1292 }
1293 u += *d;
1294 if (this->symmetric)
1295 T = !T;
1296 }
1297 if (this->debug > 0)
1298 deallog << std::endl;
1299}
1300
1301
1302
1303template <typename MatrixType, typename PreconditionerType, typename VectorType>
1304inline std::size_t
1305MGSmootherPrecondition<MatrixType, PreconditionerType, VectorType>::
1306 memory_consumption() const
1307{
1308 return sizeof(*this) + matrices.memory_consumption() +
1309 smoothers.memory_consumption() +
1310 this->vector_memory.memory_consumption();
1311}
1312
1313
1314#endif // DOXYGEN
1315
1316DEAL_II_NAMESPACE_CLOSE
1317
1318#endif
LogStream::memory_consumption
std::size_t memory_consumption() const
MGLevelObject::memory_consumption
std::size_t memory_consumption() const
Definition mg_level_object.h:340
MGLevelObject::max_level
unsigned int max_level() const
Definition mg_level_object.h:313
MGLevelObject::min_level
unsigned int min_level() const
Definition mg_level_object.h:305
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:567
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:580
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:411
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:398
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:99
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:128
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:105
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:498
DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:499
transpose
DerivativeForm< 1, spacedim, dim, Number > transpose(const DerivativeForm< 1, dim, spacedim, Number > &DF)
Definition derivative_form.h:610
level
unsigned int level
Definition grid_out.cc:4632
Assert
#define Assert(cond, exc)
Definition exceptions.h:1659
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:735
LAPACKSupport::symmetric
@ symmetric
Matrix is symmetric.
Definition lapack_support.h:112
LinearAlgebra::TpetraWrappers::TpetraTypes::VectorType
Tpetra::Vector< Number, LO, GO, NodeType< MemorySpace > > VectorType
Definition trilinos_tpetra_types.h:112
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:1639
mg
Definition mg.h:81
std::min
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
Definition vectorization.h:6960
std::max
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
Definition vectorization.h:6943
observer_pointer.h
Generated by doxygen 1.9.8