| 查看: 1792 | 回复: 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、整个过程能不能加速。 谢谢 |
回复此楼» 猜你喜欢
博后基金刷到的BUG,图片来的更直观
已经有9人回复
-
去英国的小伙伴儿都在哪儿租的房子呀？
已经有8人回复
-
太卷了
已经有22人回复
-
这个博士要读吗
已经有17人回复
-
论文已接收，但发现修改稿传的是旧版该怎么办？
已经有6人回复
-
有人中过人文社科类的博后特助吗？
已经有4人回复
-
F口信息学部拿面上，大概需要什么样的成果
已经有8人回复
-
最新消息：2024国自然下载文件名变了
已经有28人回复
-
审稿问题：为什么荧光激发波长和紫外吸收波长差的大？
已经有10人回复
-
初始合伙人来啦！（生物试剂耗材标准品）
已经有21人回复
» 本主题相关价值贴推荐，对您同样有帮助:
- 关于影响因子计算方法的疑问
已经有55人回复
- CASTEP中GGA+U计算纳米线
已经有4人回复
- VASP中怎么进行加U计算？
已经有18人回复
- Gaussian 09 计算出来的chk文件无法转化成fchk文件，求高手指教！
已经有9人回复
- Matlab三维图中如何实现x，y坐标名称与各自的坐标轴平行？
已经有4人回复
- 矩阵相乘，逆向求解，这种方法可行吗？
已经有8人回复
- 如何采用castep计算E-V曲线，请专家指示
已经有20人回复
- 固定床传质衡算
已经有4人回复
- 求matlab画一个函数的图
已经有6人回复
- 为什么我的MATLAB文件不能保存成为TXT格式，跪求大虾米们指教
已经有4人回复
- 球形扩散与平面扩散的区别
已经有9人回复
- matlab的函数调用问题
已经有3人回复
- 红外指纹图谱反映的是中药整体质量信息，体现了中药作用的整体性和模糊性。
已经有14人回复
- 关于几何矩阵在有限元中的问题
已经有6人回复
- matlab 规范切 大规模矩阵，大家帮帮忙啊。
已经有3人回复
★
感谢参与，应助指数 +1
xiegangmai: 金币+1, 谢谢参与！ 2012-07-09 00:04:16
感谢参与，应助指数 +1
xiegangmai: 金币+1, 谢谢参与！ 2012-07-09 00:04:16
|
本帖内容被屏蔽 |
2楼2012-07-08 23:39:52
baobiao007
木虫 (职业作家)
中国特色
- 应助: 201 (大学生)
- 金币: 6482.7
- 散金: 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
★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★
xiegangmai: 金币+1, 谢谢参与！ 2012-07-10 23:25:25
11hours: 金币+10, ★有帮助 2012-07-14 07:25:02
xiegangmai: 金币+1, 谢谢参与！ 2012-07-10 23:25:25
11hours: 金币+10, ★有帮助 2012-07-14 07:25:02
|
本帖内容被屏蔽 |
6楼2012-07-10 12:44:04
11hours
铜虫 (小有名气)
- 应助: 1 (幼儿园)
- 金币: 103.2
- 散金: 4
- 帖子: 66
- 在线: 40.8小时
- 虫号: 1513795
- 注册: 2011-11-28
- 性别: GG
- 专业: 固体力学
7楼2012-07-14 07:26:23
seekforunite
木虫 (著名写手)
- 应助: 6 (幼儿园)
- 金币: 12359.4
- 散金: 2074
- 红花: 3
- 帖子: 1768
- 在线: 270.9小时
- 虫号: 919266
- 注册: 2009-12-03
- 性别: GG
- 专业: 电气科学与工程
8楼2012-11-03 10:19:29
匿名
用户注销 (职业作家)
Matlab专家
- 应助: 223 (大学生)
- 金币: 2891.3
- 散金: 5302
- 红花: 59
- 沙发: 10
- 帖子: 3308
- 在线: 945小时
- 虫号: 0
- 注册: 2012-05-01
- 性别: GG
- 专业: 力学
9楼2013-07-23 09:53:34
8