Matlab计算效率1.doc

Matlab 计算效率的提高 参照文On efficient finite element assembly in Matlabby 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 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: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,:); % quadrature point % Jacobian N,dNdxi = lagrange_basis(elemType,pt); % element shape functions J0 = node(sctr,:)*dNdxi; % element Jacobian matrix % B matrix %one node can be related 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); %将耗费更多时间 %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 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. 100x100的网格一般的方法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 input 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，泊松比nu=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 : numelem sctr = element(iel,:); % element connectivity nn = length(sctr); % number of nodes per element ke = 0 ; % elementary stiffness matrix % - % Choose Gauss quadrature rules for elements if (ismember(iel,split_elem) % split element order = 2 ; phi = ls(sctr,1); W,Q = discontQ4quad(order,phi); elseif (ismember(iel,tip_elem) % tip element order = 7; phi = ls(sctr,1); nodes = node(sctr,:); W,Q = disTipQ4quad(order,phi,nodes,xTip); elseif ( any(intersect(tip_nodes,sctr) = 0)% having tip enriched nodes order = 4 ; W,Q = quadrature(order,GAUSS,2); else order = 2 ; W,Q = quadrature(order,GAUSS,2); end % - % - % Transform these Gauss points to global coords % for plotting only for igp = 1 : size(W,1) gpnt = Q(igp,:); N,dNdxi=lagrange_basis(Q4,gpnt); Gpnt = N * node(sctr,:); % global GP q = 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 of vector % u. Hence, the true nodal displacement is simply : % u_x = u(1:2:2*numnode) and u_y = u(2:2:2*numnode) % There are two kinds of element: % 1. Non-enriched one (usual element) : calculer comme dhabitude % 2. Enriched one (fully and partially as well): special traitement % Determine the position in the global matrix K sctrB = assembly(iel,enrich_node,pos); nb=length(sctrB); numl=nb*nb; r=s_i; % - % Loop on Gauss points % - for kk = 1 : size(W,1) pt = Q(kk,:); % quadrature point % B matrix B,J0 = xfemBmatrix(pt,elemType,iel,enrich_node,xCr,xTip,alpha); % Stiffness matrix %ke = ke + B*C*B*W(kk)*det(J0); K_c(r:r+numl-1)= K_c(r:r+numl-1)+reshape(B*C*B*W(kk)*det(J0),numl,1); end % end of looping on GPs IJ=repmat(sctrB,nb,1); I(r:r+numl-1)=repmat(sctrB,1,nb); J(r:r+numl-1)=reshape(IJ,numl,1); s_i=s_i+numl;end% end of looping on elementsnd=find(I=0);%找出I多余的数I(nd)=;J(nd)=;K_c(nd)=;K=sparse(I,J,K_c,total_unknown,total_unknown);Toc组集所耗时间20x40的网格旧方法Elapsed time is 1.368753 seconds新方法Elapsed time is 0.993688 seconds.4080的网格旧方法Elapsed time is 9.9833