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_INTERFACE__
34 : #define __H__UG__LIB_ALGEBRA__OPERATOR__PRECONDITIONER__PROJECTED_GAUSS_SEIDEL__PROJ_GAUSS_SEIDEL_INTERFACE__
35 :
36 : #include "obstacles/obstacle_constraint_interface.h"
37 : #include "lib_algebra/operator/preconditioner/gauss_seidel.h"
38 :
39 : namespace ug{
40 :
41 : /// Interface for Projected GaussSeidel Preconditioner
42 : /**
43 : * This class provides an interface to define a preconditioner which can be applied to solve
44 : * problems of the form
45 : *
46 : * A * u >= b (I)
47 : * c(u) >= 0 (II)
48 : * c(u)^T * [A*u - b] = 0, (III)
49 : *
50 : * where u, b are vectors and A is a matrix. '*' denotes componentwise multiplication.
51 : * c(u) denotes an obstacle-function, which depends on the solution vector u. One possible
52 : * example for such an obstacle-function could be the scalar obstacle function
53 : *
54 : * u >= 0.
55 : *
56 : * The obstacle function c(u) is defined by creating an instance of IObstacleConstraint, which is
57 : * passed to the projected preconditioner by the method 'set_obstacle_constraint'.
58 : *
59 : * Similar problems, which e.g. only differ in the sign in (I) and/or (II) can be
60 : * equivalently treated by these preconditioners.
61 : *
62 : * Note: Due to (II) the old solution needs to be stored within this method.
63 : * This is a difference to the classical smoothers/preconditioners, which usually work
64 : * on the correction and the defect only.
65 : *
66 : * Since the problem formulation (I)-(III) consists of inequalities, the projected preconditioner
67 : * performs a projection on a constraint c(u) in every preconditioner-step.
68 : *
69 : * \tparam TAlgebra Algebra type
70 : */
71 : template <typename TDomain, typename TAlgebra>
72 : class IProjGaussSeidel:
73 : public GaussSeidelBase<TAlgebra>
74 : {
75 : public:
76 : /// Base class type
77 : typedef GaussSeidelBase<TAlgebra> base_type;
78 :
79 : /// Algebra type
80 : typedef TAlgebra algebra_type;
81 :
82 : /// Matrix type
83 : typedef typename algebra_type::matrix_type matrix_type;
84 :
85 : /// Vector type
86 : typedef typename algebra_type::vector_type vector_type;
87 :
88 : /// Value type
89 : typedef typename vector_type::value_type value_type;
90 :
91 : /// Grid Function type
92 : typedef GridFunction<TDomain, TAlgebra> GF;
93 :
94 : public:
95 : /// constructor
96 0 : IProjGaussSeidel(): GaussSeidelBase<TAlgebra>(){
97 : m_spvObsConstraint.clear();
98 0 : m_bObsCons = false;
99 : };
100 :
101 :
102 : /// clone constructor
103 0 : IProjGaussSeidel( const IProjGaussSeidel<TDomain, TAlgebra> &parent )
104 0 : : base_type(parent)
105 : {
106 0 : m_spvObsConstraint = parent.m_spvObsConstraint;
107 0 : m_bObsCons = parent.m_bObsCons;
108 0 : }
109 :
110 : /// adds the obstacle constraint function c(u)
111 0 : void add_obstacle_constraint(SmartPtr<IObstacleConstraint<TDomain,TAlgebra> > spObsCons)
112 : {
113 0 : m_spvObsConstraint.push_back(spObsCons);
114 0 : m_bObsCons = true;
115 :
116 : // inits the obstacle constraint
117 0 : spObsCons->init();
118 0 : }
119 :
120 : /// Destructor
121 0 : ~IProjGaussSeidel(){};
122 :
123 : /// name
124 : virtual const char* name() const = 0;
125 :
126 : /// Prepare for Operator J(u) and linearization point u (current solution)
127 : virtual bool init(SmartPtr<ILinearOperator<vector_type> > J, const vector_type& u);
128 :
129 : /// computes a new correction c = B*d and projects on the underlying constraint
130 : /**
131 : * This method computes a new correction c = B*d. B is here the underlying matrix operator.
132 : *
133 : * \param[in] mat underlying matrix (i.e. A in A*u = b)
134 : * \param[out] c correction
135 : * \param[in] d defect
136 : * \param[in] relax relaxation parameter
137 : */
138 : virtual void step(const matrix_type& mat, vector_type& c, const vector_type& d, const number relax) = 0;
139 :
140 : /// projects the correction on the underlying constraints set by the obstacleConstraints
141 : void project_correction(value_type& c_i, const size_t i);
142 :
143 : /// Compute new correction c = B*d
144 : virtual bool apply(vector_type& c, const vector_type& d);
145 :
146 : /// Compute new correction c = B*d and return new defect d := d - A*c
147 : virtual bool apply_update_defect(vector_type& c, vector_type& d);
148 :
149 : private:
150 : /// for all indices stored in vInd:
151 : /// the entry of vec is set to zero
152 : void truncateVec(vector_type& vec, vector<DoFIndex>& vInd);
153 :
154 : /// for all indices stored in vInd:
155 : /// all rows and columns of mat are set to zero
156 : void truncateMat(matrix_type& mat, vector<DoFIndex>& vInd);
157 :
158 : protected:
159 : /// obstacle constraint
160 : vector<SmartPtr<IObstacleConstraint<TDomain,TAlgebra> > > m_spvObsConstraint;
161 :
162 : private:
163 : /// pointer to solution
164 : SmartPtr<vector_type> m_spSol;
165 :
166 : /// flag indicating if obstacle constraint has been set
167 : bool m_bObsCons;
168 :
169 : /// init flag indicating if init has been called
170 : bool m_bInit;
171 :
172 : };
173 :
174 :
175 : } // end namespace ug
176 :
177 : // include implementation
178 : #include "proj_gauss_seidel_interface_impl.h"
179 :
180 : #endif /* __H__UG__LIB_ALGEBRA__OPERATOR__PRECONDITIONER__PROJECTED_GAUSS_SEIDEL__PROJ_GAUSS_SEIDEL_INTERFACE__ */
|
