本文共 6356 字，大约阅读时间需要 21 分钟。

先从alec jacobson代码入手, 这里求得的G是一个#Face*dim by #V, 所以他求的结果是以按面来算的,然后dim分别是x,y,z上的分量.

关于梯度求法的理论请看

只是说这里把$\frac{e_1}{h_1}$根据高线性质(具体看),拆成了$\frac{e_1}{h_1} = - \frac{e_2}{h_2} - \frac{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: $\frac{e_2}{h_2}$ $-\frac{e_2}{h_2}$ $\frac{e_3}{h_3}$ $-\frac{e_3}{h_3}$

$$\begin{array}{c|cccccc} 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 & (\frac{e_2}{h_2})_x & (-\frac{e_2}{h_2})_x & (\frac{e_3}{h_3})_x & (-\frac{e_3}{h_3})_x & (\frac{e_2}{h_2})_y & (-\frac{e_2}{h_2})_y & (\frac{e_3}{h_3})_y & (-\frac{e_3}{h_3})_y&(\frac{e_2}{h_2})_z & (-\frac{e_2}{h_2})_z & (\frac{e_3}{h_3})_z & (-\frac{e_3}{h_3})_z\\ \end{array}$$

所以矩阵为

$$\begin{array}{c|ccc} & V1& V2 &V3 \\ \hline F_{1x} & (-\frac{e_2}{h_2})_x+(-\frac{e_3}{h_3})_x&(\frac{e_2}{h_2})_x&(\frac{e_3}{h_3})_x\\ F_{1y} & (-\frac{e_2}{h_2})_y+(-\frac{e_3}{h_3})_y&(\frac{e_2}{h_2})_y&(\frac{e_3}{h_3})_y\\ F_{1z} & (-\frac{e_2}{h_2})_z+(-\frac{e_3}{h_3})_z&(\frac{e_2}{h_2})_z&(\frac{e_3}{h_3})_z\\\end{array}$$

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);  endend

转载地址：http://jaxdi.baihongyu.com/