Line data Source code
1 : /*
2 : * Copyright (c) 2013-2015: G-CSC, Goethe University Frankfurt
3 : * Author: Raphael Prohl
4 : *
5 : * This file is part of UG4.
6 : *
7 : * UG4 is free software: you can redistribute it and/or modify it under the
8 : * terms of the GNU Lesser General Public License version 3 (as published by the
9 : * Free Software Foundation) with the following additional attribution
10 : * requirements (according to LGPL/GPL v3 §7):
11 : *
12 : * (1) The following notice must be displayed in the Appropriate Legal Notices
13 : * of covered and combined works: "Based on UG4 (www.ug4.org/license)".
14 : *
15 : * (2) The following notice must be displayed at a prominent place in the
16 : * terminal output of covered works: "Based on UG4 (www.ug4.org/license)".
17 : *
18 : * (3) The following bibliography is recommended for citation and must be
19 : * preserved in all covered files:
20 : * "Reiter, S., Vogel, A., Heppner, I., Rupp, M., and Wittum, G. A massively
21 : * parallel geometric multigrid solver on hierarchically distributed grids.
22 : * Computing and visualization in science 16, 4 (2013), 151-164"
23 : * "Vogel, A., Reiter, S., Rupp, M., Nägel, A., and Wittum, G. UG4 -- a novel
24 : * flexible software system for simulating pde based models on high performance
25 : * computers. Computing and visualization in science 16, 4 (2013), 165-179"
26 : *
27 : * This program is distributed in the hope that it will be useful,
28 : * but WITHOUT ANY WARRANTY; without even the implied warranty of
29 : * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
30 : * GNU Lesser General Public License for more details.
31 : */
32 :
33 : #ifndef __H__UG__LIB_ALGEBRA__OPERATOR__PRECONDITIONER__PROJECTED_GAUSS_SEIDEL__PROJ_GAUSS_SEIDEL__
34 : #define __H__UG__LIB_ALGEBRA__OPERATOR__PRECONDITIONER__PROJECTED_GAUSS_SEIDEL__PROJ_GAUSS_SEIDEL__
35 :
36 : #include "proj_gauss_seidel_interface.h"
37 :
38 : namespace ug{
39 :
40 : /// Projected GaussSeidel (SOR) -method
41 : /**
42 : * The projected GaussSeidel method can be applied to problems of the form
43 : *
44 : * A * u >= b (I)
45 : * c(u) >= 0 (II)
46 : * c(u)^T * [A*u - b] = 0, (III)
47 : *
48 : * where u, b are vectors and A is a matrix. '*' denotes componentwise multiplication.
49 : * c(u) denotes an obstacle-function, which depends on the solution vector u. One possible
50 : * example for such an obstacle-function could be the scalar obstacle function
51 : *
52 : * u >= 0.
53 : *
54 : * The obstacle function c(u) is defined by creating an instance of IObstacleConstraint, which is
55 : * passed to the projected preconditioner by the method 'IProjPreconditioner::set_obstacle_constraint'.
56 : *
57 : * Similar problems, which e.g. only differ in the sign in (I) and/or (II) can be
58 : * equivalently treated by the method.
59 : *
60 : * Note: Due to (II) the old solution needs to be stored within this method.
61 : * This is a difference to the classical smoothers/preconditioners, which usually work
62 : * on the correction and the defect only.
63 : *
64 : * By calling 'set_sor_damp(number)' one gets the successive overrelaxation-version of the
65 : * projected preconditioners of GaussSeidel type.
66 : *
67 : * References:
68 : * <ul>
69 : * <li> A. Brandt and C. W. Cryer. Multigrid algorithms for the solution of linear
70 : * <li> complementarity problems arising from free boundary problems. SIAM J. SCI. STAT. COMPUT. Vol. 4, No. 4 (1983)
71 : * </ul>
72 : *
73 : * \tparam TAlgebra Algebra type
74 : */
75 :
76 : template <typename TDomain, typename TAlgebra>
77 : class ProjGaussSeidel:
78 : public IProjGaussSeidel<TDomain,TAlgebra>
79 : {
80 : public:
81 : /// Base class type
82 : typedef IProjGaussSeidel<TDomain,TAlgebra> base_type;
83 :
84 : /// Algebra type
85 : typedef TAlgebra algebra_type;
86 :
87 : /// Matrix type
88 : typedef typename algebra_type::matrix_type matrix_type;
89 :
90 : /// Vector type
91 : typedef typename algebra_type::vector_type vector_type;
92 :
93 : protected:
94 : /// name
95 0 : virtual const char* name() const {return "Projected GaussSeidel";}
96 :
97 : public:
98 : ProjGaussSeidel() : base_type() {}
99 : /// copy constructor
100 0 : ProjGaussSeidel( const ProjGaussSeidel<TDomain, TAlgebra> &parent )
101 0 : : base_type(parent)
102 : { }
103 :
104 : /// Clone
105 0 : virtual SmartPtr<ILinearIterator<vector_type> > clone()
106 : {
107 0 : return make_sp(new ProjGaussSeidel<TDomain, algebra_type>(*this));
108 : }
109 :
110 : /// computes a new correction c = B*d and projects on the underlying constraint
111 : virtual void step(const matrix_type& mat, vector_type& c, const vector_type& d, const number relax);
112 : };
113 :
114 : template <typename TDomain, typename TAlgebra>
115 : class ProjBackwardGaussSeidel:
116 : public IProjGaussSeidel<TDomain,TAlgebra>
117 : {
118 : public:
119 : /// Base class type
120 : typedef IProjGaussSeidel<TDomain,TAlgebra> base_type;
121 :
122 : /// Algebra type
123 : typedef TAlgebra algebra_type;
124 :
125 : /// Matrix type
126 : typedef typename algebra_type::matrix_type matrix_type;
127 :
128 : /// Vector type
129 : typedef typename algebra_type::vector_type vector_type;
130 :
131 : protected:
132 : /// name
133 0 : virtual const char* name() const {return "Projected Backward GaussSeidel";}
134 :
135 : public:
136 0 : ProjBackwardGaussSeidel() : base_type() {}
137 : /// clone constructor
138 0 : ProjBackwardGaussSeidel( const ProjBackwardGaussSeidel<TDomain, TAlgebra> &parent )
139 0 : : base_type(parent)
140 : { }
141 :
142 : /// Clone
143 0 : virtual SmartPtr<ILinearIterator<vector_type> > clone()
144 : {
145 0 : return make_sp(new ProjBackwardGaussSeidel<TDomain, algebra_type>(*this));
146 : }
147 :
148 :
149 : /// computes a new correction c = B*d and projects on the underlying constraint
150 : virtual void step(const matrix_type& mat, vector_type& c, const vector_type& d, const number relax);
151 : };
152 :
153 :
154 : template <typename TDomain, typename TAlgebra>
155 : class ProjSymmetricGaussSeidel:
156 : public IProjGaussSeidel<TDomain,TAlgebra>
157 : {
158 : public:
159 : /// Base class type
160 : typedef IProjGaussSeidel<TDomain,TAlgebra> base_type;
161 :
162 : /// Algebra type
163 : typedef TAlgebra algebra_type;
164 :
165 : /// Matrix type
166 : typedef typename algebra_type::matrix_type matrix_type;
167 :
168 : /// Vector type
169 : typedef typename algebra_type::vector_type vector_type;
170 :
171 : protected:
172 : /// name
173 0 : virtual const char* name() const {return "Projected Symmetric GaussSeidel";}
174 :
175 : public:
176 0 : ProjSymmetricGaussSeidel() : base_type() {}
177 :
178 : /// copy constructor
179 0 : ProjSymmetricGaussSeidel( const ProjSymmetricGaussSeidel<TDomain, TAlgebra> &parent )
180 0 : : base_type(parent)
181 : { }
182 :
183 : /// Clone
184 0 : virtual SmartPtr<ILinearIterator<vector_type> > clone()
185 : {
186 0 : return make_sp(new ProjSymmetricGaussSeidel<TDomain, algebra_type>(*this));
187 : }
188 :
189 : /// computes a new correction c = B*d and projects on the underlying constraint
190 : virtual void step(const matrix_type& mat, vector_type& c, const vector_type& d, const number relax);
191 : };
192 :
193 : } // end namespace ug
194 :
195 : // include implementation
196 : #include "proj_gauss_seidel_impl.h"
197 :
198 : #endif /* __H__UG__LIB_ALGEBRA__OPERATOR__PRECONDITIONER__PROJECTED_GAUSS_SEIDEL__PROJ_GAUSS_SEIDEL__ */
|
