基于Fortran的结构力学位移法编程求解

本代码可用于平面桁架结构的计算。

本文是根据基于Python的结构力学位移法编程求解的另一版本的更新，内容基本相同，背景介绍此处省略。仅展示最终代码的运行过程。

由于Fortran编程语言使用比较困难一些，所以很少涉及，特别是其缺少一些相应的函数调用，过程要更复杂一些，而且调试过程也较麻烦。但就学习而言，值得去尝试一下，谁知道以后会不会用到呢。本文使用Intel Parallel Studio XE 2020和vs2019。

完整代码


program mainuse lapack95implicit nonecharacter*100 pathinteger num,point,i,j,counterreal,dimension(:),allocatable:: angle,length,uv,forces,uv_d,forces_dinteger,dimension(:),allocatable:: num1,num2,location,creal,dimension(:,:,:),allocatable:: mreal,dimension(:,:),allocatable:: mm,mm_d,a,wq,tempwrite(*,'("input file road:")')read*,pathopen(1,file=path)read(1,*)num,pointallocate(num1(num))allocate(num2(num))allocate(angle(num))allocate(length(num))allocate(m(num,4,4))allocate(uv(point*2))allocate(forces(point*2))allocate(mm(point*2,point*2))allocate(a(point*2,point*2))allocate(c(point*2))allocate(location(point*2))allocate(temp(point*2,2))!读入杆件位点，计算刚度矩阵do i=1,numread(1,*)num1(i),num2(i),angle(i),length(i)call mat(angle(i),length(i),i,num,m)enddo !读入结点的力和位移do i=1,pointread(1,*)uv(2*i-1),uv(2*i),forces(2*i-1),forces(2*i)enddo!关闭读入文件close(1)!扩阶矩阵forall(i=1:point*2,j=1:point*2)mm(i,j)=0do i=1,numcall rebig(i,m,num,num1(i),num2(i),point,mm)   enddo!获取位移向量的非0元素及其位置counter=count(uv/=0)allocate(mm_d(point*2-counter,point*2-counter))allocate(uv_d(point*2-counter))allocate(forces_d(point*2-counter))allocate(wq(point*2-counter,2))forall(i=1:2*point)location(i)=0call get_location(uv,location,counter,point)!获取缩减后的位移和力call less(forces,forces_d,location,point,counter)call less(uv,uv_d,location,point,counter)!获取缩减后的刚度矩阵call mless(mm,mm_d,location,point,counter)!刚度矩阵求逆a=mm_dcall getrf(a,c)call getri(a,c)!计算得位移矩阵forall(i=1:point*2-counter,j=1:1) wq(i,j)=forces_d(i)wq=matmul(a,wq)forall(i=1:point*2-counter,j=1:1) uv_d(i)=wq(i,j)do i=1,point*2-counteruv(location(point*2-counter-i+1))=uv_d(i)enddo!计算得力矩阵forall(i=1:point*2,j=1:1) temp(i,j)=uv(i)temp=matmul(mm,temp)forall(i=1:point*2,j=1:1) forces(i)=temp(i,j)!输出到结果文档open(2,file='result.txt')do i=1,pointwrite(2,100)i,uv(2*i-1),uv(2*i)enddodo i=1,pointwrite(2,200)i,forces(2*i-1),forces(2*i)enddo100 format(1x,'结点',i1,'的位移：',f6.2,',',f6.2)
200 format(1x,'结点',i1,'的力：',f6.2,',',f6.2)   close(2)print*,'done'pausecontains
subroutine mat(a,l,s,n,ma)integer s,n,i,jreal lamda,miu,a,l,ma(n,4,4)lamda=cosd(a)miu=sind(a)ma(s,1,:)=(/lamda**2,miu*lamda,-lamda**2,-lamda*miu/)ma(s,2,:)=(/miu*lamda,miu**2,-lamda*miu,-miu**2/)ma(s,3,:)=(/-lamda**2,-miu*lamda,lamda**2,lamda*miu/)ma(s,4,:)=(/-lamda*miu,-miu**2,lamda*miu,miu**2/)forall(i=1:4,j=1:4) ma(s,i,j)=ma(s,i,j)/l
end subroutine matsubroutine rebig(s,ma,n,n1,n2,p,mass)integer s,n,n1,n2,p,i,jreal ma(n,4,4)real mas(p*2,p*2),mass(p*2,p*2)forall(i=1:p*2,j=1:p*2)mas(i,j)=0mas(n1*2-1,n1*2-1)=ma(s,1,1);mas(n1*2-1,n1*2)=ma(s,1,2)mas(n1*2-1,n2*2-1)=ma(s,1,3);mas(n1*2-1,n2*2)=ma(s,1,4)mas(n1*2,n1*2-1)=ma(s,2,1);mas(n1*2,n1*2)=ma(s,2,2)mas(n1*2,n2*2-1)=ma(s,2,3);mas(n1*2,n2*2)=ma(s,2,4)mas(n2*2-1,n1*2-1)=ma(s,3,1);mas(n2*2-1,n1*2)=ma(s,3,2)mas(n2*2-1,n2*2-1)=ma(s,3,3);mas(n2*2-1,n2*2)=ma(s,3,4)mas(n2*2,n1*2-1)=ma(s,4,1);mas(n2*2,n1*2)=ma(s,4,2)mas(n2*2,n2*2-1)=ma(s,4,3);mas(n2*2,n2*2)=ma(s,4,4)forall(i=1:p*2,j=1:p*2)mass(i,j)=mass(i,j)+mas(i,j)
end subroutine rebigsubroutine get_location(uv_,location_,number,p)integer p,number,i,location_(2*p),nreal uv_(2*p)n=numberdo while(n>0)do i=1,p*2if(uv_(i)/=0)thenlocation(n)=in=n-1endifenddoenddo
end subroutine get_locationsubroutine less(old,new_,loca,p,counter_) integer counter_,loca(counter_),p,i,creal new_(p*2-counter_),old(2*point)c=counter_i=1do while(c>0)new_(i)=old(loca(c))c=c-1i=i+1enddo
end subroutine lesssubroutine mless(mm_,mm_d_,location_,p,count)integer count,location_(count),p,i,j,q,k,s,tfreal mm_d_t(p*2-count,p*2-count),mm_d_(p*2-count,p*2-count),mm_(2*p,2*p),temp((p*2-count)**2)q=1tf=0do i=1,2*pdo j=1,2*pdo k=1,2*p-countdo s=1,2*p-countif (i==location_(k).and.j==location_(s))thentf=1endifenddoenddoif(tf==1)thentemp(q)=mm_(i,j)q=q+1endiftf=0enddoenddomm_d_t=reshape(temp,(/p*2-count,p*2-count/))forall(i=1:p*2-count,j=1:p*2-count)mm_d_(i,j)=mm_d_t(j,i)
end subroutine mlessend

