G+Smo  25.01.0
Geometry + Simulation Modules
 
Loading...
Searching...
No Matches
gsMPBESBSplineBasis.hpp
Go to the documentation of this file.
1
18
19namespace gismo
20{
21
22template<short_t d, class T>
23gsMPBESBSplineBasis<d,T>::gsMPBESBSplineBasis (const BasisContainer & bases, gsBoxTopology const & topol,
24 index_t increaseSmoothnessLevel, index_t minEVDistance)
25{
26 m_topol = topol;
27 for(typename BasisContainer::const_iterator it = bases.begin();it!=bases.end();++it)
28 m_bases.push_back( (gsBasis<T>*)(*it)->clone().release());
29 if(increaseSmoothnessLevel==-1)
30 m_incrSmoothnessDegree=this->maxDegree()-1;
31 else
32 m_incrSmoothnessDegree=increaseSmoothnessLevel;
33 m_minDist=static_cast<unsigned>(std::max<index_t>(std::min<index_t>(m_incrSmoothnessDegree+1,maxDegree()),minEVDistance));
34 _initVertices();
35 _setDistanceOfAllVertices();
36
37 repairPatches();
38
39 if(!_checkTopologyWithBases())
40 GISMO_ERROR("topology and bases not suitable for composite basis");
41 _setMapping();
42}
43
44template<short_t d, class T>
45gsMPBESBSplineBasis<d,T>::gsMPBESBSplineBasis (const BasisContainer & bases, gsBoxTopology const & topol,std::vector<gsMatrix<T> * > coefs,
46 index_t increaseSmoothnessLevel, index_t minEVDistance)
47{
48 m_topol = topol;
49 for(typename BasisContainer::const_iterator it = bases.begin();it!=bases.end();++it)
50 m_bases.push_back( (gsBasis<T>*)(*it)->clone().release());
51 if(increaseSmoothnessLevel==-1)
52 m_incrSmoothnessDegree=this->maxDegree()-1;
53 else
54 m_incrSmoothnessDegree=increaseSmoothnessLevel;
55 m_minDist=static_cast<unsigned>(std::max<index_t>(std::min<index_t>(m_incrSmoothnessDegree+1,maxDegree()),minEVDistance));
56 _initVertices();
57 _setDistanceOfAllVertices();
58
59 repairPatches(coefs);
60
61 if(!_checkTopologyWithBases())
62 GISMO_ERROR("topology and bases not suitable for composite basis");
63 _setMapping();
64}
65
66template<short_t d, class T>
67gsMPBESBSplineBasis<d,T>::gsMPBESBSplineBasis( gsMultiPatch<T> const & mp, index_t increaseSmoothnessLevel,
68 index_t minEVDistance)
69{
70 //topol.computeAllVertices();
71 for (size_t i = 0; i < mp.nPatches(); i++)
72 {
73 GISMO_ASSERT(dynamic_cast<gsTensorBSplineBasis<2> * >(& mp.basis(i))!=NULL,"Bases is not of type gsTensorBSplineBasis<2>");
74 }
75 m_topol = mp;
76 m_bases = mp.basesCopy();
77 if(increaseSmoothnessLevel==-1)
78 m_incrSmoothnessDegree=this->maxDegree()-1;
79 else
80 m_incrSmoothnessDegree=increaseSmoothnessLevel;
81 m_minDist=static_cast<unsigned>(std::max<index_t>(std::min<index_t>(m_incrSmoothnessDegree+1,maxDegree()),minEVDistance));
82 _initVertices();
83 _setDistanceOfAllVertices();
84
85 repairPatches();
86
87 if(!_checkTopologyWithBases())
88 GISMO_ERROR("topology and bases not suitable for composite basis");
89 _setMapping();
90}
91
92template<short_t d, class T>
93gsMPBESBSplineBasis<d,T>::gsMPBESBSplineBasis( const gsMPBESBSplineBasis& other )
94{
95 m_topol = other.m_topol;
96 m_vertices = other.m_vertices;
97 m_distances=other.m_distances;
98 m_minDist=other.m_minDist;
99 // clone all geometries
100 for(ConstBasisIter it = other.m_bases.begin();it!=other.m_bases.end();++it)
101 {
102 m_bases.push_back( (gsBasis<T>*)(*it)->clone().release() );
103 }
104 m_incrSmoothnessDegree=other.m_incrSmoothnessDegree;
105 if(!this->_checkTopologyWithBases())
106 GISMO_ERROR("topology and bases not suitable for composite basis");
107 //_setMapping();
108 m_mapper=new gsWeightMapper<T>(*other.m_mapper);
109 //m_mapFactory= need a method to clone mapFactory
110}
111
112template<short_t d, class T>
113gsMPBESBSplineBasis<d,T> &gsMPBESBSplineBasis<d,T>::operator=( const gsMPBESBSplineBasis& other )
114{
115 m_topol=other.m_topol;
116 m_vertices=other.m_vertices;
117 m_distances=other.m_distances;
118 m_minDist=other.m_minDist;
119 freeAll(m_bases);
120 m_bases.clear();
121 for(ConstBasisIter it = other.m_bases.begin();it!=other.m_bases.end();++it)
122 m_bases.push_back((BasisType*)(*it)->clone().release());
123 if(m_mapper)
124 delete m_mapper;
125 m_mapper=NULL;//other.m_mapper->clone();
126 m_incrSmoothnessDegree=other.m_incrSmoothnessDegree;
127 //m_mapFactory= need a method to clone mapFactory
128 return *this;
129}
130
131template<short_t d, class T>
133{
134 // * Initializer mapper
135 gsMPBESMapB2D<d,T> maker(m_incrSmoothnessDegree,& m_topol,this);
136 if(m_mapper)
137 delete m_mapper;
138 m_mapper = maker.makeMapper();
139}
140
141template<short_t d, class T>
143{
144 if(ps.side()==1||ps.side()==2)
145 return basis(ps.patch).size(1);
146 else
147 return basis(ps.patch).size(0);
148}
149
150template<short_t d, class T>
151bool gsMPBESBSplineBasis<d,T>::isLocallyConnected(indexType i,indexType j) const
152{
153 unsigned patch_i = _getPatch(i), patch_j = _getPatch(j);
154 if( patch_i != patch_j )
155 return false;
156 indexType loc_i = _getPatchIndex(i), loc_j = _getPatchIndex(j);
157 gsVector<index_t,d> vec_i = basis(patch_i).tensorIndex(loc_i);
158 gsVector<index_t,d> vec_j = basis(patch_j).tensorIndex(loc_j);
159 unsigned distance = 0;
160 for( unsigned i2=0; i2< d;++i2 )
161 {
162 unsigned dist = std::abs(static_cast<index_t>(vec_i[i2])-static_cast<index_t>(vec_j[i2]));
163 if( dist != 0 )
164 distance+=dist;
165 }
166 return distance==1;
167}
168
169template<short_t d, class T>
171{
172 result.resize(u.cols());
173 for(index_t i = 0; i<u.cols();++i)
174 {
175 result(i)=static_cast<unsigned>((degree(patch,0)+1)*(degree(patch,1)+1));
176 }
177}
178
179template<short_t d, class T>
180void gsMPBESBSplineBasis<d,T>::refine(const index_t patch,const std::vector<T>& knots_u, const std::vector<T>& knots_v,bool updateBasis)
181{
182 std::vector<std::vector<T> >refineKnots(2);
183 refineKnots[0]=knots_u;
184 refineKnots[1]=knots_v;
185 basis(patch).insertKnots(refineKnots);
186 if(updateBasis)
187 {
188 repairPatches();
189 updateTopol();
190 }
191}
192
193template<short_t d, class T>
194void gsMPBESBSplineBasis<d,T>::refine(const index_t patch, gsMatrix<T> const & boxes, bool updateBasis)
195{
196 std::vector<std::vector<T> > refineKnots;
197 _boxToRefineKnots(patch,boxes,refineKnots);
198 refine(patch,refineKnots[0],refineKnots[1],updateBasis);
199}
200
201template<short_t d, class T>
202void gsMPBESBSplineBasis<d,T>::refineElements(const index_t patch, std::vector<index_t> const & boxes, bool updateBasis)
203{
204 gsMatrix<T> mat_boxes;
205 _boxesVectorToMatrix(boxes,mat_boxes);
206 refine(patch,mat_boxes,updateBasis);
207}
208
209template<short_t d, class T>
210void gsMPBESBSplineBasis<d,T>::refine_withCoefs(gsMatrix<T>& localCoef, const index_t patch,const std::vector<T>& knots_u, const std::vector<T>& knots_v,
211 bool updateBasis)
212{
213 std::vector<std::vector<T> >refineKnots(2);
214 refineKnots[0]=knots_u;
215 refineKnots[1]=knots_v;
216 std::vector<gsMatrix<T> *> coefs;
217 unsigned geoDim = localCoef.cols();
218 index_t start, end = -1;
219 for (size_t i = 0; i<nPatches(); ++i)
220 {
221 start=end+1;
222 end+=basis(i).size();
223 gsMatrix<T>* localMat = new gsMatrix<T>(end-start+1,geoDim);
224 *localMat << localCoef.block(start,0,end-start+1,geoDim);
225 coefs.push_back(localMat);
226 }
227 basis(patch).refine_withCoefs(*coefs[patch],refineKnots);
228 if(updateBasis)
229 repairPatches(coefs,patch);
230 unsigned totalSize=0;
231 for (size_t i = 0; i<nPatches(); ++i)
232 {
233 totalSize+=basis(i).size();
234 }
235 localCoef.resize(totalSize,geoDim);
236 end = -1;
237 for (size_t i = 0;i<nPatches();i++)
238 {
239 start=end+1;
240 end+=basis(i).size();
241 localCoef.block(start,0,end-start+1,geoDim) << *coefs[i];
242 }
243 if(updateBasis)
244 updateTopol();
245 freeAll(coefs);
246}
247
248template<short_t d, class T>
250 bool updateBasis)
251{
252 std::vector<std::vector<T> > refineKnots;
253 _boxToRefineKnots(patch,boxes,refineKnots);
254 refine_withCoefs(coefs,patch,refineKnots[0],refineKnots[1],updateBasis);
255}
256
257template<short_t d, class T>
258void gsMPBESBSplineBasis<d,T>::refineElements_withCoefs(gsMatrix<T>& coefs, const index_t patch, std::vector<index_t> const & boxes,
259 bool updateBasis)
260{
261 gsMatrix<T> mat_boxes;
262 _boxesVectorToMatrix(boxes,mat_boxes);
263 refine_withCoefs(coefs,patch,mat_boxes,updateBasis);
264}
265
266template<short_t d, class T>
267void gsMPBESBSplineBasis<d,T>::_boxesVectorToMatrix(const std::vector<index_t> & boxes, gsMatrix<T> & mat_boxes)
268{
269 GISMO_ASSERT( boxes.size() % (2 * d + 1 ) == 0, "The points did not define boxes properly. The boxes were not added to the basis.");
270 mat_boxes.resize(2,2*(boxes.size()/(2*d + 1)));
271 for( unsigned index_t i = 0; i < (boxes.size())/(2*d+1); i++)
272 {
273 mat_boxes(0,2*i)=boxes[(i*(2*d+1))+1];
274 mat_boxes(1,2*i)=boxes[(i*(2*d+1))+2];
275 mat_boxes(0,2*i+1)=boxes[(i*(2*d+1))+3];
276 mat_boxes(1,2*i+1)=boxes[(i*(2*d+1))+4];
277 }
278}
279
280template<short_t d, class T>
281void gsMPBESBSplineBasis<d,T>::_boxToRefineKnots(const index_t patch,gsMatrix<T> const & boxes,std::vector<std::vector<T> > & refineKnots)
282{
283 gsTensorBSplineBasis<d,T> & this_base=basis(patch);
284 GISMO_ASSERT( boxes.rows() == this_base.dim() , "Number of rows of refinement boxes must equal dimension of parameter space.");
285 GISMO_ASSERT( boxes.cols() % 2 == 0, "Refinement boxes must have even number of columns.");
286
287 const T tol = 0.000000001;
288 gsKnotVector<T> knots_u = this_base.component(0).knots();
289 gsKnotVector<T> knots_v = this_base.component(1).knots();
290 std::vector<T> ins_knots_u;
291 std::vector<T> ins_knots_v;
292 for(size_t i = 1;i < knots_u.size();i++)
293 if( knots_u[i]-knots_u[i-1] > tol)
294 {
295 T midpt = knots_u[i-1] + (knots_u[i]-knots_u[i-1])/2;
296 for( index_t j=0; j < boxes.cols(); j+=2 ) // loop over all boxes
297 {
298 if( boxes(0,j) < midpt && midpt < boxes(0,j+1) )
299 ins_knots_u.push_back(midpt);
300 }
301 }
302 for(size_t i = 1;i < knots_v.size();i++)
303 if( knots_v[i]-knots_v[i-1] > tol)
304 {
305 T midpt = knots_v[i-1] + (knots_v[i]-knots_v[i-1])/2;
306 for( index_t j=0; j < boxes.cols(); j+=2 ) // loop over all boxes
307 {
308 if( boxes(1,j) < midpt && midpt < boxes(1,j+1) )
309 ins_knots_v.push_back(midpt);
310 }
311 }
312 refineKnots.resize(2);
313 refineKnots[0]=ins_knots_u;
314 refineKnots[1]=ins_knots_v;
315}
316
317template<short_t d, class T>
318unsigned gsMPBESBSplineBasis<d,T>::getNrOfSpecialKnots(const gsKnotVector<T> kv,const std::vector<T>& new_knots,bool par,index_t distance)
319{
320 index_t specialKnots=0, index = 0, kv_index, nk_index;
321 for(index_t i=0;i<=distance+kv.degree();++i)
322 {
323 kv_index = par ? kv.size()-i-1 : i;
324 nk_index = par ? new_knots.size()-index-1 : index;
325 if(nk_index>=static_cast<index_t>(new_knots.size())||nk_index<0)
326 break;
327 if((!par&&new_knots[nk_index]<=kv.at(kv_index))||
328 (par&&new_knots[nk_index]>=kv.at(kv_index)))
329 {
330 specialKnots++;
331 index++;
332 }
333 }
334 return specialKnots;
335}
336
337template<short_t d, class T>
339 index_t startFromPatch)
340{
341 std::set<index_t> toCheck; //set of all indizes of patches, which have to be checked
342 if(startFromPatch==-1)
343 for(unsigned i = 0;i<m_bases.size();++i)
344 toCheck.insert(i);
345 else
346 toCheck.insert(startFromPatch);
347 unsigned curElement;
348 do
349 {
350 curElement=*(toCheck.begin());
351 toCheck.erase(curElement);
352 std::vector<std::vector<T> >refineKnots;
353 std::vector<index_t> checkPatches;
354 bool matched = _knotsMatchNeighbours(curElement,refineKnots,checkPatches);
355 if(!matched)
356 {
357 if(coefs[curElement]!=NULL)
358 basis(curElement).refine_withCoefs(*(coefs[curElement]),refineKnots);
359 else
360 basis(curElement).insertKnots(refineKnots);
361 }
362 for(unsigned i = 0;i<checkPatches.size();++i)
363 {
364 toCheck.insert(checkPatches[i]);
365 }
366 }while(!toCheck.empty());
367}
368
369template<short_t d, class T>
370bool gsMPBESBSplineBasis<d,T>::_knotsMatchNeighbours(index_t patch,std::vector<std::vector<T> >& knotsToInsert,
371 std::vector<index_t>& checkPatches,
372 T eps)
373{
374 typedef typename gsKnotVector<T>::iterator knotIter;
375 bool matched = true;
376 gsTensorBSplineBasis<d,T> & base = basis(patch);
377 patchSide ps, ps_neigh;
378 knotsToInsert.resize(2);
379 checkPatches.clear();
380 checkPatches.reserve(nPatches());
381 std::vector<bool> sideToCheck(4,false);
382 std::vector<index_t> neighbours(4,-1);
383 gsKnotVector<T> knots_neigh;
384 knotIter kv,kv_end,kvN,kvN_end;
385 for(boxSide side=boxSide::getFirst(2);side<boxSide::getEnd(2);++side)
386 {
387 ps=patchSide(patch,side);
388 if( !m_topol.getNeighbour(ps,ps_neigh) )
389 continue;
390 knots_neigh = basis(ps_neigh.patch).knots(1-(ps_neigh.direction()));
392 m_topol.getInterface(ps,bf);
393 if(!bf.dirOrientation(ps,1-ps.direction()))
394 knots_neigh.reverse();
395 kv = base.knots(1-side.direction()).begin();
396 kv_end = base.knots(1-side.direction()).end();
397 kvN = knots_neigh.begin();
398 kvN_end = knots_neigh.end();
399 neighbours[side-1]=ps_neigh.patch; //fill with the patch numbers of the neighbours
400 while(kv != kv_end || kvN != kvN_end)
401 {
402 if(math::abs(*kv-*kvN)<eps)
403 {
404 ++kv;
405 ++kvN;
406 }
407 else if(*kv < *kvN)
408 {
409 ++kv;
410 sideToCheck[side-1]=true; //set this side to true
411 }
412 else
413 {
414 knotsToInsert[1-side.direction()].push_back(*kvN);
415 matched = false;
416 ++kvN;
417 sideToCheck[side.opposite()-1]=true; //set the opposite side to true
418 }
419 }
420 }
421 //remove duplicates:
422 sort( knotsToInsert[0].begin(), knotsToInsert[0].end() );
423 knotsToInsert[0].erase( unique( knotsToInsert[0].begin(), knotsToInsert[0].end() ), knotsToInsert[0].end() );
424 sort( knotsToInsert[1].begin(), knotsToInsert[1].end() );
425 knotsToInsert[1].erase( unique( knotsToInsert[1].begin(), knotsToInsert[1].end() ), knotsToInsert[1].end() );
426 //find the patches we will have to check again:
427 for(unsigned i = 0;i<4;i++)
428 {
429 if(sideToCheck[i]&&neighbours[i]!=-1)
430 checkPatches.push_back(neighbours[i]);
431 }
432 return matched;
433}
434
435template<short_t d, class T>
436T gsMPBESBSplineBasis<d,T>::findParameter(patchSide const & ps,patchCorner const & pc,unsigned nrBasisFuncs) const
437{
438 if(nrBasisFuncs==0)
439 return 0.0;
440 gsKnotVector<T> knots = basis(ps.patch).knots(1-(ps.direction()));
441 gsVector<bool> pars;
442 pc.parameters_into(d,pars);
443 typename gsKnotVector<T>::uiterator iter;
444 if(!pars(1-ps.direction()))
445 knots.reverse();
446 iter = knots.ubegin();
447 iter+=nrBasisFuncs;
448 return *iter;
449}
450
451}
Struct which represents a certain side of a box.
Definition gsBoundary.h:85
static boxSide getEnd(short_t dim)
helper for iterating on sides of an n-dimensional box
Definition gsBoundary.h:176
static boxSide getFirst(short_t)
helper for iterating on sides of an n-dimensional box
Definition gsBoundary.h:160
short_t direction() const
Returns the parametric direction orthogonal to this side.
Definition gsBoundary.h:113
Class for representing a knot vector.
Definition gsKnotVector.h:80
uiterator ubegin() const
Returns unique iterator pointing to the beginning of the unique knots.
Definition gsKnotVector.hpp:235
T at(const size_t &i) const
Returns the value of the i - th knot (counted with repetitions).
Definition gsKnotVector.h:865
size_t size() const
Number of knots (including repetitions).
Definition gsKnotVector.h:242
iterator begin() const
Returns iterator pointing to the beginning of the repeated knots.
Definition gsKnotVector.hpp:117
std::vector< T >::const_iterator iterator
Definition gsKnotVector.h:97
void reverse()
Better directly use affineTransformTo.
Definition gsKnotVector.hpp:468
int degree() const
Returns the degree of the knot vector.
Definition gsKnotVector.hpp:700
internal::gsUKnotIterator< T > uiterator
Definition gsKnotVector.h:102
iterator end() const
Returns iterator pointing past the end of the repeated knots.
Definition gsKnotVector.hpp:124
A univariate Lagrange basis.
Definition gsMPBESBSplineBasis.h:33
void numActive_into(const index_t patch, const gsMatrix< T > &u, gsVector< index_t > &result) const
Returns the number of active (nonzero) basis functions at points u in result.
Definition gsMPBESBSplineBasis.hpp:170
gsMPBESBSplineBasis()
Default empty constructor.
Definition gsMPBESBSplineBasis.h:78
unsigned basisFunctionsOnSide(const patchSide &ps) const
Returns the amount of basis functions on a given side of a given patch.
Definition gsMPBESBSplineBasis.hpp:142
void refineElements_withCoefs(gsMatrix< T > &coefs, const index_t patch, std::vector< index_t > const &boxes, bool updateBasis=true)
Definition gsMPBESBSplineBasis.hpp:258
void _setMapping()
create a new mapping of the local basisfunctions
Definition gsMPBESBSplineBasis.hpp:132
void refineElements(const index_t patch, std::vector< index_t > const &boxes, bool updateBasis=true)
Definition gsMPBESBSplineBasis.hpp:202
T findParameter(patchSide const &ps, patchCorner const &pc, unsigned nrBasisFuncs) const
Definition gsMPBESBSplineBasis.hpp:436
void repairPatches(std::vector< gsMatrix< T > * > &coefs, index_t startFromPatch=-1)
Definition gsMPBESBSplineBasis.hpp:338
A univariate Lagrange basis.
Definition gsMPBESMapB2D.h:40
A matrix with arbitrary coefficient type and fixed or dynamic size.
Definition gsMatrix.h:41
A tensor product B-spline basis.
Definition gsTensorBSplineBasis.h:37
A vector with arbitrary coefficient type and fixed or dynamic size.
Definition gsVector.h:37
#define index_t
Definition gsConfig.h:32
#define GISMO_ERROR(message)
Definition gsDebug.h:118
#define GISMO_ASSERT(cond, message)
Definition gsDebug.h:89
Knot vector for B-splines.
Provides declaration of Basis abstract interface.
Provides declaration of Basis abstract interface.
Provides declaration of TensorBSplineBasis abstract interface.
The G+Smo namespace, containing all definitions for the library.
T distance(gsMatrix< T > const &A, gsMatrix< T > const &B, index_t i=0, index_t j=0, bool cols=false)
compute a distance between the point number in the set and the point number <j> in the set ; by def...
void freeAll(It begin, It end)
Frees all pointers in the range [begin end)
Definition gsMemory.h:312
Struct which represents an interface between two patches.
Definition gsBoundary.h:650
void parameters_into(index_t dim, gsVector< bool > &param) const
returns a vector of parameters describing the position of the corner
Definition gsBoundary.h:322
Struct which represents a certain corner of a patch.
Definition gsBoundary.h:393
Struct which represents a certain side of a patch.
Definition gsBoundary.h:232
index_t patch
The index of the patch.
Definition gsBoundary.h:234