应《网络安全法》要求，自2017年10月1日起，未进行实名认证将不得使用互联网跟帖服务。为保障您的帐号能够正常使用，请尽快对帐号进行手机号验证，感谢您的理解与支持！

24小时热门版块排行榜    

小木虫论坛-学术科研互动平台 » 计算模拟区 » 仿真模拟 » FEM基本理论 » matlab 有限元计算扩散问题，建立整体矩阵好慢。大家帮忙看看代码
91/1返回列表
查看: 2136  |  回复: 8
只看楼主 @他人 存档 新回复提醒 (忽略) 收藏    在APP中查看

11hours

铜虫 (小有名气)

[求助] matlab 有限元计算扩散问题，建立整体矩阵好慢。大家帮忙看看代码

扩散问题，其中有应力场作用。在建立整体矩阵时（不知道该叫什么名字）速度很慢。先贴段代码，就这段特别慢。

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、整个过程能不能加速。
谢谢
回复此楼

» 猜你喜欢

» 本主题相关价值贴推荐，对您同样有帮助:
1楼 2012-07-08 23:02:32
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

awenluxm

禁虫 (正式写手)
★
感谢参与，应助指数 +1
xiegangmai: 金币+1, 谢谢参与！ 2012-07-09 00:04:16
本帖内容被屏蔽
2楼2012-07-08 23:39:52
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

baobiao007

木虫 (职业作家)

中国特色
有限元的计算量很大的，完全不适合用matlab来做
赞一下 回复此楼
我同意叔本华的观点，人们投身艺术和科学领域的强烈愿望之一就是逃离痛苦、残酷和枯燥无味的现实生活，逃离自己飘忽不定的七情六欲的桎梏。--爱因斯坦
3楼2012-07-09 08:56:02
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

11hours

铜虫 (小有名气)
引用回帖:
2楼: Originally posted by awenluxm at 2012-07-08 23:39:52
我建议你用通用有限元软件！更成熟

我要计算的载荷步很多，还有些特殊结果通用的没法计算。
赞一下 回复此楼
4楼2012-07-09 10:05:31
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

11hours

铜虫 (小有名气)
引用回帖:
3楼: Originally posted by baobiao007 at 2012-07-09 08:56:02
有限元的计算量很大的，完全不适合用matlab来做

正在研究c来实现。过程基本正确吧？
赞一下 回复此楼
5楼2012-07-09 10:06:18
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

awenluxm

禁虫 (正式写手)
★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★
xiegangmai: 金币+1, 谢谢参与！ 2012-07-10 23:25:25
11hours: 金币+10, ★有帮助 2012-07-14 07:25:02
本帖内容被屏蔽
6楼2012-07-10 12:44:04
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

11hours

铜虫 (小有名气)
引用回帖:
6楼: Originally posted by awenluxm at 2012-07-10 12:44:04
你可以通过通用有限元软件的二次开发来实现
...

做过C 循环调用abaqus。结果不太满意，主要是结果文件后期超级大，拖慢计算速度。
赞一下 回复此楼
7楼2012-07-14 07:26:23
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

seekforunite

木虫 (著名写手)
★ ★
xiegangmai: 金币+2, 谢谢参与 2012-11-04 19:21:33
楼主的问题解决了没？其实采用matlab编写有限元程序时，如果采用先计算出刚度阵元素再逐个累加到总体刚度阵中，则确实是很费时。
解决的方法之一采用稀疏矩阵sparse解决，楼主可以查找sparse的相关的帖子。
赞一下 回复此楼
坚持努力！
8楼2012-11-03 10:19:29
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

匿名

用户注销 (职业作家)

Matlab专家
本帖仅楼主可见
赞一下
9楼2013-07-23 09:53:34
已阅   申请仿真EPI   回复此楼   编辑   查看我的主页
相关版块跳转 我要订阅楼主 11hours 的主题更新
91/1返回列表
普通表情 龙 兔 虎 猫
信息提示
关闭
请填处理意见
关闭