G+Smo  25.01.0
Geometry + Simulation Modules
 
Loading...
Searching...
No Matches
gsAdaptiveMeshing.hpp
Go to the documentation of this file.
1
14#pragma once
15
16
17#include <gsHSplines/gsHBSplineBasis.h>
18#include <gsHSplines/gsTHBSplineBasis.h>
19
20namespace gismo
21{
22
23template <class T>
24gsAdaptiveMeshing<T>::gsAdaptiveMeshing()
25{
26 defaultOptions();
27}
28
29template <class T>
30gsAdaptiveMeshing<T>::gsAdaptiveMeshing(gsFunctionSet<T> & input)
31:
32m_input(&input)
33{
34 defaultOptions();
35 rebuild();
36}
37
38template <class T>
39void gsAdaptiveMeshing<T>::rebuild()
40{
41 getOptions();
42
43 m_indices.clear();
44 m_boxes.clear();
45 this->_makeMap(m_input,m_indices,m_boxes);
46
47 bool check = true;
48 for (typename indexMapType::iterator it=m_indices.begin(); it!=m_indices.end(); it++)
49 check &= gsHBoxEqual<2,T>()(it->first,*m_boxes[it->second]);
50
51 for (typename boxMapType::iterator it=m_boxes.begin(); it!=m_boxes.end(); it++)
52 check &= it->first==m_indices[*it->second];
53
54 GISMO_ASSERT(check,"Something went wrong in the construction of the mappers");
55}
56
57template <class T>
58void gsAdaptiveMeshing<T>::_makeMap(const gsFunctionSet<T> * input, typename gsAdaptiveMeshing<T>::indexMapType & indexMap, typename gsAdaptiveMeshing<T>::boxMapType & boxMap)
59{
60 // typename gsAdaptiveMeshing<T>::indexMapType indexMap;
61 // typename gsAdaptiveMeshing<T>::boxMapType boxMap;
62
63 const gsBasis<T> * basis = nullptr;
64
65 // Make the container of boxes
66 // #pragma omp parallel
67 // {
68// #ifdef _OPENMP
69// const int tid = omp_get_thread_num();
70// const int nt = omp_get_num_threads();
71// index_t patch_cnt = 0;
72// #endif
73
74 index_t c = 0;
75 const gsMultiPatch<T> * mp;
76 const gsMultiBasis<T> * mb;
77 for (index_t patchInd=0; patchInd < input->nPieces(); ++patchInd)
78 {
79 // Initialize domain element iterator
80 if ( nullptr != (mp = dynamic_cast<const gsMultiPatch<T>*>(input)) ) basis = &(mp->basis(patchInd));
81 if ( nullptr != (mb = dynamic_cast<const gsMultiBasis<T>*>(input)) ) basis = &(mb->basis(patchInd));
82 GISMO_ENSURE(basis!=nullptr,"Object is not gsMultiBasis or gsMultiPatch");
83 // for all elements in patch pn
84 typename gsBasis<T>::domainIter domIt = basis->makeDomainIterator();
85 gsHDomainIterator<T,2> * domHIt = nullptr;
86 domHIt = dynamic_cast<gsHDomainIterator<T,2> *>(domIt.get());
87 GISMO_ENSURE(domHIt!=nullptr,"Domain not loaded");
88
89// #ifdef _OPENMP
90// c = patch_cnt + tid;
91// patch_cnt += domHIt->numElements();// a bit costy
92// for ( domHIt->next(tid); domHIt->good(); domHIt->next(nt) )
93// #else
94 for (; domHIt->good(); domHIt->next())
95// #endif
96 {
97 // #pragma omp critical (gsAdaptiveMeshingmakeBoxesinsert1)
98 {
99 HBox box(domHIt,patchInd);
100 std::pair<typename gsAdaptiveMeshing<T>::indexMapType::iterator,bool> mapIt = indexMap.insert({box,c});
101 if (mapIt.second)
102 {
103 // std::pair<typename gsAdaptiveMeshing<T>::boxMapType::iterator,bool> indexIt =
104 boxMap.insert({c,const_cast<gsHBox<2,real_t> *>(&(mapIt.first->first))});
105 // # ifdef _OPENMP
106 // c += nt;
107 // # else
108 c++;
109 // # endif
110
111 }
112 }
113 }
114 }
115 // }
116
117 // return std::make_pair(give(indexMap),give(boxMap));
118}
119
120
121
122template <class T>
123void gsAdaptiveMeshing<T>::_assignErrors(boxMapType & container, const std::vector<T> & elError)
124{
125 GISMO_ASSERT(elError.size()==container.size(),"The number of errors must be the same as the number of elements, but "<<elError.size()<<"!="<<container.size());
126 index_t k=0;
127
128 if (m_refRule == PBULK || m_crsRule == PBULK)
129 for (typename boxMapType::iterator it = container.begin(); it!=container.end(); it++, k++)
130 it->second->setAndProjectError(elError[k],m_alpha,m_beta);
131 else
132 for (typename boxMapType::iterator it = container.begin(); it!=container.end(); it++, k++)
133 it->second->setError(elError[k]);
134
135 m_totalError = _totalError(m_boxes);
136 m_maxError = _maxError(m_boxes);
137 if (m_alpha!=-1 && m_beta!=-1)
138 {
139 const gsBasis<T> * basis = nullptr;
140 index_t patchInd=0; // Assymes now that the minDegree is lowest on patch 0, or that the degree there is representative
141 const gsMultiPatch<T> * mp;
142 const gsMultiBasis<T> * mb;
143 if ( nullptr != (mp = dynamic_cast<const gsMultiPatch<T>*>(m_input)) ) basis = &(mp->basis(patchInd));
144 if ( nullptr != (mb = dynamic_cast<const gsMultiBasis<T>*>(m_input)) ) basis = &(mb->basis(patchInd));
145 GISMO_ASSERT(basis!=nullptr,"Object is not gsMultiBasis or gsMultiPatch");
146
147 m_uniformRefError = math::pow(1/2.,m_alpha * basis->minDegree() + m_beta) * m_totalError;
148 m_uniformCrsError = math::pow( 2.,m_alpha * basis->minDegree() + m_beta) * m_totalError;
149 }
150 else
151 {
152 m_uniformRefError = 0;
153 m_uniformCrsError = 0;
154 }
155}
156
199template <class T>
200template<bool _coarsen,bool _admissible>
201void gsAdaptiveMeshing<T>::_markElements(const std::vector<T> & elError, const index_t refCriterion, const std::vector<gsHBoxCheck<2,T> *> & predicates, HBoxContainer & elMarked) const
202{
203 GISMO_UNUSED(elError);
204 // Mark using different rules
205 switch (refCriterion)
206 {
207 case GARU:
208 _markThreshold<_coarsen,_admissible>(m_boxes,predicates,elMarked);
209 break;
210 case PUCA:
211 _markPercentage<_coarsen,_admissible>(m_boxes,predicates,elMarked);
212 break;
213 case BULK:
214 _markFraction<_coarsen,_admissible>(m_boxes,predicates,elMarked);
215 break;
216 case PBULK:
217 _markProjectedFraction<_coarsen,_admissible>(m_boxes,predicates,elMarked);
218 break;
219 default:
220 GISMO_ERROR("unknown marking strategy");
221 }
222}
223
224template <class T>
225void gsAdaptiveMeshing<T>::_crsPredicates_into(std::vector<gsHBoxCheck<2,T> *> & predicates)
226{
227 HBoxContainer empty;
228 predicates.push_back(new gsMinLvlCompare<2,T>(0));
229 predicates.push_back(new gsOverlapCompare<2,T>(empty,m_m));
230}
231
232template <class T>
233void gsAdaptiveMeshing<T>::_crsPredicates_into(const HBoxContainer & markedRef, std::vector<gsHBoxCheck<2,T> *> & predicates)
234{
235 predicates.push_back(new gsMinLvlCompare<2,T>(0));
236 predicates.push_back(new gsOverlapCompare<2,T>(markedRef,m_m));
237}
238
239template <class T>
240void gsAdaptiveMeshing<T>::_refPredicates_into(std::vector<gsHBoxCheck<2,T> *> & predicates)
241{
242 predicates.push_back(new gsMaxLvlCompare<2,T>(m_maxLvl));
243}
244
245
246template <class T>
247std::vector<index_t> gsAdaptiveMeshing<T>::_sortPermutation( const boxMapType & container)
248{
249 std::vector<index_t> idx(container.size());
250 std::iota(idx.begin(),idx.end(),0);
251 std::stable_sort(idx.begin(), idx.end(),
252 [&container](size_t i1, size_t i2) {
253 return container.at(i1)->error() < container.at(i2)->error();});
254
255 return idx;
256}
257
258template <class T>
259std::vector<index_t> gsAdaptiveMeshing<T>::_sortPermutationProjectedRef( const boxMapType & container)
260{
261 std::vector<index_t> idx(container.size());
262 std::iota(idx.begin(),idx.end(),0);
263 std::stable_sort(idx.begin(), idx.end(),
264 [&container](size_t i1, size_t i2) {
265 return container.at(i1)->projectedImprovement() < container.at(i2)->projectedImprovement();});
266
267 return idx;
268}
269
270template <class T>
271std::vector<index_t> gsAdaptiveMeshing<T>::_sortPermutationProjectedCrs( const boxMapType & container)
272{
273 std::vector<index_t> idx(container.size());
274 std::iota(idx.begin(),idx.end(),0);
275 std::stable_sort(idx.begin(), idx.end(),
276 [&container](size_t i1, size_t i2) {
277 return container.at(i1)->projectedSetBack() < container.at(i2)->projectedSetBack();});
278
279 return idx;
280}
281
282// template <class T>
283// void gsAdaptiveMeshing<T>::_sortPermutated( const std::vector<index_t> & permutation, boxContainer & container)
284// {
285// GISMO_ASSERT(permutation.size()==container.size(),"Permutation and vector should have the same size, but "<<permutation.size()<<"!="<<container.size());
286// boxContainer sorted(container.size());
287// for (size_t k=0; k!=container.size(); k++)
288// sorted.at(k) = container.at(permutation.at(k));
289
290// container.swap(sorted);
291// }
292
293template <class T>
294typename gsAdaptiveMeshing<T>::HBox * gsAdaptiveMeshing<T>::_boxPtr(const HBox & box) const
295{
296 // Fails if box is not in m_boxes
297 HBox * boxPtr = m_boxes.at(m_indices.at(box));
298 return boxPtr;
299 // return m_boxes[m_indices[box]];
300}
301
302template <class T>
303bool gsAdaptiveMeshing<T>::_checkBox( const HBox & box, const std::vector<gsHBoxCheck<2,T> *> predicates) const
304{
305 bool check = true;
306 for (typename std::vector<gsHBoxCheck<2,T>*>::const_iterator errIt = predicates.begin(); errIt!=predicates.end(); errIt++)
307 check &= (*errIt)->check(box);
308 return check;
309}
310
311template <class T>
312bool gsAdaptiveMeshing<T>::_checkBoxes( const typename HBox::Container & boxes, const std::vector<gsHBoxCheck<2,T> *> predicates) const
313{
314 bool check = true;
315 for (typename HBox::cIterator it = boxes.begin(); it!=boxes.end(); it++)
316 check &= _checkBox(*it,predicates);
317
318 return check;
319}
320
321// Use the mapTypes here!!
322
323
324template <class T>
325void gsAdaptiveMeshing<T>::_addAndMark( HBox & box, typename gsAdaptiveMeshing<T>::HBoxContainer & elMarked) const
326{
327 _boxPtr(box)->mark();
328 elMarked.add(*_boxPtr(box));
329}
330
331template <class T>
332void gsAdaptiveMeshing<T>::_addAndMark( typename HBox::Container & boxes, typename gsAdaptiveMeshing<T>::HBoxContainer & elMarked) const
333{
334 for (typename HBox::Iterator it = boxes.begin(); it!=boxes.end(); it++)
335 _addAndMark(*it,elMarked);
336}
337
338template <class T>
339void gsAdaptiveMeshing<T>::_setContainerProperties( typename HBox::Container & boxes ) const
340{
341 for (typename HBox::Iterator it = boxes.begin(); it!=boxes.end(); it++)
342 {
343 // HBox * box = _boxPtr(*it);
344 // it->setError(box->error());
345 // it->setMark(box->marked());
346 it->setAndProjectError(_boxPtr(*it)->error(),m_alpha,m_beta);
347 }
348}
349
350template <class T>
351T gsAdaptiveMeshing<T>::_totalError(const boxMapType & elements)
352{
353 // get total error
354 // Accumulation operator for boxMapType
355 auto accumulate_error_ptr = [](const T & val, const typename boxMapType::value_type & b)
356 { return val + b.second->error(); };
357 T totalError = std::accumulate(elements.begin(),elements.end(),(T)( 0 ),accumulate_error_ptr);
358 return totalError;
359}
360
361// Refinement: parameter % contributions to the total error of the highest cells are marked
362// Coarsening: parameter % contributions to the total error of the lowest cells are marked
363// In both cases, the total contribution is always lower than the threshold, i.e. if an element causes an exceed of the error, it is not taken into account
364
365template <class T>
366template<bool _coarsen,bool _admissible>
367typename std::enable_if< _coarsen && _admissible, void>::type
368gsAdaptiveMeshing<T>::_markFraction_impl( const boxMapType & elements, const std::vector<gsHBoxCheck<2,T> *> predicates, typename gsAdaptiveMeshing<T>::HBoxContainer & elMarked) const
369{
370 gsDebug<<"Fraction marking for coarsening...\n";
371 T cummulErrMarked = T(0);
372 T errorMarkSum = m_crsParam * m_totalError;
373
374 // Accumulation operator for neighborhoods (boxContainer)
375 auto accumulate_error = [](const T & val, const typename boxContainer::value_type & b)
376 { return val + b.error(); };
377
378 auto loop_action = [this,&elements,&predicates,&cummulErrMarked,&errorMarkSum,&elMarked,&accumulate_error]
379 (const index_t & index)
380 {
381 // Get the box
382 HBox_ptr box = elements.at(index);
383 // Continue if the box is already marked or if it does not satsfy all checks
384 if (_boxPtr(*box)->marked() || !_checkBox(*box,predicates))
385 return false;
386
387 // Get the neighborhoods
388 typename HBox::Container sibs = box->getSiblings();
389 _setContainerProperties(sibs);
390 if (!gsHBoxUtils<2,T>::allActive(sibs))
391 return false;
392
393 HBoxContainer siblings(sibs);
394 T siblingError = std::accumulate(sibs.begin(),sibs.end(),(T)( 0 ),accumulate_error);
395
396 // Check all siblings if they satisfy the predicates
397 if (_checkBoxes(sibs,predicates))
398 {
399 cummulErrMarked += siblingError;
400 _addAndMark(sibs,elMarked);
401 cummulErrMarked += box->error();
402 _addAndMark(*box,elMarked);
403 }
404 // return false;
405 return (cummulErrMarked > errorMarkSum);
406 };
407
408 std::vector<index_t>::const_iterator it = std::find_if(m_crsPermutation.cbegin(),m_crsPermutation.cend(),loop_action);
409 elMarked = HBoxUtils::Unique(elMarked);
410 gsDebug<<"[Mark fraction] Marked "<<elMarked.totalSize()<<" elements with marked error = "<<cummulErrMarked<<" and threshold = "<<errorMarkSum<<((it==m_crsPermutation.end()) ? " (maximum number marked)" : "")<<"\n";
411}
412
413
414template <class T>
415template<bool _coarsen,bool _admissible>
416typename std::enable_if< _coarsen && !_admissible, void>::type
417gsAdaptiveMeshing<T>::_markFraction_impl( const boxMapType & elements, const std::vector<gsHBoxCheck<2,T> *> predicates, typename gsAdaptiveMeshing<T>::HBoxContainer & elMarked) const
418{
419 gsDebug<<"Fraction marking for coarsening...\n";
420 T cummulErrMarked = T(0);
421 T errorMarkSum = m_crsParam * m_totalError;
422
423 auto loop_action = [this,&elements,&predicates,&cummulErrMarked,&errorMarkSum,&elMarked]
424 (const index_t & index)
425 {
426 // Get the box
427 HBox_ptr box = elements.at(index);
428
429 // Check box if it satisfy the predicates
430 if (_checkBox(*box,predicates))
431 {
432 cummulErrMarked += box->error();
433 _addAndMark(*box,elMarked);
434 }
435 // return false;
436 return (cummulErrMarked > errorMarkSum);
437 };
438
439 std::vector<index_t>::const_iterator it = std::find_if(m_crsPermutation.cbegin(),m_crsPermutation.cend(),loop_action);
440 elMarked = HBoxUtils::Unique(elMarked);
441 gsDebug<<"[Mark fraction] Marked "<<elMarked.totalSize()<<" elements with marked error = "<<cummulErrMarked<<" and threshold = "<<errorMarkSum<<((it==m_crsPermutation.end()) ? " (maximum number marked)" : "")<<"\n";
442}
443
444template <class T>
445template<bool _coarsen,bool _admissible>
446typename std::enable_if<!_coarsen && _admissible, void>::type
447gsAdaptiveMeshing<T>::_markFraction_impl( const boxMapType & elements, const std::vector<gsHBoxCheck<2,T> *> predicates, typename gsAdaptiveMeshing<T>::HBoxContainer & elMarked) const
448{
449 gsDebug<<"Fraction marking (admissible) for refinement...\n";
450 T cummulErrMarked = T(0);
451 T errorMarkSum = m_refParam * m_totalError;
452
453 // Accumulation operator for neighborhoods (boxContainer)
454 auto accumulate_error = [](const T & val, const typename boxContainer::value_type & b)
455 { return val + b.error(); };
456
457 auto loop_action = [this,&elements,&predicates,&cummulErrMarked,&errorMarkSum,&elMarked,&accumulate_error]
458 (const index_t & index)
459 {
460 // Get the box
461 HBox_ptr box = elements.at(index);
462 if (_boxPtr(*box)->marked())
463 {
464 return false;
465 }
466
467 // Get the neighborhoods
468 typename HBox::Container neighborhood = HBoxUtils::toContainer(HBoxUtils::markAdmissible(*box,m_m));
469 _setContainerProperties(neighborhood);
470 T neighborhoodError = std::accumulate(neighborhood.begin(),neighborhood.end(),(T)( 0 ),accumulate_error);
471
472 // if the errorMarkSum is exceeded with the current contribution, return true
473 // NOTE: this can mean that this element suddenly has a large contribution! i.e., larger then the next element in lin
474 // Maybe remove?
475 // if (cummulErrMarked + neighborhoodError > errorMarkSum)
476 // return true;
477
478 // Check all elements in the neighborhood if they satisfy the predicates
479 if (_checkBoxes(neighborhood,predicates))
480 {
481 cummulErrMarked += neighborhoodError;
482 _addAndMark(neighborhood,elMarked);
483 }
484 // return false;
485 return (cummulErrMarked > errorMarkSum);
486 };
487
488 std::vector<index_t>::const_iterator it = std::find_if(m_refPermutation.cbegin(),m_refPermutation.cend(),loop_action);
489 elMarked = HBoxUtils::Unique(elMarked);
490 gsDebug<<"[Mark fraction] Marked "<<elMarked.totalSize()<<" elements with marked error = "<<cummulErrMarked<<" and threshold = "<<errorMarkSum<<((it==m_refPermutation.end()) ? " (maximum number marked)" : "")<<"\n";
491}
492
493template <class T>
494template<bool _coarsen,bool _admissible>
495typename std::enable_if<!_coarsen && !_admissible, void>::type
496gsAdaptiveMeshing<T>::_markFraction_impl( const boxMapType & elements, const std::vector<gsHBoxCheck<2,T> *> predicates, typename gsAdaptiveMeshing<T>::HBoxContainer & elMarked) const
497{
498 gsDebug<<"Fraction marking (not admissible) for refinement...\n";
499 T cummulErrMarked = T(0);
500 T errorMarkSum = m_refParam * m_totalError;
501
502 auto loop_action = [this,&elements,&predicates,&cummulErrMarked,&errorMarkSum,&elMarked]
503 (const index_t & index)
504 {
505 HBox_ptr box = elements.at(index);
506
507 // if the errorMarkSum is exceeded with the current contribution, return true
508 // if (cummulErrMarked + box->error() > errorMarkSum)
509 // return true;
510
511 if (_checkBox(*box,predicates))
512 {
513 cummulErrMarked += box->error();
514 _addAndMark(*box,elMarked);
515 }
516 return (cummulErrMarked > errorMarkSum);
517 // return false;
518 };
519
520 std::vector<index_t>::const_iterator it = std::find_if(m_refPermutation.cbegin(),m_refPermutation.cend(),loop_action);
521 elMarked = HBoxUtils::Unique(elMarked);
522 gsDebug<<"[Mark fraction] Marked "<<elMarked.totalSize()<<" elements with marked error = "<<cummulErrMarked<<" and threshold = "<<errorMarkSum<<((it==m_refPermutation.end()) ? " (maximum number marked)" : "")<<"\n";
523}
524
525template <class T>
526template<bool _coarsen,bool _admissible>
527typename std::enable_if< _coarsen && _admissible, void>::type
528gsAdaptiveMeshing<T>::_markProjectedFraction_impl( const boxMapType & elements, const std::vector<gsHBoxCheck<2,T> *> predicates, typename gsAdaptiveMeshing<T>::HBoxContainer & elMarked) const
529{
530 gsDebug<<"Projected fraction (admissible) marking for coarsening...\n";
531 T targetError = m_crsParamExtra;
532 if (m_totalError > targetError)
533 return;
534 if (m_uniformCrsError < targetError) // then all elements should be refined
535 {
536 this->_markFraction_impl<_coarsen,_admissible>(elements,predicates,elMarked);
537 return;
538 }
539 T projectedError = m_totalError;
540
541 // Accumulation operator for neighborhoods (boxContainer)
542 auto accumulate_improvement = [](const T & val, const typename boxContainer::value_type & b)
543 { return val + b.projectedErrorCrs() - b.error(); };
544
545 auto loop_action = [this,&elements,&predicates,&projectedError,&targetError,&elMarked,&accumulate_improvement]
546 (const index_t & index)
547 {
548 // Get the box
549 HBox_ptr box = elements.at(index);
550 // Continue if the box is already marked or if it does not satsfy all checks
551 if (_boxPtr(*box)->marked() || !_checkBox(*box,predicates))
552 return false;
553
554 // Get the neighborhoods
555 typename HBox::Container sibs = box->getSiblings();
556 _setContainerProperties(sibs);
557 if (!gsHBoxUtils<2,T>::allActive(sibs))
558 return false;
559
560 HBoxContainer siblings(sibs);
561 T siblingSetBack = std::accumulate(sibs.begin(),sibs.end(),(T)( 0 ),accumulate_improvement);
562
563 // Check all siblings if they satisfy the predicates
564 if (_checkBoxes(sibs,predicates))
565 {
566 projectedError += siblingSetBack;
567 _addAndMark(sibs,elMarked);
568 projectedError += box->projectedErrorCrs() - box->error();
569 _addAndMark(*box,elMarked);
570 }
571
572 // return false;
573 return (projectedError > targetError);
574 };
575
576 std::vector<index_t>::const_iterator it = std::find_if(m_crsPermutation.cbegin(),m_crsPermutation.cend(),loop_action);
577 elMarked = HBoxUtils::Unique(elMarked);
578 gsDebug<<"[Mark projected fraction] Marked "<<elMarked.totalSize()<<" elements with projected error = "<<projectedError<<" and target error = "<<targetError<<((it==m_crsPermutation.end()) ? " (maximum number marked)" : "")<<"\n";
579}
580
581
582template <class T>
583template<bool _coarsen,bool _admissible>
584typename std::enable_if< _coarsen && !_admissible, void>::type
585//gsAdaptiveMeshing<T>::_markProjectedFraction_impl( const boxMapType & elements, const std::vector<gsHBoxCheck<2,T> *> predicates, typename gsAdaptiveMeshing<T>::HBoxContainer & elMarked) const
586gsAdaptiveMeshing<T>::_markProjectedFraction_impl( const boxMapType &, const std::vector<gsHBoxCheck<2,T> *>, typename gsAdaptiveMeshing<T>::HBoxContainer &) const
587{
588 GISMO_NO_IMPLEMENTATION
589 // gsDebug<<"Projected fraction marking for coarsening...\n";
590 // T projectedError = m_totalError;
591 // T targetError = m_crsParam;
592 // if (projectedError > targetError)
593 // return;
594
595 // auto loop_action = [this,&elements,&predicates,&projectedError,&targetError,&elMarked]
596 // (const index_t & index)
597 // {
598 // // Get the box
599 // HBox_ptr box = elements.at(index);
600
601 // // Check box if it satisfy the predicates
602 // if (_checkBox(*box,predicates))
603 // {
604 // projectedError += box->projectedImprovement();
605 // _addAndMark(*box,elMarked);
606 // }
607 // // return false;
608 // return (projectedError > targetError);
609 // };
610
611 // std::find_if(m_crsPermutation.cbegin(),m_crsPermutation.cend(),loop_action);
612 // elMarked = HBoxUtils::Unique(elMarked);
613 // gsDebug<<"[Mark fraction] Marked "<<elMarked.totalSize()<<" elements with projected error = "<<projectedError<<" and target error = "<<targetError<<"\n";
614}
615
616template <class T>
617template<bool _coarsen,bool _admissible>
618typename std::enable_if<!_coarsen && _admissible, void>::type
619gsAdaptiveMeshing<T>::_markProjectedFraction_impl( const boxMapType & elements, const std::vector<gsHBoxCheck<2,T> *> predicates, typename gsAdaptiveMeshing<T>::HBoxContainer & elMarked) const
620{
621 gsDebug<<"Projected fraction (admissible) marking for refinement...\n";
622 T targetError = m_refParamExtra;
623 if (m_totalError < targetError)
624 return;
625 if (m_uniformRefError > targetError) // then all elements should be refined
626 {
627 this->_markFraction_impl<_coarsen,_admissible>(elements,predicates,elMarked);
628 return;
629 }
630 T projectedError = m_totalError;
631
632 // Accumulation operator for neighborhoods (boxContainer)
633 auto accumulate_improvement = [](const T & val, const typename boxContainer::value_type & b)
634 {
635 return val + b.error() - b.projectedErrorRef();
636 };
637
638 auto loop_action = [this,&elements,&predicates,&projectedError,&targetError,&elMarked,&accumulate_improvement]
639 (const index_t & index)
640 {
641 // Get the box
642 HBox_ptr box = elements.at(index);
643 if (_boxPtr(*box)->marked())
644 {
645 return false;
646 }
647
648 // Get the neighborhoods
649 typename HBox::Container neighborhood = HBoxUtils::toContainer(HBoxUtils::markAdmissible(*box,m_m));
650 HBoxContainer neighborhoodtmp = HBoxContainer(neighborhood);
651 _setContainerProperties(neighborhood);
652 T neighborhoodImprovement = std::accumulate(neighborhood.begin(),neighborhood.end(),(T)( 0 ),accumulate_improvement);
653
654 // Check all elements in the neighborhood if they satisfy the predicates
655 if (_checkBoxes(neighborhood,predicates))
656 {
657 projectedError -= neighborhoodImprovement;
658 _addAndMark(neighborhood,elMarked);
659 }
660 // return false;
661 return (projectedError < targetError);
662 };
663
664 std::vector<index_t>::const_iterator it = std::find_if(m_refPermutation.cbegin(),m_refPermutation.cend(),loop_action);
665 elMarked = HBoxUtils::Unique(elMarked);
666 gsDebug<<"[Mark projected fraction] Marked "<<elMarked.totalSize()<<" elements with projected error = "<<projectedError<<" and target error = "<<targetError<<((it==m_refPermutation.end()) ? " (maximum number marked)" : "")<<"\n";
667}
668
669template <class T>
670template<bool _coarsen,bool _admissible>
671typename std::enable_if<!_coarsen && !_admissible, void>::type
672gsAdaptiveMeshing<T>::_markProjectedFraction_impl( const boxMapType & elements, const std::vector<gsHBoxCheck<2,T> *> predicates, typename gsAdaptiveMeshing<T>::HBoxContainer & elMarked) const
673{
674 gsDebug<<"Projected fraction (not admissible) marking for refinement...\n";
675 T targetError = m_refParamExtra;
676 if (m_totalError < targetError)
677 return;
678 if (m_uniformRefError > targetError) // then all elements should be refined
679 {
680 this->_markFraction_impl<_coarsen,_admissible>(elements,predicates,elMarked);
681 return;
682 }
683 T projectedError = m_totalError;
684
685 auto loop_action = [this,&elements,&predicates,&projectedError,&targetError,&elMarked]
686 (const index_t & index)
687 {
688 HBox_ptr box = elements.at(index);
689
690 // if the errorMarkSum is exceeded with the current contribution, return true
691 // if (cummulErrMarked + box->error() > errorMarkSum)
692 // return true;
693
694 if (_checkBox(*box,predicates))
695 {
696 projectedError -= box->projectedImprovement();
697 _addAndMark(*box,elMarked);
698 }
699 return (projectedError < targetError);
700 // return false;
701 };
702
703 std::vector<index_t>::const_iterator it = std::find_if(m_refPermutation.cbegin(),m_refPermutation.cend(),loop_action);
704 elMarked = HBoxUtils::Unique(elMarked);
705 gsDebug<<"[Mark projected fraction] Marked "<<elMarked.totalSize()<<" elements with projected error = "<<projectedError<<" and target error = "<<targetError<<((it==m_refPermutation.end()) ? " (maximum number marked)" : "")<<"\n";
706}
707
708template <class T>
709template<bool _coarsen,bool _admissible>
710typename std::enable_if< _coarsen && _admissible, void>::type
711gsAdaptiveMeshing<T>::_markPercentage_impl( const boxMapType & elements, const std::vector<gsHBoxCheck<2,T> *> predicates, typename gsAdaptiveMeshing<T>::HBoxContainer & elMarked) const
712{
713 gsDebug<<"Percentage marking for coarsening...\n";
714 // Total number of elements:
715 size_t NE = elements.size();
716 // Compute the index from which the refinement should start,
717 // once the vector is sorted.
718 index_t NR = cast<T,index_t>( math::floor( m_crsParam * T(NE) ) );
719
720 index_t nmarked = 0;
721 auto loop_action = [this,&elements,&predicates,&NR,&nmarked,&elMarked]
722 (const index_t & index)
723 {
724 HBox * box = elements.at(index);
725
726 // Continue if the box is already marked or if it does not satsfy all checks
727 if (_boxPtr(*box)->marked() || !_checkBox(*box,predicates))
728 return false;
729
730 // Get the neighborhoods
731 typename HBox::Container sibs = box->getSiblings();
732 _setContainerProperties(sibs);
733 if (!gsHBoxUtils<2,T>::allActive(sibs))
734 return false;
735
736 // Check all children if they satisfy the predicates
737 if (_checkBoxes(sibs,predicates))
738 {
739 nmarked += sibs.size();
740 _addAndMark(sibs,elMarked);
741 nmarked += 1;
742 _addAndMark(*box,elMarked);
743 }
744 return (nmarked > NR);
745 };
746
747 std::vector<index_t>::const_iterator it = std::find_if(m_crsPermutation.cbegin(),m_crsPermutation.cend(),loop_action);
748 elMarked = HBoxUtils::Unique(elMarked);
749 gsDebug<<"[Mark percentage] Marked "<<elMarked.totalSize()<<", ("<<nmarked<<") elements ("<<(T)nmarked/NE*100<<"%"<<" of NE "<<NE<<") and threshold = "<<NR<<" ("<<m_crsParam*100<<"%)"<<((it==m_crsPermutation.end()) ? " (maximum number marked)" : "")<<"\n";
750}
751
752template <class T>
753template<bool _coarsen,bool _admissible>
754typename std::enable_if< _coarsen && !_admissible, void>::type
755gsAdaptiveMeshing<T>::_markPercentage_impl( const boxMapType & elements, const std::vector<gsHBoxCheck<2,T> *> predicates, typename gsAdaptiveMeshing<T>::HBoxContainer & elMarked) const
756{
757 gsDebug<<"Percentage marking for coarsening...\n";
758 // Total number of elements:
759 size_t NE = elements.size();
760 // Compute the index from which the refinement should start,
761 // once the vector is sorted.
762 index_t NR = cast<T,index_t>( math::floor( m_crsParam * T(NE) ) );
763
764 index_t nmarked = 0;
765 auto loop_action = [this,&elements,&predicates,&NR,&nmarked,&elMarked]
766 (const index_t & index)
767 {
768 HBox * box = elements.at(index);
769
770 if (_checkBox(*box,predicates))
771 {
772 nmarked += 1;
773 _addAndMark(*box,elMarked);
774 }
775 return (nmarked > NR);
776 };
777
778 std::vector<index_t>::const_iterator it = std::find_if(m_crsPermutation.cbegin(),m_crsPermutation.cend(),loop_action);
779 elMarked = HBoxUtils::Unique(elMarked);
780 gsDebug<<"[Mark percentage] Marked "<<elMarked.totalSize()<<", ("<<nmarked<<") elements ("<<(T)nmarked/NE*100<<"%"<<" of NE "<<NE<<") and threshold = "<<NR<<" ("<<m_crsParam*100<<"%)"<<((it==m_crsPermutation.end()) ? " (maximum number marked)" : "")<<"\n";
781}
782
783template <class T>
784template<bool _coarsen,bool _admissible>
785typename std::enable_if<!_coarsen && _admissible, void>::type
786gsAdaptiveMeshing<T>::_markPercentage_impl( const boxMapType & elements, const std::vector<gsHBoxCheck<2,T> *> predicates, typename gsAdaptiveMeshing<T>::HBoxContainer & elMarked) const
787{
788 gsDebug<<"Percentage marking (admissible) for refinement...\n";
789 // Total number of elements:
790 size_t NE = elements.size();
791 // Compute the index from which the refinement should start,
792 // once the vector is sorted.
793 index_t NR = cast<T,index_t>( math::floor( m_refParam * T(NE) ) );
794
795 index_t nmarked = 0;
796 auto loop_action = [this,&elements,&predicates,&nmarked,&NR,&elMarked]
797 (const index_t & index)
798 {
799 // Get the box
800 HBox_ptr box = elements.at(index);
801
802 if (_boxPtr(*box)->marked())
803 return false;
804
805 typename HBox::Container neighborhood = HBoxUtils::toContainer(HBoxUtils::markAdmissible(*box,m_m));
806 _setContainerProperties(neighborhood);
807 // Check all elements in the neighborhood if they satisfy the predicates
808 if (_checkBoxes(neighborhood,predicates))
809 {
810 nmarked += neighborhood.size();
811 _addAndMark(neighborhood,elMarked);
812 }
813 return (nmarked > NR);
814 };
815
816 std::vector<index_t>::const_iterator it = std::find_if(m_refPermutation.cbegin(),m_refPermutation.cend(),loop_action);
817 elMarked = HBoxUtils::Unique(elMarked);
818 gsDebug<<"[Mark percentage] Marked "<<elMarked.totalSize()<<", ("<<nmarked<<") elements ("<<(T)nmarked/NE*100<<"%"<<" of NE "<<NE<<") and threshold = "<<NR<<" ("<<m_refParam*100<<"%)"<<((it==m_refPermutation.end()) ? " (maximum number marked)" : "")<<"\n";
819}
820
821template <class T>
822template<bool _coarsen,bool _admissible>
823typename std::enable_if<!_coarsen && !_admissible, void>::type
824gsAdaptiveMeshing<T>::_markPercentage_impl( const boxMapType & elements, const std::vector<gsHBoxCheck<2,T> *> predicates, typename gsAdaptiveMeshing<T>::HBoxContainer & elMarked) const
825{
826 gsDebug<<"Percentage marking (not admissible) for refinement...\n";
827 // Total number of elements:
828 size_t NE = elements.size();
829 // Compute the index from which the refinement should start,
830 // once the vector is sorted.
831 index_t NR = cast<T,index_t>( math::floor( m_refParam * T(NE) ) );
832
833 index_t nmarked = 0;
834 auto loop_action = [&elements,&predicates,&NR,&nmarked,&elMarked]
835 (const index_t & index)
836 {
837 HBox * box = elements.at(index);
838 bool check = true;
839 for (typename std::vector<gsHBoxCheck<2,T>*>::const_iterator errIt = predicates.begin(); errIt!=predicates.end(); errIt++)
840 check &= (*errIt)->check(*box);
841
842 if (check)
843 {
844 elMarked.add(*box);
845 nmarked++;
846 }
847
848 return (nmarked > NR);
849 };
850
851 std::vector<index_t>::const_iterator it = std::find_if(m_refPermutation.cbegin(),m_refPermutation.cend(),loop_action);
852 elMarked = HBoxUtils::Unique(elMarked);
853 gsDebug<<"[Mark percentage] Marked "<<elMarked.totalSize()<<", ("<<nmarked<<") elements ("<<(T)nmarked/NE*100<<"%"<<" of NE "<<NE<<") and threshold = "<<NR<<" ("<<m_refParam*100<<"%)"<<((it==m_refPermutation.end()) ? " (maximum number marked)" : "")<<"\n";
854}
855
856template <class T>
857T gsAdaptiveMeshing<T>::_maxError(const boxMapType & elements)
858{
859 auto larger_than = [](const typename boxMapType::value_type & a, const typename boxMapType::value_type & b)
860 {
861 return (a.second->error() < b.second->error());
862 };
863
864 // First, conduct a brutal search for the maximum local error
865 T maxErr = (std::max_element(elements.begin(), elements.end(), larger_than ))->second->error();
866
867 return maxErr;
868}
869
870template <class T>
871template<bool _coarsen,bool _admissible>
872typename std::enable_if< _coarsen && _admissible, void>::type
873gsAdaptiveMeshing<T>::_markThreshold_impl( const boxMapType & elements, const std::vector<gsHBoxCheck<2,T> *> predicates, typename gsAdaptiveMeshing<T>::HBoxContainer & elMarked) const
874{
875 gsDebug<<"Threshold marking for coarsening...\n";
876 GISMO_ASSERT(m_crsParam<=1 && m_crsParam>=0,"Coarsening parameter must be a percentage!");
877
878 T Thr = m_crsParam * m_maxError;
879 T current = 0;
880 gsHBoxCheck<2,T> * thres_predicate = new gsLargerErrCompare<2,T>(Thr);
881
882 auto loop_action = [this,&elements,&thres_predicate,&predicates,&elMarked,&current]
883 (const index_t & index)
884 {
885 HBox_ptr box = elements.at(index);
886
887 if (thres_predicate->check(*box))
888 {
889 current = std::max(current,box->error());
890 return true;
891 }
892
893 // Continue if the box is already marked or if it does not satsfy all checks
894 if (_boxPtr(*box)->marked() || !_checkBox(*box,predicates))
895 return false;
896
897
898 // Get the neighborhoods
899 typename HBox::Container sibs = box->getSiblings();
900 _setContainerProperties(sibs);
901 if (!gsHBoxUtils<2,T>::allActive(sibs))
902 return false;
903
904 // Check all siblings if they satisfy the predicates
905 if (_checkBoxes(sibs,predicates))
906 {
907 _addAndMark(sibs,elMarked);
908 _addAndMark(*box,elMarked);
909 }
910 return false;
911 };
912
913 std::vector<index_t>::const_iterator it = std::find_if(m_crsPermutation.cbegin(),m_crsPermutation.cend(),loop_action);
914 elMarked = HBoxUtils::Unique(elMarked);
915 gsDebug<<"[Mark threshold] Marked "<<elMarked.totalSize()<<" elements with largest error "<<current<<" and treshold = "<<Thr<<((it==m_crsPermutation.end()) ? " (maximum number marked)" : "")<<"\n";
916}
917
918template <class T>
919template<bool _coarsen,bool _admissible>
920typename std::enable_if< _coarsen && !_admissible, void>::type
921gsAdaptiveMeshing<T>::_markThreshold_impl( const boxMapType & elements, const std::vector<gsHBoxCheck<2,T> *> predicates, typename gsAdaptiveMeshing<T>::HBoxContainer & elMarked) const
922{
923 gsDebug<<"Threshold marking for coarsening...\n";
924 GISMO_ASSERT(m_crsParam<=1 && m_crsParam>=0,"Coarsening parameter must be a percentage!");
925
926 T Thr = m_crsParam * m_maxError;
927 T current = 0;
928 gsHBoxCheck<2,T> * thres_predicate = new gsLargerErrCompare<2,T>(Thr);
929
930 auto loop_action = [this,&elements,&thres_predicate,&predicates,&elMarked,&current]
931 (const index_t & index)
932 {
933 HBox_ptr box = elements.at(index);
934
935 if (thres_predicate->check(*box))
936 {
937 current = std::max(current,box->error());
938 return true;
939 }
940
941 // Check all children if they satisfy the predicates
942 if (_checkBox(*box,predicates))
943 _addAndMark(*box,elMarked);
944 return false;
945 };
946
947 std::vector<index_t>::const_iterator it = std::find_if(m_crsPermutation.cbegin(),m_crsPermutation.cend(),loop_action);
948 elMarked = HBoxUtils::Unique(elMarked);
949 gsDebug<<"[Mark threshold] Marked "<<elMarked.totalSize()<<" elements with largest error "<<current<<" and treshold = "<<Thr<<((it==m_crsPermutation.end()) ? " (maximum number marked)" : "")<<"\n";
950}
951
952template <class T>
953template<bool _coarsen,bool _admissible>
954typename std::enable_if<!_coarsen && _admissible, void>::type
955gsAdaptiveMeshing<T>::_markThreshold_impl( const boxMapType & elements, const std::vector<gsHBoxCheck<2,T> *> predicates, typename gsAdaptiveMeshing<T>::HBoxContainer & elMarked) const
956{
957 gsDebug<<"Threshold marking (admissible) for refinement...\n";
958 GISMO_ASSERT(m_refParam<=1 && m_refParam>=0,"Refinement parameter must be a percentage!");
959
960 T Thr = m_refParam * m_maxError;
961 T current = 0;
962 gsHBoxCheck<2,T> * thres_predicate = new gsSmallerErrCompare<2,T>(Thr);
963
964 auto loop_action = [this,&elements,&thres_predicate,&predicates,&elMarked,&current]
965 (const index_t & index)
966 {
967 HBox_ptr box = elements.at(index);
968
969 if (thres_predicate->check(*box))
970 {
971 current = std::max(current,box->error());
972 return true;
973 }
974
975 if (_boxPtr(*box)->marked())
976 {
977 return false;
978 }
979
980 typename HBox::Container neighborhood = HBoxUtils::toContainer(HBoxUtils::markAdmissible(*box,m_m));
981 _setContainerProperties(neighborhood);
982
983 // Check all elements in the neighborhood if they satisfy the predicates
984 if (_checkBoxes(neighborhood,predicates))
985 _addAndMark(neighborhood,elMarked);
986 return false;
987 };
988
989 std::vector<index_t>::const_iterator it = std::find_if(m_refPermutation.cbegin(),m_refPermutation.cend(),loop_action);
990 elMarked = HBoxUtils::Unique(elMarked);
991 gsDebug<<"[Mark threshold] Marked "<<elMarked.totalSize()<<" elements with largest error "<<current<<" and treshold = "<<Thr<<((it==m_refPermutation.end()) ? " (maximum number marked)" : "")<<"\n";
992}
993
994template <class T>
995template<bool _coarsen,bool _admissible>
996typename std::enable_if<!_coarsen && !_admissible, void>::type
997gsAdaptiveMeshing<T>::_markThreshold_impl( const boxMapType & elements, const std::vector<gsHBoxCheck<2,T> *> predicates, typename gsAdaptiveMeshing<T>::HBoxContainer & elMarked) const
998{
999 gsDebug<<"Threshold marking (not admissible) for refinement...\n";
1000 GISMO_ASSERT(m_refParam<=1 && m_refParam>=0,"Refinement parameter must be a percentage!");
1001
1002 T Thr = m_refParam * m_maxError;
1003 T current = 0;
1004 gsHBoxCheck<2,T> * thres_predicate = new gsSmallerErrCompare<2,T>(Thr);
1005
1006 auto loop_action = [this,&elements,&thres_predicate,&predicates,&elMarked,&current]
1007 (const index_t & index)
1008 {
1009 HBox_ptr box = elements.at(index);
1010
1011 if (thres_predicate->check(*box))
1012 {
1013 current = std::max(current,box->error());
1014 return true;
1015 }
1016
1017 if (_checkBox(*box,predicates))
1018 _addAndMark(*box,elMarked);
1019
1020 return false;
1021 };
1022
1023 std::vector<index_t>::const_iterator it = std::find_if(m_refPermutation.cbegin(),m_refPermutation.cend(),loop_action);
1024 elMarked = HBoxUtils::Unique(elMarked);
1025 gsDebug<<"[Mark threshold] Marked "<<elMarked.totalSize()<<" elements with largest error "<<current<<" and treshold = "<<Thr<<((it==m_refPermutation.end()) ? " (maximum number marked)" : "")<<"\n";
1026}
1027
1028template<class T>
1029void gsAdaptiveMeshing<T>::defaultOptions()
1030{
1031 m_options.addInt("CoarsenRule","Rule used for coarsening: 1=GARU, 2=PUCA, 3=BULK.",1);
1032 m_options.addInt("RefineRule","Rule used for refinement: 1=GARU, 2=PUCA, 3=BULK.",1);
1033 m_options.addReal("CoarsenParam","Parameter used for coarsening",0.1);
1034 m_options.addReal("RefineParam","Parameter used for refinement",0.1);
1035 m_options.addReal("CoarsenParamExtra","Extra parameter used for coarsening",0.1);
1036 m_options.addReal("RefineParamExtra","Extra parameter used for refinement",0.1);
1037
1038 m_options.addInt("Convergence_alpha","Estimated convergence parameter of he error, for alpha*p+beta convergence",-1);
1039 m_options.addInt("Convergence_beta","Estimated convergence parameter of he error, for alpha*p+beta convergence",-1);
1040
1041 m_options.addInt("CoarsenExtension","Extension coarsening",0);
1042 m_options.addInt("RefineExtension","Extension refinement",0);
1043
1044 m_options.addInt("MaxLevel","Maximum refinement level",3);
1045
1046 m_options.addInt("Admissibility","Admissibility region, 0=T-admissibility (default), 1=H-admissibility",0);
1047 m_options.addSwitch("Admissible","Mark the admissible region",true);
1048 m_options.addInt("Jump","Jump parameter m",2);
1049}
1050
1051template<class T>
1052void gsAdaptiveMeshing<T>::getOptions()
1053{
1054 switch (m_options.askInt("CoarsenRule",3))
1055 {
1056 case 1:
1057 m_crsRule = GARU;
1058 break;
1059 case 2:
1060 m_crsRule = PUCA;
1061 break;
1062 case 3:
1063 m_crsRule = BULK;
1064 break;
1065 case 4:
1066 m_crsRule = PBULK;
1067 break;
1068 default:
1069 GISMO_ERROR("Coarsening marking strategy unknown");
1070 break;
1071 }
1072
1073 switch (m_options.askInt("RefineRule",3))
1074 {
1075 case 1:
1076 m_refRule = GARU;
1077 break;
1078 case 2:
1079 m_refRule = PUCA;
1080 break;
1081 case 3:
1082 m_refRule = BULK;
1083 break;
1084 case 4:
1085 m_refRule = PBULK;
1086 break;
1087 default:
1088 GISMO_ERROR("Refinement marking strategy unknown");
1089 break;
1090 }
1091
1092 m_crsParam = m_options.getReal("CoarsenParam");
1093 m_crsParamExtra = m_options.getReal("CoarsenParamExtra");
1094 m_refParam = m_options.getReal("RefineParam");
1095 m_refParamExtra = m_options.getReal("RefineParamExtra");
1096
1097 m_crsExt = m_options.getInt("CoarsenExtension");
1098 m_refExt = m_options.getInt("RefineExtension");
1099
1100 m_maxLvl = m_options.getInt("MaxLevel");
1101
1102 m_admissible = m_options.getSwitch("Admissible");
1103 m_m = m_options.getInt("Jump");
1104
1105 m_alpha=m_options.askInt("Convergence_alpha",-1);
1106 m_beta=m_options.askInt("Convergence_beta",-1);
1107}
1108
1109template<class T>
1110void gsAdaptiveMeshing<T>::container_into(const std::vector<T> & elError, HBoxContainer & result)
1111{
1112 result.clear();
1113 this->_assignErrors(m_boxes,elError);
1114 for (typename boxMapType::iterator it = m_boxes.begin(); it!=m_boxes.end(); it++)
1115 result.add(*it->second);
1116}
1117
1118
1119template<class T>
1120void gsAdaptiveMeshing<T>::markRef_into(const std::vector<T> & elError, HBoxContainer & elMarked)
1121{
1122 elMarked.clear();
1123 this->_assignErrors(m_boxes,elError);
1124 if (m_refRule!=PBULK)
1125 m_refPermutation = this->_sortPermutation(m_boxes); // Index of the lowest error is first
1126 else
1127 m_refPermutation = this->_sortPermutationProjectedRef(m_boxes); // Index of the lowest error is first
1128
1129 // To do:
1130 // - sort the errors for coarsened and refined in separate functions
1131 // try again
1132 // - computeError(primalL,dualL,dualH) class in gsThinSHellAssemblerDWR which derives from gsFunctionSet
1133
1134
1135 std::reverse(m_refPermutation.begin(),m_refPermutation.end()); // Index of the highest error is first
1136
1137 std::vector<gsHBoxCheck<2,T> *> predicates;
1138 _refPredicates_into(predicates);
1139
1140 if (m_admissible)
1141 _markElements<false,true>( elError, m_refRule, predicates, elMarked);//,flag [coarse]);
1142 else
1143 _markElements<false,false>( elError, m_refRule, predicates, elMarked);//,flag [coarse]);
1144
1145 for (typename std::vector<gsHBoxCheck<2,T>*>::iterator pred=predicates.begin(); pred!=predicates.end(); pred++)
1146 delete *pred;
1147}
1148
1149template<class T>
1150void gsAdaptiveMeshing<T>::markCrs_into(const std::vector<T> & elError, const HBoxContainer & markedRef, HBoxContainer & elMarked)
1151{
1152 elMarked.clear();
1153 this->_assignErrors(m_boxes,elError);
1154
1155 if (m_crsRule!=PBULK)
1156 m_crsPermutation = this->_sortPermutation(m_boxes); // Index of the lowest error is first
1157 else
1158 m_crsPermutation = this->_sortPermutationProjectedCrs(m_boxes); // Index of the lowest error is first
1159 gsDebugVar(m_crsPermutation.size());
1160
1161 std::vector<gsHBoxCheck<2,T> *> predicates;
1162 if (markedRef.totalSize()==0 || !m_admissible)
1163 _crsPredicates_into(predicates);
1164 else
1165 _crsPredicates_into(markedRef,predicates);
1166
1167 if (m_admissible)
1168 _markElements<true,true>( elError, m_crsRule, predicates, elMarked);//,flag [coarse]);
1169 else
1170 _markElements<true,false>( elError, m_crsRule, predicates, elMarked);//,flag [coarse]);
1171
1172 gsDebugVar(m_crsPermutation.size());
1173
1174 for (typename std::vector<gsHBoxCheck<2,T>*>::iterator pred=predicates.begin(); pred!=predicates.end(); pred++)
1175 delete *pred;
1176}
1177
1178template<class T>
1179void gsAdaptiveMeshing<T>::markCrs_into(const std::vector<T> & elError, HBoxContainer & elMarked)
1180{
1181 HBoxContainer container;
1182 this->markCrs_into(elError,container,elMarked);
1183}
1184
1185template<class T>
1186void gsAdaptiveMeshing<T>::markRef(const std::vector<T> & errors)
1187{
1188 markRef_into( errors, m_markedRef);//,flag [coarse]);
1189}
1190
1191template<class T>
1192void gsAdaptiveMeshing<T>::markCrs(const std::vector<T> & errors)
1193{
1194 markCrs_into( errors, m_markedRef, m_markedCrs);//,flag [coarse]);
1195}
1196
1197template<class T>
1198bool gsAdaptiveMeshing<T>::refine(const HBoxContainer & markedRef)
1199{
1200 bool refine = (markedRef.totalSize()>0);
1201 if ( refine )
1202 _refineMarkedElements(markedRef,m_refExt);
1203 return refine;
1204}
1205
1206template<class T>
1207bool gsAdaptiveMeshing<T>::unrefine(const HBoxContainer & markedCrs)
1208{
1209 bool coarsen= (markedCrs.totalSize()>0);
1210 if (coarsen)
1211 _unrefineMarkedElements(markedCrs,m_crsExt);
1212 return coarsen;
1213}
1214
1215template<class T>
1216bool gsAdaptiveMeshing<T>::refineAll()
1217{
1218 HBoxContainer ref;
1219 for (typename boxMapType::iterator it = m_boxes.begin(); it!=m_boxes.end(); it++)
1220 ref.add(*it->second);
1221
1222 this->refine(ref);
1223
1224 return true;
1225}
1226
1227template<class T>
1228bool gsAdaptiveMeshing<T>::unrefineAll()
1229{
1230 HBoxContainer crs;
1231 for (typename boxMapType::iterator it = m_boxes.begin(); it!=m_boxes.end(); it++)
1232 crs.add(*it->second);
1233
1234 this->unrefine(crs);
1235
1236 return true;
1237}
1238
1239// template<class T>
1240// void gsAdaptiveMeshing<T>::flatten(const index_t level)
1241// {
1242// _flattenElementsToLevel(level);
1243// }
1244
1245// template<class T>
1246// void gsAdaptiveMeshing<T>::unrefineThreshold(const index_t level)
1247// {
1248// _unrefineElementsThreshold(level);
1249// }
1250
1251template<class T>
1252void gsAdaptiveMeshing<T>::_refineMarkedElements( const HBoxContainer & markedRef,
1253 index_t refExtension)
1254{
1255 gsBasis<T> * basis = nullptr;
1256
1257 gsMultiPatch<T> * mp;
1258 gsMultiBasis<T> * mb;
1259
1260 for (index_t pn=0; pn < m_input->nPieces(); ++pn )// for all patches
1261 {
1262 if ( (mp = dynamic_cast<gsMultiPatch<T>*>(m_input)) != nullptr ) basis = &(mp->basis(pn));
1263 if ( (mb = dynamic_cast<gsMultiBasis<T>*>(m_input)) != nullptr ) basis = &(mb->basis(pn));
1264 GISMO_ENSURE(basis!=nullptr,"Object is not gsMultiBasis or gsMultiPatch");
1265
1266 // if (m_options.getSwitch("Admissible"))
1267 // {
1268 // if (m_options.getInt("Admissibility")==0)
1269 // {
1270 // if ( gsTHBSplineBasis<2,T> * g = dynamic_cast<gsTHBSplineBasis<2,T> *>( basis ) )
1271 // marked.markTadmissible(m_m);
1272 // else if ( gsHBSplineBasis<2,T> * g = dynamic_cast< gsHBSplineBasis<2,T> *>( basis ) )
1273 // marked.markHadmissible(m_m);
1274 // else // if basis type unknown
1275 // marked.markHadmissible(m_m);
1276 // }
1277 // else if (m_options.getInt("Admissibility")==1)
1278 // marked.markTadmissible(m_m);
1279 // else if (m_options.getInt("Admissibility")==2)
1280 // marked.markHadmissible(m_m);
1281 // else // if basis type unknown
1282 // GISMO_ERROR("Admissibility type unknown or basis type not recognized");
1283
1284 // HBoxContainer container(marked.toUnitBoxes());
1285 // gsDebugVar(container);
1286 // std::vector<index_t> boxes = container.toRefBoxes();
1287 // if ((mp = dynamic_cast<gsMultiPatch<T>*>(m_input)))
1288 // mp->patch(pn).refineElements(container.toRefBoxes());
1289 // else if ((mb = dynamic_cast<gsMultiBasis<T>*>(m_input)))
1290 // mb->basis(pn).refineElements(container.toRefBoxes());
1291 // else
1292 // GISMO_ERROR("No gsMultiPatch or gsMultiBasis found; don't know what to refine");
1293 // }
1294 // else
1295 // {
1296 gsHBoxContainer<2,T> container = markedRef.patch(pn);
1297 container.toUnitBoxes();
1298 if (refExtension==0)
1299 if (nullptr != (mp = dynamic_cast<gsMultiPatch<T>*>(m_input)))
1300 mp->patch(pn).refineElements( container.toRefBoxes(pn) );
1301 else if (nullptr != (mb = dynamic_cast<gsMultiBasis<T>*>(m_input)))
1302 mb->basis( pn).refineElements( container.toRefBoxes(pn) );
1303 else
1304 GISMO_ERROR("No gsMultiPatch or gsMultiBasis found");
1305 else
1306 if (nullptr != (mp = dynamic_cast<gsMultiPatch<T>*>(m_input)))
1307 {
1308 // Refine all of the found refBoxes in this patch
1309 std::vector<index_t> elements = mp->patch(pn).basis().asElements(container.toCoords(pn), refExtension);
1310 mp->patch(pn).refineElements( elements );
1311 }
1312 else if (nullptr != (mb = dynamic_cast<gsMultiBasis<T>*>(m_input)))
1313 {
1314 // Refine all of the found refBoxes in this patch
1315 mb->basis( pn).refine(container.toCoords(pn), refExtension );
1316 }
1317 else
1318 GISMO_ERROR("No gsMultiPatch or gsMultiBasis found");
1319 // }
1320
1321 }
1322}
1323
1324template<class T>
1325void gsAdaptiveMeshing<T>::_unrefineMarkedElements( const HBoxContainer & markedCrs,
1326 index_t crsExtension)
1327{
1328 gsBasis<T> * basis = nullptr;
1329
1330 gsMultiPatch<T> * mp;
1331 gsMultiBasis<T> * mb;
1332
1333 for (index_t pn=0; pn < m_input->nPieces(); ++pn )// for all patches
1334 {
1335 if (nullptr != (mp = dynamic_cast<gsMultiPatch<T>*>(m_input)) ) basis = &(mp->basis(pn));
1336 if (nullptr != (mb = dynamic_cast<gsMultiBasis<T>*>(m_input)) ) basis = &(mb->basis(pn));
1337 GISMO_ENSURE(basis!=nullptr,"Object is not gsMultiBasis or gsMultiPatch");
1338
1339 // if (m_options.getSwitch("Admissible"))
1340 // {
1341 // if (m_options.getInt("Admissibility")==0)
1342 // {
1343 // if ( gsTHBSplineBasis<2,T> * g = dynamic_cast<gsTHBSplineBasis<2,T> *>( basis ) )
1344 // marked.markTadmissible(m_m);
1345 // else if ( gsHBSplineBasis<2,T> * g = dynamic_cast< gsHBSplineBasis<2,T> *>( basis ) )
1346 // marked.markHadmissible(m_m);
1347 // else // if basis type unknown
1348 // marked.markHadmissible(m_m);
1349 // }
1350 // else if (m_options.getInt("Admissibility")==1)
1351 // marked.markTadmissible(m_m);
1352 // else if (m_options.getInt("Admissibility")==2)
1353 // marked.markHadmissible(m_m);
1354 // else // if basis type unknown
1355 // GISMO_ERROR("Admissibility type unknown or basis type not recognized");
1356
1357 // HBoxContainer container(marked.toUnitBoxes());
1358 // if ((mp = dynamic_cast<gsMultiPatch<T>*>(m_input)))
1359 // mp->patch(pn).unrefineElements(container.toCrsBoxes());
1360 // else if ((mb = dynamic_cast<gsMultiBasis<T>*>(m_input)))
1361 // mb->basis(pn).unrefineElements(container.toCrsBoxes());
1362 // else
1363 // GISMO_ERROR("No gsMultiPatch or gsMultiBasis found; don't know what to refine");
1364 // }
1365 // else
1366 // {
1367 gsHBoxContainer<2,T> container = markedCrs.patch(pn);
1368 container.toUnitBoxes();
1369 if (crsExtension==0)
1370 if (nullptr != (mp = dynamic_cast<gsMultiPatch<T>*>(m_input)))
1371 mp->patch(pn).unrefineElements( container.toCrsBoxes(pn) );
1372 else if (nullptr != (mb = dynamic_cast<gsMultiBasis<T>*>(m_input)))
1373 mb->basis( pn).unrefineElements( container.toCrsBoxes(pn) );
1374 else
1375 GISMO_ERROR("No gsMultiPatch or gsMultiBasis found");
1376 else
1377 if ((mp = dynamic_cast<gsMultiPatch<T>*>(m_input))!= nullptr)
1378 {
1379 // Unrefine all of the found refBoxes in this patch
1380 std::vector<index_t> elements = mp->patch(pn).basis().asElementsUnrefine(markedCrs.toCoords(pn), crsExtension);
1381 mp->patch(pn).unrefineElements( elements );
1382 }
1383 else if ((mb = dynamic_cast<gsMultiBasis<T>*>(m_input))!= nullptr)
1384 {
1385 // Unrefine all of the found refBoxes in this patch
1386 mb->unrefine( pn, markedCrs.toCoords(pn), crsExtension );
1387 }
1388 else
1389 GISMO_ERROR("No gsMultiPatch or gsMultiBasis found");
1390 // }
1391
1392 }
1393}
1394
1395template<class T>
1396typename gsAdaptiveMeshing<T>::HBoxContainer gsAdaptiveMeshing<T>::_toContainer( const std::vector<bool> & bools) const
1397{
1398 HBoxContainer container;
1399
1400 #pragma omp parallel
1401 {
1402#ifdef _OPENMP
1403 const int tid = omp_get_thread_num();
1404 const int nt = omp_get_num_threads();
1405 index_t patch_cnt = 0;
1406#endif
1407
1408 index_t c = 0;
1409 gsBasis<T> * basis = nullptr;
1410
1411 gsMultiPatch<T> * mp;
1412 gsMultiBasis<T> * mb;
1413 typename gsBasis<T>::domainIter domIt;
1414 gsHDomainIterator<T,2> * domHIt = nullptr;
1415 for (index_t patchInd=0; patchInd < m_input->nPieces(); ++patchInd)
1416 {
1417 // Initialize domain element iterator
1418 if ((mp = dynamic_cast<gsMultiPatch<T>*>(m_input))!= nullptr ) basis = &(mp->basis(patchInd));
1419 if ((mb = dynamic_cast<gsMultiBasis<T>*>(m_input))!= nullptr ) basis = &(mb->basis(patchInd));
1420 GISMO_ASSERT(basis!=nullptr,"Object is not gsMultiBasis or gsMultiPatch");
1421 // for all elements in patch pn
1422 domIt = basis->makeDomainIterator();
1423 domHIt = dynamic_cast<gsHDomainIterator<T,2> *>(domIt.get());
1424 GISMO_ENSURE(domHIt!=nullptr,"Domain should be hierarchical");
1425
1426#ifdef _OPENMP
1427 c = patch_cnt + tid;
1428 patch_cnt += domHIt->numElements();// a bit costly
1429 for ( domHIt->next(tid); domHIt->good(); domHIt->next(nt) )
1430#else
1431 for (; domHIt->good(); domHIt->next() )
1432#endif
1433 {
1434 if (bools[c])
1435#pragma omp critical (gsAdaptiveMeshingtoContainer)
1436 container.add(HBox(domHIt,patchInd));
1437
1438# ifdef _OPENMP
1439 c += nt;
1440# else
1441 c++;
1442# endif
1443 }
1444 }
1445 }
1446
1447 return container;
1448}
1449
1450template<class T>
1451index_t gsAdaptiveMeshing<T>::numBlocked() const
1452{
1453 gsMaxLvlCompare<2,T> comp(m_maxLvl);
1454 index_t numBlocked = 0;
1455 for (typename boxMapType::const_iterator it=m_boxes.cbegin(); it!=m_boxes.cend(); it++)
1456 numBlocked += comp.check(*it->second);
1457
1458 return numBlocked;
1459}
1460
1461template<class T>
1462index_t gsAdaptiveMeshing<T>::numElements() const
1463{
1464 return m_indices.size();
1465}
1466
1467template<class T>
1468void gsAdaptiveMeshing<T>::assignErrors(const std::vector<T> & elError)
1469{
1470 this->_assignErrors(m_boxes,elError);
1471}
1472
1473template<class T>
1474T gsAdaptiveMeshing<T>::blockedError() const
1475{
1476 gsMaxLvlCompare<2,T> comp(m_maxLvl);
1477 T error = 0;
1478 for (typename boxMapType::const_iterator it=m_boxes.cbegin(); it!=m_boxes.cend(); it++)
1479 {
1480 if (!(comp.check(*it->second)))
1481 error += it->second->error();
1482 }
1483
1484 return error;
1485}
1486
1487template<class T>
1488T gsAdaptiveMeshing<T>::nonBlockedError() const
1489{
1490 gsMaxLvlCompare<2,T> comp(m_maxLvl);
1491 T error = 0;
1492 for (typename boxMapType::const_iterator it=m_boxes.cbegin(); it!=m_boxes.cend(); it++)
1493 {
1494 if (comp.check(*it->second))
1495 error += it->second->error();
1496 }
1497
1498 return error;
1499}
1500
1501// The functions below are deprecated
1502
1510// template<class T>
1511// void gsAdaptiveMeshing<T>::_flattenElementsToLevel(const index_t level)
1512// {
1513// // Get all the elements of which the level exceeds level
1514// std::vector<bool> elMarked;
1515// _markLevelThreshold(level,elMarked);
1516
1517// const index_t dim = m_input->domainDim();
1518
1519// // numMarked: Number of marked cells on current patch, also currently marked cell
1520// // poffset : offset index for the first element on a patch
1521// // globalCount: counter for the current global element index
1522// index_t numMarked, poffset = 0, globalCount = 0;
1523
1524// // crsBoxes: contains marked boxes on a given patch
1525// gsMatrix<T> crsBoxes;
1526
1527// gsBasis<T> * basis = nullptr;
1528
1529// gsMultiPatch<T> * mp;
1530// gsMultiBasis<T> * mb;
1531// for (index_t pn=0; pn < m_input->nPieces(); ++pn )// for all patches
1532// {
1533// if ( (mp = dynamic_cast<gsMultiPatch<T>*>(m_input)) ) basis = &(mp->basis(pn));
1534// if ( (mb = dynamic_cast<gsMultiBasis<T>*>(m_input)) ) basis = &(mb->basis(pn));
1535// GISMO_ASSERT(basis!=nullptr,"Object is not gsMultiBasis or gsMultiPatch");
1536// // Get number of elements to be refined on this patch
1537// const index_t numEl = basis->numElements();
1538// numMarked = std::count_if(elMarked.begin() + poffset,
1539// elMarked.begin() + poffset + numEl,
1540// GS_BIND2ND(std::equal_to<bool>(), true) );
1541
1542// poffset += numEl;
1543// crsBoxes.resize(dim, 2*numMarked);
1544// numMarked = 0;// counting current patch element to be refined
1545
1546// // for all elements in patch pn
1547// typename gsBasis<T>::domainIter domIt = basis->makeDomainIterator();
1548// for (; domIt->good(); domIt->next())
1549// {
1550// if( elMarked[ globalCount++ ] ) // refine this element ?
1551// {
1552// // Construct degenerate box by setting both
1553// // corners equal to the center
1554// crsBoxes.col(2*numMarked ) =
1555// crsBoxes.col(2*numMarked+1) = domIt->centerPoint();
1556
1557// // Advance marked cells counter
1558// numMarked++;
1559// }
1560// }
1561
1562// std::vector<index_t> elements;
1563// elements = basis->asElementsUnrefine(crsBoxes,0);
1564// const index_t offset = 2*dim+1;
1565// index_t diff = 0;
1566// GISMO_ASSERT(elements.size()%offset==0,"Element boxes have wrong size. Boxes should have size "<<offset<<" per box");
1567// for (size_t k = 0; k!=elements.size()/offset; k++)
1568// {
1569// // gsDebug<<"Old Box:\n";
1570// // for (index_t kk = 0; kk!=2*dim+1; kk++)
1571// // {
1572// // gsDebug<<elements[offset*k+kk]<<",";
1573// // }
1574// // gsDebug<<"\n";
1575
1576// if ((diff = elements[offset*k] - level) > 0)
1577// {
1578// elements[offset*k] = level;
1579// for (index_t kk = 0; kk!=2*dim; kk++)
1580// {
1581
1582// elements[offset*k+kk+1] = elements[offset*k+kk+1] >> diff;
1583// // Remove comment below
1584// // // yes, exactly: if i is index in old level a and needs to become a level b box, b<a, then the new index j is (in c++ computation):
1585// // j = i >> a-b;
1586// // // this is division by 2^{a-b} using bit-operations
1587// }
1588// }
1589// // gsDebug<<"New Box:\n";
1590// // for (index_t kk = 0; kk!=2*dim+1; kk++)
1591// // {
1592// // gsDebug<<elements[offset*k+kk]<<",";
1593// // }
1594// // gsDebug<<"\n";
1595// }
1596
1597// if ((mp = dynamic_cast<gsMultiPatch<T>*>(m_input)))
1598// {
1599// mp->patch(pn).unrefineElements( elements );
1600// }
1601// else if ((mb = dynamic_cast<gsMultiBasis<T>*>(m_input)))
1602// {
1603// // Refine all of the found refBoxes in this patch
1604// mb->unrefineElements(pn, elements );
1605// }
1606// else
1607// GISMO_ERROR("No gsMultiPatch or gsMultiBasis found");
1608// }
1609// }
1610
1620// template<class T>
1621// void gsAdaptiveMeshing<T>::_unrefineElementsThreshold(
1622// const index_t level)
1623// {
1624// // Get all the elements of which the level exceeds level
1625// std::vector<bool> elMarked;
1626// _markLevelThreshold(level,elMarked);
1627
1628// const index_t dim = m_input->domainDim();
1629
1630// // numMarked: Number of marked cells on current patch, also currently marked cell
1631// // poffset : offset index for the first element on a patch
1632// // globalCount: counter for the current global element index
1633// index_t numMarked, poffset = 0, globalCount = 0;
1634
1635// // crsBoxes: contains marked boxes on a given patch
1636// gsMatrix<T> crsBoxes;
1637
1638// gsBasis<T> * basis = nullptr;
1639
1640// gsMultiPatch<T> * mp;
1641// gsMultiBasis<T> * mb;
1642// for (index_t pn=0; pn < m_input->nPieces(); ++pn )// for all patches
1643// {
1644// if ( (mp = dynamic_cast<gsMultiPatch<T>*>(m_input)) ) basis = &(mp->basis(pn));
1645// if ( (mb = dynamic_cast<gsMultiBasis<T>*>(m_input)) ) basis = &(mb->basis(pn));
1646// GISMO_ASSERT(basis!=nullptr,"Object is not gsMultiBasis or gsMultiPatch");
1647// // Get number of elements to be refined on this patch
1648// const index_t numEl = basis->numElements();
1649// numMarked = std::count_if(elMarked.begin() + poffset,
1650// elMarked.begin() + poffset + numEl,
1651// GS_BIND2ND(std::equal_to<bool>(), true) );
1652
1653// poffset += numEl;
1654// crsBoxes.resize(dim, 2*numMarked);
1655// numMarked = 0;// counting current patch element to be refined
1656
1657// // for all elements in patch pn
1658// typename gsBasis<T>::domainIter domIt = basis->makeDomainIterator();
1659// for (; domIt->good(); domIt->next())
1660// {
1661// if( elMarked[ globalCount++ ] ) // refine this element ?
1662// {
1663// // Construct degenerate box by setting both
1664// // corners equal to the center
1665// crsBoxes.col(2*numMarked ) =
1666// crsBoxes.col(2*numMarked+1) = domIt->centerPoint();
1667
1668// // Advance marked cells counter
1669// numMarked++;
1670// }
1671// }
1672
1673// std::vector<index_t> elements;
1674// elements = basis->asElementsUnrefine(crsBoxes,0);
1675// const index_t offset = 2*dim+1;
1676// index_t diff = 0;
1677// GISMO_ASSERT(elements.size()%offset==0,"Element boxes have wrong size. Boxes should have size "<<offset<<" per box");
1678// for (size_t k = 0; k!=elements.size()/offset; k++)
1679// {
1680// // gsDebug<<"Old Box:\n";
1681// // for (index_t kk = 0; kk!=2*dim+1; kk++)
1682// // {
1683// // gsDebug<<elements[offset*k+kk]<<",";
1684// // }
1685// // gsDebug<<"\n";
1686
1687// if ((diff = elements[offset*k] - level) > 0)
1688// {
1689// elements[offset*k] = level;
1690// for (index_t kk = 0; kk!=2*dim; kk++)
1691// {
1692
1693// elements[offset*k+kk+1] = elements[offset*k+kk+1] >> 1;
1694
1695// // // yes, exactly: if i is index in old level a and needs to become a level b box, b<a, then the new index j is (in c++ computation):
1696// // j = i >> a-b;
1697// // // this is division by 2^{a-b} using bit-operations
1698// }
1699// }
1700// // gsDebug<<"New Box:\n";
1701// // for (index_t kk = 0; kk!=2*dim+1; kk++)
1702// // {
1703// // gsDebug<<elements[offset*k+kk]<<",";
1704// // }
1705// // gsDebug<<"\n";
1706// }
1707
1708// if ((mp = dynamic_cast<gsMultiPatch<T>*>(m_input)))
1709// {
1710// mp->patch(pn).unrefineElements( elements );
1711// }
1712// else if ((mb = dynamic_cast<gsMultiBasis<T>*>(m_input)))
1713// {
1714// // Refine all of the found refBoxes in this patch
1715// mb->unrefineElements(pn, elements );
1716// }
1717// else
1718// GISMO_ERROR("No gsMultiPatch or gsMultiBasis found");
1719// }
1720// }
1721
1722// template <class T>
1723// void gsAdaptiveMeshing<T>::_getElLevels( std::vector<index_t> & elLevels)
1724// {
1725// // Now just check for each element, whether the level
1726// // is above the target level or not, and mark accordingly.
1727// //
1728// index_t Nelements = 0;
1729// if ( gsMultiPatch<T> * mp = dynamic_cast<gsMultiPatch<T>*>(m_input) ) Nelements = gsMultiBasis<T>(*mp).totalElements();
1730// if ( gsMultiBasis<T> * mb = dynamic_cast<gsMultiBasis<T>*>(m_input) ) Nelements = mb->totalElements();
1731// GISMO_ASSERT(Nelements!=0,"Number of elements is zero? It might be that the m_input is not a gsMultiBasis or gsMultiPatch");
1732
1733// elLevels.resize(Nelements);
1734
1735// gsBasis<T> * basis = nullptr;
1736
1737// #pragma omp parallel
1738// {
1739// #ifdef _OPENMP
1740// const int tid = omp_get_thread_num();
1741// const int nt = omp_get_num_threads();
1742// index_t patch_cnt = 0;
1743// #endif
1744
1745// index_t c = 0;
1746// for (index_t patchInd=0; patchInd < m_input->nPieces(); ++patchInd)
1747// {
1748// // Initialize domain element iterator
1749// if ( gsMultiPatch<T> * mp = dynamic_cast<gsMultiPatch<T>*>(m_input) ) basis = &(mp->basis(patchInd));
1750// if ( gsMultiBasis<T> * mb = dynamic_cast<gsMultiBasis<T>*>(m_input) ) basis = &(mb->basis(patchInd));
1751// GISMO_ASSERT(basis!=nullptr,"Object is not gsMultiBasis or gsMultiPatch");
1752// // for all elements in patch pn
1753// typename gsBasis<T>::domainIter domIt = basis->makeDomainIterator();
1754// gsHDomainIterator<T,2> * domHIt = nullptr;
1755// domHIt = dynamic_cast<gsHDomainIterator<T,2> *>(domIt.get());
1756// GISMO_ENSURE(domHIt!=nullptr,"Domain should be 2 dimensional for flattening");
1757
1758// #ifdef _OPENMP
1759// c = patch_cnt + tid;
1760// patch_cnt += domHIt->numElements();// a bit costy
1761// for ( domHIt->next(tid); domHIt->good(); domHIt->next(nt) )
1762// #else
1763// for (; domHIt->good(); domHIt->next() )
1764// #endif
1765// {
1766// # ifdef _OPENMP
1767// elLevels[c] = domHIt->getLevel();
1768// c += nt;
1769// # else
1770// elLevels[c++] = domHIt->getLevel();
1771// # endif
1772// }
1773// }
1774// }
1775// }
1776
1777// template <class T>
1778// void gsAdaptiveMeshing<T>::_markLevelThreshold( index_t level, HBoxContainer & elMarked)
1779// {
1780// // Now just check for each element, whether the level
1781// // is above the target level or not, and mark accordingly.
1782// //
1783// index_t Nelements = 0;
1784// if ( gsMultiPatch<T> * mp = dynamic_cast<gsMultiPatch<T>*>(m_input) ) Nelements = gsMultiBasis<T>(*mp).totalElements();
1785// if ( gsMultiBasis<T> * mb = dynamic_cast<gsMultiBasis<T>*>(m_input) ) Nelements = mb->totalElements();
1786// GISMO_ASSERT(Nelements!=0,"Number of elements is zero? It might be that the m_input is not a gsMultiBasis or gsMultiPatch");
1787
1788// elMarked.resize(Nelements);
1789
1790// gsBasis<T> * basis = nullptr;
1791
1792// #pragma omp parallel
1793// {
1794// #ifdef _OPENMP
1795// const int tid = omp_get_thread_num();
1796// const int nt = omp_get_num_threads();
1797// index_t patch_cnt = 0;
1798// #endif
1799
1800// index_t c = 0;
1801// for (index_t patchInd=0; patchInd < m_input->nPieces(); ++patchInd)
1802// {
1803// // Initialize domain element iterator
1804// if ( gsMultiPatch<T> * mp = dynamic_cast<gsMultiPatch<T>*>(m_input) ) basis = &(mp->basis(patchInd));
1805// if ( gsMultiBasis<T> * mb = dynamic_cast<gsMultiBasis<T>*>(m_input) ) basis = &(mb->basis(patchInd));
1806// GISMO_ASSERT(basis!=nullptr,"Object is not gsMultiBasis or gsMultiPatch");
1807// // for all elements in patch pn
1808// typename gsBasis<T>::domainIter domIt = basis->makeDomainIterator();
1809// gsHDomainIterator<T,2> * domHIt = nullptr;
1810// domHIt = dynamic_cast<gsHDomainIterator<T,2> *>(domIt.get());
1811// GISMO_ENSURE(domHIt!=nullptr,"Domain should be 2 dimensional for flattening");
1812
1813// #ifdef _OPENMP
1814// c = patch_cnt + tid;
1815// patch_cnt += domHIt->numElements();// a bit costy
1816// for ( domIt->next(tid); domIt->good(); domIt->next(nt) )
1817// #else
1818// for (; domIt->good(); domIt->next() )
1819// #endif
1820// {
1821// # ifdef _OPENMP
1822// elMarked[c] = (domHIt->getLevel() > level);
1823// c += nt;
1824// # else
1825// elMarked[c++] = (domHIt->getLevel() > level);
1826// # endif
1827// }
1828// }
1829// }
1830// }
1831
1832
1833} // namespace gismo
gismo::gsAdaptiveMeshing
Provides adaptive meshing routines.
Definition gsAdaptiveMeshing.h:42
gismo::gsHBoxCheck
Base class for performing checks on gsHBox objects.
Definition gsAdaptiveMeshingCompare.h:32
gismo::gsMinLvlCompare
Checks if the level of a gsHBox is bigger than a minimum level.
Definition gsAdaptiveMeshingCompare.h:47
gismo::gsOverlapCompare
Checks if the coarsening neighborhood of a box is empty and if it overlaps with a refinement mask.
Definition gsAdaptiveMeshingCompare.h:135
index_t
#define index_t
Definition gsConfig.h:32
gsDebug
#define gsDebug
Definition gsDebug.h:61
GISMO_NO_IMPLEMENTATION
#define GISMO_NO_IMPLEMENTATION
Definition gsDebug.h:129
GISMO_ERROR
#define GISMO_ERROR(message)
Definition gsDebug.h:118
GISMO_UNUSED
#define GISMO_UNUSED(x)
Definition gsDebug.h:112
GISMO_ENSURE
#define GISMO_ENSURE(cond, message)
Definition gsDebug.h:102
GISMO_ASSERT
#define GISMO_ASSERT(cond, message)
Definition gsDebug.h:89
gsHBSplineBasis.h
Provides declaration of HBSplineBasis class.
gsTHBSplineBasis.h
Provides declaration of THBSplineBasis class.
gismo
The G+Smo namespace, containing all definitions for the library.