G+Smo  24.08.0
Geometry + Simulation Modules
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
gsMultiPatch.hpp
Go to the documentation of this file.
1 
14 #pragma once
15 
16 #include <gsCore/gsBasis.h>
17 #include <gsCore/gsGeometry.h>
18 #include <gsCore/gsDofMapper.h>
20 
22 
23 #include <gsMesh2/gsSurfMesh.h>
24 #include <gsTensor/gsTensorBasis.h>
25 
26 namespace gismo
27 {
28 
29 template<class T>
30 gsMultiPatch<T> gsMultiPatch<T>::coord(const index_t c) const
31 {
32  gsMultiPatch<T> result;
33  for ( const_iterator it = m_patches.begin(); it != m_patches.end(); ++it )
34  {
35  result.addPatch( (*it)->coord(c) );
36  }
37  return result;
38 }
39 
40 template<class T>
42  : BaseA( geo.parDim() )
43 {
44  m_patches.push_back(geo.clone().release());
45  //m_patches[0]->setId(0); // Note: for the single-patch constructor the id remains unchanged
46  addBox();
47  this->addAutoBoundaries();
48 }
49 
50 template<class T>
52  : BaseA( other ), BaseB( other ), m_patches( other.m_patches.size() )
53 {
54  // clone all geometries
55  cloneAll( other.m_patches.begin(), other.m_patches.end(),
56  this->m_patches.begin());
57 }
58 
59 #if EIGEN_HAS_RVALUE_REFERENCES
60 
61 template<class T>
63 {
64  if (this!=&other)
65  {
66  freeAll(m_patches);
67  BaseA::operator=(other);
68  m_patches.resize(other.m_patches.size());
69  // clone all geometries
70  cloneAll( other.m_patches.begin(), other.m_patches.end(),
71  this->m_patches.begin());
72  }
73  return *this;
74 }
75 
76 template<class T>
77 gsMultiPatch<T>& gsMultiPatch<T>::operator=( gsMultiPatch&& other )
78 {
79  freeAll(m_patches);
80  BaseA::operator=(give(other));
81  m_patches = give(other.m_patches);
82  return *this;
83 }
84 
85 
86 #endif
87 
88 template<class T>
89 gsMultiPatch<T>::gsMultiPatch(PatchContainer & patches )
90  : BaseA( patches[0]->parDim(), patches.size() )
91 {
92  m_patches.swap(patches); // patches are consumed
93  setIds();
94  this->addAutoBoundaries();
95 }
96 
97 template<class T>
98 gsMultiPatch<T>::gsMultiPatch( PatchContainer& patches,
99  const std::vector<patchSide>& boundary,
100  const std::vector<boundaryInterface>& interfaces )
101  : BaseA( patches[0]->parDim(), patches.size(), boundary, interfaces )
102 {
103  m_patches.swap(patches); // patches are consumed
104  setIds();
105 }
106 
107 template<class T>
109 {
110  freeAll(m_patches);
111 }
112 
113 template<class T>
115 {
116  size_t id = 0;
117  for ( iterator it = m_patches.begin(); it != m_patches.end(); ++it )
118  {
119  ( *it )->setId( id++ );
120  }
121 }
122 
123 template<class T>
124 std::ostream& gsMultiPatch<T>::print(std::ostream& os) const
125 {
126  if ( !this->empty() ) {
127  os << "gsMultiPatch (" << this->nPatches() << "): ";
128  os << "#Boundaries= " << nBoundary() << ", ";
129  os << "#Interfaces= " << nInterfaces() << ".\n";
130  } else {
131  os << "gsMultiPatch ( empty! ).\n";
132  }
133  return os;
134 }
135 
136 template<class T>
137 std::string gsMultiPatch<T>::detail() const
138 {
139  std::ostringstream os;
140  print(os);
141  if ( nPatches() > 0 )
142  {
143  BaseA::print( os );
144  }
145  return os.str();
146 }
147 
148 template<class T>
150 {
151  GISMO_ASSERT( m_patches.size() > 0 , "Empty multipatch object.");
152  return m_patches[0]->geoDim();
153 }
154 
155 template<class T>
157 {
158  GISMO_ASSERT( m_patches.size() > 0 , "Empty multipatch object.");
159  return m_patches[0]->geoDim() - m_dim;
160 }
161 
162 template<class T>
165 {
166  return m_patches[i]->basis().support();
167 }
168 
169 template<class T>
170 gsBasis<T> &
171 gsMultiPatch<T>::basis( size_t i ) const
172 {
173  GISMO_ASSERT( i < m_patches.size(), "Invalid patch index requested from gsMultiPatch" );
174  return m_patches[i]->basis();
175 }
176 
177 template<class T>
178 std::vector<gsBasis<T> *> gsMultiPatch<T>::basesCopy(bool NoRational) const
179 {
180  std::vector<gsBasis<T> *> bb;
181  bb.reserve(m_patches.size());
182 
183  if (NoRational)
184  {
185  for ( const_iterator it = m_patches.begin();
186  it != m_patches.end(); ++it )
187  {
188  bb.push_back( (*it)->basis().source().clone().release() );
189  }
190  }
191  else
192  {
193  for ( const_iterator it = m_patches.begin();
194  it != m_patches.end(); ++it )
195  {
196  bb.push_back( (*it)->basis().clone().release() );
197  }
198  }
199  return bb ;
200 }
201 
202 template<class T>
203 void gsMultiPatch<T>::permute(const std::vector<short_t> & perm)
204 {
205  gsAsVector<gsGeometry<T>*> a (m_patches);
206  a = gsEigen::PermutationMatrix<-1,-1,short_t>(gsAsConstVector<short_t>(perm)) * a;
207 }
208 
209 template<class T>
211 {
212  if ( m_dim == -1 )
213  {
214  m_dim = g->parDim();
215  } else
216  {
217  GISMO_ASSERT( m_dim == g->parDim(),
218  "Tried to add a patch of different dimension in a multipatch." );
219  }
220  index_t index = m_patches.size();
221  g->setId(index);
222  m_patches.push_back( g.release() ) ;
223  addBox();
224  return index;
225 }
226 
227 template<class T>
229 {
230  return addPatch(g.clone());
231 }
232 
233 template<class T>
235 {
236  const_iterator it
237  = std::find( m_patches.begin(), m_patches.end(), g );
238  GISMO_ASSERT( it != m_patches.end(), "Did not find the patch index." );
239  // note: should return g->patchId();
240  return it - m_patches.begin();
241 }
242 
243 template<class T>
245  gsGeometry<T>* g2, boxSide s2 )
246 {
247  int p1 = findPatchIndex( g1 );
248  int p2 = findPatchIndex( g2 );
249  BaseA::addInterface( p1, s1, p2, s2 );
250 }
251 
252 template<class T>
254 {
255  gsMatrix<T> coordinates;
256  m_patches[pc.patch]->eval_into(m_patches[pc.patch]->parameterCenter(pc),coordinates);
257  return coordinates;
258 }
259 
260 
261 template<class T>
263 {
264  gsMatrix<T> coordinates;
265  m_patches[ps.patch]->eval_into(m_patches[ps.patch]->parameterCenter(ps),coordinates);
266  return coordinates;
267 }
268 
269 template<class T>
270 void gsMultiPatch<T>::uniformRefine(int numKnots, int mul)
271 {
272  for ( typename PatchContainer::const_iterator it = m_patches.begin();
273  it != m_patches.end(); ++it )
274  {
275  ( *it )->uniformRefine(numKnots, mul);
276  }
277 }
278 
279 
280 template<class T>
281 void gsMultiPatch<T>::degreeElevate(short_t const elevationSteps, short_t const dir)
282 {
283  for ( typename PatchContainer::const_iterator it = m_patches.begin();
284  it != m_patches.end(); ++it )
285  {
286  ( *it )->degreeElevate(elevationSteps, dir);
287  }
288 }
289 
290 template<class T>
291 void gsMultiPatch<T>::degreeIncrease(short_t const elevationSteps, short_t const dir)
292 {
293  for ( typename PatchContainer::const_iterator it = m_patches.begin();
294  it != m_patches.end(); ++it )
295  {
296  ( *it )->degreeIncrease(elevationSteps, dir);
297  }
298 }
299 
300 template<class T>
301 void gsMultiPatch<T>::degreeReduce(int elevationSteps)
302 {
303  for ( typename PatchContainer::const_iterator it = m_patches.begin();
304  it != m_patches.end(); ++it )
305  {
306  ( *it )->degreeReduce(elevationSteps, -1);
307  }
308 }
309 
310 template<class T>
312 {
313  for ( typename PatchContainer::const_iterator it = m_patches.begin();
314  it != m_patches.end(); ++it )
315  {
316  ( *it )->uniformCoarsen(numKnots);
317  }
318 }
319 
320 template<class T>
322 {
323  result.setZero(geoDim(),2);
324  if ( m_patches.size() == 0 )
325  return;
326 
327  // to do: this is the bounding box of a gsGeometry, make a member function
328  result.col(0) = patch(0).coefs().colwise().minCoeff().transpose();
329  result.col(1) = patch(0).coefs().colwise().maxCoeff().transpose();
330 
331  for (const_iterator it = begin()+1; it != end(); ++it)
332  {
333  const gsMatrix<T> & cc = (*it)->coefs();
334  result.col(0) = result.col(0).cwiseMin( cc.colwise().minCoeff().transpose() ) ;
335  result.col(1) = result.col(1).cwiseMax( cc.colwise().maxCoeff().transpose() ) ;
336  }
337 }
338 
339 template<class T>
341 {
342  int n;
343  if (dir == -1)
344  n = math::exp2(parDim());
345  else
346  n = 2;
347  std::vector<gsGeometry<T>*> result;
348  result.reserve(nPatches() * n);
349 
350  for (size_t np = 0; np < nPatches(); ++np)
351  {
352  std::vector<gsGeometry<T>*> result_temp = m_patches[np]->uniformSplit(dir);
353  result.insert(result.end(), result_temp.begin(), result_temp.end());
354  }
355  gsMultiPatch<T> mp(result);
356  mp.computeTopology();
357  return mp;
358 }
359 
360 
361 /*
362  This is based on comparing a set of reference points of the patch
363  side and thus it implicitly assumes that the patch faces match
364 */
365 template<class T>
366 bool gsMultiPatch<T>::computeTopology( T tol, bool cornersOnly, bool)
367 {
368  BaseA::clearTopology();
369 
370  const size_t np = m_patches.size();
371  const index_t nCorP = 1 << m_dim; // corners per patch
372  const index_t nCorS = 1 << (m_dim-1); // corners per side
373 
374  gsMatrix<T> supp,
375  // Parametric coordinates of the reference points. These points
376  // are used to decide if two sides match.
377  // Currently these are the corner points and the side-centers
378  coor;
379  if (cornersOnly)
380  coor.resize(m_dim,nCorP);
381  else
382  coor.resize(m_dim,nCorP + 2*m_dim);
383 
384  gsVector<bool> boxPar(m_dim);
385 
386  // each matrix contains the physical coordinates of the reference points
387  std::vector<gsMatrix<T> > pCorners(np);
388 
389  std::vector<patchSide> pSide; // list of all candidate patchSides to compare
390  pSide.reserve(np * 2 * m_dim);
391 
392 //# pragma omp parallel for private(supp, boxPar, !coor)
393  for (size_t p=0; p<np; ++p)
394  {
395  supp = m_patches[p]->parameterRange(); // the parameter domain of patch i
396 
397  // Corners' parametric coordinates
398  for (boxCorner c=boxCorner::getFirst(m_dim); c<boxCorner::getEnd(m_dim); ++c)
399  {
400  boxPar = c.parameters(m_dim);
401  for (index_t i=0; i<m_dim;++i)
402  coor(i,c-1) = boxPar(i) ? supp(i,1) : supp(i,0);
403  }
404 
405  if (!cornersOnly)
406  {
407  // Sides' centers parametric coordinates
408  index_t l = nCorP;
409  for (boxSide c=boxSide::getFirst(m_dim); c<boxSide::getEnd(m_dim); ++c)
410  {
411  const index_t dir = c.direction();
412  const index_t s = static_cast<index_t>(c.parameter());// 0 or 1
413 
414  for (index_t i=0; i<m_dim;++i)
415  coor(i,l) = ( dir==i ? supp(i,s) :
416  (supp(i,1)+supp(i,0))/2.0 );
417  l++;
418  }
419  }
420 
421  // Evaluate the patch on the reference points
422  m_patches[p]->eval_into(coor,pCorners[p]);
423 
424  // Add all the patchSides of this patch to the candidate list
425  for (boxSide bs=boxSide::getFirst(m_dim); bs<boxSide::getEnd(m_dim); ++bs)
426  pSide.push_back(patchSide(p,bs));
427  }
428 
429  gsVector<index_t> dirMap(m_dim);
430  gsVector<bool> matched(nCorS), dirOr(m_dim);
431  std::vector<boxCorner> cId1, cId2;
432  cId1.reserve(nCorS);
433  cId2.reserve(nCorS);
434 
435  std::set<index_t> found;
436  for (size_t sideind=0; sideind<pSide.size(); ++sideind)
437  {
438  const patchSide & side = pSide[sideind];
439  for (size_t other=sideind+1; other<pSide.size(); ++other)
440  {
441  side .getContainedCorners(m_dim,cId1);
442  pSide[other].getContainedCorners(m_dim,cId2);
443  matched.setConstant(false);
444 
445  // Check whether the side center matches
446  if (!cornersOnly)
447  if ( ( pCorners[side.patch ].col(nCorP+side-1 ) -
448  pCorners[pSide[other].patch].col(nCorP+pSide[other]-1)
449  ).norm() >= tol )
450  continue;
451 
452  //t-junction
453  // check for matching vertices else
454  // invert the vertices of first side on the second and vise-versa
455  // if at least one vertex is found (at most 2^(d-1)), mark as interface
456 
457  // Check whether the vertices match and compute direction
458  // map and orientation
459  if ( matchVerticesOnSide( pCorners[side.patch] , cId1, 0,
460  pCorners[pSide[other].patch], cId2,
461  matched, dirMap, dirOr, tol ) )
462  {
463  dirMap(side.direction()) = pSide[other].direction();
464  dirOr (side.direction()) = !( side.parameter() == pSide[other].parameter() );
465  BaseA::addInterface( boundaryInterface(side, pSide[other], dirMap, dirOr));
466  found.insert(sideind);
467  found.insert(other);
468  }
469  }
470  }
471 
472  index_t k = 0;
473  found.insert(found.end(), pSide.size());
474  for (const auto & s : found)
475  {
476  for (;k<s;++k)
477  BaseA::addBoundary( pSide[k] );
478  ++k;
479  }
480 
481  return true;
482 }
483 
484 template <class T>
486 {
487  for ( typename PatchContainer::const_iterator it = m_patches.begin();
488  it != m_patches.end(); ++it )
489  if ( -1 == (*it)->orientation() )
490  (*it)->toggleOrientation();
491 
492  if (this->nInterfaces() || this->nBoundary() )
493  this->computeTopology();
494 }
495 
496 template <class T>
498  const gsMatrix<T> &cc1, const std::vector<boxCorner> &ci1, index_t start,
499  const gsMatrix<T> &cc2, const std::vector<boxCorner> &ci2,
500  const gsVector<bool> &matched, gsVector<index_t> &dirMap,
501  gsVector<bool> &dirO, T tol, index_t reference)
502 {
503  const bool computeOrientation = !(start&(start-1)) && (start != 0); // true if start is a power of 2
504  const bool setReference = start==0; // if we search for the first point then we set the reference
505 
506  const short_t dim = static_cast<short_t>(cc1.rows());
507 
508  index_t o_dir = 0, d_dir = 0;
509 
510  gsVector<bool> refPar, newPar, newMatched;
511 
512  if (computeOrientation)
513  {
514  // o_dir is the only direction in which the parameters differ between start and 0
515  const gsVector<bool> parStart = ci1[start].parameters(dim);
516  const gsVector<bool> parRef = ci1[0].parameters(dim);
517  for (; o_dir<dim && parStart(o_dir)==parRef(o_dir) ;) ++o_dir;
518  }
519 
520  if (!setReference)
521  refPar = ci2[reference].parameters(dim);
522 
523  for (size_t j=0;j<ci2.size();++j)
524  {
525  if( !matched(j) && (cc1.col(ci1[start]-1)-cc2.col(ci2[j]-1)).norm() < tol )
526  {
527  index_t newRef = (setReference) ? j : reference;
528  if (computeOrientation)
529  {
530  index_t count=0;
531  d_dir = 0;
532  newPar = ci2[j].parameters(dim);
533  for (index_t i=0; i< newPar.rows();++i)
534  {
535  if ( newPar(i)!=refPar(i) )
536  {
537  d_dir=i;
538  ++count;
539  }
540  }
541  if (count != 1)
542  {
543  // the match is wrong: we are mapping an edge to a diagonal
544  continue;
545  }
546  dirMap(o_dir) = d_dir;
547  dirO (o_dir) = (static_cast<index_t>(j) > reference);
548  }
549  if ( start + 1 == static_cast<index_t>( ci1.size() ) )
550  {
551  // we matched the last vertex, we are done
552  return true;
553  }
554  newMatched = matched;
555  newMatched(j) = true;
556  if (matchVerticesOnSide( cc1,ci1,start+1,cc2,ci2,newMatched,dirMap,dirO, tol,newRef))
557  return true;
558  }
559  }
560 
561  return false;
562 }
563 
564 
565 template<class T>
567 {
568  gsDofMapper mapper = getMapper(tol);
569 
570  gsMatrix<T> meanVal;
571  std::vector<std::pair<index_t,index_t> > dof;
572  const index_t start = mapper.freeSize() - mapper.coupledSize();
573  const index_t end = mapper.freeSize();
574 
575  for (index_t i = start; i!= end; ++i) // For all coupled DoFs
576  {
577  // Get the preimages of this global dof (as pairs (patch,index) )
578  mapper.preImage(i, dof);
579 
580  // Compute the mean value
581  meanVal = m_patches[dof.front().first]->coef(dof.front().second);
582  for (size_t k = 1; k!=dof.size(); ++k)
583  meanVal += m_patches[dof[k].first]->coef(dof[k].second);
584  meanVal.array() /= dof.size();
585 
586  // Set involved control points equal to their average value
587  for (size_t k = 0; k!=dof.size(); ++k)
588  m_patches[dof[k].first]->coef(dof[k].second) = meanVal;
589  }
590 }
591 
592 template<class T>
594 {
595  const T tol2 = tol*tol;
596  gsMatrix<index_t> bdr1, bdr2; // indices of the boundary control points
597 
598  // Create a map which assigns to all meeting patch-local indices a
599  // unique global id
600  const index_t sz = m_patches.size();
601  gsVector<index_t> patchSizes(sz);
602  for (index_t i=0; i!= sz; ++i)
603  patchSizes[i] = m_patches[i]->coefsSize();
604 
605  gsDofMapper mapper(patchSizes);
606 
607  for ( const_iiterator it = iBegin(); it != iEnd(); ++it ) // for all interfaces
608  {
609  const gsGeometry<T> & p1 = *m_patches[it->first() .patch];
610  const gsGeometry<T> & p2 = *m_patches[it->second().patch];
611 
612  // Grab boundary control points in matching configuration
613  p1.basis().matchWith(*it, p2.basis(), bdr1, bdr2);
614 
615  bool warn = true;
616  //mapper.matchDofs(it->first().patch, bdr1, it->second().patch, bdr2);
617  for (index_t i = 0; i!= bdr1.size(); ++i )
618  {
619  if ( ( p1.coef(bdr1(i)) - p2.coef(bdr2(i)) ).squaredNorm() > tol2 )
620  {
621  if (warn)
622  {
623  gsWarn<<"Big gap detected between patches "<< it->first().patch
624  <<" and "<<it->second().patch <<"\n";
625  warn = false;
626  }
627  }
628  else
629  // Match the dofs on the interface
630  mapper.matchDof(it->first().patch, bdr1(i,0), it->second().patch, bdr2(i,0) );
631  }
632  }
633 
634  // Finalize the mapper. At this point all patch-local dofs are
635  // mapped to unique global indices
636  mapper.finalize();
637  return mapper;
638 }
639 
640 template<class T>
642 {
643  GISMO_ASSERT(2==parDim(), "Works for surfaces only.");
644  gsDofMapper mapper = getMapper((T)1e-7);
645  gsSurfMesh mesh;
646  auto pid = mesh.add_vertex_property<index_t>("v:patch");
647  auto anchor = mesh.add_vertex_property<index_t>("v:anchor");
649  gsSurfMesh::Point pt(0,0,0);
650  const index_t gd = geoDim();
651  std::vector<std::pair<index_t,index_t> > pi = mapper.anyPreImages();
652  //std::pair<index_t,index_t> pi;
653 
654  for (index_t j = 0; j!= mapper.size(); ++j)
655  {
656  //pi = mapper.anyPreImage(j);
657  gsGeometry<> & pp = patch(pi[j].first);
658  pt.topRows(gd) = pp.eval( pp.basis().anchor(pi[j].second) );
659  v = mesh.add_vertex( pt );
660  pid[v] = pi[j].first;
661  anchor[v] = pi[j].second;
662  }
663 
664  size_t np = nPatches();
665  gsMatrix<> supp, coor;
666  gsVector<bool> boxPar(m_dim);
667  gsVector<index_t,2> cur, csize, strides;
668  GISMO_ENSURE( dynamic_cast<gsTensorBasis<2>*>(&patch(0).basis()), "Not a tensor basis");
669  static_cast<gsTensorBasis<2>&>(patch(0).basis()).stride_cwise(strides);
670  static_cast<gsTensorBasis<2>&>(patch(0).basis()).size_cwise (csize);
671  csize.array() -= 2;
672  gsSurfMesh::Vertex v1, v2, v3, v4;
673  for (size_t p=0; p<np; ++p)
674  {
675  // todo: basis->connectivityAtAnchors ++ basis->controlPolytope
676  gsTensorBasis<2>& pp = static_cast<gsTensorBasis<2>&>(patch(p).basis());
677  cur.setZero(2);
678  do
679  {
680  index_t ci = pp.index(cur);
681  v1 = gsSurfMesh::Vertex( mapper.index(ci, p) );
682  ci += strides[0];
683  v2 = gsSurfMesh::Vertex( mapper.index(ci, p) );
684  ci += strides[1];
685  v3 = gsSurfMesh::Vertex( mapper.index(ci, p) );
686  ci -= strides[0];
687  v4 = gsSurfMesh::Vertex( mapper.index(ci, p) );
688  mesh.add_quad(v1,v2,v3,v4);
689  } while (nextCubePoint(cur, csize));
690 
691  }
692 
693  return mesh;
694 }
695 
696 template<class T> // to do: move to boundaryInterface
698 {
699  if (scaling==0)
700  scaling=1;
701  gsMatrix<T> box1=m_patches[bi.first().patch ]->support();
702  gsMatrix<T> box2=m_patches[bi.second().patch]->support();
703  const index_t oDir1 = bi.first() .direction();
704  const index_t oDir2 = bi.second().direction();
705  const T len1=box1(oDir1,1)-box1(oDir1,0);
706 
707  if (bi.second().parameter())
708  {
709  box2(oDir2,0)=box2(oDir2,1);
710  box2(oDir2,1)+=scaling*len1;
711  }
712  else
713  {
714  box2(oDir2,1)=box2(oDir2,0);
715  box2(oDir2,0)-=scaling*len1;
716  }
717 
718  return gsAffineFunction<T>(bi.dirMap(bi.first()),bi.dirOrientation(bi.first()) ,box1,box2);
719 }
720 
721 template<class T>
723 {
724  GISMO_UNUSED(nsamples);
725  gsMultiPatch<T> result;
726  for ( typename PatchContainer::const_iterator it = m_patches.begin();
727  it != m_patches.end(); ++it )
728  {
729  result.addPatch( (*it)->approximateLinearly() );
730  }
731 
732  result.gsBoxTopology::operator=(*this);//copy the original topology
733  return result;
734 }
735 
736 template<class T>
738 {
739  gsMultiBasis<T> multiBasis(*this);
740  bool changed = false;
741 
742  std::vector<index_t> refEltsFirst;
743  std::vector<index_t> refEltsSecond;
744 
745  // Find the areas/elements that do not match...
746  switch( this->dim() )
747  {
748  case 2:
749  changed = multiBasis.template repairInterfaceFindElements<2>( bi, refEltsFirst, refEltsSecond );
750  break;
751  case 3:
752  changed = multiBasis.template repairInterfaceFindElements<3>( bi, refEltsFirst, refEltsSecond );
753  break;
754  default:
755  GISMO_ASSERT(false,"wrong dimension");
756  }
757 
758  // ...and if there are any found, refine the bases accordingly
759  if( changed )
760  {
761  if( refEltsFirst.size() > 0 )
762  {
763  int pi( bi.first().patch );
764  patch(pi).basis().refineElements_withCoefs( patch(pi).coefs(), refEltsFirst );
765  }
766  if( refEltsSecond.size() > 0 )
767  {
768  int pi( bi.second().patch );
769  patch(pi).basis().refineElements_withCoefs( patch(pi).coefs(), refEltsSecond );
770  }
771  }
772 
773  return changed;
774 }
775 
776 
777 
778 template<class T>
780  gsVector<index_t> & pids,
781  gsMatrix<T> & preim, const T accuracy) const
782 {
783  pids.resize(points.cols());
784  pids.setConstant(-1); // -1 implies not in the domain
785  preim.resize(parDim(), points.cols());//uninitialized by default
786  gsMatrix<T> pt, pr, tmp;
787 
788  for (index_t i = 0; i!=pids.size(); ++i)
789  {
790  pt = points.col(i);
791 
792  for (size_t k = 0; k!= m_patches.size(); ++k)
793  {
794  pr = m_patches[k]->parameterRange();
795  m_patches[k]->invertPoints(pt, tmp, accuracy);
796  if ( (tmp.array() >= pr.col(0).array()).all()
797  && (tmp.array() <= pr.col(1).array()).all() )
798  {
799  pids[i] = k;
800  preim.col(i) = tmp;
801  break;
802  }
803  }
804  }
805 }
806 
807 template<class T>
809  gsVector<index_t> & pid2, gsMatrix<T> & preim) const
810 {
811  // Assumes points are found on pid1 and possibly on one more patch
812  pid2.resize(points.cols());
813  pid2.setConstant(-1); // -1 implies not in the domain
814  preim.resize(parDim(), points.cols());//uninitialized by default
815  gsMatrix<T> pt, pr, tmp;
816 
817  for (index_t i = 0; i!=pid2.size(); ++i)
818  {
819  pt = points.col(i);
820 
821  for (size_t k = 0; k!= m_patches.size(); ++k)
822  {
823  if (pid1==(index_t)k) continue; // skip pid1
824 
825  pr = m_patches[k]->parameterRange();
826  m_patches[k]->invertPoints(pt, tmp);
827  if ( (tmp.array() >= pr.col(0).array()).all()
828  && (tmp.array() <= pr.col(1).array()).all() )
829  {
830  pid2[i] = k;
831  preim.col(i) = tmp;
832  break;
833  }
834  }
835  }
836 }
837 
838 namespace
839 {
840 struct __closestPointHelper
841 {
842  __closestPointHelper() : dist(math::limits::max()), pid(-1) { }
843  real_t dist;
844  index_t pid;
845  gsVector<> preim;
846 };
847 }
848 
849 template<class T> std::pair<index_t,gsVector<T> >
850 gsMultiPatch<T>::closestPointTo(const gsVector<T> & pt,
851  const T accuracy) const
852 {
853  std::pair<index_t,gsVector<T> > result;
854  this->closestDistance(pt,result,accuracy);
855  return result;
856 }
857 
858 
859 template<class T>
860 T gsMultiPatch<T>::closestDistance(const gsVector<T> & pt,
861  std::pair<index_t,gsVector<T> > & result,
862  const T accuracy) const
863 {
864  GISMO_ASSERT( pt.rows() == targetDim(), "Invalid input point." <<
865  pt.rows() <<"!="<< targetDim() );
866 
867  gsVector<T> tmp;
868 
869 #ifndef _MSC_VER
870 # pragma omp declare reduction(minimum : struct __closestPointHelper : omp_out = (omp_in.dist < omp_out.dist ? omp_in : omp_out) )
871  struct __closestPointHelper cph;
872 # pragma omp parallel for default(shared) private(tmp) reduction(minimum:cph) //OpenMP 4.0, will not work on VS2019
873 #else
874  struct __closestPointHelper cph;
875 #endif
876  for (size_t k = 0; k < m_patches.size(); ++k)
877  {
878  // possible improvement: approximate dist: eval patch on a
879  // grid. find min distance between grid and pt
880 
881  const T val = this->patch(k).closestPointTo(pt, tmp, accuracy);
882  if (cph.dist>val)
883  {
884  cph.dist = val; //need to be all in one struct for OMP
885  cph.pid = k;
886  cph.preim = tmp;
887  }
888  }
889  //gsInfo <<"--Pid="<<cph.pid<<", Dist("<<pt.transpose()<<"): "<< cph.dist <<"\n";
890  result = std::make_pair(cph.pid, give(cph.preim));
891  return cph.dist;
892 }
893 
894 template<class T>
896  const index_t nsamples,
897  const T accuracy,
898  const bool directed)
899 {
900  GISMO_ASSERT(this->nPatches()==other.nPatches(),"Number of patches should be the same, but this->nPatches()!=other.nPatches() -> "<<this->nPatches()<<"!="<<other.nPatches());
901  std::vector<T> result(this->nPatches());
902 #pragma omp parallel
903 {
904 # ifdef _OPENMP
905  const int tid = omp_get_thread_num();
906  const int nt = omp_get_num_threads();
907 # endif
908 
909 # ifdef _OPENMP
910  for ( size_t p=tid; p<this->nPatches(); p+=nt )
911 # else
912  for ( size_t p=0; p<this->nPatches(); p++ )
913 # endif
914  {
915  result.at(p) = this->patch(p).HausdorffDistance(other.patch(p),nsamples,accuracy,directed);
916  }
917 }//omp parallel
918  return result;
919 }
920 
921 template<class T>
923  const index_t nsamples,
924  const T accuracy,
925  bool directed)
926 {
927  std::vector<T> distances = HausdorffDistance(other,nsamples,accuracy,directed);
928  return std::accumulate(distances.begin(), distances.end(), (T)( 0 ) ) / distances.size();
929 }
930 
931 template<class T>
932 void gsMultiPatch<T>::constructInterfaceRep()
933 {
934  m_ifaces.clear();
935  for ( iiterator it = iBegin(); it != iEnd(); ++it ) // for all interfaces
936  {
937  const gsGeometry<T> & p1 = *m_patches[it->first() .patch];
938  const gsGeometry<T> & p2 = *m_patches[it->second().patch];
939  m_ifaces[*it] = p1.iface(*it,p2);
940  }//end for
941 }
942 
943 template<class T>
945 {
946  m_bdr.clear();
947  for ( biterator it = bBegin(); it != bEnd(); ++it ) // for all boundaries
948  {
949  const gsGeometry<T> & p1 = *m_patches[it->patch];
950  m_bdr[*it] = p1.boundary(*it);
951  }//end for
952 }
953 
954 template<class T>
955 void gsMultiPatch<T>::constructInterfaceRep(const std::string l)
956 {
957  m_ifaces.clear();
958  ifContainer ifaces = this->interfaces(l);
959  for ( iiterator it = ifaces.begin(); it != ifaces.end(); ++it ) // for all interfaces
960  {
961  const gsGeometry<T> & p1 = *m_patches[it->first() .patch];
962  const gsGeometry<T> & p2 = *m_patches[it->second().patch];
963  m_ifaces[*it] = p1.iface(*it,p2);
964  }//end for
965 }
966 
967 template<class T>
968 void gsMultiPatch<T>::constructBoundaryRep(const std::string l)
969 {
970  m_bdr.clear();
971  bContainer bdrs = this->boundaries(l);
972  for ( biterator it = bdrs.begin(); it != bdrs.end(); ++it ) // for all boundaries
973  {
974  const gsGeometry<T> & p1 = *m_patches[it->patch];
975  m_bdr[*it] = p1.boundary(*it);
976  }//end for
977 }
978 
979 template<class T>
981 {
982  for ( biterator it = bBegin(); it != bEnd(); ++it ) // for all boundaries
983  {
984  const gsGeometry<T> & p1 = *m_patches[it->patch];
985  m_sides[*it] = p1.boundary(*it);
986  }//end for
987 
988  for ( iiterator it = iBegin(); it != iEnd(); ++it ) // for all interfaces
989  {
990  const gsGeometry<T> & p1 = *m_patches[it->first() .patch];
991  const gsGeometry<T> & p2 = *m_patches[it->second().patch];
992  m_sides[it->first()] = p1.boundary(it->first());
993  m_sides[it->second()] = p2.boundary(it->second());
994  }//end for
995 }
996 
997 
998 } // namespace gismo
index_t addPatch(typename gsGeometry< T >::uPtr g)
Add a patch from a gsGeometry&lt;T&gt;::uPtr.
Definition: gsMultiPatch.hpp:210
Abstract base class representing a geometry map.
Definition: gsGeometry.h:92
gsMatrix< T >::RowXpr coef(index_t i)
Returns the i-th coefficient of the geometry as a row expression.
Definition: gsGeometry.h:346
gsMultiPatch & operator=(gsMultiPatch other)
Assignment operator (uses copy-and-swap idiom)
Definition: gsMultiPatch.h:63
void clear()
Clear (delete) all patches.
Definition: gsMultiPatch.h:319
Struct which represents a certain side of a patch.
Definition: gsBoundary.h:231
Half edge mesh structure.
bool nextCubePoint(Vec &cur, const Vec &end)
Iterate in lexigographic order through the points of the integer lattice contained in the cube [0...
Definition: gsCombinatorics.h:327
short_t geoDim() const
Dimension of the geometry (must match for all patches).
Definition: gsMultiPatch.hpp:149
#define short_t
Definition: gsConfig.h:35
static boxSide getFirst(short_t)
helper for iterating on sides of an n-dimensional box
Definition: gsBoundary.h:160
void addAutoBoundaries()
Make all patch sides which are not yet declared as interface or boundary to a boundary.
Definition: gsBoxTopology.cpp:75
gsBasis< T > & basis(const size_t i) const
Return the basis of the i-th patch.
Definition: gsMultiPatch.hpp:171
bool parameter() const
Returns the parameter value (false=0=start, true=1=end) that corresponds to this side.
Definition: gsBoundary.h:128
bool repairInterface(const boundaryInterface &bi)
Checks if the interface is fully matching, and if not, repairs it.
Definition: gsMultiPatch.hpp:737
Struct that defines the boundary sides and corners and types of a geometric object.
Definition: gsBoundary.h:55
void getContainedCorners(short_t dim, std::vector< patchCorner > &corners) const
returns the vector of the corners contained in the side
Definition: gsBoundary.cpp:35
gsSurfMesh toMesh() const
Creates a surface mesh out of this multipatch.
Definition: gsMultiPatch.hpp:641
Implements an affine function.
Provides declaration of Basis abstract interface.
S give(S &x)
Definition: gsMemory.h:266
gsMatrix< T > eval(const gsMatrix< T > &u) const
Evaluate the function,.
Definition: gsFunctionSet.hpp:120
Provides declaration of Geometry abstract interface.
#define index_t
Definition: gsConfig.h:32
#define GISMO_ENSURE(cond, message)
Definition: gsDebug.h:102
void fixOrientation()
Provides positive orientation for all patches.
Definition: gsMultiPatch.hpp:485
void uniformRefine(int numKnots=1, int mul=1)
Refine uniformly all patches by inserting numKnots in each knot-span with multipliplicity mul...
Definition: gsMultiPatch.hpp:270
void locatePoints(const gsMatrix< T > &points, gsVector< index_t > &pids, gsMatrix< T > &preim, const T accuracy=1e-6) const
For each point in points, locates the parametric coordinates of the point.
Definition: gsMultiPatch.hpp:779
Provides combinatorial unitilies.
Struct which represents a certain corner of a patch.
Definition: gsBoundary.h:392
#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
void boundingBox(gsMatrix< T > &result) const
Returns a bounding box for the multipatch domain. The output result is a matrix with two columns...
Definition: gsMultiPatch.hpp:321
index_t coupledSize() const
Returns the number of coupled (not eliminated) dofs.
Definition: gsDofMapper.cpp:601
void uniformCoarsen(int numKnots=1)
Coarsen uniformly all patches by removing numKnots in each knot-span.
Definition: gsMultiPatch.hpp:311
gsMatrix< T > parameterRange(int i=0) const
Returns the range of parameter.
Definition: gsMultiPatch.hpp:164
A halfedge data structure for polygonal meshes.
Definition: gsSurfMesh.h:45
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 size() const
size
Definition: gsMultiPatch.h:129
Provides the gsDofMapper class for re-indexing DoFs.
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
gsMultiPatch< T > approximateLinearly(index_t nsamples) const
Computes linear approximation of the patches using nsamples per direction.
Definition: gsMultiPatch.hpp:722
Creates a mapped object or data pointer to a vector without copying data.
Definition: gsLinearAlgebra.h:129
Representation of an affine function.
Definition: gsAffineFunction.h:29
void matchDof(index_t u, index_t i, index_t v, index_t j, index_t comp=0)
Couples dof i of patch u with dof j of patch v such that they refer to the same global dof at compone...
Definition: gsDofMapper.cpp:99
static boxSide getEnd(short_t dim)
helper for iterating on sides of an n-dimensional box
Definition: gsBoundary.h:176
#define gsWarn
Definition: gsDebug.h:50
Holds a set of patch-wise bases and their topology information.
Definition: gsMultiBasis.h:36
gsGeometry< T > & patch(size_t i) const
Return the i-th patch.
Definition: gsMultiPatch.h:226
size_t nPatches() const
Number of patches.
Definition: gsMultiPatch.h:208
void permute(const std::vector< short_t > &perm)
Permutes the patches according to perm.
Definition: gsMultiPatch.hpp:203
GISMO_DEPRECATED index_t direction(index_t s)
Returns the parametric direction that corresponds to side s.
Definition: gsBoundary.h:1048
std::ostream & print(std::ostream &os) const
Prints the object as a string.
Definition: gsMultiPatch.hpp:124
size_t findPatchIndex(gsGeometry< T > *g) const
Search for the given geometry and return its patch index.
Definition: gsMultiPatch.hpp:234
std::vector< gsGeometry * > boundary() const
Get boundary of this geometry as a vector of new gsGeometry instances.
Definition: gsGeometry.hpp:423
void degreeReduce(int elevationSteps=1)
Reduce the degree of all patches by elevationSteps.
Definition: gsMultiPatch.hpp:301
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
void setId(const size_t i)
Sets the patch index for this patch.
Definition: gsGeometry.h:607
Abstract base class for tensor product bases.
Definition: gsTensorBasis.h:33
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
std::vector< std::pair< index_t, index_t > > anyPreImages(index_t comp=0) const
For all global index, this function assigns a pair (patch,dof) that maps to that global index...
Definition: gsDofMapper.cpp:525
short_t parDim() const
Dimension d of the parameter domain (same as domainDim()).
Definition: gsGeometry.hpp:190
unsigned index(unsigned i, unsigned j, unsigned k=0) const
Definition: gsTensorBasis.h:882
void constructBoundaryRep()
Construct the boundary representation.
Definition: gsMultiPatch.hpp:944
uPtr clone()
Clone methode. Produceds a deep copy inside a uPtr.
memory::unique_ptr< gsGeometry > uPtr
Unique pointer for gsGeometry.
Definition: gsGeometry.h:100
virtual gsGeometry::uPtr iface(const boundaryInterface &bi, const gsGeometry &other) const
Computes and returns the interface with other as a new geometry.
Definition: gsGeometry.hpp:232
gsAffineFunction< T > getMapForInterface(const boundaryInterface &bi, T scaling=0) const
construct the affine map that places bi.first() next to bi.second() and identifies the two matching s...
Definition: gsMultiPatch.hpp:697
void closeGaps(T tol=1e-4)
Attempt to close gaps between the interfaces. Assumes that the topology is computed, ie. computeTopology() has been called.
Definition: gsMultiPatch.hpp:566
gsMatrix< T > pointOn(const patchCorner &pc)
Get coordinates of the patchCorner pc in the physical domain.
Definition: gsMultiPatch.hpp:253
std::string detail() const
Prints the object as a string with extended details.
Definition: gsMultiPatch.hpp:137
Vertex_property< T > add_vertex_property(const std::string &name, T t=T())
Definition: gsSurfMesh.h:1342
void preImage(index_t gl, std::vector< std::pair< index_t, index_t > > &result) const
For gl being a global index, this function returns a vector of pairs (patch,dof) that contains all th...
Definition: gsDofMapper.cpp:480
std::vector< T > HausdorffDistance(const gsMultiPatch< T > &other, const index_t nsamples=1000, const T accuracy=1e-6, const bool directed=false)
Construct the interface representation.
Definition: gsMultiPatch.hpp:895
#define GISMO_UNUSED(x)
Definition: gsDebug.h:112
index_t size() const
Returns the total number of dofs (free and eliminated).
Definition: gsDofMapper.h:421
Creates a mapped object or data pointer to a const vector without copying data.
Definition: gsLinearAlgebra.h:130
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
static boxCorner getFirst(short_t)
helper for iterating on corners of an n-dimensional box
Definition: gsBoundary.h:356
Struct which represents an interface between two patches.
Definition: gsBoundary.h:649
~gsMultiPatch()
Destructor.
Definition: gsMultiPatch.hpp:108
Provides declaration of TensorBasis class.
index_t freeSize() const
Returns the number of free (not eliminated) dofs.
Definition: gsDofMapper.h:436
gsMultiPatch()
Default empty constructor.
Definition: gsMultiPatch.h:55
Vertex add_vertex(const Point &p)
add a new vertex with position p
Definition: gsSurfMesh.cpp:272
virtual const gsBasis< T > & basis() const =0
Returns a const reference to the basis of the geometry.
static boxCorner getEnd(short_t dim)
helper for iterating on corners of an n-dimensional box
Definition: gsBoundary.h:372
Defines a topological arrangement of a collection of &quot;boxes&quot; (e.g., parameter domains that map to phy...
Definition: gsBoxTopology.h:38
short_t coDim() const
Co-dimension of the geometry (must match for all patches).
Definition: gsMultiPatch.hpp:156
index_t patch
The index of the patch.
Definition: gsBoundary.h:234
GISMO_DEPRECATED void addInterface(gsGeometry< T > *g1, boxSide s1, gsGeometry< T > *g2, boxSide s2)
Add an interface joint between side s1 of geometry g1 side s2 of geometry g2.
Definition: gsMultiPatch.hpp:244
void degreeIncrease(short_t const elevationSteps=1, short_t const dir=-1)
Increase the degree of all patches by elevationSteps, preserves multiplicity.
Definition: gsMultiPatch.hpp:291
Struct which represents a certain corner of a hyper-cube.
Definition: gsBoundary.h:291
short_t direction() const
Returns the parametric direction orthogonal to this side.
Definition: gsBoundary.h:113
Face add_quad(Vertex v1, Vertex v2, Vertex v3, Vertex v4)
Definition: gsSurfMesh.cpp:365
gsDofMapper getMapper(T tol) const
Used to get a mapper with unique vertices.
Definition: gsMultiPatch.hpp:593
Definition: gsSurfMesh.h:103
A basis represents a family of scalar basis functions defined over a common parameter domain...
Definition: gsBasis.h:78
bool computeTopology(T tol=1e-4, bool cornersOnly=false, bool tjunctions=false)
Attempt to compute interfaces and boundaries automatically.
Definition: gsMultiPatch.hpp:366
void degreeElevate(short_t const elevationSteps=1, short_t const dir=-1)
Elevate the degree of all patches by elevationSteps, preserves smoothness.
Definition: gsMultiPatch.hpp:281
patchSide & first()
first, returns the first patchSide of this interface
Definition: gsBoundary.h:776
gsMultiPatch< T > uniformSplit(index_t dir=-1) const
Splits each patch uniformly in each direction (if dir = -1) into two new patches, giving a total numb...
Definition: gsMultiPatch.hpp:340