27 gsMatrix<T> * innerProduct(
const gsBasis<T>& B1,
const gsBasis<T>& B2)
29 gsMatrix<T> * K =
new gsMatrix<T>(B1.size(), B2.size() ) ;
32 int nGauss = int( ceil(
double(B1.degree(0) + B2.degree(0) + 1)/2 ) );
33 if (nGauss<1) nGauss=1;
35 gsGaussRule<T> QuRule(nGauss);
38 gsMatrix<index_t> act1, act2;
39 gsMatrix<T> ev1 , ev2;
41 typename gsBasis<T>::domainIter domIt = B1.makeDomainIterator();
42 for (; domIt->good(); domIt->next())
45 QuRule.mapTo( domIt->lowerCorner(), domIt->upperCorner(), ngrid, wgrid );
47 B1.eval_into(ngrid,ev1);
48 B2.eval_into(ngrid,ev2);
49 B1.active_into(ngrid,act1);
50 B2.active_into(ngrid,act2);
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) ;
62 gsMatrix<T> * innerProduct1(
const gsBasis<T>& B1,
const gsBasis<T>& B2)
64 gsMatrix<T> * K =
new gsMatrix<T>(B1.size(), B2.size() ) ;
67 int nGauss = int( ceil(
double(B1.degree(0)-1 + B2.degree(0)-1 + 1)/2 ) );
68 if (nGauss<1) nGauss=1;
70 gsGaussRule<T> QuRule(nGauss);
73 gsMatrix<index_t> act1, act2;
74 gsMatrix<T> ev1 , ev2;
76 typename gsBasis<T>::domainIter domIt = B1.makeDomainIterator();
77 for (; domIt->good(); domIt->next())
80 QuRule.mapTo( domIt->lowerCorner(), domIt->upperCorner(), ngrid, wgrid );
82 B1.deriv_into(ngrid,ev1);
83 B2.deriv_into(ngrid,ev2);
84 B1.active_into(ngrid,act1);
85 B2.active_into(ngrid,act2);
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) ;
98 gsMatrix<T> * innerProduct2(
const gsBasis<T>& B1,
const gsBasis<T>& B2)
100 gsMatrix<T> * K =
new gsMatrix<T>(B1.size(), B2.size() ) ;
103 int nGauss = int( ceil(
double(B1.degree(0)-2 + B2.degree(0)-2 + 1)/2 ) );
104 if (nGauss<1) nGauss=1;
106 gsGaussRule<T> QuRule(nGauss);
109 gsMatrix<index_t> act1, act2;
110 gsMatrix<T> ev1 , ev2;
112 typename gsBasis<T>::domainIter domIt = B1.makeDomainIterator();
113 for (; domIt->good(); domIt->next())
116 QuRule.mapTo( domIt->lowerCorner(), domIt->upperCorner(), ngrid, wgrid );
118 B1.deriv2_into(ngrid,ev1);
119 B2.deriv2_into(ngrid,ev2);
120 B1.active_into(ngrid,act1);
121 B2.active_into(ngrid,act2);
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) ;
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) );
149 T detMatrixab = abc1(0)*abc2(1)-abc1(1)*abc2(0);
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) );
157 unknown(0) = .5*Vert1(0) + .5*Vert2(0);
158 unknown(1) = .5*Vert1(1) + .5*Vert2(1);
166 { T dotp = vec1.dot(vec2)/( vec1.norm()*vec2.norm() );
169 gsWarn<<
"gsModelingUtils: truncations done for std::acos \n";
171 if ( dotp>1 ) {dotp = 1;
172 gsWarn<<
"gsModelingUtils: truncations done for std::acos \n";
175 T ag = math::acos(dotp);
184 T ag = conditionedAngle<T>(vec1,vec2);
185 T cag = ( normal.dot( vec1.cross( vec2 ) ) >= 0 ) ? ag : (T)(2*EIGEN_PI)-ag;
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();
205 longD.head(n).setZero();
208 gsMatrix<T> result = combined.fullPivLu().solve(longD);
209 return result.block(0, 0, n, 1);
223 assert(d.cols()==1 && C.rows()==d.rows());
225 if (b.rows()!=1) bt=b.transpose();
else bt=b;
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();
234 longD.head(n) = -bt.transpose();
237 gsMatrix<T> result = combined.fullPivLu().solve(longD);
238 return result.block(0, 0, n, 1);
248 return criticalPointOfQuadratic<T>(A,b,C,d);
256 return criticalPointOfQuadratic<T>( (A.transpose())*A, (-2)*( A.transpose() )*b, C, d);
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);
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);
287 return mat.rowwise().reverse();
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();
309 for (
size_t i=0; i<= nc-1; i++)
312 rcross.col(i) = tem / tem.norm();
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;
334 d.block(0,0,nr1,1) = d1;
335 d.block(nr1,0,nr2,1) = d2;
342 T eps=2.220446049250313e-16;
347 for (
int i=0;i!=n1;i++)
349 for (
int j=0;j!=n2;j++)
351 if (
math::abs(mat(i,j))< eps) matc(i,j)=0.;
361 assert(removeEnds==1 || removeEnds==2);
362 int nPeriod = mat.cols()/nPoints;
363 assert( nPeriod*nPoints == mat.cols() );
367 for (
int i=nPeriod-1;i>=0;i--)
369 ind2 = i*nPoints + nPoints-1;
375 for (
int i=nPeriod-1;i>=0;i--)
377 ind2 = i*nPoints + nPoints-1;
396 T
const & w_reg,T
const & w_app,
410 int nip = image.cols();
411 int nn=normal.cols();
415 Nu = bs.
eval(preImage.row(0));
416 dNu = bs.
deriv(preImage.row(0));
417 Nu.transposeInPlace();
418 dNu.transposeInPlace();
419 dNu_nm = bs.
deriv(preNormal.row(0));
420 dNu_nm.transposeInPlace();
421 gsMatrix<T> AdN = normal.row(0).asDiagonal() * dNu_nm;
422 gsMatrix<T> BdN = normal.row(1).asDiagonal() * dNu_nm;
425 NuApp = bs.
eval(preImageApp.row(0));
426 NuApp.transposeInPlace();
431 int nss = dimI*ntcp + dimI*nip + nn;
434 Ass.setZero(); bss.setZero();
435 gsMatrix<T> Hess = w_reg*(T)(2)*(*Q) + w_app*(T)(2)*(NuApp.transpose())* NuApp;
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();
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();
447 Ass.block(2*ntcp,0,nip,ntcp) = Nu;
448 bss.block(2*ntcp,0,nip,1) = (image.row(0)).transpose();
450 Ass.block(2*ntcp+nip,ntcp,nip,ntcp) = Nu;
451 bss.block(2*ntcp+nip,0,nip,1) = (image.row(1)).transpose();
453 Ass.block(2*ntcp+2*nip,0,nn,ntcp) = AdN;
454 Ass.block(2*ntcp+2*nip,ntcp,nn,ntcp) = BdN;
457 tcp.col(0) = result.block(0 , 0, ntcp, 1);
458 tcp.col(1) = result.block(ntcp, 0, ntcp, 1);
480 outPointResiduals = (NuApp * tcp).transpose() - imageApp;
481 outNormalResiduals = AdN * tcp.col(0) + BdN * tcp.col(1);
509 T wEdge, T wInt, T wNormal, T wReg,
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");
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;
538 ident1.setIdentity();
539 ident2.setIdentity();
540 Nu = bs1.
evalFunc(exactPoints.row(0), ident1);
541 Nv = bs2.
evalFunc(exactPoints.row(1), ident2);
542 npts = exactPoints.cols();
545 Ecor.setZero();ecor.setZero();
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();
555 Nu = bs1.
evalFunc(appxPointsEdges.row(0), ident1);
556 Nv = bs2.
evalFunc(appxPointsEdges.row(1), ident2);
557 npts = appxPointsEdges.cols();
560 AappEdge.setZero();bappEdge.setZero();
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();
570 Nu = bs1.
evalFunc(appxPointsInt.row(0), ident1);
571 Nv = bs2.
evalFunc(appxPointsInt.row(1), ident2);
572 npts = appxPointsInt.cols();
575 AappInt.setZero();bappInt.setZero();
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();
585 Nu = bs1.
evalFunc(appxNormalPoints.row(0), ident1);
586 dNu = bs1.
derivFunc(appxNormalPoints.row(0), ident1);
587 Nv = bs2.
evalFunc(appxNormalPoints.row(1), ident2);
588 dNv = bs2.
derivFunc(appxNormalPoints.row(1), ident2);
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();}
598 dRdu.transposeInPlace();
599 dRdv.transposeInPlace();
600 int nnor = Nx.rows();
604 Anor.block(0,0,nnor,n1*n2) = Nx.asDiagonal()*dRdu;
605 Anor.block(0,n1*n2,nnor,n1*n2) = Ny.asDiagonal()*dRdu;
606 Anor.block(0,2*n1*n2,nnor,n1*n2) = Nz.asDiagonal()*dRdu;
607 Anor.block(nnor,0,nnor,n1*n2) = Nx.asDiagonal()*dRdv;
608 Anor.block(nnor,n1*n2,nnor,n1*n2) = Ny.asDiagonal()*dRdv;
609 Anor.block(nnor,2*n1*n2,nnor,n1*n2) = Nz.asDiagonal()*dRdv;
616 if (kv1==kv2) { N = M; N1 = M1; N2 = M2; }
else
618 N = innerProduct(bs2, bs2);
619 N1 = innerProduct1(bs2, bs2);
620 N2 = innerProduct2(bs2, bs2);
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;
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;
639 cp.resize( cp.rows() / 3, 3);
643 delete M;
delete M1;
delete M2;
644 if (kv1!=kv2) {
delete N;
delete N1;
delete N2;}
647 for(
index_t idxConstr = 0; idxConstr < exactPoints.cols(); idxConstr++)
650 gsMatrix<T> expectVal = exactValues.col(idxConstr);
652 GISMO_ASSERT((expectVal - surfVal).norm() < 0.0001,
"Fit surface did not satisfy exact constraints");
Provides the Gauss-Legendre quadrature rule.
A fixed-size, statically allocated 3D vector.
Definition: gsVector.h:218
gsMatrix< T > flipLR(const gsMatrix< T > &mat)
Flip columes from left to right and vice versa.
Definition: gsModelingUtils.hpp:285
#define short_t
Definition: gsConfig.h:35
gsMatrix kron(const gsEigen::MatrixBase< OtherDerived > &other) const
Returns the Kronecker product of this with other.
Definition: gsMatrix.h:484
A tensor product of d B-spline functions, with arbitrary target dimension.
Definition: gsTensorBSpline.h:44
gsMatrix< T > deriv(const gsMatrix< T > &u) const
Evaluate the derivatives,.
Definition: gsFunctionSet.hpp:129
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
gsMatrix khatriRao(const gsEigen::MatrixBase< OtherDerived > &other) const
Returns the Khatri-Rao product of this with other.
Definition: gsMatrix.h:497
Represents a tensor-product B-spline patch.
S give(S &x)
Definition: gsMemory.h:266
gsMatrix< T > eval(const gsMatrix< T > &u) const
Evaluate the function,.
Definition: gsFunctionSet.hpp:120
size_t size() const
Number of knots (including repetitions).
Definition: gsKnotVector.h:242
#define index_t
Definition: gsConfig.h:32
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
#define GISMO_ASSERT(cond, message)
Definition: gsDebug.h:89
A B-spline function of one argument, with arbitrary target dimension.
Definition: gsBSpline.h:50
A univariate B-spline basis.
Definition: gsBSplineBasis.h:694
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
void removeCol(index_t i)
Definition: gsMatrix.h:303
#define gsWarn
Definition: gsDebug.h:50
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
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
T conditionedAngle(gsVector3d< T > vec1, gsVector3d< T > vec2)
Angle between two vector: 0 <= angle <= pi.
Definition: gsModelingUtils.hpp:165
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
gsVector< T > criticalPointOfQuadratic(gsMatrix< T > &A, gsMatrix< T > &C, gsVector< T > &d)
Definition: gsModelingUtils.hpp:193
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
#define GISMO_UNUSED(x)
Definition: gsDebug.h:112
gsMatrix< T > crossNorm2Mat(gsMatrix< T > const &mat1, gsMatrix< T > const &mat2)
Definition: gsModelingUtils.hpp:301
Class for representing a knot vector.
Definition: gsKnotVector.h:79
This is the main header file that collects wrappers of Eigen for linear algebra.
EIGEN_STRONG_INLINE abs_expr< E > abs(const E &u)
Absolute value.
Definition: gsExpressions.h:4488
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
gsMatrix< T > convert2Zero(gsMatrix< T > const &mat)
convert a with abs(a) < eps=2.220446049250313e-16 into 0
Definition: gsModelingUtils.hpp:340
int degree() const
Returns the degree of the knot vector.
Definition: gsKnotVector.hpp:700
memory::shared_ptr< gsTensorBSpline > Ptr
Shared pointer for gsTensorBSpline.
Definition: gsTensorBSpline.h:59