1 /** @file gsHTensorBasis.hpp
2
3 @brief Provides implementation of HTensorBasis common operations.
4
5 This file is part of the G+Smo library.
6
7 This Source Code Form is subject to the terms of the Mozilla Public
8 License, v. 2.0. If a copy of the MPL was not distributed with this
9 file, You can obtain one at http://mozilla.org/MPL/2.0/.
10
11 Author(s): G. Kiss, A. Mantzaflaris, J. Speh
12 */
13
14 #pragma once
15
16 #include <gsUtils/gsMesh/gsMesh.h>
17
18 #include <gsIO/gsXml.h>
19 #include <gsIO/gsXmlGenericUtils.hpp>
20
21 namespace gismo
22 {
23 // Central element implementation
24 // template<short_t d, class T>
25 // gsMatrix<T> gsHTensorBasis<d,T>::elementInSupportOf(index_t i) const
26 // {
27 // index_t lvl = levelOf(i);
28 // index_t j = flatTensorIndexOf(i);
29 // return m_bases[lvl]->elementInSupportOf(j);
30 // }
31
32 template<short_t d, class T>
support() const33 gsMatrix<T> gsHTensorBasis<d,T>::support() const
34 {
35 return m_bases[0]->support();
36 }
37
38 //i in continuous numbering
39 template<short_t d, class T>
support(const index_t & i) const40 gsMatrix<T> gsHTensorBasis<d,T>::support(const index_t & i) const
41 {
42 // Get the level
43 int lvl = levelOf(i);
44 // Return the the support
45 return m_bases[lvl]->support( m_xmatrix[lvl][ i - m_xmatrix_offset[lvl] ] );
46 }
47
48 // S.K.
49 template<short_t d, class T> inline
getLevelAtPoint(const gsMatrix<T> & Pt) const50 index_t gsHTensorBasis<d,T>::getLevelAtPoint(const gsMatrix<T> & Pt) const
51 {
52 GISMO_ASSERT(Pt.cols() == 1, "Waiting for single point");
53 point loIdx;
54
55 const int maxLevel = m_tree.getMaxInsLevel();
56
57 for( int i =0; i < Dim; i++)
58 loIdx[i] = m_bases[maxLevel]->knots(i).uFind( Pt(i,0) ).uIndex();
59
60 return m_tree.levelOf( loIdx, maxLevel);
61 }
62
63 template<short_t d, class T> inline
getLevelUniqueSpanAtPoints(const gsMatrix<T> & Pt,gsVector<index_t> & lvl,gsMatrix<index_t> & loIdx) const64 void gsHTensorBasis<d,T>::getLevelUniqueSpanAtPoints(const gsMatrix<T> & Pt,
65 gsVector<index_t> & lvl,
66 gsMatrix<index_t> & loIdx ) const
67 {
68 lvl.resize( Pt.cols() );
69 loIdx.resize( Pt.rows(), Pt.cols() );
70 lvl.setZero();
71 loIdx.setZero();
72 for( index_t i = 0; i < Pt.cols(); i++)
73 {
74 lvl[i] = getLevelAtPoint( Pt.col(i) );
75 for( index_t j = 0; j < Pt.rows(); j++)
76 loIdx(j,i) = m_bases[ lvl[i] ]->knots(j).uFind( Pt(j,i) ).uIndex() ;
77 }
78 }
79
80 template<short_t d, class T> inline
numActive_into(const gsMatrix<T> & u,gsVector<index_t> & result) const81 void gsHTensorBasis<d,T>::numActive_into(const gsMatrix<T> & u, gsVector<index_t>& result) const
82 {
83 result.resize( u.cols() );
84 result.setZero();
85
86 point low, upp, cur;
87 const int maxLevel = m_tree.getMaxInsLevel();
88
89 for(index_t p = 0; p < u.cols(); p++ ) //for all input points
90 {
91 for(short_t i = 0; i != d; ++i)
92 low[i] = m_bases[maxLevel]->knots(i).uFind(u(i,p)).uIndex();
93
94 // Identify the level of the point
95 const int lvl = m_tree.levelOf(low, maxLevel);
96
97 for(int i = 0; i <= lvl; i++)
98 {
99 m_bases[i]->active_cwise(u.col(p), low, upp);
100 cur = low;
101 do //iterate over all points in [low,upp]
102 {
103 CMatrix::const_iterator it =
104 m_xmatrix[i].find_it_or_fail( m_bases[i]->index(cur) );
105
106 if( it != m_xmatrix[i].end() )// if index is found
107 result[p]++;
108 }
109 while( nextCubePoint(cur,low,upp) );
110 }
111 }
112 }
113
114 template<short_t d, class T>
addConnectivity(int lvl,gsMesh<T> & mesh) const115 void gsHTensorBasis<d, T>::addConnectivity(int lvl, gsMesh<T> & mesh) const
116 {
117 const gsVector<index_t, d> & low = gsVector<index_t, d>::Zero();
118
119 const gsTensorBSplineBasis<d, T> & bb = *m_bases[lvl];
120 const CMatrix & cmat = m_xmatrix[lvl];
121
122 // Last tensor-index in level lvl
123 gsVector<index_t, d> end(d);
124 for (index_t i = 0; i < d; ++i)
125 end(i) = bb.component(i).size() - 1;
126
127 index_t k, s;
128 gsVector<index_t, d> v, upp;
129 for (index_t i = 0; i < d; ++i) // For all axes
130 {
131 s = bb.stride(i);
132 v = low;
133 upp = end;
134 upp[i] = 0; // suppress to face v[i]==0
135
136 do // Insert all edges normal to axis i
137 {
138 k = bb.index(v);
139 for (index_t j = 0; j != end[i]; ++j)
140 {
141 if (cmat.bContains(k) && cmat.bContains(k + s))
142 {
143 // inefficient for now
144 const index_t kInd = m_xmatrix_offset[lvl] +
145 (std::lower_bound(cmat.begin(), cmat.end(), k)
146 - cmat.begin());
147
148 // inefficient for now
149 const index_t kNextInd = m_xmatrix_offset[lvl] +
150 (std::lower_bound(cmat.begin(), cmat.end(), k + s)
151 - cmat.begin());
152
153 mesh.addEdge(kInd, kNextInd);
154 }
155 k += s;
156 }
157 } while (nextCubePoint(v, low, upp));
158 }
159 }
160
161 template<short_t d, class T>
connectivity(const gsMatrix<T> & nodes,int level,gsMesh<T> & mesh) const162 void gsHTensorBasis<d, T>::connectivity(const gsMatrix<T> & nodes, int level, gsMesh<T> & mesh) const
163 {
164 const index_t sz = size();
165 GISMO_ASSERT(nodes.rows() == sz, "Invalid input.");
166
167 // Add all vertices
168 for (index_t i = 0; i< sz; ++i)
169 mesh.addVertex(nodes.row(i).transpose());
170
171 addConnectivity(level, mesh);
172 }
173
174 template<short_t d, class T>
connectivity(const gsMatrix<T> & nodes,gsMesh<T> & mesh) const175 void gsHTensorBasis<d,T>::connectivity(const gsMatrix<T> & nodes, gsMesh<T> & mesh) const
176 {
177 const index_t sz = size();
178 GISMO_ASSERT( nodes.rows() == sz, "Invalid input.");
179
180 // Add vertices
181 for(index_t i = 0; i< sz; ++i )
182 mesh.addVertex( nodes.row(i).transpose() );
183
184 // For all levels
185 for(size_t lvl = 0; lvl <= maxLevel(); lvl++)
186 {
187 addConnectivity(lvl, mesh);
188 }
189 }
190
191 template<short_t d, class T>
size() const192 index_t gsHTensorBasis<d,T>::size() const
193 {
194 return m_xmatrix_offset.back();
195 }
196 template<short_t d, class T>
refine_withCoefs(gsMatrix<T> & coefs,gsMatrix<T> const & boxes)197 void gsHTensorBasis<d,T>::refine_withCoefs(gsMatrix<T> & coefs, gsMatrix<T> const & boxes)
198 {
199 std::vector<gsSortedVector<index_t> > OX = m_xmatrix;
200 refine(boxes);
201 gsSparseMatrix<> transf;
202 this->transfer(OX, transf);
203 gsDebug<<"tranf orig:\n"<<transf<<std::endl;
204 coefs = transf*coefs;
205 }
206
207 template<short_t d, class T>
refineElements_withCoefs(gsMatrix<T> & coefs,std::vector<index_t> const & boxes)208 void gsHTensorBasis<d,T>::refineElements_withCoefs(gsMatrix<T> & coefs,std::vector<index_t> const & boxes)
209 {
210 std::vector<gsSortedVector<index_t> > OX = m_xmatrix;
211 refineElements(boxes);
212 gsSparseMatrix<> transf;
213 this->transfer(OX, transf);
214 //gsDebug<<"tranf orig:\n"<<transf<<std::endl;
215 coefs = transf*coefs;
216 }
217
218 template<short_t d, class T>
refineElements_withTransfer(std::vector<index_t> const & boxes,gsSparseMatrix<T> & tran)219 void gsHTensorBasis<d,T>::refineElements_withTransfer(std::vector<index_t> const & boxes, gsSparseMatrix<T> & tran)
220 {
221 std::vector<gsSortedVector<index_t> > OX = m_xmatrix;
222 this->refineElements(boxes);
223 this->transfer(OX, tran);
224 }
225
226
227 template<short_t d, class T>
refineElements_withCoefs2(gsMatrix<T> & coefs,std::vector<index_t> const & boxes)228 void gsHTensorBasis<d,T>::refineElements_withCoefs2(gsMatrix<T> & coefs,std::vector<index_t> const & boxes)
229 {
230 std::vector<gsSortedVector<index_t> > OX = m_xmatrix;
231 refineElements(boxes);
232 gsSparseMatrix<> transf;
233 this->transfer2(OX, transf);
234 //gsDebug<<"tranf 2:\n"<<transf<<std::endl;
235 coefs = transf*coefs;
236 }
237
238 template<short_t d, class T>
uniformRefine_withCoefs(gsMatrix<T> & coefs,int numKnots,int mul)239 void gsHTensorBasis<d,T>::uniformRefine_withCoefs(gsMatrix<T>& coefs, int numKnots, int mul)
240 {
241 std::vector<gsSortedVector<index_t> > OX = m_xmatrix;
242 //uniformRefine(numKnots);
243 //gsMatrix<> transf;
244 //CMatrix temp;
245 //this->m_xmatrix.insert(this->m_xmatrix.begin(), temp);
246 //gsVector<index_t> level;
247 //gsMatrix<index_t> p1, p2;
248 //this->m_tree.getBoxes(p1,p2,level);
249 std::vector<index_t> boxes;
250 index_t lvl;
251 for ( typename hdomain_type::literator it = m_tree.beginLeafIterator(); it.good(); it.next() )
252 {
253 // gsDebug <<" level : "<< it.level() <<"\n";
254 // gsDebug <<" lower : "<< it.lowerCorner() <<"\n";
255 // gsDebug <<" upper : "<< it.upperCorner() <<"\n";
256
257 lvl = it.level() + 1;
258 const point & l = it.lowerCorner();
259 const point & u = it.upperCorner();
260
261 boxes.push_back(lvl);
262 for( short_t i = 0; i < d; i++)
263 boxes.push_back( l(i) * 2);
264 for( short_t i = 0; i < d; i++)
265 boxes.push_back( u(i) * 2);
266 }
267
268 this->clone()->refineElements_withCoefs(coefs, boxes);
269 this->uniformRefine(numKnots, mul);
270 //this->m_xmatrix.erase(this->m_xmatrix.begin(),this->m_xmatrix.begin()+1);
271 //coefs = transf*coefs;
272 }
273
274 template<short_t d, class T>
refine(gsMatrix<T> const & boxes,int refExt)275 void gsHTensorBasis<d,T>::refine(gsMatrix<T> const & boxes, int refExt)
276 {
277 GISMO_ASSERT(boxes.rows() == d, "refine() needs d rows of boxes.");
278 GISMO_ASSERT(boxes.cols()%2 == 0, "Each box needs two corners but you don't provide refine() with them.");
279
280 #ifndef NDEBUG
281 gsMatrix<T> para = support();
282 for(int i = 0; i < boxes.cols()/2; i++)
283 {
284 for( short_t j = 0; j < d; j++ )
285 {
286 GISMO_ASSERT( para(j,0) <= boxes(j, 2*i) ,
287 "In refine() the first corner is outside the computational domain.");
288 GISMO_ASSERT( para(j,1) >= boxes(j, 2*i+1),
289 "In refine() the second corner is outside the computational domain." );
290 }
291 }
292 #endif
293
294 if( refExt == 0 )
295 {
296 // If there is no refinement-extension, just use the
297 // "regular" refinement function refine( gsMatrix )
298 this->refine( boxes );
299
300 // Make sure there are enough levels
301 needLevel( m_tree.getMaxInsLevel() );
302 }
303 else
304 {
305 // Make an element vector
306 std::vector<index_t> refVector = this->asElements(boxes, refExt);//std::vector<unsigned> refVector = this->asElements(boxes, refExt);
307
308 // ...and refine
309 this->refineElements( refVector );
310 }
311
312 // Update the basis (already done by now)
313 //update_structure();
314 }
315
316 template<short_t d, class T>
asElements(gsMatrix<T> const & boxes,int refExt) const317 std::vector<index_t> gsHTensorBasis<d,T>::asElements(gsMatrix<T> const & boxes, int refExt) const
318 {
319 // If there is a refinement-extension, we will have to use
320 // refineElements( std::vector )
321 //
322 // Each box will be represented by 2*d+1 entries specifying
323 // <level to be refined to>,<lower corner>,<upper corner>
324 const int offset = 2*d+1;
325
326 // Initialize vector of size
327 // "entries per box" times "number of boxes":
328 std::vector<index_t> refVector( offset * boxes.cols()/2 );
329 gsMatrix<T> ctr(d,1);
330
331 // Loop over all boxes:
332 for(index_t i = 0; i < boxes.cols()/2; i++)
333 {
334 ctr = ( boxes.col( 2*i ) + boxes.col( 2*i+1) )/2;
335
336 // Compute the level we want to refine to.
337 // Note that, if the box extends over several elements,
338 // the level at the centerpoint will be taken for reference
339 const int refLevel = getLevelAtPoint( ctr ) + 1;
340
341 // Make sure there are enough levels
342 needLevel( refLevel );
343
344 for(index_t j = 0; j < boxes.rows();j++)
345 {
346 // Convert the parameter coordinates to (unique) knot indices
347 const gsKnotVector<T> & kv = m_bases[refLevel]->knots(j);
348 int k1 = (std::upper_bound(kv.domainUBegin(), kv.domainUEnd(),
349 boxes(j,2*i ) ) - 1).uIndex();
350 int k2 = (std::upper_bound(kv.domainUBegin(), kv.domainUEnd()+1,
351 boxes(j,2*i+1) ) - 1).uIndex();
352
353 // Trivial boxes trigger some refinement
354 if ( k1 == k2)
355 {
356 if (0!=k1) {--k1;}
357 ++k2;
358 }
359
360 // If applicable, add the refinement extension.
361 // Note that extending by one cell on level L means
362 // extending by two cells in level L+1
363 ( k1 < 2*refExt ? k1=0 : k1-=2*refExt );
364 const index_t maxKtIndex = kv.size();
365 ( k2 + 2*refExt >= maxKtIndex ? k2=maxKtIndex-1 : k2+=2*refExt);
366
367 // Store the data...
368 refVector[i*offset] = refLevel;
369 refVector[i*offset+1+j] = k1;
370 refVector[i*offset+1+j+d] = k2;
371 }
372 }
373 // gsDebug<<"begin\n";
374 // for (std::vector<unsigned>::const_iterator i = refVector.begin(); i != refVector.end(); ++i)
375 // std::cout << *i << ' ';
376 // gsDebug<<"end\n";
377
378 return refVector;
379 }
380
381 template<short_t d, class T>
refine(gsMatrix<T> const & boxes)382 void gsHTensorBasis<d,T>::refine(gsMatrix<T> const & boxes)
383 {
384 GISMO_ASSERT(boxes.rows() == d, "refine() needs d rows of boxes.");
385 GISMO_ASSERT(boxes.cols()%2 == 0, "Each box needs two corners but you don't provide refine() with them.");
386
387 #ifndef NDEBUG
388 gsMatrix<T> para = support();
389 for(int i = 0; i < boxes.cols()/2; i++)
390 {
391 for( short_t j = 0; j < d; j++ )
392 {
393 GISMO_ASSERT( para(j,0) <= boxes(j, 2*i) ,
394 "In refine() the first corner is outside the computational domain.");
395 GISMO_ASSERT( para(j,1) >= boxes(j, 2*i+1),
396 "In refine() the second corner is outside the computational domain." );
397 }
398 }
399 #endif
400
401 gsVector<index_t,d> k1, k2;
402 for(index_t i = 0; i < boxes.cols()/2; i++)
403 {
404 // 1. Get a small cell containing the box
405 const int fLevel = m_bases.size()-1;
406
407 for(index_t j = 0; j < k1.size();j++)
408 {
409 const gsKnotVector<T> & kv = m_bases.back()->knots(j);
410 k1[j] = (std::upper_bound(kv.domainUBegin(), kv.domainUEnd(),
411 boxes(j,2*i ) ) - 1).uIndex();
412 k2[j] = (std::upper_bound(kv.domainUBegin(), kv.domainUEnd()+1,
413 boxes(j,2*i+1) ) - 1).uIndex();
414
415 // Trivial boxes trigger some refinement
416 if ( k1[j] == k2[j])
417 {
418 if (0!=k1[j]) {--k1[j];}
419 ++k2[j];
420 }
421 }
422
423 // 2. Find the smallest level in which the box is completely contained
424 //const int level = m_tree.query3(k1,k2,fLevel) + 1;
425 // make sure that the grid is computed ( needLevel(level) )
426 //const tensorBasis & tb = tensorLevel(level);
427 //GISMO_UNUSED(tb);
428
429 // Sink box
430 m_tree.sinkBox(k1, k2, fLevel);
431 // Make sure we have enough levels
432 needLevel( m_tree.getMaxInsLevel() );
433 }
434
435 // Update the basis
436 // update_structure();
437 }
438
439 template<short_t d, class T>
refineBasisFunction(const index_t i)440 void gsHTensorBasis<d,T>::refineBasisFunction(const index_t i)
441 {
442 // Get current level
443 const index_t lvl = this->levelOf(i);
444 // Get the support endpoints
445 gsMatrix<index_t, d, 2> elements;
446 m_bases[lvl]->elementSupport_into(m_xmatrix[lvl][ i - m_xmatrix_offset[lvl] ],
447 elements);
448 point low = elements.col(0);
449 point upp = elements.col(1);
450 // Advance the indices to one level deeper
451 for ( short_t k = 0; k!=d; ++k )
452 {
453 low[k] = low[k] << 1;
454 upp[k] = upp[k] << 1;
455 }
456 // Insert the domain to the lvl+1 nested domain
457 m_tree.insertBox(low,upp,lvl+1);
458 // Update the basis
459 update_structure();
460 }
461
462
463 /*
464 template<short_t d, class T>
465 void gsHTensorBasis<d,T>::refine(gsDomainIterator<T> const & boxes)
466 {
467
468 }
469 */
470
471 template<short_t d, class T>
refineElements(std::vector<index_t> const & boxes)472 void gsHTensorBasis<d,T>::refineElements(std::vector<index_t> const & boxes)
473 {
474 point i1;
475 point i2;
476
477 GISMO_ASSERT( (boxes.size()%(2*d + 1))==0,
478 "The points did not define boxes properly. The boxes were not added to the basis.");
479 for(size_t i = 0; i < (boxes.size())/(2*d+1); i++)
480 {
481 for( short_t j = 0; j < d; j++ )
482 {
483 i1[j] = boxes[(i*(2*d+1))+j+1];
484 i2[j] = boxes[(i*(2*d+1))+d+j+1];
485 }
486 insert_box(i1,i2,boxes[i*(2*d+1)]);
487 }
488
489 update_structure();
490 }
491
492 template<short_t d, class T>
refineSide(const boxSide side,index_t lvl)493 void gsHTensorBasis<d,T>::refineSide(const boxSide side, index_t lvl)
494 {
495 const index_t dir = side.direction();
496 const index_t par = side.parameter();
497 gsMatrix<T> rf = this->support();
498 rf(dir,!par) = rf(dir,par);
499 for (index_t i = 0; i!=lvl; ++i) // lazy impl., this can be more efficient
500 this->refine(rf);
501 }
502
503 template<short_t d, class T>
matchWith(const boundaryInterface & bi,const gsBasis<T> & other,gsMatrix<index_t> & bndThis,gsMatrix<index_t> & bndOther) const504 void gsHTensorBasis<d,T>::matchWith(const boundaryInterface & bi,
505 const gsBasis<T> & other,
506 gsMatrix<index_t> & bndThis,
507 gsMatrix<index_t> & bndOther) const
508 {
509 if( const Self_t * _other = dynamic_cast<const Self_t*>( &other) )
510 {
511 gsVector<index_t> N(d);
512
513 // tens1 will store the tensor-index on side second(),...
514 gsVector<index_t> tens0(d), tens1(d);
515
516 // see if the orientation is preserved on side second()
517 const gsVector<bool> dirOrient = bi.dirOrientation();
518
519 const gsVector<index_t> dirMap = bi.dirMap();
520
521 // get the global indices of the basis functions which are
522 // active on the interface
523 bndThis = this->boundary( bi.first().side() );
524
525 // this is only for checking whether, at least, both involved
526 // bases have the same number of DOF on the interface.
527 bndOther= _other->boundary( bi.second().side() );
528 GISMO_ASSERT( bndThis.rows() == bndOther.rows(),
529 "Input error, sizes do not match: "
530 <<bndThis.rows()<<"!="<<bndOther.rows() );
531 // bndOther gets overwritten completely, so here is the setZero():
532 bndOther.setZero();
533
534 for( index_t i=0; i < bndThis.rows(); i++)
535 {
536 // get the level of the basis function on side first()
537 index_t L = this->levelOf( bndThis(i,0) );
538 // get the flat tensor index
539 // (i.e., the single-number-index on level L)...
540 index_t flat0 = this->flatTensorIndexOf( bndThis(i,0) );
541 // ... and change it to the tensor-index.
542 tens0 = this->tensorLevel(L).tensorIndex( flat0 );
543
544 // ...flat1 the corresponding flat index
545 // (single-number on level)...
546 index_t flat1 = 0;
547 // ...and cont1 the corresponding continued (global) index.
548 index_t cont1 = 0;
549
550 // get the sizes of the components of the tensor-basis on this level,
551 // i.e., the sizes of the univariate bases corresponding
552 // to the respective coordinate directions
553 for( short_t j=0; j < d; j++)
554 N[j] = _other->tensorLevel(L).component(j).size();
555
556 // get the tensor-index of the basis function on level L on
557 // second() that should be matched with flatp/tens0
558 for( short_t j=0; j<d; j++)
559 {
560 // coordinate direction j on first() gets
561 // mapped to direction jj on second()
562 index_t jj = dirMap[j];
563 // store the respective component of the tensor-index
564 tens1[jj] = tens0[j];
565
566 if( jj == bi.second().direction() )
567 {
568 // if jj is the direction() of the interface,
569 // however, we need either the first
570 // or last basis function
571 if( bi.second().parameter() ) // true = 1 = end
572 tens1[jj] = N[jj]-1;
573 else
574 tens1[jj] = 0;
575 }
576 else
577 {
578 // otherwise, check if the orientation is
579 // preserved. If necessary, flip it.
580 if( !dirOrient[j] )
581 tens1[jj] = N[jj]-1 - tens1[jj];
582 }
583
584 }
585
586 flat1 = _other->tensorLevel(L).index( tens1 );
587
588 // compute the "continuous" index on second(), i.e., the index
589 // in the numbering which is global over all levels.
590 cont1 = _other->flatTensorIndexToHierachicalIndex( flat1, L );
591 // this is the index that has to be matched with bndThis(i,0)
592 bndOther( i, 0 ) = cont1;
593 }
594 return;
595 }
596 gsWarn<<"Cannot match with "<< other <<"\n";
597 }
598
599 //protected functions
600
601 // Construct the characteristic matrix of level \a level ; i.e., set
602 // all the matrix entries corresponding to active functions to one and
603 // the rest to zero.
604 template<short_t d, class T>
set_activ1(int level)605 void gsHTensorBasis<d,T>::set_activ1(int level)
606 {
607 typedef typename gsKnotVector<T>::smart_iterator knotIter;
608
609 //gsDebug<<" Setting level "<< level <<"\n";
610 point low, upp;
611
612 CMatrix & cmat = m_xmatrix[level];
613
614 // Clear previous entries
615 cmat.clear();
616
617 // If a level is to be checked which is larger than
618 // the maximum inserted level, nothing need so to be done
619 if ( level > static_cast<int>(m_tree.getMaxInsLevel() ) )
620 return;
621
622
623 gsVector<knotIter,d> starts, ends, curr;
624 point ind;
625 ind[0] = 0; // for d==1: warning: may be used uninitialized in this function (snap-ci)
626
627 for(short_t i = 0; i != d; ++i)
628 {
629 // beginning of the iteration in i-th direction
630 starts[i] = m_bases[level]->knots(i).sbegin() ;
631 // end of the iteration in i-th direction
632 ends [i] = m_bases[level]->knots(i).send()-m_deg[i]-1;
633 }
634
635 curr = starts;// start iteration
636 do
637 {
638 for(short_t i = 0; i != d; ++i)
639 {
640 low[i] = curr[i].uIndex(); // lower left corner of the support of the function
641 upp[i] = (curr[i]+m_deg[i]+1).uIndex(); // upper right corner of the support
642 ind[i] = curr[i].index(); // index of the function in the matrix
643 }
644
645 if ( m_tree.query3(low, upp,level) == level) //if active
646 cmat.push_unsorted( m_bases[level]->index( ind ) );
647
648 }
649 while ( nextLexicographic(curr,starts,ends) ); // while there are some functions (i.e., some combinations of iterators) left
650
651 cmat.sort();
652
653 }
654
655 template<short_t d, class T>
functionOverlap(const point & boxLow,const point & boxUpp,const int level,point & actLow,point & actUpp)656 void gsHTensorBasis<d,T>::functionOverlap(const point & boxLow, const point & boxUpp,
657 const int level, point & actLow, point & actUpp)
658 {
659 const tensorBasis & tb = *m_bases[level];
660 for(short_t i = 0; i != d; ++i)
661 {
662 actLow[i] = tb.knots(i).lastKnotIndex (boxLow[i]) - m_deg[i];
663 actUpp[i] = tb.knots(i).firstKnotIndex(boxUpp[i]) - 1 ;
664
665 // Note aao:
666 //actLow[i] = firstKnotIndex(boxLow[i]);
667 //actUpp[i] = tb.knots(i).lastKnotIndex(boxUpp[i])-m_deg[i]-1;
668 }
669 }
670
671 template<short_t d, class T>
setActive()672 void gsHTensorBasis<d,T>::setActive()
673 {
674 // for(size_t lvl = 0; lvl != m_xmatrix.size(); lvl++)
675 // set_activ1(lvl);
676 // return;
677
678 // iterate over leaf-boxes
679 // for all overlapping supports with the box
680 // set obvious to active
681 // for the rest candidates (supp. not fully contained in box ~ !query2)
682 // (equiv: actives on the boundary cells of the box)
683 // query3(supp,box.level) == level (min. is level: no coarser)
684 // take care: duplicates from different leaves or adj. cells
685 point curr, actUpp;
686 gsMatrix<index_t,d,2> elSupp;
687
688 // try: iteration per level
689 for ( typename hdomain_type::literator it = m_tree.beginLeafIterator();
690 it.good(); it.next() )
691 {
692 const int lvl = it.level();
693 CMatrix & cmat = m_xmatrix[lvl];
694
695 // Get candidate functions
696 functionOverlap(it.lowerCorner(), it.upperCorner(), lvl, curr, actUpp);
697
698 do
699 {
700 const index_t gi = m_bases[lvl]->index( curr );
701
702 // Get element support
703 m_bases[lvl]->elementSupport_into(gi, elSupp);
704
705 if ( (elSupp.col(0).array() >= it.lowerCorner().array()).all() &&
706 (elSupp.col(1).array() <= it.upperCorner().array()).all() )
707 {
708 // to do: all-at-once
709 cmat.push_unsorted( gi );
710 }
711 else
712 {
713 // Check if active (iff no overlap with level less than lvl)
714 if ( m_tree.query3(elSupp.col(0), elSupp.col(1), lvl) == lvl)
715 cmat.push_unsorted( gi );
716 }
717 }
718 while( nextCubePoint(curr, actUpp) );
719 }
720
721 for(size_t lvl = 0; lvl != m_xmatrix.size(); ++lvl)
722 {
723 m_xmatrix[lvl].sort();
724 m_xmatrix[lvl].erase( std::unique( m_xmatrix[lvl].begin(), m_xmatrix[lvl].end() ),
725 m_xmatrix[lvl].end() );
726 }
727 }
728
729 template<short_t d, class T>
setActiveToLvl(int level,std::vector<CMatrix> & x_matrix_lvl) const730 void gsHTensorBasis<d,T>::setActiveToLvl(int level,
731 std::vector<CMatrix> & x_matrix_lvl) const
732 {
733 x_matrix_lvl.resize(level+1);
734
735 // vectors of iterators. one iterator for each dimension
736 gsVector<typename gsKnotVector<T>::smart_iterator,d> starts, ends, curr;
737 gsVector<index_t,d> ind;
738 ind[0] = 0; // for d==1: warning: may be used uninitialized in this function (snap-ci)
739 gsVector<index_t,d> low, upp;
740 for(int j =0; j < level+1; j++)
741 {
742 // Clear previous entries
743 x_matrix_lvl[j].clear();
744
745 for(index_t i = 0; i != d; ++i)
746 {
747 // beginning of the iteration in i-th direction
748 starts[i] = m_bases[j]->knots(i).sbegin() ;
749 // end of the iteration in i-th direction
750 ends [i] = m_bases[j]->knots(i).send()-m_deg[i]-1;
751 }
752
753 curr = starts; // set start of iteration
754 do
755 {
756 for(index_t i = 0; i != d; ++i)
757 {
758 low[i] = curr[i].uIndex(); // lower left corner of the support of the function
759 upp[i] = (curr[i]+m_deg[i]+1).uIndex(); // upper right corner of the support
760 ind[i] = curr[i].index(); // index of the function in the matrix
761 }
762 if(j < level)
763 {
764 if ( m_tree.query3(low, upp,j) == j)
765 { //if active
766 x_matrix_lvl[j].push_unsorted( m_bases[j]->index( ind ) );
767 }
768 }else{
769 if ( m_tree.query3(low, upp,j) >= j)
770 { //if active
771 x_matrix_lvl[j].push_unsorted( m_bases[j]->index( ind ) );
772 }
773 }
774 }
775 while ( nextLexicographic(curr,starts,ends) ); // while there are some functions (i.e., some combinations of iterators) left
776
777 x_matrix_lvl[j].sort();
778 }
779
780 }
781
782 ///private functions
783 template<short_t d, class T> inline
insert_box(typename gsHTensorBasis<d,T>::point const & k1,typename gsHTensorBasis<d,T>::point const & k2,int lvl)784 void gsHTensorBasis<d,T>::insert_box(typename gsHTensorBasis<d,T>::point const & k1,
785 typename gsHTensorBasis<d,T>::point const & k2,
786 int lvl )
787 {
788 // Remember box in History (for debugging)
789 // m_boxHistory.push_back( box(k1,k2,lvl) );
790
791 m_tree.insertBox(k1,k2, lvl);
792 needLevel( m_tree.getMaxInsLevel() );
793 }
794
795 template<short_t d, class T>
makeCompressed()796 void gsHTensorBasis<d,T>::makeCompressed()
797 {
798 // Compress the tree
799 // m_tree.makeCompressed();
800
801 while ( ! m_xmatrix_offset[1] )
802 {
803 delete m_bases.front();
804 m_bases.erase( m_bases.begin() );
805 m_tree.decrementLevel();
806 m_xmatrix.erase( m_xmatrix.begin() );
807 m_xmatrix_offset.erase( m_xmatrix_offset.begin() );
808 }
809 // Note/to do: cleaning up empty levels at the end as well.
810 }
811
812 template<short_t d, class T>
flatTensorIndexesToHierachicalIndexes(gsSortedVector<int> & indexes,const int level) const813 void gsHTensorBasis<d,T>::flatTensorIndexesToHierachicalIndexes(gsSortedVector< int > & indexes,const int level) const
814 {
815 GISMO_ASSERT( static_cast<size_t>(level) < m_xmatrix.size(), "Requested level does not exist.\n");
816
817 CMatrix::const_iterator xmat_pointer = m_xmatrix[level].begin();
818 CMatrix::const_iterator xmat_end = m_xmatrix[level].end();
819 gsSortedVector< int >::iterator ind_pointer = indexes.begin();
820 gsSortedVector< int >::iterator ind_end = indexes.end();
821 index_t index = 0;
822 while(ind_pointer!=ind_end&&xmat_pointer!=xmat_end)
823 {
824 if(*ind_pointer<static_cast<int>(*xmat_pointer))
825 {
826 (*ind_pointer)=-1;
827 ++ind_pointer;
828 }
829 else if(*ind_pointer==static_cast<int>(*xmat_pointer))
830 {
831 (*ind_pointer)=m_xmatrix_offset[level]+index;
832 ++xmat_pointer;
833 ++index;
834 ++ind_pointer;
835 }
836 else
837 {
838 ++xmat_pointer;
839 ++index;
840 }
841 }
842 while(ind_pointer!=ind_end)
843 {
844 (*ind_pointer)=-1;
845 ++ind_pointer;
846 }
847 }
848
849 template<short_t d, class T>
flatTensorIndexToHierachicalIndex(index_t index,const int level) const850 int gsHTensorBasis<d,T>::flatTensorIndexToHierachicalIndex(index_t index,const int level) const
851 {
852 if( m_xmatrix.size()<=static_cast<size_t>(level) )
853 return -1;
854 CMatrix::const_iterator it = std::lower_bound(m_xmatrix[level].begin(), m_xmatrix[level].end(), index );
855 if(it == m_xmatrix[level].end() || *it != index)
856 return -1;
857 else
858 return m_xmatrix_offset[level] + ( it - m_xmatrix[level].begin() );
859 }
860
861 template<short_t d, class T>
activeBoundaryFunctionsOfLevel(const unsigned level,const boxSide & s,std::vector<bool> & actives) const862 void gsHTensorBasis<d,T>::activeBoundaryFunctionsOfLevel(const unsigned level,const boxSide & s,std::vector<bool>& actives) const
863 {
864 needLevel( level );
865
866 const gsMatrix<index_t> bound = m_bases[level]->boundary(s);
867 const index_t sz = bound.rows();
868 //gsSortedVector< int > indexes(bound->data(),bound->data()+sz);
869 gsSortedVector< int > indexes;
870 indexes.resize(sz,-1);
871 if(level<=maxLevel())
872 {
873 for(index_t i = 0;i<sz;++i)
874 indexes[i]=(bound)(i,0);
875 flatTensorIndexesToHierachicalIndexes(indexes,level);
876 }
877 actives.resize(indexes.size(),false);
878 std::fill (actives.begin(),actives.end(),false);
879 for(size_t i = 0;i<indexes.size();i++)
880 if(indexes[i]!=-1)
881 actives[i]=true;
882 }
883
884 template<short_t d, class T>
update_structure()885 void gsHTensorBasis<d,T>::update_structure() // to do: rename as updateHook
886 {
887 // Make sure we have computed enough levels
888 needLevel( m_tree.getMaxInsLevel() );
889
890 // Setup the characteristic matrices
891 m_xmatrix.clear();
892 m_xmatrix.resize( m_bases.size() );
893
894 // Compress the tree
895 m_tree.makeCompressed();
896
897 for(size_t i = 0; i != m_xmatrix.size(); i ++)
898 set_activ1(i);
899
900 // Store all indices of active basis functions to m_matrix
901 //setActive();
902
903 // Compute offsets
904 m_xmatrix_offset.clear();
905 m_xmatrix_offset.reserve(m_xmatrix.size()+1);
906 m_xmatrix_offset.push_back(0);
907 for (size_t i = 0; i != m_xmatrix.size(); i++)
908 {
909 m_xmatrix_offset.push_back(
910 m_xmatrix_offset.back() + m_xmatrix[i].size() );
911 }
912 }
913
914 template<short_t d, class T>
needLevel(int maxLevel) const915 void gsHTensorBasis<d,T>::needLevel(int maxLevel) const
916 {
917 // +1 for the initial basis in m_bases
918 const int extraLevels = maxLevel + 1 - m_bases.size();
919
920 for ( int i = 0; i < extraLevels; ++i )
921 {
922 tensorBasis * next_basis = m_bases.back()->clone().release();
923 next_basis->uniformRefine(1);
924 m_bases.push_back (next_basis); //note: m_bases is mutable
925 }
926 }
927
928 template<short_t d, class T>
initialize_class(gsBasis<T> const & tbasis)929 void gsHTensorBasis<d,T>::initialize_class(gsBasis<T> const& tbasis)
930 {
931 // Degrees
932 //m_deg = tbasis.cwiseDegree();
933 m_deg.resize(d);
934 for( index_t i = 0; i < d; i++)
935 m_deg[i] = tbasis.degree(i);
936
937 // Construct the initial basis
938 if ( const tensorBasis * tb2 =
939 dynamic_cast<const tensorBasis*>(&tbasis) )
940 {
941 m_bases.push_back(tb2->clone().release());
942 }
943 else
944 {
945 GISMO_ERROR("Cannot construct a Hierarchical basis from "<< tbasis );
946 }
947
948 // Initialize the binary tree
949 point upp;
950 for ( index_t i = 0; i!=d; ++i )
951 upp[i] = m_bases[0]->knots(i).numElements();
952
953 m_tree.init(upp);
954
955 // Produce a couple of tensor-product spaces by dyadic refinement
956 m_bases.reserve(3);
957 for(index_t i = 1; i <= 2; i++)
958 {
959 tensorBasis* next_basis = m_bases[i-1]->clone().release();
960 next_basis->uniformRefine(1);
961 m_bases.push_back( next_basis );
962 }
963
964 }
965
966
967 template<short_t d, class T>
active_into(const gsMatrix<T> & u,gsMatrix<index_t> & result) const968 void gsHTensorBasis<d,T>::active_into(const gsMatrix<T> & u, gsMatrix<index_t>& result) const
969 {
970 gsMatrix<T> currPoint;
971 point low, upp, cur;
972 const int maxLevel = m_tree.getMaxInsLevel();
973
974 std::vector<std::vector<index_t> > temp_output;//collects the outputs
975 temp_output.resize( u.cols() );
976 size_t sz = 0;
977
978 //gsMatrix<index_t> activesLvl;
979
980 for(index_t p = 0; p < u.cols(); p++) //for all input points
981 {
982 currPoint = u.col(p);
983 for(short_t i = 0; i != d; ++i)
984 low[i] = m_bases[maxLevel]->knots(i).uFind( currPoint(i,0) ).uIndex();
985
986 // Identify the level of the point
987 const int lvl = m_tree.levelOf(low, maxLevel);
988
989 for(int i = 0; i <= lvl; i++)
990 {
991 /*
992 m_bases[i]->active_into(currPoint, activesLvl);
993
994 std::set_intersection(m_xmatrix[i].begin(), m_xmatrix[i].end(),
995 activesLvl.data(), activesLvl.data() + activesLvl.size(),
996 std::back_inserter( temp_output[p] ) );
997 +++ Renumbering to H-basis indexing
998 // */
999
1000 // /*
1001 m_bases[i]->active_cwise(currPoint, low, upp);
1002 cur = low;
1003 do
1004 {
1005 CMatrix::const_iterator it =
1006 m_xmatrix[i].find_it_or_fail( m_bases[i]->index(cur) );
1007
1008 if( it != m_xmatrix[i].end() )// if index is found
1009 {
1010 temp_output[p].push_back(
1011 this->m_xmatrix_offset[i] + (it - m_xmatrix[i].begin() )
1012 );
1013 }
1014 }
1015 while( nextCubePoint(cur,low,upp) );
1016 //*/
1017 }
1018
1019 // update result size
1020 if ( temp_output[p].size() > sz )
1021 sz = temp_output[p].size();
1022 }
1023
1024 result.resize(sz, u.cols() );
1025 for(index_t i = 0; i < result.cols(); i++)
1026 {
1027 result.col(i).topRows(temp_output[i].size())
1028 = gsAsConstVector<index_t>(temp_output[i]);
1029 result.col(i).bottomRows(sz-temp_output[i].size()).setZero();
1030 }
1031 }
1032
1033 template<short_t d, class T>
allBoundary() const1034 gsMatrix<index_t> gsHTensorBasis<d,T>::allBoundary( ) const
1035 {
1036 std::vector<index_t> temp;
1037 gsVector<index_t, d> ind;
1038 for(unsigned i = 0; i <= this->maxLevel(); i++)
1039 for (CMatrix::const_iterator it = m_xmatrix[i].begin();
1040 it != m_xmatrix[i].end(); it++)
1041 {
1042 ind = this->m_bases[i]->tensorIndex(*it);
1043 for (unsigned j=0; j!=d; ++j )
1044 if ( (ind[j]==0) || (ind[j]==(this->m_bases[i]->size(j)-1)) )
1045 {
1046 temp.push_back(m_xmatrix_offset[i] + (it-m_xmatrix[i].begin()) );
1047 break;
1048 }
1049 }
1050 return makeMatrix<index_t>(temp.begin(),temp.size(),1 );
1051 }
1052
1053 template<short_t d, class T>
1054 gsMatrix<index_t> gsHTensorBasis<d,T>::
boundaryOffset(boxSide const & s,index_t offset) const1055 boundaryOffset(boxSide const & s,index_t offset) const
1056 {
1057 //get information on the side
1058 const index_t k = s.direction();
1059 const bool par = s.parameter();
1060
1061 std::vector<index_t> temp;
1062 gsVector<index_t,d> ind;
1063 // i goes through all levels of the hierarchical basis
1064 for(unsigned i = 0; i <= this->maxLevel(); i++)
1065 {
1066 GISMO_ASSERT(static_cast<int>(offset)<this->m_bases[i]->size(k),
1067 "Offset cannot be bigger than the amount of basis"
1068 "functions orthogonal to Boxside s!");
1069
1070 index_t r = ( par ? this->m_bases[i]->size(k) - 1 -offset : offset);
1071 for (CMatrix::const_iterator it = m_xmatrix[i].begin();
1072 it != m_xmatrix[i].end(); it++)
1073 {
1074 ind = this->m_bases[i]->tensorIndex(*it);
1075 if ( ind[k]==r )
1076 temp.push_back(
1077 m_xmatrix_offset[i] + (it - m_xmatrix[i].begin() )
1078 );
1079 }
1080 }
1081 return makeMatrix<index_t>(temp.begin(),temp.size(),1 );
1082 }
1083
1084 /*
1085 template<short_t d, class T>
1086 void gsHTensorBasis<d,T>::evalAllDers_into(const gsMatrix<T> & u, int n,
1087 std::vector<gsMatrix<T> >& result) const;
1088 {
1089 result.resize(n+1);
1090
1091 }
1092 */
1093
1094 template<short_t d, class T>
uniformRefine(int numKnots,int mul)1095 void gsHTensorBasis<d,T>::uniformRefine(int numKnots, int mul)
1096 {
1097 GISMO_UNUSED(numKnots);
1098 GISMO_ASSERT(numKnots == 1, "Only implemented for numKnots = 1");
1099
1100 GISMO_ASSERT( m_tree.getMaxInsLevel() < static_cast<unsigned>(m_bases.size()),
1101 "Problem with max inserted levels: "<< m_tree.getMaxInsLevel()
1102 <<"<" << m_bases.size() <<"\n");
1103
1104 // Delete the first level
1105 delete m_bases.front();
1106 m_bases.erase( m_bases.begin() );
1107
1108 // Keep consistency of finest level
1109 tensorBasis * last_basis = m_bases.back()->clone().release();
1110 last_basis->uniformRefine(1,mul);
1111 m_bases.push_back( last_basis );
1112
1113 // Lift all indices in the tree by one level
1114 m_tree.multiplyByTwo();
1115
1116 update_structure();
1117 }
1118
1119 template<short_t d, class T>
domainBoundariesParams(std::vector<std::vector<std::vector<std::vector<T>>>> & result) const1120 std::vector< std::vector< std::vector<index_t > > > gsHTensorBasis<d,T>::domainBoundariesParams( std::vector< std::vector< std::vector< std::vector< T > > > >& result) const
1121 {
1122 std::vector< std::vector< std::vector< std::vector<index_t > > > > dummy;
1123 return domainBoundariesGeneric( dummy, result, false );
1124 }
1125
1126 template<short_t d, class T>
domainBoundariesIndices(std::vector<std::vector<std::vector<std::vector<index_t>>>> & result) const1127 std::vector< std::vector< std::vector<index_t > > > gsHTensorBasis<d,T>::domainBoundariesIndices( std::vector< std::vector< std::vector< std::vector<index_t > > > >& result ) const
1128 {
1129 std::vector< std::vector< std::vector< std::vector< T > > > > dummy;
1130 return domainBoundariesGeneric( result, dummy, true );
1131 }
1132
1133
1134 template<short_t d, class T>
domainBoundariesGeneric(std::vector<std::vector<std::vector<std::vector<index_t>>>> & indices,std::vector<std::vector<std::vector<std::vector<T>>>> & params,bool indicesFlag) const1135 std::vector< std::vector< std::vector<index_t > > > gsHTensorBasis<d,T>::domainBoundariesGeneric(std::vector< std::vector< std::vector< std::vector<index_t > > > >& indices,
1136 std::vector< std::vector< std::vector< std::vector< T > > > >& params,
1137 bool indicesFlag ) const
1138 {
1139 indices.clear();
1140 params.clear();
1141 std::vector< std::vector< std::vector< int > > > res_aabb;
1142 std::vector< std::vector< std::vector<index_t > > > res_aabb_unsigned;
1143 std::vector< std::vector< std::vector< std::vector< index_t > > > > polylines;
1144
1145 polylines = this->m_tree.getPolylines();
1146 res_aabb.resize( polylines.size() );
1147 // We cannot simply say
1148 // result = this->m_tree.getPolylines();
1149 // because the return value of getPolylines() are vectors of ints and we have no implicit conversion of int to T.
1150
1151 // We want indices/params to be of the same size as polylines. We achieve this here and in the for cycles.
1152 if( indicesFlag )
1153 indices.resize( polylines.size() );
1154
1155 else
1156 params.resize(polylines.size());
1157
1158 int maxLevel = static_cast<int>( this->maxLevel() );
1159 // We precompute the parameter values corresponding to indices of m_maxInsLevel
1160 // although we don't need them if indicesFlag == true.
1161 std::vector<T> x_dir(m_bases[maxLevel]->knots(0).unique());
1162 std::vector<T> y_dir(m_bases[maxLevel]->knots(1).unique());
1163
1164 for(unsigned int i0 = 0; i0 < polylines.size(); i0++)
1165 {
1166 if( indicesFlag )
1167 indices[i0].resize( polylines[i0].size() );
1168 else
1169 params[i0].resize( polylines[i0].size() );
1170
1171 res_aabb[i0].resize( polylines[i0].size() );
1172 for(unsigned int i1 = 0; i1 < polylines[i0].size(); i1++)
1173 {
1174 if( indicesFlag )
1175 indices[i0][i1].resize( polylines[i0][i1].size() );
1176 else
1177 params[i0][i1].resize( polylines[i0][i1].size() );
1178
1179 res_aabb[i0][i1].resize( 4 );
1180 res_aabb[i0][i1][0] = 1000000;
1181 res_aabb[i0][i1][1] = 1000000;
1182 res_aabb[i0][i1][2] = -10000000;
1183 res_aabb[i0][i1][3] = -10000000;
1184 for(unsigned int i2 = 0; i2 < polylines[i0][i1].size(); i2++)
1185 {
1186 if( indicesFlag )
1187 {
1188 indices[i0][i1][i2].resize(4);
1189 indices[i0][i1][i2][0] = polylines[i0][i1][i2][0];
1190 indices[i0][i1][i2][1] = polylines[i0][i1][i2][1];
1191 indices[i0][i1][i2][2] = polylines[i0][i1][i2][2];
1192 indices[i0][i1][i2][3] = polylines[i0][i1][i2][3];
1193
1194 }
1195 else
1196 {
1197 params[i0][i1][i2].resize(4);
1198 // We could as well successively push_back() them but this should be slightly more efficient.
1199 params[i0][i1][i2][0] = (x_dir[polylines[i0][i1][i2][0]]);
1200 params[i0][i1][i2][1] = (y_dir[polylines[i0][i1][i2][1]]);
1201 params[i0][i1][i2][2] = (x_dir[polylines[i0][i1][i2][2]]);
1202 params[i0][i1][i2][3] = (y_dir[polylines[i0][i1][i2][3]]);
1203 }
1204 if(res_aabb[i0][i1][0]>signed(polylines[i0][i1][i2][0]))
1205 {
1206 res_aabb[i0][i1][0] = polylines[i0][i1][i2][0];
1207 }
1208 if(res_aabb[i0][i1][1]>signed(polylines[i0][i1][i2][1]))
1209 {
1210 res_aabb[i0][i1][1] = polylines[i0][i1][i2][1];
1211 }
1212 if(res_aabb[i0][i1][2]<signed(polylines[i0][i1][i2][2]))
1213 {
1214 res_aabb[i0][i1][2] = polylines[i0][i1][i2][2];
1215 }
1216 if(res_aabb[i0][i1][3]<signed(polylines[i0][i1][i2][3]))
1217 {
1218 res_aabb[i0][i1][3] = polylines[i0][i1][i2][3];
1219 }
1220 }
1221 }
1222 }
1223
1224 res_aabb_unsigned.resize(res_aabb.size());
1225 for (unsigned int i = 0; i < res_aabb.size(); i++)
1226 {
1227 res_aabb_unsigned[i].resize(res_aabb[i].size());
1228 for (unsigned int j = 0; j < res_aabb[i].size(); j++)
1229 {
1230 res_aabb_unsigned[i][j].resize(res_aabb[i][j].size());
1231 for (unsigned int k = 0; k < res_aabb[i][j].size(); k++)
1232 {
1233 if(res_aabb[i][j][k]<0)
1234 gsWarn << "conversion form signed to unsigned\n";
1235 res_aabb_unsigned[i][j][k] = res_aabb[i][j][k];
1236 }
1237 }
1238 }
1239 return res_aabb_unsigned;
1240 }
1241
1242
1243 template<short_t d, class T>
transfer(const std::vector<gsSortedVector<index_t>> & old,gsSparseMatrix<T> & result)1244 void gsHTensorBasis<d,T>::transfer(const std::vector<gsSortedVector<index_t> >& old, gsSparseMatrix<T>& result)
1245 {
1246 // Note: implementation assumes number of old + 1 m_bases exists in this basis
1247 needLevel( old.size() );
1248
1249 tensorBasis T_0_copy = this->tensorLevel(0);
1250
1251 std::vector< gsSparseMatrix<T,RowMajor> > transfer(m_bases.size()-1);
1252 std::vector<std::vector<T> > knots(d);
1253
1254 for(size_t i = 1; i < m_bases.size(); ++i)
1255 {
1256 //T_0_copy.uniformRefine_withTransfer(transfer[i-1], 1);
1257 for(unsigned dim = 0; dim != d; ++dim)
1258 {
1259 const gsKnotVector<T> & ckv = m_bases[i-1]->knots(dim);
1260 const gsKnotVector<T> & fkv = m_bases[i ]->knots(dim);
1261 ckv.symDifference(fkv, knots[dim]);
1262 // equivalent (dyadic ref.):
1263 // ckv.getUniformRefinementKnots(1, knots[dim]);
1264
1265 //gsDebug << "level: " << i << "\n"
1266 // << "direction: " << dim << "\n";
1267 //gsDebugVar(gsAsMatrix<T>(dirKnots));
1268 }
1269
1270 T_0_copy.refine_withTransfer(transfer[i-1], knots);
1271 }
1272
1273 // Add missing empty char. matrices
1274 while ( old.size() >= m_xmatrix.size() )
1275 m_xmatrix.push_back( gsSortedVector<index_t>() );
1276
1277 result = this->coarsening_direct(old, m_xmatrix, transfer);
1278
1279 // This function automatically adds additional characteristic matrices,
1280 // even if they are not needed.
1281 // Check whether the characteristic matrices corresponding to the finest
1282 // levels are actually used. If they are empty, i.e., if there are no
1283 // active functions on that level, drop them...
1284 while( m_xmatrix.back().size() == 0 )
1285 m_xmatrix.pop_back();
1286
1287 // ...similarly, erase all those fine bases which are actually not used.
1288 const int sizeDiff = static_cast<int>( m_bases.size() - m_xmatrix.size() );
1289 if( sizeDiff > 0 )
1290 {
1291 freeAll(m_bases.end() - sizeDiff, m_bases.end());
1292 m_bases.resize(m_xmatrix.size());
1293 }
1294 }
1295
1296 template<short_t d, class T>
transfer2(const std::vector<gsSortedVector<index_t>> & old,gsSparseMatrix<T> & result)1297 void gsHTensorBasis<d,T>::transfer2(const std::vector<gsSortedVector<index_t> >& old, gsSparseMatrix<T>& result)
1298 {
1299 // Note: implementation assumes number of old + 1 m_bases exists in this basis
1300 needLevel( old.size() );
1301
1302 tensorBasis T_0_copy = this->tensorLevel(0);
1303 std::vector< gsSparseMatrix<T,RowMajor> > transfer( m_bases.size()-1 );
1304 std::vector<std::vector<T> > knots(d);
1305
1306 for(size_t i = 1; i < m_bases.size(); ++i)
1307 {
1308 //T_0_copy.uniformRefine_withTransfer(transfer[i], 1);
1309 for(short_t dim = 0; dim != d; ++dim)
1310 {
1311 const gsKnotVector<T> & ckv = m_bases[i-1]->knots(dim);
1312 const gsKnotVector<T> & fkv = m_bases[i ]->knots(dim);
1313 ckv.symDifference(fkv, knots[dim]);
1314 // equivalent (dyadic ref.):
1315 // ckv.getUniformRefinementKnots(1, knots[dim]);
1316
1317 //gsDebug << "level: " << i << "\n"
1318 // << "direction: " << dim << "\n";
1319 //gsDebugVar(gsAsMatrix<T>(dirKnots));
1320 }
1321 T_0_copy.refine_withTransfer(transfer[i-1], knots);
1322 }
1323
1324 // Add missing empty char. matrices
1325 while ( old.size() >= m_xmatrix.size())
1326 m_xmatrix.push_back( gsSortedVector<index_t>() );
1327
1328 result = this->coarsening_direct2(old, m_xmatrix, transfer);
1329 }
1330
1331
1332 template<short_t d, class T>
increaseMultiplicity(index_t lvl,int dir,T knotValue,int mult)1333 void gsHTensorBasis<d,T>::increaseMultiplicity(index_t lvl, int dir, T knotValue, int mult)
1334 {
1335 GISMO_ASSERT( static_cast<size_t>(lvl) < m_xmatrix.size(),
1336 "Requested level does not exist.\n");
1337
1338 if (m_bases[lvl]->knots(dir).has(knotValue))
1339 {
1340 for(unsigned int i =lvl;i < m_bases.size();i++)
1341 m_bases[i]->component(dir).insertKnot(knotValue,mult);
1342 }
1343 else
1344 {
1345 gsWarn<<"Knot value not in the given knot vector."<<std::endl;
1346 }
1347
1348 update_structure();
1349 }
1350
1351
1352 template<short_t d, class T>
increaseMultiplicity(index_t lvl,int dir,const std::vector<T> & knotValue,int mult)1353 void gsHTensorBasis<d,T>::increaseMultiplicity(index_t lvl, int dir, const std::vector<T> & knotValue, int mult)
1354 {
1355 for(unsigned int k =0; k < knotValue.size(); ++k)
1356 {
1357 if (m_bases[lvl]->knots(dir).has(knotValue[k]))
1358 {
1359 for(unsigned int i =lvl;i < m_bases.size(); ++i)
1360 m_bases[i]->component(dir).insertKnot(knotValue[k],mult);
1361 }
1362 else
1363 {
1364 gsWarn<<"Knot value not in the given knot vector."<<std::endl;
1365 }
1366 }
1367 update_structure();
1368 }
1369
1370 template<short_t d, class T>
getBoxesAlongSlice(int dir,T par,std::vector<index_t> & boxes) const1371 void gsHTensorBasis<d,T>::getBoxesAlongSlice( int dir, T par,std::vector<index_t>& boxes ) const
1372 {
1373 gsMatrix<index_t> b1,b2;
1374 gsVector<index_t> level;
1375 m_tree.getBoxesInLevelIndex(b1,b2,level);
1376 gsVector<index_t> min,max;
1377 for(index_t i = 0;i<level.rows();i++)
1378 {
1379 min = b1.row(i);
1380 max = b2.row(i);
1381 const index_t l = level(i);
1382 const index_t par_index = m_bases[l]->knots(dir).uFind(par).uIndex();
1383 if(l>0 && (par_index>=min(dir)) && (par_index<=max(dir)))
1384 {
1385 boxes.push_back(l);
1386 for(index_t j=0;j<min.rows();++j)
1387 if(j!=dir)
1388 boxes.push_back(min(j));
1389 for(index_t j=0;j<max.rows();++j)
1390 if(j!=dir)
1391 boxes.push_back(max(j));
1392 }
1393 }
1394 }
1395
1396 template<short_t d, class T>
degreeElevate(int const & i,int const dir)1397 void gsHTensorBasis<d,T>::degreeElevate(int const & i, int const dir)
1398 {
1399 for (size_t level=0;level<m_bases.size();++level)
1400 m_bases[level]->degreeElevate(i,dir);
1401
1402 for(unsigned c=0; c<d; ++c)
1403 m_deg[c]=m_bases[0]->degree(c);
1404
1405 this->update_structure();
1406 }
1407
1408 template<short_t d, class T>
degreeReduce(int const & i,int const dir)1409 void gsHTensorBasis<d,T>::degreeReduce(int const & i, int const dir)
1410 {
1411 for (size_t level=0;level<m_bases.size();++level)
1412 m_bases[level]->degreeReduce(i,dir);
1413
1414 for(unsigned c=0; c<d; ++c)
1415 m_deg[c]=m_bases[0]->degree(c);
1416
1417 this->update_structure();
1418 }
1419
1420 template<short_t d, class T>
degreeIncrease(int const & i,int const dir)1421 void gsHTensorBasis<d,T>::degreeIncrease(int const & i, int const dir)
1422 {
1423 for (size_t level=0;level<m_bases.size();++level)
1424 m_bases[level]->degreeIncrease(i,dir);
1425
1426 for(unsigned c=0; c<d; ++c)
1427 m_deg[c]=m_bases[0]->degree(c);
1428
1429 this->update_structure();
1430 }
1431
1432 template<short_t d, class T>
degreeDecrease(int const & i,int const dir)1433 void gsHTensorBasis<d,T>::degreeDecrease(int const & i, int const dir)
1434 {
1435 for (size_t level=0;level<m_bases.size();++level)
1436 m_bases[level]->degreeDecrease(i,dir);
1437
1438 for(unsigned c=0; c<d; ++c)
1439 m_deg[c]=m_bases[0]->degree(c);
1440
1441 this->update_structure();
1442 }
1443
1444 namespace internal
1445 {
1446
1447 /// Get a HTensorBasis from XML data
1448 template<short_t d, class T>
1449 class gsXml< gsHTensorBasis<d,T> >
1450 {
1451 private:
gsXml()1452 gsXml() { }
1453 public:
1454 GSXML_COMMON_FUNCTIONS(gsHTensorBasis<TMPLA2(d,T)>);
tag()1455 static std::string tag () { return "Basis"; }
type()1456 static std::string type () { return ""; } // tag ?
1457
get(gsXmlNode * node)1458 static gsHTensorBasis<d,T> * get (gsXmlNode * node)
1459 {
1460 gsXmlAttribute * btype = node->first_attribute("type");
1461 if ( ! btype )
1462 {
1463 gsWarn<< "Basis without a type in the xml file.\n";
1464 return NULL;
1465 }
1466 std::string s = btype->value() ;
1467 if ( s.compare(0, 9, "HBSplineB" , 9 ) == 0 ) // needs correct d as well
1468 return gsXml< gsHBSplineBasis<d,T> >::get(node);
1469 if ( s.compare(0, 10,"THBSplineB", 10) == 0 )
1470 return gsXml< gsTHBSplineBasis<d,T> >::get(node);
1471
1472 gsWarn<<"gsXmlUtils: gsHTensorBasis: No known basis \""<<s<<"\". Error.\n";
1473 return NULL;
1474 }
1475
put(const gsHTensorBasis<d,T> & obj,gsXmlTree & data)1476 static gsXmlNode * put (const gsHTensorBasis<d,T> & obj,
1477 gsXmlTree & data )
1478 {
1479 const gsBasis<T> * ptr = & obj;
1480
1481 // Hier. B-splines
1482 if ( const gsHBSplineBasis<d,T> * g =
1483 dynamic_cast<const gsHBSplineBasis<d,T> *>( ptr ) )
1484 return gsXml< gsHBSplineBasis<d,T> >::put(*g,data);
1485
1486 // Truncated hier. B-splines
1487 if ( const gsTHBSplineBasis<d,T> * g =
1488 dynamic_cast<const gsTHBSplineBasis<d,T> *>( ptr ) )
1489 return gsXml< gsTHBSplineBasis<d,T> >::put(*g,data);
1490
1491 gsWarn<<"gsXmlUtils put: getBasis: No known basis \""<<obj<<"\". Error.\n";
1492 return NULL;
1493 }
1494 };
1495
1496 } // namespace internal
1497
1498
1499 } // namespace gismo
1500