扩散问题,其中有应力场作用。在建立整体矩阵时(不知道该叫什么名字)速度很慢。先贴段代码,就这段特别慢。
UK1=zeros(NodeNum,NodeNum);%NodeNum=443 是节点数量
UK2=zeros(NodeNum,NodeNum);
%不知道组装整体矩阵的方法对不对,我是把单元矩阵的元素放到整体矩阵对应的节点位置上,累加。这是不是叫贡献法?
for EN=1:EleNum%EleNum=401单元数量
Bi=subs(B,{x1,x2,x3,x4,y1,y2,y3,y4},...%{C,x C,y}T=[B]{C},{C}是单元的四个节点浓度,B就是求C梯度的矩阵。B中变量有,单元局部坐标xi 和 eta,整体坐标{x1,x2,x3,x4,y1,y2,y3,y4}
{Node(Elem(EN,2),2),Node(Elem(EN,3),2),Node(Elem(EN,4),2),Node(Elem(EN,5),2),...
Node(Elem(EN,2),3),Node(Elem(EN,3),3),Node(Elem(EN,4),3),Node(Elem(EN,5),3)});
Press=[Pressure(Elem(EN,2),2),Pressure(Elem(EN,3),2),Pressure(Elem(EN,4),2),Pressure(Elem(EN,5),2)]';%Pressure是节点应力矩阵。第一列是节点号,第二列是应力值,
%Node()是节点坐标矩阵,第一列是编号,第二列和三列是x y坐标。Elem()是单元矩阵,第一列是单元编号,第2至5列是四个节点的编号。
%这段我是想把每个单元节点的坐标代入B,下来就可以对xi 和 eta积分了。
dK1=Bi'*Bi;
dK2=Bi'*VH/R/T*Bi*Press*A; %VH/R/T都是常数,Press每个载荷步要变——!
%积分K1
[r,c]=size(dK1);
K1=zeros(r,c);
for m=1:r
for n=1:c
dk1=matlabFunction(dK1(m,n));
K1(m,n)=dblquad(dk1,-1,1,-1,1);
end
end
%积分K2
[r,c]=size(dK2);
K2=zeros(r,c);
for m=1:r
for n=1:c
dk2=matlabFunction(dK2(m,n));
K2(m,n)=dblquad(dk2,-1,1,-1,1);
end
end
%组装UK1,UK2
UK1Loc=Elem(EN,2:5); %确定总体矩阵中的存放位置
for i=1:4
for j=1:4
UK1(UK1Loc(1,i),UK1Loc(1,j))=UK1(UK1Loc(1,i),UK1Loc(1,j))+K1(i,j);
UK2(UK1Loc(1,i),UK1Loc(1,j))=UK2(UK1Loc(1,i),UK1Loc(1,j))+K2(i,j);
end
end
%UK组装完毕
end
%%%%%%%%%结束
UK1也就罢了,组装一次也就行了,以后不会变。但UK2中有个Press是每个载荷步要变。
dK2=Bi'*VH/R/T*Bi*Press*A;
[r,c]=size(dK2);
K2=zeros(r,c);
for m=1:r
for n=1:c
dk2=matlabFunction(dK2(m,n));
K2(m,n)=dblquad(dk2,-1,1,-1,1);
end
end
这样每个载荷步都要求一次K2,而每次求都会很慢。
初次自己写有限元程序,大家帮我看看:
1、整体矩阵组装的方法对不对
2、积分有没有更快的方法。
3、整个过程能不能加速。
谢谢 |