G+Smo  25.01.0
Geometry + Simulation Modules
 
Loading...
Searching...
No Matches
gsModelingUtils.hpp
Go to the documentation of this file.
1
14#pragma once
15#include <iostream>
16#include <gsCore/gsLinearAlgebra.h>
17#include <gsNurbs/gsTensorBSpline.h>
18#include <gsIO/gsFileData.h>
19#include <gsAssembler/gsGaussRule.h>
20
21
22namespace gismo
23{
24
25
26template<class T>
27gsMatrix<T> * innerProduct( const gsBasis<T>& B1, const gsBasis<T>& B2)
28{
29 gsMatrix<T> * K = new gsMatrix<T>(B1.size(), B2.size() ) ;
30 K->setZero();
31
32 int nGauss = int( ceil( double(B1.degree(0) + B2.degree(0) + 1)/2 ) );
33 if (nGauss<1) nGauss=1;
34
35 gsGaussRule<T> QuRule(nGauss); // Reference Quadrature rule
36 gsMatrix<T> ngrid; // tensor Gauss nodes
37 gsVector<T> wgrid; // tensor Gauss weights
38 gsMatrix<index_t> act1, act2;
39 gsMatrix<T> ev1 , ev2;
40
41 typename gsBasis<T>::domainIter domIt = B1.makeDomainIterator();
42 for (; domIt->good(); domIt->next())
43 {
44 // Map the Quadrature rule to the element
45 QuRule.mapTo( domIt->lowerCorner(), domIt->upperCorner(), ngrid, wgrid );
46
47 B1.eval_into(ngrid,ev1);
48 B2.eval_into(ngrid,ev2);
49 B1.active_into(ngrid,act1);
50 B2.active_into(ngrid,act2);
51
52 for (index_t k=0; k!= ngrid.cols(); ++k)
53 for (index_t i=0; i!=act1.rows(); ++i)
54 for (index_t j=0; j!=act2.rows(); ++j)
55 (*K)( act1(i,k) , act2(j,k) ) += wgrid(k) * ev1(i,k) * ev2(j,k) ;
56 }
57
58 return K;
59}
60
61template<class T>
62gsMatrix<T> * innerProduct1( const gsBasis<T>& B1, const gsBasis<T>& B2)
63{
64 gsMatrix<T> * K = new gsMatrix<T>(B1.size(), B2.size() ) ;
65 K->setZero();
66
67 int nGauss = int( ceil( double(B1.degree(0)-1 + B2.degree(0)-1 + 1)/2 ) );
68 if (nGauss<1) nGauss=1;
69
70 gsGaussRule<T> QuRule(nGauss); // Reference Quadrature rule
71 gsMatrix<T> ngrid; // tensor Gauss nodes
72 gsVector<T> wgrid; // tensor Gauss weights
73 gsMatrix<index_t> act1, act2;
74 gsMatrix<T> ev1 , ev2;
75
76 typename gsBasis<T>::domainIter domIt = B1.makeDomainIterator();
77 for (; domIt->good(); domIt->next())
78 {
79 // Map the Quadrature rule to the element
80 QuRule.mapTo( domIt->lowerCorner(), domIt->upperCorner(), ngrid, wgrid );
81
82 B1.deriv_into(ngrid,ev1);
83 B2.deriv_into(ngrid,ev2);
84 B1.active_into(ngrid,act1);
85 B2.active_into(ngrid,act2);
86
87 for (index_t k=0; k!= ngrid.cols(); ++k)
88 for (index_t i=0; i!=act1.rows(); ++i)
89 for (index_t j=0; j!=act2.rows(); ++j)
90 (*K)( act1(i,k) , act2(j,k) ) += wgrid(k) * ev1(i,k) * ev2(j,k) ;
91
92 }
93
94 return K;
95}
96
97template<class T>
98gsMatrix<T> * innerProduct2( const gsBasis<T>& B1, const gsBasis<T>& B2)
99{
100 gsMatrix<T> * K = new gsMatrix<T>(B1.size(), B2.size() ) ;
101 K->setZero();
102
103 int nGauss = int( ceil( double(B1.degree(0)-2 + B2.degree(0)-2 + 1)/2 ) );
104 if (nGauss<1) nGauss=1;
105
106 gsGaussRule<T> QuRule(nGauss); // Reference Quadrature rule
107 gsMatrix<T> ngrid; // tensor Gauss nodes
108 gsVector<T> wgrid; // tensor Gauss weights
109 gsMatrix<index_t> act1, act2;
110 gsMatrix<T> ev1 , ev2;
111
112 typename gsBasis<T>::domainIter domIt = B1.makeDomainIterator();
113 for (; domIt->good(); domIt->next())
114 {
115 // Map the Quadrature rule to the element
116 QuRule.mapTo( domIt->lowerCorner(), domIt->upperCorner(), ngrid, wgrid );
117
118 B1.deriv2_into(ngrid,ev1);
119 B2.deriv2_into(ngrid,ev2);
120 B1.active_into(ngrid,act1);
121 B2.active_into(ngrid,act2);
122
123 for (index_t k=0; k!= ngrid.cols(); ++k)
124 for (index_t i=0; i!=act1.rows(); ++i)
125 for (index_t j=0; j!=act2.rows(); ++j)
126 (*K)( act1(i,k) , act2(j,k) ) += wgrid(k) * ev1(i,k) * ev2(j,k) ;
127 }
128
129 return K;
130}
131
132
134template <class T>
135gsVector<T> vectorIntersect(gsVector<T> const & tangent1,
136 gsVector<T> const & tangent2,
137 gsMatrix<T> const & Vert1,
138 gsMatrix<T> const & Vert2)
139{
140 gsVector3d<T> abc1;
141 gsVector3d<T> abc2;
142 abc1(0) = tangent1(1);
143 abc1(1) = -tangent1(0);
144 abc1(2) = -( abc1(0)*Vert1(0,0) + abc1(1)*Vert1(0,1) );
145 abc2(0) = tangent2(1);
146 abc2(1) = -tangent2(0);
147 abc2(2) = -( abc2(0)*Vert2(0,0) + abc2(1)*Vert2(0,1) );
148 gsVector<T> unknown(2);
149 T detMatrixab = abc1(0)*abc2(1)-abc1(1)*abc2(0);
150 if (detMatrixab!=0)
151 {
152 unknown(0) = -(1/detMatrixab)*( abc2(1)*abc1(2)-abc1(1)*abc2(2) );
153 unknown(1) = -(1/detMatrixab)*( -abc2(0)*abc1(2)+abc1(0)*abc2(2) );
154 }
155 else
156 {
157 unknown(0) = .5*Vert1(0) + .5*Vert2(0);
158 unknown(1) = .5*Vert1(1) + .5*Vert2(1);
159 };
160 return unknown;
161}
162
164template <class T>
165T conditionedAngle(gsVector3d<T> vec1, gsVector3d<T> vec2)
166{ T dotp = vec1.dot(vec2)/( vec1.norm()*vec2.norm() );
167 if ( dotp<-1 ) {
168 dotp = -1;
169 gsWarn<<"gsModelingUtils: truncations done for std::acos \n";
170 return EIGEN_PI;}
171 if ( dotp>1 ) {dotp = 1;
172 gsWarn<<"gsModelingUtils: truncations done for std::acos \n";
173 return 0;}
174
175 T ag = math::acos(dotp);
176 return ag;
177}
178
181template <class T>
182T conditionedAngle(gsVector3d<T> vec1, gsVector3d<T> vec2, gsVector3d<T> normal)
183{
184 T ag = conditionedAngle<T>(vec1,vec2);
185 T cag = ( normal.dot( vec1.cross( vec2 ) ) >= 0 ) ? ag : (T)(2*EIGEN_PI)-ag;
186 return cag;
187}
188
189
192template<class T>
193gsVector<T> criticalPointOfQuadratic(gsMatrix<T> & A, gsMatrix<T> & C, gsVector<T> & d)
194{
195 index_t n = A.rows();
196 index_t m = d.rows();
197
198 gsMatrix<T> combined(n + m, n + m);
199 combined.block(0, 0, n, n) = 2*A;
200 combined.block(0, n, n, m) = C.transpose();
201 combined.block(n, 0, m, n) = C;
202 combined.block(n, n, m, m).setZero();
203
204 gsVector<T> longD(n + m);
205 longD.head(n).setZero();
206 longD.tail(m) = d;
207
208 gsMatrix<T> result = combined.fullPivLu().solve(longD);
209 return result.block(0, 0, n, 1);
210}
211
214template<class T>
215gsMatrix<T> criticalPointOfQuadratic(gsMatrix<T> const & A, gsMatrix<T> const & b,
216 gsMatrix<T> const & C, gsMatrix<T> const & d)
217{
218 index_t n = A.rows(); // dimension of X
219 index_t m = d.rows(); // number of exact constraints
220 assert(m<=n); // if not, the problem is ill defined
221 assert(A.cols()==n); // A must be a square matrix
222 assert(C.cols()==n);
223 assert(d.cols()==1 && C.rows()==d.rows());
224 gsMatrix<T> bt;
225 if (b.rows()!=1) bt=b.transpose(); else bt=b;
226
227 gsMatrix<T> combined(n + m, n + m);
228 combined.block(0, 0, n, n) = 2*A;
229 combined.block(0, n, n, m) = C.transpose();
230 combined.block(n, 0, m, n) = C;
231 combined.block(n, n, m, m).setZero();
232
233 gsVector<T> longD(n + m);
234 longD.head(n) = -bt.transpose();
235 longD.tail(m) = d;
236
237 gsMatrix<T> result = combined.fullPivLu().solve(longD);
238 return result.block(0, 0, n, 1);
239}
240
243template<class T>
244gsMatrix<T> criticalPointOfQuadratic(gsMatrix<T> const & A, gsMatrix<T> const & C,
245 gsMatrix<T> const & d)
246{
247 gsMatrix<T> b(1,A.cols()); b.setZero();
248 return criticalPointOfQuadratic<T>(A,b,C,d);
249}
250
252template<class T>
253gsMatrix<T> optQuadratic(gsMatrix<T> const & A, gsMatrix<T> const & b,
254 gsMatrix<T> const & C, gsMatrix<T> const & d)
255{
256 return criticalPointOfQuadratic<T>( (A.transpose())*A, (-2)*( A.transpose() )*b, C, d);
257}
258
260template<class T>
261gsMatrix<T> optQuadratic(gsMatrix<T> const & A1, gsMatrix<T> const & b1,
262 T const & w1, gsMatrix<T> const & A2,
263 gsMatrix<T> const & b2, T const & w2,
264 gsMatrix<T> const & C, gsMatrix<T> const & d)
265{
266 return criticalPointOfQuadratic<T>( w1*(A1.transpose())*A1 + w2*(A2.transpose())*A2,
267 (-2)*w1*( A1.transpose() )*b1 + (-2)*w2*( A2.transpose() )*b2, C, d);
268}
269
272template<class T>
273gsMatrix<T> optQuadratic(gsMatrix<T> const & A1, gsMatrix<T> const & b1,
274 T const & w1, gsMatrix<T> const & A2,
275 gsMatrix<T> const & b2, T const & w2,
276 gsMatrix<T> const & C, gsMatrix<T> const & d,
277 T const & w3, gsMatrix<T> const & Q)
278{
279 return criticalPointOfQuadratic<T>( w1*(A1.transpose())*A1 + w2*(A2.transpose())*A2 + w3*Q,
280 (-2)*w1*( A1.transpose() )*b1 + (-2)*w2*( A2.transpose() )*b2, C, d);
281}
282
284template <class T>
285gsMatrix<T> flipLR(const gsMatrix<T> & mat)
286{
287 return mat.rowwise().reverse();
288
289 // size_t const ncol = mat.cols();
290 // gsMatrix<T> nMat(mat.rows(),ncol);
291 // for (size_t i=0; i<= ncol-1; i++)
292 // {
293 // nMat.col(i) = mat.col(ncol-1-i);
294 // };
295 // return nMat;
296}
297
300template <class T>
301gsMatrix<T> crossNorm2Mat(gsMatrix<T> const & mat1,gsMatrix<T> const & mat2)
302{
303 assert(mat1.rows()==3 && mat2.rows()==3);
304 assert(mat1.cols()==mat2.cols());
305 size_t const nr = mat1.rows();
306 size_t const nc = mat1.cols();
307 gsMatrix<T> rcross(nr,mat1.cols()); rcross.setZero();
308 gsMatrix<T> tem(3,1);
309 for (size_t i=0; i<= nc-1; i++)
310 {
311 tem = ( gsVector3d<T>( mat1.col(i) ) ).cross( gsVector3d<T>( mat2.col(i) ) );
312 rcross.col(i) = tem / tem.norm();
313 };
314 return rcross;
315}
316
318template <class T>
319void addConstraints(gsMatrix<T> const & C1, gsMatrix<T> const & d1,
320 gsMatrix<T> const & C2, gsMatrix<T> const & d2,
321 gsMatrix<T> & C, gsMatrix<T> & d)
322{
323 int nr1 = C1.rows();
324 int nc1 = C1.cols();
325 int nr2 = C2.rows();
326
327 assert(nc1== C2.cols() );
328 assert(nr1==d1.rows());
329 assert(nr2==d2.rows());
330 C.resize(nr1+nr2,nc1);
331 C.block(0,0,nr1,nc1) = C1;
332 C.block(nr1,0,nr2,nc1) = C2;
333 d.resize(nr1+nr2,1);
334 d.block(0,0,nr1,1) = d1;
335 d.block(nr1,0,nr2,1) = d2;
336}
337
339template <class T>
340gsMatrix<T> convert2Zero(gsMatrix<T> const & mat)
341{
342 T eps=2.220446049250313e-16;
343
344 gsMatrix<T> matc = mat;
345 int n1 = mat.rows();
346 int n2 = mat.cols();
347 for (int i=0;i!=n1;i++)
348 {
349 for (int j=0;j!=n2;j++)
350 {
351 if (math::abs(mat(i,j))< eps) matc(i,j)=0.;
352 }
353 }
354 return matc;
355}
356
358template <class T>
359void removeCol(gsMatrix<T> & mat, int const & removeEnds, int const & nPoints)
360{
361 assert(removeEnds==1 || removeEnds==2);
362 int nPeriod = mat.cols()/nPoints;
363 assert( nPeriod*nPoints == mat.cols() );
364 int ind1,ind2;
365 if (removeEnds==1)
366 {
367 for (int i=nPeriod-1;i>=0;i--)
368 {
369 ind2 = i*nPoints + nPoints-1;
370 mat.removeCol(ind2);
371 };
372 };
373 if (removeEnds==2)
374 {
375 for (int i=nPeriod-1;i>=0;i--)
376 {
377 ind2 = i*nPoints + nPoints-1;
378 ind1 = i*nPoints ; //+0
379 mat.removeCol(ind2);
380 mat.removeCol(ind1);
381 };
382 };
383}
384
391template <class T>
392gsBSpline<T> gsInterpolate(gsKnotVector<T> & kv,const gsMatrix<T> & preImage,
393 const gsMatrix<T> & image,
394 const gsMatrix<T> & preNormal,const gsMatrix<T> & normal,
395 const gsMatrix<T> & preImageApp,const gsMatrix<T> & imageApp,
396 T const & w_reg,T const & w_app,
397 gsMatrix<T> &outPointResiduals, gsMatrix<T> &outNormalResiduals)
398{
399 const int ntcp = kv.size()-kv.degree()-1;
400 gsMatrix<T> tcp (ntcp, 2);
401
402 // Quadratic forms which (approximately) constitute the beam strain energy
403 gsBSplineBasis<T> bs(kv);
404 gsMatrix<T> *Q = innerProduct2(bs, bs);
405
406 // Exact constraints: point interpolation
407 short_t dimPI = 1; // dimension of space of preImage, TODO: put dimPI, dimI to template<dimPI,...
408 short_t dimI = 2; // dimension of space of image
409 GISMO_UNUSED(dimPI);
410 int nip = image.cols(); // number of interpolating points
411 int nn=normal.cols(); // number of prescribed normals
412 gsMatrix<T> Nu, dNu, dNu_nm, NuApp;
413 //--
414 GISMO_ASSERT(dimPI==1 && dimI==2," "); // can be easily extended for other dimensions
415 Nu = bs.eval(preImage.row(0)); // u-BSplines
416 dNu = bs.deriv(preImage.row(0)); // u-BSplines
417 Nu.transposeInPlace();
418 dNu.transposeInPlace();
419 dNu_nm = bs.deriv(preNormal.row(0)); // u-BSplines for normals
420 dNu_nm.transposeInPlace();
421 gsMatrix<T> AdN = normal.row(0).asDiagonal() * dNu_nm;
422 gsMatrix<T> BdN = normal.row(1).asDiagonal() * dNu_nm;
423
424 // Approximate constraints
425 NuApp = bs.eval(preImageApp.row(0));
426 NuApp.transposeInPlace();
427 gsMatrix<T> X0 = imageApp.row(0);
428 gsMatrix<T> Y0 = imageApp.row(1);
429
430 //-- resulting Saddle point linear System
431 int nss = dimI*ntcp + dimI*nip + nn;
432 gsMatrix<T> Ass(nss,nss);
433 gsMatrix<T> bss(nss,1);
434 Ass.setZero(); bss.setZero();
435 gsMatrix<T> Hess = w_reg*(T)(2)*(*Q) + w_app*(T)(2)*(NuApp.transpose())* NuApp;
436 //--- row 0
437 Ass.block(0,0,ntcp,ntcp) = Hess;
438 Ass.block(0,2*ntcp,ntcp,nip) = Nu.transpose();
439 Ass.block(0,2*ntcp+2*nip,ntcp,nn) = AdN.transpose();
440 bss.block(0,0,ntcp,1) = w_app*(T)(2)*NuApp.transpose()*X0.transpose();
441 //--- row 1
442 Ass.block(ntcp,ntcp,ntcp,ntcp) = Hess;
443 Ass.block(ntcp,2*ntcp+nip,ntcp,nip) = Nu.transpose();
444 Ass.block(ntcp,2*ntcp+2*nip,ntcp,nn) = BdN.transpose();
445 bss.block(ntcp,0,ntcp,1) = w_app*(T)(2)*NuApp.transpose()*Y0.transpose();
446 //--- row 2
447 Ass.block(2*ntcp,0,nip,ntcp) = Nu;
448 bss.block(2*ntcp,0,nip,1) = (image.row(0)).transpose();
449 //--- row 3
450 Ass.block(2*ntcp+nip,ntcp,nip,ntcp) = Nu;
451 bss.block(2*ntcp+nip,0,nip,1) = (image.row(1)).transpose();
452 //--- row 4
453 Ass.block(2*ntcp+2*nip,0,nn,ntcp) = AdN;
454 Ass.block(2*ntcp+2*nip,ntcp,nn,ntcp) = BdN;
455
456 gsMatrix<T> result = Ass.fullPivLu().solve(bss);
457 tcp.col(0) = result.block(0 , 0, ntcp, 1);
458 tcp.col(1) = result.block(ntcp, 0, ntcp, 1);
459
460// gsDebug<< " parameterRange: \n"<< *trimLoop[sourceID]->basis().parameterRange()<<"\n";
461// gsDebug<<" Ass: \n"<<Ass<<"\n";
462// gsDebug<<" bss: \n"<<bss<<"\n";
463// gsDebug<<" result: \n"<<result<<"\n";
464// gsDebug<<" Q: \n"<< *Q<<"\n";
465// gsDebug<<" preimage: \n"<< preImage<<"\n";
466// gsDebug<<" prenormal: \n"<< preNormal<<"\n";
467// gsDebug<<" image: \n"<< image<<"\n";
468// gsDebug<<" normal: \n"<< normal<<"\n";
469// gsDebug<<" Nu: \n"<< *Nu<<"\n";
470// gsDebug<<" dNu: \n"<< *dNu<<"\n";
471// gsDebug<<" AdN: \n"<< AdN<<"\n";
472// gsDebug<<" BdN: \n"<< BdN<<"\n";
473// gsDebug<<" tcp: \n"<< tcp<<"\n";
474// gsDebug<<" preimageApp: \n"<< preImageApp<<"\n";
475// gsDebug<<" imageApp: \n"<< imageApp<<"\n";
476// gsDebug<<" residual of app x constraints: \n"<< *NuApp*tcp.col(0)-imageApp.row(0).transpose()<<std::endl;
477// gsDebug<<" residual of app y constraints: \n"<< *NuApp*tcp.col(1)-imageApp.row(1).transpose()<<std::endl;
478// gsDebug<<" residual of normal constraints: \n"<< AdN*tcp.col(0)+BdN*tcp.col(1)<<std::endl;
479
480 outPointResiduals = (NuApp * tcp).transpose() - imageApp;
481 outNormalResiduals = AdN * tcp.col(0) + BdN * tcp.col(1);
482 //gsDebug << std::flush;
483
484 delete Q;
485
486 gsBSpline<T> tcurve(kv, give(tcp));
487
488 return tcurve;
489}
490
491
503template<class T>
504typename gsTensorBSpline<2,T>::Ptr gsInterpolateSurface(
505 const gsMatrix<T> &exactPoints, const gsMatrix<T> &exactValues,
506 const gsMatrix<T> &appxPointsEdges, const gsMatrix<T> &appxValuesEdges,
507 const gsMatrix<T> &appxPointsInt, const gsMatrix<T> &appxValuesInt,
508 const gsMatrix<T> &appxNormalPoints, const gsMatrix<T> &appxNormals,
509 T wEdge, T wInt, T wNormal, T wReg,
510 const gsKnotVector<T> &kv1, const gsKnotVector<T> &kv2,
511 bool force_normal
512 )
513{
514 GISMO_ASSERT(exactPoints.rows() == 2 && exactValues.rows() == 3, "Matrix input has incorrect dimension");
515 GISMO_ASSERT(exactPoints.cols() == exactValues.cols(), "Matrix input has incorrect dimension");
516 GISMO_ASSERT(appxPointsEdges.rows() == 2 && appxValuesEdges.rows() == 3, "Matrix input has incorrect dimension");
517 GISMO_ASSERT(appxPointsEdges.cols() == appxValuesEdges.cols(), "Matrix input has incorrect dimension");
518 GISMO_ASSERT(appxPointsInt.rows() == 2 && appxValuesInt.rows() == 3, "Matrix input has incorrect dimension");
519 GISMO_ASSERT(appxPointsInt.cols() == appxValuesInt.cols(), "Matrix input has incorrect dimension");
520 GISMO_ASSERT(appxNormalPoints.rows() == 2 && appxNormals.rows() == 3, "Matrix input has incorrect dimension");
521 GISMO_ASSERT(appxNormalPoints.cols() == appxNormals.cols(), "Matrix input has incorrect dimension");
522
523
524 int const patchDeg1 = kv1.degree();
525 int const patchDeg2 = kv2.degree();
526 int const n1 = kv1.size() - patchDeg1 - 1;
527 int const n2 = kv2.size() - patchDeg2 - 1;
528
529 gsBSplineBasis<T> bs1(kv1);
530 gsBSplineBasis<T> bs2(kv2);
531 gsMatrix<T> Nu, Nv, dNu, dNv;
532 gsMatrix<T> R;
533 int npts;
534
535 //Assemble Exact constraints for corners
536 gsMatrix<T> ident1(n1, n1);
537 gsMatrix<T> ident2(n2, n2);
538 ident1.setIdentity();
539 ident2.setIdentity();
540 Nu = bs1.evalFunc(exactPoints.row(0), ident1); // u-BSplines
541 Nv = bs2.evalFunc(exactPoints.row(1), ident2); // v-BSplines
542 npts = exactPoints.cols();
543 gsMatrix<T> Ecor(3*npts,3*n1*n2);
544 gsMatrix<T> ecor(3*npts,1);
545 Ecor.setZero();ecor.setZero();
546 R = Nv.khatriRao(Nu); // R = M tensors N
547 Ecor.block(0,0,npts,n1*n2) = R.transpose();
548 Ecor.block(npts,n1*n2,npts,n1*n2) = R.transpose();
549 Ecor.block(2*npts,2*n1*n2,npts,n1*n2) = R.transpose();
550 ecor.block(0,0,npts,1) = (exactValues.row(0)).transpose();
551 ecor.block(npts,0,npts,1) = (exactValues.row(1)).transpose();
552 ecor.block(2*npts,0,npts,1) = (exactValues.row(2)).transpose();
553
554 //Assemble Approximative constraints for inner points on edges
555 Nu = bs1.evalFunc(appxPointsEdges.row(0), ident1); // u-BSplines
556 Nv = bs2.evalFunc(appxPointsEdges.row(1), ident2); // v-BSplines
557 npts = appxPointsEdges.cols();
558 gsMatrix<T> AappEdge(3*npts,3*n1*n2);
559 gsMatrix<T> bappEdge(3*npts,1);
560 AappEdge.setZero();bappEdge.setZero();
561 R = Nv.khatriRao(Nu); // R = M tensors N
562 AappEdge.block(0,0,npts,n1*n2) = R.transpose();
563 AappEdge.block(npts,n1*n2,npts,n1*n2) = R.transpose();
564 AappEdge.block(2*npts,2*n1*n2,npts,n1*n2) = R.transpose();
565 bappEdge.block(0,0,npts,1) = (appxValuesEdges.row(0)).transpose();
566 bappEdge.block(npts,0,npts,1) = (appxValuesEdges.row(1)).transpose();
567 bappEdge.block(2*npts,0,npts,1) = (appxValuesEdges.row(2)).transpose();
568
569 //Assemble Approximate constraints for interior
570 Nu = bs1.evalFunc(appxPointsInt.row(0), ident1); // u-BSplines
571 Nv = bs2.evalFunc(appxPointsInt.row(1), ident2); // v-BSplines
572 npts = appxPointsInt.cols();
573 gsMatrix<T> AappInt(3*npts,3*n1*n2);
574 gsMatrix<T> bappInt(3*npts,1);
575 AappInt.setZero();bappInt.setZero();
576 R = Nv.khatriRao(Nu); // R = M tensors N
577 AappInt.block(0,0,npts,n1*n2) = R.transpose();
578 AappInt.block(npts,n1*n2,npts,n1*n2) = R.transpose();
579 AappInt.block(2*npts,2*n1*n2,npts,n1*n2) = R.transpose();
580 bappInt.block(0,0,npts,1) = (appxValuesInt.row(0)).transpose();
581 bappInt.block(npts,0,npts,1) = (appxValuesInt.row(1)).transpose();
582 bappInt.block(2*npts,0,npts,1) = (appxValuesInt.row(2)).transpose();
583
584 //Assemble Approximative constraints for normals
585 Nu = bs1.evalFunc(appxNormalPoints.row(0), ident1); // u-BSplines
586 dNu = bs1.derivFunc(appxNormalPoints.row(0), ident1);
587 Nv = bs2.evalFunc(appxNormalPoints.row(1), ident2); // v-BSplines
588 dNv = bs2.derivFunc(appxNormalPoints.row(1), ident2);
589 gsMatrix<T> Nx,Ny,Nz; // x, y, z components of normals
590 gsMatrix<T> trNormals = appxNormals.transpose();
591 Nx = trNormals.col(0);
592 Ny = trNormals.col(1);
593 Nz = trNormals.col(2);
594 if (force_normal==true) { Nx.setZero(); Ny.setZero();Nz.setOnes();} //TODO: do this automatically
595
596 gsMatrix<T> dRdu = Nv.khatriRao(dNu); // R = M tensors N
597 gsMatrix<T> dRdv = dNv.khatriRao(Nu); // R = M tensors N
598 dRdu.transposeInPlace();
599 dRdv.transposeInPlace();
600 int nnor = Nx.rows(); // number of normals
601 gsMatrix<T> Anor(2*nnor,3*n1*n2),bnor(2*nnor,1);
602 Anor.setZero();
603 bnor.setZero();
604 Anor.block(0,0,nnor,n1*n2) = Nx.asDiagonal()*dRdu; // sigma_u . normal=0, x part
605 Anor.block(0,n1*n2,nnor,n1*n2) = Ny.asDiagonal()*dRdu; // sigma_u . normal=0, y part
606 Anor.block(0,2*n1*n2,nnor,n1*n2) = Nz.asDiagonal()*dRdu; // sigma_u . normal=0, z part
607 Anor.block(nnor,0,nnor,n1*n2) = Nx.asDiagonal()*dRdv; // sigma_v . normal=0, x part
608 Anor.block(nnor,n1*n2,nnor,n1*n2) = Ny.asDiagonal()*dRdv; // sigma_v . normal=0, y part
609 Anor.block(nnor,2*n1*n2,nnor,n1*n2) = Nz.asDiagonal()*dRdv; // sigma_v . normal=0, z part
610
611 // Quadratic forms which constitute the plate bending energy
612 gsMatrix<T> * M = innerProduct(bs1, bs1);
613 gsMatrix<T> * M1 = innerProduct1(bs1, bs1);
614 gsMatrix<T> * M2 = innerProduct2(bs1, bs1);
615 gsMatrix<T> *N, *N1, *N2;
616 if (kv1==kv2) { N = M; N1 = M1; N2 = M2; } else
617 {
618 N = innerProduct(bs2, bs2);
619 N1 = innerProduct1(bs2, bs2);
620 N2 = innerProduct2(bs2, bs2);
621 };
622 gsMatrix<T> Q1D = N->kron(*M2) + 2*N1->kron(*M1)+N2->kron(*M);
623 gsMatrix<T> Q(3*n1*n2,3*n1*n2);
624 Q.setZero();
625 Q.block(0,0,n1*n2,n1*n2) = Q1D;
626 Q.block(n1*n2,n1*n2,n1*n2,n1*n2) = Q1D;
627 Q.block(2*n1*n2,2*n1*n2,n1*n2,n1*n2) = Q1D;
628
629 // now solve the cps
630 gsMatrix<T> coefA = wEdge*(AappEdge.transpose())*AappEdge + wInt*(AappInt.transpose())*AappInt + wNormal*(Anor.transpose())*Anor + wReg*Q;
631 gsMatrix<T> coefb = (T)(-2)*wEdge*(AappEdge.transpose())*bappEdge + (T)(-2)*wInt*(AappInt.transpose())*bappInt + (T)(-2)*wNormal*(Anor.transpose())*bnor;
632
633 gsMatrix<T> cp = criticalPointOfQuadratic(
634 coefA,
635 coefb,
636 Ecor,
637 ecor);
638
639 cp.resize( cp.rows() / 3, 3);
640
641 typename gsTensorBSpline<2,T>::Ptr master(new gsTensorBSpline<2,T>( kv1, kv2, give(cp) ));
642
643 delete M; delete M1; delete M2;
644 if (kv1!=kv2) {delete N; delete N1; delete N2;}
645
646 // check that the spline actually satisfies the exact constraints
647 for(index_t idxConstr = 0; idxConstr < exactPoints.cols(); idxConstr++)
648 {
649 gsMatrix<T> pt = exactPoints.col(idxConstr);
650 gsMatrix<T> expectVal = exactValues.col(idxConstr);
651 gsMatrix<T> surfVal = master->eval(pt);
652 GISMO_ASSERT((expectVal - surfVal).norm() < 0.0001, "Fit surface did not satisfy exact constraints");
653 }
654
655 return master;
656}
657
663 template<class T>
664 void threeOnDiag(const gsSparseMatrix<T>& block,
665 gsSparseMatrix<T>& result) //const
666 {
667
668 index_t brows = block.rows();
669 index_t bcols = block.cols();
670
671 result = gsSparseMatrix<T>(3*brows,3*bcols);
672 result.reservePerColumn(block.nonZeros() / bcols);
673
674 for(index_t j = 0; j != bcols; j++)
675 {
676 for (auto it = block.begin(j); it; ++it)
677 {
678 result.insertTo(it.row(), j, it.value());
679 result.insertTo(brows+it.row(), bcols+j, it.value());
680 result.insertTo(2*brows+it.row(), 2*bcols+j, it.value());
681 }
682 }
683 result.makeCompressed();
684 }
685
687 template<class T>
688 void scalePoints(const gsMatrix<T> & xyz,
689 gsMatrix<T> & points)
690 {
691 T p_min = xyz.minCoeff(),
692 p_max = xyz.maxCoeff();
693 T den = p_max - p_min;
694
695 points.resize(xyz.rows(), xyz.cols());
696 points = (1/den)*(xyz - p_min * gsMatrix<T>::Ones(xyz.rows(), xyz.cols()));
697 }
698
699
701 template<class T>
702 T scaleTo01(T tMin, T t, T tMax)
703 {
704 return (t - tMin) / (tMax - tMin);
705 }
706
708 template <class T>
709 void scaleTo01(real_t tMin, gsMatrix<T>& mT, real_t tMax)
710 {
711 for(index_t i=0; i<mT.rows(); i++)
712 for(index_t j=0; j<mT.cols(); j++)
713 mT(i, j) = scaleTo01(tMin, mT(i, j), tMax);
714 }
715
717 template <class T>
718 void scaleTo01(gsMatrix<T>& xyz, bool verbose)
719 {
720 T xyzMin = xyz.minCoeff();
721 T xyzMax = xyz.maxCoeff();
722 scaleTo01(xyzMin, xyz, xyzMax);
723
724 if(verbose)
725 gsInfo << std::setprecision(15)
726 << "Scaling TO [0, 1]^3 as follows.\n"
727 << "min: " << xyzMin << "\n"
728 << "max: " << xyzMax << "\n"
729 << "scale: " << T(1.0) / (xyzMax - xyzMin) << std::endl;
730 }
731
732
734 template <class T>
735 T scaleFrom01(T tMin, T t, T tMax)
736 {
737 return (tMax - tMin) * t + tMin;
738 }
739
740
742 template <class T>
743 void scaleFrom01(T tMin, gsMatrix<T>& mT, T tMax, bool verbose)
744 {
745 for(index_t i=0; i<mT.rows(); i++)
746 for(index_t j=0; j<mT.cols(); j++)
747 mT(i, j) = scaleFrom01(tMin, mT(i, j), tMax);
748
749 if(verbose)
750 gsInfo << "Scaling FROM [0, 1]^3.\n"
751 << "inverted scale: " << T(1.0) / (tMax - tMin) << std::scientific << std::endl;
752 }
753
754
756 template <class T>
757 void scaleFrom01(T tMin, gsGeometry<T>& geo, T tMax, bool verbose)
758 {
759 gsMatrix<T> coefs = geo.coefs();
760 scaleFrom01(tMin, coefs, tMax, verbose);
761 geo.coefs() = coefs;
762 }
763
764
766 template <class T>
767 void scaleGeo(const std::string& fin,
768 const std::string& fout,
769 T tMin,
770 T tMax,
771 bool verbose)
772 {
773 gsFileData<T> fdIn(fin);
774 typename gsGeometry<T>::uPtr geo = fdIn.template getFirst< gsGeometry<T> >();
775 scaleFrom01(tMin, *geo.get(), tMax, verbose);
776 gsFileData<T> fdOut;
777 fdOut << *geo.get();
778 fdOut.dump(fout);
779 }
780
781
783 template <class T>
784 void scalePts(const std::string& fin,
785 const std::string& fout,
786 index_t uvIdIn,
787 index_t uvIdOut,
788 index_t xyzIdIn,
789 index_t xyzIdOut,
790 T tMin,
791 T tMax,
792 bool verbose)
793 {
794 // Reading the inputs.
795 gsFileData<T> fdIn(fin);
796 gsMatrix<T> xyz, uv;
797 fdIn.template getId<gsMatrix<T>>(uvIdIn, uv);
798 fdIn.template getId<gsMatrix<T>>(xyzIdIn, xyz);
799
800 // Scaling the data.
801 if(tMax < tMin) // defaults, i.e., scaling down
802 scaleTo01(xyz, verbose);
803 else
804 scaleFrom01(tMin, xyz, tMax, verbose);
805
806 // Writing the outputs
807 gsFileData<T> fdOut;
808
809 if(uvIdOut < xyzIdOut)
810 {
811 fdOut << uv;
812 fdOut << xyz;
813 }
814 else
815 {
816 fdOut << xyz;
817 fdOut << uv;
818 }
819 fdOut.dump(fout);
820 }
821
822
831 template<class T>
832 void sortPointCloud(gsMatrix<T> & parameters,
833 gsMatrix<T> & points,
834 std::vector<index_t> & corners)
835 {
836 // The following matrices and vectors store the parameters and points values and indeces.
837 // There is no need to store these information, we could also use only one matrix and 1 std::vector and overwirte them each time.
838 gsMatrix<T> uv_interiors, uv_south, uv_east, uv_north, uv_west;
839 gsMatrix<T> p_interiors, p_south, p_east, p_north, p_west;
840 std::vector<index_t> interiors, b_west, b_east, b_south, b_north;
841
842 // Determine the parameter domain by mi/max of parameter values
843 T u_min = parameters.row(0).minCoeff(),
844 u_max = parameters.row(0).maxCoeff(),
845 v_min = parameters.row(1).minCoeff(),
846 v_max = parameters.row(1).maxCoeff();
847
848 gsVector<T> curr_point(2,1);
849 for(index_t i=0; i < parameters.cols(); i++)
850 {
851 curr_point = parameters.col(i);
852 if( (u_min < curr_point(0)) && (curr_point(0) < u_max) && (v_min < curr_point(1)) && (curr_point(1) < v_max) )
853 interiors.push_back(i);
854 else // not interior point
855 {
856 if( (math::abs(curr_point(0) - u_min) < 1e-15) && (curr_point(1) > v_min) )
857 b_west.push_back(i);//west edge
858 else if( (math::abs(curr_point(0) - u_max) < 1e-15) && curr_point(1) < v_max)
859 b_east.push_back(i);// east edge
860 else if( (math::abs(curr_point(1) - v_min) < 1e-15) && (curr_point(0) < u_max) )
861 b_south.push_back(i);// south edge
862 else
863 b_north.push_back(i);// north edge
864 }
865 }
866
867 corners.push_back(interiors.size()); // c1
868 corners.push_back(interiors.size() + b_south.size()); // c2
869 corners.push_back(interiors.size() + b_south.size() + b_east.size()); // c3
870 corners.push_back(interiors.size() + b_south.size() + b_east.size() + b_north.size()); // c4
871
872 uv_interiors.resize(2, interiors.size());
873 p_interiors.resize(3, interiors.size());
874 for( size_t i = 0; i < interiors.size(); i++ )
875 {
876 uv_interiors.col(i) = parameters.col(interiors[i]);
877 p_interiors.col(i) = points.col(interiors[i]);
878 }
879
880 uv_west.resize(2, b_west.size());
881 gsMatrix<T> tmp_west(3, b_west.size());
882 for( size_t i = 0; i < b_west.size(); i++ )
883 {
884 uv_west.col(i) = parameters.col(b_west[i]);
885 tmp_west.col(i) = points.col(b_west[i]);
886 }
887
888 uv_east.resize(2, b_east.size());
889 gsMatrix<T> tmp_east(3, b_east.size());
890 for( size_t i = 0; i < b_east.size(); i++ )
891 {
892 uv_east.col(i) = parameters.col(b_east[i]);
893 tmp_east.col(i) = points.col(b_east[i]);
894 }
895
896 uv_south.resize(2, b_south.size());
897 gsMatrix<T> tmp_south(3, b_south.size());
898 for( size_t i = 0; i < b_south.size(); i++ )
899 {
900 uv_south.col(i) = parameters.col(b_south[i]);
901 tmp_south.col(i) = points.col(b_south[i]);
902 }
903
904 uv_north.resize(2, b_north.size());
905 gsMatrix<T> tmp_north(3, b_north.size());
906 for( size_t i = 0; i < b_north.size(); i++ )
907 {
908 uv_north.col(i) = parameters.col(b_north[i]);
909 tmp_north.col(i) = points.col(b_north[i]);
910 }
911
912 uv_south.transposeInPlace();
913 uv_east.transposeInPlace();
914 uv_north.transposeInPlace();
915 uv_west.transposeInPlace();
916
917
918 std::vector<index_t> tmp = uv_south.idxByColumn(0);
919 p_south.resize(tmp_south.rows(), tmp_south.cols());
920 for(size_t i = 0; i<tmp.size(); i++)
921 {
922 p_south.col(i) = tmp_south.col(tmp[i]);
923 }
924 uv_south.transposeInPlace();
925
926
927 tmp.clear();
928 tmp = uv_east.idxByColumn(1);
929 p_east.resize(tmp_east.rows(), tmp_east.cols());
930 for(size_t i = 0; i<tmp.size(); i++)
931 {
932 p_east.col(i) = tmp_east.col(tmp[i]);
933 }
934 uv_east.transposeInPlace();
935
936
937 tmp.clear();
938 tmp = uv_north.idxByColumn(0);
939 std::reverse(tmp.begin(),tmp.end());
940 gsVector<T> tcol = uv_north.col(0).reverse();
941 uv_north.col(0) = tcol;
942 tcol = uv_north.col(1).reverse();
943 uv_north.col(1) = tcol;
944
945 p_north.resize(tmp_north.rows(), tmp_north.cols());
946 for(size_t i = 0; i<tmp.size(); i++)
947 {
948 p_north.col(i) = tmp_north.col(tmp[i]);
949 }
950 uv_north.transposeInPlace();
951
952 tmp.clear();
953 tmp = uv_west.idxByColumn(1);
954 tcol = uv_west.col(0).reverse();
955 uv_west.col(0) = tcol;
956 tcol = uv_west.col(1).reverse();
957 uv_west.col(1) = tcol;
958 std::reverse(tmp.begin(),tmp.end());
959
960 p_west.resize(tmp_west.rows(), tmp_west.cols());
961 for(size_t i = 0; i<tmp.size(); i++)
962 {
963 p_west.col(i) = tmp_west.col(tmp[i]);
964 }
965 uv_west.transposeInPlace();
966
967
968 // reordering of the input point cloud (parameters and points)
969 parameters.resize(uv_interiors.rows(), points.cols());
970 parameters << uv_interiors.row(0), uv_south.row(0), uv_east.row(0), uv_north.row(0), uv_west.row(0),
971 uv_interiors.row(1), uv_south.row(1), uv_east.row(1), uv_north.row(1), uv_west.row(1);
972
973 points.resize(p_interiors.rows(), parameters.cols());
974 points << p_interiors.row(0), p_south.row(0), p_east.row(0), p_north.row(0), p_west.row(0),
975 p_interiors.row(1), p_south.row(1), p_east.row(1), p_north.row(1), p_west.row(1),
976 p_interiors.row(2), p_south.row(2), p_east.row(2), p_north.row(2), p_west.row(2);
977
978 } // end sortPointCloud
979
980
981
989 template<class T>
990 void sampleGridGeometry(const gsMultiPatch<T> & mp,
991 const index_t & numPatch,
992 const index_t & numSamples,
993 gsMatrix<T> & params,
994 gsMatrix<T> & points)
995 {
996 GISMO_ASSERT( numPatch <= mp.nPatches()-1 , "Patch number not found, quitting.");
997
998 const gsGeometry<T> & geometry = mp.patch(numPatch);
999
1000 gsVector<unsigned> numPtsVec(2);
1001 numPtsVec<<numSamples,numSamples;
1002 gsVector<T> a = geometry.support().col(0);
1003 gsVector<T> b = geometry.support().col(1);
1004
1005 params = gsPointGrid(a,b, numPtsVec);
1006 geometry.eval_into(params, points);
1007 }
1008
1009
1018 template<class T>
1019 void sampleScatteredGeometry(const gsMultiPatch<T> & mp,
1020 const index_t & numPatch,
1021 const index_t & numSamples,
1022 index_t & numBdr,
1023 gsMatrix<T> & params,
1024 gsMatrix<T> & points)
1025 {
1026 GISMO_ASSERT( numPatch <= mp.nPatches()-1 , "Patch number not found, quitting.");
1027
1028 const gsGeometry<T> & geometry = mp.patch(numPatch);
1029
1030
1031 // Sample the interior parameters
1032 gsVector<unsigned> numPtsVec(2);
1033 numPtsVec<<numSamples,numSamples;
1034 gsVector<T> a = geometry.support().col(0);
1035 gsVector<T> b = geometry.support().col(1);
1036
1037 T urange= b(0)-a(0); // umax - umin
1038
1039 gsMatrix<T> mu = gsMatrix<T>::Random(1,numSamples); // 3x3 Matrix filled with random numbers between (-1,1)
1040 mu = (mu + gsMatrix<T>::Constant(1,numSamples,1))*urange/2.; // add 1 to the matrix to have values between 0 and 2; multiply with range/2
1041 mu = (mu + gsMatrix<T>::Constant(1,numSamples,a(0))); //set LO as the lower bound (offset)
1042
1043 T vrange= b(1)-a(1); // vmax - vmin
1044 gsMatrix<T> mv= gsMatrix<T>::Random(1,numSamples); // 3x3 Matrix filled with random numbers between (-1,1)
1045 mv = (mv + gsMatrix<T>::Constant(1,numSamples,1))*vrange/2.; // add 1 to the matrix to have values between 0 and 2; multiply with range/2
1046 mv = (mv + gsMatrix<T>::Constant(1,numSamples,a(1))); //set LO as the lower bound (offset)
1047
1048 gsMatrix<T> uv_interiors(2, numSamples);
1049 uv_interiors << mu, mv; //interior parameters
1050
1051 // Sample the boundary and the corner parameters
1052 if (numBdr < 2)
1053 numBdr = math::ceil(math::sqrt(numSamples)) + 2; // number of boundary points
1054
1055 gsMatrix<T> uv_boundary(2, numBdr*4-4);
1056 gsMatrix<T> b_0(1, numBdr-1);
1057 gsMatrix<T> b_1(1, numBdr-1);
1058 gsMatrix<T> b_2(1, numBdr-1);
1059 gsMatrix<T> b_3(1, numBdr-1);
1060
1061 for (index_t pace=0; pace < 4; pace++)
1062 {
1063 if (pace == 0){
1064 gsMatrix<T> mu = gsMatrix<T>::Random(numBdr-2, 1); // 1xnumBdr-1 Matrix filled with random numbers between (-1,1)
1065 mu = (mu + gsMatrix<T>::Constant(numBdr-2, 1, 1))*urange/2.; // add 1 to the matrix to have values between 0 and 2; multiply with range/2
1066 mu = (mu + gsMatrix<T>::Constant(numBdr-2,1,a(0))); //set LO as the lower bound (offset)
1067 mu.sortByColumn(0);
1068 mu = mu.reshape(1, numBdr-2);
1069 b_0 << a(0), mu;
1070 }
1071 if (pace == 1){
1072 gsMatrix<T> mu = gsMatrix<T>::Random(numBdr-2, 1); // 1xnumBdr-1 Matrix filled with random numbers between (-1,1)
1073 mu = (mu + gsMatrix<T>::Constant(numBdr-2, 1, 1))*vrange/2.; // add 1 to the matrix to have values between 0 and 2; multiply with range/2
1074 mu = (mu + gsMatrix<T>::Constant(numBdr-2,1, a(1))); //set LO as the lower bound (offset)
1075 mu.sortByColumn(0);
1076 mu = mu.reshape(1, numBdr-2);
1077 b_1 << a(1), mu;
1078 }
1079 if (pace == 2){
1080 gsMatrix<T> mu = gsMatrix<T>::Random(numBdr-2, 1); // 1xnumBdr-1 Matrix filled with random numbers between (-1,1)
1081 mu = (mu + gsMatrix<T>::Constant(numBdr-2, 1, 1))*urange/2.; // add 1 to the matrix to have values between 0 and 2; multiply with range/2
1082 mu = (mu + gsMatrix<T>::Constant(numBdr-2,1,a(0))); //set LO as the lower bound (offset)
1083 mu.sortByColumn(0);
1084 mu = mu.reshape(1, numBdr-2);
1085 b_2 << b(0), mu.reverse();
1086 }
1087 if (pace == 3){
1088 gsMatrix<T> mu = gsMatrix<T>::Random(numBdr-2, 1); // 1xnumBdr-1 Matrix filled with random numbers between (-1,1)
1089 mu = (mu + gsMatrix<T>::Constant(numBdr-2, 1, 1))*vrange/2.; // add 1 to the matrix to have values between 0 and 2; multiply with range/2
1090 mu = (mu + gsMatrix<T>::Constant(numBdr-2,1,a(1))); //set LO as the lower bound (offset)
1091 mu.sortByColumn(0);
1092 mu = mu.reshape(1, numBdr-2);
1093 b_3 << b(1), mu.reverse();
1094 }
1095 }
1096
1097 gsMatrix<T> u_zeros = gsMatrix<T>::Constant(1, numBdr-1, a(1));
1098 gsMatrix<T> v_zeros = gsMatrix<T>::Constant(1, numBdr-1, a(0));
1099 gsMatrix<T> u_ones = gsMatrix<T>::Constant(1, numBdr-1, b(1));
1100 gsMatrix<T> v_ones = gsMatrix<T>::Constant(1, numBdr-1, b(0));
1101 gsMatrix<T> zeros = gsMatrix<T>::Zero(1, numBdr-1);
1102 gsMatrix<T> ones = gsMatrix<T>::Ones(1, numBdr-1);
1103
1104
1105 uv_boundary << b_0, v_ones, b_2, v_zeros, u_zeros, b_1, u_ones, b_3;
1106
1107 params.resize(2, numSamples + numBdr*4-4);
1108 params << uv_interiors, uv_boundary;
1109
1110 geometry.eval_into(params, points);
1111
1112 }
1113
1114
1115
1116} // namespace gismo
gismo::gsBSplineBasis
A univariate B-spline basis.
Definition gsBSplineBasis.h:700
gismo::gsBSpline
A B-spline function of one argument, with arbitrary target dimension.
Definition gsBSpline.h:51
gismo::gsBasis::derivFunc
gsMatrix< T > derivFunc(const gsMatrix< T > &u, const gsMatrix< T > &coefs) const
Evaluate the derivatives of the function described by coefs at points u.
Definition gsBasis.h:215
gismo::gsBasis::evalFunc
gsMatrix< T > evalFunc(const gsMatrix< T > &u, const gsMatrix< T > &coefs) const
Evaluate the function described by coefs at points u.
Definition gsBasis.h:184
gismo::gsFileData
This class represents an XML data tree which can be read from or written to a (file) stream.
Definition gsFileData.h:34
gismo::gsFileData::dump
void dump(String const &fname="dump") const
Dump file contents to an xml file.
Definition gsFileData.hpp:86
gismo::gsFunctionSet::eval
gsMatrix< T > eval(const gsMatrix< T > &u) const
Evaluate the function,.
Definition gsFunctionSet.hpp:120
gismo::gsFunctionSet::deriv
gsMatrix< T > deriv(const gsMatrix< T > &u) const
Evaluate the derivatives,.
Definition gsFunctionSet.hpp:129
gismo::gsGeometry
Abstract base class representing a geometry map.
Definition gsGeometry.h:93
gismo::gsGeometry::coefs
gsMatrix< T > & coefs()
Definition gsGeometry.h:340
gismo::gsGeometry::uPtr
memory::unique_ptr< gsGeometry > uPtr
Unique pointer for gsGeometry.
Definition gsGeometry.h:100
gismo::gsGeometry::eval_into
void eval_into(const gsMatrix< T > &u, gsMatrix< T > &result) const
Evaluate the function at points u into result.
Definition gsGeometry.hpp:166
gismo::gsGeometry::support
gsMatrix< T > support() const
Returns the range of parameters (same as parameterRange())
Definition gsGeometry.hpp:193
gismo::gsKnotVector
Class for representing a knot vector.
Definition gsKnotVector.h:80
gismo::gsKnotVector::size
size_t size() const
Number of knots (including repetitions).
Definition gsKnotVector.h:242
gismo::gsKnotVector::degree
int degree() const
Returns the degree of the knot vector.
Definition gsKnotVector.hpp:700
gismo::gsMatrix
A matrix with arbitrary coefficient type and fixed or dynamic size.
Definition gsMatrix.h:41
gismo::gsMatrix::reshape
gsAsMatrix< T, Dynamic, Dynamic > reshape(index_t n, index_t m)
Returns the matrix resized to n x m matrix (data is not copied) This function assumes that the matrix...
Definition gsMatrix.h:221
gismo::gsMatrix::khatriRao
gsMatrix khatriRao(const gsEigen::MatrixBase< OtherDerived > &other) const
Returns the Khatri-Rao product of this with other.
Definition gsMatrix.h:537
gismo::gsMatrix::sortByColumn
void sortByColumn(const index_t j)
Sorts rows of matrix by column j.
Definition gsMatrix.h:388
gismo::gsMatrix::idxByColumn
std::vector< index_t > idxByColumn(const index_t j)
Returns the vector permutation of the rows of the matrix by column j.
Definition gsMatrix.h:415
gismo::gsMatrix::kron
gsMatrix kron(const gsEigen::MatrixBase< OtherDerived > &other) const
Returns the Kronecker product of this with other.
Definition gsMatrix.h:524
gismo::gsMultiPatch
Container class for a set of geometry patches and their topology, that is, the interface connections ...
Definition gsMultiPatch.h:100
gismo::gsSparseMatrix
Sparse matrix class, based on gsEigen::SparseMatrix.
Definition gsSparseMatrix.h:139
gismo::gsSparseMatrix::begin
iterator begin(const index_t outer) const
Returns an iterator to the first non-zero elemenent of column \ a outer (or row outer if the matrix i...
Definition gsSparseMatrix.h:256
gismo::gsSparseMatrix::insertTo
void insertTo(_Index i, _Index j, const T val)
Definition gsSparseMatrix.h:289
gismo::gsTensorBSpline
A tensor product of d B-spline functions, with arbitrary target dimension.
Definition gsTensorBSpline.h:45
gismo::gsTensorBSpline::Ptr
memory::shared_ptr< gsTensorBSpline > Ptr
Shared pointer for gsTensorBSpline.
Definition gsTensorBSpline.h:59
gismo::gsVector3d
A fixed-size, statically allocated 3D vector.
Definition gsVector.h:219
gismo::gsVector
A vector with arbitrary coefficient type and fixed or dynamic size.
Definition gsVector.h:37
gismo::gsPointGrid
gsMatrix< T > gsPointGrid(gsVector< T > const &a, gsVector< T > const &b, gsVector< unsigned > const &np)
Construct a Cartesian grid of uniform points in a hypercube, using np[i] points in direction i.
Definition gsPointGrid.hpp:82
short_t
#define short_t
Definition gsConfig.h:35
index_t
#define index_t
Definition gsConfig.h:32
gsWarn
#define gsWarn
Definition gsDebug.h:50
GISMO_UNUSED
#define GISMO_UNUSED(x)
Definition gsDebug.h:112
gsInfo
#define gsInfo
Definition gsDebug.h:43
GISMO_ASSERT
#define GISMO_ASSERT(cond, message)
Definition gsDebug.h:89
gsFileData.h
Utility class which holds I/O XML data to read/write to/from files.
gsGaussRule.h
Provides the Gauss-Legendre quadrature rule.
gsLinearAlgebra.h
This is the main header file that collects wrappers of Eigen for linear algebra.
gsTensorBSpline.h
Represents a tensor-product B-spline patch.
gismo
The G+Smo namespace, containing all definitions for the library.
gismo::sampleScatteredGeometry
void sampleScatteredGeometry(const gsMultiPatch< T > &mp, const index_t &numPatch, const index_t &numSamples, index_t &numBdr, gsMatrix< T > &params, gsMatrix< T > &points)
sampleScatteredGeometry: samples a scattered point cloud from a given geometry
Definition gsModelingUtils.hpp:1019
gismo::removeCol
void removeCol(gsMatrix< T > &mat, int const &removeEnds, int const &nPoints)
remove columes 0, nPoints, 2*nPoints,.. of a given matrix
Definition gsModelingUtils.hpp:359
gismo::scalePoints
void scalePoints(const gsMatrix< T > &xyz, gsMatrix< T > &points)
Function to scale the input points xyz in [0,1]^D and saves it to points.
Definition gsModelingUtils.hpp:688
gismo::flipLR
gsMatrix< T > flipLR(const gsMatrix< T > &mat)
Flip columes from left to right and vice versa.
Definition gsModelingUtils.hpp:285
gismo::crossNorm2Mat
gsMatrix< T > crossNorm2Mat(gsMatrix< T > const &mat1, gsMatrix< T > const &mat2)
Definition gsModelingUtils.hpp:301
gismo::addConstraints
void addConstraints(gsMatrix< T > const &C1, gsMatrix< T > const &d1, gsMatrix< T > const &C2, gsMatrix< T > const &d2, gsMatrix< T > &C, gsMatrix< T > &d)
addConstraints
Definition gsModelingUtils.hpp:319
gismo::scalePts
void scalePts(const std::string &fin, const std::string &fout, index_t uvIdIn, index_t uvIdOut, index_t xyzIdIn, index_t xyzIdOut, T tMin, T tMax, bool verbose)
Scale the points contained in file fin from [0, 1]^D to [tMin, tMax]^D and save it to fout.
Definition gsModelingUtils.hpp:784
gismo::threeOnDiag
void threeOnDiag(const gsSparseMatrix< T > &block, gsSparseMatrix< T > &result)
Definition gsModelingUtils.hpp:664
gismo::gsInterpolateSurface
gsTensorBSpline< 2, T >::Ptr gsInterpolateSurface(const gsMatrix< T > &exactPoints, const gsMatrix< T > &exactValues, const gsMatrix< T > &appxPointsEdges, const gsMatrix< T > &appxValuesEdges, const gsMatrix< T > &appxPointsInt, const gsMatrix< T > &appxValuesInt, const gsMatrix< T > &appxNormalPoints, const gsMatrix< T > &appxNormals, T wEdge, T wInt, T wNormal, T wReg, const gsKnotVector< T > &kv1, const gsKnotVector< T > &kv2, bool force_normal)
Definition gsModelingUtils.hpp:504
gismo::conditionedAngle
T conditionedAngle(gsVector3d< T > vec1, gsVector3d< T > vec2)
Angle between two vector: 0 <= angle <= pi.
Definition gsModelingUtils.hpp:165
gismo::gsInterpolate
gsBSpline< T > gsInterpolate(gsKnotVector< T > &kv, const gsMatrix< T > &preImage, const gsMatrix< T > &image, const gsMatrix< T > &preNormal, const gsMatrix< T > &normal, const gsMatrix< T > &preImageApp, const gsMatrix< T > &imageApp, T const &w_reg, T const &w_app, gsMatrix< T > &outPointResiduals, gsMatrix< T > &outNormalResiduals)
Definition gsModelingUtils.hpp:392
gismo::scaleTo01
T scaleTo01(T tMin, T t, T tMax)
Scale the interval [tMin, tMax] to [0, 1].
Definition gsModelingUtils.hpp:702
gismo::criticalPointOfQuadratic
gsVector< T > criticalPointOfQuadratic(gsMatrix< T > &A, gsMatrix< T > &C, gsVector< T > &d)
Definition gsModelingUtils.hpp:193
gismo::sortPointCloud
void sortPointCloud(gsMatrix< T > &parameters, gsMatrix< T > &points, std::vector< index_t > &corners)
sortPointCloud: sorts the point cloud into interior and boundary points. parameters and points ordere...
Definition gsModelingUtils.hpp:832
gismo::convert2Zero
gsMatrix< T > convert2Zero(gsMatrix< T > const &mat)
convert a with abs(a) < eps=2.220446049250313e-16 into 0
Definition gsModelingUtils.hpp:340
gismo::vectorIntersect
gsVector< T > vectorIntersect(gsVector< T > const &tangent1, gsVector< T > const &tangent2, gsMatrix< T > const &Vert1, gsMatrix< T > const &Vert2)
intersection of two vectors
Definition gsModelingUtils.hpp:135
gismo::scaleFrom01
T scaleFrom01(T tMin, T t, T tMax)
Scale the inteval [0,1] to [tMin, tMax].
Definition gsModelingUtils.hpp:735
gismo::give
S give(S &x)
Definition gsMemory.h:266
gismo::optQuadratic
gsMatrix< T > optQuadratic(gsMatrix< T > const &A, gsMatrix< T > const &b, gsMatrix< T > const &C, gsMatrix< T > const &d)
Find X which solves: min (AX-b)^T (AX-b), s.t. CX=d.
Definition gsModelingUtils.hpp:253
gismo::sampleGridGeometry
void sampleGridGeometry(const gsMultiPatch< T > &mp, const index_t &numPatch, const index_t &numSamples, gsMatrix< T > &params, gsMatrix< T > &points)
sampleGridGeometry: samples a grid point cloud from a given geometry
Definition gsModelingUtils.hpp:990
gismo::scaleGeo
void scaleGeo(const std::string &fin, const std::string &fout, T tMin, T tMax, bool verbose)
Scale the geometry contained in file fin from [0, 1]^D to [tMin, tMax]^D and save it to fout.
Definition gsModelingUtils.hpp:767