G+Smo  25.01.0
Geometry + Simulation Modules
 
Loading...
Searching...
No Matches
gsHBox.hpp
Go to the documentation of this file.
1
14#pragma once
15
16#include <gsUtils/gsCombinatorics.h>
17#include <gsHSplines/gsHBSplineBasis.h>
18#include <gsHSplines/gsTHBSplineBasis.h>
19#include <gsHSplines/gsAABB.h>
20
21#include <gsIO/gsXml.h>
22
23namespace gismo
24{
25
26template <short_t d, class T>
27gsHBox<d, T>::gsHBox()
28:
29m_pid(-1),
30m_error(0),
31m_error_ref(0),
32m_error_crs(0),
33m_index(-1),
34m_marked(false)
35{
36 m_basis = nullptr;
37}
38
39template <short_t d, class T>
40gsHBox<d, T>::gsHBox(const gsHDomainIterator<T,d> * domHIt)
41:
42gsHBox(domHIt,-1)
43{ }
44
45template <short_t d, class T>
46gsHBox<d, T>::gsHBox(const gsHDomainIterator<T,d> * domHIt, const index_t pid)
47:
48m_pid(pid),
49m_error(0),
50m_error_ref(0),
51m_error_crs(0),
52m_index(-1),
53m_marked(false)
54{
55 m_basis = nullptr;
56 m_basis = static_cast<const gsHTensorBasis<d,T> *>(domHIt->m_basis);
57 GISMO_ASSERT(m_basis!=nullptr,"basis is not a gsHTensorBasis");
58
59 m_coords.resize(d,2);
60 m_coords.col(0) = domHIt->lowerCorner();
61 m_coords.col(1) = domHIt->upperCorner();
62 m_center = (m_coords.col(1) + m_coords.col(0))/2;
63 m_indices = _computeIndices(m_coords,m_center);
64}
65
66template <short_t d, class T>
67gsHBox<d, T>::gsHBox(const typename gsHBox<d,T>::point & low,const typename gsHBox<d,T>::point & upp, index_t level, const gsHTensorBasis<d,T> * basis, const index_t pid)
68:
69m_indices(low,upp,level),
70m_pid(pid),
71m_error(0),
72m_error_ref(0),
73m_error_crs(0),
74m_index(-1),
75m_marked(false)
76{
77 m_basis = basis;
78}
79
80template <short_t d, class T>
81gsHBox<d, T>::gsHBox(const gsAABB<d, index_t> & box, const gsHTensorBasis<d,T> * basis, const index_t pid)
82:
83m_indices(box),
84m_pid(pid),
85m_error(0),
86m_error_ref(0),
87m_error_crs(0),
88m_index(-1),
89m_marked(false)
90{
91 m_basis = basis;
92}
93
94template <short_t d, class T>
95gsHBox<d, T>::gsHBox(const std::vector<index_t> & indices, const gsHTensorBasis<d,T> * basis, const index_t pid)
96:
97m_pid(pid),
98m_error(0),
99m_error_ref(0),
100m_error_crs(0),
101m_index(-1),
102m_marked(false)
103{
104 GISMO_ENSURE(indices.size()==2*d+1,"Index size is wrong");
105 typename gsHBox<d,T>::point low, upp;
106 for (index_t k=0; k!=d; k++)
107 {
108 low[k] = indices[k+1];
109 upp[k] = indices[k+d+1];
110 }
111
112 m_indices = gsAABB<d, index_t>(low,upp,indices[0]);
113 // = gsAsVector<index_t,d>()
114 // upp;
115
116 m_basis = basis;
117}
118
119template <short_t d, class T>
120gsHBox<d, T>::gsHBox( const gsHBox<d,T> & other )
121{
122 operator=(other);
123}
124
125template <short_t d, class T>
126gsHBox<d, T>::gsHBox( gsHBox<d,T> && other )
127{
128 operator=(give(other));
129}
130
131template <short_t d, class T>
132gsHBox<d,T> & gsHBox<d, T>::operator= ( const gsHBox<d,T> & other )
133{
134 if (this!=&other)
135 {
136 m_indices = other.m_indices;
137 m_pid = other.m_pid;
138 m_coords = other.m_coords;
139 m_center = other.m_center;
140 m_basis = other.m_basis;
141 m_error = other.m_error;
142 m_error_ref = other.m_error_ref;
143 m_error_crs = other.m_error_crs;
144 m_marked = other.m_marked;
145 m_index = other.m_index;
146 }
147 return *this;
148}
149
150template <short_t d, class T>
151gsHBox<d,T> & gsHBox<d, T>::operator= ( gsHBox<d,T> && other )
152{
153 m_indices = give(other.m_indices);
154 m_pid = give(other.m_pid);
155 m_coords = give(other.m_coords);
156 m_center = give(other.m_center);
157 m_basis = give(other.m_basis);
158 m_error = give(other.m_error);
159 m_error_ref = give(other.m_error_ref);
160 m_error_crs = give(other.m_error_crs);
161 m_marked = give(other.m_marked);
162 m_index = give(other.m_index);
163 return *this;
164}
165
166template <short_t d, class T>
167bool gsHBox<d, T>::isContained(const gsHBox<d,T> & other) const
168{
169 // bool res = true;
170 // res &= this->level() == other.level();
171 // for (index_t i=0; i!=d && res; i++)
172 // {
173 // res &= this->lowerIndex().at(i) <= other.lowerIndex().at(i);
174 // res &= this->upperIndex().at(i) >= other.upperIndex().at(i);
175 // }
176 // return res;
177
178
179 return gsHBoxIsContained<d,T>()(*this,other);
180
181}
182
183template <short_t d, class T>
184bool gsHBox<d, T>::contains(const gsHBox<d,T> & other) const
185{
186 // bool res = true;
187 // res &= this->level() == other.level();
188 // for (index_t i=0; i!=d && res; i++)
189 // {
190 // res &= this->lowerIndex().at(i) >= other.lowerIndex().at(i);
191 // res &= this->upperIndex().at(i) <= other.upperIndex().at(i);
192 // }
193 // return res;
194
195 return gsHBoxContains<d,T>()(*this,other);
196}
197
198template <short_t d, class T>
199bool gsHBox<d, T>::isSame(const gsHBox<d,T> & other) const
200{
201 bool res = true;
202 res &= this->patch() == other.patch();
203 res &= this->level() == other.level();
204 for (index_t i=0; i!=d && res; i++)
205 {
206 res &= this->lowerIndex().at(i) == other.lowerIndex().at(i);
207 res &= this->upperIndex().at(i) == other.upperIndex().at(i);
208 }
209 return res;
210}
211
212template <short_t d, class T>
213bool gsHBox<d, T>::isActive() const
214{
215 return m_basis->getLevelAtPoint(this->getCenter()) == this->level();
216}
217
218template <short_t d, class T>
219bool gsHBox<d, T>::isActiveOrContained() const
220{
221 return m_basis->getLevelAtPoint(this->getCenter()) >= this->level();
222}
223
224template <short_t d, class T>
225index_t gsHBox<d, T>::levelInCenter() const
226{
227 return m_basis->getLevelAtPoint(this->getCenter());
228}
229
230template <short_t d, class T>
231const gsMatrix<T> & gsHBox<d, T>::getCoordinates() const
232{
233 GISMO_ASSERT(m_coords.size()!=0,"Coordinates have not been computed! Please call computeCoordinates()");
234 return m_coords;
235}
236
237template <short_t d, class T>
238const gsMatrix<T> & gsHBox<d, T>::getCenter() const
239{
240 // GISMO_ASSERT(m_center.size()!=0,"Center has not been computed! Please call computeCenter()");
241 if (m_center.size()==0)
242 this->computeCenter();
243 return m_center;
244}
245
246template <short_t d, class T>
247gsVector<T,d> gsHBox<d, T>::lowerCorner() const
248{
249 GISMO_ASSERT(m_coords.size()!=0,"Coordinates have not been computed! Please call computeCoordinates()");
250 return m_coords.col(0);
251}
252
253template <short_t d, class T>
254gsVector<T,d> gsHBox<d, T>::upperCorner() const
255{
256 GISMO_ASSERT(m_coords.size()!=0,"Coordinates have not been computed! Please call computeCoordinates()");
257 return m_coords.col(1);
258}
259
260template <short_t d, class T>
261T gsHBox<d, T>::getCellSize() const
262{
263 return (upperCorner() - lowerCorner()).norm();
264}
265
266template <short_t d, class T>
267T gsHBox<d, T>::getMinCellLength() const
268{
269 return (upperCorner() - lowerCorner()).minCoeff();
270}
271
272template <short_t d, class T>
273T gsHBox<d, T>::getMaxCellLength() const
274{
275 return (upperCorner() - lowerCorner()).maxCoeff();
276}
277
278template <short_t d, class T>
279const typename gsHBox<d,T>::point & gsHBox<d, T>::lowerIndex() const { return m_indices.first; }
280template <short_t d, class T>
281const typename gsHBox<d,T>::point & gsHBox<d, T>::upperIndex() const { return m_indices.second; }
282
283template <short_t d, class T>
284index_t gsHBox<d, T>::patch() const { return m_pid; }
285
286template <short_t d, class T>
287index_t gsHBox<d, T>::level() const { return m_indices.level; }
288
289template <short_t d, class T>
290void gsHBox<d, T>::setError(T error)
291{
292 m_error = error;
293}
294
295template <short_t d, class T>
296void gsHBox<d, T>::setAndProjectError(T error, index_t alpha, index_t beta)
297{
298 this->setError(error);
299 m_error_ref= math::pow(1./2.,alpha*m_basis->maxDegree()+beta) * error ;
300 m_error_crs= math::pow(2.,alpha*m_basis->maxDegree()+beta) * error ;
301}
302
303template <short_t d, class T>
304T gsHBox<d, T>::error() const { return m_error; }
305
306template <short_t d, class T>
307T gsHBox<d, T>::projectedErrorRef() const { return m_error_ref; }
308
309template <short_t d, class T>
310T gsHBox<d, T>::projectedErrorCrs() const { return m_error_crs; }
311
312template <short_t d, class T>
313T gsHBox<d, T>::projectedImprovement() const
314{
315 return m_error-m_error_ref;
316}
317
318template <short_t d, class T>
319T gsHBox<d, T>::projectedSetBack() const
320{
321 return m_error_crs-m_error;
322}
323
324template <short_t d, class T>
325void gsHBox<d, T>::setIndex(index_t index) { m_index = index; }
326
327template <short_t d, class T>
328index_t gsHBox<d, T>::index() const { return m_index; }
329
330template <short_t d, class T>
331void gsHBox<d, T>::mark() { m_marked = true; }
332
333template <short_t d, class T>
334void gsHBox<d, T>::unmark() { m_marked = false; }
335
336template <short_t d, class T>
337bool gsHBox<d, T>::marked() const { return m_marked; }
338
339template <short_t d, class T>
340void gsHBox<d, T>::setMark(bool mark) { m_marked = mark; }
341
342template <short_t d, class T>
343gsHBox<d,T> gsHBox<d, T>::getParent() const
344{
345 GISMO_ENSURE(this->level()>0,"Box is at ground level and has no parent");
346 gsAABB<d, index_t> box = _elevateBox(m_indices);
347 return gsHBox<d,T>(box,m_basis,m_pid);
348}
349
350template <short_t d, class T>
351gsHBox<d,T> gsHBox<d, T>::getAncestor(index_t k) const
352{
353 index_t lvl = this->level();
354 // GISMO_ASSERT(lvl > k,"Current level should be larger than requested level, l = "<<lvl<<", k = "<<k);
355 GISMO_ASSERT(k >= 0,"Requested ancestor k = "<<k<<" should be greater than or equal to 0");
356 GISMO_ASSERT(lvl >= 0,"Level lvl = "<<lvl<<" should be larger than 0");
357 // GISMO_ASSERT(k > 0,"Current level should be larger than requested level, l = "<<lvl<<", k = "<<k);
358 gsHBox<d,T> parent = this->getParent();
359 gsHBox<d,T> ancestor;
360 if (k < lvl - 1)
361 {
362 ancestor = parent.getAncestor(k);
363 return ancestor;
364 }
365 else if (k==lvl-1)
366 return parent;
367 else if (k==lvl)
368 return *this;
369 else
370 GISMO_ERROR("Cannot find ancestor with index k="<<k<<" on level l="<<lvl<<". Something went wrong?");
371}
372
373template <short_t d, class T>
374typename gsHBox<d,T>::Container gsHBox<d, T>::getChildren() const
375{
376 gsAABB<d, index_t> box = _lowerBox(m_indices);
377 gsHBox<d,T> childRegion(box,m_basis,m_pid);
378 return childRegion.toUnitBoxes();
379}
380
381template <short_t d, class T>
382typename gsHBox<d,T>::Container gsHBox<d, T>::getDescendants(index_t k) const
383{
384 index_t lvl = this->level();
385 // GISMO_ASSERT(lvl > k,"Current level should be larger than requested level, l = "<<lvl<<", k = "<<k);
386 GISMO_ASSERT(k >= 0,"Requested descendant k = "<<k<<" should be larger than or equal to 0");
387 GISMO_ASSERT(lvl >= 0,"Level lvl = "<<lvl<<" should be larger than 0");
388 Container descendants(std::pow(std::pow(2,d),k-lvl));
389 if (k==lvl)
390 descendants.push_back(*this);
391 if (k==lvl+1)
392 descendants = this->getChildren();
393 else
394 {
395 Container tmp = this->getChildren(); // children have level lvl+1
396 for (index_t k_tmp = lvl+1; k_tmp!=k; k_tmp++)
397 {
398 descendants.clear();
399 // descendants.reserve()
400 for (Iterator it=tmp.begin(); it!=tmp.end(); it++)
401 {
402 Container tmp2 = it->getChildren();
403 for (Iterator it2=tmp2.begin(); it2!=tmp2.end(); it2++)
404 descendants.push_back(*it2);
405 }
406 tmp = descendants;
407 }
408 }
409 return descendants;
410}
411
412template <short_t d, class T>
413typename gsHBox<d,T>::Container gsHBox<d, T>::getSiblings() const
414{
415 auto isthis = [this](const gsHBox<d,T> & other) { return gsHBoxEqual<d,T>()(other,*this); };
416 gsHBox<d,T> parent = this->getParent();
417 typename gsHBox<d,T>::Container children = parent.getChildren();
418 Iterator toErase = std::find_if(children.begin(),children.end(),isthis);
419 if (toErase!=children.end())
420 children.erase(toErase);
421 else
422 GISMO_ERROR("Something went wrong");
423 return children;
424}
425
426template <short_t d, class T>
427typename gsHBox<d,T>::Container gsHBox<d, T>::toContainer()
428{
429 Container container;
430 container.push_back(*this);
431 return container;
432}
433
434template <short_t d, class T>
435typename gsHBox<d,T>::HContainer gsHBox<d, T>::toHContainer()
436{
437 HContainer container;
438 container.resize(this->level()+1);
439 container.at(this->level()).push_back(*this);
440 return container;
441}
442
443template <short_t d, class T>
444typename gsHBox<d, T>::Container gsHBox<d, T>::getSupportExtension()
445{
446 this->computeCenter();
447 index_t lvl = this->level();
448 // Compute actives
449 gsMatrix<index_t> acts;
450 m_basis->tensorLevel(lvl).active_into(m_center,acts);
451
452 // Support extension
453 gsMatrix<T> cells(d,2*acts.rows());
454 gsMatrix<T> support;
455 gsHBox<d,T> supportBox;
456 gsAABB<d, index_t> aabb;
457 Container container;
458 Container tmpContainer;
459 for (index_t act = 0; act!=acts.rows(); act++)
460 {
461 support = m_basis->tensorLevel(lvl).support(acts(act,0));
462 aabb = _computeIndices(support,lvl);
463
464 supportBox = gsHBox<d,T>(aabb,m_basis,m_pid);
465
466 // Split the boxes into index interval with coordinate delta 1
467 tmpContainer = supportBox.toUnitBoxes();
468 for (cIterator it = tmpContainer.begin(); it!=tmpContainer.end(); it++)
469 container.push_back(*it);
470 }
471 // Remove duplicates
472 Container container2 = _makeUnique(container);
473 return container2;
474}
475
476template <short_t d, class T>
477typename gsHBox<d,T>::Container gsHBox<d, T>::getMultiLevelSupportExtension(index_t k)
478{
479 index_t lvl = this->level();
480 GISMO_ASSERT(k <= lvl,"Ancestor must be requested from a level lower than the current level (k <= l), but k="<<k<<" and lvl="<<lvl);
481 if (k==lvl)
482 {
483 return this->getSupportExtension();
484 }
485 else
486 {
487 gsHBox<d,T> ancestor = this->getAncestor(k);
488 return ancestor.getSupportExtension();
489 }
490}
491
492template <short_t d, class T>
493typename gsHBox<d,T>::Container gsHBox<d, T>::getHneighborhood(index_t m)
494{
495 Container extension;
496 Container neighborhood;
497
498 index_t lvl = this->level();
499 index_t k = lvl - m + 1;
500 if (k>=0)
501 {
502 // Get multi level support extension on level k
503 extension = this->getMultiLevelSupportExtension(k);
504
505 // Eliminate elements which are too low
506 for (Iterator it = extension.begin(); it!=extension.end(); it++)
507 {
508 it->computeCenter(); // needed to check active
509 if (it->isActive())
510 neighborhood.push_back(*it);
511 }
512 }
513 return neighborhood;
514}
515
516template <short_t d, class T>
517typename gsHBox<d,T>::Container gsHBox<d, T>::getTneighborhood(index_t m)
518{
519 // Everything is in the same level so we can use normal Container
520 Container neighborhood;
521 Container parents, extension;
522
523 index_t lvl = this->level();
524 index_t k = lvl - m + 2;
525 if (k-1>=0)
526 {
527 // Get multi level support extension on level k
528 extension = this->getMultiLevelSupportExtension(k);
529 // Eliminate elements which are too low
530 parents = _getParents(extension);
531
532 for (Iterator it = parents.begin(); it!=parents.end(); it++)
533 {
534 it->computeCenter(); // needed to check active
535 if (it->isActive())
536 neighborhood.push_back(*it);
537 }
538 // neighborhood = parents;
539 }
540 return neighborhood;
541}
542
543template <short_t d, class T>
544typename gsHBox<d,T>::Container gsHBox<d, T>::getNeighborhood(index_t m)
545{
546 Container result;
547 if (dynamic_cast<const gsTHBSplineBasis<d,T>*>(m_basis))
548 result = this->getNeighborhood<gsHNeighborhood::T>(m);
549 else if (dynamic_cast<const gsHBSplineBasis<d,T>*>(m_basis))
550 result = this->getNeighborhood<gsHNeighborhood::H>(m);
551 else
552 GISMO_ERROR("Basis type should be gsTHBSplineBasis or gsHBSplineBasis");
553 return result;
554}
555
556template <short_t d, class T>
557typename gsHBox<d,T>::Container gsHBox<d, T>::getCextension(index_t m)
558{
559 // See Definition 3.5 from [1]
560 // [1] Carraturo, M., Giannelli, C., Reali, A. & Vázquez, R.
561 // Suitably graded THB-spline refinement and coarsening: Towards an adaptive isogeometric analysis of additive manufacturing processes.
562 // Comput. Methods Appl. Mech. Eng. 348, 660–679 (2019).
563 //
564 // This function returns the coarsening neighborhood of \a this, which is an element to be reactivated
565 Container extension, children, childExtension, descendants, result;
566 index_t targetLvl = this->level() + m;
567
568 children = this->getChildren();
569 for (typename gsHBox<d,T>::Iterator child=children.begin(); child!=children.end(); child++)
570 {
571 childExtension = child->getSupportExtension();
572 extension = gsHBoxUtils<d,T>::Union(extension,childExtension);
573 }
574
575 // All elements in the support extension of the children (so at level l+1)
576 extension = gsHBoxUtils<d,T>::Union(extension,children);
577 // question: needed MINUS the siblings??
578 // extension = gsHBoxUtils<d,T>::Difference(extension,children);
579
580 // result = extension;
581 // gsDebugVar(targetLvl);
582 for (typename gsHBox<d,T>::Iterator eIt = extension.begin(); eIt!=extension.end(); eIt++)
583 {
584 descendants = eIt->getDescendants(targetLvl);
585 result.insert(result.end(), descendants.begin(), descendants.end());
586 }
587 return result;
588}
589
590template <short_t d, class T>
591typename gsHBox<d,T>::Container gsHBox<d, T>::getCneighborhood(index_t m)
592{
593 // See Definition 3.5 from [1]
594 // [1] Carraturo, M., Giannelli, C., Reali, A. & Vázquez, R.
595 // Suitably graded THB-spline refinement and coarsening: Towards an adaptive isogeometric analysis of additive manufacturing processes.
596 // Comput. Methods Appl. Mech. Eng. 348, 660–679 (2019).
597 //
598 // This function returns the coarsening neighborhood of \a this, which is an element to be reactivated
599 Container descendants, result;
600 descendants = this->getCextension(m);
601
602 for (typename gsHBox<d,T>::Iterator dIt = descendants.begin(); dIt!=descendants.end(); dIt++)
603 {
604 dIt->computeCenter(); // needed to check active
605 if (dIt->levelInCenter()>=dIt->level()) // the level is even larger (i.e. even higher decendant)!!
606 result.push_back(*dIt);
607 }
608 return result;
609}
610
611template <short_t d, class T>
612std::ostream& gsHBox<d, T>::print( std::ostream& os ) const
613{
614 os <<"Cell on patch "<<m_pid
615 <<" on level "<<m_indices.level<<". "
616 <<"\nIndices:\n"
617 <<"("<<m_indices.first.transpose()<<")"
618 <<" -- "
619 <<"("<<m_indices.second.transpose()<<")";
620 if (m_coords.cols()!=0)
621 {
622 os <<"\nKnot values:\n"
623 <<"("<<m_coords.col(0).transpose()<<")"
624 <<" -- "
625 <<"("<<m_coords.col(1).transpose()<<")";
626 }
627 os <<"\nmarked = "<<m_marked<<"";
628 os <<"\nerror = "<<m_error<<"";
629 return os;
630}
631
632template <short_t d, class T>
633void gsHBox<d, T>::computeCoordinates() const
634{
635 m_coords.resize(d,2);
636 gsVector<T> low(d), upp(d);
637
638 for (index_t i=0; i!=d; i++)
639 {
640 const gsKnotVector<T> & kv = m_basis->tensorLevel(m_indices.level).knots(i);
641 low[i] = kv.uValue(m_indices.first.at(i));
642 upp[i] = kv.uValue(m_indices.second.at(i));
643 }
644 m_coords.col(0) = low;
645 m_coords.col(1) = upp;
646}
647
648template <short_t d, class T>
649void gsHBox<d, T>::computeCenter() const
650{
651 this->computeCoordinates();
652 m_center = (m_coords.col(0) + m_coords.col(1))/2;
653}
654
655template <short_t d, class T>
656void gsHBox<d, T>::_computeIndices()
657{
658 m_indices = _computeIndices(m_coords);
659}
660
661template <short_t d, class T>
662gsAABB<d, index_t> gsHBox<d, T>::_computeIndices(const gsMatrix<T> & coords, index_t level)
663{
664 typename gsHBox<d,T>::point low,upp;
665 for(index_t j = 0; j < d;j++)
666 {
667 // Convert the parameter coordinates to (unique) knot indices
668 const gsKnotVector<T> & kv = m_basis->tensorLevel(level).knots(j);
669 int k1 = (std::upper_bound(kv.domainUBegin(), kv.domainUEnd(),
670 coords(j,0) ) - 1).uIndex();
671 int k2 = (std::upper_bound(kv.domainUBegin(), kv.domainUEnd()+1,
672 coords(j,1) ) - 1).uIndex();
673
674 // Trivial cells trigger some refinement
675 if ( k1 == k2)
676 {
677 if (0!=k1) {--k1;}
678 ++k2;
679 }
680
681 // Store the data...
682 low.at(j) = k1;
683 upp.at(j) = k2;
684 }
685 return gsAABB<d, index_t>(low,upp,level);
686}
687
688template <short_t d, class T>
689gsAABB<d, index_t> gsHBox<d, T>::_computeIndices(const gsMatrix<T> & coords)
690{
691 gsMatrix<T> center = (coords.col(0) + coords.col(1))/2;
692 index_t level = m_basis->getLevelAtPoint(center);
693 return _computeIndices(coords,level);
694}
695
696template <short_t d, class T>
697gsAABB<d, index_t> gsHBox<d, T>::_computeIndices(const gsMatrix<T> & coords, const gsMatrix<T> & center)
698{
699 index_t level = m_basis->getLevelAtPoint(center);
700 return _computeIndices(coords,level);
701}
702
703template <short_t d, class T>
704gsAABB<d, index_t> gsHBox<d, T>::_elevateBox(const gsAABB<d, index_t> & box) const
705{
706 typename gsHBox<d,T>::point low,upp;
707 for (index_t i = 0; i!=d; i++)
708 {
709 low.at(i) = box.first .at(i) / 2;
710 upp.at(i) = box.second.at(i) / 2 + (index_t)(box.second.at(i) % 2 != 0);
711 }
712 return gsAABB<d, index_t>(low,upp,box.level-1);
713}
714
715template <short_t d, class T>
716gsAABB<d, index_t> gsHBox<d, T>::_lowerBox(const gsAABB<d, index_t> & box) const
717{
718 typename gsHBox<d,T>::point low,upp;
719 for (index_t i = 0; i!=d; i++)
720 {
721 low.at(i) = box.first .at(i) * 2;
722 upp.at(i) = box.second.at(i) * 2;
723 }
724 return gsAABB<d, index_t>(low,upp,box.level+1);
725}
726
727template <short_t d, class T>
728typename gsHBox<d,T>::HContainer gsHBox<d, T>::_getParents(HContainer & container) const
729{
730 HContainer result;
731 result.resize(container.size()-1);
732
733 // Handle level 0 separately (operation cannot be performed on level 0)
734 GISMO_ASSERT(container[0].size()==0,"Boxes at level 0 cannot have a parent. Did something go wrong? You can run check() to see if the boxes are allocated coorectly");
735
736 HIterator resIt = result.begin();
737 for (HIterator hit = std::next(container.begin()); hit!=container.end(); hit++, resIt++)
738 for (Iterator it=hit->begin(); it!=hit->end(); it++)
739 resIt->push_back(it->getParent());
740
741 return result;
742}
743
744template <short_t d, class T>
745typename gsHBox<d,T>::Container gsHBox<d, T>::_getParents(Container & container) const
746{
747 Container result;
748 for (Iterator it=container.begin(); it!=container.end(); it++)
749 result.push_back(it->getParent());
750
751 return result;
752}
753
754
755template <short_t d, class T>
756typename gsHBox<d,T>::Container gsHBox<d, T>::toUnitBoxes() const
757{
758 point low, upp, cur, curupp, ones;
759 ones = gsVector<index_t,d>::Ones();
760 cur = low = this->lowerIndex();
761 upp = this->upperIndex();
762 // index_t nboxes = (upp-low).prod();
763 Container result;
764
765 bool next = true;
766 while (next)
767 {
768 curupp = cur + ones;
769 result.push_back(gsHBox<d,T>(cur,curupp,this->level(),m_basis,m_pid));
770 next = nextLexicographic(cur,low,upp);
771 }
772 return result;
773}
774
775template <short_t d, class T>
776typename gsHBox<d,T>::RefBox gsHBox<d, T>::toBox() const
777{
778 std::vector<index_t> result(5);
779 result[0] = this->level();
780 result[1] = this->lowerIndex()[0];
781 result[2] = this->lowerIndex()[1];
782 result[3] = this->upperIndex()[0];
783 result[4] = this->upperIndex()[1];
784 return result;
785}
786
787template <short_t d, class T>
788typename gsHBox<d,T>::RefBox gsHBox<d, T>::toRefBox(index_t targetLevel) const
789{
790 GISMO_ASSERT(targetLevel > this->level(),"Target level must be larger than the current level, but "<<targetLevel<<"=<"<<this->level());
791 std::vector<index_t> result(2*d+1);
792 index_t diff = targetLevel - this->level();
793 result[0] = this->level() + diff;
794 index_t lowerIndex, upperIndex, degree; //, maxKtIndex
795 for (index_t i = 0; i!=d; i++)
796 {
797 degree = m_basis->degree(i);
798 // maxKtIndex = m_basis->tensorLevel(this->level()+1).knots(i).size();
799
800 if (degree % 2 == 1 && degree>1)
801 {
802 lowerIndex = this->lowerIndex()[i]*std::pow(2,diff);
803 ( lowerIndex < (degree-1)/2-1 ? lowerIndex=0 : lowerIndex-=(degree-1)/2-1 );
804 result[i+1] = lowerIndex;
805 upperIndex = this->upperIndex()[i]*std::pow(2,diff);
806 // ( upperIndex + (degree)/2+1 >= maxKtIndex ? upperIndex=maxKtIndex-1 : upperIndex+=(degree)/2+1);
807 result[d+i+1] = upperIndex;
808 }
809 else
810 {
811 lowerIndex = this->lowerIndex()[i]*std::pow(2,diff);
812 ( lowerIndex < (degree-1)/2 ? lowerIndex=0 : lowerIndex-=(degree-1)/2 );
813 result[i+1] = lowerIndex;
814 upperIndex = this->upperIndex()[i]*std::pow(2,diff);
815 // ( upperIndex + (degree)/2 >= maxKtIndex ? upperIndex=maxKtIndex-1 : upperIndex+=(degree)/2);
816 result[d+i+1] = upperIndex;
817 }
818 }
819 return result;
820}
821
822template <short_t d, class T>
823typename gsHBox<d,T>::RefBox gsHBox<d, T>::toCrsBox(index_t targetLevel) const
824{
825 GISMO_ASSERT(targetLevel < this->level(),"Target level must be smaller than the current level, but "<<targetLevel<<">="<<this->level());
826 std::vector<index_t> result(2*d+1);
827 index_t diff = this->level() - targetLevel;
828 result[0] = this->level() - diff;
829 for (index_t i = 0; i!=d; i++)
830 {
831 result[i+1] = this->lowerIndex()[i]/std::pow(2,diff);
832 result[d+i+1] = this->upperIndex()[i]/std::pow(2,diff);
833 }
834 return result;
835}
836
837
838template <short_t d, class T>
839typename gsHBox<d,T>::HContainer gsHBox<d, T>::boxUnion(const HContainer & container1, const HContainer & container2) const
840{
841 HContainer result, region1, region2;
842
843 region1 = container1;
844 region2 = container2;
845
846 index_t lmax = std::max(region1.size(),region2.size());
847 region1.resize(lmax);
848 region2.resize(lmax);
849 result.resize(lmax);
850
851 for (index_t l = 0; l!=lmax; l++)
852 result[l] = _boxUnion(region1[l],region2[l]);
853
854 return result;
855}
856
857template <short_t d, class T>
858typename gsHBox<d, T>::Container gsHBox<d, T>::_boxUnion(const Container & container1, const Container & container2) const
859{
860 // SortedContainer sortedResult;
861
862 // SortedContainer scontainer1 = gsHBoxSort<d,T>(container1);
863 // SortedContainer scontainer2 = gsHBoxSort<d,T>(container2);
864
865
866 // struct
867 // {
868 // bool operator()(const gsHBox<d,T> & a, const gsHBox<d,T> & b) const
869 // {
870 // return
871 // (a.patch() < b.patch())
872 // ||
873 // ((a.patch() == b.patch()) &&
874 // (a.level() < b.level()) )
875 // ||
876 // ((a.patch() == b.patch()) &&
877 // (a.level() == b.level()) &&
878 // std::lexicographical_compare( a.lowerIndex().begin(), a.lowerIndex().end(),
879 // b.lowerIndex().begin(), b.lowerIndex().end()) )
880 // ||
881 // ((a.patch() == b.patch()) &&
882 // (a.level() == b.level()) &&
883 // (a.lowerIndex() == b.lowerIndex()) &&
884 // std::lexicographical_compare( a.upperIndex().begin(), a.upperIndex().end(),
885 // b.upperIndex().begin(), b.upperIndex().end()) );
886 // };
887 // }
888 // comp;
889
890 // sortedResult.reserve(scontainer1.size() + scontainer2.size());
891 // if (scontainer1.size()!=0 && scontainer2.size()!=0)
892 // {
893 // // First sort (otherwise union is wrong)
894 // std::sort(scontainer1.begin(),scontainer1.end(),comp);
895 // std::sort(scontainer2.begin(),scontainer2.end(),comp);
896
897 // std::set_union( scontainer1.begin(),scontainer1.end(),
898 // scontainer2.begin(),scontainer2.end(),
899 // std::inserter(sortedResult,sortedResult.begin()),
900 // comp);
901 // }
902 // else if (scontainer1.size()!=0 && container2.size()==0)
903 // {
904 // scontainer1 = SortedContainer (container1.begin(), container1.end());
905 // scontainer2 = SortedContainer (container2.begin(), container2.end());
906 // sortedResult.insert(sortedResult.end(),scontainer1.begin(),scontainer1.end());
907 // }
908 // else if (scontainer1.size()==0 && container2.size()!=0)
909 // {
910
911 // sortedResult.insert(sortedResult.end(),scontainer2.begin(),scontainer2.end());
912 // }
913 // else { /* Do nothing */ }
914
915 return gsHBoxUtils<d,T>::Union(container1,container2);
916}
917
918template <short_t d, class T>
919typename gsHBox<d, T>::Container gsHBox<d, T>::_makeUnique(const Container & container) const
920{
921 // SortedContainer scontainer = gsHBoxSort<d,T>(container);
922
923 // struct
924 // {
925 // bool operator()(const gsHBox<d,T> & a, const gsHBox<d,T> & b) const
926 // {
927 // return a.isSame(b);
928 // };
929 // }
930 // pred;
931
932 // // Get unique entries
933 // typename SortedContainer::iterator it = std::unique(scontainer.begin(),scontainer.end(),pred);
934 // scontainer.resize(distance(scontainer.begin(), it));
935 // Container result(scontainer.begin(),scontainer.end());
936 // return result;
937
938 return gsHBoxUtils<d,T>::Unique(container);
939}
940
941template <short_t d, class T>
942bool gsHBox<d, T>::good() const
943{
944 return (m_indices.first.array() >=0).any() && (m_indices.second.array() >=0).any();
945}
946
947template <short_t d, class T>
948void gsHBox<d, T>::clean(Container & container) const
949{
950 std::function<bool(const gsHBox<d,T> &)> pred = [](const gsHBox<d,T> & box) { return !(box.good()); };
951 Iterator beg = std::remove_if(container.begin(),container.end(),pred);
952 container.erase(beg,container.end());
953}
954
955namespace internal
956{
957
959template<short_t d, class T>
960class gsXml< gsHBox<d,T> >
961{
962private:
963 gsXml() { }
964 typedef gsHBox<d,T> Object;
965public:
966 GSXML_COMMON_FUNCTIONS(Object);
967 static std::string tag () { return "HBox"; }
968 static std::string type () { return "HBox"+std::to_string(d); }
969
970 GSXML_GET_POINTER(Object);
971
972 static void get_into (gsXmlNode * node, Object & obj)
973 {
974 gsXmlNode * tmp = node->first_node("box");
975 GISMO_ASSERT( tmp , "Expected to find a box node.");
976
977 std::vector<index_t> box;
978 box.push_back(atoi( tmp->first_attribute("level")->value() ));
979
980 std::istringstream str;
981 str.clear();
982 str.str( tmp->value() );
983 unsigned c;
984 for( unsigned i = 0; i < 2*d; i++)
985 {
986 str>> c;
987 box.push_back(c);
988 }
989
990 if (tmp->first_attribute("patch"))
991 {
992 index_t patch = atoi(tmp->first_attribute("patch")->value());
993 obj = gsHBox<d,T>(box,nullptr,patch);
994 }
995 else
996 obj = gsHBox<d,T>(box,nullptr);
997 }
998
999 static gsXmlNode * put (const Object & obj,
1000 gsXmlTree & data )
1001 {
1002 gsXmlNode * tmp = makeNode("HBox", data);
1003 tmp->append_attribute( makeAttribute("type",internal::gsXml<Object>::type().c_str(), data) );
1004 // tmp->append_attribute(makeAttribute("level", 0, data));
1005
1006 gsMatrix<index_t> box(1,2*d);
1007 box.leftCols(d) = obj.lowerIndex().transpose();
1008 box.rightCols(d) = obj.upperIndex().transpose();
1009 gsXmlNode * mat = putMatrixToXml( box, data, "box" );
1010 mat->append_attribute(makeAttribute("level", obj.level(), data));
1011 if (obj.patch()!=-1)
1012 mat->append_attribute(makeAttribute("patch", obj.patch(), data));
1013 tmp->append_node(mat);
1014 return tmp;
1015 }
1016};
1017
1018} // internal
1019
1020} // namespace gismo
gismo::gsDomainIterator::m_basis
const gsBasis< T > * m_basis
The basis on which the domain iterator is defined.
Definition gsDomainIterator.h:219
gismo::gsHBSplineBasis
A hierarchical B-spline basis of parametric dimension d.
Definition gsHBSplineBasis.h:36
gismo::gsHBox
This class provides a Hierarchical Box (gsHBox)
Definition gsHBox.h:55
gismo::gsHBox::getMultiLevelSupportExtension
Container getMultiLevelSupportExtension(index_t k)
Gets the multi-level support extension.
Definition gsHBox.hpp:477
gismo::gsHBox::toContainer
Container toContainer()
Returns a container representation of the object.
Definition gsHBox.hpp:427
gismo::gsHBox::getMaxCellLength
T getMaxCellLength() const
Return the length of the largest edge of the element.
Definition gsHBox.hpp:273
gismo::gsHBox::gsHBox
gsHBox()
Default constructor.
Definition gsHBox.hpp:27
gismo::gsHBox::marked
bool marked() const
Returns whether the element is marked or not.
Definition gsHBox.hpp:337
gismo::gsHBox::isActive
bool isActive() const
Determines if the box is active on its current level.
Definition gsHBox.hpp:213
gismo::gsHBox::levelInCenter
index_t levelInCenter() const
Gets the level in the center of the object.
Definition gsHBox.hpp:225
gismo::gsHBox::getSupportExtension
Container getSupportExtension()
Gets the support extension.
Definition gsHBox.hpp:444
gismo::gsHBox::projectedErrorRef
T projectedErrorRef() const
The error contribution of *this when it is refined.
Definition gsHBox.hpp:307
gismo::gsHBox::getNeighborhood
gsHBox< d, T >::Container getNeighborhood(index_t m)
Gets the neighborhood.
Definition gsHBox.h:422
gismo::gsHBox::computeCenter
void computeCenter() const
Computes the center of this.
Definition gsHBox.hpp:649
gismo::gsHBox::index
index_t index() const
Gets the index stored in the object.
Definition gsHBox.hpp:328
gismo::gsHBox::getAncestor
gsHBox< d, T > getAncestor(index_t k) const
Gets the ancestor of the object on level k.
Definition gsHBox.hpp:351
gismo::gsHBox::contains
bool contains(const gsHBox< d, T > &other) const
Checks if the other cell is contains in this cell.
Definition gsHBox.hpp:184
gismo::gsHBox::toBox
RefBox toBox() const
Returns a box representation of the object.
Definition gsHBox.hpp:776
gismo::gsHBox::getTneighborhood
Container getTneighborhood(index_t m)
Gets the T-neighborhood.
Definition gsHBox.hpp:517
gismo::gsHBox::getHneighborhood
Container getHneighborhood(index_t m)
Gets the H-neighborhood.
Definition gsHBox.hpp:493
gismo::gsHBox::level
index_t level() const
Gets the level of the object.
Definition gsHBox.hpp:287
gismo::gsHBox::setIndex
void setIndex(index_t index)
Assigns an index to the object.
Definition gsHBox.hpp:325
gismo::gsHBox::computeCoordinates
void computeCoordinates() const
Computes the parametric coordinates of this.
Definition gsHBox.hpp:633
gismo::gsHBox::getParent
gsHBox< d, T > getParent() const
Gets the parent of the object.
Definition gsHBox.hpp:343
gismo::gsHBox::lowerIndex
const point & lowerIndex() const
Gets the lower index of the box.
Definition gsHBox.hpp:279
gismo::gsHBox::mark
void mark()
Marks this element for refinement.
Definition gsHBox.hpp:331
gismo::gsHBox::projectedImprovement
T projectedImprovement() const
Gives the projected improvement that can be expected for this.
Definition gsHBox.hpp:313
gismo::gsHBox::toHContainer
HContainer toHContainer()
Returns a hierarchical container representation of the object.
Definition gsHBox.hpp:435
gismo::gsHBox::projectedSetBack
T projectedSetBack() const
Gives the projected set-back that can be expected for this.
Definition gsHBox.hpp:319
gismo::gsHBox::getChildren
Container getChildren() const
Gets the children of the object.
Definition gsHBox.hpp:374
gismo::gsHBox::upperIndex
const point & upperIndex() const
Gets the upper index of the box.
Definition gsHBox.hpp:281
gismo::gsHBox::setAndProjectError
void setAndProjectError(T error, index_t alpha=2, index_t beta=0)
Sets the error of the object and compute the projection of the error on a finer mesh....
Definition gsHBox.hpp:296
gismo::gsHBox::operator=
gsHBox< d, T > & operator=(const gsHBox< d, T > &other)
Assignment operator.
Definition gsHBox.hpp:132
gismo::gsHBox::clean
void clean(Container &container) const
Cleans the container from bad elements (see good())
Definition gsHBox.hpp:948
gismo::gsHBox::isContained
bool isContained(const gsHBox< d, T > &other) const
Checks if the other cell is contained in this cell.
Definition gsHBox.hpp:167
gismo::gsHBox::getCextension
Container getCextension(index_t m)
Gets the Coarsening extension, which is the coarsening neighborhood before checking for active elemen...
Definition gsHBox.hpp:557
gismo::gsHBox::good
bool good() const
Checks if the box has positive indices.
Definition gsHBox.hpp:942
gismo::gsHBox::patch
index_t patch() const
Gets the patch ID of the object.
Definition gsHBox.hpp:284
gismo::gsHBox::getMinCellLength
T getMinCellLength() const
Return the length of the smallest edge of the element.
Definition gsHBox.hpp:267
gismo::gsHBox::getCenter
const gsMatrix< T > & getCenter() const
Gets the center of the box.
Definition gsHBox.hpp:238
gismo::gsHBox::getCellSize
T getCellSize() const
Return the diagonal of the element.
Definition gsHBox.hpp:261
gismo::gsHBox::setMark
void setMark(bool mark)
Sets the mark.
Definition gsHBox.hpp:340
gismo::gsHBox::setError
void setError(T error)
Sets the error of the object.
Definition gsHBox.hpp:290
gismo::gsHBox::isActiveOrContained
bool isActiveOrContained() const
Determines if active or contained in the active elemebt. In other words; checks if the level of this ...
Definition gsHBox.hpp:219
gismo::gsHBox::error
T error() const
Gets the error stored in the object.
Definition gsHBox.hpp:304
gismo::gsHBox::print
std::ostream & print(std::ostream &os) const
Prints the object.
Definition gsHBox.hpp:612
gismo::gsHBox::projectedErrorCrs
T projectedErrorCrs() const
The error contribution of *this when it is coarsened.
Definition gsHBox.hpp:310
gismo::gsHBox::getCoordinates
const gsMatrix< T > & getCoordinates() const
Gets the coordinates of the box (first column lower corner, second column higher corner).
Definition gsHBox.hpp:231
gismo::gsHBox::getDescendants
Container getDescendants(index_t k) const
Gets the descendants of the object on level k.
Definition gsHBox.hpp:382
gismo::gsHBox::getCneighborhood
Container getCneighborhood(index_t m)
Gets the Coarsening neighborhood.
Definition gsHBox.hpp:591
gismo::gsHBox::upperCorner
gsVector< T, d > upperCorner() const
Gets the upper corner of the box.
Definition gsHBox.hpp:254
gismo::gsHBox::unmark
void unmark()
Unmarks this element for refinement.
Definition gsHBox.hpp:334
gismo::gsHBox::lowerCorner
gsVector< T, d > lowerCorner() const
Gets the lower corner of the box.
Definition gsHBox.hpp:247
gismo::gsHBox::isSame
bool isSame(const gsHBox< d, T > &other) const
Determines whether the this box is the same as the other box.
Definition gsHBox.hpp:199
gismo::gsHBox::toUnitBoxes
Container toUnitBoxes() const
Returns unit boxes representation of the object.
Definition gsHBox.hpp:756
gismo::gsHDomainIterator
Re-implements gsDomainIterator for iteration over all boundary elements of a hierarchical parameter d...
Definition gsHDomainIterator.h:39
gismo::gsHDomainIterator::upperCorner
const gsVector< T > & upperCorner() const
Returns the upper corner of the current element.
Definition gsHDomainIterator.h:122
gismo::gsHDomainIterator::lowerCorner
const gsVector< T > & lowerCorner() const
Returns the lower corner of the current element.
Definition gsHDomainIterator.h:120
gismo::gsHTensorBasis
Class representing a (scalar) hierarchical tensor basis of functions .
Definition gsHTensorBasis.h:75
gismo::gsKnotVector
Class for representing a knot vector.
Definition gsKnotVector.h:80
gismo::gsKnotVector::domainUBegin
uiterator domainUBegin() const
Returns a unique iterator pointing to the starting knot of the domain.
Definition gsKnotVector.h:359
gismo::gsKnotVector::domainUEnd
uiterator domainUEnd() const
Returns a unique iterator pointing to the ending knot of the domain.
Definition gsKnotVector.h:367
gismo::gsKnotVector::uValue
T uValue(const size_t &i) const
Provides the knot with unique index i.
Definition gsKnotVector.h:264
gismo::gsMatrix
A matrix with arbitrary coefficient type and fixed or dynamic size.
Definition gsMatrix.h:41
gismo::gsTHBSplineBasis
Truncated hierarchical B-spline basis.
Definition gsTHBSplineBasis.h:36
gismo::gsVector::at
T at(index_t i) const
Returns the i-th element of the vector.
Definition gsVector.h:177
gismo::nextLexicographic
bool nextLexicographic(Vec &cur, const Vec &size)
Iterates through a tensor lattice with the given size. Updates cur and returns true if another entry ...
Definition gsCombinatorics.h:196
gsAABB.h
Declarations of the class gsAABB, i.e., the axis-aligned bounding box.
gsCombinatorics.h
Provides combinatorial unitilies.
index_t
#define index_t
Definition gsConfig.h:32
GISMO_ERROR
#define GISMO_ERROR(message)
Definition gsDebug.h:118
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.
gsXml.h
Provides declaration of input/output XML utilities struct.
gismo::internal::putMatrixToXml
gsXmlNode * putMatrixToXml(gsMatrix< T > const &mat, gsXmlTree &data, std::string name="Matrix")
Helper to insert matrices into XML.
Definition gsXml.hpp:103
gismo::internal::makeNode
gsXmlNode * makeNode(const std::string &name, gsXmlTree &data)
Helper to allocate XML node.
Definition gsXml.cpp:54
gismo::internal::makeAttribute
gsXmlAttribute * makeAttribute(const std::string &name, const std::string &value, gsXmlTree &data)
Helper to allocate XML attribute.
Definition gsXml.cpp:37
gismo
The G+Smo namespace, containing all definitions for the library.
gismo::give
S give(S &x)
Definition gsMemory.h:266
gismo::gsAABB
Struct of for an Axis-aligned bounding box.
Definition gsAABB.h:31
gismo::gsHBoxUtils
The gsHBoxUtils provide basic utilities to modify HBoxes.
Definition gsHBoxUtils.h:46