G+Smo  24.08.0
Geometry + Simulation Modules
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 
20 namespace gismo
21 {
22 
23 template <class T>
24 gsAdaptiveMeshing<T>::gsAdaptiveMeshing()
25 {
26  defaultOptions();
27 }
28 
29 template <class T>
30 gsAdaptiveMeshing<T>::gsAdaptiveMeshing(gsFunctionSet<T> & input)
31 :
32 m_input(&input)
33 {
34  defaultOptions();
35  rebuild();
36 }
37 
38 template <class T>
39 void 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 
57 template <class T>
58 void 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 
122 template <class T>
123 void 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 
199 template <class T>
200 template<bool _coarsen,bool _admissible>
201 void 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 
224 template <class T>
225 void 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 
232 template <class T>
233 void 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 
239 template <class T>
240 void gsAdaptiveMeshing<T>::_refPredicates_into(std::vector<gsHBoxCheck<2,T> *> & predicates)
241 {
242  predicates.push_back(new gsMaxLvlCompare<2,T>(m_maxLvl));
243 }
244 
245 
246 template <class T>
247 std::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 
258 template <class T>
259 std::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 
270 template <class T>
271 std::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 
293 template <class T>
294 typename 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 
302 template <class T>
303 bool 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 
311 template <class T>
312 bool 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 
324 template <class T>
325 void gsAdaptiveMeshing<T>::_addAndMark( HBox & box, typename gsAdaptiveMeshing<T>::HBoxContainer & elMarked) const
326 {
327  _boxPtr(box)->mark();
328  elMarked.add(*_boxPtr(box));
329 }
330 
331 template <class T>
332 void 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 
338 template <class T>
339 void 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 
350 template <class T>
351 T 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 
365 template <class T>
366 template<bool _coarsen,bool _admissible>
367 typename std::enable_if< _coarsen && _admissible, void>::type
368 gsAdaptiveMeshing<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 
414 template <class T>
415 template<bool _coarsen,bool _admissible>
416 typename std::enable_if< _coarsen && !_admissible, void>::type
417 gsAdaptiveMeshing<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 
444 template <class T>
445 template<bool _coarsen,bool _admissible>
446 typename std::enable_if<!_coarsen && _admissible, void>::type
447 gsAdaptiveMeshing<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 
493 template <class T>
494 template<bool _coarsen,bool _admissible>
495 typename std::enable_if<!_coarsen && !_admissible, void>::type
496 gsAdaptiveMeshing<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 
525 template <class T>
526 template<bool _coarsen,bool _admissible>
527 typename std::enable_if< _coarsen && _admissible, void>::type
528 gsAdaptiveMeshing<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 
582 template <class T>
583 template<bool _coarsen,bool _admissible>
584 typename 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
586 gsAdaptiveMeshing<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 
616 template <class T>
617 template<bool _coarsen,bool _admissible>
618 typename std::enable_if<!_coarsen && _admissible, void>::type
619 gsAdaptiveMeshing<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 
669 template <class T>
670 template<bool _coarsen,bool _admissible>
671 typename std::enable_if<!_coarsen && !_admissible, void>::type
672 gsAdaptiveMeshing<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 
708 template <class T>
709 template<bool _coarsen,bool _admissible>
710 typename std::enable_if< _coarsen && _admissible, void>::type
711 gsAdaptiveMeshing<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 
752 template <class T>
753 template<bool _coarsen,bool _admissible>
754 typename std::enable_if< _coarsen && !_admissible, void>::type
755 gsAdaptiveMeshing<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 
783 template <class T>
784 template<bool _coarsen,bool _admissible>
785 typename std::enable_if<!_coarsen && _admissible, void>::type
786 gsAdaptiveMeshing<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 
821 template <class T>
822 template<bool _coarsen,bool _admissible>
823 typename std::enable_if<!_coarsen && !_admissible, void>::type
824 gsAdaptiveMeshing<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 
856 template <class T>
857 T 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 
870 template <class T>
871 template<bool _coarsen,bool _admissible>
872 typename std::enable_if< _coarsen && _admissible, void>::type
873 gsAdaptiveMeshing<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 
918 template <class T>
919 template<bool _coarsen,bool _admissible>
920 typename std::enable_if< _coarsen && !_admissible, void>::type
921 gsAdaptiveMeshing<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 
952 template <class T>
953 template<bool _coarsen,bool _admissible>
954 typename std::enable_if<!_coarsen && _admissible, void>::type
955 gsAdaptiveMeshing<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 
994 template <class T>
995 template<bool _coarsen,bool _admissible>
996 typename std::enable_if<!_coarsen && !_admissible, void>::type
997 gsAdaptiveMeshing<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 
1028 template<class T>
1029 void 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 
1051 template<class T>
1052 void 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 
1109 template<class T>
1110 void 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 
1119 template<class T>
1120 void 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 
1149 template<class T>
1150 void 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 
1178 template<class T>
1179 void gsAdaptiveMeshing<T>::markCrs_into(const std::vector<T> & elError, HBoxContainer & elMarked)
1180 {
1181  HBoxContainer container;
1182  this->markCrs_into(elError,container,elMarked);
1183 }
1184 
1185 template<class T>
1186 void gsAdaptiveMeshing<T>::markRef(const std::vector<T> & errors)
1187 {
1188  markRef_into( errors, m_markedRef);//,flag [coarse]);
1189 }
1190 
1191 template<class T>
1192 void gsAdaptiveMeshing<T>::markCrs(const std::vector<T> & errors)
1193 {
1194  markCrs_into( errors, m_markedRef, m_markedCrs);//,flag [coarse]);
1195 }
1196 
1197 template<class T>
1198 bool 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 
1206 template<class T>
1207 bool 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 
1215 template<class T>
1216 bool 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 
1227 template<class T>
1228 bool 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 
1251 template<class T>
1252 void 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 
1324 template<class T>
1325 void 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 
1395 template<class T>
1396 typename 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 
1450 template<class T>
1451 index_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 
1461 template<class T>
1462 index_t gsAdaptiveMeshing<T>::numElements() const
1463 {
1464  return m_indices.size();
1465 }
1466 
1467 template<class T>
1468 void gsAdaptiveMeshing<T>::assignErrors(const std::vector<T> & elError)
1469 {
1470  this->_assignErrors(m_boxes,elError);
1471 }
1472 
1473 template<class T>
1474 T 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 
1487 template<class T>
1488 T 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_NO_IMPLEMENTATION
#define GISMO_NO_IMPLEMENTATION
Definition: gsDebug.h:129
gsDebug
#define gsDebug
Definition: gsDebug.h:61
index_t
#define index_t
Definition: gsConfig.h:32
GISMO_ENSURE
#define GISMO_ENSURE(cond, message)
Definition: gsDebug.h:102
GISMO_ASSERT
#define GISMO_ASSERT(cond, message)
Definition: gsDebug.h:89
gsTHBSplineBasis.h
Provides declaration of THBSplineBasis class.
gismo::gsHBoxCheck
Base class for performing checks on gsHBox objects.
Definition: gsAdaptiveMeshingCompare.h:31
gismo::gsOverlapCompare
Checks if the coarsening neighborhood of a box is empty and if it overlaps with a refinement mask...
Definition: gsAdaptiveMeshingCompare.h:134
GISMO_UNUSED
#define GISMO_UNUSED(x)
Definition: gsDebug.h:112
GISMO_ERROR
#define GISMO_ERROR(message)
Definition: gsDebug.h:118
gismo::gsAdaptiveMeshing::_markElements
void _markElements(const std::vector< T > &elError, const index_t refCriterion, const std::vector< gsHBoxCheck< 2, T > * > &predicates, HBoxContainer &elMarked) const
Marks elements/cells for refinement.
Definition: gsAdaptiveMeshing.hpp:201
gismo::gsAdaptiveMeshing
Provides adaptive meshing routines.
Definition: gsAdaptiveMeshing.h:41
gsHBSplineBasis.h
Provides declaration of HBSplineBasis class.
gismo::gsMinLvlCompare
Checks if the level of a gsHBox is bigger than a minimum level.
Definition: gsAdaptiveMeshingCompare.h:46