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__ELEM_DISC__ELEM_DISC_INTERFACE__
34 : #define __H__UG__LIB_DISC__SPATIAL_DISC__ELEM_DISC__ELEM_DISC_INTERFACE__
35 :
36 : // extern headers
37 : #include <vector>
38 : #include <string>
39 :
40 : // intern headers
41 : #include "lib_disc/common/local_algebra.h"
42 : #include "lib_disc/time_disc/solution_time_series.h"
43 : #include "lib_disc/function_spaces/approximation_space.h"
44 : #include "lib_disc/local_finite_element/local_finite_element_id.h"
45 : #include "lib_disc/reference_element/reference_element_traits.h"
46 : #include "lib_disc/spatial_disc/user_data/data_import.h"
47 : #include "common/util/provider.h"
48 : #include "lib_disc/domain_util.h"
49 : #include "lib_disc/domain_traits.h"
50 : #include "elem_modifier.h"
51 : #include "lib_disc/spatial_disc/elem_disc/err_est_data.h"
52 : #include "bridge/util_algebra_dependent.h"
53 : #include "lib_disc/common/multi_index.h"
54 :
55 : namespace ug{
56 :
57 : /// Types of elem disc
58 : enum ElemDiscType
59 : {
60 : EDT_NONE = 0,
61 : EDT_ELEM = 1 << 0,
62 : EDT_SIDE = 1 << 1,
63 : EDT_BND = 1 << 2,
64 : EDT_ALL = EDT_NONE | EDT_SIDE | EDT_ELEM | EDT_BND
65 : };
66 :
67 :
68 : /// Proxy struct for generic passing of any vector type
69 : struct VectorProxyBase
70 : {
71 : virtual ~VectorProxyBase() {};
72 : virtual number evaluate(const DoFIndex& di) const = 0;
73 : };
74 :
75 : template <typename TVector>
76 0 : struct VectorProxy : public VectorProxyBase
77 : {
78 0 : VectorProxy(const TVector& v) : m_v(v) {}
79 :
80 0 : virtual number evaluate(const DoFIndex& di) const {return DoFRef(m_v, di);}
81 :
82 : const TVector& m_v;
83 :
84 : };
85 :
86 : /**
87 : * Element Discretizations
88 : *
89 : * \defgroup lib_disc_elem_disc Elem Disc
90 : * \ingroup lib_discretization
91 : */
92 :
93 : /// \ingroup lib_disc_elem_disc
94 : /// @{
95 : /*
96 : template <typename TDomain>
97 : class IElemDiscBaseData
98 : {
99 : public:
100 : /// Domain type
101 : typedef TDomain domain_type;
102 :
103 : /// World dimension
104 : static const int dim = TDomain::dim;
105 : };
106 : */
107 : /// This class encapsulates all functions related to error estimation
108 : template <typename TLeaf, typename TDomain>
109 : class IElemAssembleFuncs
110 : {
111 : public:
112 : /// constructor
113 0 : IElemAssembleFuncs() { set_default_add_fct(); }
114 :
115 : /// Virtual destructor
116 0 : virtual ~IElemAssembleFuncs(){}
117 :
118 : /// Barton Nackman trick (TODO: needed?)
119 : typedef TLeaf leaf_type;
120 :
121 0 : TLeaf& asLeaf()
122 0 : { return static_cast<TLeaf&>(*this); }
123 :
124 : /// Domain type
125 : typedef TDomain domain_type;
126 :
127 : /// World dimension
128 : static const int dim = TDomain::dim;
129 :
130 : ////////////////////////////
131 : // assembling functions
132 : ////////////////////////////
133 : public:
134 : /// virtual prepares the loop over all elements of one type
135 0 : virtual void prep_assemble_loop() {}
136 :
137 : /// virtual prepares the loop over all elements of one type
138 0 : virtual void post_assemble_loop() {}
139 :
140 : /// prepare the time step
141 : virtual void prep_timestep(number future_time, number time, VectorProxyBase* u);
142 :
143 : /// prepare the time step element-wise
144 : virtual void prep_timestep_elem(const number time, const LocalVector& u, GridObject* elem, const MathVector<dim> vCornerCoords[]);
145 :
146 : /// virtual prepares the loop over all elements of one type
147 : virtual void prep_elem_loop(const ReferenceObjectID roid, const int si);
148 :
149 : /// virtual prepare one elements for assembling
150 : virtual void prep_elem(const LocalVector& u, GridObject* elem, const ReferenceObjectID roid, const MathVector<dim> vCornerCoords[]);
151 :
152 : /// virtual postprocesses the loop over all elements of one type
153 : virtual void fsh_elem_loop();
154 :
155 : /// finish the time step
156 : virtual void fsh_timestep(number time, VectorProxyBase* u);
157 :
158 : /// virtual finish the time step element-wise
159 : virtual void fsh_timestep_elem(const number time, const LocalVector& u, GridObject* elem, const MathVector<dim> vCornerCoords[]);
160 :
161 : /// Assembling of Jacobian (Stiffness part)
162 : virtual void add_jac_A_elem(LocalMatrix& J, const LocalVector& u, GridObject* elem, const MathVector<dim> vCornerCoords[]);
163 :
164 : /// Assembling of Jacobian (Mass part)
165 : virtual void add_jac_M_elem(LocalMatrix& J, const LocalVector& u, GridObject* elem, const MathVector<dim> vCornerCoords[]);
166 :
167 : /// virtual Assembling of Defect (Stiffness part)
168 : virtual void add_def_A_elem(LocalVector& d, const LocalVector& u, GridObject* elem, const MathVector<dim> vCornerCoords[]);
169 :
170 : /// defect for explicit terms
171 : virtual void add_def_A_expl_elem(LocalVector& d, const LocalVector& u, GridObject* elem, const MathVector<dim> vCornerCoords[]);
172 :
173 : /// virtual Assembling of Defect (Mass part)
174 : virtual void add_def_M_elem(LocalVector& d, const LocalVector& u, GridObject* elem, const MathVector<dim> vCornerCoords[]);
175 :
176 : /// virtual Assembling of Right-Hand Side
177 : virtual void add_rhs_elem(LocalVector& rhs, GridObject* elem, const MathVector<dim> vCornerCoords[]);
178 :
179 :
180 : /// function dispatching call to implementation
181 : /// \{
182 : void do_prep_timestep(number future_time, const number time, VectorProxyBase* u, size_t algebra_id);
183 : void do_prep_timestep_elem(const number time, LocalVector& u, GridObject* elem, const MathVector<dim> vCornerCoords[]);
184 : void do_prep_elem_loop(const ReferenceObjectID roid, const int si);
185 : void do_prep_elem(LocalVector& u, GridObject* elem, const ReferenceObjectID roid, const MathVector<dim> vCornerCoords[]);
186 : void do_fsh_elem_loop();
187 : void do_fsh_timestep(const number time, VectorProxyBase* u, size_t algebra_id);
188 : void do_fsh_timestep_elem(const number time, LocalVector& u, GridObject* elem, const MathVector<dim> vCornerCoords[]);
189 : void do_add_jac_A_elem(LocalMatrix& J, LocalVector& u, GridObject* elem, const MathVector<dim> vCornerCoords[]);
190 : void do_add_jac_M_elem(LocalMatrix& J, LocalVector& u, GridObject* elem, const MathVector<dim> vCornerCoords[]);
191 : void do_add_def_A_elem(LocalVector& d, LocalVector& u, GridObject* elem, const MathVector<dim> vCornerCoords[]);
192 : void do_add_def_A_expl_elem(LocalVector& d, LocalVector& u, GridObject* elem, const MathVector<dim> vCornerCoords[]);
193 : void do_add_def_M_elem(LocalVector& d, LocalVector& u, GridObject* elem, const MathVector<dim> vCornerCoords[]);
194 : void do_add_rhs_elem(LocalVector& rhs, GridObject* elem, const MathVector<dim> vCornerCoords[]);
195 :
196 :
197 :
198 :
199 : protected:
200 : // register the functions
201 : template <typename TAssFunc> void set_prep_timestep_fct(size_t algebra_id, TAssFunc func);
202 : template <typename TAssFunc> void set_prep_timestep_elem_fct(ReferenceObjectID id, TAssFunc func);
203 : template <typename TAssFunc> void set_fsh_timestep_fct(size_t algebra_id, TAssFunc func);
204 : template <typename TAssFunc> void set_fsh_timestep_elem_fct(ReferenceObjectID id, TAssFunc func);
205 :
206 : template <typename TAssFunc> void set_prep_elem_loop_fct(ReferenceObjectID id, TAssFunc func);
207 : template <typename TAssFunc> void set_prep_elem_fct(ReferenceObjectID id, TAssFunc func);
208 : template <typename TAssFunc> void set_fsh_elem_loop_fct(ReferenceObjectID id, TAssFunc func);
209 :
210 : template <typename TAssFunc> void set_add_jac_A_elem_fct(ReferenceObjectID id, TAssFunc func);
211 : template <typename TAssFunc> void set_add_jac_M_elem_fct(ReferenceObjectID id, TAssFunc func);
212 : template <typename TAssFunc> void set_add_def_A_elem_fct(ReferenceObjectID id, TAssFunc func);
213 : template <typename TAssFunc> void set_add_def_A_expl_elem_fct(ReferenceObjectID id, TAssFunc func);
214 : template <typename TAssFunc> void set_add_def_M_elem_fct(ReferenceObjectID id, TAssFunc func);
215 : template <typename TAssFunc> void set_add_rhs_elem_fct(ReferenceObjectID id, TAssFunc func);
216 :
217 :
218 :
219 : // unregister functions
220 : void remove_prep_timestep_fct(size_t algebra_id);
221 : void remove_prep_timestep_elem_fct(ReferenceObjectID id);
222 : void remove_fsh_timestep_fct(size_t algebra_id);
223 : void remove_fsh_timestep_elem_fct(ReferenceObjectID id);
224 :
225 : void remove_prep_elem_loop_fct(ReferenceObjectID id);
226 : void remove_prep_elem_fct(ReferenceObjectID id);
227 : void remove_fsh_elem_loop_fct(ReferenceObjectID id);
228 :
229 : void remove_add_jac_A_elem_fct(ReferenceObjectID id);
230 : void remove_add_jac_M_elem_fct(ReferenceObjectID id);
231 : void remove_add_def_A_elem_fct(ReferenceObjectID id);
232 : void remove_add_def_A_expl_elem_fct(ReferenceObjectID id);
233 : void remove_add_def_M_elem_fct(ReferenceObjectID id);
234 : void remove_add_rhs_elem_fct(ReferenceObjectID id);
235 :
236 : protected:
237 : /// sets all assemble functions to the corresponding virtual ones
238 : void set_default_add_fct();
239 :
240 : /// sets all assemble functions to NULL for a given ReferenceObjectID
241 : void clear_add_fct(ReferenceObjectID id);
242 :
243 : /// sets all assemble functions to NULL (for all ReferenceObjectID's)
244 : void clear_add_fct();
245 :
246 : private:
247 : // abbreviation for own type
248 : typedef IElemAssembleFuncs<TLeaf, TDomain> T;
249 :
250 : // types of timestep function pointers
251 : typedef void (T::*PrepareTimestepFct)(number, number, VectorProxyBase*);
252 : typedef void (T::*PrepareTimestepElemFct)(number, const LocalVector& u, GridObject* elem, const MathVector<dim> vCornerCoords[]);
253 : typedef void (T::*FinishTimestepFct)(number, VectorProxyBase*);
254 : typedef void (T::*FinishTimestepElemFct)(number, const LocalVector& u, GridObject* elem, const MathVector<dim> vCornerCoords[]);
255 :
256 : // types of loop function pointers
257 : typedef void (T::*PrepareElemLoopFct)(ReferenceObjectID roid, int si);
258 : typedef void (T::*PrepareElemFct)(const LocalVector& u, GridObject* elem, const ReferenceObjectID roid, const MathVector<dim> vCornerCoords[]);
259 : typedef void (T::*FinishElemLoopFct)();
260 :
261 : // types of Jacobian assemble functions
262 : typedef void (T::*ElemJAFct)(LocalMatrix& J, const LocalVector& u, GridObject* elem, const MathVector<dim> vCornerCoords[]);
263 : typedef void (T::*ElemJMFct)(LocalMatrix& J, const LocalVector& u, GridObject* elem, const MathVector<dim> vCornerCoords[]);
264 :
265 : // types of Defect assemble functions
266 : typedef void (T::*ElemdAFct)(LocalVector& d, const LocalVector& u, GridObject* elem, const MathVector<dim> vCornerCoords[]);
267 : typedef void (T::*ElemdMFct)(LocalVector& d, const LocalVector& u, GridObject* elem, const MathVector<dim> vCornerCoords[]);
268 :
269 : // types of right hand side assemble functions
270 : typedef void (T::*ElemRHSFct)(LocalVector& rhs, GridObject* elem, const MathVector<dim> vCornerCoords[]);
271 :
272 :
273 : private:
274 : // timestep function pointers
275 : PrepareTimestepFct m_vPrepareTimestepFct[bridge::NUM_ALGEBRA_TYPES];
276 : PrepareTimestepElemFct m_vPrepareTimestepElemFct[NUM_REFERENCE_OBJECTS];
277 : FinishTimestepFct m_vFinishTimestepFct[bridge::NUM_ALGEBRA_TYPES];
278 : FinishTimestepElemFct m_vFinishTimestepElemFct[NUM_REFERENCE_OBJECTS];
279 :
280 : // loop function pointers
281 : PrepareElemLoopFct m_vPrepareElemLoopFct[NUM_REFERENCE_OBJECTS];
282 : PrepareElemFct m_vPrepareElemFct[NUM_REFERENCE_OBJECTS];
283 : FinishElemLoopFct m_vFinishElemLoopFct[NUM_REFERENCE_OBJECTS];
284 :
285 : // Jacobian function pointers
286 : ElemJAFct m_vElemJAFct[NUM_REFERENCE_OBJECTS];
287 : ElemJMFct m_vElemJMFct[NUM_REFERENCE_OBJECTS];
288 :
289 : // Defect function pointers
290 : ElemdAFct m_vElemdAFct[NUM_REFERENCE_OBJECTS];
291 : ElemdAFct m_vElemdAExplFct[NUM_REFERENCE_OBJECTS];
292 : ElemdMFct m_vElemdMFct[NUM_REFERENCE_OBJECTS];
293 :
294 : // Rhs function pointers
295 : ElemRHSFct m_vElemRHSFct[NUM_REFERENCE_OBJECTS];
296 :
297 : public:
298 : /// sets the geometric object type
299 : /**
300 : * This functions set the geometric object type of the object, that is
301 : * assembled next. The user has to call this function before most of the
302 : * assembling routines can be called. Keep in mind, that the elements are
303 : * looped type by type, thus this function has to be called very few times.
304 : */
305 : void set_roid(ReferenceObjectID id, int discType);
306 :
307 : /// check, if all inputs have been set
308 : void check_roid(ReferenceObjectID roid, int discType);
309 :
310 : protected:
311 : /// current Geometric Object
312 : ReferenceObjectID m_roid;
313 : };
314 :
315 :
316 :
317 :
318 : /// This class encapsulates all functions related to error estimation
319 : template <typename TLeaf, typename TDomain>
320 : class IElemEstimatorFuncs
321 : {
322 : public:
323 : /// constructor
324 0 : IElemEstimatorFuncs() : m_bDoErrEst(false), m_spErrEstData(SPNULL)
325 0 : { set_default_add_fct(); }
326 :
327 : /// Virtual destructor
328 0 : virtual ~IElemEstimatorFuncs(){}
329 :
330 : /// Barton Nackman trick (TODO: needed?)
331 : typedef TLeaf leaf_type;
332 :
333 0 : TLeaf& asLeaf()
334 0 : { return static_cast<TLeaf&>(*this); }
335 :
336 : /// Domain type
337 : typedef TDomain domain_type;
338 :
339 : /// World dimension
340 : static const int dim = TDomain::dim;
341 :
342 : void do_prep_err_est_elem_loop(const ReferenceObjectID roid, const int si);
343 : void do_prep_err_est_elem(LocalVector& u, GridObject* elem, const MathVector<dim> vCornerCoords[]);
344 : void do_compute_err_est_A_elem(LocalVector& u, GridObject* elem, const MathVector<dim> vCornerCoords[], const number& scale);
345 : void do_compute_err_est_M_elem(LocalVector& u, GridObject* elem, const MathVector<dim> vCornerCoords[], const number& scale);
346 : void do_compute_err_est_rhs_elem(GridObject* elem, const MathVector<dim> vCornerCoords[], const number& scale);
347 : void do_fsh_err_est_elem_loop();
348 :
349 : public:
350 : /// virtual prepares the loop over all elements of one type for the computation of the error estimator
351 : virtual void prep_err_est_elem_loop(const ReferenceObjectID roid, const int si);
352 :
353 : /// virtual prepares the loop over all elements of one type for the computation of the error estimator
354 : virtual void prep_err_est_elem(const LocalVector& u, GridObject* elem, const MathVector<dim> vCornerCoords[]);
355 :
356 : /// virtual compute the error estimator (stiffness part) contribution for one element
357 : virtual void compute_err_est_A_elem(const LocalVector& u, GridObject* elem, const MathVector<dim> vCornerCoords[], const number& scale);
358 :
359 : /// virtual compute the error estimator (mass part) contribution for one element
360 : virtual void compute_err_est_M_elem(const LocalVector& u, GridObject* elem, const MathVector<dim> vCornerCoords[], const number& scale);
361 :
362 : /// virtual compute the error estimator (rhs part) contribution for one element
363 : virtual void compute_err_est_rhs_elem(GridObject* elem, const MathVector<dim> vCornerCoords[], const number& scale);
364 :
365 : /// virtual postprocesses the loop over all elements of one type in the computation of the error estimator
366 : virtual void fsh_err_est_elem_loop();
367 :
368 : protected:
369 : template <typename TAssFunc> void set_prep_err_est_elem_loop(ReferenceObjectID id, TAssFunc func);
370 : template <typename TAssFunc> void set_prep_err_est_elem(ReferenceObjectID id, TAssFunc func);
371 : template <typename TAssFunc> void set_compute_err_est_A_elem(ReferenceObjectID id, TAssFunc func);
372 : template <typename TAssFunc> void set_compute_err_est_M_elem(ReferenceObjectID id, TAssFunc func);
373 : template <typename TAssFunc> void set_compute_err_est_rhs_elem(ReferenceObjectID id, TAssFunc func);
374 : template <typename TAssFunc> void set_fsh_err_est_elem_loop(ReferenceObjectID id, TAssFunc func);
375 :
376 : void remove_prep_err_est_elem_loop(ReferenceObjectID id);
377 : void remove_prep_err_est_elem(ReferenceObjectID id);
378 : void remove_compute_err_est_A_elem(ReferenceObjectID id);
379 : void remove_compute_err_est_M_elem(ReferenceObjectID id);
380 : void remove_compute_err_est_rhs_elem(ReferenceObjectID id);
381 : void remove_fsh_err_est_elem_loop(ReferenceObjectID id);
382 :
383 : /// sets all assemble functions to NULL for a given ReferenceObjectID
384 : void clear_add_fct(ReferenceObjectID id);
385 :
386 : /// sets all assemble functions to NULL (for all ReferenceObjectID's)
387 : void clear_add_fct();
388 :
389 : /// sets all assemble functions to the corresponding virtual ones
390 : void set_default_add_fct();
391 :
392 :
393 :
394 : private:
395 : // abbreviation for own type
396 : typedef IElemEstimatorFuncs<TLeaf, TDomain> T;
397 :
398 : // types of the error estimator assembler
399 : typedef void (T::*PrepareErrEstElemLoopFct)(ReferenceObjectID roid, int si);
400 : typedef void (T::*PrepareErrEstElemFct)(const LocalVector& u, GridObject* elem, const MathVector<dim> vCornerCoords[]);
401 : typedef void (T::*ElemComputeErrEstAFct)(const LocalVector& u, GridObject* elem, const MathVector<dim> vCornerCoords[], const number&);
402 : typedef void (T::*ElemComputeErrEstMFct)(const LocalVector& u, GridObject* elem, const MathVector<dim> vCornerCoords[], const number&);
403 : typedef void (T::*ElemComputeErrEstRhsFct)(GridObject* elem, const MathVector<dim> vCornerCoords[], const number&);
404 : typedef void (T::*FinishErrEstElemLoopFct)();
405 :
406 : // Error estimator functions
407 : PrepareErrEstElemLoopFct m_vPrepareErrEstElemLoopFct[NUM_REFERENCE_OBJECTS];
408 : PrepareErrEstElemFct m_vPrepareErrEstElemFct[NUM_REFERENCE_OBJECTS];
409 : ElemComputeErrEstAFct m_vElemComputeErrEstAFct[NUM_REFERENCE_OBJECTS];
410 : ElemComputeErrEstMFct m_vElemComputeErrEstMFct[NUM_REFERENCE_OBJECTS];
411 : ElemComputeErrEstRhsFct m_vElemComputeErrEstRhsFct[NUM_REFERENCE_OBJECTS];
412 : FinishErrEstElemLoopFct m_vFinishErrEstElemLoopFct[NUM_REFERENCE_OBJECTS];
413 :
414 :
415 : // //////////////////////////
416 : // Error estimator
417 : // //////////////////////////
418 : public:
419 : /// sets the pointer to an error estimator data object (or NULL)
420 : /**
421 : * This function sets the pointer to an error estimator data object
422 : * that should be used for this discretization. Note that the ElemDisc
423 : * object must use RTTI to try to convert this pointer to the type
424 : * of the objects accepted by it for this purpose. If the conversion
425 : * fails than an exception must be thrown since this situation is not
426 : * allowed.
427 : */
428 0 : void set_error_estimator(SmartPtr<IErrEstData<TDomain> > ee) {m_spErrEstData = ee; m_bDoErrEst = true;}
429 :
430 : /// find out whether or not a posteriori error estimation is to be performed for this disc
431 0 : bool err_est_enabled() const {return m_bDoErrEst;}
432 :
433 : /// returns the pointer to the error estimator data object (or NULL)
434 0 : virtual SmartPtr<IErrEstData<TDomain> > err_est_data() {return m_spErrEstData;}
435 :
436 : private:
437 : /// flag indicating whether or not a posteriori error estimation is to be performed for this disc
438 : bool m_bDoErrEst;
439 :
440 : protected:
441 : /// error estimation object associated to the element discretization
442 : SmartPtr<IErrEstData<TDomain> > m_spErrEstData;
443 :
444 : public:
445 : /// sets the geometric object type
446 : /**
447 : * This functions set the geometric object type of the object, that is
448 : * assembled next. The user has to call this function before most of the
449 : * assembling routines can be called. Keep in mind, that the elements are
450 : * looped type by type, thus this function has to be called very few times.
451 : */
452 : void set_roid(ReferenceObjectID id, int discType);
453 :
454 : /// check, if all inputs have been set
455 : void check_roid(ReferenceObjectID roid, int discType);
456 :
457 : protected:
458 : /// current Geometric Object
459 : ReferenceObjectID m_roid;
460 : };
461 :
462 : /// base class for all element-wise discretizations
463 : /**
464 : * This class is the base class for element-wise discretizations. An
465 : * implementation of this class must provide local stiffness/mass-matrix
466 : * contribution of one element to the global jacobian and local contributions
467 : * of one element to the local defect.
468 : */
469 :
470 : template <typename TDomain>
471 : class IElemDiscBase
472 :
473 : {
474 : public:
475 : /// Domain type
476 : typedef TDomain domain_type;
477 :
478 : /// Position type
479 : typedef typename TDomain::position_type position_type;
480 :
481 : /// World dimension
482 : static const int dim = TDomain::dim;
483 :
484 : public:
485 : /// Constructor
486 : IElemDiscBase(const char* functions = "", const char* subsets = "");
487 :
488 : /// Constructor
489 : IElemDiscBase(const std::vector<std::string>& vFct, const std::vector<std::string>& vSubset);
490 :
491 : /// Virtual destructor
492 0 : virtual ~IElemDiscBase(){}
493 :
494 : public:
495 : /// sets the approximation space
496 : /** Calls protected virtual 'approximation_space_changed', when a new approximation space
497 : * has been set. Note that 'approximation_space_changed' is only called once if the
498 : * same approximation space is set multiple times.*/
499 : void set_approximation_space(SmartPtr<ApproximationSpace<TDomain> > approxSpace);
500 :
501 : /// returns approximation space
502 0 : SmartPtr<ApproximationSpace<TDomain> > approx_space() {return m_spApproxSpace;}
503 :
504 : /// returns approximation space
505 0 : ConstSmartPtr<ApproximationSpace<TDomain> > approx_space() const {return m_spApproxSpace;}
506 :
507 : /// returns the domain
508 0 : SmartPtr<TDomain> domain(){return m_spApproxSpace->domain();}
509 :
510 : /// returns the domain
511 0 : ConstSmartPtr<TDomain> domain() const{return m_spApproxSpace->domain();}
512 :
513 : /// returns the subset handler
514 0 : typename TDomain::subset_handler_type& subset_handler()
515 : {
516 : UG_ASSERT(m_spApproxSpace.valid(), "ApproxSpace not set.");
517 0 : return *m_spApproxSpace->domain()->subset_handler();
518 : }
519 :
520 : /// returns the subset handler
521 0 : const typename TDomain::subset_handler_type& subset_handler() const
522 : {
523 : UG_ASSERT(m_spApproxSpace.valid(), "ApproxSpace not set.");
524 0 : return *m_spApproxSpace->domain()->subset_handler();
525 : }
526 : /*
527 : void add_elem_modifier(SmartPtr<IElemDiscModifier<TDomain> > elemModifier )
528 : {
529 : m_spElemModifier.push_back(elemModifier);
530 : elemModifier->set_elem_disc(this);
531 : }
532 : std::vector<SmartPtr<IElemDiscModifier<TDomain> > >& get_elem_modifier(){ return m_spElemModifier;}
533 : */
534 :
535 : protected:
536 : /// callback invoked, when approximation space is changed
537 0 : virtual void approximation_space_changed() {}
538 :
539 : /// Approximation Space
540 : SmartPtr<ApproximationSpace<TDomain> > m_spApproxSpace;
541 :
542 : /// Approximation Space
543 : // std::vector<SmartPtr<IElemDiscModifier<TDomain> > > m_spElemModifier;
544 :
545 : ////////////////////////////
546 : // Functions and Subsets
547 : ////////////////////////////
548 : public:
549 : /// sets functions by name list, divided by ','
550 : void set_functions(const std::string& functions);
551 :
552 : /// sets functions by vector of names
553 : void set_functions(const std::vector<std::string>& functions);
554 :
555 : /// sets subset(s) by name list, divided by ','
556 : void set_subsets(const std::string& subsets);
557 :
558 : /// sets subset(s) by name list, divided by ','
559 : void set_subsets(const std::vector<std::string>& subsets);
560 :
561 : /// number of functions this discretization handles
562 0 : size_t num_fct() const {return m_vFct.size();}
563 :
564 : /// returns the symbolic functions
565 0 : const std::vector<std::string>& symb_fcts() const {return m_vFct;}
566 :
567 : /// number of subsets this discretization handles
568 0 : size_t num_subsets() const {return m_vSubset.size();}
569 :
570 : /// returns the symbolic subsets
571 0 : const std::vector<std::string>& symb_subsets() const {return m_vSubset;}
572 :
573 : /// returns the current function pattern
574 0 : ConstSmartPtr<FunctionPattern> function_pattern() const {return m_spFctPattern;}
575 :
576 : /// returns the current function group
577 0 : const FunctionGroup& function_group() const {return m_fctGrp;}
578 :
579 : /// returns the current function index mapping
580 0 : const FunctionIndexMapping& map() const {return m_fctIndexMap;}
581 :
582 : /// checks the setup of the elem disc
583 : void check_setup(bool bNonRegularGrid);
584 :
585 : protected:
586 : /// vector holding name of all symbolic functions
587 : std::vector<std::string> m_vFct;
588 :
589 : /// vector holding name of all symbolic subsets
590 : std::vector<std::string> m_vSubset;
591 :
592 : /// current function pattern
593 : ConstSmartPtr<FunctionPattern> m_spFctPattern;
594 :
595 : /// current function group
596 : FunctionGroup m_fctGrp;
597 :
598 : /// current function index mapping
599 : FunctionIndexMapping m_fctIndexMap;
600 :
601 : /// sets current function pattern
602 : void set_function_pattern(ConstSmartPtr<FunctionPattern> fctPatt);
603 :
604 : /// updates the function index mapping
605 : void update_function_index_mapping();
606 :
607 : ////////////////////////////
608 : // UserData and Coupling
609 : ////////////////////////////
610 : public:
611 : /// registers a data import
612 : void register_import(IDataImport<dim>& Imp);
613 :
614 : /// returns number of imports
615 0 : size_t num_imports() const {return m_vIImport.size();}
616 :
617 : /// returns an import
618 0 : IDataImport<dim>& get_import(size_t i)
619 : {
620 : UG_ASSERT(i < num_imports(), "Invalid index");
621 0 : return *m_vIImport[i];
622 : }
623 :
624 : /// removes all imports
625 0 : void clear_imports() {m_vIImport.clear();}
626 :
627 : protected:
628 : /// data imports
629 : std::vector<IDataImport<dim>*> m_vIImport;
630 :
631 : ////////////////////////////
632 : // time handling
633 : ////////////////////////////
634 : public:
635 : /// sets if assembling should be time-dependent and the local time series
636 : /**
637 : * This function specifies if the assembling is time-dependent. If NULL is
638 : * passed, the assembling is assumed to be time-independent. If a local
639 : * time series is passed, this series is used as previous solution.
640 : *
641 : * \param[in] locTimeSeries Time series of previous solutions
642 : */
643 : void set_time_dependent(LocalVectorTimeSeries& locTimeSeries,
644 : const std::vector<number>& vScaleMass,
645 : const std::vector<number>& vScaleStiff);
646 :
647 : /// sets that the assembling is time independent
648 : void set_time_independent();
649 :
650 : /// returns if assembling is time-dependent
651 0 : bool is_time_dependent() const {return (m_pLocalVectorTimeSeries != NULL) && !m_bStationaryForced;}
652 :
653 : /// sets that the assembling is always stationary (even in instationary case)
654 0 : void set_stationary(bool bStationaryForced = true) {m_bStationaryForced = bStationaryForced;}
655 0 : void set_stationary() {set_stationary(true);}
656 :
657 : /// returns if local time series needed by assembling
658 : /**
659 : * This callback must be implemented by a derived Elem Disc in order to handle
660 : * time-dependent data. As return the derived Elem Disc can specify, if
661 : * it really needs data from previous time steps for the (spatial) disc. The
662 : * default is false.
663 : *
664 : * \returns if elem disc needs time series local solutions
665 : */
666 0 : virtual bool requests_local_time_series() {return false;}
667 :
668 0 : bool local_time_series_needed() {return is_time_dependent() && requests_local_time_series();}
669 :
670 : /// sets the current time point
671 0 : void set_time_point(const size_t timePoint) {m_timePoint = timePoint;}
672 :
673 : /// returns the currently considered time point of the time-disc scheme
674 0 : size_t time_point() const {return m_timePoint;}
675 :
676 : /// returns currently set timepoint
677 0 : number time() const
678 : {
679 0 : if(m_pLocalVectorTimeSeries) return m_pLocalVectorTimeSeries->time(m_timePoint);
680 : else return 0.0;
681 : }
682 :
683 : /// returns the local time solutions
684 : /**
685 : * This function returns the local time solutions. This a type of vector,
686 : * that holds the local unknowns for each time point of a time series.
687 : * Note, that the first sol is always the current (iterated, unknown)
688 : * time point, while all remaining sols are the already computed time steps
689 : * used e.g. in a time stepping scheme.
690 : *
691 : * \returns vLocalTimeSol vector of local time Solutions
692 : */
693 0 : const LocalVectorTimeSeries* local_time_solutions() const
694 0 : {return m_pLocalVectorTimeSeries;}
695 :
696 : /// returns the weight factors of the time-disc scheme
697 : /// \{
698 0 : const std::vector<number>& mass_scales() const {return m_vScaleMass;}
699 0 : const std::vector<number>& stiff_scales() const {return m_vScaleStiff;}
700 :
701 0 : number mass_scale(const size_t timePoint) const {return m_vScaleMass[timePoint];}
702 0 : number stiff_scale(const size_t timePoint) const {return m_vScaleStiff[timePoint];}
703 :
704 0 : number mass_scale() const {return m_vScaleMass[m_timePoint];}
705 0 : number stiff_scale() const {return m_vScaleStiff[m_timePoint];}
706 : /// \}
707 :
708 : protected:
709 : /// time point
710 : size_t m_timePoint;
711 :
712 : /// list of local vectors for all solutions of the time series
713 : LocalVectorTimeSeries* m_pLocalVectorTimeSeries;
714 :
715 : /// weight factors for time dependent assembling
716 : /// \{
717 : std::vector<number> m_vScaleMass;
718 : std::vector<number> m_vScaleStiff;
719 : /// \}
720 :
721 : /// flag if stationary assembling is to be used even in instationary assembling
722 : bool m_bStationaryForced;
723 :
724 :
725 :
726 : ////////////////////////////
727 : // general info
728 : ////////////////////////////
729 : public:
730 : /// returns the type of elem disc
731 0 : virtual int type() const {return EDT_ELEM | EDT_SIDE;}
732 :
733 : /// requests assembling for trial spaces and grid type
734 : /**
735 : * This function is called before the assembling starts. The
736 : * IElemDisc-Implementation is supposed to checks if it can assemble the set
737 : * of LFEID and the grid type. It may register corresponding assembling
738 : * functions or perform other initialization.
739 : * If the ElemDisc does not support the setting it should throw an exception.
740 : *
741 : * \param[in] vLfeID vector of Local Finite Element IDs
742 : * \param[in] bNonRegularGrid regular grid type
743 : */
744 : virtual void prepare_setting(const std::vector<LFEID>& vLfeID, bool bNonRegularGrid) = 0;
745 :
746 : /// returns if discretization acts on hanging nodes if present
747 : /**
748 : * This function returns if a discretization really needs the hanging nodes
749 : * in case of non-regular grid. This may not be the case for e.g. finite
750 : * element assemblings but is needed for finite volumes
751 : */
752 0 : virtual bool use_hanging() const {return false;}
753 : };
754 :
755 :
756 : template <typename TDomain>
757 : class IElemError :
758 : public IElemDiscBase<TDomain>,
759 : public IElemEstimatorFuncs<IElemDisc<TDomain>, TDomain>
760 : {
761 : public:
762 : typedef TDomain domain_type;
763 : static const int dim = TDomain::dim;
764 :
765 : friend class IElemEstimatorFuncs<IElemDisc<TDomain>, TDomain>;
766 : typedef IElemEstimatorFuncs<IElemDisc<TDomain>, TDomain> estimator_base_type;
767 :
768 0 : IElemError(const char* functions, const char* subsets)
769 0 : : IElemDiscBase<TDomain>(functions, subsets), estimator_base_type() {}
770 :
771 0 : IElemError(const std::vector<std::string>& vFct, const std::vector<std::string>& vSubset)
772 0 : : IElemDiscBase<TDomain>(vFct, vSubset), estimator_base_type() {}
773 :
774 :
775 : protected:
776 :
777 : /// sets all assemble functions to NULL for a given ReferenceObjectID
778 0 : void clear_add_fct(ReferenceObjectID id)
779 0 : { estimator_base_type::clear_add_fct(id); }
780 :
781 : /// sets all assemble functions to NULL (for all ReferenceObjectID's)
782 0 : void clear_add_fct()
783 0 : { estimator_base_type::clear_add_fct(); }
784 :
785 : /// sets all assemble functions to the corresponding virtual ones
786 : //void set_default_add_fct();
787 : using estimator_base_type::set_default_add_fct;
788 :
789 : };
790 :
791 :
792 : /**
793 : * Element discretization (including error indicator)
794 : * TODO: Should be separated!!!
795 : * */
796 : template <typename TDomain>
797 : class IElemDisc :
798 : public IElemAssembleFuncs<IElemDisc<TDomain>, TDomain>,
799 : public IElemError<TDomain>
800 : {
801 : public:
802 : typedef TDomain domain_type;
803 : static const int dim = TDomain::dim;
804 :
805 : /// real base class
806 : typedef IElemError<TDomain> base_type;
807 : typedef IElemEstimatorFuncs<IElemDisc<TDomain>, TDomain> estimator_base_type;
808 : typedef IElemAssembleFuncs<IElemDisc<TDomain>, TDomain> assemble_base_type;
809 :
810 : friend class IElemEstimatorFuncs<IElemDisc<TDomain>, TDomain>;
811 : friend class IElemAssembleFuncs<IElemDisc<TDomain>, TDomain>;
812 :
813 :
814 0 : IElemDisc(const char* functions, const char* subsets)
815 0 : : assemble_base_type(), IElemError<TDomain>(functions, subsets) {}
816 :
817 0 : IElemDisc(const std::vector<std::string>& vFct, const std::vector<std::string>& vSubset)
818 0 : : assemble_base_type(), IElemError<TDomain>(vFct, vSubset) {}
819 :
820 : protected:
821 :
822 : /// sets all assemble functions to NULL for a given ReferenceObjectID
823 0 : void clear_add_fct(ReferenceObjectID id)
824 : {
825 : base_type::clear_add_fct(id);
826 0 : assemble_base_type::clear_add_fct(id);
827 0 : }
828 :
829 : /// sets all assemble functions to NULL (for all ReferenceObjectID's)
830 0 : void clear_add_fct()
831 : {
832 : base_type::clear_add_fct();
833 0 : assemble_base_type::clear_add_fct();
834 0 : }
835 :
836 : /// sets all assemble functions to the corresponding virtual ones
837 0 : void set_default_add_fct()
838 : {
839 0 : base_type::set_default_add_fct();
840 0 : assemble_base_type::set_default_add_fct();
841 0 : }
842 :
843 :
844 : public:
845 0 : void add_elem_modifier(SmartPtr<IElemDiscModifier<TDomain> > elemModifier )
846 : {
847 0 : m_spElemModifier.push_back(elemModifier);
848 : elemModifier->set_elem_disc(this);
849 0 : }
850 :
851 0 : std::vector<SmartPtr<IElemDiscModifier<TDomain> > >& get_elem_modifier()
852 0 : { return m_spElemModifier;}
853 :
854 : protected:
855 : /// Approximation Space
856 : std::vector<SmartPtr<IElemDiscModifier<TDomain> > > m_spElemModifier;
857 :
858 :
859 : };
860 :
861 :
862 :
863 : /// @}
864 :
865 : } // end namespace ug
866 :
867 : #include "elem_disc_interface_impl.h"
868 :
869 : #endif /* __H__UG__LIB_DISC__SPATIAL_DISC__ELEM_DISC__ELEM_DISC_INTERFACE__ */
|
