mass矩阵入门，matlab

function M = massmatrix(V,F, type)
% MASSMATRIX mass matrix for the mesh given by V and F
%
% M = massmatrix(V,F, type)
% M = massmatrix(V,T, type)
%
% Inputs:
% V #V x 3 matrix of vertex coordinates
% F #F x simplex-size matrix of indices of triangle corners
% type string containing type of mass matrix to compute
% ‘full’: full mass matrix for p.w. linear fem
% ‘barycentric’: diagonal lumped mass matrix obtained by summing 1/3
% ‘voronoi’: true voronoi area, except in cases where triangle is obtuse
% then uses 1/2, 1/4, 1/4 {simplex size 3 only}
% Output:
% M #V by #V sparse mass matrix
%
% Copyright 2011, Alec Jacobson (jacobson@inf.ethz.ch)
%

function r = cross2(a,b)
% Optimizes r = cross(a,b,2), that is it computes cross products per row
% Faster than cross if I know that I’m calling it correctly
r =[a(:,2).*b(:,3)-a(:,3).*b(:,2), …
a(:,3).*b(:,1)-a(:,1).*b(:,3), …
a(:,1).*b(:,2)-a(:,2).*b(:,1)];
end

% simplex size
ss = size(F,2);
if nargin<3
switch ss
case 3
type = ‘voronoi’;
case {2,4}
type = ‘barycentric’;
end
end

