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 ACTIVE_SET_H_
34 : #define ACTIVE_SET_H_
35 :
36 : #include "lib_algebra/active_set/lagrange_multiplier_disc_interface.h"
37 : #include "lib_disc/function_spaces/grid_function.h"
38 :
39 : using namespace std;
40 :
41 : namespace ug {
42 :
43 : /// Active Set method
44 : /**
45 : * The active Set method is a well-known method in constrained optimization theory.
46 : * A general formulation for these problems reads
47 : * \f{eqnarray*}{
48 : * min_{x \in \mathbb{R}^n} f(x)
49 : * \f}
50 : * s.t.
51 : * \f{eqnarray*}{
52 : * c_i(x) = 0, \qquad i \in E \\
53 : * c_i(x) \ge 0 \qquad i \in I,
54 : * \f}
55 : *
56 : * where \f$ f \f$, \f$ c_i \f$ are smooth, real-valued functions on a subset of \f$ \mathbb{R}^n \f$.
57 : * \f$ I \f$ (set of inequality constraints) and \f$ E \f$ (set of equality constraints)
58 : * are two finite sets of indices. \f$ f \f$ is called objective function.
59 : *
60 : * The active Set \f$ A \f$ is defined as:
61 : * \f{eqnarray*}{
62 : * A(x) := E \cup \{ i \in I | c_i(x) = 0 \},
63 : * \f}
64 : * i.e. for \f$ i \in I \f$ the inequality constraint is said to be active, if \f$ c_i(x) = 0 \f$.
65 : * Otherwise (\f$ c_i(x) > 0 \f$) it is called inactive.
66 : *
67 : * A common approach to treat the inequality constraints is its reformulation as
68 : * equations by using so called complementarity functions, see e.g. C.Hager und B. I. Wohlmuth:
69 : * "Hindernis- und Kontaktprobleme" for a simple introduction into this topic.
70 : * By means of complementarity functions, constraints of the form
71 : * \f{eqnarray*}{
72 : * a \ge 0, \, b \ge 0, \, a b = 0
73 : * \f}
74 : * with \f$ a, b \in \mathbb{R}^n \f$ can be reformulated as
75 : * \f{eqnarray*}{
76 : * C(a,b) = 0,
77 : * \f}
78 : * with \f$ C: \mathbb{R}^n x \mathbb{R}^n \to \mathbb{R}^n \f$ being an appropriate complementarity function.
79 : * The value of \f$ C \f$ indicates, whether the index is active or inactive (see method 'active_index').
80 : * Thus, it determines in every step the set of active indices, for which the original system
81 : * of equations needs to be adapted. The influence of these active indices on the original system
82 : * of equations ( \f$ K u = f \f$, with \f$ K \f$: system-matrix; \f$ u, f \f$ vectors) can be modelled by means of a
83 : * lagrange multiplier '\f$ \lambda \f$'. \f$ \lambda \f$ can either be defined by the residual ( \f$ \lambda := f - K u \f$,
84 : * see method 'residual_lagrange_mult') or by a problem-dependent computation (see method 'lagrange_multiplier').
85 : *
86 : * In every Active Set step a linear or linearized system is solved. The algorithm stops when the active
87 : * and inactive Set remains unchanged (see method 'check_conv').
88 : *
89 : *
90 : * References:
91 : * <ul>
92 : * <li> J. Nocedal and S. J. Wright. Numerical optimization.(2000)
93 : * </ul>
94 : *
95 : * \tparam TDomain Domain type
96 : * \tparam TAlgebra Algebra type
97 : */
98 : template <typename TDomain, typename TAlgebra>
99 : class ActiveSet
100 : {
101 : public:
102 : /// Type of algebra
103 : typedef TAlgebra algebra_type;
104 :
105 : /// Type of algebra matrix
106 : typedef typename algebra_type::matrix_type matrix_type;
107 :
108 : /// Type of algebra vector
109 : typedef typename algebra_type::vector_type vector_type;
110 :
111 : /// Type of algebra value
112 : typedef typename vector_type::value_type value_type;
113 :
114 : /// Type of grid function
115 : typedef GridFunction<TDomain, TAlgebra> function_type;
116 :
117 : /// base element type of associated domain
118 : typedef typename domain_traits<TDomain::dim>::grid_base_object TBaseElem;
119 :
120 : /// domain dimension
121 : static const int dim = TDomain::dim;
122 :
123 : public:
124 : /// constructor
125 0 : ActiveSet() : m_bObs(false), m_spLagMultDisc(NULL) {
126 : // specifies the number of fcts
127 : //value_type u_val;
128 : //m_nrFcts = GetSize(u_val); //ToDo: This field is only used in check_dist_to_obs which is commented out. Remove it completely?
129 : };
130 :
131 : /// sets obstacle gridfunction, which limits the solution u
132 0 : void set_obstacle(ConstSmartPtr<function_type> obs) {
133 0 : m_spObs = obs; m_bObs = true;
134 : // if 'obs'-gridfunction is defined on a subset,
135 : // which is not a boundary-subset -> UG_LOG
136 0 : }
137 :
138 : /// sets a discretization in order to compute the lagrange multiplier
139 0 : void set_lagrange_multiplier_disc(
140 : SmartPtr<ILagrangeMultiplierDisc<TDomain, function_type> > lagMultDisc)
141 0 : {m_spLagMultDisc = lagMultDisc;};
142 :
143 : void prepare(function_type& u);
144 :
145 : /// checks the distance to the prescribed obstacle/constraint
146 : //bool check_dist_to_obs(vector_type& u);
147 :
148 : template <typename TElem, typename TIterator>
149 : void active_index_elem(TIterator iterBegin,
150 : TIterator iterEnd, function_type& u,
151 : function_type& rhs, function_type& lagrangeMult);
152 :
153 : /// determines the active indices, stores them in a vector and sets dirichlet values in rhs for active indices
154 : bool active_index(function_type& u, function_type& rhs, function_type& lagrangeMult,
155 : function_type& gap);
156 :
157 : void set_dirichlet_rows(matrix_type& mat);
158 :
159 : /// computes the lagrange multiplier for a given disc
160 : void lagrange_multiplier(function_type& lagrangeMult, const function_type& u);
161 :
162 : /// computes the lagrange multiplier by means of
163 : /// the residuum (lagMult = rhs - mat * u)
164 : void residual_lagrange_mult(vector_type& lagMult, const matrix_type& mat,
165 : const vector_type& u, vector_type& rhs);
166 :
167 : template <typename TElem, typename TIterator>
168 : bool check_conv_elem(TIterator iterBegin,
169 : TIterator iterEnd, function_type& u, const function_type& lambda);
170 :
171 : /// checks if all constraints are fulfilled & the activeSet remained unchanged
172 : bool check_conv(function_type& u, const function_type& lambda, const size_t step);
173 :
174 : /// checks if all inequalities are fulfilled
175 : bool check_inequ(const matrix_type& mat, const vector_type& u,
176 : const vector_type& lagrangeMult, const vector_type& rhs);
177 :
178 : template <typename TElem, typename TIterator>
179 : void lagrange_mat_inv_elem(TIterator iterBegin,
180 : TIterator iterEnd, matrix_type& lagrangeMatInv);
181 :
182 : void lagrange_mat_inv(matrix_type& lagrangeMatInv);
183 :
184 : private:
185 : /// pointer to the DofDistribution on the whole domain
186 : SmartPtr<DoFDistribution> m_spDD;
187 : SmartPtr<TDomain> m_spDom;
188 :
189 : /// number of functions
190 : //size_t m_nrFcts; //ToDo: This field is only used in check_dist_to_obs which is commented out. Remove it completely?
191 :
192 : /// pointer to a gridfunction describing an obstacle-constraint
193 : ConstSmartPtr<function_type> m_spObs;
194 : bool m_bObs;
195 :
196 : /// pointer to a lagrangeMultiplier-Disc
197 : SmartPtr<ILagrangeMultiplierDisc<TDomain, function_type> > m_spLagMultDisc;
198 :
199 : /// vector of possible active subsets
200 : vector<int> m_vActiveSubsets;
201 :
202 : /*template <typename TElem>
203 : struct activeElemAndLocInd{
204 : TElem* pElem; // pointer to active elem
205 : vector<vector<size_t> > vlocInd; // vector of local active indices
206 : };*/
207 :
208 : /// vector of the current active set of global DoFIndices
209 : vector<DoFIndex> m_vActiveSetGlob;
210 : /// vector remembering the active set of global DoFIndices
211 : vector<DoFIndex> m_vActiveSetGlobOld;
212 : };
213 :
214 : } // namespace ug
215 :
216 : #include "active_set_impl.h"
217 :
218 : #endif /* ACTIVE_SET_H_ */
|
