note: this passage serves for the analysis of Alec Jacobson’s thesis

1. what’s mass matrix

According to (2.50), mass matrix is given by:

Mij=∫ΩϕiϕjdA

M_{ij}=\int_\Omega\phi_i\phi_jdA

2. matlab code analysis

function M = massmatrix(V,F, type)% MASSMATRIX mass matrix for the mesh given by V and F%% M = massmatrix(V,F, type)%% Inputs:%  V  #V x 3 matrix of vertex coordinates%  F  #F x 3  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% Output:%  M  #V by #V sparse mass matrix%% Copyright 2011, Alec Jacobson (jacobson@inf.ethz.ch)%% See also: massmatrix3%% should change code below, so we don't need this transposeif(size(F,1) == 3)warning('F seems to be 3 by #F, it should be #F by 3');endF = F';% renaming indices of vertices of triangles for conveniencei1 = F(1,:); i2 = F(2,:); i3 = F(3,:); %#F x 3 matrices of triangle edge vectors, named after opposite verticesv1 = V(i3,:) - V(i2,:);  v2 = V(i1,:) - V(i3,:); v3 = V(i2,:) - V(i1,:);    % computing the areasif size(V,2) == 2% 2d vertex datadblA = 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 rowsdblA = (sqrt(sum((n').^2)))';else error('unsupported vertex dimension %d', size(V,2))endif 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));%seamanj: 根据Quadrature rules, 对角线上为A_T/6,非对角线上为A_T/12①%注意这里dblA为双倍的三角形面积elseif strcmp(type,'barycentric')% only diagonal elementsi = [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));%the entry M^d_i is the one third the sum of the areas of incident%triangles on vertex i.②elseif strcmp(type,'voronoi')% just ported version of intrinsic code% edges numbered same as opposite verticesFT = 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');%The voronoi mass matrix entry M^d_i for vertex i is the sum of its%corresponding quadrilaterals from all incident triangles.%具体请参照下一个文件③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
end

①

注意这里它采用的是second set of quadrature rules for triangular elements

