Line data Source code
1 : /*
2 : * Copyright (c) 2010-2015: G-CSC, Goethe University Frankfurt
3 : * Author: Andreas Vogel
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_DISC__SPATIAL_DISC__CONSTRAINTS__LAGRANGE_DIRICHLET_BOUNDARY__
34 : #define __H__UG__LIB_DISC__SPATIAL_DISC__CONSTRAINTS__LAGRANGE_DIRICHLET_BOUNDARY__
35 :
36 : #include "lib_disc/common/function_group.h"
37 : #include "lib_disc/common/groups_util.h"
38 : #include "lib_disc/spatial_disc/domain_disc_interface.h"
39 : #include "lib_disc/function_spaces/approximation_space.h"
40 : #include "lib_disc/spatial_disc/user_data/const_user_data.h"
41 : #include "lib_grid/tools/subset_handler_interface.h"
42 :
43 : #include "lib_disc/spatial_disc/constraints/constraint_interface.h"
44 :
45 : #include <map>
46 : #include <vector>
47 :
48 :
49 : // #define LAGRANGE_DIRICHLET_ADJ_TRANSFER_FIX
50 :
51 : namespace ug{
52 :
53 : template < typename TDomain, typename TAlgebra>
54 : class DirichletBoundary
55 : : public IDomainConstraint<TDomain, TAlgebra>
56 : {
57 : public:
58 : /// Base Type
59 : typedef IDomainConstraint<TDomain, TAlgebra> base_type;
60 :
61 : /// Type of domain
62 : typedef TDomain domain_type;
63 :
64 : /// world Dimension
65 : static const int dim = domain_type::dim;
66 :
67 : /// Type of position coordinates (e.g. position_type)
68 : typedef typename domain_type::position_type position_type;
69 :
70 : /// Type of algebra
71 : typedef TAlgebra algebra_type;
72 :
73 : /// Type of algebra matrix
74 : typedef typename algebra_type::matrix_type matrix_type;
75 :
76 : /// Type of algebra vector
77 : typedef typename algebra_type::vector_type vector_type;
78 :
79 : /// Type of value type
80 : typedef typename matrix_type::value_type value_type;
81 :
82 : /// error estimator type
83 : typedef MultipleSideAndElemErrEstData<TDomain> err_est_type;
84 :
85 :
86 : public:
87 : /// If you want to have Dirichlet columns than use the constructor with the flag
88 : /// The default ist without Dirichlet columns.
89 :
90 : /// constructor
91 0 : DirichletBoundary()
92 0 : : m_bInvertSubsetSelection(false),
93 0 : m_bDirichletColumns(false),
94 0 : m_A(NULL)
95 : #ifdef LAGRANGE_DIRICHLET_ADJ_TRANSFER_FIX
96 : , m_bAdjustTransfers(true)
97 : #endif
98 0 : {clear();}
99 :
100 : /// constructor with flag for Dirichlet-Columns.
101 0 : DirichletBoundary(bool DirichletColumns)
102 0 : : m_bInvertSubsetSelection(false),
103 0 : m_bDirichletColumns(DirichletColumns),
104 0 : m_A(NULL)
105 : #ifdef LAGRANGE_DIRICHLET_ADJ_TRANSFER_FIX
106 : , m_bAdjustTransfers(true)
107 : #endif
108 0 : {clear();}
109 :
110 : #ifdef LAGRANGE_DIRICHLET_ADJ_TRANSFER_FIX
111 : /// constructor with flag for Dirichlet-Columns.
112 : DirichletBoundary(bool DirichletColumns, bool bAdjustTransfers)
113 : : m_bInvertSubsetSelection(false),
114 : m_bDirichletColumns(DirichletColumns),
115 : m_A(NULL),
116 : m_bAdjustTransfers(bAdjustTransfers)
117 :
118 : {clear();}
119 : #endif
120 :
121 : /// destructor
122 0 : ~DirichletBoundary() {}
123 :
124 : /// adds a lua callback (cond and non-cond)
125 : #ifdef UG_FOR_LUA
126 : void add(const char* name, const char* function, const char* subsets);
127 : void add(const char* name, const std::vector<std::string>& Fcts, const std::vector<std::string>& Subsets);
128 : void add(LuaFunctionHandle fct, const char* function, const char* subsets);
129 : void add(LuaFunctionHandle fct, const std::vector<std::string>& Fcts, const std::vector<std::string>& Subsets);
130 : #endif
131 :
132 : /// adds a conditional user-defined value as dirichlet condition for a function on subsets
133 : void add(SmartPtr<UserData<number, dim, bool> > func, const char* function, const char* subsets);
134 : void add(SmartPtr<UserData<number, dim, bool> > func, const std::vector<std::string>& Fcts, const std::vector<std::string>& Subsets);
135 :
136 : /// adds a user-defined value as dirichlet condition for a function on subsets
137 : void add(SmartPtr<UserData<number, dim> > func, const char* function, const char* subsets);
138 : void add(SmartPtr<UserData<number, dim> > func, const std::vector<std::string>& Fcts, const std::vector<std::string>& Subsets);
139 :
140 : /// adds a constant value as dirichlet condition for a function on subsets
141 : void add(number value, const char* function, const char* subsets);
142 : void add(number value, const std::vector<std::string>& Fcts, const std::vector<std::string>& Subsets);
143 :
144 : /// adds a user-defined vector as dirichlet condition for a vector-function on subsets
145 : void add(SmartPtr<UserData<MathVector<dim>, dim> > func, const char* functions, const char* subsets);
146 : void add(SmartPtr<UserData<MathVector<dim>, dim> > func, const std::vector<std::string>& Fcts, const std::vector<std::string>& Subsets);
147 :
148 : /// adds a number Dirichlet condition whose value is taken from the solution vector itself
149 : void add(const char* functions, const char* subsets);
150 : void add(const std::vector<std::string>& Fcts, const std::vector<std::string>& Subsets);
151 :
152 : /// inverts the subset selection making the conditions be imposed on the rest of the domain
153 0 : void invert_subset_selection() {m_bInvertSubsetSelection = true;};
154 :
155 : /// sets the approximation space to work on
156 : void set_approximation_space(SmartPtr<ApproximationSpace<TDomain> > approxSpace);
157 :
158 : /// removes all scheduled dirichlet data.
159 : void clear();
160 :
161 : /// Sets dirichlet rows for all registered dirichlet values
162 : /** (implemented by Mr. Xylouris and Mr. Reiter)
163 : *
164 : * This method is just an attempt to allow to set dirichlet rows in a matrix.
165 : * It should probably be a virtual method derived from IDirichletPostProcess.
166 : *
167 : * Note that adjust_jacobian does the same (...!!!)
168 : * You should thus use adjust_jacobian.
169 : *
170 : * This method is probably removed in the near future!
171 : *
172 : * It could make sense to keep it but implement it as an overload of
173 : * adjust_jacobian...
174 : *
175 : * If Mr. Vogel decides that this is nonsense, he may of course remove it!!!
176 : */
177 : void assemble_dirichlet_rows(matrix_type& mat, ConstSmartPtr<DoFDistribution> dd, number time = 0.0);
178 :
179 : public:
180 : ///////////////////////////////
181 : // Implement Interface
182 : ///////////////////////////////
183 :
184 : /// sets a unity row for all dirichlet indices
185 : void adjust_jacobian(matrix_type& J, const vector_type& u,
186 : ConstSmartPtr<DoFDistribution> dd, int type, number time = 0.0,
187 : ConstSmartPtr<VectorTimeSeries<vector_type> > vSol = NULL,
188 : const number s_a0 = 1.0);
189 :
190 : /// sets a zero value in the defect for all dirichlet indices
191 : void adjust_defect(vector_type& d, const vector_type& u,
192 : ConstSmartPtr<DoFDistribution> dd, int type, number time = 0.0,
193 : ConstSmartPtr<VectorTimeSeries<vector_type> > vSol = NULL,
194 : const std::vector<number>* vScaleMass = NULL,
195 : const std::vector<number>* vScaleStiff = NULL);
196 :
197 : /// sets the dirichlet value in the solution for all dirichlet indices
198 : void adjust_solution(vector_type& u,
199 : ConstSmartPtr<DoFDistribution> dd, int type, number time = 0.0);
200 :
201 : /// sets zero to correction
202 : virtual void adjust_correction(vector_type& c, ConstSmartPtr<DoFDistribution> dd, int type,
203 : number time = 0.0);
204 :
205 : /// sets unity rows in A and dirichlet values in right-hand side b
206 : void adjust_linear(matrix_type& A, vector_type& b,
207 : ConstSmartPtr<DoFDistribution> dd, int type, number time = 0.0);
208 :
209 : /// sets the dirichlet value in the right-hand side
210 : void adjust_rhs(vector_type& b, const vector_type& u,
211 : ConstSmartPtr<DoFDistribution> dd, int type, number time = 0.0);
212 :
213 : /// @copydoc IConstraint::adjust_error()
214 : virtual void adjust_error(const vector_type& u, ConstSmartPtr<DoFDistribution> dd, int type,
215 : number time = 0.0,
216 : ConstSmartPtr<VectorTimeSeries<vector_type> > vSol = SPNULL,
217 : const std::vector<number>* vScaleMass = NULL,
218 : const std::vector<number>* vScaleStiff = NULL);
219 :
220 : /// sets constraints in prolongation
221 : virtual void adjust_prolongation(matrix_type& P,
222 : ConstSmartPtr<DoFDistribution> ddFine,
223 : ConstSmartPtr<DoFDistribution> ddCoarse,
224 : int type,
225 : number time = 0.0);
226 :
227 : /// sets constraints in restriction
228 : virtual void adjust_restriction(matrix_type& R,
229 : ConstSmartPtr<DoFDistribution> ddCoarse,
230 : ConstSmartPtr<DoFDistribution> ddFine,
231 : int type,
232 : number time = 0.0);
233 :
234 : /// returns the type of the constraints
235 0 : virtual int type() const {return CT_DIRICHLET;}
236 :
237 : protected:
238 : void check_functions_and_subsets(FunctionGroup& functionGroup,
239 : SubsetGroup& subsetGroup, size_t numFct) const;
240 :
241 : void extract_data();
242 :
243 : template <typename TUserData, typename TScheduledUserData>
244 : void extract_data(std::map<int, std::vector<TUserData*> >& mvUserDataBndSegment,
245 : std::vector<TScheduledUserData>& vUserData);
246 :
247 : template <typename TUserData>
248 : void adjust_jacobian(const std::map<int, std::vector<TUserData*> >& mvUserData,
249 : matrix_type& J, const vector_type& u,
250 : ConstSmartPtr<DoFDistribution> dd, number time);
251 :
252 : template <typename TBaseElem, typename TUserData>
253 : void adjust_jacobian(const std::vector<TUserData*>& vUserData, int si,
254 : matrix_type& J, const vector_type& u,
255 : ConstSmartPtr<DoFDistribution> dd, number time);
256 :
257 : template <typename TUserData>
258 : void adjust_defect(const std::map<int, std::vector<TUserData*> >& mvUserData,
259 : vector_type& d, const vector_type& u,
260 : ConstSmartPtr<DoFDistribution> dd, number time);
261 :
262 : template <typename TBaseElem, typename TUserData>
263 : void adjust_defect(const std::vector<TUserData*>& vUserData, int si,
264 : vector_type& d, const vector_type& u,
265 : ConstSmartPtr<DoFDistribution> dd, number time);
266 :
267 : template <typename TUserData>
268 : void adjust_solution(const std::map<int, std::vector<TUserData*> >& mvUserData,
269 : vector_type& u, ConstSmartPtr<DoFDistribution> dd, number time);
270 :
271 : template <typename TUserData>
272 : void adjust_correction(const std::map<int, std::vector<TUserData*> >& mvUserData,
273 : vector_type& c, ConstSmartPtr<DoFDistribution> dd, number time);
274 :
275 : template <typename TBaseElem, typename TUserData>
276 : void adjust_solution(const std::vector<TUserData*>& vUserData, int si,
277 : vector_type& u, ConstSmartPtr<DoFDistribution> dd, number time);
278 :
279 : template <typename TBaseElem, typename TUserData>
280 : void adjust_correction(const std::vector<TUserData*>& vUserData, int si,
281 : vector_type& c, ConstSmartPtr<DoFDistribution> dd, number time);
282 :
283 : template <typename TUserData>
284 : void adjust_linear(const std::map<int, std::vector<TUserData*> >& mvUserData,
285 : matrix_type& A, vector_type& b,
286 : ConstSmartPtr<DoFDistribution> dd, number time);
287 :
288 : template <typename TBaseElem, typename TUserData>
289 : void adjust_linear(const std::vector<TUserData*>& vUserData, int si,
290 : matrix_type& A, vector_type& b,
291 : ConstSmartPtr<DoFDistribution> dd, number time);
292 :
293 : template <typename TUserData>
294 : void adjust_rhs(const std::map<int, std::vector<TUserData*> >& mvUserData,
295 : vector_type& b, const vector_type& u,
296 : ConstSmartPtr<DoFDistribution> dd, number time);
297 :
298 : template <typename TBaseElem, typename TUserData>
299 : void adjust_rhs(const std::vector<TUserData*>& vUserData, int si,
300 : vector_type& b, const vector_type& u,
301 : ConstSmartPtr<DoFDistribution> dd, number time);
302 :
303 : template <typename TUserData>
304 : void adjust_error(const std::map<int, std::vector<TUserData*> >& mvUserData,
305 : const vector_type& u, ConstSmartPtr<DoFDistribution> dd, number time);
306 :
307 : template <typename TUserData>
308 : void adjust_prolongation(const std::map<int, std::vector<TUserData*> >& mvUserData,
309 : matrix_type& P,
310 : ConstSmartPtr<DoFDistribution> ddFine,
311 : ConstSmartPtr<DoFDistribution> ddCoarse,
312 : number time);
313 :
314 : template <typename TBaseElem, typename TUserData>
315 : void adjust_prolongation(const std::vector<TUserData*>& vUserData, int si,
316 : matrix_type& P,
317 : ConstSmartPtr<DoFDistribution> ddFine,
318 : ConstSmartPtr<DoFDistribution> ddCoarse,
319 : number time);
320 :
321 : template <typename TUserData>
322 : void adjust_restriction(const std::map<int, std::vector<TUserData*> >& mvUserData,
323 : matrix_type& R,
324 : ConstSmartPtr<DoFDistribution> ddCoarse,
325 : ConstSmartPtr<DoFDistribution> ddFine,
326 : number time);
327 :
328 : template <typename TBaseElem, typename TUserData>
329 : void adjust_restriction(const std::vector<TUserData*>& vUserData, int si,
330 : matrix_type& R,
331 : ConstSmartPtr<DoFDistribution> ddCoarse,
332 : ConstSmartPtr<DoFDistribution> ddFine,
333 : number time);
334 :
335 : protected:
336 : /// grouping for subset and non-conditional data
337 : struct NumberData
338 : {
339 : const static bool isConditional = false;
340 : const static bool setSolValue = true;
341 : const static size_t numFct = 1;
342 : typedef MathVector<1> value_type;
343 0 : NumberData(SmartPtr<UserData<number, dim> > functor_,
344 : std::string fctName_, std::string ssName_)
345 0 : : spFunctor(functor_), fctName(fctName_), ssName(ssName_)
346 0 : {}
347 :
348 : bool operator()(MathVector<1>& val, const MathVector<dim> x,
349 : number time, int si) const
350 : {
351 0 : (*spFunctor)(val[0], x, time, si); return true;
352 : }
353 :
354 : SmartPtr<UserData<number, dim> > spFunctor;
355 : std::string fctName;
356 : std::string ssName;
357 : size_t fct[numFct];
358 : SubsetGroup ssGrp;
359 : };
360 :
361 : /// grouping for subset and conditional data
362 : struct CondNumberData
363 : {
364 : const static bool isConditional = true;
365 : const static bool setSolValue = true;
366 : const static size_t numFct = 1;
367 : typedef MathVector<1> value_type;
368 0 : CondNumberData(SmartPtr<UserData<number, dim, bool> > functor_,
369 : std::string fctName_, std::string ssName_)
370 0 : : spFunctor(functor_), fctName(fctName_), ssName(ssName_)
371 0 : {}
372 : bool operator()(MathVector<1>& val, const MathVector<dim> x,
373 : number time, int si) const
374 : {
375 0 : return (*spFunctor)(val[0], x, time, si);
376 : }
377 :
378 : SmartPtr<UserData<number, dim, bool> > spFunctor;
379 : std::string fctName;
380 : std::string ssName;
381 : size_t fct[numFct];
382 : SubsetGroup ssGrp;
383 : };
384 :
385 : /// grouping for subset and constant data
386 : struct ConstNumberData
387 : {
388 : const static bool isConditional = false;
389 : const static bool setSolValue = true;
390 : const static size_t numFct = 1;
391 : typedef MathVector<1> value_type;
392 0 : ConstNumberData(number value_,
393 : std::string fctName_, std::string ssName_)
394 0 : : functor(value_), fctName(fctName_), ssName(ssName_)
395 0 : {}
396 : inline bool operator()(MathVector<1>& val, const MathVector<dim> x,
397 : number time, int si) const
398 : {
399 0 : val[0] = functor; return true;
400 : }
401 :
402 : number functor;
403 : std::string fctName;
404 : std::string ssName;
405 : size_t fct[numFct];
406 : SubsetGroup ssGrp;
407 : };
408 :
409 : /// grouping for subset and non-conditional vector data
410 : struct VectorData
411 : {
412 : const static bool isConditional = false;
413 : const static bool setSolValue = true;
414 : const static size_t numFct = dim;
415 : typedef MathVector<dim> value_type;
416 0 : VectorData(SmartPtr<UserData<MathVector<dim>, dim> > value_,
417 : std::string fctName_, std::string ssName_)
418 0 : : spFunctor(value_), fctName(fctName_), ssName(ssName_)
419 0 : {}
420 : bool operator()(MathVector<dim>& val, const MathVector<dim> x,
421 : number time, int si) const
422 : {
423 0 : (*spFunctor)(val, x, time, si); return true;
424 : }
425 :
426 : SmartPtr<UserData<MathVector<dim>, dim> > spFunctor;
427 : std::string fctName;
428 : std::string ssName;
429 : size_t fct[numFct];
430 : SubsetGroup ssGrp;
431 : };
432 :
433 : /// grouping for subset and the data already stored in the solution
434 : struct OldNumberData
435 : {
436 : const static bool isConditional = false;
437 : const static bool setSolValue = false;
438 : const static size_t numFct = 1;
439 : typedef MathVector<1> value_type;
440 0 : OldNumberData(std::string fctName_, std::string ssName_)
441 0 : : fctName(fctName_), ssName(ssName_)
442 0 : {}
443 : inline bool operator()(MathVector<1>& val, const MathVector<dim> x,
444 : number time, int si) const
445 : {
446 : return true; // note that we do not set val because setSolValue == false
447 : }
448 :
449 : number functor;
450 : std::string fctName;
451 : std::string ssName;
452 : size_t fct[numFct];
453 : SubsetGroup ssGrp;
454 : };
455 :
456 : std::vector<CondNumberData> m_vBNDNumberData;
457 : std::vector<NumberData> m_vNumberData;
458 : std::vector<ConstNumberData> m_vConstNumberData;
459 :
460 : std::vector<VectorData> m_vVectorData;
461 :
462 : std::vector<OldNumberData> m_vOldNumberData;
463 :
464 : /// non-conditional boundary values for all subsets
465 : std::map<int, std::vector<NumberData*> > m_mNumberBndSegment;
466 :
467 : /// constant boundary values for all subsets
468 : std::map<int, std::vector<ConstNumberData*> > m_mConstNumberBndSegment;
469 :
470 : /// conditional boundary values for all subsets
471 : std::map<int, std::vector<CondNumberData*> > m_mBNDNumberBndSegment;
472 :
473 : /// non-conditional boundary values for all subsets
474 : std::map<int, std::vector<VectorData*> > m_mVectorBndSegment;
475 :
476 : /// non-conditional boundary values for all subsets
477 : std::map<int, std::vector<OldNumberData*> > m_mOldNumberBndSegment;
478 :
479 : protected:
480 : /// flag for inverting the subset selection: use Dirichlet throughout except for the given subsets
481 : bool m_bInvertSubsetSelection;
482 :
483 : /// flag for setting dirichlet columns
484 : bool m_bDirichletColumns;
485 :
486 : /// maps a column dirichlet index to the
487 : /// row and its corresponding matrix entry.
488 : std::map<int, std::map<int, std::map<int, value_type> > > m_dirichletMap;
489 : ///
490 : matrix_type* m_A;
491 :
492 : /// current ApproxSpace
493 : SmartPtr<ApproximationSpace<TDomain> > m_spApproxSpace;
494 :
495 : /// current domain
496 : SmartPtr<TDomain> m_spDomain;
497 :
498 : /// current position accessor
499 : typename domain_type::position_accessor_type m_aaPos;
500 : #ifdef LAGRANGE_DIRICHLET_ADJ_TRANSFER_FIX
501 : /// flag for setting dirichlet columns
502 : bool m_bAdjustTransfers;
503 : #endif
504 : };
505 :
506 : } // end namespace ug
507 :
508 : #include "lagrange_dirichlet_boundary_impl.h"
509 :
510 : #endif /* __H__UG__LIB_DISC__SPATIAL_DISC__CONSTRAINTS__LAGRANGE_DIRICHLET_BOUNDARY__ */
|