基于Fortran的结构力学位移法编程求解相关推荐

基于Python的结构力学位移法编程求解
如果想用Fortran编,可以参考基于Fortran的结构力学位移法编程求解目录 1.背景介绍 2.编程思路 2.1 结点编号.元素编号,建立总体坐标系xoy 2.2 建立各元素在总体坐标系xoy ...
matlab 矩阵位移法编程结构力学,matlab 矩阵位移法编程结构力学
矩阵位移法编程大作业 (091210211) 一.编制原理本程序的原理是基于结构力学矩阵位移法原理,以结构结点位移作基本未知量,将要分析的结构拆成已知节点力-结点力位移关系的单跨梁集合,通过强令结构 ...
matlab 矩阵位移法编程结构力学,matlab 矩阵位移法编程结构力学.doc
matlab 矩阵位移法编程结构力学.doc 矩阵位移法编程大作业(091210211)一.编制原理本程序的原理是基于结构力学矩阵位移法原理,以结构结点位移作基本未知量,将要分析的结构拆成已知节点力 ...
C语言练习题：三色球分组,编程计算三色球问题。若一个口袋中放有12个球，其中有3个红色的，3个白色的，6个黑色的，从中任取8个球，问共有多少种不同的颜色搭配？请用穷举法编程求解。
编程计算三色球问题.若一个口袋中放有12个球,其中有3个红色的,3个白色的,6个黑色的,从中任取8个球,问共有多少种不同的颜色搭配?请用穷举法编程求解. **输入格式:无 **输出格式:"i ...
矩阵位移法是用于求解杆系结构的计算机方法,结构力学的教学思路
结构力学的教学思路 2019-10-14 版权声明举报文章一.经典位移法(以下简称位移法)和矩阵位移法都是求解杆系结构的基本方法,是结构力学课程中两个十分重要的内容,两种方法都是结构力学课程中讲授 ...
C语言编程实现求解应力向量,基于MATLAB与Visual C 实现主应力的编程求解.doc
您所在位置:网站首页 > 海量文档 &nbsp>&nbsp计算机&nbsp>&nbspmatlab 基于MATLAB与Visual C 实现主应力的编 ...
矩阵位移法matlab编程,矩阵位移法_MATLAB_GUI.doc
Matrix_Displacement_Method--by MATLAB GUI PAGE58 / NUMPAGES64 yanfeng39@zju.edu.cn <结构力学>课程设计之 ...
用gulp计算合金弹性常数_考研结构力学用3种方法解东南大学位移法大题
位移法中含无穷刚度杆的结构是考研结构力学的一大难点,很多热门院校都喜欢出这类型的题目,下面以两道有复杂牵连位移的含无穷刚度杆位移法题目为例,对三种解法进行讲解,题目取自东南大学真题和烟台大学真题. 1 ...
结构力学计算机矩阵位移法,结构力学-矩阵位移法.ppt
<结构力学-矩阵位移法.ppt>由会员分享,可在线阅读,更多相关<结构力学-矩阵位移法.ppt(42页珍藏版)>请在装配图网上搜索. 1.结构力学,2 / 42,第七节平面刚 ...

基于Fortran的结构力学位移法编程求解

基于Fortran的结构力学位移法编程求解相关推荐

最新文章

热门文章