[http://math2.uncc.edu/~shaodeng/TEACHING/math5172/Lectures/Lect_15.PDF]

②

③

function [M] = massmatrix_intrinsic(l,F,nvert,masstype)% MASSMATRIX_INTRINSIC compute the mass matrix from edge lengths only%% [M] = massmatrix_intrinsic(l,F)%% Inputs:%  l: #F by 3, array of edge lengths of edges opposite each face in F%  F: #F by 3, list of indices of triangle corners%  nvert: number of vertices, only needed to set size%  masstype: full, barycentric, or voronoi% TODO: this is almost identical to massmatrix, % only the area computation is different, need to refactor%% here's a handy line to view mass matrix entries on plot:% text(UV(:,1), UV(:,2),zeros(size(UV,1),1),num2str(M(M>0)))%% Copyright 2011, Alec Jacobson (jacobson@inf.ethz.ch)%% See also: massmatrix%% should change code below, so we don't need this transposeif(size(F,1) == 3)warning('F seems to be 3 by #F, it should be #F by 3');endF = F';% renaming indices of vertices of triangles for conveniencel1 = l(:,1); l2 = l(:,2); l3 = l(:,3);% semiperimeterss = (l1 + l2 + l3)*0.5;% Heron's formula for areadblA = 2*sqrt( s.*(s-l1).*(s-l2).*(s-l3));% renaming indices of vertices of triangles for conveniencei1 = F(1,:); i2 = F(2,:); i3 = F(3,:); if strcmp(masstype,'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];  elseif strcmp(masstype,'barycentric')% only diagonal elementsi = [i1 i2 i3];j = [i1 i2 i3];diag_v = dblA/6.;v = [diag_v,diag_v,diag_v];elseif strcmp(masstype,'voronoi')cosines = [ ...(l(:,3).^2+l(:,2).^2-l(:,1).^2)./(2*l(:,2).*l(:,3)), ...(l(:,1).^2+l(:,3).^2-l(:,2).^2)./(2*l(:,1).*l(:,3)), ...(l(:,1).^2+l(:,2).^2-l(:,3).^2)./(2*l(:,1).*l(:,2))];%seamanj:求cosinebarycentric = cosines.*l;%seamanj:求质心,质心坐标为  a^2(b^2+c^2-a^2),b^2(c^2+a^2-b^2),c^2(a^2+b^2-c^2)参见④normalized_barycentric = barycentric./[sum(barycentric')' sum(barycentric')' sum(barycentric')'];%seamanj:Barycentric coordinates are homogeneous, so(t_1,t_2,t_3)=(ut_1,ut_2,ut_3) 引用④areas = 0.25*sqrt( ...(l(:,1) + l(:,2) - l(:,3)).* ...(l(:,1) - l(:,2) + l(:,3)).* ...(-l(:,1) + l(:,2) + l(:,3)).* ...(l(:,1) + l(:,2) + l(:,3)));partial_triangle_areas = normalized_barycentric.*[areas areas areas];%seamanj:the areas of the triangles ΔA_1A_2P, ΔA_1A_3P, and ΔA_2A_3P are proportional to the barycentric coordinates t_3, t_2, and t_1 of P  引用④  quads = [ (partial_triangle_areas(:,2)+ partial_triangle_areas(:,3))*0.5 ...(partial_triangle_areas(:,1)+ partial_triangle_areas(:,3))*0.5 ...(partial_triangle_areas(:,1)+ partial_triangle_areas(:,2))*0.5];%seamanj:这里条件cosines(:,1)<0当筛选器，注意左右两边都筛选
quads(cosines(:,1)<0,:) = [areas(cosines(:,1)<0,:)*0.5, ...areas(cosines(:,1)<0,:)*0.25, areas(cosines(:,1)<0,:)*0.25];
quads(cosines(:,2)<0,:) = [areas(cosines(:,2)<0,:)*0.25, ...areas(cosines(:,2)<0,:)*0.5, areas(cosines(:,2)<0,:)*0.25];
quads(cosines(:,3)<0,:) = [areas(cosines(:,3)<0,:)*0.25, ...areas(cosines(:,3)<0,:)*0.25, areas(cosines(:,3)<0,:)*0.5];i = [i1 i2 i3];j = [i1 i2 i3];v = reshape(quads,size(quads,1)*3,1);else error('bad mass matrix type')endM = sparse(i,j,v,nvert, nvert);  %这里会做叠加的事
end

④

[http://mathworld.wolfram.com/BarycentricCoordinates.html]
大部分的内容图片已经给出,这里我想说明下为什么质心坐标的公式为
barycentric = cosines.*l;
注意质心到三角形三个顶点的距离相等,所以这里质心也是外接圆的球心
即这里我们看到的

而matlab里面为cosines.*l, 这里我想推导下
将cosines写开来:[(b2+c2−a2)/2bc,(c2+a2−b2)/2ac,(a2+b2−c2)/2ab][(b^2+c^2-a^2)/2bc,(c^2+a^2-b^2)/2ac,(a^2+b^2-c^2)/2ab]
将l写开来:[a,b,c][a,b,c]
那么cosines.*l则为[(b2+c2−a2)∗a/2bc,(c2+a2−b2)∗b/2ac,(a2+b2−c2)∗c/2ab][(b^2+c^2-a^2)*a/2bc,(c^2+a^2-b^2)*b/2ac,(a^2+b^2-c^2)*c/2ab]
再根据

我们将cosines.*l里面的三个分量同时乘以2abc2abc,则得到[(b2+c2−a2)∗a2,(c2+a2−b2)∗b2,(a2+b2−c2)∗c2][(b^2+c^2-a^2)*a^2,(c^2+a^2-b^2)*b^2,(a^2+b^2-c^2)*c^2]

solve mass matrix in matlab相关推荐

solve stiffness matrix in matlab
note: this passage serves for the analysis of Alec Jacobson's thesis 1.what's stiffness matrix Accor ...
matlab读取txt到矩阵,如何在MATLAB中将文本文件中的数据读入矩阵(How to read data from a text file into a matrix in MATLAB)...
如何在MATLAB中将文本文件中的数据读入矩阵(How to read data from a text file into a matrix in MATLAB) 我在将.txt文件读入单个矩阵时遇 ...
多分类问题中混淆矩阵(Confusion Matrix)的Matlab画法
在多分类问题中,有一种很实用的分类问题结果统计图. 比如说多类别文类问题,那么每一个类别分到其他类别都有一些数据,但是分到自己类别的毕竟多,这样计算百分比之后就形成了一个矩阵,如果分类正确率高的话,那 ...
matlab solve 解的范围,matlab怎么解方程，如何规定解的范围？
方法/步骤: 一.解一元方程 1.先举一例,解方程"x^2+100*x+99=0". 在Matlab "Command Window"中输入如下命令:x=sol ...
matlab solve 保留正实数解,matlab solve命令中解的顺序问题
满意答案 ____阿牛 2013.04.06 采纳率:57% 等级:12 已帮助:12175人在matlab里面solve命令主要是用来求解代数方程(即多项式)的解,但是也不是说其它方程一个也 ...
matlab获得solve得到的值,matlab solve函数赋值方程组
matlab函数赋值如何实现? 定义全局变量就可以的,你可以在命令窗口中先分别定义如:a=2,b=5,c=7,d=11,e=12;x=[abcde];f=f(x)回车就OK啦再如:a=30;b=45; ...
MATLAB | solve函数求解析解时不支持分段函数的解决方案
MATLAB符号求解功能居然不能求分段函数??这么离谱的事情你敢信? 离谱的问题遇到一个很神奇的问题,这两天逛CSDN的时候发现了一个提问: 这个人在求解多元方程组的时候,遇到了以下问题,即求解时遇 ...
matlab中ode指令,matlab中ode5函数编写.doc
matlab中ode5函数编写 function varargout = ode45(ode,tspan,y0,options,varargin) %ODE45 Solve non-stiff dif ...
matlab初值微分方程,常微分方程初值问题的MATLAB解法
大连圣亚海洋世界官网-2021年2月7日发(作者:转身之后还是你) 用 Matlab 求常微分方程 < br>(ODE) 的初值问题 (IVP) 本节考虑一阶常微分方程  u   f ...

solve mass matrix in matlab

1. what’s mass matrix

2. matlab code analysis

solve mass matrix in matlab相关推荐

最新文章

热门文章