《matlab计算效率1》由会员分享，可在线阅读，更多相关《matlab计算效率1（7页珍藏版）》请在金锄头文库上搜索。
1、Matlab 计算效率的提高参照文On efficient finite element assembly in Matlabby Alejandro A. Ortiz. 用 Matlab进行有限元计算，往往采用稀疏矩阵节约内存，一般是一开始就定义整体刚度矩阵为稀疏矩阵，然后再组集，但需要耗费大量时间（见 Alejandro A.Ortiz） ，而作者先使用全矩阵方式求得每个单元矩阵，然后利用 K=sparse（I,J，X，freedom，freedom ）组集整体刚度矩阵，随着网格的细分，组集效率的提高变得明显。但是作者只举了一个常规有限元三角形单元分析温度场的例子，而我们常常用常规有限元四
2、节点等参单元（更高阶单元）来进行分析问题。需要对该程序进行一些修改。首先说一下 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
3、，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 的块体在重力作用下组集整体刚度矩阵所需时间。新的组集的程序：% C
4、ompute KI=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=;ticfor iel=1:numelemsctr = element(iel,:) ;nn = length(sctr);ke = 0;sctr = element(iel,:); % element connectivity%choose Gauss quadrature rules
5、for elementsW,Q=quadrature(2,GAUSS,2);%Transfrom these Gauss points to global coords for plotting ONLYfor 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 Ksctr=element(iel,:);for i=1:length
6、(sctr)sctrB(2*i-1)=2*sctr(i)-1;sctrB(2*i)=2*sctr(i);endr=s;% -% Loop on Gauss points% -for kk = 1 : size(W,1)B=;pt = Q(kk,:); % quadrature point% JacobianN,dNdxi = lagrange_basis(elemType,pt); % element shape functionsJ0 = node(sctr,:)*dNdxi; % element Jacobian matrix% B matrix%one node can be relat
7、ed to more than one crack!B=B BbmatQ4(pt,elemType,iel);% Stiffness matrix K_c(r:r+63)= K_c(r:r+63)+reshape(B*C*B*W(kk)*det(J0),64,1);end % end of looping on GPs nb=length(sctrB);numl=nb*nb;IJ=repmat(sctrB,nb,1);I(r:r+numl-1)=repmat(sctrB,1,nb);J(r:r+numl-1)=reshape(IJ,numl,1);%I=I;repmat(sctrB,1,nb)
8、; %将耗费更多时间%J=J;reshape(IJ,numl,1);s=s+64;end % end of looping on elementsK=sparse(I,J,K_c,2*numnode,2*numnode);toc相同网格下所耗时间网格 50x50 时一般方法please input the x direction divide:50please input the y direction divide:50Elapsed time is 4.989717 seconds.新方法please input the x direction divide:50please input
9、the y direction divide:50Elapsed time is 2.463251 seconds. 80x80 的网格一般方法please input the x direction divide:80please input the y direction divide:80Elapsed time is 37.003733 seconds.新方法please input the x direction divide:80please input the y direction divide:80Elapsed time is 9.246477 seconds. 100x1
10、00 的网格一般的方法please input the x direction divide:100please input the y direction divide:100Elapsed time is 91.080788 seconds.新方法please input the x direction divide:100please input the y direction divide:100Elapsed time is 18.478487 seconds.200x200 的网格一般方法please input the x direction divide:200please i
11、nput the y direction divide:200Elapsed time is 1750.599756 seconds.新方法please input the x direction divide:200please input the y direction divide:200Elapsed time is 231.757482 seconds.可以看出效果十分明显。如果能够加上并行，效果将大增。这种方法其实和谢老的方法十分相似，过程相对复杂，却赢得了时间。算例 2：边界裂纹扩展有限元分析（平面板的尺寸 1mX2m,裂纹长度 a=0.45m，弹性模量E=1e3Pa，泊松比 n
12、u=0.3）% -% Loop on elements% -ticI=zeros(64*numelem+1600*8,1);J=zeros(64*numelem+1600*8,1);K_c=zeros(64*numelem+1600*8,1);%由于需要考虑加强自由度（见上图），将预设矩阵放大，再计算完所有单元刚度矩阵后，再将多余的数据去掉。即64*numelem+1600*8需根据网格细分来调整s_i=1;for iel = 1 : numelemsctr = element(iel,:); % element connectivitynn = length(sctr); % number
13、of nodes per elementke = 0 ; % elementary stiffness matrix % -% Choose Gauss quadrature rules for elementsif (ismember(iel,split_elem) % split elementorder = 2 ;phi = ls(sctr,1);W,Q = discontQ4quad(order,phi);elseif (ismember(iel,tip_elem) % tip elementorder = 7;phi = ls(sctr,1);nodes = node(sctr,:)
14、;W,Q = disTipQ4quad(order,phi,nodes,xTip);elseif ( any(intersect(tip_nodes,sctr) = 0)% having tip enriched nodesorder = 4 ;W,Q = quadrature(order,GAUSS,2);elseorder = 2 ;W,Q = quadrature(order,GAUSS,2);end% -% -% Transform these Gauss points to global coords% for plotting onlyfor igp = 1 : size(W,1)
15、gpnt = Q(igp,:);N,dNdxi=lagrange_basis(Q4,gpnt);Gpnt = N * node(sctr,:); % global GPq = q;Gpnt;end% -% Stiffness matrix Ke = BT C B% B = Bfem | Bxfem % The nodal parameters are stored in the following manner:% u = u1 v1 u2 v2 . un vn | a1 b1 . an bn;% It means that the additional unknowns are placed at the end
