先从alec jacobson代码入手, 这里求得的G是一个#Face*dim by #V, 所以他求的结果是以按面来算的,然后dim分别是x,y,z上的分量.
关于梯度求法的理论请看
http://blog.csdn.net/seamanj/article/details/52070611
只是说这里把
e1h1
e
1
h
1
根据高线性质(具体看http://blog.csdn.net/seamanj/article/details/52075447),拆成了
e1h1=−e2h2−e3h3
e
1
h
1
=
−
e
2
h
2
−
e
3
h
3
分析一下矩阵构成的代码, 最简单假设只有一个三角形, 那么
I : 1 1 1 1 2 2 2 2 3 3 3 3
J : 2 1 3 1 2 1 3 1 2 1 3 1
V:
e2h2
e
2
h
2
−e2h2
−
e
2
h
2
e3h3
e
3
h
3
−e3h3
−
e
3
h
3
所以矩阵为
function [G] = grad(V,F)
% GRAD Compute the numerical gradient operator for triangle or tet meshes
%
% G = grad(V,F)
%
% Inputs:
% V #vertices by dim list of mesh vertex positions
% F #faces by simplex-size list of mesh face indices
% Outputs:
% G #faces*dim by #V Gradient operator
%
% Example:
% L = cotmatrix(V,F);
% G = grad(V,F);
% dblA = doublearea(V,F);
% GMG = -G'*repdiag(diag(sparse(dblA)/2),size(V,2))*G;
%
% % Columns of W are scalar fields
% G = grad(V,F);
% % Compute gradient magnitude for each column in W
% GM = squeeze(sqrt(sum(reshape(G*W,size(F,1),size(V,2),size(W,2)).^2,2)));
%
dim = size(V,2);
ss = size(F,2);
switch ss
case 2
% Edge lengths
len = abs(V(F(:,2))-V(F(:,1)));
% Gradient is just staggered grid finite difference
G = sparse(repmat(1:size(F,1),2,1)',F,[1 -1]./len, size(F,1),size(V,1));
case 3
% append with 0s for convenience
if size(V,2) == 2
V = [V zeros(size(V,1),1)];
end
% Gradient of a scalar function defined on piecewise linear elements (mesh)
% is constant on each triangle i,j,k:
% grad(Xijk) = (Xj-Xi) * (Vi - Vk)^R90 / 2A + (Xk-Xi) * (Vj - Vi)^R90 / 2A
% grad(Xijk) = Xj * (Vi - Vk)^R90 / 2A + Xk * (Vj - Vi)^R90 / 2A +
% -Xi * (Vi - Vk)^R90 / 2A - Xi * (Vj - Vi)^R90 / 2A
% where Xi is the scalar value at vertex i, Vi is the 3D position of vertex
% i, and A is the area of triangle (i,j,k). ^R90 represent a rotation of
% 90 degrees
%
% renaming indices of vertices of triangles for convenience
i1 = F(:,1); i2 = F(:,2); i3 = F(:,3);
% #F x 3 matrices of triangle edge vectors, named after opposite vertices
v32 = V(i3,:) - V(i2,:); v13 = V(i1,:) - V(i3,:); v21 = V(i2,:) - V(i1,:);
% area of parallelogram is twice area of triangle
% area of parallelogram is || v1 x v2 ||
n = cross(v32,v13,2);
% This does correct l2 norm of rows, so that it contains #F list of twice
% triangle areas
dblA = normrow(n);
% now normalize normals to get unit normals
u = normalizerow(n);
% rotate each vector 90 degrees around normal
%eperp21 = bsxfun(@times,normalizerow(cross(u,v21)),normrow(v21)./dblA);
%eperp13 = bsxfun(@times,normalizerow(cross(u,v13)),normrow(v13)./dblA);
eperp21 = bsxfun(@times,cross(u,v21),1./dblA);
eperp13 = bsxfun(@times,cross(u,v13),1./dblA);
%g = ...
% ( ...
% repmat(X(F(:,2)) - X(F(:,1)),1,3).*eperp13 + ...
% repmat(X(F(:,3)) - X(F(:,1)),1,3).*eperp21 ...
% );
G = sparse( ...
[0*size(F,1)+repmat(1:size(F,1),1,4) ...
1*size(F,1)+repmat(1:size(F,1),1,4) ...
2*size(F,1)+repmat(1:size(F,1),1,4)]', ...
repmat([F(:,2);F(:,1);F(:,3);F(:,1)],3,1), ...
[eperp13(:,1);-eperp13(:,1);eperp21(:,1);-eperp21(:,1); ...
eperp13(:,2);-eperp13(:,2);eperp21(:,2);-eperp21(:,2); ...
eperp13(:,3);-eperp13(:,3);eperp21(:,3);-eperp21(:,3)], ...
3*size(F,1), size(V,1));
%% Alternatively
%%
%% f(x) is piecewise-linear function:
%%
%% f(x) = ∑ φi(x) fi, f(x ∈ T) = φi(x) fi + φj(x) fj + φk(x) fk
%% ∇f(x) = ... = ∇φi(x) fi + ∇φj(x) fj + ∇φk(x) fk
%% = ∇φi fi + ∇φj fj + ∇φk) fk
%%
%% ∇φi = 1/hjk ((Vj-Vk)/||Vj-Vk||)^perp =
%% = ||Vj-Vk|| /(2 Aijk) * ((Vj-Vk)/||Vj-Vk||)^perp
%% = 1/(2 Aijk) * (Vj-Vk)^perp
%%
%m = size(F,1);
%eperp32 = bsxfun(@times,cross(u,v32),1./dblA);
%G = sparse( ...
% [0*m + repmat(1:m,1,3) ...
% 1*m + repmat(1:m,1,3) ...
% 2*m + repmat(1:m,1,3)]', ...
% repmat([F(:,1);F(:,2);F(:,3)],3,1), ...
% [eperp32(:,1);eperp13(:,1);eperp21(:,1); ...
% eperp32(:,2);eperp13(:,2);eperp21(:,2); ...
% eperp32(:,3);eperp13(:,3);eperp21(:,3)], ...
% 3*m,size(V,1));
if dim == 2
G = G(1:(size(F,1)*dim),:);
end
% Should be the same as:
% g = ...
% bsxfun(@times,X(F(:,1)),cross(u,v32)) + ...
% bsxfun(@times,X(F(:,2)),cross(u,v13)) + ...
% bsxfun(@times,X(F(:,3)),cross(u,v21));
% g = bsxfun(@rdivide,g,dblA);
case 4
% really dealing with tets
T = F;
% number of dimensions
assert(dim == 3);
% number of vertices
n = size(V,1);
% number of elements
m = size(T,1);
% simplex size
assert(size(T,2) == 4);
% f(x) is piecewise-linear function:
%
% f(x) = ∑ φi(x) fi, f(x ∈ T) = φi(x) fi + φj(x) fj + φk(x) fk + φl(x) fl
% ∇f(x) = ... = ∇φi(x) fi + ∇φj(x) fj + ∇φk(x) fk + ∇φl(x) fl
% = ∇φi fi + ∇φj fj + ∇φk fk + ∇φl fl
%
% ∇φi = 1/hjk = Ajkl / 3V * (Facejkl)^perp
% = Ajkl / 3V * (Vj-Vk)x(Vl-Vk)
% = Ajkl / 3V * Njkl / ||Njkl||
%
% get all faces
F = [ ...
T(:,1) T(:,2) T(:,3); ...
T(:,1) T(:,3) T(:,4); ...
T(:,1) T(:,4) T(:,2); ...
T(:,2) T(:,4) T(:,3)];
% compute areas of each face
A = doublearea(V,F)/2;
N = normalizerow(normals(V,F));
% compute volume of each tet
vol = volume(V,T);
G = sparse( ...
[0*m + repmat(1:m,1,4) ...
1*m + repmat(1:m,1,4) ...
2*m + repmat(1:m,1,4)], ...
repmat([T(:,4);T(:,2);T(:,3);T(:,1)],3,1), ...
repmat(A./(3*repmat(vol,4,1)),3,1).*N(:), ...
3*m,n);
end
end