http://blog.163.com/zpfzcjndx@126/blog/static/6354568120135196735577/

最近一直在计算向量式有限元的膜单元，笔者真心给Matlab的计算速度给跪了，之前分析梁杆单元的时候还能接受，但现在真心无语了，电脑常跑了一夜，第二天早上却发现计算结果是错的。因此我不得不投向数值计算最经典的Fortan语言。

本文以一个简单的平面桁架算例，如图1，分别用Fortran和Matlab编写相同的分析程序进行计算，它们的计算速度可以说有天壤之别，如图2，分别增大计算循环步数两种语言的计算时间比较，图3为个人计算机的性能配置。这也印证了网上盛传的一句话“MATLAB只是个高级计算器”，当然，决不能因此一棒子拍死MATLAB，它在很多方面的功能还是很强大的，是其他一些软件是无法比拟的。

图1

计算模型

图2 计算速度对比

图3 计算机配置

MATLAB程序~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

clc

clear

Node=xlsread('Data.xls','Node');

Element=xlsread('Data.xls','Element');

Boundary=xlsread('Data.xls','Boundary');

Section=xlsread('Data.xls','Section');

Force=xlsread('Data.xls','Force');

%%读取数据

N_Element=size(Element,1);

N_Node=size(Node,1);

N_Force=size(Force,1);

N_Boundary=size(Boundary,1);

h=1e-4;

alltime=20;

Zeta =0.1;

EA=Section(1,1)*Section(1,2);

%%阻尼系数

C1 = 1+Zeta*h/2;

C2 = 1-Zeta*h/2;

Squre_h=h^2;

Mass=0.1;

%%定义位移

ExtForce=zeros(N_Node,2);

for i=1:N_Force

ExtForce(Force(i,1),:)=Force(i,2:3);

end

ExpB=ones(N_Node,2);

for i=1:N_Boundary

ExpB(Boundary(i,1),:)=Boundary(i,2:3);

end

First=Node;

Second=Node;

%约束标识

BD=ones(N_Node,2)-ExpB;

InternalForce=zeros(N_Node,2);

EleForce=zeros(N_Element,1);

tic

for time=h:h: alltime

UnbalancedForce=InternalForce+time*ExtForce/alltime;

Third=

(2*Second - C2*First + Squre_h*UnbalancedForce/Mass)/C1;

Third=Third.*ExpB+Node.*BD;

InternalForce=zeros(N_Node,2);

for

i=1:N_Element

n1=Element(i,1);

n2=Element(i,2);

Storeforce=EleForce(i,1);

Curforce=Truss2D(Third(n1,:),Third(n2,:),Second(n1,:),Second(n2,:),EA,Storeforce);

EleForce(i)=Curforce(1,5);

InternalForce(n1,:)=InternalForce(n1,:)+Curforce(1,1:2);

InternalForce(n2,:)=InternalForce(n2,:)+Curforce(1,3:4);

end

First=Second;

Second=Third;

end

toc

~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

function

Eleforce=Truss2D(Th1,Th2,Sec1,Sec2,Mat,F)

Eleforce=zeros(1,5);

L=sqrt((Th1(1)-Th2(1))^2+(Th1(2)-Th2(2))^2);

L0=sqrt((Sec1(1)-Sec2(1))^2+(Sec1(2)-Sec2(2))^2);

Eleforce(1,5)=F+Mat*(L-L0)/L0;

Eleforce(1,1)=Eleforce(1,5)*(Th2(1)-Th1(1))/L;

Eleforce(1,2)=Eleforce(1,5)*(Th2(2)-Th1(2))/L;

Eleforce(1,3)=-Eleforce(1,1);

Eleforce(1,4)=-Eleforce(1,2);

end

Fortran程序~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

program main

implicit none

integer,parameter::Nnode=42,Nelement=81,NBoundary=2,Nforce=1

integer i,j,n1,n2

real*8 Node(Nnode,2),Third(Nnode,2),First(Nnode,2),Second(Nnode,2),Timeend,Timebegin

real*8 Section(1,2),EleForce(Nelement),Force(Nforce,3),EA,c1,c2,t,Squre_h

integer

Element(Nelement,2),Boundary(NBoundary,3),Oboundary(Nnode,2)

integer::Gboundary(Nnode,2)=1

real*8::

Gforce(Nnode,2)=0,InternalForce(Nnode,2)=0,UnbalancedForce(Nnode,2)=0,Curforce(5),Storeforce

real*8,parameter:: h=1e-4,zeta=0.1,alltime=10,mass=0.1

interface

function Truss2D(Th1,Th2,Sec1,Sec2,Mat,F)

implicit none

real*8 Truss2D(5),Th1(2),Th2(2),Sec1(2),Sec2(2),Mat,F,L0,L

end function

end interface

!定义函数接口

OPEN(UNIT=7,FILE='Node.dat',STATUS='OLD')

OPEN(UNIT=8,FILE='Element.dat',STATUS='OLD')

OPEN(UNIT=9,FILE='Boundary.dat',STATUS='OLD')

OPEN(UNIT=10,FILE='Section.dat',STATUS='OLD')

OPEN(UNIT=11,FILE='Force.dat',STATUS='OLD')

OPEN(UNIT=12,FILE='Result.dat',STATUS='UNKNOWN') !避开1,2,5,6,

