您所在位置：主题：有个程序 写的没问题 可是结果不对 请大家帮忙看看 怎么回事
程序为
program main
implicit none
integer iek,eno,i,j,m
real xi,yi,xj,yj,xm,ym
real E,NU,t,p
real betai,betaj,betam
real gammai,gammaj,gammam
real B,D,k,area
c real, allocatable,dimension(:):: k,area,B,D
open(unit=11,file='setup.txt',status='old')
read(11,*)
read(11,*)t
read(11,*)
read(11,*)E
read(11,*)
read(11,*)NU
read(11,*)
read(11,*)eno
write(*,*)t,E,NU,eno
call sspp(E,NU,t,xi,yi,xj,yj,xm,ym,
& p,area,betai,betaj,betam,gammai,gammaj,gammam,
& B,D,k,iek,eno)
close(11)
end
c------------------------------------------------------------------
subroutine sspp(E,NU,t,xi,yi,xj,yj,xm,ym,
& p,area, betai,betaj,betam,gammai,gammaj,gammam,
& B,D,k,iek,eno)
implicit none
integer iek,eno,lll
real xi,yi,xj,yj,xm,ym
real area(iek)
real E,NU,t,p
real betai,betaj,betam
real gammai,gammaj,gammam
real B(iek,3,6),D(iek,3,3),k(iek,6,6),BT(iek,6,3)
real c1(6,3),c2(6,6)
open(unit=11,file='coordinates.nas',status='old')
read(11,*)
c do lll=1,1000
read(11,*)iek,xi,yi,xj,yj,xm,ym;
write(*,*)iek,xi,yi,xj,yj,xm,ym;
write(*,*)'area(',iek,')=',area(iek);
call LinearTriangleElementStiffness(E,NU,t,xi,yi,xj,yj,xm,ym,
& p,area, betai,betaj,betam,gammai,gammaj,gammam,
& B,D,k,iek,eno)
c101 end do
close(11)
end
c----------------------------------------------------------------------e
c LinearTriangleElementStiffness This function returns the element
c stiffness matrix for a linear
c triangular element with modulus of
c elasticity E,Poisson's ratio NU,
c thickness t,coordinates of the first
c node(xi,yi),coordinates of the second
c node(xj,yj) and coordinates of the
c third node(xm,ym). Use p=1 for case of
c plane stress, and p=2 for cases for
c plane strain.
c the size of the element stiffness
c matrix is 6*6.
subroutine LinearTriangleElementStiffness(E,NU,t,xi,yi,xj,yj,xm,
& ym,p,area,betai,betaj,betam,gammai,gammaj,gammam,
& B,D,k,iek,eno)
implicit none
integer iek,eno
real E,NU,t
real xi,yi,xj,yj,xm,ym,p
real betai,betaj,betam
real gammai,gammaj,gammam
integer i,j,m
real B(iek,3,6),D(iek,3,3),k(iek,6,6),BT(iek,6,3),
& area(iek),c1(iek,6,3),c2(iek,6,6)
B=0;
D=0;
k=0;
BT=0;
c1=0;
c2=0;
area(iek)=(xi*(yj-ym)+xj*(ym-yi)+xm*(yi-yj))/2;
write(*,*)'area(',iek,')=',area(iek);
betai=yj-ym;
betaj=ym-yi;
betam=yi-yj;
gammai=xm-xj;
gammaj=xi-xm;
gammam=xj-xi;
B(iek,1,1)=betai/(2*area(iek));
B(iek,1,2)=0/(2*area(iek));
B(iek,1,3)=betaj/(2*area(iek));
B(iek,1,4)=0/(2*area(iek));
B(iek,1,5)=betam/(2*area(iek));
B(iek,1,6)=0/(2*area(iek));
B(iek,2,1)=0/(2*area(iek));
B(iek,2,2)=gammai/(2*area(iek));
B(iek,2,3)=0/(2*area(iek));
B(iek,2,4)=gammaj/(2*area(iek));
B(iek,2,5)=0/(2*area(iek));
B(iek,2,6)=gammam/(2*area(iek));
B(iek,3,1)=gammai/(2*area(iek));
B(iek,3,2)=betai/(2*area(iek));
B(iek,3,3)=gammaj/(2*area(iek));
B(iek,3,4)=betaj/(2*area(iek));
B(iek,3,5)=gammam/(2*area(iek));
B(iek,3,6)=betam/(2*area(iek));
write(*,*)'B(',iek,')=',B(iek,1,:);
write(*,*)B(iek,2,:);
write(*,*)B(iek,3,:);
c---------------------has:B(iek,3,6)---------------------------------
do i=1,3
do j=1,6
BT(iek,j,i)=B(iek,i,j)
enddo
enddo
write(*,*)'BT(',iek,')=',BT(iek,1,:);
write(*,*)BT(iek,2,:);
write(*,*)BT(iek,3,:);
write(*,*)BT(iek,4,:);
write(*,*)BT(iek,5,:);
write(*,*)BT(iek,6,:);
c-----------------has:BT(iek,6,3)---------------------------------
open(unit=11,file='p_number.txt',status='old')
read(11,*)
read(11,*)p
IF (p==1) THEN
D(iek,1,1)=(E/(1-NU**2))*1;
D(iek,1,2)=(E/(1-NU**2))*NU;
D(iek,1,3)=(E/(1-NU**2))*0;
D(iek,2,1)=(E/(1-NU**2))*NU;
D(iek,2,2)=(E/(1-NU**2))*1;
D(iek,2,3)=(E/(1-NU**2))*0;
D(iek,3,1)=(E/(1-NU**2))*0;
D(iek,3,2)=(E/(1-NU**2))*0;
D(iek,3,3)=(E/(1-NU**2))*((1-NU)/2);
ELSEIF (p==2) THEN
D(iek,1,1)=(E/(1+NU)/(1-2*NU))*(1-NU);
D(iek,1,2)=(E/(1+NU)/(1-2*NU))*NU;
D(iek,1,3)=(E/(1+NU)/(1-2*NU))*0;
D(iek,2,1)=(E/(1+NU)/(1-2*NU))*NU;
D(iek,2,2)=(E/(1+NU)/(1-2*NU))*(1-NU);
D(iek,2,3)=(E/(1+NU)/(1-2*NU))*0;
D(iek,3,1)=(E/(1+NU)/(1-2*NU))*0;
D(iek,3,2)=(E/(1+NU)/(1-2*NU))*0;
D(iek,3,3)=(E/(1+NU)/(1-2*NU))*((1-2*NU)/2);
ENDIF
write(*,*)'p=',p;
write(*,*)'D(',iek,')=',D(iek,1,:);
write(*,*)D(iek,2,:);
write(*,*)D(iek,3,:);
c------------------has:D(iek,3,3)----------------------------------
c1=0
do i=1,6
do j=1,3
c1(iek,i,j)=0
do m=1,3
c1(iek,i,j)=c1(iek,i,j)+BT(iek,i,m)*D(iek,m,j)
enddo
enddo
enddo
write(*,*)'c1(',iek,')=',c1(iek,1,:);
write(*,*)c1(iek,2,:);
write(*,*)c1(iek,3,:);
write(*,*)c1(iek,4,:);
write(*,*)c1(iek,5,:);
write(*,*)c1(iek,6,:);
c------------------has:c1(iek,6,3)----------------------------------------
write(*,*)'c2(',iek,')=',c2(iek,1,:);
write(*,*)'c1(',iek,')=',c1(iek,1,:);
write(*,*)'B(',iek,')=',B(iek,1,:);
do i=1,6
do j=1,6
c2(iek,i,j)=0
do m=1,3
c2(iek,i,j)=c2(iek,i,j)+c1(iek,i,m)*B(iek,m,j)
enddo
end do
end do
write(*,*)'c2(',iek,')=',c2(iek,1,:);
write(*,*)'c1(',iek,')=',c1(iek,1,:);
write(*,*)'B(',iek,')=',B(iek,1,:);
c write(*,*)c2(iek,2,:);
c write(*,*)c2(iek,3,:);
c write(*,*)c2(iek,4,:);
c write(*,*)c2(iek,5,:);
c write(*,*)c2(iek,6,:);
c-------------------has:c2(iek,6,6)--------------------------------------
k=0
do i=1,6
do j=1,6
k(iek,i,j)=0
k(iek,i,j)=k(iek,i,j)+t*area(iek)*c2(iek,i,j)
enddo
enddo
write(*,*)'k(',iek,')=',k(iek,1,:);
write(*,*)k(iek,2,:);
write(*,*)k(iek,3,:);
write(*,*)k(iek,4,:);
write(*,*)k(iek,5,:);
write(*,*)k(iek,6,:);
write(*,*)'------------------------------------'
c-------------------has:k(iek,6,6)--------------------------------------
close(11)
endsubroutine
可是结果出现到 数组k 有值，且能输出，但是从数组c2始就不对了，错的离谱，结果还有下面的语句
forrtl: severe (174): SIGSEGV, segmentation fault occurred
Image PC Routine Line Source
a.out 00000000004031A9 Unknown Unknown Unknown
Stack trace terminated abnormally.
不知道问题在哪里 ？