Line data Source code
1 : /*
2 : * Copyright (c) 2014-2015: G-CSC, Goethe University Frankfurt
3 : * Authors: Dmitry Logashenko, Markus Breit
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 : /*
34 : * Data shared by element discretizations for a-posteriori error estimation
35 : */
36 :
37 : #ifndef __H__UG__LIB_DISC__SPATIAL_DISC__ELEM_DISC__ERR_EST_DATA__
38 : #define __H__UG__LIB_DISC__SPATIAL_DISC__ELEM_DISC__ERR_EST_DATA__
39 :
40 : // extern headers
41 : #include <vector>
42 : #include <string>
43 : #include <limits>
44 :
45 : // intern headers
46 : #include "lib_grid/tools/surface_view.h"
47 : #include "lib_grid/algorithms/multi_grid_util.h"
48 : #include "lib_disc/function_spaces/integrate.h"
49 :
50 : #ifdef UG_PARALLEL
51 : #include "lib_grid/parallelization/util/compol_attachment_reduce.h"
52 : #include "lib_grid/parallelization/util/compol_copy_attachment.h"
53 : #endif
54 :
55 : #include <boost/mpl/for_each.hpp>
56 :
57 :
58 : namespace ug{
59 :
60 : /// Base class for error estimator data
61 : /**
62 : * This virtual class should be the base of any particular error estimator
63 : * implemented in the elem_disc's. Every elem_disc class (not object!) should
64 : * declare its own derived class for keeping intermediate information that
65 : * should be accumulated in the computation of the local error estimators.
66 : * Several objects of the elem_disc class may share the same object of the
67 : * derived class for a consistent computation of the error estimator.
68 : */
69 : template <typename TDomain>
70 : class IErrEstData
71 : {
72 : public:
73 : /// world dimension
74 : static const int dim = TDomain::dim;
75 :
76 : /// class constructor
77 0 : IErrEstData () : m_consider(true), m_scale(1.0) {};
78 :
79 : /// virtual class destructor
80 : virtual ~IErrEstData () {};
81 :
82 : /// virtual function to allocate data structures for the error estimator
83 : virtual void alloc_err_est_data (ConstSmartPtr<SurfaceView> spSV, const GridLevel& gl) = 0;
84 :
85 : /// virtual function called after the computation of the error estimator data in all the elements
86 : virtual void summarize_err_est_data (SmartPtr<TDomain> spDomain) = 0;
87 :
88 : /// calculate L2 integrals
89 : virtual number get_elem_error_indicator(GridObject* elem, const MathVector<dim> vCornerCoords[]) = 0;
90 :
91 : /// virtual function to release data structures for the error estimator
92 : virtual void release_err_est_data () = 0;
93 :
94 : /// virtual function granting get access to the m_consider member
95 0 : bool consider_me() const {return m_consider;};
96 :
97 : /// whether or not this instance is to be considered by domainDisc
98 : /**
99 : * The domainDisc calls alloc_err_est_data(), summarize_err_est_data(),
100 : * get_elem_error_indicator() etc. only for ErrEstData objects that have
101 : * consider_me() == true.
102 : * This is useful when using a MultipleErrEstData object combining ErrEstData objects
103 : * which are already set to some ElemDisc. In this case, set_consider_me(false).
104 : * The default value is true.
105 : */
106 0 : void set_consider_me(bool b) {m_consider = b;};
107 :
108 : /// set scaling factor for final error calculation
109 0 : void set_scaling_factor(number scale) {m_scale = scale;};
110 :
111 : /// get scaling factor
112 : /**
113 : * This factor is used in the element-wise calculation of error indicators.
114 : * After calculation of indicators for each IErrEstData object in the domain
115 : * discretization (typically one per equation), the final overall indicator
116 : * is calculated as weighted sum (with the scaling factors as weights).
117 : * The default value (if not set) is 1.0.
118 : *
119 : * @return scale scaling factor
120 : */
121 0 : number scaling_factor() {return m_scale;}
122 :
123 : private:
124 : bool m_consider;
125 : number m_scale;
126 : };
127 :
128 : /// Error estimator data class storing one scalar number per side
129 : /**
130 : * This class allocates an attachment keeping one number per full-dimensional
131 : * element side. Furthermore, the data are collected at the boundaries of the
132 : * patches (in the case of the adaptive refinement).
133 : *
134 : * \tparam TDomain domain type
135 : */
136 : template <typename TDomain>
137 : class SideFluxErrEstData : public IErrEstData<TDomain>
138 : {
139 : public:
140 : /// domain type
141 : typedef TDomain domain_type;
142 :
143 : /// world dimension
144 : static const int dim = TDomain::dim;
145 :
146 : /// type of the sides (face, edge) and the elems (volume, face)
147 : typedef typename domain_traits<dim>::side_type side_type;
148 :
149 : public:
150 : /// constructor
151 0 : SideFluxErrEstData() : IErrEstData<TDomain>() {};
152 :
153 : /// virtual class destructor
154 0 : virtual ~SideFluxErrEstData() {};
155 :
156 : // Functions to access data
157 :
158 : /// get the data reference for a given side
159 : number& operator()
160 : (
161 : side_type* pSide ///< pointer to the side
162 : )
163 : {
164 : return m_aaFluxJump[pSide];
165 : };
166 :
167 : /// get the surface view
168 : ConstSmartPtr<SurfaceView>& surface_view() {return m_spSV;};
169 :
170 : // Interface virtual functions inherited from IErrEstData
171 :
172 : /// virtual function to allocate data structures for the error estimator
173 : virtual void alloc_err_est_data (ConstSmartPtr<SurfaceView> spSV, const GridLevel& gl);
174 :
175 : /// virtual function called after the computation of the error estimator data in all the elements
176 : virtual void summarize_err_est_data (SmartPtr<TDomain> spDomain);
177 :
178 : /// calculate L2 integrals
179 0 : virtual number get_elem_error_indicator(GridObject* elem, const MathVector<dim> vCornerCoords[]) {return 0;};
180 :
181 : /// virtual function to release data structures of the error estimator
182 : virtual void release_err_est_data ();
183 :
184 : private:
185 : /// Flux jumps for the error estimator
186 : ANumber m_aFluxJumpOverSide;
187 :
188 : /// Attachment accessor
189 : MultiGrid::AttachmentAccessor<side_type, ANumber> m_aaFluxJump;
190 :
191 : /// Grid for the attachment
192 : ConstSmartPtr<SurfaceView> m_spSV;
193 :
194 : /// Finest grid level
195 : GridLevel m_errEstGL;
196 : };
197 :
198 :
199 :
200 : /// Error estimator data class storing a number vector per side and per element.
201 : /**
202 : * This class represents an H1 error estimator.
203 : * It can integrate expressions on elements and their sides with arbitrary order.
204 : * A vector of values at defined integration points is attached to any element and
205 : * side to that end.
206 : *
207 : * RECOMMENDED (INTENDED) USAGE
208 : * The data will typically consist of the values of certain functions at integration
209 : * points (IP) on the sides and the element.
210 : * A pointer to an object of this class can be handed to any element discretization
211 : * involved in a discretization. They will access the attachments in their method
212 : * compute_err_est_elem and add their respective parts of the function to be
213 : * integrated for the error estimator. Exactly one of them (or maybe some other object)
214 : * then has to do the actual integration using the given values at the IPs and add up
215 : * side and element terms according to the error estimator formula used.
216 : *
217 : * \tparam TDomain domain type
218 : */
219 : template <typename TDomain>
220 : class SideAndElemErrEstData : public IErrEstData<TDomain>
221 : {
222 : public:
223 : /// domain type
224 : typedef TDomain domain_type;
225 :
226 : /// world dimension
227 : static const int dim = TDomain::dim;
228 :
229 : /// type of the sides (face, edge) and the elems (volume, face)
230 : typedef typename domain_traits<dim>::side_type side_type;
231 : typedef typename domain_traits<dim>::element_type elem_type;
232 :
233 : /// attachment type
234 : typedef Attachment<std::vector<number> > attachment_type;
235 :
236 : /// this class
237 : typedef SideAndElemErrEstData<TDomain> this_type;
238 :
239 : /// maximal number of sides of any element
240 : static const int MAX_NUM_SIDES = 8;
241 :
242 : public:
243 : /// constructors
244 : SideAndElemErrEstData(size_t _sideOrder, size_t _elemOrder, const char* subsets);
245 : SideAndElemErrEstData(size_t _sideOrder, size_t _elemOrder,
246 : std::vector<std::string> subsets = std::vector<std::string>(0));
247 :
248 : /// virtual class destructor
249 0 : virtual ~SideAndElemErrEstData() {};
250 :
251 : // Functions to access data
252 :
253 : /// getting the side integration order
254 0 : size_t side_order() const {return sideOrder;}
255 :
256 : /// getting the elem integration order
257 0 : size_t elem_order() const {return elemOrder;}
258 :
259 : /// get the data reference for a given side and ip
260 : number& operator()
261 : (
262 : side_type* pSide, ///< pointer to the side
263 : size_t ip ///< integration point id on the side
264 : );
265 :
266 : /// get the data reference for a given elem and ip
267 : number& operator()
268 : (
269 : elem_type* pElem, ///< pointer to the elem
270 : size_t ip ///< integration point id on the elem
271 : );
272 :
273 : /// get the local side integration points for a specific roid
274 : template <int refDim>
275 : const MathVector<refDim>* side_local_ips(const ReferenceObjectID roid);
276 :
277 : /// get the local elem integration points for a specific roid
278 : template <int refDim>
279 : const MathVector<refDim>* elem_local_ips(const ReferenceObjectID roid);
280 :
281 : /// get all global side integration points
282 : MathVector<TDomain::dim>* all_side_global_ips(GridObject* elem, const MathVector<dim> vCornerCoords[]);
283 :
284 : /// get the global side integration points for a specific side roid
285 : MathVector<TDomain::dim>* side_global_ips(GridObject* elem, const MathVector<dim> vCornerCoords[]);
286 :
287 : /// get the global elem integration points for a specific roid
288 : MathVector<TDomain::dim>* elem_global_ips(GridObject* elem, const MathVector<dim> vCornerCoords[]);
289 :
290 : /// get number of side IPs of a specific side
291 : size_t num_side_ips(const side_type* pSide);
292 :
293 : /// get number of side IPs of a specific side type
294 : size_t num_side_ips(const ReferenceObjectID roid);
295 :
296 : /// get number of first IP belonging to a specific side
297 : size_t first_side_ips(const ReferenceObjectID roid, const size_t side);
298 :
299 : /// get number of side IPs
300 : size_t num_all_side_ips(const ReferenceObjectID roid);
301 :
302 : /// get number of elem IPs
303 : size_t num_elem_ips(const ReferenceObjectID roid);
304 :
305 : /// get index of specific side IP in sideIP array returned by side_local_ips
306 : size_t side_ip_index(const ReferenceObjectID roid, const size_t side, const size_t ip);
307 :
308 : /// get the surface view
309 : ConstSmartPtr<SurfaceView>& surface_view () {return m_spSV;};
310 :
311 : // virtual functions inherited from IErrEstData
312 : /// virtual function to allocate data structures for the error estimator
313 : virtual void alloc_err_est_data (ConstSmartPtr<SurfaceView> spSV, const GridLevel& gl);
314 :
315 : /// virtual function called after the computation of the error estimator data in all the elements
316 : virtual void summarize_err_est_data (SmartPtr<TDomain> spDomain);
317 :
318 : /// calculate L2 integrals
319 : virtual number get_elem_error_indicator(GridObject* elem, const MathVector<dim> vCornerCoords[]);
320 :
321 : /// virtual function to release data structures of the error estimator
322 : virtual void release_err_est_data ();
323 :
324 : /// select L2/H1 Estimator
325 0 : void set_type(int type) {m_type = (type<=H1_ERROR_TYPE) ? type : H1_ERROR_TYPE;}
326 :
327 : protected:
328 : /// initialization of quadrature (to be called during construction)
329 : void init_quadrature();
330 :
331 : /// helper struct for getting quadrature rules by use of mpl::lists
332 : template<int refDim>
333 : struct GetQuadRules
334 : {
335 : GetQuadRules(QuadratureRule<refDim>** ppQuadRule, size_t quadOrder) :
336 : m_ppQuadRule(ppQuadRule), m_quadOrder(quadOrder) {}
337 : QuadratureRule<refDim>** m_ppQuadRule;
338 : size_t m_quadOrder;
339 : template< typename TElem > void operator()(TElem&)
340 : {
341 : const ReferenceObjectID roid = geometry_traits<TElem>::REFERENCE_OBJECT_ID;
342 0 : m_ppQuadRule[roid] =
343 0 : const_cast<QuadratureRule<refDim>*>(&QuadratureRuleProvider<refDim>::get(roid, m_quadOrder));
344 : }
345 : };
346 :
347 : private:
348 : /// order of side and elem function approximations for integrating
349 : size_t sideOrder;
350 : size_t elemOrder;
351 :
352 : /// the subsets this error estimator will produce values for
353 : std::vector<std::string> m_vSs;
354 : SubsetGroup m_ssg;
355 :
356 : /// storage for integration rules
357 : QuadratureRule<dim-1>* quadRuleSide[NUM_REFERENCE_OBJECTS];
358 : QuadratureRule<dim>* quadRuleElem[NUM_REFERENCE_OBJECTS];
359 :
360 : /// extra storage for local side IPs (elem IPs are contained in elem quad rules)
361 : std::vector<MathVector<TDomain::dim> > m_SideIPcoords[NUM_REFERENCE_OBJECTS];
362 :
363 : /// storage for global elem and side IPs
364 : std::vector<MathVector<TDomain::dim> > m_sideGlobalIPcoords;
365 : std::vector<MathVector<TDomain::dim> > m_singleSideGlobalIPcoords; // not the most elegant solution...
366 : std::vector<MathVector<TDomain::dim> > m_elemGlobalIPcoords;
367 :
368 : /// the first index for IPs of a specific side in the sideIP series for a roid
369 : size_t m_sideIPsStartIndex[NUM_REFERENCE_OBJECTS][MAX_NUM_SIDES];
370 :
371 : /// vector of attachments for sides
372 : attachment_type m_aSide;
373 :
374 : /// vector of attachments for elems
375 : attachment_type m_aElem;
376 :
377 : /// vector of side attachment accessors
378 : MultiGrid::AttachmentAccessor<side_type, attachment_type > m_aaSide;
379 :
380 : /// vector of elem attachment accessors
381 : MultiGrid::AttachmentAccessor<elem_type, attachment_type > m_aaElem;
382 :
383 : /// Grid for the attachment
384 : ConstSmartPtr<SurfaceView> m_spSV;
385 :
386 : /// Finest grid level
387 : GridLevel m_errEstGL;
388 :
389 : /// set type
390 : enum type {L2_ERROR_TYPE=0, H1_ERROR_TYPE};
391 : int m_type;
392 : };
393 :
394 :
395 : /// Error estimator data class for discretizations with more than one unknown.
396 : /**
397 : * This class is a kind of wrapper for a bundle of error estimator objects.
398 : * It can be useful if a discretization depends on more than one unknown and needs to
399 : * compute error estimators for both of them:
400 : * One can only pass one error estimator to any ElemDisc, but at the same time
401 : * it is necessary to calculate the errors for different unknowns separately!
402 : *
403 : * This class will not actually do anything, but pass any request on to the underlying
404 : * objects.
405 : * As they most probably figure in some other ElemDisc of their corresponding unknown,
406 : * the considerMe property is set to false by default.
407 : *
408 : * Make sure that the error estimation routines of any elem disc that is given this
409 : * MultipleErrEstData object write to the correct sub-objects (i.e. ErrEstData objects)!
410 : * There has to be some kind of mapping between the unknowns of the elem disc and
411 : * the order in which the single-function ErrEstData objects are added to this object.
412 : * This will typically be the exact same order. So try not to add the same object of
413 : * MultipleErrEstData to elem discs with unknowns defined in a different order!
414 : *
415 : * The template parameter TErrEstData must implement the IErrEstData interface.
416 : *
417 : * \todo Maybe find a better way to deal with the orders of ErrEstData objects here
418 : * and unknowns in the elem discs.
419 : * \tparam TDomain domain type
420 : */
421 : template <typename TDomain, typename TErrEstData>
422 : class MultipleErrEstData : public IErrEstData<TDomain>
423 : {
424 : public:
425 : /// world dimension
426 : static const int dim = TDomain::dim;
427 :
428 : /// class constructor
429 0 : MultipleErrEstData(ConstSmartPtr<ApproximationSpace<TDomain> > approx)
430 : : IErrEstData<TDomain>(), m_spApprox(approx),
431 0 : m_fctGrp(m_spApprox->dof_distribution_info())
432 : {
433 : this->set_consider_me(false);
434 0 : }
435 :
436 : /// virtual class destructor
437 0 : virtual ~MultipleErrEstData() {};
438 :
439 : /// adding error estimator data objects
440 0 : virtual void add(SmartPtr<TErrEstData> spEed, const char* fct)
441 : {
442 : // check that fct is not already contained in fctGrp
443 : size_t uid = m_spApprox->fct_id_by_name(fct);
444 0 : if (m_fctGrp.contains(uid))
445 : {
446 0 : UG_THROW("Error estimator for function '" << fct << "' can not be added\n"
447 : "as another error estimator object for the same function is already\n"
448 : "contained here.")
449 : }
450 :
451 : // add to function group
452 : try
453 : {
454 0 : m_fctGrp.add(fct);
455 : }
456 0 : UG_CATCH_THROW("Error estimator data object for function '"
457 : << fct << "' could not be added.");
458 :
459 0 : m_vEed.push_back(spEed.get());
460 0 : }
461 :
462 : /// getting the number of underlying error estimator data objects
463 : size_t num() const {return m_vEed.size();};
464 :
465 : /// accessing the underlying error estimator data objects via function id
466 0 : TErrEstData* get(size_t uid)
467 : {
468 0 : if (!m_fctGrp.contains(uid))
469 : {
470 : std::string name = m_spApprox->name(uid);
471 0 : UG_THROW("Trying to access error estimator data object "
472 : "for unique function index " << uid << "(aka '" << name << "')\n"
473 : "which is not present in this collection.")
474 : }
475 :
476 0 : return m_vEed[m_fctGrp.local_index(uid)];
477 : }
478 :
479 : // inherited from IErrEstData
480 : /// virtual function to allocate data structures for the error estimator
481 : virtual void alloc_err_est_data(ConstSmartPtr<SurfaceView> spSV, const GridLevel& gl);
482 :
483 : /// virtual function called after the computation of the error estimator data in all the elements
484 : virtual void summarize_err_est_data(SmartPtr<TDomain> spDomain);
485 :
486 : /// calculate L2 integrals
487 : virtual number get_elem_error_indicator(GridObject* elem, const MathVector<dim> vCornerCoords[]);
488 :
489 : /// virtual function to release data structures for the error estimator
490 : virtual void release_err_est_data();
491 :
492 : protected:
493 : std::vector<TErrEstData*> m_vEed;
494 :
495 : /// approx space
496 : ConstSmartPtr<ApproximationSpace<TDomain> > m_spApprox;
497 :
498 : /// function group (in order to map fcts to error estimator objects)
499 : FunctionGroup m_fctGrp;
500 : };
501 :
502 :
503 : template <typename TDomain>
504 : class MultipleSideAndElemErrEstData
505 : : public MultipleErrEstData<TDomain, SideAndElemErrEstData<TDomain> >
506 : {
507 : public:
508 : /// world dimension
509 : static const int dim = TDomain::dim;
510 :
511 : /// type of the sides (face, edge) and the elems (volume, face)
512 : typedef typename SideAndElemErrEstData<TDomain>::side_type side_type;
513 : typedef typename SideAndElemErrEstData<TDomain>::elem_type elem_type;
514 :
515 : /// constructor
516 0 : MultipleSideAndElemErrEstData(ConstSmartPtr<ApproximationSpace<TDomain> > approx)
517 : : MultipleErrEstData<TDomain, SideAndElemErrEstData<TDomain> >(approx),
518 0 : m_bEqSideOrder(false), m_bEqElemOrder(false) {};
519 :
520 : /// destructor
521 0 : virtual ~MultipleSideAndElemErrEstData() {};
522 :
523 : /// adding error estimator data objects
524 : /// overrides parent add method; performs check for equal order after adding
525 : virtual void add(SmartPtr<SideAndElemErrEstData<TDomain> > spEed, const char* fct);
526 :
527 : /// returns whether all underlying err ests have the same elem and side integration orders
528 : bool equal_side_order() const {return m_bEqSideOrder;}
529 :
530 : /// returns whether all underlying err ests have the same elem and side integration orders
531 : bool equal_elem_order() const {return m_bEqElemOrder;}
532 :
533 : protected:
534 : /// find out whether all underlying err_ests have the same integration orders
535 : /// (makes assembling easier)
536 : void check_equal_order();
537 :
538 : /// find out whether all underlying err_ests have the same side integration orders
539 : /// (makes assembling easier)
540 : void check_equal_side_order();
541 :
542 : /// find out whether all underlying err_ests have the same elem integration orders
543 : /// (makes assembling easier)
544 : void check_equal_elem_order();
545 :
546 : private:
547 : bool m_bEqSideOrder;
548 : bool m_bEqElemOrder;
549 :
550 : };
551 :
552 : } // end of namespace ug
553 :
554 : #include "err_est_data_impl.h"
555 :
556 : #endif /* __H__UG__LIB_DISC__SPATIAL_DISC__ELEM_DISC__ERR_EST_DATA__ */
557 :
558 : /* End of File */
|
