Line data Source code
1 : /*
2 : * Copyright (c) 2015: G-CSC, Goethe University Frankfurt
3 : * Author: Arne Nägel
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 : #ifndef __H__UG_DISC__ERROR_ELEM_MARKING_STRATEGY__
35 : #define __H__UG_DISC__ERROR_ELEM_MARKING_STRATEGY__
36 :
37 : #include "lib_grid/algorithms/debug_util.h" // for ElementDebugInfo
38 : #include "lib_grid/multi_grid.h"
39 : #include "lib_grid/refinement/refiner_interface.h"
40 : #include "lib_disc/dof_manager/dof_distribution.h"
41 : #include "lib_disc/function_spaces/grid_function.h"
42 : #include "error_indicator_util.h"
43 :
44 : namespace ug{
45 :
46 : /// Abstract base class for element marking (in adaptive refinement)
47 : template <typename TDomain> class IElementMarkingStrategy;
48 :
49 :
50 : /// This class encapsulates the multi-grid attachments for error estimation
51 : /** Purpose: replaces direct access to 'm_aaError' etc. */
52 : template <typename TDomain>
53 : class IMultigridElementIndicators
54 : {
55 : public:
56 : /// world dimension
57 : static const int dim = TDomain::dim;
58 :
59 : typedef typename domain_traits<dim>::element_type elem_type;
60 : typedef Attachment<number> error_attachment_type;
61 : typedef MultiGrid::AttachmentAccessor<elem_type, error_attachment_type > attachment_accessor_type;
62 :
63 : /// CTOR
64 0 : IMultigridElementIndicators() {}
65 :
66 : /// DTOR
67 : ~IMultigridElementIndicators()
68 : {
69 0 : detach_indicators();
70 0 : }
71 :
72 : /// Attach error indicator to multigrid
73 0 : void attach_indicators(SmartPtr<MultiGrid> pMG)
74 : {
75 : typedef typename domain_traits<dim>::element_type elem_type;
76 :
77 : // default value negative in order to distinguish between newly added elements (e.g. after refinement)
78 : // and elements which an error indicator is known for
79 0 : if (!pMG->has_attachment<elem_type>(m_aError))
80 0 : pMG->template attach_to_dv<elem_type>(m_aError, -1.0); // attach with default value
81 0 : m_pMG = pMG;
82 0 : m_aaError = attachment_accessor_type(*m_pMG, m_aError);
83 0 : }
84 :
85 : /// Detach error indicator from multigrid
86 0 : void detach_indicators()
87 : {
88 0 : if (m_pMG.invalid()) return; // no elements attached
89 0 : if (m_pMG->has_attachment<elem_type>(m_aError))
90 0 : m_pMG->template detach_from<elem_type>(m_aError);
91 : }
92 :
93 :
94 : /// returns error indicator value
95 : number& error(typename attachment_accessor_type::atraits::ConstElemPtr pElem)
96 : { return this->m_aaError[pElem]; }
97 :
98 : /// returns error indicator value
99 : const number& error(typename attachment_accessor_type::atraits::ConstElemPtr pElem) const
100 : { return this->m_aaError[pElem]; }
101 :
102 :
103 : friend class IElementMarkingStrategy<TDomain>;
104 :
105 :
106 : /// TODO: remove this function
107 : /// (mbreit: no, please leave it, it is very useful, at least with const access)
108 : attachment_accessor_type& errors()
109 0 : { return m_aaError; }
110 :
111 : protected:
112 : SmartPtr<MultiGrid> m_pMG; // multigrid
113 : error_attachment_type m_aError; // holding 'number' attachments
114 : attachment_accessor_type m_aaError;
115 : };
116 :
117 :
118 :
119 :
120 : /**
121 : *
122 : */
123 : template <typename TDomain>
124 : class IElementMarkingStrategy
125 : {
126 : public:
127 : /// world dimension
128 : static const int dim = TDomain::dim;
129 :
130 : /// element type to be marked
131 : typedef typename domain_traits<dim>::element_type elem_type;
132 : typedef typename Grid::AttachmentAccessor<elem_type, ug::Attachment<number> > elem_accessor_type;
133 :
134 0 : IElementMarkingStrategy()
135 0 : : m_latest_error(-1), m_latest_error_per_elem_max(-1), m_latest_error_per_elem_min(-1)
136 : {}
137 : virtual ~IElementMarkingStrategy(){};
138 :
139 :
140 : /// This function marks all elements
141 0 : void mark(IMultigridElementIndicators<TDomain>& mgElemIndicators,
142 : IRefiner& refiner,
143 : ConstSmartPtr<DoFDistribution> dd)
144 : {
145 : // forward call
146 0 : mark(mgElemIndicators.errors(), refiner, dd);
147 0 : };
148 :
149 : protected:
150 : /// DEPRECATED:
151 : virtual void mark(elem_accessor_type& aaError,
152 : IRefiner& refiner,
153 : ConstSmartPtr<DoFDistribution> dd) = 0;
154 :
155 : number m_latest_error;
156 : number m_latest_error_per_elem_max;
157 : number m_latest_error_per_elem_min;
158 :
159 : public:
160 0 : number global_estimated_error() const {return m_latest_error;}
161 0 : number global_estimated_error_per_elem_max() const {return m_latest_error_per_elem_max;}
162 0 : number global_estimated_error_per_elem_min() const {return m_latest_error_per_elem_min;}
163 : };
164 :
165 : /// M. Breit's standard refinement strategy
166 : template <typename TDomain>
167 : class StdRefinementMarkingStrategy : public IElementMarkingStrategy<TDomain>
168 : {
169 :
170 : public:
171 : typedef IElementMarkingStrategy<TDomain> base_type;
172 :
173 0 : StdRefinementMarkingStrategy(number tol, int max_level)
174 0 : : m_tol(tol), m_max_level(max_level) {};
175 :
176 0 : void set_tolerance(number tol) {m_tol = tol;}
177 0 : void set_max_level(int max_level) {m_max_level = max_level;}
178 :
179 : protected:
180 : void mark(typename base_type::elem_accessor_type& aaError,
181 : IRefiner& refiner,
182 : ConstSmartPtr<DoFDistribution> dd);
183 :
184 : number m_tol;
185 : int m_max_level;
186 : };
187 :
188 : template <typename TDomain>
189 0 : void StdRefinementMarkingStrategy<TDomain>::mark(typename base_type::elem_accessor_type& aaError,
190 : IRefiner& refiner,
191 : ConstSmartPtr<DoFDistribution> dd)
192 : {
193 0 : MarkElementsForRefinement<typename base_type::elem_type>(aaError, refiner, dd, m_tol, m_max_level);
194 0 : }
195 :
196 :
197 :
198 : /// mark everything if error too high and refinement allowed
199 : template <typename TDomain>
200 : class GlobalMarking : public IElementMarkingStrategy<TDomain>
201 : {
202 :
203 : public:
204 : typedef IElementMarkingStrategy<TDomain> base_type;
205 :
206 0 : GlobalMarking(number tol, size_t max_level)
207 0 : : m_tol(tol), m_max_level(max_level) {};
208 :
209 0 : void set_tolerance(number tol) {m_tol = tol;}
210 0 : void set_max_level(size_t max_level) {m_max_level = max_level;}
211 :
212 : protected:
213 : void mark(typename base_type::elem_accessor_type& aaError,
214 : IRefiner& refiner,
215 : ConstSmartPtr<DoFDistribution> dd);
216 :
217 : protected:
218 :
219 : number m_tol;
220 : size_t m_max_level;
221 : };
222 :
223 : template <typename TDomain>
224 0 : void GlobalMarking<TDomain>::mark(typename base_type::elem_accessor_type& aaError,
225 : IRefiner& refiner,
226 : ConstSmartPtr<DoFDistribution> dd)
227 : {
228 : number minElemErr;
229 : number maxElemErr;
230 : number errTotal;
231 : size_t numElem;
232 :
233 0 : ComputeMinMaxTotal(aaError, dd, minElemErr, maxElemErr, errTotal, numElem);
234 0 : if (errTotal <= m_tol || dd->multi_grid()->num_levels() > m_max_level)
235 0 : return;
236 :
237 : typedef typename DoFDistribution::traits<typename base_type::elem_type>::const_iterator const_iterator;
238 :
239 0 : const_iterator iter = dd->template begin<typename base_type::elem_type>();
240 0 : const_iterator iterEnd = dd->template end<typename base_type::elem_type>();
241 :
242 : // loop elements for marking
243 0 : for (; iter != iterEnd; ++iter)
244 0 : refiner.mark(*iter, RM_REFINE);
245 : }
246 :
247 :
248 :
249 : /// M. Breit's standard coarsening strategy
250 : template <typename TDomain>
251 : class StdCoarseningMarkingStrategy : public IElementMarkingStrategy<TDomain>
252 : {
253 :
254 : public:
255 : typedef IElementMarkingStrategy<TDomain> base_type;
256 :
257 0 : StdCoarseningMarkingStrategy(number tol)
258 0 : : m_tol(tol), m_safety(8.0), m_minLvl(0) {}
259 :
260 0 : StdCoarseningMarkingStrategy(number tol, number safety)
261 0 : : m_tol(tol), m_safety(safety), m_minLvl(0) {}
262 :
263 0 : StdCoarseningMarkingStrategy(number tol, number safety, int minLvl)
264 0 : : m_tol(tol), m_safety(safety), m_minLvl(minLvl) {}
265 :
266 0 : void set_tolerance(number tol) {m_tol = tol;}
267 0 : void set_safety_factor(number safety) {m_safety = safety;}
268 :
269 : protected:
270 : void mark(typename base_type::elem_accessor_type& aaError,
271 : IRefiner& refiner,
272 : ConstSmartPtr<DoFDistribution> dd);
273 :
274 : protected:
275 : number m_tol;
276 : number m_safety;
277 : int m_minLvl;
278 : };
279 :
280 : template <typename TDomain>
281 0 : void StdCoarseningMarkingStrategy<TDomain>::mark(typename base_type::elem_accessor_type& aaError,
282 : IRefiner& refiner,
283 : ConstSmartPtr<DoFDistribution> dd)
284 : {
285 0 : MarkElementsForCoarsening<typename base_type::elem_type>(aaError, refiner, dd, m_tol, m_safety, m_minLvl);
286 0 : }
287 :
288 :
289 : /* Generate an ordered list of \eta_k^2 (in descending order, ie. largest first) */
290 : template<class TElem>
291 0 : number CreateListOfElemWeights(
292 : Grid::AttachmentAccessor<TElem, ug::Attachment<number> > &aaError,
293 : typename DoFDistribution::traits<TElem>::const_iterator iterBegin,
294 : const typename DoFDistribution::traits<TElem>::const_iterator iterEnd,
295 : std::vector<double> &etaSq)
296 : {
297 : number localErrSq=0;
298 : typename DoFDistribution::traits<TElem>::const_iterator iter;
299 : size_t i=0;
300 0 : for (iter = iterBegin; iter != iterEnd; ++iter)
301 : {
302 0 : const double elemErrSq = aaError[*iter];
303 :
304 : // If no error value exists: ignore (might be newly added by refinement);
305 : // newly added elements are supposed to have a negative error estimator
306 0 : if (elemErrSq < 0) continue;
307 :
308 0 : etaSq[i++] = elemErrSq;
309 0 : localErrSq += elemErrSq;
310 : }
311 :
312 : // Sort in descending order (using default comparison).
313 0 : std::sort (etaSq.begin(), etaSq.end(), std::greater<double>());
314 0 : return localErrSq;
315 : };
316 :
317 :
318 :
319 : /* Generate an ordered list of \eta_k^2 (in descending order, ie. largest first) */
320 : template<class TElem>
321 0 : number CreateSortedListOfElems(
322 : Grid::AttachmentAccessor<TElem, ug::Attachment<number> > &aaError,
323 : typename DoFDistribution::traits<TElem>::const_iterator iterBegin,
324 : const typename DoFDistribution::traits<TElem>::const_iterator iterEnd,
325 : std::vector< std::pair<double, TElem*> > &etaSqList)
326 : {
327 :
328 : //typedef typename std::pair<double, TElem*> TPair;
329 : typename DoFDistribution::traits<TElem>::const_iterator iter;
330 :
331 : number localErrSq=0;
332 : size_t i=0;
333 0 : for (iter = iterBegin; iter != iterEnd; ++iter)
334 : {
335 0 : const double elemErrSq = aaError[*iter];
336 :
337 : // If no error value exists: ignore (might be newly added by refinement);
338 : // newly added elements are supposed to have a negative error estimator
339 0 : if (elemErrSq < 0) continue;
340 :
341 0 : etaSqList[i++] = std::make_pair<>(elemErrSq, *iter);
342 0 : localErrSq += elemErrSq;
343 : }
344 :
345 : // Sort in descending order (using default comparison).
346 0 : std::sort (etaSqList.begin(), etaSqList.end());
347 0 : return localErrSq;
348 : };
349 :
350 :
351 :
352 : /// Refine as many largest-error elements as necessary to reach tolerance.
353 : /// The expected tolerance is calculated using a user-given expected error reduction factor.
354 : template <typename TDomain>
355 : class ExpectedErrorMarkingStrategy : public IElementMarkingStrategy<TDomain>
356 : {
357 :
358 : public:
359 : typedef IElementMarkingStrategy<TDomain> base_type;
360 :
361 0 : ExpectedErrorMarkingStrategy(number tol, int max_level, number safetyFactor, number expectedReductionFactor)
362 0 : : m_tol(tol), m_max_level(max_level), m_safety(safetyFactor), m_expRedFac(expectedReductionFactor) {};
363 :
364 0 : void set_tolerance(number tol) {m_tol = tol;}
365 0 : void set_max_level(int max_level) {m_max_level = max_level;}
366 0 : void set_safety_factor(number safetyFactor) {m_safety = safetyFactor;}
367 0 : void set_expected_reduction_factor(number expectedReductionFactor) {m_expRedFac = expectedReductionFactor;}
368 :
369 : void mark(typename base_type::elem_accessor_type& aaError,
370 : IRefiner& refiner,
371 : ConstSmartPtr<DoFDistribution> dd);
372 :
373 : protected:
374 : void merge_sorted_lists
375 : (
376 : std::vector<std::pair<number, int> >::iterator beginFirst,
377 : std::vector<std::pair<number, int> >::iterator beginSecond,
378 : std::vector<std::pair<number, int> >::iterator end
379 : );
380 :
381 : protected:
382 : number m_tol;
383 : int m_max_level;
384 : number m_safety;
385 : number m_expRedFac;
386 : };
387 :
388 :
389 : template <typename TDomain>
390 : struct ElemErrorSortDesc
391 : {
392 : typedef typename domain_traits<TDomain::dim>::element_type elem_type;
393 : typedef typename Grid::AttachmentAccessor<elem_type, ug::Attachment<number> > error_accessor_type;
394 : ElemErrorSortDesc(const error_accessor_type& aaErr)
395 : : m_aaErr(aaErr) {}
396 :
397 : bool operator()(const elem_type* elem1, const elem_type* elem2)
398 : {
399 0 : return m_aaErr[elem1] > m_aaErr[elem2];
400 : }
401 :
402 : const error_accessor_type& m_aaErr;
403 : };
404 :
405 :
406 : template <typename TDomain>
407 : void ExpectedErrorMarkingStrategy<TDomain>::merge_sorted_lists
408 : (
409 : std::vector<std::pair<number, int> >::iterator beginFirst,
410 : std::vector<std::pair<number, int> >::iterator beginSecond,
411 : std::vector<std::pair<number, int> >::iterator end
412 : )
413 : {
414 : const size_t nVal = std::distance(beginFirst, end);
415 : std::vector<std::pair<number, int> > sorted;
416 : sorted.reserve(nVal);
417 :
418 : std::vector<std::pair<number, int> >::iterator it1 = beginFirst;
419 : std::vector<std::pair<number, int> >::iterator it2 = beginSecond;
420 : while (it1 != beginSecond && it2 != end)
421 : {
422 : if (it1->first > it2->first)
423 : {
424 : sorted.push_back(*it1);
425 : ++it1;
426 : }
427 : else
428 : {
429 : sorted.push_back(*it2);
430 : ++it2;
431 : }
432 : }
433 :
434 : for (; it1 != beginSecond; ++it1)
435 : sorted.push_back(*it1);
436 : for (; it2 != end; ++it2)
437 : sorted.push_back(*it2);
438 :
439 : size_t i = 0;
440 : for (it1 = beginFirst; it1 != end; ++it1, ++i)
441 : *it1 = sorted[i];
442 : }
443 :
444 :
445 :
446 : template <typename TDomain>
447 0 : void ExpectedErrorMarkingStrategy<TDomain>::mark
448 : (
449 : typename base_type::elem_accessor_type& aaErrorSq,
450 : IRefiner& refiner,
451 : ConstSmartPtr<DoFDistribution> dd
452 : )
453 : {
454 : typedef typename base_type::elem_type TElem;
455 : typedef typename DoFDistribution::traits<TElem>::const_iterator const_iterator;
456 :
457 : // create vector of local element (squared) errors
458 0 : const_iterator iter = dd->template begin<TElem>();
459 0 : const const_iterator iterEnd = dd->template end<TElem>();
460 : std::vector<TElem*> elemVec;
461 : number locError = 0.0;
462 0 : for (; iter != iterEnd; ++iter)
463 : {
464 0 : TElem* elem = *iter;
465 :
466 0 : if (aaErrorSq[elem] < 0)
467 0 : continue;
468 :
469 0 : if (aaErrorSq[elem] != aaErrorSq[elem])
470 : {
471 0 : UG_LOG_ALL_PROCS("NaN error indicator on " << ElementDebugInfo(*refiner.grid(), elem));
472 0 : continue;
473 0 : }
474 :
475 0 : locError += aaErrorSq[elem];
476 :
477 0 : if (dd->multi_grid()->get_level(elem) < m_max_level)
478 0 : elemVec.push_back(elem);
479 : }
480 : const size_t nLocalElem = elemVec.size();
481 :
482 : // sort vector of elements locally (descending)
483 : ElemErrorSortDesc<TDomain> eeSort(aaErrorSq);
484 0 : std::sort(elemVec.begin(), elemVec.end(), eeSort);
485 :
486 : // create sorted vector of errors
487 0 : std::vector<number> etaSq(nLocalElem);
488 0 : for (size_t i = 0; i < nLocalElem; ++i)
489 0 : etaSq[i] = aaErrorSq[elemVec[i]];
490 :
491 : #ifdef UG_PARALLEL
492 : const int nProcs = pcl::NumProcs();
493 : if (nProcs > 1)
494 : {
495 : // calculate global error
496 : pcl::ProcessCommunicator pc;
497 : const int rootProc = pcl::NumProcs() - 1;
498 : const number globError = pc.allreduce(locError, PCL_RO_SUM);
499 : UG_LOGN(" +++ Element errors: sumEtaSq = " << globError << ".");
500 :
501 : // construct global sorted vector of elem errors
502 : // gather on root proc
503 : // choose root as the highest-rank proc (proc 0 has least free memory)
504 : std::vector<number> rcvBuf;
505 : std::vector<int> vSizes;
506 : pc.gatherv(rcvBuf, etaSq, rootProc, &vSizes);
507 :
508 : // merge of pre-sorted local vectors
509 : std::vector<int> nRefineElems;
510 : size_t globNumRefineElems = 0;
511 : if (pcl::ProcRank() == rootProc)
512 : {
513 : // associate each error value with the proc it belongs to
514 : // and calculate offsets at the same time
515 : std::vector<size_t> vOffsets(nProcs, 0);
516 : const size_t nGlobElem = rcvBuf.size();
517 : std::vector<std::pair<number, int> > vGlobErrors(nGlobElem);
518 : size_t e = 0;
519 : for (int p = 0; p < nProcs; ++p)
520 : {
521 : const size_t szp = vSizes[p];
522 : for (size_t i = 0; i < szp; ++i, ++e)
523 : {
524 : vGlobErrors[e].first = rcvBuf[e];
525 : vGlobErrors[e].second = p;
526 : }
527 : if (p < nProcs-1)
528 : vOffsets[p+1] = e;
529 : }
530 :
531 : // free rcv buffer
532 : {
533 : std::vector<number> tmp;
534 : tmp.swap(rcvBuf);
535 : }
536 :
537 : // merge pre-sorted proc error lists ((log p) * N)
538 : size_t nLists = 2;
539 : while (true)
540 : {
541 : const size_t nMerge = nProcs / std::min(nLists, (size_t) nProcs);
542 : for (size_t m = 0; m < nMerge; ++m)
543 : {
544 : std::vector<std::pair<number, int> >::iterator beginFirst =
545 : vGlobErrors.begin() + vOffsets[m*nLists];
546 : std::vector<std::pair<number, int> >::iterator beginSecond =
547 : vGlobErrors.begin() + vOffsets[m*nLists + nLists/2];
548 : std::vector<std::pair<number, int> >::iterator endSecond =
549 : (m+1)*nLists < (size_t) nProcs ? vGlobErrors.begin() + vOffsets[(m+1)*nLists] : vGlobErrors.end();
550 :
551 : merge_sorted_lists(beginFirst, beginSecond, endSecond);
552 : }
553 :
554 : if (nLists >= (size_t) nProcs)
555 : break;
556 :
557 : nLists = nLists << 1;
558 : }
559 :
560 : // decide how many elements on each proc have to be refined
561 : const number requiredReduction = globError > m_tol ? globError - m_safety*m_tol : 0.0;
562 : number red = 0.0;
563 : nRefineElems.resize(nProcs, 0);
564 : for (; red < requiredReduction && globNumRefineElems < nGlobElem; ++globNumRefineElems)
565 : {
566 : red += (1.0-m_expRedFac) * vGlobErrors[globNumRefineElems].first;
567 : ++nRefineElems[vGlobErrors[globNumRefineElems].second];
568 : }
569 : }
570 :
571 : // tell each proc how many of their elements are to be marked
572 : int nRefElems = 0;
573 : pc.scatter(GetDataPtr(nRefineElems), 1, PCL_DT_INT, &nRefElems, 1, PCL_DT_INT, rootProc);
574 :
575 : pc.broadcast(globNumRefineElems, rootProc);
576 :
577 : // mark for refinement
578 : UG_COND_THROW((size_t) nRefElems > nLocalElem, "More elements supposedly need refinement here ("
579 : << nRefElems << ") than are refineable (" << nLocalElem << ").");
580 : for (size_t i = 0; i < (size_t) nRefElems; ++i)
581 : refiner.mark(elemVec[i], RM_REFINE);
582 :
583 : if (globNumRefineElems)
584 : {
585 : UG_LOGN(" +++ Marked for refinement: " << globNumRefineElems << " elements");
586 : }
587 : else
588 : {
589 : UG_LOGN(" +++ No refinement necessary.");
590 : }
591 : }
592 : else
593 : {
594 : #endif
595 :
596 : // plan refinement of elements until expected error reduction is enough
597 : // to get below tolerance (or all elements are marked)
598 0 : const number requiredReduction = locError - m_safety*m_tol;
599 : UG_LOGN(" +++ Element errors: sumEtaSq = " << locError << ".");
600 : number red = 0.0;
601 : size_t i = 0;
602 0 : for (; red < requiredReduction && i < nLocalElem; ++i)
603 : {
604 0 : red += m_expRedFac * etaSq[i];
605 0 : refiner.mark(elemVec[i], RM_REFINE);
606 : }
607 0 : if (i)
608 : {
609 : UG_LOGN(" +++ Marked for refinement: " << i << " elements.");
610 : }
611 : else
612 : {
613 : UG_LOGN(" +++ No refinement necessary.");
614 : }
615 :
616 : #ifdef UG_PARALLEL
617 : }
618 : #endif
619 0 : }
620 :
621 :
622 :
623 : /// Marks elements with \eta_K >= \theta \max_{K'} \eta_{K'} for refinement
624 : //! (cf. Verfuerth script)
625 : template <typename TDomain>
626 : class MaximumMarking : public IElementMarkingStrategy<TDomain>{
627 :
628 : public:
629 : typedef IElementMarkingStrategy<TDomain> base_type;
630 0 : MaximumMarking(number theta=1.0)
631 0 : : m_theta(theta), m_theta_min(0.0), m_eps(0.01), m_max_level(100), m_min_level(0) {};
632 0 : MaximumMarking(number theta, number eps)
633 0 : : m_theta(theta), m_theta_min(0.0), m_eps (eps), m_max_level(100), m_min_level(0) {};
634 0 : MaximumMarking(number theta_max, number theta_min, number eps)
635 0 : : m_theta(theta_max), m_theta_min(theta_min), m_eps (eps), m_max_level(100), m_min_level(0) {};
636 :
637 : protected:
638 : void mark(typename base_type::elem_accessor_type& aaErrorSq, IRefiner& refiner, ConstSmartPtr<DoFDistribution> dd);
639 :
640 : public:
641 0 : void set_max_level(int lvl) {m_max_level = lvl;}
642 0 : void set_min_level(int lvl) {m_min_level = lvl;}
643 :
644 : protected:
645 :
646 : number m_theta, m_theta_min;
647 : number m_eps;
648 : int m_max_level, m_min_level;
649 : };
650 :
651 :
652 :
653 : template <typename TDomain>
654 0 : void MaximumMarking<TDomain>::mark(typename base_type::elem_accessor_type& aaErrorSq,
655 : IRefiner& refiner,
656 : ConstSmartPtr<DoFDistribution> dd)
657 : {
658 : typedef typename base_type::elem_type TElem;
659 : typedef typename DoFDistribution::traits<TElem>::const_iterator const_iterator;
660 :
661 : // compute minimal/maximal/ total error and number of elements
662 :
663 : number minElemErr, minElemErrLocal;
664 : number maxElemErr, maxElemErrLocal;
665 : number errTotal, errLocal;
666 : size_t numElem, numElemLocal;
667 :
668 0 : ComputeMinMax(aaErrorSq, dd, minElemErr, maxElemErr, errTotal, numElem,
669 : minElemErrLocal, maxElemErrLocal, errLocal, numElemLocal);
670 :
671 0 : this->m_latest_error = sqrt(errTotal);
672 0 : this->m_latest_error_per_elem_max = maxElemErr;
673 0 : this->m_latest_error_per_elem_min = minElemErr;
674 :
675 : // init iterators
676 : const_iterator iter;
677 0 : const const_iterator iterEnd = dd->template end<TElem>();
678 :
679 : // determine (local) number of excess elements
680 0 : const size_t ndiscard = (size_t) (numElemLocal*m_eps); // TODO: on every process?
681 : UG_LOG(" +++ MaximumMarking: Found max "<< maxElemErr << ", ndiscard="<<ndiscard<<".\n");
682 :
683 0 : if (numElemLocal > 0)
684 : {
685 : // create sorted array of elem weights
686 : std::vector<double> etaSq;
687 0 : etaSq.resize(numElemLocal);
688 0 : CreateListOfElemWeights<TElem>(aaErrorSq,dd->template begin<TElem>(), iterEnd, etaSq);
689 : UG_ASSERT(numElemLocal==etaSq.size(), "Huhh: number of elements does not match!");
690 : UG_ASSERT(numElemLocal > ndiscard, "Huhh: number of elements does not match!");
691 0 : maxElemErr = etaSq[ndiscard];
692 0 : }
693 :
694 : // compute parallel threshold
695 : #ifdef UG_PARALLEL
696 : if (pcl::NumProcs() > 1)
697 : {
698 : pcl::ProcessCommunicator com;
699 : number maxElemErrLocal = maxElemErr;
700 : maxElemErr = com.allreduce(maxElemErrLocal, PCL_RO_MAX);
701 : UG_LOG(" +++ Maximum for refinement: " << maxElemErr << ".\n");
702 : }
703 : #else
704 0 : UG_LOG(" +++ Skipping " << ndiscard << " elements; new (parallel) max." << maxElemErr << ".\n");
705 : #endif
706 :
707 : // refine all element above threshold
708 0 : const number top_threshold = maxElemErr*m_theta;
709 0 : const number bot_threshold = maxElemErr*m_theta_min;
710 0 : UG_LOG(" +++ Refining elements, if error > " << maxElemErr << "*" << m_theta <<
711 : " = " << top_threshold << ".\n");
712 0 : UG_LOG(" +++ Coarsening elements, if error < " << maxElemErr << "*" << m_theta_min <<
713 : " = " << bot_threshold << ".\n");
714 :
715 : // reset counter
716 : std::size_t numMarkedRefine = 0;
717 : std::size_t numMarkedCoarse = 0;
718 :
719 : // loop elements for marking
720 0 : for (iter = dd->template begin<TElem>(); iter != iterEnd; ++iter)
721 : {
722 : // get element
723 : TElem* elem = *iter;
724 :
725 : // if no error value exists: ignore (might be newly added by refinement);
726 : // newly added elements are supposed to have a negative error estimator
727 0 : if (aaErrorSq[elem] < 0) continue;
728 :
729 : // mark for refinement
730 0 : if ((aaErrorSq[elem] > top_threshold) && (dd->multi_grid()->get_level(elem) <= m_max_level))
731 : {
732 0 : refiner.mark(elem, RM_REFINE);
733 0 : numMarkedRefine++;
734 : }
735 :
736 : // mark for coarsening
737 0 : if ((aaErrorSq[elem] < bot_threshold) && (dd->multi_grid()->get_level(elem) >= m_min_level))
738 : {
739 0 : refiner.mark(elem, RM_COARSEN);
740 0 : numMarkedCoarse++;
741 : }
742 : }
743 :
744 : #ifdef UG_PARALLEL
745 : if (pcl::NumProcs() > 1)
746 : {
747 : pcl::ProcessCommunicator com;
748 : std::size_t numMarkedRefineLocal = numMarkedRefine;
749 : numMarkedRefine = com.allreduce(numMarkedRefineLocal, PCL_RO_SUM);
750 : UG_LOG(" +++ MaximumMarking: Marked for refinement: " << numMarkedRefine << " ("<< numMarkedRefineLocal << ") elements.\n");
751 : }
752 : else
753 : #endif
754 : { UG_LOG(" +++ MaximumMarking: refinement: " << numMarkedRefine << ", coarsening" << numMarkedCoarse << " elements.\n") }
755 :
756 :
757 0 : }
758 :
759 :
760 :
761 : /// Marks element with smallest \eta_i^2
762 : //! (cf. Verfuerth script)
763 : template <typename TDomain>
764 : class APosterioriCoarsening : public IElementMarkingStrategy<TDomain>{
765 :
766 : public:
767 : typedef IElementMarkingStrategy<TDomain> base_type;
768 0 : APosterioriCoarsening(number theta=0.1)
769 0 : : m_theta(theta), m_max_level(100), m_min_level(0) {};
770 : protected:
771 : void mark(typename base_type::elem_accessor_type& aaErrorSq, IRefiner& refiner, ConstSmartPtr<DoFDistribution> dd);
772 : public:
773 : void set_max_level(int lvl) {m_max_level = lvl;}
774 : void set_min_level(int lvl) {m_min_level = lvl;}
775 :
776 : protected:
777 :
778 : number m_theta;
779 : int m_max_level, m_min_level;
780 : };
781 :
782 :
783 :
784 : template <typename TDomain>
785 0 : void APosterioriCoarsening<TDomain>::mark(typename base_type::elem_accessor_type& aaErrorSq,
786 : IRefiner& refiner,
787 : ConstSmartPtr<DoFDistribution> dd)
788 : {
789 : typedef typename base_type::elem_type TElem;
790 : //typedef typename DoFDistribution::traits<TElem>::const_iterator const_iterator;
791 : typedef typename std::pair<double, TElem*> TPair;
792 : typedef typename std::vector<TPair> TPairVector;
793 :
794 : // Compute minimal/maximal/ total error and number of elements.
795 : number minElemErrSq, minElemErrSqLocal;
796 : number maxElemErrSq, maxElemErrSqLocal;
797 : number errSqTotal, errSqLocal;
798 : size_t numElem, numElemLocal;
799 :
800 : size_t numCoarsened=0;
801 :
802 :
803 0 : ComputeMinMax(aaErrorSq, dd, minElemErrSq, maxElemErrSq, errSqTotal, numElem,
804 : minElemErrSqLocal, maxElemErrSqLocal, errSqLocal, numElemLocal);
805 :
806 0 : this->m_latest_error = sqrt(errSqTotal);
807 0 : this->m_latest_error_per_elem_max = maxElemErrSq;
808 0 : this->m_latest_error_per_elem_min = minElemErrSq;
809 :
810 :
811 :
812 0 : if (numElemLocal > 0)
813 : {
814 : // Create sorted array of elem weights.
815 : TPairVector etaSqVec;
816 0 : etaSqVec.resize(numElemLocal);
817 0 : CreateSortedListOfElems<TElem>(aaErrorSq,dd->template begin<TElem>(), dd->template end<TElem>(), etaSqVec);
818 :
819 :
820 : UG_ASSERT(numElemLocal==etaSqVec.size(), "Huhh: number of elements does not match!");
821 :
822 : {
823 :
824 0 : const double mySumSq = errSqLocal*m_theta;
825 : double localSumSq = 0.0;
826 :
827 : typename TPairVector::const_iterator iterEnd = etaSqVec.end();
828 :
829 0 : for (typename TPairVector::iterator iter = etaSqVec.begin();
830 0 : (iter != iterEnd) && (localSumSq < mySumSq); ++iter)
831 : {
832 :
833 0 : localSumSq += iter->first;
834 0 : if (localSumSq < mySumSq) {
835 0 : refiner.mark(iter->second, RM_COARSEN);
836 0 : numCoarsened++;
837 : }
838 :
839 : }
840 : UG_LOG(" +++ APosterioriCoarsening: localSumSq = "<< localSumSq << " >" << mySumSq << std::endl);
841 :
842 : }
843 :
844 :
845 :
846 0 : }
847 :
848 :
849 : #ifdef UG_PARALLEL
850 : if (pcl::NumProcs() > 1)
851 : {
852 : pcl::ProcessCommunicator com;
853 : std::size_t numCoarsenedLocal = numCoarsened;
854 : numCoarsened = com.allreduce(numCoarsenedLocal, PCL_RO_SUM);
855 : UG_LOG(" +++ APosterioriCoarsening: Marked for coarsening: " << numCoarsened << " ("<< numCoarsenedLocal << ") elements.\n");
856 : }
857 : else
858 : #endif
859 : { UG_LOG(" +++ APosterioriCoarsening: coarsening" << numCoarsened << " elements.\n") }
860 :
861 :
862 0 : }
863 :
864 :
865 : //! marks elements above a certain fraction of the maximum
866 : /*! Cf. Verfuerth' scriptum */
867 : template <typename TDomain>
868 : class EquilibrationMarkingStrategy : public IElementMarkingStrategy<TDomain>{
869 :
870 : public:
871 : typedef IElementMarkingStrategy<TDomain> base_type;
872 0 : EquilibrationMarkingStrategy(number theta_top=0.9)
873 0 : : m_theta_top(theta_top), m_theta_bot(0.0) {} //, m_max_level(100) {};
874 0 : EquilibrationMarkingStrategy(number theta_top, number theta_bot)
875 0 : : m_theta_top(theta_top), m_theta_bot(theta_bot) {} ;// , m_max_level(100) {};
876 :
877 :
878 : protected:
879 : void mark(typename base_type::elem_accessor_type& aaErrorSq,
880 : IRefiner& refiner,
881 : ConstSmartPtr<DoFDistribution> dd);
882 :
883 : number m_theta_top; // refine at least a certain fraction, 0.0 <= m_theta_top <= 1.0
884 : number m_theta_bot; // refine at least a certain fraction, 0.0 <= m_theta_bot <= 1.0
885 : // int m_max_level;
886 : };
887 :
888 :
889 :
890 :
891 :
892 : template <typename TDomain>
893 0 : void EquilibrationMarkingStrategy<TDomain>::mark(typename base_type::elem_accessor_type& aaErrorSq,
894 : IRefiner& refiner,
895 : ConstSmartPtr<DoFDistribution> dd)
896 : {
897 : typedef typename base_type::elem_type TElem;
898 : typedef typename DoFDistribution::traits<TElem>::const_iterator const_iterator;
899 :
900 : // compute minimal/maximal/total error and number of elements
901 : number minElemErr, minElemErrLocal;
902 : number maxElemErr, maxElemErrLocal;
903 0 : number errTotal= 0.0, errLocal= 0.0;
904 0 : size_t numElem= 0, numElemLocal;
905 :
906 : // compute weights
907 0 : ComputeMinMax(aaErrorSq, dd, minElemErr, maxElemErr, errTotal, numElem,
908 : minElemErrLocal, maxElemErrLocal, errLocal, numElemLocal);
909 :
910 : // init iterators
911 0 : const const_iterator iterEnd = dd->template end<TElem>();
912 : const_iterator iter;
913 :
914 : // create and fill sorted array of $\etaSq^2_i$ for all (local) elements
915 : std::vector<double> etaSq;
916 0 : etaSq.resize(numElemLocal);
917 0 : CreateListOfElemWeights<TElem>(aaErrorSq,dd->template begin<TElem>(), iterEnd, etaSq);
918 : UG_ASSERT(numElemLocal==etaSq.size(), "Huhh: number of elements does not match!");
919 :
920 : // compute thresholds
921 : UG_ASSERT( ((m_theta_top>=0.0) && (m_theta_top<=1.0)), "Huhh: m_theta_top invalid!");
922 : UG_ASSERT( ((m_theta_bot>=0.0) && (m_theta_bot<=1.0)), "Huhh: m_theta_top invalid!");
923 : UG_ASSERT( (m_theta_top>m_theta_bot), "Huhh: m_theta_top invalid!");
924 :
925 : // discard a fraction of elements
926 : // a) largest elements
927 : typename std::vector<double>::const_iterator top = etaSq.begin();
928 0 : for (number sumSq=0.0;
929 0 : (sumSq<m_theta_top*errLocal) && (top !=(etaSq.end()-1));
930 0 : ++top) { sumSq += *top; }
931 0 : number top_threshold = (top != etaSq.begin()) ? (*top + *(top-1))/2.0 : *top;
932 :
933 : // a) smallest elements
934 : typename std::vector<double>::const_iterator bot = etaSq.end()-1;
935 0 : for (number sumSq=0.0;
936 0 : (sumSq<m_theta_bot*errLocal) && (bot !=etaSq.begin() );
937 0 : --bot) { sumSq += *bot; }
938 0 : number bot_threshold = (bot != (etaSq.end()-1)) ? (*bot + *(bot+1))/2.0 : 0.0;
939 :
940 : UG_LOG(" +++ error = "<< errLocal << std::endl);
941 : UG_LOG(" +++ top_threshold= "<< top_threshold <<"( "<< top - etaSq.begin() << " cells)" << std::endl);
942 : UG_LOG(" +++ bot_threshold= "<< bot_threshold <<"( "<< etaSq.end()-bot << " cells)" << std::endl);
943 :
944 : // mark elements with maximal contribution
945 : size_t numMarkedRefine = 0;
946 : size_t numMarkedCoarse = 0;
947 0 : for (iter = dd->template begin<TElem>(); iter != iterEnd; ++iter)
948 : {
949 :
950 0 : const double elemErr = aaErrorSq[*iter]; // get element
951 0 : if (elemErr < 0) continue; // skip invalid
952 :
953 0 : if (elemErr > top_threshold)
954 : {
955 0 : refiner.mark(*iter, RM_REFINE);
956 0 : numMarkedRefine++;
957 : }
958 :
959 0 : if (elemErr < bot_threshold)
960 : {
961 0 : refiner.mark(*iter, RM_COARSEN);
962 0 : numMarkedCoarse++;
963 : }
964 : }
965 :
966 :
967 :
968 : UG_LOG(" +++ EquilibrationMarkingStrategy: Marked for refinement: "<< numMarkedRefine << ", " << numMarkedCoarse << " elements.\n");
969 :
970 0 : }
971 :
972 :
973 : /// Marks elements above \f$ \theta * (\mu + width * \sigma) \f$ for refinement
974 : //! where \f$ \mu = E[\eta^2], \sigma^2 = Var[\eta^2] \f$
975 : template <typename TDomain>
976 : class VarianceMarking : public IElementMarkingStrategy<TDomain>{
977 :
978 : public:
979 : typedef IElementMarkingStrategy<TDomain> base_type;
980 0 : VarianceMarking(number theta) : m_theta(theta), m_width(3.0), m_max_level(100) {};
981 0 : VarianceMarking(number theta, number width) : m_theta(theta), m_width (width), m_max_level(100) {};
982 :
983 :
984 : protected:
985 : void mark(typename base_type::elem_accessor_type& aaError2,
986 : IRefiner& refiner,
987 : ConstSmartPtr<DoFDistribution> dd);
988 :
989 : number m_theta;
990 : number m_width;
991 : int m_max_level;
992 : };
993 :
994 : template <typename TDomain>
995 0 : void VarianceMarking<TDomain>::mark(typename base_type::elem_accessor_type& aaError2,
996 : IRefiner& refiner,
997 : ConstSmartPtr<DoFDistribution> dd)
998 : {
999 : typedef typename base_type::elem_type TElem;
1000 : typedef typename DoFDistribution::traits<TElem>::const_iterator const_iterator;
1001 :
1002 : // compute minimal/maximal/ total error and number of elements
1003 :
1004 : number minElemErr, minElemErrLocal;
1005 : number maxElemErr, maxElemErrLocal;
1006 : number errTotal, errLocal;
1007 : size_t numElem, numElemLocal;
1008 :
1009 0 : ComputeMinMax(aaError2, dd, minElemErr, maxElemErr, errTotal, numElem,
1010 : minElemErrLocal, maxElemErrLocal, errLocal, numElemLocal);
1011 :
1012 0 : this->m_latest_error = sqrt(errTotal);
1013 0 : this->m_latest_error_per_elem_max = maxElemErr;
1014 0 : this->m_latest_error_per_elem_min = minElemErr;
1015 :
1016 : // number elemMean = sqrt(errTotal) / numElem;
1017 0 : number elemMean = errTotal / numElem;
1018 :
1019 : UG_LOG(" +++ VarianceMarking: Mean error : " << elemMean << " on "<< numElem << " elements.\n");
1020 :
1021 : // init iterators
1022 : const_iterator iter;
1023 0 : const const_iterator iterEnd = dd->template end<TElem>();
1024 :
1025 : number elemVar = 0.0;
1026 0 : for (iter = dd->template begin<TElem>(); iter != iterEnd; ++iter)
1027 : {
1028 : TElem* elem = *iter;
1029 0 : number elemError = aaError2[elem]; // eta_i^2
1030 :
1031 0 : if (elemError < 0) continue;
1032 0 : elemVar += (elemMean-elemError) * (elemMean-elemError);
1033 :
1034 :
1035 : }
1036 :
1037 : #ifdef UG_PARALLEL
1038 : if (pcl::NumProcs() > 1)
1039 : {
1040 : pcl::ProcessCommunicator com;
1041 : number elemVarLocal = elemVar;
1042 : elemVar = com.allreduce(elemVarLocal, PCL_RO_SUM);
1043 :
1044 : }
1045 : #endif
1046 : UG_LOG(" +++ VarianceMarking: Est. variance (1) : " << elemVar << " on "<< numElem << " elements.\n");
1047 :
1048 :
1049 0 : elemVar /= (numElem-1.0);
1050 : UG_LOG(" +++ VarianceMarking: Est. variance (2): " << elemVar << " on "<< numElem << " elements.\n");
1051 :
1052 :
1053 : // refine all element above threshold
1054 0 : const number sigma = sqrt(elemVar);
1055 0 : const number maxError = elemMean + sigma*m_width;
1056 0 : UG_LOG(" +++ Refining elements if error greater " << sigma << "*" << m_width <<" + "<< elemMean <<
1057 : " = " << maxError << ".\n");
1058 :
1059 :
1060 :
1061 0 : const number minErrToRefine = maxError*m_theta;
1062 0 : UG_LOG(" +++ Refining elements if error greater " << maxError << "*" << m_theta <<
1063 : " = " << minErrToRefine << ".\n");
1064 :
1065 : // reset counter
1066 : std::size_t numMarkedRefine = 0;
1067 :
1068 : // loop elements for marking
1069 0 : for (iter = dd->template begin<TElem>(); iter != iterEnd; ++iter)
1070 : {
1071 : // get element
1072 : TElem* elem = *iter;
1073 :
1074 : // if no error value exists: ignore (might be newly added by refinement);
1075 : // newly added elements are supposed to have a negative error estimator
1076 0 : if (aaError2[elem] < 0) continue;
1077 :
1078 : // marks for refinement
1079 0 : if (aaError2[elem] >= minErrToRefine)
1080 0 : if (dd->multi_grid()->get_level(elem) <= m_max_level)
1081 : {
1082 0 : refiner.mark(elem, RM_REFINE);
1083 0 : numMarkedRefine++;
1084 : }
1085 : }
1086 :
1087 : #ifdef UG_PARALLEL
1088 : if (pcl::NumProcs() > 1)
1089 : {
1090 : pcl::ProcessCommunicator com;
1091 : std::size_t numMarkedRefineLocal = numMarkedRefine;
1092 : numMarkedRefine = com.allreduce(numMarkedRefineLocal, PCL_RO_SUM);
1093 : UG_LOG(" +++ MaximumMarking: Marked for refinement: " << numMarkedRefine << " ("<< numMarkedRefineLocal << ") elements.\n");
1094 : }
1095 : #else
1096 : UG_LOG(" +++ MaximumMarking: Marked for refinement: " << numMarkedRefine << " elements.\n");
1097 : #endif
1098 :
1099 :
1100 0 : }
1101 :
1102 :
1103 : /// Marks elements above a certain threshold for refinement
1104 : /*! Threshold given by
1105 : * \f$ \theta * (\mu + width * \sigma) \f$
1106 : * where
1107 : * \f$ \mu = E[\eta], \sigma^2 = Var[\eta] \f$*
1108 : */
1109 : template <typename TDomain>
1110 : class VarianceMarkingEta : public IElementMarkingStrategy<TDomain>{
1111 :
1112 : public:
1113 : typedef IElementMarkingStrategy<TDomain> base_type;
1114 0 : VarianceMarkingEta(number theta) :
1115 0 : m_theta(theta), m_width(3.0), m_max_level(100), m_theta_coarse(0.0), m_min_level(0)
1116 : {};
1117 0 : VarianceMarkingEta(number theta, number width) :
1118 0 : m_theta(theta), m_width (width), m_max_level(100), m_theta_coarse(0.0), m_min_level(0)
1119 : {};
1120 0 : VarianceMarkingEta(number theta, number width, number theta_coarse) :
1121 0 : m_theta(theta), m_width (width), m_max_level(100), m_theta_coarse(theta_coarse), m_min_level(0)
1122 : {};
1123 :
1124 0 : void init_refinement(number theta, int max_level)
1125 0 : {m_theta = theta; m_max_level=max_level;}
1126 :
1127 0 : void init_coarsening(number theta, int min_level)
1128 0 : {m_theta_coarse = theta; m_min_level=min_level;}
1129 :
1130 : protected:
1131 : void mark(typename base_type::elem_accessor_type& aaError2,
1132 : IRefiner& refiner,
1133 : ConstSmartPtr<DoFDistribution> dd);
1134 : protected:
1135 :
1136 : number m_theta;
1137 : number m_width;
1138 : int m_max_level;
1139 :
1140 : number m_theta_coarse;
1141 : int m_min_level;
1142 : };
1143 :
1144 : template <typename TDomain>
1145 0 : void VarianceMarkingEta<TDomain>::mark(typename base_type::elem_accessor_type& aaError2,
1146 : IRefiner& refiner,
1147 : ConstSmartPtr<DoFDistribution> dd)
1148 : {
1149 : typedef typename base_type::elem_type TElem;
1150 : typedef typename DoFDistribution::traits<TElem>::const_iterator const_iterator;
1151 :
1152 : // compute minimal/maximal/ total error and number of elements
1153 :
1154 : number minElemErr, minElemErrLocal;
1155 : number maxElemErr, maxElemErrLocal;
1156 : number errSum, errTotalSq, errLocal;
1157 : size_t numElem, numElemLocal;
1158 :
1159 0 : const number elemMean = ComputeAvg(aaError2, dd, minElemErr, maxElemErr, errSum, errTotalSq, numElem,
1160 : minElemErrLocal, maxElemErrLocal, errLocal, numElemLocal);
1161 :
1162 :
1163 0 : this->m_latest_error = sqrt(errTotalSq);
1164 0 : this->m_latest_error_per_elem_max = maxElemErr;
1165 0 : this->m_latest_error_per_elem_min = minElemErr;
1166 :
1167 0 : UG_LOG(" +++ VarianceMarkingEta: error : "<< this->m_latest_error << " (meanEta : " << elemMean << " on "<< numElem << " elements).\n");
1168 :
1169 : // init iterators
1170 : const_iterator iter;
1171 0 : const const_iterator iterEnd = dd->template end<TElem>();
1172 :
1173 : number elemVar = 0.0;
1174 0 : for (iter = dd->template begin<TElem>(); iter != iterEnd; ++iter)
1175 : {
1176 : TElem* elem = *iter;
1177 :
1178 0 : if (aaError2[elem] < 0) continue;
1179 :
1180 0 : number elemError = sqrt(aaError2[elem]); // eta_i
1181 :
1182 0 : elemVar += (elemMean-elemError) * (elemMean-elemError);
1183 : }
1184 :
1185 0 : UG_LOG(" +++ VarianceMarkingEta: Est. variance (1) : " << (elemVar / (numElem -1.0)) << " on "<< numElem << " elements.\n");
1186 : #ifdef UG_PARALLEL
1187 : if (pcl::NumProcs() > 1)
1188 : {
1189 : pcl::ProcessCommunicator com;
1190 : number elemVarLocal = elemVar;
1191 : elemVar = com.allreduce(elemVarLocal, PCL_RO_SUM);
1192 :
1193 : }
1194 : #endif
1195 : elemVar /= (numElem-1.0);
1196 : UG_LOG(" +++ VarianceMarkingEta: Est. variance (2): " << elemVar << " on "<< numElem << " elements.\n");
1197 :
1198 :
1199 : // refine all element above threshold
1200 0 : const number sigma = sqrt(elemVar);
1201 0 : const number maxError = elemMean + sigma*m_width;
1202 0 : UG_LOG(" +++ Refining elements if error > " << elemMean << " + "<< sigma << "*" << m_width <<
1203 : " = " << maxError << ".\n");
1204 :
1205 0 : const number elemErrorRefine = maxError*m_theta;
1206 0 : UG_LOG(" +++ Refining elements if error > " << maxError << "*" << m_theta <<
1207 : " = " << elemErrorRefine << ".\n");
1208 :
1209 : // coarsen all element not contributing significantly
1210 : //const number minError = std::max(elemMean /*- sigma*m_width*/, 0.0);
1211 : //UG_LOG(" +++ Coarsening elements if error greater " << elemMean << " - "<< sigma << "*" << m_width <<
1212 : // " = " << minError << ".\n");
1213 :
1214 0 : const number elemErrorCoarsen = elemMean*m_theta_coarse;
1215 0 : UG_LOG(" +++ Coarsening elements if error < " << elemMean << "*" << m_theta_coarse <<
1216 : " = " << elemErrorCoarsen << ".\n");
1217 :
1218 :
1219 :
1220 : // reset counters
1221 : std::size_t numMarkedRefine = 0;
1222 : std::size_t numMarkedCoarsen = 0;
1223 :
1224 : // loop elements for marking
1225 0 : for (iter = dd->template begin<TElem>(); iter != iterEnd; ++iter)
1226 : {
1227 : // get element
1228 : TElem* elem = *iter;
1229 :
1230 : // if no error value exists: ignore (might be newly added by refinement);
1231 : // newly added elements are supposed to have a negative error estimator
1232 0 : if (aaError2[elem] < 0) continue;
1233 :
1234 0 : const number elemError = sqrt(aaError2[elem]); // eta_i
1235 :
1236 : // mark for refinement
1237 0 : if ((elemError >= elemErrorRefine) && (dd->multi_grid()->get_level(elem) <= m_max_level))
1238 : {
1239 0 : refiner.mark(elem, RM_REFINE);
1240 0 : numMarkedRefine++;
1241 : }
1242 :
1243 : // mark for coarsening
1244 0 : if ((elemError < elemErrorCoarsen) && (dd->multi_grid()->get_level(elem) >= m_min_level))
1245 : {
1246 0 : refiner.mark(elem, RM_COARSEN);
1247 0 : numMarkedCoarsen++;
1248 : }
1249 :
1250 : }
1251 :
1252 : #ifdef UG_PARALLEL
1253 : if (pcl::NumProcs() > 1)
1254 : {
1255 : pcl::ProcessCommunicator com;
1256 : std::size_t numMarkedRefineLocal = numMarkedRefine;
1257 : numMarkedRefine = com.allreduce(numMarkedRefineLocal, PCL_RO_SUM);
1258 : UG_LOG(" +++ VarianceMarkingEta: Marked for refinement: " << numMarkedRefine << " ("<< numMarkedRefineLocal << ") elements.\n");
1259 : } else
1260 : #endif
1261 : {
1262 : UG_LOG(" +++ VarianceMarkingEta: Marked for refinement: " << numMarkedRefine << " elements.\n");
1263 : UG_LOG(" +++ VarianceMarkingEta: Marked for coarsening: " << numMarkedCoarsen << " elements.\n");
1264 : }
1265 0 : }
1266 :
1267 :
1268 : /// marks elements above a certain fraction of the maximum
1269 : //!
1270 : template <typename TDomain>
1271 : class MeanValueMarking : public IElementMarkingStrategy<TDomain>{
1272 :
1273 : public:
1274 : typedef IElementMarkingStrategy<TDomain> base_type;
1275 0 : MeanValueMarking(number theta, number factor) : m_theta(theta), m_factor (factor) {};
1276 :
1277 : protected:
1278 : void mark(typename base_type::elem_accessor_type& aaError,
1279 : IRefiner& refiner,
1280 : ConstSmartPtr<DoFDistribution> dd);
1281 : protected:
1282 :
1283 : number m_theta;
1284 : number m_factor;
1285 : };
1286 :
1287 : template <typename TDomain>
1288 0 : void MeanValueMarking<TDomain>::mark(typename base_type::elem_accessor_type& aaError,
1289 : IRefiner& refiner,
1290 : ConstSmartPtr<DoFDistribution> dd)
1291 : {
1292 : typedef typename base_type::elem_type TElem;
1293 : typedef typename DoFDistribution::traits<TElem>::const_iterator const_iterator;
1294 :
1295 : // compute minimal/maximal/ total error and number of elements
1296 :
1297 : number minElemErr, minElemErrLocal;
1298 : number maxElemErrSq, maxElemErrLocal;
1299 : number errTotalSq, errLocal;
1300 : size_t numElem, numElemLocal;
1301 :
1302 0 : ComputeMinMax(aaError, dd, minElemErr, maxElemErrSq, errTotalSq, numElem,
1303 : minElemErrLocal, maxElemErrLocal, errLocal, numElemLocal);
1304 :
1305 :
1306 : // init iterators
1307 0 : number avgElemErr = errTotalSq / numElem;
1308 0 : UG_LOG(" +++ Global max: "<< maxElemErrSq << ", Global min: "<< minElemErr <<".\n");
1309 : UG_LOG(" +++ MeanValueMarking: Mean value : " << avgElemErr << " (Global error sum: "<<errTotalSq<<" on "<< numElem <<" elements).\n");
1310 :
1311 : // refine all element above threshold
1312 0 : const number minThetaErrToRefine = maxElemErrSq*m_theta;
1313 0 : const number minFactorAvgErrToRefine = avgElemErr*m_factor;
1314 0 : UG_LOG(" +++ MeanValueMarking: Min theta error : "<<minThetaErrToRefine<<" (Global Max Error: " << maxElemErrSq<< " for theta: "<<m_theta<<").\n");
1315 0 : UG_LOG(" +++ MeanValueMarking: Min factor avg error : "<<minFactorAvgErrToRefine<<" (Global Avg Error: " << avgElemErr<< " for factor: "<<m_factor<<").\n");
1316 :
1317 0 : const number minErrToRefine = std::min(minThetaErrToRefine, minFactorAvgErrToRefine);
1318 : UG_LOG(" +++ MeanValueMarking: Refining if error >= : "<<minErrToRefine<<".\n");
1319 :
1320 :
1321 : // reset counter
1322 : std::size_t numMarkedRefine = 0;
1323 :
1324 : // loop elements for marking
1325 : const_iterator iter;
1326 0 : const const_iterator iterEnd = dd->template end<TElem>();
1327 0 : for (iter = dd->template begin<TElem>(); iter != iterEnd; ++iter)
1328 : {
1329 : // get element
1330 : TElem* elem = *iter;
1331 :
1332 : // if no error value exists: ignore (might be newly added by refinement);
1333 : // newly added elements are supposed to have a negative error estimator
1334 0 : if (aaError[elem] < 0) continue;
1335 :
1336 : // marks for refinement
1337 0 : if (aaError[elem] >= minErrToRefine)
1338 : {
1339 0 : refiner.mark(elem, RM_REFINE);
1340 0 : numMarkedRefine++;
1341 : }
1342 : }
1343 :
1344 : #ifdef UG_PARALLEL
1345 : if (pcl::NumProcs() > 1)
1346 : {
1347 : pcl::ProcessCommunicator com;
1348 : std::size_t numMarkedRefineLocal = numMarkedRefine;
1349 : numMarkedRefine = com.allreduce(numMarkedRefineLocal, PCL_RO_SUM);
1350 : UG_LOG(" +++ MeanValueMarking: Marked for refinement: " << numMarkedRefine << " ("<< numMarkedRefineLocal << ") elements.\n");
1351 : }
1352 : #else
1353 : UG_LOG(" +++ MeanValueMarking: Marked for refinement: " << numMarkedRefine << " elements.\n");
1354 : #endif
1355 :
1356 :
1357 0 : }
1358 :
1359 :
1360 :
1361 :
1362 :
1363 :
1364 : /// marks elements above an absolute threshold (based on S. Reiter's idea)
1365 : template <typename TDomain>
1366 : class AbsoluteMarking : public IElementMarkingStrategy<TDomain>{
1367 :
1368 : public:
1369 : typedef IElementMarkingStrategy<TDomain> base_type;
1370 0 : AbsoluteMarking(number eta) : m_eta(eta), m_max_level(100) {};
1371 :
1372 : protected:
1373 : void mark(typename base_type::elem_accessor_type& aaError,
1374 : IRefiner& refiner,
1375 : ConstSmartPtr<DoFDistribution> dd);
1376 : protected:
1377 :
1378 : number m_eta;
1379 : int m_max_level;
1380 : };
1381 :
1382 : template <typename TDomain>
1383 0 : void AbsoluteMarking<TDomain>::mark(typename base_type::elem_accessor_type& aaError,
1384 : IRefiner& refiner,
1385 : ConstSmartPtr<DoFDistribution> dd)
1386 : {
1387 : typedef typename base_type::elem_type TElem;
1388 : typedef typename DoFDistribution::traits<TElem>::const_iterator const_iterator;
1389 :
1390 : // loop elements for marking
1391 0 : const const_iterator iterEnd = dd->template end<TElem>();
1392 :
1393 0 : number lowerBound = m_eta*m_eta;
1394 :
1395 : number errTotal=0.0;
1396 :
1397 : size_t numMarkedRefine = 0;
1398 0 : for (const_iterator iter = dd->template begin<TElem>(); iter != iterEnd; ++iter)
1399 : {
1400 : // get element error
1401 : TElem* elem = *iter;
1402 0 : const number elemError = aaError[elem];
1403 :
1404 : // skip newly added
1405 0 : if (elemError < 0) continue;
1406 0 : errTotal += elemError;
1407 :
1408 : // mark for refinement
1409 0 : if ((elemError >= lowerBound) && (dd->multi_grid()->get_level(elem) <= m_max_level))
1410 : {
1411 0 : refiner.mark(elem, RM_REFINE);
1412 0 : numMarkedRefine++;
1413 : }
1414 : }
1415 :
1416 :
1417 :
1418 : #ifdef UG_PARALLEL
1419 : if (pcl::NumProcs() > 1)
1420 : {
1421 : pcl::ProcessCommunicator com;
1422 : std::size_t numMarkedRefineLocal = numMarkedRefine;
1423 : number errLocal = errTotal;
1424 : numMarkedRefine = com.allreduce(numMarkedRefineLocal, PCL_RO_SUM);
1425 : errTotal = com.allreduce(errLocal, PCL_RO_SUM);
1426 : }
1427 : #endif
1428 :
1429 0 : this->m_latest_error = sqrt(errTotal);
1430 :
1431 : UG_LOG(" +++ AbsoluteMarking: Error**2 = "<<errTotal <<"; marked "<< numMarkedRefine << " elements for refinement "<<std::endl);
1432 :
1433 0 : }
1434 :
1435 :
1436 :
1437 : /// Mark surface layer for coarsening.
1438 : template <typename TDomain, typename TAlgebra>
1439 0 : void MarkForCoarsenening_SurfaceLayer(const GridFunction<TDomain, TAlgebra> &u, IRefiner& refiner){
1440 : typedef typename domain_traits<TDomain::dim>::element_type TElem;
1441 : ConstSmartPtr<DoFDistribution> spDD=u.dof_distribution();
1442 0 : refiner.mark(spDD->begin<TElem>(), spDD->end<TElem>(), RM_COARSEN);
1443 0 : refiner.coarsen();
1444 0 : }
1445 :
1446 :
1447 : }// end of namespace
1448 :
1449 : #endif
1450 :
|
