| 查看: 2528 | 回复: 8 | ||
| 当前只显示满足指定条件的回帖，点击这里查看本话题的所有回帖 | ||
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、整个过程能不能加速。 谢谢 |
回复此楼» 猜你喜欢
售SCI-T0P文章，我:8O.5.5.1.O.54,科目齐全,可+急
已经有3人回复
-
售SCI文章，我:8O5.5.1.O.54,科目齐全,可+急
已经有3人回复
-
售SCI文章，我:8O.5.5.1O.54,科目全,可十急
已经有3人回复
-
售SCI一区T0P文章，我:8O.55.1.O.5.4,科目齐全,可+急
已经有3人回复
-
售SCI一区T0P文章，我:8.O.55.1.O.5.4,科目全,可+急
已经有3人回复
-
售SCI文章，我:8O.5.5.1O.54,科目全,可十急
已经有3人回复
-
售SCI文章，我:8O5.5.1.O.54,科目齐全,可+急
已经有3人回复
-
售SCI一区文章，我:8O5.5.1.O5.4,科目全,可伽急
已经有3人回复
-
售SCI-T0P文章，我:8O.5.5.1.O.54,科目齐全,可+急
已经有3人回复
-
售SCI一区T0P文章，我:8.O55.1.O.54,科目全,可十急
已经有4人回复
» 本主题相关价值贴推荐，对您同样有帮助:
- 关于影响因子计算方法的疑问
已经有55人回复
- CASTEP中GGA+U计算纳米线
已经有4人回复
- Gaussian 09 计算出来的chk文件无法转化成fchk文件，求高手指教！
已经有9人回复
- Matlab三维图中如何实现x，y坐标名称与各自的坐标轴平行？
已经有4人回复
- 矩阵相乘，逆向求解，这种方法可行吗？
已经有8人回复
- 如何采用castep计算E-V曲线，请专家指示
已经有20人回复
- 固定床传质衡算
已经有4人回复
- 求matlab画一个函数的图
已经有6人回复
- 球形扩散与平面扩散的区别
已经有9人回复
- matlab的函数调用问题
已经有3人回复
- 红外指纹图谱反映的是中药整体质量信息，体现了中药作用的整体性和模糊性。
已经有14人回复
- 关于几何矩阵在有限元中的问题
已经有6人回复
- matlab 规范切 大规模矩阵，大家帮帮忙啊。
已经有3人回复
匿名
用户注销 (职业作家)
Matlab专家
- 应助: 223 (大学生)
- 金币: 2891.3
- 散金: 5302
- 红花: 59
- 沙发: 10
- 帖子: 3308
- 在线: 945小时
- 虫号: 0
- 注册: 2012-05-01
- 性别: GG
- 专业: 力学
9楼2013-07-23 09:53:34
baobiao007
木虫 (职业作家)
中国特色
- 应助: 201 (大学生)
- 金币: 6484.2
- 散金: 557
- 红花: 40
- 帖子: 3050
- 在线: 1009.9小时
- 虫号: 505962
- 注册: 2008-02-18
- 专业: 应用地球物理学
3楼2012-07-09 08:56:02
11hours
铜虫 (小有名气)
- 应助: 1 (幼儿园)
- 金币: 103.2
- 散金: 4
- 帖子: 66
- 在线: 40.8小时
- 虫号: 1513795
- 注册: 2011-11-28
- 性别: GG
- 专业: 固体力学
4楼2012-07-09 10:05:31
11hours
铜虫 (小有名气)
- 应助: 1 (幼儿园)
- 金币: 103.2
- 散金: 4
- 帖子: 66
- 在线: 40.8小时
- 虫号: 1513795
- 注册: 2011-11-28
- 性别: GG
- 专业: 固体力学
5楼2012-07-09 10:06:18
20