diff --git a/.gitignore b/.gitignore index 9874cac..21cc2e5 100644 --- a/.gitignore +++ b/.gitignore @@ -3,4 +3,5 @@ ToDo-Liste.txt kernel/matrixfree/.svn *.o temp -*. +*.m~ +*.m~ diff --git a/FAIRstartup.m b/FAIRstartup.m index e1a9041..7d22e5f 100644 --- a/FAIRstartup.m +++ b/FAIRstartup.m @@ -42,6 +42,7 @@ addpath(genpath(fullfile(FAIRpath,'kernel'))); addpath(fullfile(FAIRpath,'temp')); addpath(genpath(fullfile(FAIRpath,'add-ons'))); + message({ '[new to FAIR? try the numerous examples in kernel/examples]' '% ==============================================================================' diff --git a/add-ons/FAIRFEM/FEMPIREobjFctn.m b/add-ons/FAIRFEM/FEMPIREobjFctn.m index 92614c8..83f8ad4 100755 --- a/add-ons/FAIRFEM/FEMPIREobjFctn.m +++ b/add-ons/FAIRFEM/FEMPIREobjFctn.m @@ -75,6 +75,7 @@ yc = reshape(yc,[],dim); yRef = reshape(yRef,[],dim); doDerivative = (nargout>2); % flag for necessity of derivatives +matrixFree = regularizer('get','matrixFree'); % do the work ------------------------------------------------------------ % define barycentric interpolation @@ -83,7 +84,7 @@ vol = Mesh.vol; % compute interpolated image and derivative -[Tc,dT] = imgModel(T,omega,Mesh.mfPi(yc,'C'),'doDerivative',doDerivative); +[Tc,dT] = imgModel(T,omega,Mesh.mfPi(yc,'C'),'doDerivative',doDerivative,'matrixFree',matrixFree); % apply intensity modulation [det,dDet] = getDeterminant(Mesh,yc,'doDerivative',doDerivative); Tcmod = Tc .* det; @@ -121,29 +122,137 @@ H = dr'*d2psi*dr + d2S; else % derivatives rather explicit - dr = dres*dT; - dD = dD*dT; - dJ = reshape(Mesh.mfPi(reshape(dD',[],Mesh.dim),'C'),1,[]) + dS; + %P = @(x) reshape(Mesh.mfPi(reshape(x,[],dim),'C'),[],1); + %dTcmod = P((sdiag(det)*dT)')' + sdiag(Tc)*dDet; + p = (dD(:).*det).*dT; + dD = reshape(Mesh.mfPi(p,'C'),1,[]) + (dDet.dDetadj(dD'.*Tc))'; + + dJ = dD + dS; % approximation to d2D in matrix free mode % d2D = P'*dr'*d2psi*dr*P % P and P' are operators matrix free - H.Mesh = Mesh; + H.Mesh = Mesh; + H.omega = omega; + H.m = m; H.d2D.how = 'P''*dr''*d2psi*dr*P'; - H.d2D.P = @(x) Mesh.mfPi(reshape(x(:),[],Mesh.dim),'C'); - H.d2D.dr = dr; + H.d2D.P = @(x) reshape(Mesh.mfPi(x,'C'),[],1); + H.d2D.Tc = Tc; + H.d2D.dT = dT; + H.d2D.Jac = det; + H.d2D.dJac = dDet.dDet; + H.d2D.dJacadj = dDet.dDetadj; + H.d2D.diag = @(Mesh,yc) getDiag(Mesh,det,Tc,dT,yc); + H.d2D.dres = dres; H.d2D.d2psi = d2psi; + H.solver = d2S.solver; H.d2S = d2S; end; - % shortcut for sparse diagonal matrix function A = sdiag(v) A = spdiags(v(:),0,numel(v),numel(v)); +% diagonal of hessian matrix +% sum(sdiag(vol)dTcmod.^2,1) +% where dTcmod = sdiag(det)*dT*P + sdiag(Tc)*dDet; +% implemention: sum(sdiag(vol)*(sdiag(det)*dT*P(.^2 + +% sdiag(vol)*(sdiag(Tc)*dDet).^2,1) +function D = getDiag(Mesh,det,Tc,dT,yc) + +vol = Mesh.vol; +dim = Mesh.dim; + +if dim == 2 + + % first term + D1 = Mesh.mfPi(vol.*(det.*dT).^2,'C'); + D1 = D1(:)./3; + + % second term + dphi = Mesh.dphi; + dphi1 = dphi{1}; dphi2 = dphi{2}; dphi3 = dphi{3}; + + % compute diagonal of volume regularizer + % MB : dDet = [sdiag(By(:,4)), sdiag(-By(:,3)) , sdiag(-By(:,2)), sdiag(By(:,1))] * B; + [dx1, dx2] = getGradientMatrixFEM(Mesh,1); + yc = reshape(yc,[],2); + By = [dx1.D(yc(:,1)) dx2.D(yc(:,1)) dx1.D(yc(:,2)) dx2.D(yc(:,2))]; + + % get boundaries + Dxi = @(i,j) Mesh.mfPi(vol.*(Tc.*By(:,j).*dphi1(:,i)).^2,1) ... + + Mesh.mfPi(vol.*(Tc.*By(:,j).*dphi2(:,i)).^2,2) ... + + Mesh.mfPi(vol.*(Tc.*By(:,j).*dphi3(:,i)).^2,3); % diagonal of Dxi'*Dxi + + Dxy = @(i,j,k,l) Mesh.mfPi(vol.*(Tc.*By(:,j).*dphi1(:,i)).*(Tc.*By(:,l).*dphi1(:,k)),1) ... % byproduct terms for verifications + + Mesh.mfPi(vol.*(Tc.*By(:,j).*dphi2(:,i)).*(Tc.*By(:,l).*dphi2(:,k)),2) ... + + Mesh.mfPi(vol.*(Tc.*By(:,j).*dphi3(:,i)).*(Tc.*By(:,l).*dphi3(:,k)),3); % diagonal of Dxi'*Dxi + + D2 = [Dxi(1,4)+ Dxi(2,3) - 2*Dxy(1,4,2,3) ;Dxi(1,2)+ Dxi(2,1) - 2*Dxy(1,2,2,1)]; + + % third byproduct term + D3 = [(Mesh.mfPi(vol.*((1/3)*det.*dT(:,1)).*Tc.*(By(:,4).*dphi1(:,1) - By(:,3).*dphi1(:,2)),1)) + ... + (Mesh.mfPi(vol.*((1/3)*det.*dT(:,1)).*Tc.*(By(:,4).*dphi2(:,1) - By(:,3).*dphi2(:,2)),2)) + ... + (Mesh.mfPi(vol.*((1/3)*det.*dT(:,1)).*Tc.*(By(:,4).*dphi3(:,1) - By(:,3).*dphi3(:,2)),3)) ;... + + (Mesh.mfPi(vol.*((1/3)*det.*dT(:,2)).*Tc.*(By(:,1).*dphi1(:,2) - By(:,2).*dphi1(:,1)),1)) + ... + (Mesh.mfPi(vol.*((1/3)*det.*dT(:,2)).*Tc.*(By(:,1).*dphi2(:,2) - By(:,2).*dphi2(:,1)),2)) + ... + (Mesh.mfPi(vol.*((1/3)*det.*dT(:,2)).*Tc.*(By(:,1).*dphi3(:,2) - By(:,2).*dphi3(:,1)),3))]; + D = D1 + D2 + 2*D3; + +else + + % first term + D1 = Mesh.mfPi(vol.*(det.*dT).^2,'C'); + D1 = D1(:)./4; + + dphi = Mesh.dphi; + dphi1 = dphi{1}; dphi2 = dphi{2}; dphi3 = dphi{3}; dphi4 = dphi{4}; + + [cof, dCof] = cofactor3D([],Mesh,0,1); + yc = reshape(yc,[],3); + dx1 = Mesh.mfdx1; dx2 = Mesh.mfdx2; dx3 = Mesh.mfdx3; + %det = dx1.D(yc(:,1)).*cof{1}(yc) + dx2.D(yc(:,1)).*cof{2}(yc) + dx3.D(yc(:,1)).*cof{3}(yc); + + % get boundaries + Dxi = @(i,j) Mesh.mfPi(vol.*(Tc.*cof{j}(yc).*dphi1(:,i)).^2,1) ... + + Mesh.mfPi(vol.*(Tc.*cof{j}(yc).*dphi2(:,i)).^2,2) ... + + Mesh.mfPi(vol.*(Tc.*cof{j}(yc).*dphi3(:,i)).^2,3) ... + + Mesh.mfPi(vol.*(Tc.*cof{j}(yc).*dphi4(:,i)).^2,4); % diagonal of Dxi'*Dxi + + Dxy = @(i,j,k,l) Mesh.mfPi(vol.*(Tc.*cof{j}(yc).*dphi1(:,i)).*(Tc.*cof{l}(yc).*dphi1(:,k)),1) ... % byproduct terms for verifications + + Mesh.mfPi(vol.*(Tc.*cof{j}(yc).*dphi2(:,i)).*(Tc.*cof{l}(yc).*dphi2(:,k)),2) ... + + Mesh.mfPi(vol.*(Tc.*cof{j}(yc).*dphi3(:,i)).*(Tc.*cof{l}(yc).*dphi3(:,k)),3) ... + + Mesh.mfPi(vol.*(Tc.*cof{j}(yc).*dphi4(:,i)).*(Tc.*cof{l}(yc).*dphi4(:,k)),4); % diagonal of Dxi'*Dxi + + D21 = Dxi(1,1)+ Dxi(2,2) + Dxi(3,3) + 2*Dxy(1,1,2,2) + 2*Dxy(1,1,3,3) + 2*Dxy(2,2,3,3); + D22 = Dxi(1,4)+ Dxi(2,5) + Dxi(3,6) + 2*Dxy(1,4,2,5) + 2*Dxy(1,4,3,6) + 2*Dxy(2,5,3,6); + D23 = Dxi(1,7)+ Dxi(2,8) + Dxi(3,9) + 2*Dxy(1,7,2,8) + 2*Dxy(1,7,3,9) + 2*Dxy(2,8,3,9); + + D2 = [D21;D22;D23]; + + % third term + D3 = [(Mesh.mfPi(vol.*((1/4)*det.*dT(:,1)).*Tc.*(cof{1}(yc).*dphi1(:,1) + cof{2}(yc).*dphi1(:,2) + cof{3}(yc).*dphi1(:,3)),1)) + ... + (Mesh.mfPi(vol.*((1/4)*det.*dT(:,1)).*Tc.*(cof{1}(yc).*dphi2(:,1) + cof{2}(yc).*dphi2(:,2) + cof{3}(yc).*dphi2(:,3)),2)) + ... + (Mesh.mfPi(vol.*((1/4)*det.*dT(:,1)).*Tc.*(cof{1}(yc).*dphi3(:,1) + cof{2}(yc).*dphi3(:,2) + cof{3}(yc).*dphi3(:,3)),3)) + ... + (Mesh.mfPi(vol.*((1/4)*det.*dT(:,1)).*Tc.*(cof{1}(yc).*dphi4(:,1) + cof{2}(yc).*dphi4(:,2) + cof{3}(yc).*dphi4(:,3)),4));... + + (Mesh.mfPi(vol.*((1/4)*det.*dT(:,2)).*Tc.*(cof{4}(yc).*dphi1(:,1) + cof{5}(yc).*dphi1(:,2) + cof{6}(yc).*dphi1(:,3)),1)) + ... + (Mesh.mfPi(vol.*((1/4)*det.*dT(:,2)).*Tc.*(cof{4}(yc).*dphi2(:,1) + cof{5}(yc).*dphi2(:,2) + cof{6}(yc).*dphi2(:,3)),2)) + ... + (Mesh.mfPi(vol.*((1/4)*det.*dT(:,2)).*Tc.*(cof{4}(yc).*dphi3(:,1) + cof{5}(yc).*dphi3(:,2) + cof{6}(yc).*dphi3(:,3)),3)) + ... + (Mesh.mfPi(vol.*((1/4)*det.*dT(:,2)).*Tc.*(cof{4}(yc).*dphi4(:,1) + cof{5}(yc).*dphi4(:,2) + cof{6}(yc).*dphi4(:,3)),4));... + + (Mesh.mfPi(vol.*((1/4)*det.*dT(:,3)).*Tc.*(cof{7}(yc).*dphi1(:,1) + cof{8}(yc).*dphi1(:,2) + cof{9}(yc).*dphi1(:,3)),1)) + ... + (Mesh.mfPi(vol.*((1/4)*det.*dT(:,3)).*Tc.*(cof{7}(yc).*dphi2(:,1) + cof{8}(yc).*dphi2(:,2) + cof{9}(yc).*dphi2(:,3)),2)) + ... + (Mesh.mfPi(vol.*((1/4)*det.*dT(:,3)).*Tc.*(cof{7}(yc).*dphi3(:,1) + cof{8}(yc).*dphi3(:,2) + cof{9}(yc).*dphi3(:,3)),3)) + ... + (Mesh.mfPi(vol.*((1/4)*det.*dT(:,3)).*Tc.*(cof{7}(yc).*dphi4(:,1) + cof{8}(yc).*dphi4(:,2) + cof{9}(yc).*dphi4(:,3)),4))]; + + D = D1 + D2 + 2*D3; +end + function runMinimalExample setup2DGaussianData diff --git a/add-ons/FAIRFEM/FEMPIREsolveGN_PCG.m b/add-ons/FAIRFEM/FEMPIREsolveGN_PCG.m new file mode 100644 index 0000000..c72ba79 --- /dev/null +++ b/add-ons/FAIRFEM/FEMPIREsolveGN_PCG.m @@ -0,0 +1,61 @@ +%============================================================================== +% This code is part of the VAMPIRE app for the Matlab-based toolbox +% FAIR - Flexible Algorithms for Image Registration. +% For details see +% - https://github.com/C4IR/VAMPIRE.m +%============================================================================== +% PCG solver for FEMPIRE +% + +function dy = FEMPIREsolveGN_PCG(rhs, H, maxIterCG, tolCG) + +h = (H.omega(2:2:end)-H.omega(1:2:end))./H.m; +hd = prod(h); +% set pointers to simplify notation +Mesh = H.Mesh; +P = H.d2D.P; +Jac = H.d2D.Jac; +dJac = H.d2D.dJac; +dJacadj = H.d2D.dJacadj; +Tc = H.d2D.Tc; +dT = H.d2D.dT; +yc = H.d2S.yc; +dres = H.d2D.dres; +dim = Mesh.dim; + +% Tcmod(y) = T(P*y) * Jac(y) +% ==> dTcmod(y) = Jac(y) * (dT(y)*P) + T(P*y) * dJac(y) +% ==> dTcmod'(y) = (P'* dT') Jac(y) + (dJac(y))' * Tc +dTcmod = @(x) vecmatprod1(Jac,dT,H.d2D.P(reshape(x,[],dim))) + Tc .* dJac(reshape(x,[],dim)); +dTcmodadj = @(x) H.d2D.P(vecmatprod1adj(dT,Jac,x)) + dJacadj(Tc.*x); + +M = @(x) dTcmodadj(dres'*H.d2D.d2psi*(dres * dTcmod(x))); +%Hoperator = @(x) M(x) + H.d2S.d2S(x,H.omega,H.m,H.d2S.yc); +Hoperator = @(x) M(x) + H.d2S.d2S(x); + +D = H.d2D.diag(Mesh,yc) + H.d2S.diag(H.d2S.yc); +Preconditioner = @(x) D.\x; % Jacobi preconditioner +[dy,flag,relres,iter] = pcg(Hoperator,rhs,tolCG,maxIterCG,Preconditioner); + + +function p= vecmatprod1(Jac,dI,x) +% We aim to compute diag(Jac) * [dI] * x +dim = size(dI,2); +n = length(Jac); + +x = reshape(x,[],dim); +p = zeros(n,1); +for d=1:dim + p = p + dI(:,d) .* x(:,d); +end +p = p.*Jac; + +function p = vecmatprod1adj(dI,Jac,x) +% We aim to compute dI' * [Jac] * x +dim = size(dI,2); +n = length(Jac); + +p = zeros(n,dim); +for d=1:dim + p(:,d) = dI(:,d) .* Jac.*x; +end diff --git a/add-ons/FAIRFEM/FEMSSDsolveGN_PCG.m b/add-ons/FAIRFEM/FEMSSDsolveGN_PCG.m new file mode 100644 index 0000000..7720d5f --- /dev/null +++ b/add-ons/FAIRFEM/FEMSSDsolveGN_PCG.m @@ -0,0 +1,44 @@ +%============================================================================== +% This code is part of the VAMPIRE app for the Matlab-based toolbox +% FAIR - Flexible Algorithms for Image Registration. +% For details see +% - https://github.com/C4IR/VAMPIRE.m +%============================================================================== +% PCG solver for FEMSSD +% + +function dy = FEMSSDsolveGN_PCG(rhs, H, maxIterCG, tolCG) + +% set pointers to simplify notation +Mesh = H.Mesh; +P = H.d2D.P; +Tc = H.d2D.Tc; +dT = H.d2D.dT; +yc = H.d2S.yc; +dres = H.d2D.dres; +dim = Mesh.dim; + +% Tcmod(y) = T(P*y) +% ==> dTcmod(y) = (dT(y)*P) +% ==> dTcmod'(y) = (P'* dT')(y) +%dTcmod = @(x) vecmatprod1(dT,H.d2D.P(reshape(x,[],dim))); +%dTcmodadj = @(x) H.d2D.P(vecmatprod1adj(dT,x)); + +dTcmod = @(x) vecmatprod1(dT,H.d2D.P(reshape(x,[],dim))); +dTcmodadj = @(x) H.d2D.P(dT.*x); + +M = @(x) dTcmodadj(dres'*H.d2D.d2psi*(dres * dTcmod(x))); +Hoperator = @(x) M(x) + H.d2S.d2S(x); + +D = H.d2D.diag() + H.d2S.diag(H.d2S.yc); +Preconditioner = @(x) D.\x; % Jacobi preconditioner +[dy,flag,relres,iter] = pcg(Hoperator,rhs,tolCG,maxIterCG,Preconditioner); + + +function p= vecmatprod1(dI,x) +% We aim to compute [dI] * x +dim = size(dI,2); + +x = reshape(x,[],dim); +p = sum(dI.*x,2); + diff --git a/add-ons/FAIRFEM/FEMobjFctn.m b/add-ons/FAIRFEM/FEMobjFctn.m index ec1b23f..feaa822 100755 --- a/add-ons/FAIRFEM/FEMobjFctn.m +++ b/add-ons/FAIRFEM/FEMobjFctn.m @@ -68,6 +68,7 @@ yc = reshape(yc,[],dim); yRef = reshape(yRef,[],dim); doDerivative = (nargout>2); % flag for necessity of derivatives +matrixFree = regularizer('get','matrixFree'); % do the work ------------------------------------------------------------ @@ -75,7 +76,7 @@ vol = Mesh.vol; % compute interpolated image and derivative -[Tc,dT] = imgModel(T,omega,Mesh.mfPi(yc,'C'),'doDerivative',doDerivative); +[Tc,dT] = imgModel(T,omega,Mesh.mfPi(yc,'C'),'doDerivative',doDerivative,'matrixFree',matrixFree); % compute SSD distance rc = Tc-Rc; % the residual @@ -109,22 +110,28 @@ H = dr'*d2psi*dr + d2S; else % derivatives rather explicit - P = @(x) reshape(Mesh.mfPi(reshape(x,[],dim),'C'),[],1); - dr = dres*dT; - dD = dD*dT; - dJ = P(dD')' + dS; - - % approximation to d2D in matrix free mode - % d2D = P'*dr'*d2psi*dr*P - % P and P' are operators matrix free - H.Mesh = Mesh; - H.d2D.how = 'P''*dr''*d2psi*dr*P'; - H.d2D.P = P; - H.d2D.dr = dr; - H.d2D.d2psi = d2psi; - H.solver = d2S.solver; - - H.d2S = d2S; + p = dD(:).*dT; + dD = reshape(Mesh.mfPi(p,'C'),1,[]); + + dJ = dD + dS; + + % approximation to d2D in matrix free mode + % d2D = P'*dr'*d2psi*dr*P + % P and P' are operators matrix free + H.Mesh = Mesh; + H.omega = omega; + H.m = m; + H.d2D.how = 'P''*dr''*d2psi*dr*P'; + H.d2D.P = @(x) reshape(Mesh.mfPi(x,'C'),[],1); + H.d2D.Tc = Tc; + H.d2D.dT = dT; + H.d2D.diag = @() getDiag(Mesh,dT); + H.d2D.dres = dres; + H.d2D.d2psi = d2psi; + H.solver = d2S.solver; + + H.d2S = d2S; + end; @@ -134,6 +141,28 @@ A = spdiags(v(:),0,numel(v),numel(v)); +% diagonal of hessian matrix +% sum(sdiag(vol)dTcmod.^2,1) +% where dTcmod = dT*P; +function D = getDiag(Mesh,dT) + +vol = Mesh.vol; +dim = Mesh.dim; + +if dim == 2 + + % first term + D = Mesh.mfPi(vol.*(dT).^2,'C'); + D = D(:)./3; + +else + + % first term + D = Mesh.mfPi(vol.*(dT).^2,'C'); + D = D(:)./4; + +end + function runMinimalExample setup2DhandData diff --git a/add-ons/FAIRFEM/MLIRFEM.m b/add-ons/FAIRFEM/MLIRFEM.m index a8b8fff..eb0cc19 100755 --- a/add-ons/FAIRFEM/MLIRFEM.m +++ b/add-ons/FAIRFEM/MLIRFEM.m @@ -179,7 +179,7 @@ yn = trafo(wc,xc(:)); yRef = reshape(yn,[],dim); end; - Mesh.xn = reshape(yn,[],dim); + Mesh.xn = yRef;%reshape(yn,[],dim); % compute starting guess y0 if level == minLevel, diff --git a/add-ons/FAIRFEM/cofactor3D.m b/add-ons/FAIRFEM/cofactor3D.m new file mode 100644 index 0000000..e44f7cb --- /dev/null +++ b/add-ons/FAIRFEM/cofactor3D.m @@ -0,0 +1,107 @@ +% For details see +% - https://github.com/C4IR/FAIRFEM +%============================================================================== +% +% function [cof,dCof] = cofactor3D(By,Mesh,doDerivative,matrixFree) +% +% cofactor matrix, and derivative +% +% Input: +% +% By - partial derivative of y wrt x +% Mesh - description of mesh, struct +% +% Output: +% +% [cof, dcof] - cofactor and partial derivative operators +% +%============================================================================== + +function [cof,dCof] = cofactor3D(By,Mesh,doDerivative,matrixFree) + +if ~matrixFree + + dCof = []; + D = @(i,j) By(:,(j-1)*3+i); % gets partial_i y^j + cof = [ + D(2,2).*D(3,3)-D(3,2).*D(2,3),...%ok + -D(1,2).*D(3,3)+D(3,2).*D(1,3),...%trouble + D(1,2).*D(2,3)-D(2,2).*D(1,3),...%ok + -D(2,1).*D(3,3)+D(3,1).*D(2,3),...%ok + D(1,1).*D(3,3)-D(3,1).*D(1,3),...%ok + -D(1,1).*D(2,3)+D(2,1).*D(1,3),...%ok + D(2,1).*D(3,2)-D(3,1).*D(2,2),...%ok + -D(1,1).*D(3,2)+D(3,1).*D(1,2),...%ok + D(1,1).*D(2,2)-D(2,1).*D(1,2)... %ok + ]; + + if doDerivative + D1 = Mesh.dx1; D2 = Mesh.dx2; D3 = Mesh.dx3; + Z = sparse(size(D1,1), size(D2,2)); + D = @(i,j) sdiag(By(:,(j-1)*3+i)); + + dCof{1} = [ Z , D(3,3)*D2-D(2,3)*D3, D(2,2)*D3-D(3,2)*D2]; + dCof{2} = [ Z , D(1,3)*D3-D(3,3)*D1, D(3,2)*D1-D(1,2)*D3]; + dCof{3} = [ Z , D(2,3)*D1-D(1,3)*D2, D(1,2)*D2-D(2,2)*D1]; + dCof{4} = [ D(2,3)*D3-D(3,3)*D2, Z , D(3,1)*D2-D(2,1)*D3]; + dCof{5} = [ D(3,3)*D1-D(1,3)*D3, Z , D(1,1)*D3-D(3,1)*D1]; + dCof{6} = [ D(1,3)*D2-D(2,3)*D1, Z , D(2,1)*D1-D(1,1)*D2]; + dCof{7} = [ D(3,2)*D2-D(2,2)*D3, D(2,1)*D3-D(3,1)*D2, Z]; + dCof{8} = [ D(1,2)*D3-D(3,2)*D1, D(3,1)*D1-D(1,1)*D3, Z]; + dCof{9} = [ D(2,2)*D1-D(1,2)*D2, D(1,1)*D2-D(2,1)*D1, Z]; + end + +else + + dCof = []; + dx1 = Mesh.mfdx1; dx2 = Mesh.mfdx2; dx3 = Mesh.mfdx3; + + cof{1} = @(yc) dx2.D(yc(:,2)).*dx3.D(yc(:,3)) - dx3.D(yc(:,2)).*dx2.D(yc(:,3)); + cof{2} = @(yc) -dx1.D(yc(:,2)).*dx3.D(yc(:,3)) + dx3.D(yc(:,2)).*dx1.D(yc(:,3)); + cof{3} = @(yc) dx1.D(yc(:,2)).*dx2.D(yc(:,3)) - dx2.D(yc(:,2)).*dx1.D(yc(:,3)); + cof{4} = @(yc) -dx2.D(yc(:,1)).*dx3.D(yc(:,3)) + dx3.D(yc(:,1)).*dx2.D(yc(:,3)); + cof{5} = @(yc) dx1.D(yc(:,1)).*dx3.D(yc(:,3)) - dx3.D(yc(:,1)).*dx1.D(yc(:,3)); + cof{6} = @(yc) -dx1.D(yc(:,1)).*dx2.D(yc(:,3)) + dx2.D(yc(:,1)).*dx1.D(yc(:,3)); + cof{7} = @(yc) dx2.D(yc(:,1)).*dx3.D(yc(:,2)) - dx3.D(yc(:,1)).*dx2.D(yc(:,2)); + cof{8} = @(yc) -dx1.D(yc(:,1)).*dx3.D(yc(:,2)) + dx3.D(yc(:,1)).*dx1.D(yc(:,2)); + cof{9} = @(yc) dx1.D(yc(:,1)).*dx2.D(yc(:,2)) - dx2.D(yc(:,1)).*dx1.D(yc(:,2)); + + if doDerivative + + dx{1}.D = Mesh.mfdx1.D; dx{1}.Dadj = Mesh.mfdx1.Dadj; + dx{2}.D = Mesh.mfdx2.D; dx{2}.Dadj = Mesh.mfdx2.Dadj; + dx{3}.D = Mesh.mfdx3.D; dx{3}.Dadj = Mesh.mfdx3.Dadj; + + % adjoint/transpose of an element of dCof matrix + dCofadjele = @(i,j,k,x,yc) dx{i}.Dadj(dx{j}.D(yc(:,k)).*x); + + dCof.dCofadj{1} = @(x,yc) [zeros(Mesh.nnodes,1); dCofadjele(2,3,3,x,yc)-dCofadjele(3,2,3,x,yc); dCofadjele(3,2,2,x,yc)-dCofadjele(2,3,2,x,yc)]; + dCof.dCofadj{2} = @(x,yc) [zeros(Mesh.nnodes,1); dCofadjele(3,1,3,x,yc)-dCofadjele(1,3,3,x,yc); dCofadjele(1,3,2,x,yc)-dCofadjele(3,1,2,x,yc)]; + dCof.dCofadj{3} = @(x,yc) [zeros(Mesh.nnodes,1); dCofadjele(1,2,3,x,yc)-dCofadjele(2,1,3,x,yc); dCofadjele(2,1,2,x,yc)-dCofadjele(1,2,2,x,yc)]; + + dCof.dCofadj{4} = @(x,yc) [dCofadjele(3,2,3,x,yc)-dCofadjele(2,3,3,x,yc); zeros(Mesh.nnodes,1); dCofadjele(2,3,1,x,yc)-dCofadjele(3,2,1,x,yc)]; + dCof.dCofadj{5} = @(x,yc) [dCofadjele(1,3,3,x,yc)-dCofadjele(3,1,3,x,yc); zeros(Mesh.nnodes,1); dCofadjele(3,1,1,x,yc)-dCofadjele(1,3,1,x,yc)]; + dCof.dCofadj{6} = @(x,yc) [dCofadjele(2,1,3,x,yc)-dCofadjele(1,2,3,x,yc); zeros(Mesh.nnodes,1); dCofadjele(1,2,1,x,yc)-dCofadjele(2,1,1,x,yc)]; + + dCof.dCofadj{7} = @(x,yc) [dCofadjele(2,3,2,x,yc)-dCofadjele(3,2,2,x,yc); dCofadjele(3,2,1,x,yc)-dCofadjele(2,3,1,x,yc); zeros(Mesh.nnodes,1)]; + dCof.dCofadj{8} = @(x,yc) [dCofadjele(3,1,2,x,yc)-dCofadjele(1,3,2,x,yc); dCofadjele(1,3,1,x,yc)-dCofadjele(3,1,1,x,yc); zeros(Mesh.nnodes,1)]; + dCof.dCofadj{9} = @(x,yc) [dCofadjele(1,2,2,x,yc)-dCofadjele(2,1,2,x,yc); dCofadjele(2,1,1,x,yc)-dCofadjele(1,2,1,x,yc); zeros(Mesh.nnodes,1)]; + + % element of dCof matrix + dCofele = @(i,j,k,l,x,yc) dx{i}.D(yc(:,j)).*dx{k}.D(x(:,l)); + + dCof.dCof{1} = @(x,yc) dCofele(3,3,2,2,x,yc)-dCofele(2,3,3,2,x,yc) + dCofele(2,2,3,3,x,yc) - dCofele(3,2,2,3,x,yc); + dCof.dCof{2} = @(x,yc) dCofele(1,3,3,2,x,yc)-dCofele(3,3,1,2,x,yc) + dCofele(3,2,1,3,x,yc) - dCofele(1,2,3,3,x,yc); + dCof.dCof{3} = @(x,yc) dCofele(2,3,1,2,x,yc)-dCofele(1,3,2,2,x,yc) + dCofele(1,2,2,3,x,yc) - dCofele(2,2,1,3,x,yc); + + dCof.dCof{4} = @(x,yc) dCofele(2,3,3,1,x,yc)-dCofele(3,3,2,1,x,yc) + dCofele(3,1,2,3,x,yc) - dCofele(2,1,3,3,x,yc); + dCof.dCof{5} = @(x,yc) dCofele(3,3,1,1,x,yc)-dCofele(1,3,3,1,x,yc) + dCofele(1,1,3,3,x,yc) - dCofele(3,1,1,3,x,yc); + dCof.dCof{6} = @(x,yc) dCofele(1,3,2,1,x,yc)-dCofele(2,3,1,1,x,yc) + dCofele(2,1,1,3,x,yc) - dCofele(1,1,2,3,x,yc); + + dCof.dCof{7} = @(x,yc) dCofele(3,2,2,1,x,yc)-dCofele(2,2,3,1,x,yc) + dCofele(2,1,3,2,x,yc) - dCofele(3,1,2,2,x,yc); + dCof.dCof{8} = @(x,yc) dCofele(1,2,3,1,x,yc)-dCofele(3,2,1,1,x,yc) + dCofele(3,1,1,2,x,yc) - dCofele(1,1,3,2,x,yc); + dCof.dCof{9} = @(x,yc) dCofele(2,2,1,1,x,yc)-dCofele(1,2,2,1,x,yc) + dCofele(1,1,2,2,x,yc) - dCofele(2,1,1,2,x,yc); + + end + +end \ No newline at end of file diff --git a/add-ons/FAIRFEM/examples/EFEM_FEMPIRE_Mice3D.m b/add-ons/FAIRFEM/examples/EFEM_FEMPIRE_Mice3D.m index 87ec57a..02f86fd 100755 --- a/add-ons/FAIRFEM/examples/EFEM_FEMPIRE_Mice3D.m +++ b/add-ons/FAIRFEM/examples/EFEM_FEMPIRE_Mice3D.m @@ -6,12 +6,15 @@ % % ================================================================================== close all; clc; clear; -setup3DmiceData; level = 4; +%setup3DmiceData; +setup3DmouseData; + +level = 4; omega = ML{level}.omega; m = ML{level}.m; viewImage('reset','viewImage','imgmontage','colormap','gray'); imgModel('reset','imgModel','splineInterMex','regularizer','moments','theta',1e-1); -regularizer('reset','regularizer','mbHyperElasticFEM','alpha',1e1); +regularizer('reset','regularizer','mbHyperElasticFEM','alpha',1e2); [T,R] = imgModel('coefficients',ML{level}.T,ML{level}.R,omega,'out',0); Mesh = TetraMesh1(omega,m); @@ -24,6 +27,11 @@ FAIRplotsFEM('reset','mode','FEM','fig',level,'plots',1); FAIRplotsFEM('init',struct('Tc',T,'Rc',R,'Mesh',Mesh)); -yc = GaussNewton(fctn,xc(:),'Plots',@FAIRplotsFEM,'solver','pcg'); +yc = GaussNewton(fctn,xc(:),'Plots',@FAIRplotsFEM,'solver','mbPCG-Jacobi'); + +% matrixfree version +regularizer('reset','regularizer','mfHyperElasticFEM','alpha',1e2); +ycMF = GaussNewton(fctn,xc(:),'Plots',@FAIRplotsFEM,'solver',@FEMPIREsolveGN_PCG); -showResults(ML,yc,'level',level); +relerr = norm(yc-ycMF)/norm(yc) +%showResults(ML,yc,'level',level); diff --git a/add-ons/FAIRFEM/examples/EFEM_FEMSSD_Mice3D.m b/add-ons/FAIRFEM/examples/EFEM_FEMSSD_Mice3D.m new file mode 100644 index 0000000..08a3891 --- /dev/null +++ b/add-ons/FAIRFEM/examples/EFEM_FEMSSD_Mice3D.m @@ -0,0 +1,37 @@ +% ================================================================================== +% (c) Lars Ruthotto 2012/07/26, see FAIR.2 and FAIRcopyright.m. +% http://www.mic.uni-luebeck.de/people/lars-ruthotto.html +% +% Registration of 3D PET of a mice heart using FEM +% +% ================================================================================== +close all; clc; clear; +%setup3DmiceData; +setup3DmouseData; + +level = 4; +omega = ML{level}.omega; m = ML{level}.m; + +viewImage('reset','viewImage','imgmontage','colormap','gray'); +imgModel('reset','imgModel','splineInterMex','regularizer','moments','theta',1e-1); +regularizer('reset','regularizer','mbHyperElasticFEM','alpha',1e2); +[T,R] = imgModel('coefficients',ML{level}.T,ML{level}.R,omega,'out',0); + +Mesh = TetraMesh1(omega,m); +xc = Mesh.xn; +Rc = imgModel(R,omega,Mesh.mfPi(xc,'C')); +Tc = imgModel(T,omega,Mesh.mfPi(xc,'C')); +fctn = @(yc) FEMobjFctn(T,Rc,Mesh,xc,yc(:)); +fctn([]); + +FAIRplotsFEM('reset','mode','FEM','fig',level,'plots',1); +FAIRplotsFEM('init',struct('Tc',T,'Rc',R,'Mesh',Mesh)); + +yc = GaussNewton(fctn,xc(:),'Plots',@FAIRplotsFEM,'solver','mbPCG-Jacobi'); + +% matrixfree version +regularizer('reset','regularizer','mfHyperElasticFEM','alpha',1e2); +ycMF = GaussNewton(fctn,xc(:),'Plots',@FAIRplotsFEM,'solver',@FEMSSDsolveGN_PCG); + +relerr = norm(yc-ycMF)/norm(yc) +%showResults(ML,yc,'level',level); diff --git a/add-ons/FAIRFEM/examples/EFEM_MBvsMF.m b/add-ons/FAIRFEM/examples/EFEM_MBvsMF.m new file mode 100644 index 0000000..eec07de --- /dev/null +++ b/add-ons/FAIRFEM/examples/EFEM_MBvsMF.m @@ -0,0 +1,39 @@ +% ================================================================================== +% (c) Lars Ruthotto 2012/07/26, see FAIR.2 and FAIRcopyright.m. +% http://www.mic.uni-luebeck.de/people/lars-ruthotto.html +% +% Simplified Registration of Cardiac PET +% +% data Gaussian Blobs +% +% ================================================================================== +close all; clc; clear; +setup2DGaussianData; + +distance('reset','distance','SSD'); +imgModel('reset','imgModel','splineInterMex','regularizer','moments','theta',4e-1); +regularizer('reset','regularizer','mbHyperElasticFEM','alpha',1e2); + +minLevel = 3; +maxLevel = 3; + +NPIRpara = optPara('NPIR-GN'); +NPIRpara.Plots = @FAIRplotsFEM; +NPIRpara.lineSearch = @ArmijoDiffeomorphicFEM; +NPIRpara.solver = 'mbPCG-Jacobi'; + +[ySSD_MB,~,hisSSD_MB] = MLIRFEM(ML,'parametric',0,... + 'NPIRobj',@FEMobjFctn,'NPIRpara',NPIRpara,... + 'minLevel',minLevel,'maxLevel',maxLevel); + +regularizer('reset','regularizer','mfHyperElasticFEM','alpha',5e1); +NPIRpara.solver = @FEMSSDsolveGN_PCG; + +[ySSD_MF,~,hisSSD_MF] = MLIRFEM(ML,'parametric',0,... + 'NPIRobj',@FEMobjFctn,'NPIRpara',NPIRpara,... + 'minLevel',minLevel,'maxLevel',maxLevel); + +relerr = norm(ySSD_MB-ySSD_MF)/norm(ySSD_MB) + +%% +%save('EFEM_SSDvsMP','ML','yFEMPIRE','ySSD','hisFEMPIRE','hisSSD'); \ No newline at end of file diff --git a/add-ons/FAIRFEM/examples/EFEM_MP_MBvsMF.m b/add-ons/FAIRFEM/examples/EFEM_MP_MBvsMF.m new file mode 100644 index 0000000..7820ab0 --- /dev/null +++ b/add-ons/FAIRFEM/examples/EFEM_MP_MBvsMF.m @@ -0,0 +1,40 @@ +% ================================================================================== +% (c) Lars Ruthotto 2012/07/26, see FAIR.2 and FAIRcopyright.m. +% http://www.mic.uni-luebeck.de/people/lars-ruthotto.html +% +% Verification code for matrixfree implementation +% +% 'relerr' shows the relative difference between matrixfree and matrix +% based implementations +% +% ================================================================================== +close all; clc; clear; +setup2DGaussianData; + +distance('reset','distance','SSD'); +imgModel('reset','imgModel','splineInterMex','regularizer','moments','theta',4e-1); +regularizer('reset','regularizer','mbHyperElasticFEM','alpha',5e1); + +minLevel = 5; +maxLevel = 5; + +NPIRpara = optPara('NPIR-GN'); +NPIRpara.Plots = @FAIRplotsFEM; +NPIRpara.lineSearch = @ArmijoDiffeomorphicFEM; +NPIRpara.solver = 'mbPCG-Jacobi'; + +[yMP_MB,~,hisMP_MB] = MLIRFEM(ML,'parametric',0,... + 'NPIRobj',@FEMPIREobjFctn,'NPIRpara',NPIRpara,... + 'minLevel',minLevel,'maxLevel',maxLevel); + +regularizer('reset','regularizer','mfHyperElasticFEM','alpha',5e1); +NPIRpara.solver = @FEMPIREsolveGN_PCG;%'PCG-hyperElastic'; + +[yMP_MF,~,hisMP_MF] = MLIRFEM(ML,'parametric',0,... + 'NPIRobj',@FEMPIREobjFctn,'NPIRpara',NPIRpara,... + 'minLevel',minLevel,'maxLevel',maxLevel); + +relerr = norm(yMP_MB-yMP_MF)/norm(yMP_MB) + +%% +%save('EFEM_SSDvsMP','ML','yFEMPIRE','ySSD','hisFEMPIRE','hisSSD'); \ No newline at end of file diff --git a/add-ons/FAIRFEM/examples/EFEM_SSDvsMP.m b/add-ons/FAIRFEM/examples/EFEM_SSDvsMP.m index 842347d..bd108a3 100755 --- a/add-ons/FAIRFEM/examples/EFEM_SSDvsMP.m +++ b/add-ons/FAIRFEM/examples/EFEM_SSDvsMP.m @@ -30,5 +30,6 @@ [ySSD,~,hisSSD] = MLIRFEM(ML,'parametric',0,... 'NPIRobj',@FEMobjFctn,'NPIRpara',NPIRpara,... 'minLevel',minLevel,'maxLevel',maxLevel); + %% save('EFEM_SSDvsMP','ML','yFEMPIRE','ySSD','hisFEMPIRE','hisSSD'); \ No newline at end of file diff --git a/add-ons/FAIRFEM/examples/verifyMF_getdeterminant.m b/add-ons/FAIRFEM/examples/verifyMF_getdeterminant.m new file mode 100644 index 0000000..b4165f3 --- /dev/null +++ b/add-ons/FAIRFEM/examples/verifyMF_getdeterminant.m @@ -0,0 +1,29 @@ +% verify getdeterminant matrix free implementation + +%omega = [0 3 0 5]; m = [2 4]; +%Mesh = TriMesh1(omega,m); +omega = [0 3 0 5 0 2]; m = [2 4 3]; +Mesh = TetraMesh1(omega,m); +xc = Mesh.xn; +% rotate x-y plane by 45 degree +yc(:,1) = (1/sqrt(2))*(xc(:,1) - xc(:,2)); +yc(:,2) = (1/sqrt(2))*(xc(:,1) + xc(:,2)); +yc(:,3) = xc(:,3); +xc = yc; + +[det,dD] = getDeterminant(Mesh,xc,'matrixFree',0); + +[detMF,dDMF] = getDeterminant(Mesh,xc,'matrixFree',1); + +% error in determinant value +err_det = norm(det-detMF); + +% error in dDet +x = rand(size(xc)); +err_dD = norm(dD*x(:) - dDMF.dDet(x)); + +% error in adjoint of dDet +z = rand(Mesh.ntri,1); +err_dDT = norm(dD'*z - dDMF.dDetadj(z)); + + diff --git a/add-ons/FAIRFEM/getBasisGradient.m b/add-ons/FAIRFEM/getBasisGradient.m new file mode 100644 index 0000000..8c0c999 --- /dev/null +++ b/add-ons/FAIRFEM/getBasisGradient.m @@ -0,0 +1,67 @@ +% For details see +% - https://github.com/C4IR/FAIRFEM +%============================================================================== +% +% function [dphi] = getBasisGradient(Mesh) +% +% builds gradient operator for piecewise linear basis functions +% +% Input: +% +% Mesh - description of mesh, struct +% +% Output: +% +% [dphi] - discrete partial derivative operators +% +%============================================================================== + +function [dphi] = getBasisGradient(Mesh) + +xn = Mesh.xn; +vol = Mesh.vol; +if Mesh.dim == 3 +% compute edges + v1 = Mesh.mfPi(xn,1); + v2 = Mesh.mfPi(xn,2); + v3 = Mesh.mfPi(xn,3); + v4 = Mesh.mfPi(xn,4); + e1 = v1 - v4; + e2 = v2 - v4; + e3 = v3 - v4; + % compute inverse transformation to reference element + cofA = [ + e2(:,2).*e3(:,3)-e2(:,3).*e3(:,2), ... + -(e1(:,2).*e3(:,3)-e1(:,3).*e3(:,2)), ... + e1(:,2).*e2(:,3)-e1(:,3).*e2(:,2), ... + -(e2(:,1).*e3(:,3)-e2(:,3).*e3(:,1)), ... + e1(:,1).*e3(:,3)-e1(:,3).*e3(:,1), ... + -(e1(:,1).*e2(:,3)-e1(:,3).*e2(:,1)), ... + e2(:,1).*e3(:,2)-e2(:,2).*e3(:,1), ... + -(e1(:,1).*e3(:,2)-e1(:,2).*e3(:,1)), ... + e1(:,1).*e2(:,2)-e1(:,2).*e2(:,1), ... + ]; + detA = e1(:,1).*cofA(:,1) + e2(:,1).*cofA(:,2)+ e3(:,1).*cofA(:,3); + + + % compute gradients of basis functions + dphi{1} = cofA(:,[1 4 7])./repmat(detA,[1 3]); + dphi{2} = cofA(:,[2 5 8])./repmat(detA,[1 3]); + dphi{3} = cofA(:,[3 6 9])./repmat(detA,[1 3]); + dphi{4} = -dphi{1} - dphi{2} - dphi{3}; + +else + + % compute edges + x1 = Mesh.mfPi(xn,1); + x2 = Mesh.mfPi(xn,2); + x3 = Mesh.mfPi(xn,3); + e1 = x1 - x3; + e2 = x2 - x3; + + % compute gradients of basis functions + dphi{1} = [ e2(:,2) -e2(:,1)]./[2*vol,2*vol]; + dphi{2} = [-e1(:,2) e1(:,1)]./[2*vol,2*vol]; + dphi{3} = -dphi{1} - dphi{2}; + +end \ No newline at end of file diff --git a/add-ons/FAIRFEM/getDeterminant.m b/add-ons/FAIRFEM/getDeterminant.m index 5842232..e6f0323 100755 --- a/add-ons/FAIRFEM/getDeterminant.m +++ b/add-ons/FAIRFEM/getDeterminant.m @@ -23,57 +23,105 @@ dDet = []; if nargin==0, help(mfilename); runMinimalExample; return; end; +matrixFree = regularizer('get','matrixFree'); +dim = Mesh.dim; +yc = reshape(yc,[],dim); + doDerivative = (nargout>1); -for k=1:2:length(varargin), % overwrites default parameter +for k=1:2:length(varargin) % overwrites default parameter eval([varargin{k},'=varargin{',int2str(k+1),'};']); end; -dim = Mesh.dim; -yc = reshape(yc,[],dim); - switch dim case 2 - dx1 = Mesh.dx1; dx2 = Mesh.dx2; + if ~matrixFree + dx1 = Mesh.dx1; dx2 = Mesh.dx2; + + D1y = dx1*yc; + D2y = dx2*yc; + % compute volume + det = D1y(:,1).*D2y(:,2) - D2y(:,1).*D1y(:,2); - D1y = dx1*yc; - D2y = dx2*yc; - % compute volume - det = D1y(:,1).*D2y(:,2) - D2y(:,1).*D1y(:,2); - if doDerivative, - w = [D2y(:,2), D1y(:,2), D1y(:,1), D2y(:,1)]; + if doDerivative + w = [D2y(:,2), D1y(:,2), D1y(:,1), D2y(:,1)]; + dDet =[ + sdiag(w(:,1))*(dx1)-sdiag(w(:,2))*(dx2),... + sdiag(w(:,3))*(dx2)-sdiag(w(:,4))*(dx1)]; + end + else + dx1 = Mesh.mfdx1; dx2 = Mesh.mfdx2; + By = [dx1.D(yc(:,1)) dx2.D(yc(:,1)) dx1.D(yc(:,2)) dx2.D(yc(:,2))]; + det = By(:,1).*By(:,4) - By(:,3) .* By(:,2); - dDet =[ - sdiag(w(:,1))*(dx1) ... - - sdiag(w(:,2))*(dx2),... - sdiag(w(:,3))*(dx2) ... - - sdiag(w(:,4))*(dx1) - ]; + dDet.dDet = @(x) By(:,4).*dx1.D(x(:,1)) - By(:,3).*dx2.D(x(:,1)) + By(:,1).*dx2.D(x(:,2)) - By(:,2).*dx1.D(x(:,2)); + + dDet.dDetadj = @(x) [ ... + dx1.Dadj(By(:,4).*x)-dx2.Dadj(By(:,3).*x);... + -dx1.Dadj(By(:,2).*x)+dx2.Dadj(By(:,1).*x)]; end case 3 - % compute derivative - GRAD = Mesh.GRAD; - By = reshape(GRAD*reshape(yc,[],3),[],9); - D = @(i,j) By(:,(i-1)*3+j); - cof = [ - D(2,2).*D(3,3)-D(2,3).*D(3,2),... - D(2,1).*D(3,3)-D(2,3).*D(3,2),... - D(2,1).*D(3,2)-D(2,2).*D(3,1),... - ]; - - det = sum(By(:,1:3).*cof(:,1:3),2); - if doDerivative, + + if ~matrixFree + % compute derivative + GRAD = Mesh.GRAD; + By = reshape(GRAD*reshape(yc,[],3),[],9); + D = @(i,j) By(:,(j-1)*3+i); % gets partial_i y^j + cof = [ + D(2,2).*D(3,3)-D(3,2).*D(2,3),...%C(1,1) + -D(1,2).*D(3,3)+D(3,2).*D(1,3),...%C(1,2) + D(1,2).*D(2,3)-D(2,2).*D(1,3),...%C(1,3) + -D(2,1).*D(3,3)+D(3,1).*D(2,3),...%C(2,1) + D(1,1).*D(3,3)-D(3,1).*D(1,3),...%C(2,2) + -D(1,1).*D(2,3)+D(2,1).*D(1,3),...%C(2,3) + D(2,1).*D(3,2)-D(3,1).*D(2,2),...%C(3,1) + -D(1,1).*D(3,2)+D(3,1).*D(1,2),...%C(3,2) + D(1,1).*D(2,2)-D(2,1).*D(1,2)... %C(3,3) + ]; + + det = sum(By(:,1:3).*cof(:,1:3),2); + if doDerivative + + D1 = Mesh.dx1; D2 = Mesh.dx2; D3 = Mesh.dx3; + + dDet = [sdiag(cof(:,1))*D1 + sdiag(cof(:,2))*D2 + sdiag(cof(:,3))*D3,... + sdiag(cof(:,4))*D1 + sdiag(cof(:,5))*D2 + sdiag(cof(:,6))*D3,... + sdiag(cof(:,7))*D1 + sdiag(cof(:,8))*D2 + sdiag(cof(:,9))*D3]; + + end + + else - D1 = Mesh.dx1; D2 = Mesh.dx2; D3 = Mesh.dx3; - Z = sparse(size(D1,1), size(D2,2)); - D = @(i,j) sdiag(By(:,(i-1)*3+j)); + dx1 = Mesh.mfdx1; dx2 = Mesh.mfdx2; dx3 = Mesh.mfdx3; - dCof{1} = [ Z , D(3,3)*D2-D(3,2)*D3, D(2,2)*D3-D(2,3)*D2]; - dCof{2} = [ Z , D(3,3)*D1-D(3,2)*D3, D(2,1)*D3-D(2,3)*D2]; - dCof{3} = [ Z , D(3,2)*D1-D(3,1)*D2, D(2,1)*D2-D(2,2)*D1]; + cof11 = @(yc) dx2.D(yc(:,2)).*dx3.D(yc(:,3)) - dx3.D(yc(:,2)).*dx2.D(yc(:,3)); + cof12 = @(yc) -dx1.D(yc(:,2)).*dx3.D(yc(:,3)) + dx3.D(yc(:,2)).*dx1.D(yc(:,3)); + cof13 = @(yc) dx1.D(yc(:,2)).*dx2.D(yc(:,3)) - dx2.D(yc(:,2)).*dx1.D(yc(:,3)); + cof21 = @(yc) -dx2.D(yc(:,1)).*dx3.D(yc(:,3)) + dx3.D(yc(:,1)).*dx2.D(yc(:,3)); + cof22 = @(yc) dx1.D(yc(:,1)).*dx3.D(yc(:,3)) - dx3.D(yc(:,1)).*dx1.D(yc(:,3)); + cof23 = @(yc) -dx1.D(yc(:,1)).*dx2.D(yc(:,3)) + dx2.D(yc(:,1)).*dx1.D(yc(:,3)); + cof31 = @(yc) dx2.D(yc(:,1)).*dx3.D(yc(:,2)) - dx3.D(yc(:,1)).*dx2.D(yc(:,2)); + cof32 = @(yc) -dx1.D(yc(:,1)).*dx3.D(yc(:,2)) + dx3.D(yc(:,1)).*dx1.D(yc(:,2)); + cof33 = @(yc) dx1.D(yc(:,1)).*dx2.D(yc(:,2)) - dx2.D(yc(:,1)).*dx1.D(yc(:,2)); + + det = dx1.D(yc(:,1)).*cof11(yc) + dx2.D(yc(:,1)).*cof12(yc) + dx3.D(yc(:,1)).*cof13(yc); + if doDerivative + + dDet.dDet = @(x) cof11(yc).*dx1.D(x(:,1)) + cof12(yc).*dx2.D(x(:,1)) + cof13(yc).*dx3.D(x(:,1)) +... + cof21(yc).*dx1.D(x(:,2)) + cof22(yc).*dx2.D(x(:,2)) + cof23(yc).*dx3.D(x(:,2)) +... + cof31(yc).*dx1.D(x(:,3)) + cof32(yc).*dx2.D(x(:,3)) + cof33(yc).*dx3.D(x(:,3)); + + dDet.dDetadj = @(x) [ ... + dx1.Dadj(cof11(yc).*x) + dx2.Dadj(cof12(yc).*x) + dx3.Dadj(cof13(yc).*x);... + dx1.Dadj(cof21(yc).*x) + dx2.Dadj(cof22(yc).*x) + dx3.Dadj(cof23(yc).*x);... + dx1.Dadj(cof31(yc).*x) + dx2.Dadj(cof32(yc).*x) + dx3.Dadj(cof33(yc).*x) + ]; + + + end - dDet = [sdiag(cof(:,1))*D1 + sdiag(cof(:,2))*D2 + sdiag(cof(:,3))*D3,Z,Z]... - + sdiag(By(:,1))*dCof{1} + sdiag(By(:,2))*dCof{2} + sdiag(By(:,3))*dCof{3}; end + + otherwise error('only dim=2,3.'); end diff --git a/add-ons/FAIRFEM/hyperElasticFEM.m b/add-ons/FAIRFEM/hyperElasticFEM.m index 7944069..4450e58 100755 --- a/add-ons/FAIRFEM/hyperElasticFEM.m +++ b/add-ons/FAIRFEM/hyperElasticFEM.m @@ -79,7 +79,7 @@ % ===================================================================== % area regularization if dim==3, - [cof, dCof] = cofactor3D(By,Mesh,doDerivative); + [cof, dCof] = cofactor3D(By,Mesh,doDerivative,matrixFree); % compute areas area = [ sum(cof(:,[1 4 7]).^2,2); @@ -103,7 +103,6 @@ d2Sarea = d2Sarea + dA'*sdiag(alphaArea*vol.*d2H(:,i))*dA; end clear dA dH d2H -% dCof{4:9} = []; end else @@ -127,10 +126,15 @@ det = By(:,1).*cof(:,1)+By(:,2).*cof(:,2)+By(:,3).*cof(:,3); if doDerivative, - Z = sparse(size(D1,1),size(D1,2)); + %Z = sparse(size(D1,1),size(D1,2)); % simple product rule - dDet = [sdiag(cof(:,1))*D1 + sdiag(cof(:,2))*D2 + sdiag(cof(:,3))*D3,Z,Z]... - + sdiag(By(:,1))*dCof{1} + sdiag(By(:,2))*dCof{2} + sdiag(By(:,3))*dCof{3}; + %dDet = [sdiag(cof(:,1))*D1 + sdiag(cof(:,2))*D2 + sdiag(cof(:,3))*D3,Z,Z]... + % + sdiag(By(:,1))*dCof{1} + sdiag(By(:,2))*dCof{2} + sdiag(By(:,3))*dCof{3}; + + dDet = [sdiag(cof(:,1))*D1 + sdiag(cof(:,2))*D2 + sdiag(cof(:,3))*D3,... + sdiag(cof(:,4))*D1 + sdiag(cof(:,5))*D2 + sdiag(cof(:,6))*D3,... + sdiag(cof(:,7))*D1 + sdiag(cof(:,8))*D2 + sdiag(cof(:,9))*D3]; + end end [G,dG,d2G] = psi(det,doDerivative); @@ -152,78 +156,128 @@ else % matrix-free d2S.regularizer = regularizer; - d2S.alpha = alpha; + d2S.alpha = alpha; + d2S.yc = yc; + d2S.solver = 'PCG-hyperElastic';%@FEMMultiGridSolveHyper; + % code only diffusion part for B d2S.By = @(u,Mesh) Mesh.mfGRAD.D(u); d2S.BTy = @(u,Mesh) Mesh.mfGRAD.Dadj(u); d2S.B = @(Mesh) Mesh.GRAD; + % give seperate handles for diagonals of length, area and volume d2S.diagLength = @(Mesh) getDiagLength(Mesh,alphaLength); - d2S.diagArea = [];%@(Mesh) getDiagArea(Mesh,alphaLength); + d2S.diagArea = @(Mesh) getDiagArea(Mesh,yc,alphaArea); d2S.diagVol = @(Mesh) getDiagVolume(Mesh,yc,alphaVolume); - d2S.diag = @(Mesh) d2S.diagLength(Mesh) + d2S.diagVol(Mesh); - d2S.solver = @FEMMultiGridSolveHyper; - d2S.d2S = @(uc,Mesh) ... - alphaLength * ... - diffusionOperator(Mesh.vol *Mesh.mfB(uc,Mesh,'By'),Mesh,'BTy') ... - + alphaVolume * ... - volumeOperator(yc,Mesh,uc, Mesh.vol); % length regularizer dSlength = transpose( alphaLength* d2S.BTy(repmat(Mesh.vol,[dim^2,1]).*d2S.By(uc,Mesh),Mesh) ); Slength = .5*dSlength*uc; - % volume term - dx1 = Mesh.mfdx1; dx2 = Mesh.mfdx2; - yc = reshape(yc,[],2); - By = [dx1.D(yc(:,1)) dx2.D(yc(:,1)) dx1.D(yc(:,2)) dx2.D(yc(:,2))]; - det = By(:,1).*By(:,4) - By(:,3) .* By(:,2); - [G,dG,d2G] = psi(det,doDerivative); - Svol = alphaVolume*sum(Mesh.vol.*G); + if doDerivative + d2Slength = @(uc) alphaLength * d2S.BTy(repmat(Mesh.vol,[dim^2,1]).*d2S.By(uc,Mesh),Mesh); + end - dDet = @(x) [ ... - dx1.Dadj(By(:,4).*x)-dx2.Dadj(By(:,3).*x);... - -dx1.Dadj(By(:,2).*x)+dx2.Dadj(By(:,1).*x)]; - dSvol = alphaVolume * transpose(dDet(dG.*Mesh.vol) ); + % area regularizer + if dim==3 + yc = reshape(yc,[],dim); + [cof, dCof] = cofactor3D([],Mesh,doDerivative,matrixFree); + % compute areas + area = [ + cof{1}(yc).^2+cof{4}(yc).^2+cof{7}(yc).^2; + cof{2}(yc).^2+cof{5}(yc).^2+cof{8}(yc).^2; + cof{3}(yc).^2+cof{6}(yc).^2+cof{9}(yc).^2; + ]; + % compute penalty + [H,dH,d2H] = phiDW(area,doDerivative); + % compute area regularizer + Sarea = alphaArea * repmat(vol,[3 1])'*H; + + if doDerivative + dH = reshape(dH,[],3); + d2H = reshape(d2H,[],3); + dAadj = @(x,i) 2*(dCof.dCofadj{i}(x.*cof{i}(yc),yc) + dCof.dCofadj{i+3}(x.*cof{i+3}(yc),yc) + dCof.dCofadj{i+6}(x.*cof{i+6}(yc),yc)); + dSarea = zeros(1,Mesh.nnodes*dim); + for i=1:3 + dSarea = dSarea + alphaArea*dAadj(vol.*dH(:,i),i)'; + end + + dA = @(x,i) 2*(cof{i}(yc).*dCof.dCof{i}(x,yc) + cof{i+3}(yc).*dCof.dCof{i+3}(x,yc) + cof{i+6}(yc).*dCof.dCof{i+6}(x,yc)); + d2Sarea = @(x) dAadj(alphaArea*vol.*d2H(:,1).*dA(reshape(x,[],3),1),1) + dAadj(alphaArea*vol.*d2H(:,2).*dA(reshape(x,[],3),2),2) + dAadj(alphaArea*vol.*d2H(:,3).*dA(reshape(x,[],3),3),3); + + end + + else + Sarea = 0; + dSarea = 0; + d2Sarea = 0; + end + % volume regularizer + if dim==2 + % volume term + dx1 = Mesh.mfdx1; dx2 = Mesh.mfdx2; + yc = reshape(yc,[],2); + By = [dx1.D(yc(:,1)) dx2.D(yc(:,1)) dx1.D(yc(:,2)) dx2.D(yc(:,2))]; + det = By(:,1).*By(:,4) - By(:,3) .* By(:,2); + [G,dG,d2G] = psi(det,doDerivative); + Svolume = alphaVolume*sum(Mesh.vol.*G); - Sc = Slength + Svol; - dS = dSlength + dSvol; + if doDerivative + + dDet = @(x) By(:,4).*dx1.D(x(:,1))- By(:,3).*dx2.D(x(:,1))... + -By(:,2).*dx1.D(x(:,2))+By(:,1).*dx2.D(x(:,2)); + + dDetadj = @(x) [ ... + dx1.Dadj(By(:,4).*x)-dx2.Dadj(By(:,3).*x);... + -dx1.Dadj(By(:,2).*x)+dx2.Dadj(By(:,1).*x)]; + + dSvolume = alphaVolume * transpose(dDetadj(dG.*Mesh.vol)); + + d2Svolume = @(x) dDetadj(vol.*d2G.*dDet(reshape(x,[],2))); + + end + + else + dx1 = Mesh.mfdx1; dx2 = Mesh.mfdx2; dx3 = Mesh.mfdx3; + det = dx1.D(yc(:,1)).*cof{1}(yc) + dx2.D(yc(:,1)).*cof{2}(yc) + dx3.D(yc(:,1)).*cof{3}(yc); + [G,dG,d2G] = psi(det,doDerivative); + Svolume = alphaVolume*sum(Mesh.vol.*G); + + if doDerivative + + dDet = @(x) cof{1}(yc).*dx1.D(x(:,1)) + cof{2}(yc).*dx2.D(x(:,1)) + cof{3}(yc).*dx3.D(x(:,1)) +... + cof{4}(yc).*dx1.D(x(:,2)) + cof{5}(yc).*dx2.D(x(:,2)) + cof{6}(yc).*dx3.D(x(:,2)) +... + cof{7}(yc).*dx1.D(x(:,3)) + cof{8}(yc).*dx2.D(x(:,3)) + cof{9}(yc).*dx3.D(x(:,3)); + + dDetadj = @(x) [ ... + dx1.Dadj(cof{1}(yc).*x) + dx2.Dadj(cof{2}(yc).*x) + dx3.Dadj(cof{3}(yc).*x);... + dx1.Dadj(cof{4}(yc).*x) + dx2.Dadj(cof{5}(yc).*x) + dx3.Dadj(cof{6}(yc).*x);... + dx1.Dadj(cof{7}(yc).*x) + dx2.Dadj(cof{8}(yc).*x) + dx3.Dadj(cof{9}(yc).*x) + ]; + + dSvolume = alphaVolume * transpose(dDetadj(dG.*Mesh.vol)); + + d2Svolume = @(x) alphaVolume * dDetadj(vol.*d2G.*dDet(reshape(x,[],3))); + + end + + end + Sc = Slength +Sarea + Svolume; - d2S.Dy = By; - d2S.d2G = d2G; + if doDerivative + if dim==3 + dS = dSlength + dSarea + dSvolume; + d2S.d2S = @(x) d2Slength(x) + d2Sarea(x) + d2Svolume(x); + d2S.diag = @(yc) d2S.diagLength(Mesh) + d2S.diagArea(Mesh) + d2S.diagVol(Mesh); + elseif dim==2 + dS = dSlength + dSvolume; + d2S.d2S = @(x) d2Slength(x) + d2Svolume(x); + d2S.diag = @(yc) d2S.diagLength(Mesh) + d2S.diagVol(Mesh); + end + end -end - - function [cof,dCof] = cofactor3D(By,Mesh,doDerivative) -dCof = []; -D = @(i,j) By(:,(j-1)*3+i); % gets partial_i y^j -cof = [ - D(2,2).*D(3,3)-D(3,2).*D(2,3),...%ok - -D(1,2).*D(3,3)+D(3,2).*D(1,3),...%trouble - D(1,2).*D(2,3)-D(2,2).*D(1,3),...%ok - -D(2,1).*D(3,3)+D(3,1).*D(2,3),...%ok - D(1,1).*D(3,3)-D(3,1).*D(1,3),...%ok - -D(1,1).*D(2,3)+D(2,1).*D(1,3),...%ok - D(2,1).*D(3,2)-D(3,1).*D(2,2),...%ok - -D(1,1).*D(3,2)+D(3,1).*D(1,2),...%ok - D(1,1).*D(2,2)-D(2,1).*D(1,2)... %ok - ]; - -if doDerivative - D1 = Mesh.dx1; D2 = Mesh.dx2; D3 = Mesh.dx3; - Z = sparse(size(D1,1), size(D2,2)); - D = @(i,j) sdiag(By(:,(j-1)*3+i)); - dCof{1} = [ Z , D(3,3)*D2-D(2,3)*D3, D(2,2)*D3-D(3,2)*D2]; - dCof{2} = [ Z , D(1,3)*D3-D(3,3)*D1, D(3,2)*D1-D(1,2)*D3]; - dCof{3} = [ Z , D(2,3)*D1-D(1,3)*D2, D(1,2)*D2-D(2,2)*D1]; - dCof{4} = [ D(2,3)*D3-D(3,3)*D2, Z , D(3,1)*D2-D(2,1)*D3]; - dCof{5} = [ D(3,3)*D1-D(1,3)*D3, Z , D(1,1)*D3-D(3,1)*D1]; - dCof{6} = [ D(1,3)*D2-D(2,3)*D1, Z , D(2,1)*D1-D(1,1)*D2]; - dCof{7} = [ D(3,2)*D2-D(2,2)*D3 , D(2,1)*D3-D(3,1)*D2, Z]; - dCof{8} = [ D(1,2)*D3-D(3,2)*D1 , D(3,1)*D1-D(1,1)*D3, Z]; - dCof{9} = [ D(2,2)*D1-D(1,2)*D2 , D(1,1)*D2-D(2,1)*D1, Z]; end @@ -231,56 +285,152 @@ dim = Mesh.dim; vol = Mesh.vol; if dim==2 - xn = Mesh.xn; - % get edges - e1 = Mesh.mfP1(xn) - Mesh.mfP3(xn); - e2 = Mesh.mfP2(xn) - Mesh.mfP3(xn); - e3 = Mesh.mfP1(xn) - Mesh.mfP2(xn); - % compute gradients of basis functions - dphi1 = ([ e2(:,2) -e2(:,1)]./vol); - dphi2 = ([-e1(:,2) e1(:,1)]./vol); - dphi3 = ([ e3(:,2) -e3(:,1)]./vol); + dphi = Mesh.dphi; + dphi1 = dphi{1}; dphi2 = dphi{2}; dphi3 = dphi{3}; % compute diagonal of volume regularizer % MB : dDet = [sdiag(By(:,4)), sdiag(-By(:,3)) , sdiag(-By(:,2)), sdiag(By(:,1))] * B; - [dx1, dx2] = getGradientMatrixFEM(Mesh); + [dx1, dx2] = getGradientMatrixFEM(Mesh,1); yc = reshape(yc,[],2); By = [dx1.D(yc(:,1)) dx2.D(yc(:,1)) dx1.D(yc(:,2)) dx2.D(yc(:,2))]; det = By(:,1).*By(:,4) - By(:,3) .* By(:,2); [~,~,d2G] = psi(det,1); % get boundaries - Dxi = @(i,j) Mesh.mfP1(d2G.*(By(:,j).*dphi1(:,i)).^2) ... - + Mesh.mfP2(d2G.*(By(:,j).*dphi2(:,i)).^2) ... - + Mesh.mfP3(d2G.*(By(:,j).*dphi3(:,i)).^2); % diagonal of Dxi'*Dxi + Dxi = @(i,j) Mesh.mfPi(vol.*d2G.*(By(:,j).*dphi1(:,i)).^2,1) ... + + Mesh.mfPi(vol.*d2G.*(By(:,j).*dphi2(:,i)).^2,2) ... + + Mesh.mfPi(vol.*d2G.*(By(:,j).*dphi3(:,i)).^2,3); % diagonal of Dxi'*Dxi + + Dxy = @(i,j,k,l) Mesh.mfPi(vol.*d2G.*(By(:,j).*dphi1(:,i)).*(By(:,l).*dphi1(:,k)),1) ... % byproduct terms for verifications + + Mesh.mfPi(vol.*d2G.*(By(:,j).*dphi2(:,i)).*(By(:,l).*dphi2(:,k)),2) ... + + Mesh.mfPi(vol.*d2G.*(By(:,j).*dphi3(:,i)).*(By(:,l).*dphi3(:,k)),3); % diagonal of Dxi'*Dxi - D = [Dxi(1,4)+ Dxi(2,3) ;Dxi(1,2)+ Dxi(2,1)]; + D = [Dxi(1,4)+ Dxi(2,3) - 2*Dxy(1,4,2,3) ;Dxi(1,2)+ Dxi(2,1) - 2*Dxy(1,2,2,1)]; + else - error('nyi'); + dphi = Mesh.dphi; + dphi1 = dphi{1}; dphi2 = dphi{2}; dphi3 = dphi{3}; dphi4 = dphi{4}; + + [cof, dCof] = cofactor3D([],Mesh,0,1); + yc = reshape(yc,[],3); + dx1 = Mesh.mfdx1; dx2 = Mesh.mfdx2; dx3 = Mesh.mfdx3; + det = dx1.D(yc(:,1)).*cof{1}(yc) + dx2.D(yc(:,1)).*cof{2}(yc) + dx3.D(yc(:,1)).*cof{3}(yc); + [~,~,d2G] = psi(det,1); + + % get boundaries + Dxi = @(i,j) Mesh.mfPi(vol.*d2G.*(cof{j}(yc).*dphi1(:,i)).^2,1) ... + + Mesh.mfPi(vol.*d2G.*(cof{j}(yc).*dphi2(:,i)).^2,2) ... + + Mesh.mfPi(vol.*d2G.*(cof{j}(yc).*dphi3(:,i)).^2,3) ... + + Mesh.mfPi(vol.*d2G.*(cof{j}(yc).*dphi4(:,i)).^2,4); % diagonal of Dxi'*Dxi + + Dxy = @(i,j,k,l) Mesh.mfPi(vol.*d2G.*(cof{j}(yc).*dphi1(:,i)).*(cof{l}(yc).*dphi1(:,k)),1) ... % byproduct terms for verifications + + Mesh.mfPi(vol.*d2G.*(cof{j}(yc).*dphi2(:,i)).*(cof{l}(yc).*dphi2(:,k)),2) ... + + Mesh.mfPi(vol.*d2G.*(cof{j}(yc).*dphi3(:,i)).*(cof{l}(yc).*dphi3(:,k)),3) ... + + Mesh.mfPi(vol.*d2G.*(cof{j}(yc).*dphi4(:,i)).*(cof{l}(yc).*dphi4(:,k)),4); % diagonal of Dxi'*Dxi + + D1 = Dxi(1,1)+ Dxi(2,2) + Dxi(3,3) + 2*Dxy(1,1,2,2) + 2*Dxy(1,1,3,3) + 2*Dxy(2,2,3,3); + D2 = Dxi(1,4)+ Dxi(2,5) + Dxi(3,6) + 2*Dxy(1,4,2,5) + 2*Dxy(1,4,3,6) + 2*Dxy(2,5,3,6); + D3 = Dxi(1,7)+ Dxi(2,8) + Dxi(3,9) + 2*Dxy(1,7,2,8) + 2*Dxy(1,7,3,9) + 2*Dxy(2,8,3,9); + + D = [D1;D2;D3]; + end -D = alphaVolume*vol*D; +D = alphaVolume*D; -% compute d2Svol*x -function By = volumeOperator(yc,Mesh,x,vol) -[dx1, dx2] = getGradientMatrixFEM(Mesh); -% volume regularization -yc = reshape(yc,[],2); -By = [dx1.D(yc(:,1)) dx2.D(yc(:,1)) dx1.D(yc(:,2)) dx2.D(yc(:,2))]; -det = By(:,1).*By(:,4) - By(:,3) .* By(:,2); -[~,~,d2G] = psi(det,1); +function D = getDiagArea(Mesh,yc,alphaArea) +%dim = Mesh.dim; +vol = Mesh.vol; + +dphi = Mesh.dphi; +%dphi1 = dphi{1}; dphi2 = dphi{2}; dphi3 = dphi{3}; dphi4 = dphi{4}; + +[cof, ~] = cofactor3D([],Mesh,0,1); +yc = reshape(yc,[],3); +% compute areas +area = [ + cof{1}(yc).^2+cof{4}(yc).^2+cof{7}(yc).^2; + cof{2}(yc).^2+cof{5}(yc).^2+cof{8}(yc).^2; + cof{3}(yc).^2+cof{6}(yc).^2+cof{9}(yc).^2; + ]; +% compute penalty +[~,~,d2H] = phiDW(area,1); +d2H = reshape(d2H,[],3); + +dx{1}.D = Mesh.mfdx1.D; dx{1}.Dadj = Mesh.mfdx1.Dadj; +dx{2}.D = Mesh.mfdx2.D; dx{2}.Dadj = Mesh.mfdx2.Dadj; +dx{3}.D = Mesh.mfdx3.D; dx{3}.Dadj = Mesh.mfdx3.Dadj; + +D = @(i,j) dx{i}.D(yc(:,j)); +%dA = 2*(cof{i}(yc)*dCof{i}+sdiag(cof(:,i+3))*dCof{i+3}+sdiag(cof(:,i+6))*dCof{i+6}); +%vol.*d2H(:,i))*dA; + +% dCof{1} = [ Z , D(3,3)*D2-D(2,3)*D3, D(2,2)*D3-D(3,2)*D2]; +% dCof{2} = [ Z , D(1,3)*D3-D(3,3)*D1, D(3,2)*D1-D(1,2)*D3]; +% dCof{3} = [ Z , D(2,3)*D1-D(1,3)*D2, D(1,2)*D2-D(2,2)*D1]; +% dCof{4} = [ D(2,3)*D3-D(3,3)*D2, Z , D(3,1)*D2-D(2,1)*D3]; +% dCof{5} = [ D(3,3)*D1-D(1,3)*D3, Z , D(1,1)*D3-D(3,1)*D1]; +% dCof{6} = [ D(1,3)*D2-D(2,3)*D1, Z , D(2,1)*D1-D(1,1)*D2]; +% dCof{7} = [ D(3,2)*D2-D(2,2)*D3, D(2,1)*D3-D(3,1)*D2, Z]; +% dCof{8} = [ D(1,2)*D3-D(3,2)*D1, D(3,1)*D1-D(1,1)*D3, Z]; +% dCof{9} = [ D(2,2)*D1-D(1,2)*D2, D(1,1)*D2-D(2,1)*D1, Z]; +Z = zeros(Mesh.ntri,1); +dCof{1} = @(j) [ Z , D(3,3).*dphi{j}(:,2)-D(2,3).*dphi{j}(:,3), D(2,2).*dphi{j}(:,3)-D(3,2).*dphi{j}(:,2)]; +dCof{2} = @(j) [ Z , D(1,3).*dphi{j}(:,3)-D(3,3).*dphi{j}(:,1), D(3,2).*dphi{j}(:,1)-D(1,2).*dphi{j}(:,3)]; +dCof{3} = @(j) [ Z , D(2,3).*dphi{j}(:,1)-D(1,3).*dphi{j}(:,2), D(1,2).*dphi{j}(:,2)-D(2,2).*dphi{j}(:,1)]; +dCof{4} = @(j) [ D(2,3).*dphi{j}(:,3)-D(3,3).*dphi{j}(:,2) , Z , D(3,1).*dphi{j}(:,2)-D(2,1).*dphi{j}(:,3)]; +dCof{5} = @(j) [ D(3,3).*dphi{j}(:,1)-D(1,3).*dphi{j}(:,3) , Z , D(1,1).*dphi{j}(:,3)-D(3,1).*dphi{j}(:,1)]; +dCof{6} = @(j) [ D(1,3).*dphi{j}(:,2)-D(2,3).*dphi{j}(:,1) , Z , D(2,1).*dphi{j}(:,1)-D(1,1).*dphi{j}(:,2)]; +dCof{7} = @(j) [ D(3,2).*dphi{j}(:,2)-D(2,2).*dphi{j}(:,3) , D(2,1).*dphi{j}(:,3)-D(3,1).*dphi{j}(:,2), Z]; +dCof{8} = @(j) [ D(1,2).*dphi{j}(:,3)-D(3,2).*dphi{j}(:,1) , D(3,1).*dphi{j}(:,1)-D(1,1).*dphi{j}(:,3), Z]; +dCof{9} = @(j) [ D(2,2).*dphi{j}(:,1)-D(1,2).*dphi{j}(:,2) , D(1,1).*dphi{j}(:,2)-D(2,1).*dphi{j}(:,1), Z]; + +% vol.*d2H(:,i))*2*cof{i}(yc)*dCof{i} -% MB : dDet = [sdiag(By(:,4)), sdiag(-By(:,3)) , sdiag(-By(:,2)), sdiag(By(:,1))] * B; +coeff = @(i,j) 4.*vol.*d2H(:,i).*cof{j}(yc).^2; -dDet = @(x) By(:,4).*dx1.D(x(:,1))- By(:,3).*dx2.D(x(:,1))... - -By(:,2).*dx1.D(x(:,2))+By(:,1).*dx2.D(x(:,2)); +Dxi = @(i,j) Mesh.mfPi(coeff(i,j).*dCof{j}(1).^2,1) ... + + Mesh.mfPi(coeff(i,j).*dCof{j}(2).^2,2) ... + + Mesh.mfPi(coeff(i,j).*dCof{j}(3).^2,3) ... + + Mesh.mfPi(coeff(i,j).*dCof{j}(4).^2,4); % diagonal of Dxi'*Dxi + +coeff2 = @(i,j,k) 4.*vol.*d2H(:,i).*cof{j}(yc).*cof{k}(yc); + +Dxy = @(i,j,k) Mesh.mfPi(coeff2(i,j,k).*dCof{j}(1).*dCof{k}(1),1) ... + + Mesh.mfPi(coeff2(i,j,k).*dCof{j}(2).*dCof{k}(2),2) ... + + Mesh.mfPi(coeff2(i,j,k).*dCof{j}(3).*dCof{k}(3),3) ... + + Mesh.mfPi(coeff2(i,j,k).*dCof{j}(4).*dCof{k}(4),4); % diagonal of Dxi'*Dxi + +D = zeros(Mesh.nnodes,3); +for i=1:3 + D = D + Dxi(i,i) + Dxi(i,i+3) + Dxi(i,i+6) + 2*Dxy(i,i,i+3) + 2*Dxy(i,i,i+6) + 2*Dxy(i,i+3,i+6); +end -dDetAdj = @(x) [ ... - dx1.Dadj(By(:,4).*x)-dx2.Dadj(By(:,3).*x);... - -dx1.Dadj(By(:,2).*x)+dx2.Dadj(By(:,1).*x) - ]; -By = dDetAdj(vol.*d2G.*dDet(reshape(x,[],2))); +D = alphaArea*D(:); + + + + +% compute d2Svol*x +% function By = volumeOperator(yc,Mesh,x,vol) +% [dx1, dx2] = getGradientMatrixFEM(Mesh,1); +% % volume regularization +% yc = reshape(yc,[],2); +% By = [dx1.D(yc(:,1)) dx2.D(yc(:,1)) dx1.D(yc(:,2)) dx2.D(yc(:,2))]; +% det = By(:,1).*By(:,4) - By(:,3) .* By(:,2); +% [~,~,d2G] = psi(det,1); +% +% % MB : dDet = [sdiag(By(:,4)), sdiag(-By(:,3)) , sdiag(-By(:,2)), sdiag(By(:,1))] * B; +% +% dDet = @(x) By(:,4).*dx1.D(x(:,1))- By(:,3).*dx2.D(x(:,1))... +% -By(:,2).*dx1.D(x(:,2))+By(:,1).*dx2.D(x(:,2)); +% +% dDetAdj = @(x) [ ... +% dx1.Dadj(By(:,4).*x)-dx2.Dadj(By(:,3).*x);... +% -dx1.Dadj(By(:,2).*x)+dx2.Dadj(By(:,1).*x) +% ]; +% By = dDetAdj(vol.*d2G.*dDet(reshape(x,[],2))); @@ -288,27 +438,30 @@ dim = Mesh.dim; vol = Mesh.vol; if dim==2 - xn = Mesh.xn; - % get edges - e1 = Mesh.mfP1(xn) - Mesh.mfP3(xn); - e2 = Mesh.mfP2(xn) - Mesh.mfP3(xn); - e3 = Mesh.mfP1(xn) - Mesh.mfP2(xn); - % compute gradients of basis functions - dphi1 = ([ e2(:,2) -e2(:,1)]./vol); - dphi2 = ([-e1(:,2) e1(:,1)]./vol); - dphi3 = ([ e3(:,2) -e3(:,1)]./vol); + dphi = Mesh.dphi; + dphi1 = dphi{1}; dphi2 = dphi{2}; dphi3 = dphi{3}; % get boundaries - Dxi = @(i) Mesh.mfP1(dphi1(:,i).^2) + Mesh.mfP2(dphi2(:,i).^2) + Mesh.mfP3(dphi3(:,i).^2); % diagonal of Dx1'*Dx1 - + Dxi = @(i) Mesh.mfPi(vol.*dphi1(:,i).^2,1) + Mesh.mfPi(vol.*dphi2(:,i).^2,2) + Mesh.mfPi(vol.*dphi3(:,i).^2,3); % diagonal of Dx1'*Dx1 + D = Dxi(1)+ Dxi(2); D = [D;D]; else - error('nyi'); -end -D = alphaLength*vol*D(:); + dphi = Mesh.dphi; + dphi1 = dphi{1}; dphi2 = dphi{2}; dphi3 = dphi{3}; dphi4 = dphi{4}; + + % get boundaries + Dxi = @(i) Mesh.mfPi(vol.*dphi1(:,i).^2,1) + ... + Mesh.mfPi(vol.*dphi2(:,i).^2,2) + ... + Mesh.mfPi(vol.*dphi3(:,i).^2,3) + ... + Mesh.mfPi(vol.*dphi4(:,i).^2,4); % diagonal of Dx1'*Dx1 + + D = Dxi(1)+ Dxi(2) + Dxi(3); + D = [D;D;D]; +end +D = alphaLength*D(:); function [G dG d2G] = psi(x,doDerivative) @@ -393,12 +546,30 @@ %========== 3D omega = [0,10,0,8,0 4]; m = [5,6,8]; type = 1; +%omega = [0,5,0,5,0 5]; m = [5,5,5]; type = 1; Mesh = TetraMesh1(omega,m); xn = Mesh.xn; yn = xn + 1e-1* randn(size(xn)); -regOptn = {'alpha',1,'alphaLength',0,'alphaArea',1,'alphaVolume',0}; +regOptn = {'alpha',1,'alphaLength',0,'alphaArea',0,'alphaVolume',1}; fctn = @(y) feval(mfilename,y-xn(:),xn(:),Mesh,regOptn{:},'matrixFree',0); [Sc, dS, d2S] = fctn(yn(:)); checkDerivative(fctn,yn(:)); + +fctn = @(y) feval(mfilename,y-xn(:),xn(:),Mesh,regOptn{:},'matrixFree',1); +[ScMF, dSMF, d2SMF] = fctn(yn(:)); + +% error MB - MF +errdS = norm(dS - dSMF) + +z = rand(3*Mesh.nnodes,1); +errd2S = norm(d2S*z - d2SMF.d2S(z)) + +errDiag = norm(diag(d2S) - d2SMF.diag(yn)) + +ans + + + + diff --git a/add-ons/FAIRFEM/meshes/TetraMesh1.m b/add-ons/FAIRFEM/meshes/TetraMesh1.m index a2a117b..8400498 100755 --- a/add-ons/FAIRFEM/meshes/TetraMesh1.m +++ b/add-ons/FAIRFEM/meshes/TetraMesh1.m @@ -171,6 +171,10 @@ mfdx2 mfdx3 mfGRAD + % =================================== + % Partial derivatives of Basis functions + % =================================== + dphi end properties (Access = private) @@ -193,6 +197,7 @@ mfGRAD_ boundaryIdx_ boundaryProj_ + dphi_ end methods @@ -764,6 +769,13 @@ function runMinimalExample(~) mfGRAD = this.mfGRAD_; end + function dphi = get.dphi(this) + if isempty(this.dphi_), + [this.dphi_] = getBasisGradient(this); + end + dphi = this.dphi_; + end + function idx = get.boundaryIdx(this) if isempty(this.boundaryIdx_), mkvc = @(v) v(:); diff --git a/add-ons/FAIRFEM/meshes/TriMesh1.m b/add-ons/FAIRFEM/meshes/TriMesh1.m index 7b5c935..8185556 100755 --- a/add-ons/FAIRFEM/meshes/TriMesh1.m +++ b/add-ons/FAIRFEM/meshes/TriMesh1.m @@ -172,6 +172,10 @@ mfdx1 mfdx2 mfGRAD + % =================================== + % Partial derivatives of Basis functions + % =================================== + dphi end properties (Access = private) @@ -191,6 +195,7 @@ mfGRAD_ boundaryIdx_ boundaryProj_ + dphi_ end methods @@ -439,6 +444,14 @@ function runMinimalExample(~) end mfGRAD = this.mfGRAD_; end + + function dphi = get.dphi(this) + if isempty(this.dphi_), + [this.dphi_] = getBasisGradient(this); + end + dphi = this.dphi_; + end + function idx = get.boundaryIdx(this) if isempty(this.boundaryIdx_), id = reshape(1:prod(this.m+1),this.m+1); diff --git a/add-ons/FAIRFEM/meshes/TriMesh2.m b/add-ons/FAIRFEM/meshes/TriMesh2.m index 25ee2d7..f913583 100755 --- a/add-ons/FAIRFEM/meshes/TriMesh2.m +++ b/add-ons/FAIRFEM/meshes/TriMesh2.m @@ -172,6 +172,10 @@ mfdx1 mfdx2 mfGRAD + % =================================== + % Partial derivatives of Basis functions + % =================================== + dphi end properties (Access = private) @@ -191,6 +195,7 @@ mfGRAD_ boundaryIdx_ boundaryProj_ + dphi_ end methods @@ -446,6 +451,13 @@ function runMinimalExample(~) mfGRAD = this.mfGRAD_; end + function dphi = get.dphi(this) + if isempty(this.dphi_), + [this.dphi_] = getBasisGradient(this); + end + dphi = this.dphi_; + end + function idx = get.boundaryIdx(this) if isempty(this.boundaryIdx_), id = reshape(1:prod(this.m+1),this.m+1); diff --git a/add-ons/FAIRFEM/meshes/TriMesh3.m b/add-ons/FAIRFEM/meshes/TriMesh3.m index f5bc136..6355be3 100755 --- a/add-ons/FAIRFEM/meshes/TriMesh3.m +++ b/add-ons/FAIRFEM/meshes/TriMesh3.m @@ -176,6 +176,10 @@ mfdx1 mfdx2 mfGRAD + % =================================== + % Partial derivatives of Basis functions + % =================================== + dphi end properties (Access = private) @@ -196,6 +200,7 @@ boundaryIdx_ boundaryProj_ AvN_ + dphi_ end @@ -535,6 +540,15 @@ function runMinimalExample(~) end mfGRAD = this.mfGRAD_; end + + function dphi = get.dphi(this) + if isempty(this.dphi_), + [this.dphi_] = getBasisGradient(this); + end + dphi = this.dphi_; + end + + function idx = get.boundaryIdx(this) if isempty(this.boundaryIdx_), id = reshape(1:prod(this.m+1),this.m+1); diff --git a/kernel/numerics/solveLinearSystem.m b/kernel/numerics/solveLinearSystem.m index a08cbcc..4da85ff 100644 --- a/kernel/numerics/solveLinearSystem.m +++ b/kernel/numerics/solveLinearSystem.m @@ -158,14 +158,14 @@ [dy,flag,relres,iter] = pcg(H,rhs,tolCG,maxIterCG); end - case {'PCG-hyperElastic'} + case {'PCG-hyperElastic'} M = @(x) H.d2D.P((H.d2D.dr'*H.d2D.d2psi*H.d2D.dr)*H.d2D.P(x)); Hoperator = @(x) M(x) + H.d2S.d2S(x,H.omega,H.m,H.d2S.yc); Ddiag = diag(H.d2D.dr'*H.d2D.d2psi*H.d2D.dr); D = H.d2D.P(full(Ddiag)) + H.d2S.diag(H.d2S.yc); Preconditioner = @(x) D.\x; % Jacobi preconditioner [dy,flag,relres,iter] = pcg(Hoperator,rhs,tolCG,maxIterCG,Preconditioner); - + otherwise, keyboard error(1) diff --git a/kernel/regularizers/regularizer.m b/kernel/regularizers/regularizer.m index 29da886..1bed112 100644 --- a/kernel/regularizers/regularizer.m +++ b/kernel/regularizers/regularizer.m @@ -54,7 +54,6 @@ % setup default solver for Gauss-Newton systems if strcmp(task,'set') || strcmp(task,'reset'), - switch method, case 'mfElastic', scheme = 'elastic'; @@ -105,6 +104,12 @@ grid = 'FEM'; solver = 'backslash'; + case 'mfHyperElasticFEM', + scheme = 'hyperElasticFEM'; + matrixFree = 1; + grid = 'FEM'; + solver = 'PCG-hyperElastic'; + case 'mbElasticFEM', scheme = 'elasticFEM'; matrixFree = 0;