Loading [MathJax]/extensions/tex2jax.js
ug4
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages
error_indicator_util.h
Go to the documentation of this file.
1/*
2 * Copyright (c) 2014-2015: G-CSC, Goethe University Frankfurt
3 * Author: Sebastian Reiter
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_DISC__ERROR_INDICATOR_UTIL__
34#define __H__UG_DISC__ERROR_INDICATOR_UTIL__
35
36#include "lib_grid/multi_grid.h"
39
40namespace ug{
41
42
44
52template<typename TElem>
56 number& min, number& max, number& sum, number& errSq, size_t& numElem,
57 number& minLocal, number& maxLocal, number& sumLocal, size_t& numElemLocal
58)
59{
60 typedef typename DoFDistribution::traits<TElem>::const_iterator const_iterator;
61
62// reset maximum of error
63 max = 0.0, min = std::numeric_limits<number>::max();
64 sum = 0.0;
65 errSq = 0.0;
66 numElem = 0;
67
68// get element iterator
69 const_iterator iter = dd->template begin<TElem>();
70 const_iterator iterEnd = dd->template end<TElem>();
71
72// loop all elements to find the maximum of the error
73 for (; iter != iterEnd; ++iter)
74 {
75 // get element
76 TElem* elem = *iter;
77
78 const number elemEta2 = aaError2[elem];
79
80 // if no error value exists: ignore (might be newly added by refinement);
81 // newly added elements are supposed to have a negative error estimator
82 if (aaError2[elem] < 0) UG_THROW("error value invalid!");//continue;
83
84
85 const number elemEta = sqrt(aaError2[elem]);
86
87 // search for maximum and minimum
88 if (elemEta > max) max = elemEta;
89 if (elemEta < min) min = elemEta;
90
91 // sum up total error
92 sum += elemEta;
93 errSq += elemEta2;
94 ++numElem;
95 }
96
97 // set local variables
98 maxLocal = max;
99 minLocal = min;
100 sumLocal = sum;
101#ifdef UG_PARALLEL
102 number errLocal = errSq;
103#endif
104 numElemLocal = numElem;
105
106#ifdef UG_PARALLEL
107 if (pcl::NumProcs() > 1)
108 {
110 max = com.allreduce(maxLocal, PCL_RO_MAX);
111 min = com.allreduce(minLocal, PCL_RO_MIN);
112 sum = com.allreduce(sumLocal, PCL_RO_SUM);
113 errSq = com.allreduce(errLocal, PCL_RO_SUM);
114 numElem = com.allreduce(numElemLocal, PCL_RO_SUM);
115 }
116#endif
117 UG_LOG(" +++++ Error indicator on " << numElem << " elements +++++\n");
118 UG_LOG(" +++ Element errors: maxEta=" << max << ", minEta="
119 << min << ", sumEta=" << sum << ", sumEtaSq=" << errSq << ".\n");
120
121 return (sum/numElem);
122}
123
124
126
133template<typename TElem>
137 number& min, number& max, number& totalErr, size_t& numElem,
138 number& minLocal, number& maxLocal, number& totalErrLocal, size_t& numElemLocal
139)
140{
141 typedef typename DoFDistribution::traits<TElem>::const_iterator const_iterator;
142
143// reset maximum of error
144 max = 0.0, min = std::numeric_limits<number>::max();
145 totalErr = 0.0;
146 numElem = 0;
147
148// get element iterator
149 const_iterator iter = dd->template begin<TElem>();
150 const_iterator iterEnd = dd->template end<TElem>();
151
152// loop all elements to find the maximum of the error
153 for (; iter != iterEnd; ++iter)
154 {
155 // get element
156 TElem* elem = *iter;
157
158 // if no error value exists: ignore (might be newly added by refinement);
159 // newly added elements are supposed to have a negative error estimator
160 if (aaError[elem] < 0) continue;
161
162 // search for maximum and minimum
163 if (aaError[elem] > max) max = aaError[elem];
164 if (aaError[elem] < min) min = aaError[elem];
165
166 // sum up total error
167 totalErr += aaError[elem];
168 ++numElem;
169 }
170
171 // set local variables
172 maxLocal = max;
173 minLocal = min;
174 totalErrLocal = totalErr;
175 numElemLocal = numElem;
176
177#ifdef UG_PARALLEL
178 if (pcl::NumProcs() > 1)
179 {
181 max = com.allreduce(maxLocal, PCL_RO_MAX);
182 min = com.allreduce(minLocal, PCL_RO_MIN);
183 totalErr = com.allreduce(totalErrLocal, PCL_RO_SUM);
184 numElem = com.allreduce(numElemLocal, PCL_RO_SUM);
185 }
186#endif
187 UG_LOG(" +++++ Error indicator on " << numElem << " elements +++++\n");
188 UG_LOG(" +++ Element errors: maxEtaSq=" << max << ", minEtaSq="
189 << min << ", sumEtaSq=" << totalErr << ".\n");
190}
191
192
194template<typename TElem>
198 number& min, number& max, number& totalErr, size_t& numElem
199
200)
201{
202 number minLocal, maxLocal, totalErrLocal;
203 size_t numElemLocal;
204 ComputeMinMax(aaError2, dd, min, max, totalErr, numElem, minLocal, maxLocal, totalErrLocal, numElemLocal);
205}
207
221template<typename TElem>
223 IRefiner& refiner,
225 number TOL,
226 number refineFrac, number coarseFrac,
227 int maxLevel)
228{
229 typedef typename DoFDistribution::traits<TElem>::const_iterator const_iterator;
230
231// compute minimal/maximal/ total error and number of elements
232 number min, max, totalErr;
233 size_t numElem;
234 ComputeMinMaxTotal(aaError, dd, min, max, totalErr, numElem);
235
236// check if total error is smaller than tolerance. If that is the case we're done
237 if(totalErr < TOL)
238 {
239 UG_LOG(" +++ Total error "<<totalErr<<" smaller than TOL ("<<TOL<<"). done.");
240 return;
241 }
242
243// Compute minimum
244 number minErrToRefine = max * refineFrac;
245 UG_LOG(" +++ Refining elements if error greater " << refineFrac << "*" << max <<
246 " = " << minErrToRefine << ".\n");
247 number maxErrToCoarse = min * (1 + coarseFrac);
248 if(maxErrToCoarse < TOL/numElem) maxErrToCoarse = TOL/numElem;
249 UG_LOG(" +++ Coarsening elements if error smaller " << maxErrToCoarse << ".\n");
250
251// reset counter
252 int numMarkedRefine = 0, numMarkedCoarse = 0;
253
254 const_iterator iter = dd->template begin<TElem>();
255 const_iterator iterEnd = dd->template end<TElem>();
256
257// loop elements for marking
258 for(; iter != iterEnd; ++iter)
259 {
260 // get element
261 TElem* elem = *iter;
262
263 // marks for refinement
264 if(aaError[elem] >= minErrToRefine)
265 if(dd->multi_grid()->get_level(elem) <= maxLevel)
266 {
267 refiner.mark(elem, RM_REFINE);
268 numMarkedRefine++;
269 }
270
271 // marks for coarsening
272 if(aaError[elem] <= maxErrToCoarse)
273 {
274 refiner.mark(elem, RM_COARSEN);
275 numMarkedCoarse++;
276 }
277 }
278
279#ifdef UG_PARALLEL
280 if(pcl::NumProcs() > 1)
281 {
282 UG_LOG(" +++ Marked for refinement on Proc "<<pcl::ProcRank()<<": " << numMarkedRefine << " Elements.\n");
283 UG_LOG(" +++ Marked for coarsening on Proc "<<pcl::ProcRank()<<": " << numMarkedCoarse << " Elements.\n");
285 int numMarkedRefineLocal = numMarkedRefine, numMarkedCoarseLocal = numMarkedCoarse;
286 numMarkedRefine = com.allreduce(numMarkedRefineLocal, PCL_RO_SUM);
287 numMarkedCoarse = com.allreduce(numMarkedCoarseLocal, PCL_RO_SUM);
288 }
289#endif
290
291 UG_LOG(" +++ Marked for refinement: " << numMarkedRefine << " Elements.\n");
292 UG_LOG(" +++ Marked for coarsening: " << numMarkedCoarse << " Elements.\n");
293}
294
296
311template<typename TElem>
314 IRefiner& refiner,
316 number tol,
317 int maxLevel
318)
319{
320 typedef typename DoFDistribution::traits<TElem>::const_iterator const_iterator;
321
322// compute minimal/maximal/ total error and number of elements
323 number min, max, totalErr;
324 size_t numElem;
325 ComputeMinMaxTotal(aaError, dd, min, max, totalErr, numElem);
326
327// check if total error is smaller than tolerance; if that is the case we're done
328 if (totalErr < tol)
329 {
330 UG_LOG(" +++ Total error "<<totalErr<<" smaller than TOL (" << tol << "). "
331 "No refinement necessary." << std::endl);
332 return;
333 }
334
335// compute minimum
336 //number minErrToRefine = max * refineFrac;
337 number minErrToRefine = tol / numElem;
338 UG_LOG(" +++ Refining elements if error greater " << tol << "/" << numElem <<
339 " = " << minErrToRefine << ".\n");
340
341// reset counter
342 size_t numMarkedRefine = 0;
343
344 const_iterator iter = dd->template begin<TElem>();
345 const_iterator iterEnd = dd->template end<TElem>();
346
347// loop elements for marking
348 for (; iter != iterEnd; ++iter)
349 {
350 // get element
351 TElem* elem = *iter;
352
353 // if no error value exists: ignore (might be newly added by refinement);
354 // newly added elements are supposed to have a negative error estimator
355 if (aaError[elem] < 0) continue;
356
357 // marks for refinement
358 if (aaError[elem] >= minErrToRefine)
359 if (dd->multi_grid()->get_level(elem) < maxLevel)
360 {
361 refiner.mark(elem, RM_REFINE);
362 numMarkedRefine++;
363 }
364 }
365
366#ifdef UG_PARALLEL
367 if (pcl::NumProcs() > 1)
368 {
370 size_t numMarkedRefineLocal = numMarkedRefine;
371 numMarkedRefine = com.allreduce(numMarkedRefineLocal, PCL_RO_SUM);
372 }
373#endif
374
375 UG_LOG(" +++ Marked for refinement: " << numMarkedRefine << " elements.\n");
376}
377
379
393template<typename TElem>
396 ug::Attachment<number> >& aaError,
397 IRefiner& refiner,
399 number TOL,
400 number safety,
401 int minLevel = 0
402)
403{
404 typedef typename DoFDistribution::traits<TElem>::const_iterator const_iterator;
405
406// compute minimal/maximal/ total error and number of elements
407 number min, max, totalErr;
408 size_t numElem;
409 ComputeMinMaxTotal(aaError, dd, min, max, totalErr, numElem);
410
411// compute maximum
412 //number maxErrToCoarse = min * (1+coarseFrac);
413 //if (maxErrToCoarse < TOL/numElem) maxErrToCoarse = TOL/numElem;
414 number maxErrToCoarse = TOL / numElem / 4.0 / safety; // = tol_per_elem / (2^2*safety_factor)
415 UG_LOG(" +++ Coarsening elements if avg child error smaller than "<< maxErrToCoarse << ".\n");
416
417// reset counter
418 size_t numMarkedCoarse = 0;
419
420 const_iterator iter = dd->template begin<TElem>();
421 const_iterator iterEnd = dd->template end<TElem>();
422
423// loop elements for marking
424 for (; iter != iterEnd; ++iter)
425 {
426 // get element
427 TElem* elem = *iter;
428
429 // ignore if already marked for coarsening
430 if (refiner.get_mark(elem) & RM_COARSEN) continue;
431
432 // ignore if level too low
433 if (dd->multi_grid()->get_level(elem) <= minLevel)
434 continue;
435
436 // get parent
437 TElem* parent = dynamic_cast<TElem*>(dd->multi_grid()->get_parent(elem));
438 if (!parent) continue; // either dynamic casting failed or parent does not exist
439
440 // sum up error values over all children of parent
441 number sum = 0.0;
442 size_t cnt = 0;
443 size_t nCh = dd->multi_grid()->num_children<TElem>(parent);
444 for (size_t ch = 0; ch < nCh; ch++)
445 {
446 TElem* child = dd->multi_grid()->get_child<TElem>(parent, ch);
447
448 // if no error value exists: ignore (might be newly added by refinement);
449 // newly added elements are supposed to have a negative error estimator
450 if (aaError[child] < 0) continue;
451
452 sum += aaError[child];
453 cnt++;
454 }
455
456 // marks for coarsening
457 if (sum/cnt <= maxErrToCoarse)
458 {
459 for (size_t ch = 0; ch < nCh; ch++)
460 {
461 TElem* child = dd->multi_grid()->get_child<TElem>(parent, ch);
462 if (refiner.get_mark(child) & RM_COARSEN) continue;
463 refiner.mark(child, RM_COARSEN);
464 numMarkedCoarse++;
465 }
466 }
467 }
468
469#ifdef UG_PARALLEL
470 if (pcl::NumProcs() > 1)
471 {
473 size_t numMarkedCoarseLocal = numMarkedCoarse;
474 numMarkedCoarse = com.allreduce(numMarkedCoarseLocal, PCL_RO_SUM);
475 }
476#endif
477
478 UG_LOG(" +++ Marked for coarsening: " << numMarkedCoarse << " Elements.\n");
479}
480
482
496template<typename TElem>
498 IRefiner& refiner,
500 number refTol,
501 number coarsenTol,
502 int minLevel,
503 int maxLevel,
504 bool refTopLvlOnly = false)
505{
506 typedef typename DoFDistribution::traits<TElem>::const_iterator const_iterator;
507
508 int numMarkedRefine = 0, numMarkedCoarse = 0;
509 const_iterator iter = dd->template begin<TElem>();
510 const_iterator iterEnd = dd->template end<TElem>();
511 const MultiGrid* mg = dd->multi_grid().get();
512 int topLvl = 0;
513 if(mg)
514 topLvl = (int)mg->top_level();
515 else
516 refTopLvlOnly = false;
517
518// loop elements for marking
519 for(; iter != iterEnd; ++iter)
520 {
521 TElem* elem = *iter;
522
523 // marks for refinement
524 if((refTol >= 0)
525 && (aaError[elem] > refTol)
526 && (dd->multi_grid()->get_level(elem) < maxLevel)
527 && ((!refTopLvlOnly) || (mg->get_level(elem) == topLvl)))
528 {
529 refiner.mark(elem, RM_REFINE);
530 numMarkedRefine++;
531 }
532
533 // marks for coarsening
534 if((coarsenTol >= 0)
535 && (aaError[elem] < coarsenTol)
536 && (dd->multi_grid()->get_level(elem) > minLevel))
537 {
538 refiner.mark(elem, RM_COARSEN);
539 numMarkedCoarse++;
540 }
541 }
542
543#ifdef UG_PARALLEL
544 if(pcl::NumProcs() > 1)
545 {
546 UG_LOG(" +++ Marked for refinement on Proc "<<pcl::ProcRank()<<": " << numMarkedRefine << " Elements.\n");
547 UG_LOG(" +++ Marked for coarsening on Proc "<<pcl::ProcRank()<<": " << numMarkedCoarse << " Elements.\n");
549 int numMarkedRefineLocal = numMarkedRefine, numMarkedCoarseLocal = numMarkedCoarse;
550 numMarkedRefine = com.allreduce(numMarkedRefineLocal, PCL_RO_SUM);
551 numMarkedCoarse = com.allreduce(numMarkedCoarseLocal, PCL_RO_SUM);
552 }
553#endif
554
555 UG_LOG(" +++ Marked for refinement: " << numMarkedRefine << " Elements.\n");
556 UG_LOG(" +++ Marked for coarsening: " << numMarkedCoarse << " Elements.\n");
557}
558
559}// end of namespace
560
561#endif
Definition smart_pointer.h:296
const T * get() const
Definition smart_pointer.h:409
Definition pcl_process_communicator.h:70
void allreduce(const void *sendBuf, void *recBuf, int count, DataType type, ReduceOperation op) const
performs MPI_Allreduce on the processes of the communicator.
Definition pcl_process_communicator.cpp:318
A generic specialization of IAttachment.
Definition attachment_pipe.h:263
the generic attachment-accessor for access to grids attachment pipes.
Definition grid.h:182
The refiner interface allows to mark elements for refinement and to call refine.
Definition refiner_interface.h:67
virtual bool mark(Vertex *v, RefinementMark refMark=RM_REFINE)
Marks an element for refinement. Default implementation is empty.
Definition refiner_interface.h:103
virtual RefinementMark get_mark(Vertex *v) const
Returns the mark of a given element. Default returns RM_REFINE.
Definition refiner_interface.h:198
Definition multi_grid.h:72
int get_level(TElem *elem) const
Definition multi_grid.h:206
size_t top_level() const
index of the highest level.
Definition multi_grid_impl.hpp:86
@ RM_COARSEN
the element is coarsened (only valid for adaptive multi-grid refinement)
Definition refiner_interface.h:56
@ RM_REFINE
DEPRECATED. Use RM_FULL instead.
Definition refiner_interface.h:55
#define PCL_RO_SUM
Definition pcl_methods.h:63
int ProcRank()
returns the rank of the process
Definition pcl_base.cpp:83
#define PCL_RO_MAX
Definition pcl_methods.h:61
#define PCL_RO_MIN
Definition pcl_methods.h:62
int NumProcs()
returns the number of processes
Definition pcl_base.cpp:91
#define UG_THROW(msg)
Definition error.h:57
#define UG_LOG(msg)
Definition log.h:367
double number
Definition types.h:124
the ug namespace
void ComputeMinMax(MultiGrid::AttachmentAccessor< TElem, ug::Attachment< number > > &aaError, ConstSmartPtr< DoFDistribution > dd, number &min, number &max, number &totalErr, size_t &numElem, number &minLocal, number &maxLocal, number &totalErrLocal, size_t &numElemLocal)
helper function that computes min/max and total of error indicators
Definition error_indicator_util.h:135
void MarkElements(MultiGrid::AttachmentAccessor< TElem, ug::Attachment< number > > &aaError, IRefiner &refiner, ConstSmartPtr< DoFDistribution > dd, number TOL, number refineFrac, number coarseFrac, int maxLevel)
marks elements according to an attached error value field
Definition error_indicator_util.h:222
void MarkElementsAbsolute(MultiGrid::AttachmentAccessor< TElem, ug::Attachment< number > > &aaError, IRefiner &refiner, ConstSmartPtr< DoFDistribution > dd, number refTol, number coarsenTol, int minLevel, int maxLevel, bool refTopLvlOnly=false)
marks elements according to an attached error value field
Definition error_indicator_util.h:497
void MarkElementsForCoarsening(MultiGrid::AttachmentAccessor< TElem, ug::Attachment< number > > &aaError, IRefiner &refiner, ConstSmartPtr< DoFDistribution > dd, number TOL, number safety, int minLevel=0)
marks elements for coarsening according to an attached error value field
Definition error_indicator_util.h:395
number ComputeAvg(MultiGrid::AttachmentAccessor< TElem, ug::Attachment< number > > &aaError2, ConstSmartPtr< DoFDistribution > dd, number &min, number &max, number &sum, number &errSq, size_t &numElem, number &minLocal, number &maxLocal, number &sumLocal, size_t &numElemLocal)
helper function that computes min/max and total of error indicators
Definition error_indicator_util.h:54
void ComputeMinMaxTotal(MultiGrid::AttachmentAccessor< TElem, ug::Attachment< number > > &aaError2, ConstSmartPtr< DoFDistribution > dd, number &min, number &max, number &totalErr, size_t &numElem)
helper function that computes min/max and total of error indicators
Definition error_indicator_util.h:196
void MarkElementsForRefinement(MultiGrid::AttachmentAccessor< TElem, ug::Attachment< number > > &aaError, IRefiner &refiner, ConstSmartPtr< DoFDistribution > dd, number tol, int maxLevel)
marks elements according for refinement to an attached error value field
Definition error_indicator_util.h:313
SurfaceView::traits< TElem >::const_iterator const_iterator
Definition dof_distribution.h:82