Продолжение приложения 1

if((diag1/diag2).gt.1.02) goto41

j1=nr(1)

j2=nr(ij+1)

j3=nr(ij+2)

goto40

41 j1=nr(ij)

j2=nr(ij+1)

j3=nr(4)

40 lb(1)=iabs(ne(j1)-ne(j2))+1

lb(2)=iabs(ne(j2)-ne(j3))+1

lb(3)=iabs(ne(j1)-ne(j3))+1

do 107 ik=1,3

if(lb(ik).le.nbw) goto 107

nbw=lb(ik)

nelbw=nel

107 continue

write(10,301) nel,ne(j1),ne(j2),ne(j3),xe(j1),ye(j1),xe(j2),ye(j2),xe(j3),ye(j3)

nop(nel,1)=ne(j1)

nop(nel,2)=ne(j2)

nop(nel,3)=ne(j3)

cord(ne(j1),1)=xe(j1)

cord(ne(j1),2)=ye(j1)

cord(ne(j2),1)=xe(j2)

cord(ne(j2),2)=ye(j2)

cord(ne(j3),1)=xe(j3)

cord(ne(j3),2)=ye(j3)

301 format(1x,4i5,3x,6f12.4)

if(ipch.eq.0) goto15

write(10,303) nel,ne(j1),ne(j2),ne(j3),xe(j1),ye(j1),xe(j2),ye(j2),xe(j3),ye(j3)

303 format(4i3,6f10.4)

Продолжение приложения 1

15 continue

16 continue

777 continue

endfile(11)

close(11)

write (10,51) nbw,nelbw

51 format(///1x,22h bandwidth quantity is,i4,22hcalculated in element ,i4)

write(10,*)' cord(500,3)'

do 998 i=1,500

998 write(10,*)' i=',i,'x=',cord(i,1),'y=',cord(i,2),'z=',cord(i,3)

write(10,*)' nop(500,3)'

do 997 i=1,500

997 write(10,*)' i=',i,'nop(i,1-3)=',nop(i,1),nop(i,2),nop(i,3)

DO 71 I=1,496

DO 71 J=1,3

71 NOP1(I,J)=NOP(I,J)

DO 72 I=1,290

DO 72 J=1,3

72 CRD(I,J)=CORD(I,J)

OPEN(31,FILE='NOP',FORM='UNFORMATTED')

WRITE(31)NOP1,CRD

CLOSE(31)

endfile(10)

close(10)

stop

end

Приложение 2

c main program mesh

common/total/nelem

dimension a(15000),nco(6),xx(6),yy(6),zz(6)

open(6,file='F6')

open(5,file='F5')

open(7,file='F7')

c set up dynamic allocation

limit=500

nelem=0

n1=1

n2=n1+limit

n3=n2+limit

n4=n3+limit

n5=n4+limit

n6=n5+limit*8

n7=n6+limit

n8=n7+limit

n9=n8+limit

if(n9.gt.15000)write(6,10)

10 format(' error .. increase the size of array a(20000)')

c read the element type and number of division on each side

11 write(6,*)' input element type(1,2), number of division'

itype=1

m=10

c terminate if itype is zero, no more generation

Продолжение приложения 2

if(itype.eq.0)go to 100

c read first 3 corner key points

c 3

c / \

c / \

c 1––––-2

write(6,*)' input key point node number and x,y,z coordinates'

do 20 i =1,6

read(5,*)nco(i),xx(i),yy(i),zz(i)

if(nco(i).lt.0)go to 30

20 continue

c generate the mesh

30 call surf(itype,m,nco,xx,yy,zz,limit,a(n1),a(n2),a(n3)

*,a(n4),a(n5))

go to 11

c write out the nodal point coordinates and element connections

c on tape number 7.

100 call output(limit,a(n1),a(n2),a(n3),a(n4),a(n5)

*,a(n6),a(n7),a(n8))

close(5)

close(6)

close(7)

stop

end

subroutine surf(itype,m,nco,xx,yy,zz,limit,node,x,y,z,icon)

c generate surface mesh

c itype=element type

c =1 for 3 node triangular element

c =2 for 6 node

c m =number of division on each side

c nco = array of 6 with key point node numbers

Продолжение приложения 2

c xx =arrary of 6 with global x coordinate of the key pointes

c yy =

c zz =

c limit=total number of nodes in the mesh

c node =nodal point numbers

c x,y,z= " " coordinates (global)

c icon = element connection

common/tutal/ne

dimension nco(6),xx(6),yy(6),zz(6),node(limit),x(limit)

*,y(limit),z(limit),icon(limit,8),an(6)

c first generate the coordinates by using area-coordinate system

ndi=m

print *' m=',m

if(ndi.eq.0)ndi=nco(2)-nco(1)

if(ndi.eq.0)go to 999

inc=(nco(2)-nco(1))/ndi

if(inc.le.0)go to 999

c check middle node

do 30 i=4,6

ii=i-2

if(ii.eq.4)ii=1

if(nco(i).ne.0)go to 30

xx(i)=(xx(i-3)+xx(ii))/2.

yy(i)=(yy(i-3)+yy(ii))/2.

zz(i)=(zz(i-3)+zz(ii))/2.

30 continue

m=ndi

nod=nco(1)

lll=m+1

do 40 j =1,lll

kk=m+2-j

do 40 i=1,kk

al3=j-1