| 查看: 1794 | 回复: 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、整个过程能不能加速。 谢谢 |
回复此楼» 猜你喜欢
博后网上办公系统维护
已经有21人回复
-
自然基金变化？？
已经有25人回复
-
做什么表征可以检测塑料中碳氢氧的含量
已经有7人回复
-
基金评审时，评审专家还回去看申请人代表作之外的文章吗？
已经有11人回复
-
今年什么时候会评啊
已经有8人回复
-
回家两天，不想呆了
已经有13人回复
-
国基在研影响申请结果吗
已经有10人回复
-
对象没有，还非常想读博，难以抉择
已经有51人回复
-
有机物的酸酐如何准确表征
已经有7人回复
-
大佬们，打的核磁氢谱与chemdraw预测的有偏差
已经有5人回复
» 本主题相关价值贴推荐，对您同样有帮助:
- 关于影响因子计算方法的疑问
已经有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人回复
baobiao007
木虫 (职业作家)
中国特色
- 应助: 201 (大学生)
- 金币: 6482.7
- 散金: 557
- 红花: 40
- 帖子: 3050
- 在线: 1009.9小时
- 虫号: 505962
- 注册: 2008-02-18
- 专业: 应用地球物理学
3楼2012-07-09 08:56:02
10