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

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

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

program main use lapack95 implicit none character*100 path integer num,point,i,j,counter real,dimension(:),allocatable:: angle,length,uv,forces,uv_d,forces_d integer,dimension(:),allocatable:: num1,num2,location,c real,dimension(:,:,:),allocatable:: m real,dimension(:,:),allocatable:: mm,mm_d,a,wq,temp write(*,'("input file road:")') read*,path open(1,file=path) read(1,*)num,point allocate(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,num read(1,*)num1(i),num2(i),angle(i),length(i) call mat(angle(i),length(i),i,num,m) enddo !读入结点的力和位移 do i=1,point read(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)=0 do i=1,num call 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)=0 call 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_d call 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-counter uv(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,point write(2,100)i,uv(2*i-1),uv(2*i) enddo do i=1,point write(2,200)i,forces(2*i-1),forces(2*i) enddo 100 format(1x,'结点',i1,'的位移：',f6.2,',',f6.2) 200 format(1x,'结点',i1,'的力：',f6.2,',',f6.2) close(2) print*,'done' pause contains subroutine mat(a,l,s,n,ma) integer s,n,i,j real 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 mat subroutine rebig(s,ma,n,n1,n2,p,mass) integer s,n,n1,n2,p,i,j real 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)=0 mas(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 rebig subroutine get_location(uv_,location_,number,p) integer p,number,i,location_(2*p),n real uv_(2*p) n=number do while(n>0) do i=1,p*2 if(uv_(i)/=0)then location(n)=i n=n-1 endif enddo enddo end subroutine get_location subroutine less(old,new_,loca,p,counter_) integer counter_,loca(counter_),p,i,c real new_(p*2-counter_),old(2*point) c=counter_ do while(c>0) new_(i)=old(loca(c)) c=c-1 i=i+1 enddo end subroutine less subroutine mless(mm_,mm_d_,location_,p,count) integer count,location_(count),p,i,j,q,k,s,tf real 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) do i=1,2*p do j=1,2*p do k=1,2*p-count do s=1,2*p-count if (i==location_(k).and.j==location_(s))then endif enddo enddo if(tf==1)then temp(q)=mm_(i,j) q=q+1 endif enddo enddo mm_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 mless def circshift(matrix,shiftnum1,shiftnum2): h,w=matrix.shape matrix=np.vstack((matrix[(h-shiftnum1):,:],matrix[:(h-shiftnum1),:])) matrix=np.hstack((matrix[:,(w-shiftnum2):],matrix[:,:(w-shiftnum2)])) return matrix 上述代码中的matrix代