G+Smo  24.08.0
Geometry + Simulation Modules
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
gsMultiBasis.hpp
Go to the documentation of this file.
1 
14 #pragma once
15 
16 #include <gsCore/gsMultiPatch.h>
19 #include <gsIO/gsOptionList.h>
20 
21 namespace gismo
22 {
23 
24 template<class T>
25 gsMultiBasis<T>::gsMultiBasis( const gsBasis<T> & bb, bool NoRational)
26 : m_topology( bb.dim() )
27 {
28  if (NoRational)
29  m_bases.push_back( bb.source().clone().release() );
30  else
31  m_bases.push_back( bb.clone().release() );
32  m_topology.addBox();
33  m_topology.addAutoBoundaries();
34 }
35 
36 template<class T>
37 gsMultiBasis<T>::gsMultiBasis( const gsMultiPatch<T> & mpatch, bool NoRational)
38 : m_topology( mpatch.topology() )
39 {
40  m_bases = mpatch.basesCopy(NoRational);
41 }
42 
43 template<class T>
45 : Base(other),
46  m_bases ( other.m_bases.size() ),
47  m_topology ( other.m_topology )
48 {
49  cloneAll( other.m_bases.begin(), other.m_bases.end(), this->m_bases.begin() );
50 }
51 
52 #if EIGEN_HAS_RVALUE_REFERENCES
53 template<class T>
55 {
56  freeAll(m_bases);
57  m_bases.resize( other.m_bases.size() );
58  cloneAll( other.m_bases.begin(), other.m_bases.end(), this->m_bases.begin() );
59  m_topology = other.m_topology;
60  return *this;
61 }
62 #endif
63 
64 template<class T>
66 {
67  freeAll(m_bases);
68 }
69 
70 template<class T>
71 std::ostream& gsMultiBasis<T>::print( std::ostream& os ) const
72 {
73  os<<"Topology: "<< m_topology <<"\n";
74  return os;
75 }
76 
77 
78 template<class T>
80 {
81  //gsDebug<< "TO DO\n";
82  if ( m_topology.dim() == -1 )
83  {
84  m_topology.setDim( g->dim() );
85  }
86  else
87  {
88  assert( g->dim() == m_topology.dim() );
89  }
90  m_bases.push_back( g ) ;
91  m_topology.addBox();
92  g = NULL;
93 }
94 
95 template<class T>
97 {
98  if ( m_topology.dim() == -1 )
99  {
100  m_topology.setDim( g->dim() );
101  }
102  else
103  {
104  GISMO_ASSERT( g->dim() == m_topology.dim(), "Dimensions do not match.");
105  }
106 
107  m_bases.push_back( g.release() ) ;
108  m_topology.addBox();
109 }
110 
111 template<class T>
113 {
114  typename BasisContainer::const_iterator it
115  = std::find( m_bases.begin(), m_bases.end(), g );
116  assert( it != m_bases.end() );
117  return it - m_bases.begin();
118 }
119 
120 template<class T>
122  gsBasis<T>* g2, boxSide s2 )
123 {
124  int p1 = findBasisIndex( g1 );
125  int p2 = findBasisIndex( g2 );
126  m_topology.addInterface( p1, s1, p2, s2);
127 }
128 
129 namespace {
130 template <typename T>
131 struct take_first {
132  T operator() (const T& a, const T&) { return a; }
133 };
134 }
135 
136 template <typename T>
138  const std::vector< gsSparseMatrix<T, RowMajor> >& localTransferMatrices,
139  const gsDofMapper& coarseMapper,
140  const gsDofMapper& fineMapper,
141  gsSparseMatrix<T, RowMajor>& transferMatrix
142  )
143 {
144  const index_t nBases = localTransferMatrices.size();
145 
146  index_t nonzeros = 0;
147  index_t trRows = 0;
148  index_t trCols = 0;
149 
150  for (index_t j=0; j<nBases; ++j)
151  {
152  nonzeros += localTransferMatrices[j].nonZeros();
153  trRows += localTransferMatrices[j].rows();
154  trCols += localTransferMatrices[j].cols();
155  }
156 
157  gsSparseEntries<T> entries;
158  entries.reserve( nonzeros );
159 
160  for (index_t j=0; j<nBases; ++j)
161  {
162  for (index_t k=0; k < localTransferMatrices[j].outerSize(); ++k)
163  {
164  for (typename gsSparseMatrix<T, RowMajor>::iterator it(localTransferMatrices[j],k); it; ++it)
165  {
166  const index_t coarse_dof_idx = coarseMapper.index(it.col(),j);
167  const index_t fine_dof_idx = fineMapper.index(it.row(),j);
168 
169  if (coarseMapper.is_free_index(coarse_dof_idx) && fineMapper.is_free_index(fine_dof_idx))
170  entries.add(fine_dof_idx, coarse_dof_idx, it.value());
171  }
172  }
173  }
174 
175  transferMatrix.resize(fineMapper.freeSize(), coarseMapper.freeSize());
176  transferMatrix.setFromTriplets(entries.begin(), entries.end(), take_first<T>());
177  transferMatrix.makeCompressed();
178 }
179 
180 template <typename T>
182  gsSparseMatrix<T, RowMajor>& transferMatrix,
183  const gsBoundaryConditions<T>& boundaryConditions,
184  const gsOptionList& assemblerOptions,
185  int numKnots,
186  int mul,
187  index_t unk)
188 {
189  // Get coarse mapper
190  gsDofMapper coarseMapper;
191  this->getMapper(
192  (dirichlet::strategy)assemblerOptions.askInt("DirichletStrategy",11),
193  (iFace ::strategy)assemblerOptions.askInt("InterfaceStrategy", 1),
194  boundaryConditions,
195  coarseMapper,
196  unk
197  );
198 
199  // Refine
200  std::vector< gsSparseMatrix<T, RowMajor> > localTransferMatrices(nBases());
201  for (size_t k = 0; k < m_bases.size(); ++k)
202  {
203  m_bases[k]->uniformRefine_withTransfer(localTransferMatrices[k],numKnots,mul);
204  }
205 
206  // Get fine mapper
207  gsDofMapper fineMapper;
208  this->getMapper(
209  (dirichlet::strategy)assemblerOptions.askInt("DirichletStrategy",11),
210  (iFace ::strategy)assemblerOptions.askInt("InterfaceStrategy", 1),
211  boundaryConditions,
212  fineMapper,
213  unk
214  );
215 
216  // restrict to free dofs
217  combineTransferMatrices( localTransferMatrices, coarseMapper, fineMapper, transferMatrix );
218 
219 }
220 
221 template <typename T>
223  gsSparseMatrix<T, RowMajor>& transferMatrix,
224  const gsBoundaryConditions<T>& boundaryConditions,
225  const gsOptionList& assemblerOptions,
226  int numKnots,
227  index_t unk)
228 {
229  // Get fine mapper
230  gsDofMapper fineMapper;
231  this->getMapper(
232  (dirichlet::strategy)assemblerOptions.askInt("DirichletStrategy",11),
233  (iFace ::strategy)assemblerOptions.askInt("InterfaceStrategy", 1),
234  boundaryConditions,
235  fineMapper,
236  unk
237  );
238 
239  // Refine
240  std::vector< gsSparseMatrix<T, RowMajor> > localTransferMatrices(nBases());
241  for (size_t k = 0; k < m_bases.size(); ++k)
242  {
243  m_bases[k]->uniformCoarsen_withTransfer(localTransferMatrices[k],numKnots);
244  }
245 
246  // Get coarse mapper
247  gsDofMapper coarseMapper;
248  this->getMapper(
249  (dirichlet::strategy)assemblerOptions.askInt("DirichletStrategy",11),
250  (iFace ::strategy)assemblerOptions.askInt("InterfaceStrategy", 1),
251  boundaryConditions,
252  coarseMapper,
253  unk
254  );
255 
256  // restrict to free dofs
257  combineTransferMatrices( localTransferMatrices, coarseMapper, fineMapper, transferMatrix );
258 
259 }
260 
261 template <typename T>
263  patchComponent pc,
264  const gsDofMapper& dm,
265  gsMatrix<index_t>& indices,
266  bool no_lower
267  ) const
268 {
269  typename gsBasis<T>::uPtr result = m_bases[pc.patch()]->componentBasis_withIndices(pc, indices, no_lower);
270 
271  const index_t sz = indices.rows();
272  index_t j = 0;
273  for (index_t i=0; i<sz; ++i)
274  {
275  const index_t loc = indices(i,0);
276  if (dm.is_free(loc, pc.patch()))
277  {
278  indices(j,0) = dm.index(loc, pc.patch());
279  ++j;
280  }
281  }
282  indices.conservativeResize(j,1);
283  return result;
284 }
285 
286 template <typename T>
287 std::vector<typename gsBasis<T>::uPtr> gsMultiBasis<T>::componentBasis_withIndices(
288  const std::vector<patchComponent>& pc,
289  const gsDofMapper& dm,
290  gsMatrix<index_t>& indices,
291  bool no_lower
292  ) const
293 {
294  const index_t nrpc = pc.size();
295  std::vector<typename gsBasis<T>::uPtr> bases;
296  bases.reserve(nrpc);
297  std::vector< gsMatrix<index_t> > local_indices(nrpc);
298 
299  index_t sz = 0;
300 
301  for (index_t i=0; i<nrpc; ++i)
302  {
303  bases.push_back(
304  this->componentBasis_withIndices(pc[i], dm, local_indices[i], no_lower)
305  );
306  sz += local_indices[i].rows();
307  }
308 
309  std::vector<index_t> global_indices;
310  global_indices.reserve(sz);
311 
312  for (index_t i=0; i<nrpc; ++i)
313  {
314  const index_t nr_local_indices = local_indices[i].rows();
315  for (index_t j=0; j<nr_local_indices; ++j)
316  if ( i==0 || std::find(global_indices.begin(), global_indices.end(), local_indices[i](j,0) )
317  == global_indices.end()
318  )
319  global_indices.push_back( local_indices[i](j,0) );
320  }
321 
322  const index_t final_size = global_indices.size();
323  indices.resize(final_size,1);
324  for (index_t i=0; i<final_size; ++i)
325  indices(i,0) = global_indices[i];
326 
327  return bases;
328 }
329 
330 template<class T>
332 {
333  GISMO_ASSERT(m_bases.size(), "Empty multibasis.");
334  short_t result = m_bases[0]->degree(k);
335  for (size_t i = 0; i < m_bases.size(); ++i)
336  if (m_bases[i]->degree(k) > result )
337  result = m_bases[i]->degree(k);
338  return result;
339 }
340 
341 template<class T>
343 {
344  GISMO_ASSERT(m_bases.size(), "Empty multibasis.");
345  short_t result = m_bases[0]->maxDegree();
346  for (size_t i = 0; i < m_bases.size(); ++i)
347  result = math::max(m_bases[i]->maxDegree(), result);
348  return result;
349 }
350 
351 template<class T>
353 {
354  GISMO_ASSERT(m_bases.size(), "Empty multibasis.");
355  short_t result = m_bases[0]->minDegree();
356  for (size_t i = 0; i < m_bases.size(); ++i)
357  result = math::min(m_bases[i]->minDegree(), result);
358  return result;
359 }
360 
361 template<class T>
363 {
364  GISMO_ASSERT(m_bases.size(), "Empty multibasis.");
365  short_t result = m_bases[0]->degree(k);
366  for (size_t i = 0; i < m_bases.size(); ++i)
367  if (m_bases[i]->degree(k) < result )
368  result = m_bases[i]->degree(k);
369  return result;
370 }
371 
372 template<class T>
373 void gsMultiBasis<T>::getMapper(bool conforming,
374  gsDofMapper & mapper,
375  bool finalize) const
376 {
377  mapper = gsDofMapper(*this);//.init(*this);
378 
379  if ( conforming ) // Conforming boundaries ?
380  {
381  for ( gsBoxTopology::const_iiterator it = m_topology.iBegin();
382  it != m_topology.iEnd(); ++it )
383  {
384  matchInterface(*it,mapper);
385  }
386  }
387 
388  if (finalize)
389  mapper.finalize();
390 }
391 
392 template<class T>
393 void gsMultiBasis<T>::getMapper(bool conforming,
394  const gsBoundaryConditions<T> & bc,
395  int unk,
396  gsDofMapper & mapper,
397  bool finalize) const
398 {
399  mapper = gsDofMapper(*this, bc, unk); //.init(*this, bc, unk);
400 
401  if ( conforming ) // Conforming boundaries ?
402  {
403  for ( gsBoxTopology::const_iiterator it = m_topology.iBegin();
404  it != m_topology.iEnd(); ++it )
405  {
406  matchInterface(*it,mapper);
407  }
408  }
409 
410  if (finalize)
411  mapper.finalize();
412 }
413 
414 template<class T>
416 {
417  // should work for all basis which have matchWith() implementeds
418  gsMatrix<index_t> b1, b2;
419  m_bases[bi.first().patch]->matchWith(bi, *m_bases[bi.second().patch],
420  b1, b2);
421 
422  // Match the dofs on the interface
423  for (size_t i = 0; i!=mapper.componentsSize(); ++i)
424  mapper.matchDofs(bi.first().patch, b1, bi.second().patch, b2, i );
425 }
426 
427 template<class T>
429 {
430  bool changed = false;
431 
432  std::vector<index_t> refEltsFirst;
433  std::vector<index_t> refEltsSecond;
434 
435  // Find the areas/elements that do not match...
436  switch( this->dim() )
437  {
438  case 2:
439  changed = this->template repairInterfaceFindElements<2>( bi, refEltsFirst, refEltsSecond );
440  break;
441  case 3:
442  changed = this->template repairInterfaceFindElements<3>( bi, refEltsFirst, refEltsSecond );
443  break;
444  default:
445  GISMO_ASSERT(false,"wrong dimension");
446  }
447 
448  // ...and if there are any found, refine the bases accordingly
449  if( changed )
450  {
451  if( refEltsFirst.size() > 0 )
452  m_bases[ bi.first().patch ]->refineElements( refEltsFirst );
453  if( refEltsSecond.size() > 0 )
454  m_bases[ bi.second().patch ]->refineElements( refEltsSecond );
455  }
456 
457  return changed;
458 }
459 
460 template<class T>
461 template<short_t d>
463  const boundaryInterface & bi,
464  std::vector<index_t> & refEltsFirst,
465  std::vector<index_t> & refEltsSecond )
466 {
467  GISMO_ASSERT( d == 2 || d == 3, "Dimension must be 2 or 3.");
468 
469  refEltsFirst.clear();
470  refEltsSecond.clear();
471 
472  // get direction and orientation maps
473  const gsVector<bool> dirOrient = bi.dirOrientation();
474  const gsVector<index_t> dirMap= bi.dirMap();
475 
476  // get the bases of both sides as gsHTensorBasis
477  const gsHTensorBasis<d,T> * bas0 = dynamic_cast< const gsHTensorBasis<d,T> * >( m_bases[ bi.first().patch ] );
478  const gsHTensorBasis<d,T> * bas1 = dynamic_cast< const gsHTensorBasis<d,T> * >( m_bases[ bi.second().patch ] );
479 
480  GISMO_ASSERT( bas0 != 0 && bas1 != 0, "Cannot cast basis as needed.");
481 
482  gsMatrix<index_t> lo0;
483  gsMatrix<index_t> up0;
484  gsVector<index_t> level0;
485  gsMatrix<index_t> lo1;
486  gsMatrix<index_t> up1;
487  gsVector<index_t> level1;
488 
489  unsigned idxExponent;
490 
491  // get the higher one of both indexLevels
492  unsigned indexLevelUse = ( bas0->tree().getIndexLevel() > bas1->tree().getIndexLevel() ? bas0->tree().getIndexLevel() : bas1->tree().getIndexLevel() );
493  unsigned indexLevelDiff0 = indexLevelUse - bas0->tree().getIndexLevel();
494  unsigned indexLevelDiff1 = indexLevelUse - bas1->tree().getIndexLevel();
495 
496  // get upper corners, but w.r.t. level "indexLevelUse"
497  gsVector<index_t> upperCorn0 = bas0->tree().upperCorner();
498  gsVector<index_t> upperCorn1 = bas1->tree().upperCorner();
499  for( short_t i=0; i < d; i++)
500  {
501  upperCorn0[i] = upperCorn0[i] << indexLevelDiff0;
502  upperCorn1[i] = upperCorn1[i] << indexLevelDiff1;
503  }
504 
505 // GISMO_ASSERT( upperCorn0[0] == upperCorn1[0] &&
506 // upperCorn0[1] == upperCorn1[1], "The meshes are not matching as they should be!");
507 // GISMO_ASSERT( (d<3) || (upperCorn0[2] == upperCorn1[2]),
508 // "The meshes are not matching as they should be!");
509 
510  // get the box-representation of the gsHDomain on the interface
511  bas0->tree().getBoxesOnSide( bi.first().side(), lo0, up0, level0);
512  bas1->tree().getBoxesOnSide( bi.second().side(), lo1, up1, level1);
513 
514  // Compute the indices on the same level (indexLevelUse)
515  idxExponent = ( indexLevelUse - bas0->tree().getMaxInsLevel());
516  for( index_t i=0; i < lo0.rows(); i++)
517  for( short_t j=0; j < d; j++)
518  {
519  lo0(i,j) = lo0(i,j) << idxExponent;
520  up0(i,j) = up0(i,j) << idxExponent;
521  }
522  idxExponent = ( indexLevelUse - bas1->tree().getMaxInsLevel());
523  for( index_t i=0; i < lo1.rows(); i++)
524  for( short_t jj=0; jj < d; jj++)
525  {
526  // Computation done via dirMap, because...
527  unsigned j = dirMap[jj];
528  lo1(i,j) = lo1(i,j) << idxExponent;
529  up1(i,j) = up1(i,j) << idxExponent;
530 
531  //... we also have to check whether the orientation
532  // is preserved or not.
533  if( !dirOrient[jj] )
534  {
535  unsigned tmp = upperCorn1[j] - lo1(i,j);
536  lo1(i,j) = upperCorn1[j] - up1(i,j);
537  up1(i,j) = tmp;
538  }
539  }
540 
541  // Find the merged interface mesh with
542  // the respective levels.
543  // Not efficient, but simple to implement.
544 
545  // a, b will correspond to the coordinate
546  // directions which "span" the interface
547  unsigned a0, b0, a1, b1;
548  // c corresponds to the coordinate direction which
549  // defines the interface-side by being set to 0 or 1
550  unsigned c0, c1;
551  switch( bi.first().direction() )
552  {
553  case 0:
554  a0 = 1; b0 = 2; c0 = 0;
555  break;
556  case 1:
557  a0 = 0; b0 = 2; c0 = 1;
558  break;
559  case 2:
560  a0 = 0; b0 = 1; c0 = 2;
561  break;
562  default:
563  a0 = 99; b0 = 99; c0 = 99;
564  // Initialize for CDash, let it crash if this happens.
565  }
566 
567  // If d == 2, the "b"'s are not needed.
568  // Setting them to the "a"'s will
569  // result in some steps and tests being repeated,
570  // but the implementation not optimized w.r.t. efficiency.
571  if( d == 2 )
572  b0 = a0;
573 
574  a1 = dirMap[a0];
575  b1 = dirMap[b0];
576  c1 = dirMap[c0];
577 
578  // Run through all possible pairings of
579  // boxes and see if they overlap.
580  // If so, their overlap is a box of the merged
581  // interface mesh.
582  std::vector< std::vector<unsigned> > iU;
583  for( index_t i0 = 0; i0 < lo0.rows(); i0++)
584  for( index_t i1 = 0; i1 < lo1.rows(); i1++)
585  {
586  if( lo0(i0,a0) < up1(i1,a1) &&
587  lo0(i0,b0) < up1(i1,b1) &&
588  lo1(i1,a1) < up0(i0,a0) &&
589  lo1(i1,b1) < up0(i0,b0) )// overlap
590  {
591  std::vector<unsigned> tmp;
592  tmp.push_back( std::max( lo0(i0,a0), lo1(i1,a1) ) );
593  tmp.push_back( std::max( lo0(i0,b0), lo1(i1,b1) ) ); // duplicate in 2D
594  tmp.push_back( std::min( up0(i0,a0), up1(i1,a1) ) );
595  tmp.push_back( std::min( up0(i0,b0), up1(i1,b1) ) ); // duplicate in 2D
596  tmp.push_back( level0[i0] );
597  tmp.push_back( level1[i1] );
598  iU.push_back( tmp );
599  }
600  }
601 
602  std::vector<unsigned> tmpvec(1+2*d);
603  for( size_t i = 0; i < size_t( iU.size() ); i++)
604  {
605  // the levels on both sides of the
606  // box of the interface
607  unsigned L0 = iU[i][4];
608  unsigned L1 = iU[i][5];
609 
610  if( L0 != L1 ) // one side has to be refined
611  {
612  unsigned a, b, c, Luse;
613  unsigned refSideIndex;
614  unsigned upperCornOnLevel;
615 
616  if( L0 < L1 ) // refine first()
617  {
618  Luse = L1;
619  a = a0;
620  b = b0;
621  c = c0;
622  refSideIndex = bi.first().side().index();
623  upperCornOnLevel = ( upperCorn0[c] >> ( indexLevelUse - Luse ) );
624  }
625  else // refine second()
626  {
627  Luse = L0;
628  a = a1;
629  b = b1;
630  c = c1;
631  refSideIndex = bi.second().side().index();
632  upperCornOnLevel = ( upperCorn1[c] >> ( indexLevelUse - Luse ) );
633  }
634 
635  // store the new level
636  tmpvec[0] = Luse;
637  // store the box on the interface that
638  // has to be refined to that new level
639  tmpvec[1+a] = iU[i][0] >> ( indexLevelUse - Luse );
640  tmpvec[1+d+a] = iU[i][2] >> ( indexLevelUse - Luse );
641  if( d == 3 )
642  {
643  tmpvec[1+b] = iU[i][1] >> ( indexLevelUse - Luse );
644  tmpvec[1+d+b] = iU[i][3] >> ( indexLevelUse - Luse );
645  }
646 
647  if( refSideIndex % 2 == 1 )
648  {
649  // west, south, front:
650  tmpvec[1+c] = 0;
651  tmpvec[1+d+c] = 1;
652  }
653  else
654  {
655  // east, north, back:
656  tmpvec[1+c] = upperCornOnLevel-1;
657  tmpvec[1+d+c] = upperCornOnLevel;
658  }
659 
660  if( Luse == L1 ) // refine first
661  {
662  // no messing around with orientation and
663  // maps needed.
664  for( unsigned j=0; j < tmpvec.size(); j++)
665  refEltsFirst.push_back( tmpvec[j] );
666  }
667  else // refine second
668  {
669  // if the orientation is changed, flip where necessary
670  for( unsigned jj = 0; jj < d; jj++)
671  {
672  unsigned j = dirMap[jj];
673  if( j != c && !dirOrient[ jj ] )
674  {
675  upperCornOnLevel = ( upperCorn1[j] >> ( indexLevelUse - Luse ) );
676  unsigned tmp = tmpvec[1+j];
677  tmpvec[1+j] = upperCornOnLevel - tmpvec[1+d+j];
678  tmpvec[1+d+j] = upperCornOnLevel - tmp;
679  }
680  }
681 
682  for( unsigned j=0; j < (1+2*d); j++)
683  refEltsSecond.push_back( tmpvec[j] );
684  }
685  }
686  }
687 
688  return ( ( refEltsFirst.size() > 0 ) || ( refEltsSecond.size() > 0 ) );
689 }
690 
691 template<class T>
693 {
694  // get direction and orientation maps
695  const gsVector<bool> dirOrient = bi.dirOrientation();
696 
697  // get the bases of both sides as gsHTensorBasis
698  const gsHTensorBasis<2,T> * bas0 = dynamic_cast< const gsHTensorBasis<2,T> * >( m_bases[ bi.first().patch ] );
699  const gsHTensorBasis<2,T> * bas1 = dynamic_cast< const gsHTensorBasis<2,T> * >( m_bases[ bi.second().patch ] );
700 
701  GISMO_ASSERT( bas0 != 0 && bas1 != 0, "Cannot cast basis as needed.");
702 
705  gsVector<index_t> level;
706 
707  unsigned idxExponent;
708 
709  // get the higher one of both indexLevels
710  unsigned indexLevelUse = ( bas0->tree().getIndexLevel() > bas1->tree().getIndexLevel() ? bas0->tree().getIndexLevel() : bas1->tree().getIndexLevel() );
711  unsigned indexLevelDiff0 = indexLevelUse - bas0->tree().getIndexLevel();
712  unsigned indexLevelDiff1 = indexLevelUse - bas1->tree().getIndexLevel();
713 
714  // get the box-representation of the gsHDomain on the interface
715  bas0->tree().getBoxesOnSide( bi.first().side(), lo, up, level);
716 
717  int dir0 = ( bi.first().direction() + 1 ) % 2;
718  bool orientPreserv = dirOrient[ dir0 ];
719  // for mapping the indices to the same
720  idxExponent = ( indexLevelUse - bas0->tree().getMaxInsLevel());
721  gsMatrix<unsigned> intfc0( lo.rows(), 3 );
722  for( index_t i=0; i < lo.rows(); i++)
723  {
724  intfc0(i,0) = lo(i,dir0) << idxExponent;
725  intfc0(i,1) = up(i,dir0) << idxExponent;
726  intfc0(i,2) = level[i];
727  }
728  intfc0.sortByColumn(0);
729 
730  // get the box-representation of the gsHDomain on the interface
731  bas1->tree().getBoxesOnSide( bi.second().side(), lo, up, level);
732  int dir1 = ( bi.second().direction() + 1 ) % 2;
733  idxExponent = ( indexLevelUse - bas1->tree().getMaxInsLevel());
734  gsMatrix<unsigned> intfc1( lo.rows(), 3 );
735  for( index_t i=0; i < lo.rows(); i++)
736  {
737  intfc1(i,0) = lo(i,dir1) << idxExponent;
738  intfc1(i,1) = up(i,dir1) << idxExponent;
739  intfc1(i,2) = level[i];
740  }
741 
742  // now the knot indices in intfc0 and intfc1 both correspond to
743  // numbering on level "indexLevelUse"
744 
745  // get upper corners, but w.r.t. level "indexLevelUse"
746  gsVector<index_t,2> upperCorn0 = bas0->tree().upperCorner();
747  upperCorn0[0] = upperCorn0[0] << indexLevelDiff0;
748  upperCorn0[1] = upperCorn0[1] << indexLevelDiff0;
749 
750  gsVector<index_t,2> upperCorn1 = bas1->tree().upperCorner();
751  upperCorn1[0] = upperCorn1[0] << indexLevelDiff1;
752  upperCorn1[1] = upperCorn1[1] << indexLevelDiff1;
753 
754  if( !orientPreserv )
755  {
756  // flip the knot indices
757  for( index_t i=0; i < lo.rows(); i++)
758  {
759  const index_t tmp = upperCorn1[dir1] - intfc1(i, 1);
760  intfc1(i,1) = upperCorn1[dir1] - intfc1(i, 0);
761  intfc1(i,0) = tmp;
762  }
763  }
764  intfc1.sortByColumn(0);
765 
766  GISMO_ASSERT(intfc0( intfc0.rows()-1, 1) == intfc1( intfc1.rows()-1, 1)," Something wrong with interfaces! Mark 264");
767 
768  // Merge the knot spans from both sides into intfcU
769  // intfcU[i][0]: end-knot-index
770  // intfcU[i][1]: level on first()
771  // intfcU[i][2]: level on second()
772  int i0 = 0; int i1 = 0;
773  std::vector< std::vector< unsigned > > intfcU;
774  while( i0 < intfc0.rows() && i1 < intfc1.rows() )
775  {
776  std::vector<unsigned> tmp(3);
777 
778  if( intfc0( i0, 1 ) == intfc1( i1, 1 ) )
779  {
780  tmp[0] = intfc0(i0,1);
781  tmp[1] = intfc0(i0,2);
782  tmp[2] = intfc1(i1,2);
783  intfcU.push_back( tmp );
784  i0++;
785  i1++;
786  }
787  else if( intfc0( i0, 1 ) > intfc1( i1, 1 ) )
788  {
789  tmp[0] = intfc1(i1,1);
790  tmp[1] = intfc0(i0,2);
791  tmp[2] = intfc1(i1,2);
792  intfcU.push_back( tmp );
793  i1++;
794  }
795  else
796  {
797  tmp[0] = intfc0(i0,1);
798  tmp[1] = intfc0(i0,2);
799  tmp[2] = intfc1(i1,2);
800  intfcU.push_back( tmp );
801  i0++;
802  }
803  }
804 
805  // create the refineboxes needed for
806  // reparing the interface
807  unsigned knot0;
808  unsigned knot1 = 0;
809  std::vector<index_t> refElts0;
810  std::vector<index_t> refElts1;
811 
812  for( unsigned i=0; i < intfcU.size(); i++)
813  {
814  knot0 = knot1;
815  knot1 = intfcU[i][0];
816  unsigned L0 = intfcU[i][1];
817  unsigned L1 = intfcU[i][2];
818 
819  if( L0 < L1 ) // refine first()
820  {
821  refElts0.push_back( L1 );
822 
823  // knot indices on level L1:
824  unsigned knot0L = knot0 >> ( indexLevelUse - L1 );
825  unsigned knot1L = knot1 >> ( indexLevelUse - L1 );
826 
827  unsigned upperCornOnLevel;
828  switch( bi.first().side().index() )
829  {
830  case 1: // west
831  refElts0.push_back( 0 );
832  refElts0.push_back( knot0L );
833  refElts0.push_back( 1 );
834  refElts0.push_back( knot1L );
835  break;
836  case 2: // east
837  upperCornOnLevel = ( upperCorn0[0] >> ( indexLevelUse - L1 ) );
838 
839  refElts0.push_back( upperCornOnLevel-1 );
840  refElts0.push_back( knot0L );
841  refElts0.push_back( upperCornOnLevel );
842  refElts0.push_back( knot1L );
843  break;
844  case 3: // south
845  refElts0.push_back( knot0L );
846  refElts0.push_back( 0 );
847  refElts0.push_back( knot1L );
848  refElts0.push_back( 1 );
849  break;
850  case 4: // north
851  upperCornOnLevel = ( upperCorn0[1] >> ( indexLevelUse - L1 ) );
852 
853  refElts0.push_back( knot0L );
854  refElts0.push_back( upperCornOnLevel-1 );
855  refElts0.push_back( knot1L );
856  refElts0.push_back( upperCornOnLevel );
857  break;
858  default:
859  GISMO_ASSERT(false,"3D not implemented yet. You can do it, if you want.");
860  break;
861  }
862  }
863  else if( L0 > L1 ) // refine second()
864  {
865  refElts1.push_back( L0 );
866 
867  // knot indices on level "indexLevelUse":
868  unsigned knot0L = knot0;
869  unsigned knot1L = knot1;
870  // flip, if necessary
871  if( !orientPreserv )
872  {
873  unsigned tmp = knot0L;
874  knot0L = ( upperCorn1[dir1] - knot1L );
875  knot1L = ( upperCorn1[dir1] - tmp );
876  }
877  // push to level L0
878  knot0L = knot0L >> ( indexLevelUse - L0 );
879  knot1L = knot1L >> ( indexLevelUse - L0 );
880 
881  gsVector<unsigned> upperCornOnLevel(2);
882  switch( bi.second().side().index() )
883  {
884  case 1: // west
885  refElts1.push_back( 0 );
886  refElts1.push_back( knot0L );
887  refElts1.push_back( 1 );
888  refElts1.push_back( knot1L );
889  break;
890  case 2: // east
891  upperCornOnLevel[0] = ( upperCorn1[0] >> ( indexLevelUse - L0 ) );
892  upperCornOnLevel[1] = ( upperCorn1[1] >> ( indexLevelUse - L0 ) );
893 
894  refElts1.push_back( upperCornOnLevel[0]-1 );
895  refElts1.push_back( knot0L );
896  refElts1.push_back( upperCornOnLevel[0] );
897  refElts1.push_back( knot1L );
898  break;
899  case 3: // south
900  refElts1.push_back( knot0L );
901  refElts1.push_back( 0 );
902  refElts1.push_back( knot1L );
903  refElts1.push_back( 1 );
904  break;
905  case 4: // north
906  upperCornOnLevel[0] = ( upperCorn1[0] >> ( indexLevelUse - L0 ) );
907  upperCornOnLevel[1] = ( upperCorn1[1] >> ( indexLevelUse - L0 ) );
908 
909  refElts1.push_back( knot0L );
910  refElts1.push_back( upperCornOnLevel[1]-1 );
911  refElts1.push_back( knot1L );
912  refElts1.push_back( upperCornOnLevel[1] );
913  break;
914  default:
915  GISMO_ASSERT(false,"3D not implemented yet. You can do it, if you want.");
916  break;
917  }
918  }
919  }
920 
921  if( refElts0.size() > 0 )
922  m_bases[ bi.first().patch ]->refineElements( refElts0 );
923  if( refElts1.size() > 0 )
924  m_bases[ bi.second().patch ]->refineElements( refElts1 );
925 
926  return ( ( refElts0.size() > 0 ) || ( refElts1.size() > 0 ) );
927 
928 }
929 
930 } // namespace gismo
gsBasis< T >::uPtr componentBasis_withIndices(patchComponent pc, const gsDofMapper &dm, gsMatrix< index_t > &indices, bool no_lower=true) const
Returns the basis that corresponds to the component.
Definition: gsMultiBasis.hpp:262
Provides definition of HTensorBasis abstract interface.
Class that provides a container for triplets (i,j,value) to be filled in a sparse matrix...
Definition: gsSparseMatrix.h:33
int findBasisIndex(gsBasis< T > *g) const
Search for the given basis and return its index.
Definition: gsMultiBasis.hpp:112
void addInterface(gsBasis< T > *g1, boxSide s1, gsBasis< T > *g2, boxSide s2)
Add an interface joint between side s1 of geometry g1 side s2 of geometry g2.
Definition: gsMultiBasis.hpp:121
index_t patch() const
Returns the patch number.
Definition: gsBoundary.h:626
gsMultiBasis()
Default empty constructor.
Definition: gsMultiBasis.h:58
gsMultiBasis & operator=(gsMultiBasis other)
Assignment operator (uses copy-and-swap idiom)
Definition: gsMultiBasis.h:114
void matchDofs(index_t u, const gsMatrix< index_t > &b1, index_t v, const gsMatrix< index_t > &b2, index_t comp=0)
Couples dofs b1 of patch u with dofs b2 of patch v one by one such that they refer to the same global...
Definition: gsDofMapper.cpp:159
#define short_t
Definition: gsConfig.h:35
void addAutoBoundaries()
Make all patch sides which are not yet declared as interface or boundary to a boundary.
Definition: gsBoxTopology.cpp:75
short_t maxDegree(short_t k) const
Maximum degree with respect to variable k.
Definition: gsMultiBasis.hpp:331
bool repairInterface2d(const boundaryInterface &bi)
Checks if the 2D-interface is fully matching, and if not, repairs it.
Definition: gsMultiBasis.hpp:692
short_t minCwiseDegree() const
Minimum degree with respect to all variables.
Definition: gsMultiBasis.hpp:352
bool is_free(index_t i, index_t k=0, index_t c=0) const
Returns true if local dof i of patch k is not eliminated.
Definition: gsDofMapper.h:382
#define index_t
Definition: gsConfig.h:32
const gsHDomain< d > & tree() const
Returns a reference to m_tree.
Definition: gsHTensorBasis.h:601
void matchInterface(const boundaryInterface &bi, gsDofMapper &mapper) const
Matches the degrees of freedom along an interface.
Definition: gsMultiBasis.hpp:415
std::ostream & print(std::ostream &os) const
Prints the object as a string.
Definition: gsMultiBasis.hpp:71
Provides combinatorial unitilies.
#define GISMO_ASSERT(cond, message)
Definition: gsDebug.h:89
void cloneAll(It start, It end, ItOut out)
Clones all pointers in the range [start end) and stores new raw pointers in iterator out...
Definition: gsMemory.h:295
Maintains a mapping from patch-local dofs to global dof indices and allows the elimination of individ...
Definition: gsDofMapper.h:68
bool repairInterface(const boundaryInterface &bi)
Checks if the interface is fully matching, and if not, repairs it.
Definition: gsMultiBasis.hpp:428
void finalize()
Must be called after all boundaries and interfaces have been marked to set up the dof numbering...
Definition: gsDofMapper.cpp:240
index_t index(index_t i, index_t k=0, index_t c=0) const
Returns the global dof index associated to local dof i of patch k.
Definition: gsDofMapper.h:325
Provides a list of labeled parameters/options that can be set and accessed easily.
bool is_free_index(index_t gl) const
Returns true if global dof gl is not eliminated.
Definition: gsDofMapper.h:376
Class representing a (scalar) hierarchical tensor basis of functions .
Definition: gsHTensorBasis.h:74
Holds a set of patch-wise bases and their topology information.
Definition: gsMultiBasis.h:36
short_t maxCwiseDegree() const
Maximum degree with respect to all variables.
Definition: gsMultiBasis.hpp:342
Provides declaration of the MultiPatch class.
Interface for the set of functions defined on a domain (the total number of functions in the set equa...
Definition: gsFuncData.h:23
void freeAll(It begin, It end)
Frees all pointers in the range [begin end)
Definition: gsMemory.h:312
std::vector< gsBasis< T > * > basesCopy(bool numeratorOnly=false) const
Makes a deep copy of all bases and puts them in a vector.
Definition: gsMultiPatch.hpp:178
Container class for a set of geometry patches and their topology, that is, the interface connections ...
Definition: gsMultiPatch.h:33
Struct which represents a certain side of a box.
Definition: gsBoundary.h:84
virtual const gsBasis & source() const
Definition: gsBasis.h:704
uPtr clone()
Clone methode. Produceds a deep copy inside a uPtr.
Struct which represents a certain component (interior, face, egde, corner) of a particular patch...
Definition: gsBoundary.h:565
memory::unique_ptr< gsBasis > uPtr
Unique pointer for gsBasis.
Definition: gsBasis.h:89
const point & upperCorner() const
Accessor for gsHDomain::m_upperIndex.
Definition: gsHDomain.h:249
short_t minDegree(short_t k) const
Minimum degree with respect to variable k.
Definition: gsMultiBasis.hpp:362
Class containing a set of boundary conditions.
Definition: gsBoundaryConditions.h:341
index_t askInt(const std::string &label, const index_t &value=0) const
Reads value for option label from options.
Definition: gsOptionList.cpp:117
void uniformRefine_withTransfer(gsSparseMatrix< T, RowMajor > &transfer, const gsBoundaryConditions< T > &boundaryConditions, const gsOptionList &assemblerOptions, int numKnots=1, int mul=1, index_t unk=0)
Refine every basis uniformly.
Definition: gsMultiBasis.hpp:181
bool repairInterfaceFindElements(const boundaryInterface &bi, std::vector< index_t > &refEltsFirst, std::vector< index_t > &refEltsSecond)
Finds the elements that need to be refined in order to repair an interface.
Definition: gsMultiBasis.hpp:462
~gsMultiBasis()
Destructor.
Definition: gsMultiBasis.hpp:65
patchSide & second()
second, returns the second patchSide of this interface
Definition: gsBoundary.h:782
void addBox(index_t i=1)
Add i new boxes.
Definition: gsBoxTopology.h:198
Struct which represents an interface between two patches.
Definition: gsBoundary.h:649
short_t index() const
Returns the index (as specified in boundary::side) of the box side.
Definition: gsBoundary.h:140
Class which holds a list of parameters/options, and provides easy access to them. ...
Definition: gsOptionList.h:32
index_t freeSize() const
Returns the number of free (not eliminated) dofs.
Definition: gsDofMapper.h:436
void sortByColumn(const index_t j)
Sorts rows of matrix by column j.
Definition: gsMatrix.h:388
index_t patch
The index of the patch.
Definition: gsBoundary.h:234
short_t direction() const
Returns the parametric direction orthogonal to this side.
Definition: gsBoundary.h:113
void getBoxesOnSide(boundary::side s, gsMatrix< Z > &b1, gsMatrix< Z > &b2, gsVector< Z > &level) const
Returns the boxes which make up the hierarchical domain and the respective levels touching side s...
Definition: gsHDomain.hpp:739
A basis represents a family of scalar basis functions defined over a common parameter domain...
Definition: gsBasis.h:78
void uniformCoarsen_withTransfer(gsSparseMatrix< T, RowMajor > &transfer, const gsBoundaryConditions< T > &boundaryConditions, const gsOptionList &assemblerOptions, int numKnots=1, index_t unk=0)
Coarsen every basis uniformly.
Definition: gsMultiBasis.hpp:222
static void combineTransferMatrices(const std::vector< gsSparseMatrix< T, RowMajor > > &localTransferMatrices, const gsDofMapper &coarseMapper, const gsDofMapper &fineMapper, gsSparseMatrix< T, RowMajor > &transferMatrix)
This function takes local transfer matrices (per patch) and combines them using given DofMappers to a...
Definition: gsMultiBasis.hpp:137
void addBasis(gsBasis< T > *g)
Add a basis (ownership of the pointer is also acquired)
Definition: gsMultiBasis.hpp:79
patchSide & first()
first, returns the first patchSide of this interface
Definition: gsBoundary.h:776