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__OBSTACLE_CONSTRAINT_INTERFACE__
34 : #define __H__UG__LIB_ALGEBRA__OPERATOR__PRECONDITIONER__PROJECTED_GAUSS_SEIDEL__OBSTACLE_CONSTRAINT_INTERFACE__
35 :
36 : #include "lib_disc/common/multi_index.h"
37 : #include "lib_disc/function_spaces/grid_function.h"
38 : #include "lib_disc/spatial_disc/user_data/const_user_data.h"
39 :
40 : #include "lib_disc/spatial_disc/constraints/constraint_interface.h"
41 :
42 : using namespace std;
43 :
44 : namespace ug{
45 :
46 : /// Interface for Obstacle Constraints
47 : /**
48 : * The Interface for Obstacle Constraints provides the framework to define obstacle constraints
49 : * of the form
50 : *
51 : * c(u) <= upObs (cf. 'set_upper_obstacle')
52 : *
53 : * and
54 : *
55 : * c(u) >= lowObs (cf. 'set_lower_obstacle')
56 : *
57 : * where u is the solution vector. Here, 'upObs' and 'lowObs' are user-defined functions,
58 : * which need to be of the same size as the function of unknowns u. The obstacle function
59 : * c is defined by creating a derived class of this interface and by specializing the way
60 : * the correction should be computed (cf. correction_for_lower_obs, correction_for_upper_obs,
61 : * correction_for_lower_and_upper_obs).
62 : *
63 : * One simple example is the ScalarObstacle-class. Such an obstacle functions can be used in
64 : * combination with projected preconditioners. They should be passed to the preconditioner
65 : * by 'IProjPreconditioner::set_obstacle_constraint'.
66 : */
67 : template <typename TDomain, typename TAlgebra>
68 : class IObstacleConstraint:
69 : public IDomainConstraint<TDomain, TAlgebra>
70 : // TODO: think about, restructuring the IConstraint-Interface in order to distinguish between
71 : // a constraint having an effect on the assembling process and a constraint, which is considered in the
72 : // solver (e.g. by means of 'adjust_restriction') only!
73 : {
74 : public:
75 : /// Base Type
76 : typedef IDomainConstraint<TDomain, TAlgebra> base_type;
77 :
78 : /// Algebra type
79 : typedef TAlgebra algebra_type;
80 :
81 : /// Matrix type
82 : typedef typename algebra_type::matrix_type matrix_type;
83 :
84 : /// Vector type
85 : typedef typename algebra_type::vector_type vector_type;
86 :
87 : /// Value type
88 : typedef typename vector_type::value_type value_type;
89 :
90 : /// Type of domain
91 : typedef TDomain domain_type;
92 :
93 : /// world Dimension
94 : static const int dim = domain_type::dim;
95 :
96 : /// Type of position coordinates (e.g. position_type)
97 : typedef typename domain_type::position_type position_type;
98 :
99 : public:
100 : /// constructor for an obstacle defined on some subset(s)
101 0 : IObstacleConstraint(const GridFunction<TDomain, TAlgebra>& u)
102 0 : {
103 0 : clear();
104 :
105 0 : m_spDD = u.dof_distribution();
106 0 : m_spDomain = u.domain();
107 0 : };
108 :
109 : /// constructor
110 0 : IObstacleConstraint(){
111 : //clear();
112 0 : UG_THROW("In 'IObstacleConstraint()': A constructor with a GridFunction as parameter"
113 : "is needed here!");
114 0 : };
115 :
116 : /// adds a lua callback (cond and non-cond)
117 : #ifdef UG_FOR_LUA
118 : void add(const char* name, const char* function);
119 : void add(const char* name, const char* function, const char* subsets);
120 : #endif
121 :
122 : /// adds a conditional user-defined value as dirichlet condition for a function on subsets and on whole domain
123 : void add(SmartPtr<UserData<number, dim, bool> > func, const char* function);
124 : void add(SmartPtr<UserData<number, dim, bool> > func, const char* function, const char* subsets);
125 :
126 : /// adds a user-defined value as dirichlet condition for a function on subsets and on whole domain
127 : void add(SmartPtr<UserData<number, dim> > func, const char* function);
128 : void add(SmartPtr<UserData<number, dim> > func, const char* function, const char* subsets);
129 :
130 : /// adds a constant value as dirichlet condition for a function on subsets and on whole domain
131 : void add(number value, const char* function);
132 : void add(number value, const char* function, const char* subsets);
133 :
134 : /// adds a user-defined vector as dirichlet condition for a vector-function on subsets and on whole domain
135 : void add(SmartPtr<UserData<MathVector<dim>, dim> > func, const char* functions);
136 : void add(SmartPtr<UserData<MathVector<dim>, dim> > func, const char* functions, const char* subsets);
137 :
138 :
139 : /// init function calls 'extract_data' and 'obstacle_value'-methods in order to
140 : /// store the obstacle values set by UserDatas
141 : void init();
142 :
143 : /// checks if a given dof is in an obstacle subset
144 : bool is_obs_dof(const DoFIndex& dof);
145 :
146 : /// resets the vector storing the active dofs
147 0 : void reset_active_dofs(){m_vActiveDofs.resize(0);}
148 :
149 : /// returns the vector storing the active dofs
150 : void active_dofs(vector<DoFIndex>& vActiveDoFs)
151 0 : {vActiveDoFs = m_vActiveDofs;}
152 :
153 :
154 : /// this preprocess-method is called in every init-call of the underlying proj. preconditioner
155 : /// it is useful to attach e.g. additional data to the obstacle DoFs, like the e.g. the normal vector at
156 : /// this DoF
157 0 : virtual void preprocess(){};
158 :
159 : /// projects the i-th index of the solution onto the admissible set and adjusts the correction
160 : virtual void adjust_sol_and_cor(value_type& sol_i, value_type& c_i, bool& dofIsAdmissible,
161 : const DoFIndex& dof) = 0;
162 :
163 : /// the defect needs to be adjusted for the active indices (those indices, which are in contact)
164 : virtual void adjust_defect_to_constraint(vector_type& d) = 0;
165 :
166 :
167 : /// restricts the obstacle values to a coarser grid in a multigrid hierarchy
168 : virtual void restrict_obs_values() = 0;
169 :
170 :
171 : /// Destructor
172 0 : virtual ~IObstacleConstraint(){};
173 :
174 : public:
175 : ///////////////////////////////
176 : // Implement Interface
177 : ///////////////////////////////
178 :
179 : /// sets a unity row for all dirichlet indices
180 0 : void adjust_jacobian(matrix_type& J, const vector_type& u,
181 : ConstSmartPtr<DoFDistribution> dd, int type, number time = 0.0,
182 : ConstSmartPtr<VectorTimeSeries<vector_type> > vSol = NULL,
183 : const number s_a0 = 1.0){
184 : UG_LOG("IObstacleConstraint::adjust_jacobian() \n");
185 0 : };
186 :
187 : /// sets a zero value in the defect for all dirichlet indices
188 0 : void adjust_defect(vector_type& d, const vector_type& u,
189 : ConstSmartPtr<DoFDistribution> dd, int type, number time = 0.0,
190 : ConstSmartPtr<VectorTimeSeries<vector_type> > vSol = NULL,
191 : const std::vector<number>* vScaleMass = NULL,
192 : const std::vector<number>* vScaleStiff = NULL){
193 : UG_LOG("IObstacleConstraint::adjust_defect() \n");
194 0 : };
195 :
196 : /// sets the dirichlet value in the solution for all dirichlet indices
197 0 : void adjust_solution(vector_type& u,
198 : ConstSmartPtr<DoFDistribution> dd, int type, number time = 0.0){
199 : UG_LOG("IObstacleConstraint::adjust_solution() \n");
200 0 : };
201 :
202 : /// sets unity rows in A and dirichlet values in right-hand side b
203 0 : void adjust_linear(matrix_type& A, vector_type& b,
204 : ConstSmartPtr<DoFDistribution> dd, int type, number time = 0.0){
205 : UG_LOG("IObstacleConstraint::adjust_linear() \n");
206 0 : };
207 :
208 : /// sets the dirichlet value in the right-hand side
209 0 : void adjust_rhs(vector_type& b, const vector_type& u,
210 : ConstSmartPtr<DoFDistribution> dd, int type, number time = 0.0){
211 : UG_LOG("IObstacleConstraint::adjust_rhs() \n");
212 0 : };
213 :
214 : /// sets constraints in prolongation
215 0 : virtual void adjust_prolongation(matrix_type& P,
216 : ConstSmartPtr<DoFDistribution> ddFine,
217 : ConstSmartPtr<DoFDistribution> ddCoarse,
218 : int type,
219 : number time = 0.0){
220 : UG_LOG("IObstacleConstraint::adjust_prolongationP() \n");
221 0 : };
222 :
223 : /// sets constraints in restriction
224 : virtual void adjust_restriction(matrix_type& R,
225 : ConstSmartPtr<DoFDistribution> ddCoarse,
226 : ConstSmartPtr<DoFDistribution> ddFine,
227 : int type,
228 : number time = 0.0);
229 :
230 : /// returns the type of the constraints
231 : //TODO: does another type make more sense?
232 0 : virtual int type() const {return CT_CONSTRAINTS;}
233 :
234 : private:
235 : /// extract the UserDatas
236 : void extract_data();
237 :
238 : template <typename TUserData, typename TScheduledUserData>
239 : void extract_data(std::map<int, std::vector<TUserData*> >& mvUserDataObsSegment,
240 : std::vector<TScheduledUserData>& vUserData);
241 :
242 : /// clear all UserData-member variables
243 : void clear();
244 :
245 : void check_functions_and_subsets(FunctionGroup& functionGroup, SubsetGroup& subsetGroup,
246 : size_t numFct) const;
247 :
248 : /// store the obstacle value set by means of UserDatas
249 : void init_obstacle_dofs_with_values(number time);
250 :
251 : template <typename TUserData>
252 : void init_obstacle_dofs_with_values(const std::map<int, std::vector<TUserData*> >& mvUserData, number time);
253 :
254 : template <typename TBaseElem, typename TUserData>
255 : void init_obstacle_dofs_with_values(const std::vector<TUserData*>& vUserData, int si, number time);
256 :
257 : protected:
258 : /// map to store obstacle values with its corresponding DoFs
259 : map<DoFIndex, number> m_mObstacleValues;
260 :
261 : /// stores the subset-indices of the obstacle subsets
262 : vector<int> m_vObsSubsets;
263 :
264 : /// stores the dofs, which satisfy the constraints with equality
265 : vector<DoFIndex> m_vActiveDofs;
266 :
267 : private:
268 : /// pointer to the domain
269 : ConstSmartPtr<TDomain> m_spDomain;
270 :
271 : /// pointer to the DofDistribution on the whole domain
272 : ConstSmartPtr<DoFDistribution> m_spDD;
273 :
274 : /// grouping for subset and non-conditional data
275 : struct NumberData
276 : {
277 : const static bool isConditional = false;
278 : const static size_t numFct = 1;
279 : typedef MathVector<1> value_type;
280 0 : NumberData(SmartPtr<UserData<number, dim> > functor_,
281 : std::string fctName_)
282 0 : : spFunctor(functor_), fctName(fctName_), bWholeDomain(true)
283 0 : {}
284 0 : NumberData(SmartPtr<UserData<number, dim> > functor_,
285 : std::string fctName_, std::string ssName_)
286 0 : : spFunctor(functor_), fctName(fctName_), ssName(ssName_),
287 0 : bWholeDomain(false)
288 0 : {}
289 : bool operator()(MathVector<1>& val, const MathVector<dim> x,
290 : number time, int si) const
291 : {
292 0 : (*spFunctor)(val[0], x, time, si); return true;
293 : }
294 :
295 : SmartPtr<UserData<number, dim> > spFunctor;
296 : std::string fctName;
297 : std::string ssName;
298 : size_t fct[numFct];
299 : SubsetGroup ssGrp;
300 : bool bWholeDomain;
301 : };
302 :
303 : /// grouping for subset and conditional data
304 : struct CondNumberData
305 : {
306 : const static bool isConditional = true;
307 : const static size_t numFct = 1;
308 : typedef MathVector<1> value_type;
309 0 : CondNumberData(SmartPtr<UserData<number, dim, bool> > functor_,
310 : std::string fctName_)
311 0 : : spFunctor(functor_), fctName(fctName_), bWholeDomain(true)
312 0 : {}
313 0 : CondNumberData(SmartPtr<UserData<number, dim, bool> > functor_,
314 : std::string fctName_, std::string ssName_)
315 0 : : spFunctor(functor_), fctName(fctName_), ssName(ssName_),
316 0 : bWholeDomain(true)
317 0 : {}
318 : bool operator()(MathVector<1>& val, const MathVector<dim> x,
319 : number time, int si) const
320 : {
321 0 : return (*spFunctor)(val[0], x, time, si);
322 : }
323 :
324 : SmartPtr<UserData<number, dim, bool> > spFunctor;
325 : std::string fctName;
326 : std::string ssName;
327 : size_t fct[numFct];
328 : SubsetGroup ssGrp;
329 : bool bWholeDomain;
330 : };
331 :
332 : /// grouping for subset and conditional data
333 : struct ConstNumberData
334 : {
335 : const static bool isConditional = false;
336 : const static size_t numFct = 1;
337 : typedef MathVector<1> value_type;
338 0 : ConstNumberData(number value_,
339 : std::string fctName_)
340 0 : : functor(value_), fctName(fctName_), bWholeDomain(true)
341 0 : {}
342 0 : ConstNumberData(number value_,
343 : std::string fctName_, std::string ssName_)
344 0 : : functor(value_), fctName(fctName_), ssName(ssName_),
345 0 : bWholeDomain(false)
346 0 : {}
347 : inline bool operator()(MathVector<1>& val, const MathVector<dim> x,
348 : number time, int si) const
349 : {
350 0 : val[0] = functor; return true;
351 : }
352 :
353 : number functor;
354 : std::string fctName;
355 : std::string ssName;
356 : size_t fct[numFct];
357 : SubsetGroup ssGrp;
358 : bool bWholeDomain;
359 : };
360 :
361 : /// grouping for subset and non-conditional data
362 : struct VectorData
363 : {
364 : const static bool isConditional = false;
365 : const static size_t numFct = dim;
366 : typedef MathVector<dim> value_type;
367 0 : VectorData(SmartPtr<UserData<MathVector<dim>, dim> > value_,
368 : std::string fctName_)
369 0 : : spFunctor(value_), fctName(fctName_), bWholeDomain(true)
370 0 : {}
371 0 : VectorData(SmartPtr<UserData<MathVector<dim>, dim> > value_,
372 : std::string fctName_, std::string ssName_)
373 0 : : spFunctor(value_), fctName(fctName_), ssName(ssName_),
374 0 : bWholeDomain(false)
375 0 : {}
376 : bool operator()(MathVector<dim>& val, const MathVector<dim> x,
377 : number time, int si) const
378 : {
379 0 : (*spFunctor)(val, x, time, si); return true;
380 : }
381 :
382 : SmartPtr<UserData<MathVector<dim>, dim> > spFunctor;
383 : std::string fctName;
384 : std::string ssName;
385 : size_t fct[numFct];
386 : SubsetGroup ssGrp;
387 : bool bWholeDomain;
388 : };
389 :
390 : std::vector<CondNumberData> m_vCondNumberData;
391 : std::vector<NumberData> m_vNumberData;
392 : std::vector<ConstNumberData> m_vConstNumberData;
393 :
394 : std::vector<VectorData> m_vVectorData;
395 :
396 : /// non-conditional obstacle values for all subsets
397 : std::map<int, std::vector<NumberData*> > m_mNumberObsSegment;
398 :
399 : /// constant obstacle values for all subsets
400 : std::map<int, std::vector<ConstNumberData*> > m_mConstNumberObsSegment;
401 :
402 : /// conditional obstacle values for all subsets
403 : std::map<int, std::vector<CondNumberData*> > m_mCondNumberObsSegment;
404 :
405 : /// non-conditional obstacle values for all subsets
406 : std::map<int, std::vector<VectorData*> > m_mVectorObsSegment;
407 : };
408 :
409 :
410 : } // end namespace ug
411 :
412 : // include implementation
413 : #include "obstacle_constraint_interface_impl.h"
414 :
415 : #endif /* __H__UG__LIB_ALGEBRA__OPERATOR__PRECONDITIONER__PROJECTED_GAUSS_SEIDEL__OBSTACLE_CONSTRAINT_INTERFACE__ */
416 :
|
