G+Smo  24.08.0
Geometry + Simulation Modules
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 
19 #include <gsAssembler/gsGaussRule.h>
20 
21 
22 namespace gismo
23 {
24 
25 
26 template<class T>
27 gsMatrix<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 
61 template<class T>
62 gsMatrix<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 
97 template<class T>
98 gsMatrix<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 
134 template <class T>
135 gsVector<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 
164 template <class T>
165 T 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 
181 template <class T>
182 T 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 
192 template<class T>
193 gsVector<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 
214 template<class T>
215 gsMatrix<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 
243 template<class T>
244 gsMatrix<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 
252 template<class T>
253 gsMatrix<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 
260 template<class T>
261 gsMatrix<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 
272 template<class T>
273 gsMatrix<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 
284 template <class T>
285 gsMatrix<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 
300 template <class T>
301 gsMatrix<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 
318 template <class T>
319 void 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 
339 template <class T>
340 gsMatrix<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 
358 template <class T>
359 void 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 
391 template <class T>
392 gsBSpline<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 
503 template<class T>
504 typename 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 
658 
659 } // namespace gismo
gsGaussRule.h
Provides the Gauss-Legendre quadrature rule.
gismo::gsVector3d
A fixed-size, statically allocated 3D vector.
Definition: gsVector.h:218
gismo::flipLR
gsMatrix< T > flipLR(const gsMatrix< T > &mat)
Flip columes from left to right and vice versa.
Definition: gsModelingUtils.hpp:285
short_t
#define short_t
Definition: gsConfig.h:35
gismo::gsMatrix::kron
gsMatrix kron(const gsEigen::MatrixBase< OtherDerived > &other) const
Returns the Kronecker product of this with other.
Definition: gsMatrix.h:484
gismo::gsTensorBSpline
A tensor product of d B-spline functions, with arbitrary target dimension.
Definition: gsTensorBSpline.h:44
gismo::gsFunctionSet::deriv
gsMatrix< T > deriv(const gsMatrix< T > &u) const
Evaluate the derivatives,.
Definition: gsFunctionSet.hpp:129
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::gsMatrix::khatriRao
gsMatrix khatriRao(const gsEigen::MatrixBase< OtherDerived > &other) const
Returns the Khatri-Rao product of this with other.
Definition: gsMatrix.h:497
gsTensorBSpline.h
Represents a tensor-product B-spline patch.
gismo::give
S give(S &x)
Definition: gsMemory.h:266
gismo::gsFunctionSet::eval
gsMatrix< T > eval(const gsMatrix< T > &u) const
Evaluate the function,.
Definition: gsFunctionSet.hpp:120
gismo::gsKnotVector::size
size_t size() const
Number of knots (including repetitions).
Definition: gsKnotVector.h:242
index_t
#define index_t
Definition: gsConfig.h:32
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_ASSERT
#define GISMO_ASSERT(cond, message)
Definition: gsDebug.h:89
gismo::gsBSpline
A B-spline function of one argument, with arbitrary target dimension.
Definition: gsBSpline.h:50
gismo::gsBSplineBasis
A univariate B-spline basis.
Definition: gsBSplineBasis.h:694
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::gsMatrix::removeCol
void removeCol(index_t i)
Definition: gsMatrix.h:303
gsWarn
#define gsWarn
Definition: gsDebug.h:50
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::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::conditionedAngle
T conditionedAngle(gsVector3d< T > vec1, gsVector3d< T > vec2)
Angle between two vector: 0 &lt;= angle &lt;= pi.
Definition: gsModelingUtils.hpp:165
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::criticalPointOfQuadratic
gsVector< T > criticalPointOfQuadratic(gsMatrix< T > &A, gsMatrix< T > &C, gsVector< T > &d)
Definition: gsModelingUtils.hpp:193
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_UNUSED
#define GISMO_UNUSED(x)
Definition: gsDebug.h:112
gismo::crossNorm2Mat
gsMatrix< T > crossNorm2Mat(gsMatrix< T > const &mat1, gsMatrix< T > const &mat2)
Definition: gsModelingUtils.hpp:301
gismo::gsKnotVector
Class for representing a knot vector.
Definition: gsKnotVector.h:79
gsLinearAlgebra.h
This is the main header file that collects wrappers of Eigen for linear algebra.
gismo::expr::abs
EIGEN_STRONG_INLINE abs_expr< E > abs(const E &u)
Absolute value.
Definition: gsExpressions.h:4486
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::convert2Zero
gsMatrix< T > convert2Zero(gsMatrix< T > const &mat)
convert a with abs(a) &lt; eps=2.220446049250313e-16 into 0
Definition: gsModelingUtils.hpp:340
gismo::gsKnotVector::degree
int degree() const
Returns the degree of the knot vector.
Definition: gsKnotVector.hpp:700
gismo::gsTensorBSpline::Ptr
memory::shared_ptr< gsTensorBSpline > Ptr
Shared pointer for gsTensorBSpline.
Definition: gsTensorBSpline.h:59