Matlab计算效率1.docVIP

Matlab 计算效率的提高 参照文《On efficient finite element assembly in Matlab》by Alejandro A. Ortiz. 用Matlab进行有限元计算，往往采用稀疏矩阵节约内存，一般是一开始就定义整体刚度矩阵为稀疏矩阵，然后再组集，但需要耗费大量时间（见Alejandro A.Ortiz），而作者先使用全矩阵方式求得每个单元矩阵，然后利用K=sparse（I,J，X，freedom，freedom）组集整体刚度矩阵，随着网格的细分，组集效率的提高变得明显。但是作者只举了一个常规有限元三角形单元分析温度场的例子，而我们常常用常规有限元四节点等参单元（更高阶单元）来进行分析问题。需要对该程序进行一些修改。 首先说一下 sparse的用法。 一般我们将full矩阵转换为稀疏矩阵直接采用K=sparse（K）。但是采用K=sparse（I,J，X，freedom，freedom）使得生成稀疏矩阵变得灵活。其中I为freedom X freedom的稀疏矩阵K的非零元素的行标，J为列标，X则为K（I，J）的值。这样只需要知道非零元素在整个矩阵K的位置，即可生成稀疏矩阵。 比如 I=[1 3 2 3];J=[1 1 2 4]; X=[1 3 2 -1]; Freedom=5; 则 输入 K=sparse（I,J，X，freedom，freedom） 得到： K = (1,1) 1 (3,1) 3 (2,2) 2 (3,4) -1 有限元中： 一般组集方法： 先设置整体刚度矩阵 K=sparse（freedom，freedom），然后组集整体刚度矩阵 新的组集方法： 根据Alejandro A.Ortiz的组集方法是先确定单元full矩阵在整体矩阵的位置，最后K=sparse（I,J，X，freedom，freedom）组集整体刚度矩阵。 算例1. 下面分别使用两种组集方法测试尺寸为1m x 1m，弹性模量为1，泊松比取0.3，密度取1的块体在重力作用下组集整体刚度矩阵所需时间。 新的组集的程序： % Compute K I=zeros(64*numelem,1); J=zeros(64*numelem,1); K_c=zeros(64*numelem,1); %这里若使用I=[],J=[],K_c=[];在下面循环计算反而要花更多的时间，因为四节点等参单元生成的单元刚度矩阵为8X8,故I维数为64numelem. s=1; q=[]; tic for iel=1:numelem sctr = element(iel,:) ; nn = length(sctr); ke = 0; sctr = element(iel,:); % element connectivity %choose Gauss quadrature rules for elements [W,Q]=quadrature(2,GAUSS,2); %Transfrom these Gauss points to global coords for plotting ONLY for igp = 1:size(W,1) gpnt = Q(igp,:) ; [N,dNdxi] = lagrange_basis(elemType,gpnt) ; Gpnt = N*node(sctr,:) ; q = [q;Gpnt] ; end % Determine the position in the global matrix K sctr=element(iel,:); for i=1:length(sctr) sctrB(2*i-1)=2*sctr(i)-1; sctrB(2*i)=2*sctr(i); end r=s; % --------------------- % Loop on Gauss points % --------------------- for kk = 1 : size(W,1) B=[]; pt = Q(kk,:);