switch ss
case 2
l = edge_lengths(V,F);
switch type
case {‘voronoi’,‘barycentric’}
M = sparse(F(😃,F(😃,[l l]/size(V,2),size(V,1),size(V,1));
case ‘full’
% renaming indices of vertices of triangles for convenience
i1 = F(:,1); i2 = F(:,2);
i = [i1 i2 i1 i2];
j = [i2 i1 i1 i2];
v = [l/4 l/4 l/2 l/2];
M = sparse(i,j,v,size(V,1),size(V,1));
end
case 3
% should change code below, so we don’t need this transpose
if(size(F,1) == 3)
warning(‘F seems to be 3 by #F, it should be #F by 3’);
end
F = F’;

% 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
v1 = V(i3,:) - V(i2,:);  v2 = V(i1,:) - V(i3,:); v3 = V(i2,:) - V(i1,:);
% computing the areas
if size(V,2) == 2
% 2d vertex data
  dblA = v1(:,1).*v2(:,2)-v1(:,2).*v2(:,1);
elseif size(V,2) == 3
  %n  = cross(v1,v2,2);  dblA  = multinorm(n,2);
  n  = cross(v1,v2,2);  
  
  % dblA  = norm(n,2);
  % This does correct l2 norm of rows
  dblA = (sqrt(sum((n').^2)))';
else 
  error('unsupported vertex dimension %d', size(V,2))
end
if strcmp(type,'full')
    % arrays for matrix assembly using 'sparse'
    % indices and values of the element mass matrix entries in the order 
    % (1,2), (2,1),(2,3), (3,2), (3,1), (1,3) (1,1), (2,2), (3,3);
    i = [i1 i2 i2 i3 i3 i1  i1 i2 i3];
    j = [i2 i1 i3 i2 i1 i3  i1 i2 i3];
    offd_v = dblA/24.;
    diag_v = dblA/12.;
    v = [offd_v,offd_v, offd_v,offd_v, offd_v,offd_v, diag_v,diag_v,diag_v];  
    M = sparse(i,j,v,size(V,1), size(V,1));
elseif strcmp(type,'barycentric')
    % only diagonal elements
    i = [i1 i2 i3];
    j = [i1 i2 i3];
    diag_v = dblA/6.;
    v = [diag_v,diag_v,diag_v];
    M = sparse(i,j,v,size(V,1), size(V,1));
elseif strcmp(type,'voronoi')

  % just ported version of intrinsic code

  % edges numbered same as opposite vertices
  FT = F';
  l = [ ...
    sqrt(sum((V(FT(:,2),:)-V(FT(:,3),:)).^2,2)) ...
    sqrt(sum((V(FT(:,3),:)-V(FT(:,1),:)).^2,2)) ...
    sqrt(sum((V(FT(:,1),:)-V(FT(:,2),:)).^2,2)) ...
    ];
  M = massmatrix_intrinsic(l,F',size(V,1),'voronoi');
else 
    error('bad mass matrix type')
end

% warn if any rows are all zero (probably unreferenced vertices)
if(any(sum(M,2) == 0))
  warning('Some rows have all zeros... probably unreferenced vertices..');
end

case 4
% vertices must be defined in 3D
assert(size(V,2)==3);

% should change code below, so we don't need this transpose
if(size(F,1) == 4)
  warning('F seems to be 4 by #F, it should be #F by 4');
end

switch type
case 'full'
  % e.g., from _Finite Elements for Analysis and Design_,  J. E. Akin, pp
  % 198. Or Zienkiewicz-4
  %
  % Earliest reference I could find is "Numerical Integration over
  % Simplexes and Cones", [Hammer et al. 1956]
  % 
  I = [F(:,[2 3 4 3 4 1 4 1 2 1 2 3 1 2 3 4])];
  J = [F(:,[1 1 1 2 2 2 3 3 3 4 4 4 1 2 3 4])];
  vol = abs(volume(V,F));
  VV = [repmat(vol/20,1,3*4) repmat(vol/10,1,4)];
  M = sparse(I,J,VV,size(V,1),size(V,1));
  % Uniform mid-face quadrature rules don't seem to be quadratically
  % precise
  % I = [F(:,[2 3 4 3 4 1 4 1 2 1 2 3 1 2 3 4])];
  % J = [F(:,[1 1 1 2 2 2 3 3 3 4 4 4 1 2 3 4])];
  % vol = abs(volume(V,F));
  % VV = [repmat(vol/18,1,3*4) repmat(vol/12,1,4)];
  % M = sparse(I,J,VV,size(V,1),size(V,1));
case 'barycentric'
  %a = V(F(:,1),:);
  %b = V(F(:,2),:);
  %c = V(F(:,3),:);
  %d = V(F(:,4),:);
  %% http://en.wikipedia.org/wiki/Tetrahedron#Volume
  %% volume for each tetrahedron
  %v = repmat(abs(dot((a-d),cross2(b-d,c-d),2))./6./4,1,4);
  v = repmat(abs(volume(V,F)),1,4);

  % only diagonal elements
  M = sparse(F(:),F(:),v/4,size(V,1),size(V,1));
case 'voronoi'
  pa = V(F(:,1),:);
  pb = V(F(:,2),:);
  pc = V(F(:,3),:);
  pd = V(F(:,4),:);
  % as if pd is the origin
  a = pa-pd;
  b = pb-pd;
  c = pc-pd;
  % circumcenter:
  % http://en.wikipedia.org/wiki/Tetrahedron#More_vector_formulas_in_a_general_tetrahedron
  cc = pd + bsxfun(@rdivide, ...
    bsxfun(@times,sum(a.*a,2),cross2(b,c)) + ...
    bsxfun(@times,sum(b.*b,2),cross2(c,a)) + ...
    bsxfun(@times,sum(c.*c,2),cross2(a,b)), ...
    2*sum(a.*cross2(b,c),2));
  % get correct sign
  sa = sign(sum((pa-pb).*cross2(pc-pb,pd-pb),2)/6);
  sb = sign(sum((pb-pc).*cross2(pd-pc,pa-pc),2)/6);
  sc = sign(sum((pc-pd).*cross2(pa-pd,pb-pd),2)/6);
  sd = sign(sum((pd-pa).*cross2(pb-pa,pc-pa),2)/6);
  % get area of sub-tet and correct sign
  la = sa.*sum((cc-pb).*cross2(pc-pb,pd-pb),2)/6;
  lb = sb.*sum((cc-pc).*cross2(pd-pc,pa-pc),2)/6;
  lc = sc.*sum((cc-pd).*cross2(pa-pd,pb-pd),2)/6;
  ld = sd.*sum((cc-pa).*cross2(pb-pa,pc-pa),2)/6;
  v = abs(sum(a.*cross2(b,c),2)/6);
  %max(abs((la+lb+lc+ld - v)))
  assert(all((la+lb+lc+ld - v)<10*eps));
  % partial volumes attached to each corner
  pv = [(lb+lc+ld)/3 (lc+ld+la)/3 (ld+la+lb)/3 (la+lb+lc)/3];
  M = sparse(F,F,pv,size(V,1),size(V,1));
otherwise
  error('bad mass matrix type')
end

end
end

clear all;
mesh_source = ‘C:\Users\lixiaohu\Desktop\101028\101028_l\TOOTH_0.obj’;
[v,F] = load_mesh(mesh_source);

massmatrix(v,F,‘barycentric’)

(4592,4592) 0.033451
(4593,4593) 0.018297
(4594,4594) 0.059208
(4595,4595) 0.126
(4596,4596) 0.15465
(4597,4597) 0.24073
(4598,4598) 0.20689
(4599,4599) 0.12835
(4600,4600) 0.1287
(4601,4601) 0.089461
(4602,4602) 0.13076
(4603,4603) 0.12632
(4604,4604) 0.15629
(4605,4605) 0.18538
(4606,4606) 0.045072
(4607,4607) 0.044621
(4608,4608) 0.01189
(4609,4609) 0.081306
(4610,4610) 0.13137
(4611,4611) 0.066875
(4612,4612) 0.070747
(4613,4613) 0.13602
(4614,4614) 0.035903
(4615,4615) 0.089015
(4616,4616) 0.18846
(4617,4617) 0.20401
(4618,4618) 0.14487
(4619,4619) 0.16409
(4620,4620) 0.18081
(4621,4621) 0.26814
(4622,4622) 0.31501
(4623,4623) 0.36441
(4624,4624) 0.20949
(4625,4625) 0.2264
(4626,4626) 0.15011
(4627,4627) 0.19432
(4628,4628) 0.19907
(4629,4629) 0.13658
(4630,4630) 0.31603
(4631,4631) 0.16798
(4632,4632) 0.21317
(4633,4633) 0.12932
(4634,4634) 0.25533
(4635,4635) 0.27766
(4636,4636) 0.12725
(4637,4637) 0.23106
(4638,4638) 0.17924
(4639,4639) 0.14798
(4640,4640) 0.11679
(4641,4641) 0.24793
(4642,4642) 0.21979
(4643,4643) 0.20105
(4644,4644) 0.15137
(4645,4645) 0.064404
(4646,4646) 0.14623
(4647,4647) 0.21603
(4648,4648) 0.17891
(4649,4649) 0.17528
(4650,4650) 0.072994
(4651,4651) 0.19146
(4652,4652) 0.15995
(4653,4653) 0.15465
(4654,4654) 0.076551
(4655,4655) 0.023715
(4656,4656) 0.014295
(4657,4657) 0.025876
(4658,4658) 0.091234
(4659,4659) 0.063501
(4660,4660) 0.11813
(4661,4661) 0.03099
(4662,4662) 0.0047765
(4663,4663) 0.014873
(4664,4664) 0.004479
(4665,4665) 0.010133
(4666,4666) 0.010514
(4667,4667) 0.015758
(4668,4668) 0.035992
(4669,4669) 0.0063828
(4670,4670) 0.016875
(4671,4671) 0.013173
(4672,4672) 0.0073152
(4673,4673) 0.015478
(4674,4674) 0.015032
(4675,4675) 0.0061446
(4676,4676) 0.013652
(4677,4677) 0.0092215
(4678,4678) 0.015346
(4679,4679) 0.01324
(4680,4680) 0.027328
(4681,4681) 0.022611
(4682,4682) 0.037459
(4683,4683) 0.02517
(4684,4684) 0.026819
(4685,4685) 0.02
(4686,4686) 0.016032
(4687,4687) 0.033259
(4688,4688) 0.025533
(4689,4689) 0.031606
(4690,4690) 0.03342