call cpu_time(Timebegin)

do i=1, Nnode

read(7,*) Node(i,:)

end do

close(unit=7)

!放入坐标

do i=1, Nelement

read(8,*) Element(i,:)

end do

close(unit=8)

!放入单元

do i=1, NBoundary

read(9,*) Boundary(i,:)

end do

close(unit=9)

!放入边界条件

read(10,*) Section

close(unit=10)

!放入截面参数

do i=1, Nforce

read(11,*) Force(i,:)

end do

!放入荷载

close(unit=11)

do i=1,NBoundary

j=Boundary(i,1)

Gboundary(j,:)=Boundary(i,2:3)

end do

Oboundary=1-Gboundary

!定义整体边界约束矩阵和逆向约束矩阵

do i=1,Nforce

j=int(Force(i,1))

Gforce(j,:)=Force(i,2:3)

end do

!定义整体力矩阵

Squre_h=h*h

c1 = 1+zeta*h/2

c2 = 1-zeta*h/2

EA=Section(1,1)*Section(1,2)

Second=Node

First=Second

!x(n)

!%x(n-1)

do t=h,alltime,h

UnbalancedForce=InternalForce+t*Gforce/alltime

Third= (2*Second - c2*First +

Squre_h*UnbalancedForce/mass)/c1

Third=Third*Gboundary+Node*Oboundary

InternalForce=0;

do i=1,Nelement

n1=Element(i,1)

n2=Element(i,2)

Storeforce=EleForce(i)

Curforce(:)=Truss2D(Third(n1,:),Third(n2,:),Second(n1,:),Second(n2,:),EA,Storeforce)

EleForce(i)=Curforce(5)

InternalForce(n1,:)=InternalForce(n1,:)+Curforce(1:2)

InternalForce(n2,:)=InternalForce(n2,:)+Curforce(3:4)

end do

First=Second

Second=Third

end do

call cpu_time(Timeend)

write(12,*) Timeend-Timebegin

end program

!定义子程序

function Truss2D(Th1,Th2,Sec1,Sec2,Mat,F)

implicit none

real*8 Truss2D(5),Th1(2),Th2(2),Sec1(2),Sec2(2),Mat,F,L0,L

L=sqrt((Th1(1)-Th2(1))**2+(Th1(2)-Th2(2))**2)

L0=sqrt((Sec1(1)-Sec2(1))**2+(Sec1(2)-Sec2(2))**2)

Truss2D(5)=F+Mat*(L-L0)/L0

Truss2D(1)=Truss2D(5)*(Th2(1)-Th1(1))/L

Truss2D(2)=Truss2D(5)*(Th2(2)-Th1(2))/L

Truss2D(3)=-Truss2D(1)

Truss2D(4)=-Truss2D(2)

end function

转发至微博

转发至微博

阅读(0)|

评论(0)

|

用微信 “扫一扫”

将文章分享到朋友圈。

用易信 “扫一扫”

将文章分享到朋友圈。

喜欢 推荐 0人|

转载

历史上的今天

最近读者

热度

关闭

玩LOFTER，免费冲印20张照片，人人有奖！

评论

this.p={ m:2,

b:2,

loftPermalink:'',

id:'fks_087070085080082064084086095064072084080064081080080065083087',

blogTitle:'【转载】Fortran与Matlab的计算速度对比(Code by myself)',

blogAbstract:'\n

http://blog.163.com/zpfzcjndx@126/blog/static/6354568120135196735577/

最近一直在计算向量式有限元的膜单元，笔者真心给Matlab的计算速度给跪了，之前分析梁杆单元的时候还能接受，但现在真心无语了，电脑常跑了一夜，第二天早上却发现计算结果是错的。因此我不得不投向数值计算最经典的Fortan语言。

\n

本文以一个简单的平面桁架算例，如图1，分别用Fortran和Matlab编写相同的分析程序进行计算，它们的计算速度可以说有天壤之别，如图2，分别增大计算循环步数两种语言的计算时间比较，图3为个

', blogTag:'', blogUrl:'blog/static/213566271201492792140850',

isPublished:1, istop:false, type:0, modifyTime:0,

publishTime:1414416100850,

permalink:'blog/static/213566271201492792140850', commentCount:0,

mainCommentCount:0, recommendCount:0, bsrk:-100, publisherId:0,

recomBlogHome:false, currentRecomBlog:false, attachmentsFileIds:[],

vote:{}, groupInfo:{}, friendstatus:'none',

followstatus:'unFollow', pubSucc:'', visitorProvince:'',

visitorCity:'', visitorNewUser:false, postAddInfo:{}, mset:'000',

mcon:'', srk:-100, remindgoodnightblog:false, isBlackVisitor:false,

isShowYodaoAd:false, hostIntro:'', hmcon:'0',

selfRecomBlogCount:'0', lofter_single:'' }

{list a as x}

{if !!x}

{if

x.visitorName==visitor.userName}

{else}

{/if}

{if x.moveFrom=='wap'} {elseif

x.moveFrom=='iphone'} {elseif

x.moveFrom=='android'} {elseif

x.moveFrom=='mobile'} {/if}

${fn(x.visitorNickname,8)|escape}

{/if} {/list}

{if !!a}

${fn(a.nickname,8)|escape}

${a.selfIntro|escape}{if

great260}${suplement}{/if}

{/if}