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

1.what’s stiffness matrix

According to (2.16), stiffness matrix is given by:

Lij=∫Ω−∇ϕi∇ϕjdA

L_{ij}=\int\limits_{\Omega}-\nabla\phi_i\nabla\phi_jdA

2. matlab code analysis

function L = cotmatrix(V,F)% COTMATRIX computes cotangent matrix (laplacian mesh operator), (mass/area% terms already cancelled out)%% L = cotmatrix(V,F)%% Inputs:%   V  #V x 3 matrix of vertex coordinates%   F  #F x 3  matrix of indices of triangle corners% Outputs:%   L  sparse #V x #V matrix of cot weights %% Copyright 2011, Alec Jacobson (jacobson@inf.ethz.ch), Denis Zorin%% See also: cotmatrix3%% 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,:); %i1 = i1'; i2 = i2'; i3 = i3';% #F x 3 matrices of triangle edge vectors, named after opposite vertices%seamanj:这里i1,i2,i3是索引，应该为一个向量，无论列向量，行向量都行v1 = V(i3,:) - V(i2,:);  v2 = V(i1,:) - V(i3,:); v3 = V(i2,:) - V(i1,:);%seamanj:这里v1,v2,v3,对于某个索引来说是个行向量，然后多个索引从垂直方向延伸% computing areas if size(V,2) == 2% 2d vertex datadblA = v1(:,1).*v2(:,2)-v1(:,2).*v2(:,1);% seamanj added: dbl for doubleelseif size(V,2) == 3%n  = cross(v1,v2,2);  dblA  = multinorm(n,2);% area of parallelogram is twice area of triangle% area of parallelogram is || v1 x v2 || n  = cross(v1,v2,2); % THIS DOES MATRIX NORM!!! don't use it!!% dblA  = norm(n,2);% This does correct l2 norm of rowsdblA = (sqrt(sum((n').^2)))';else error('unsupported vertex dimension %d', size(V,2))end% cotangents and diagonal entries for element matricescot12 = -dot(v1,v2,2)./dblA; cot23 = -dot(v2,v3,2)./dblA; cot31 = -dot(v3,v1,2)./dblA;%注意v1和v2向量的起点不在一个点上,算余弦要把一个向量反转, 如果为锐角这里的cot都为正% diag entries computed from the condition that rows of the matrix sum up% to 0% (follows from  the element matrix formula E_{ij} = (v_i dot v_j)/4/A )diag1 = -cot12/2-cot31/2; diag2 = -cot12/2-cot23/2; diag3 = -cot31/2-cot23/2;%diag1 = -diag1; diag2 = -diag2; diag3 = -diag3;% indices of nonzero elements in the matrix for sparse() constructori = [i1 i2 i2 i3 i3 i1  i1 i2 i3];j = [i2 i1 i3 i2 i1 i3  i1 i2 i3];% values corresponding to pairs form (i,j)v = [cot12/2 cot12/2 cot23/2 cot23/2 cot31/2 cot31/2 diag1 diag2 diag3];% for repeated indices (i,j) sparse automatically sums up elements, as we% wantL = sparse(i,j,v,size(V,1),size(V,1));
end

I have already added some comments in the code. Therefore, I would like to pick some important lines to explain.

For one edge in each triangle, say edge e12e_{12}, the corresponding element in stiffness matrix L12L_{12} is like:

L12=∫Ω−∇ϕ1∇ϕ2dA=−A∇ϕ1⋅∇ϕ2=−A∗(−cotα122A)=cotα122

L_{12}=\int\limits_{\Omega}-\nabla\phi_1\nabla\phi _2dA=-A\nabla\phi_1\cdot\nabla\phi_2=-A*(-\frac{cot\alpha_{12}}{2A})=\frac{cot\alpha_{12}}{2}

note: there are some little errors in Jacobson’s original thesis.

At line 52 to 54, diag entries computed from the condition that rows of the matrix sum up to 0(follows from the element matrix formula Eij=(Vi⋅Vj)/4/A)E_{ij} = (V_i \cdot V_j)/4/A ).

I would like to explain why the sum equals 0. As the comment said, the element matrix formula is Eij=(Vi⋅Vj)/4AE_{ij}=(V_i\cdot V_j)/4A, the sum of a row is like: ∑jEij=∑j(Vi⋅Vj)/4A=(Vi⋅∑jVj)/4A\sum\limits_{j}E_{ij}=\sum\limits_{j}(V_i\cdot V_j)/4A=(V_i\cdot \sum\limits_{j}V_j)/4A. As we all know, in a triangle the sum of three edge vectors is zero which reveals ∑jVj=0\sum\limits_{j}V_j=0, proved.

solve stiffness matrix in matlab相关推荐

solve mass matrix in matlab
note: this passage serves for the analysis of Alec Jacobson's thesis 1. what's mass matrix According ...
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的时候发现了一个提问: 这个人在求解多元方程组的时候,遇到了以下问题,即求解时遇 ...
abaqus inp扫盲与提高 *MATRIX GENERATE,STIFFNESS的验证
可以通过inp文件可以直接定义一个job: 通过改变source实现,输入input file的路径即可. 不过报错了.报错的原因是我输入的一个集合名称不存在.实际上我定义了这个集合FIBRE-0,但 ...
MATLAB未定义变量example,小虫求助“ 未定义函数或变量 'Beam_InputData547'。
小虫刚学习MATLAB,现正在学习<Matlab有限元结构动力学分析与工程应用>,在5.4.3瞬态问题分析.例5.7中按照书本附带的源程序运行提示未定义函数或变量 'Beam_Input ...

solve stiffness matrix in matlab

1.what’s stiffness matrix

2. matlab code analysis

solve stiffness matrix in matlab相关推荐

最新文章

热门文章