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

al2=i-1

al1=m-al2-al3

am=m

al1=al1/am

al2=al2/am

al3=al3/am

an(1)=(2.*al1-1.)*al1

an(2)=(2.*al2-1.)*al2

an(3)=(2.*al3-1.)*al3

an(4)=4.*al1*al2

an(5)=4.*al2*al3

an(6)=4.*al3*al1

ax=0.

ay=0.

az=0.

do 35 k=1,6

ax=ax+an(k)*xx(k)

ay=ay+an(k)*yy(k)

az=az+an(k)*zz(k)

35 continue

x(nod)=ax

y(nod)=ay

z(nod)=az

node(nod)=nod

nod=nod+inc

40 continue

go to 199

999 write(6,*)' error,triangulat surface generation'

write(6,*)' will be terminated because: either number'

write(6,*)' of divisions or node increment are wrong'

stop

c now generate element connectivity

199 if(itype.ne.1)go to 300

c triangular element with 3 nodes

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

n1t=nco(1)

n4t=n1t+(nco(2)-nco(1))/m

if(n4t.eq.n1t)go to 999

do 210 i=1,m

k1=m-i+1

k2=k1-1

n1=n1t

do 200 j=1,k1

ne=ne+1

n2=n1+1

n3=n2+k1

icon(ne,1)=ne

icon(ne,2)=n1

icon(ne,3)=n2

icon(ne,4)=n3

n1=n1+1

200 continue

n1t=n1t+k1+1

if(k2.eq.0)go to 210

n4=n4t

do 230 j=1,k2

ne=ne+1

n5=n4+k1+1

n6=n5-1

icon(ne,1)=ne

icon(ne,2)=n4

icon(ne,3)=n5

icon(ne,4)=n6

n4=n4+1

230 continue

n4t=n4t+k1+1

210 continue

return

300 continue

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

c triangular element with 6 nodes

mt=(m/2)*2

if(mt.ne.m)write(6,*)' error, for 6 nodetriangular element

*the number of division must be even'

n1t=nco(1)

n7t=n1t+2

mt=m/2

ne=0

do 411 i=1,mt

inc=m-2*i+2

k1=mt-i+1

k2=k1-1

n1=n1t

do 412 j=1,k1

ne=ne+1

n2=n1+2

n3=n2+2*inc-1

n4=n1+1

n5=n2+inc

n6=n1+inc+1

icon(ne,1)=ne

icon(ne,2)=n1

icon(ne,3)=n2

icon(ne,4)=n3

icon(ne,5)=n4

icon(ne,6)=n5

icon(ne,7)=n6

n1=n1+2

412 continue

if(k2.eq.0)go to 411

n1t=n1t+2*inc+1

n7=n7t

do 413 j=1,k2

ne=ne+1

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

n8=n7+2*inc+1

n9=n8-2

n10=n7+inc+1

n11=n10+inc-1

n12=n10-1

icon(ne,1)=ne

icon(ne,2)=n7

icon(ne,3)=n8

icon(ne,4)=n9

icon(ne,5)=n10

icon(ne,6)=n11

icon(ne,7)=n12

n7=n7+2

413 continue

n7t=n7t+2*inc+1

411 continue

return

end

subroutine OUTPUT(limit,nstor,x,y,z,icon,iold,new,nodold)

c write out the nodal point coordinate and its element connectivity

c to tape 7

dimension nstor(limit),x(limit),y(limit),z(limit),icon(limit,8)

*,iold(limit),new(limit),nodold(limit)

c merge algorithm

c first merge the coincide nodes and renumber the mesh

c find coincident nodes

ico=0

zone=0.001

itotn=0

ll=limit-1

if(ll.le.0) return

do 21 i=1,ll

if(nstor(i).eq.0)go to 20

j=i+1

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

do 19 k=j,limit

if(nstor(k).eq.0) go to 18

a=abs(x(i)-x(k))

b=abs(y(i)-y(k))

c=abs(z(i)-z(k))

if(a.le.zone.and.b.le.zone.and.c.le.zone)go to 22

go to 18

c i,k are coincide

22 ico=ico+1

itotn=itotn+1

new(ico)=nstor(i)

iold(ico)=nstor(k)

nstor(k)=-1

18 continue

19 continue

20 continue

21 continue

if(ico.ne.0)go to 30

write(6,*)' ++w++there were no coincident nodes with zone=',zone

go to 500

c change the node numbering sequence

30 isub=0

ki=1

do 40 i=1,limit

if(nstor(i).ne.-1)go to 35

c it is a coincident node

isub=isub+1

nodold(i)=0

nstor(i)=0

go to 40

35 if(isub.eq.0)go to 36

c change the old number to new number

ico=ico+1

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

iold(ico)=nstor(i)

new(ico)=nstor(i)-isub

36 nstor(ki)=nstor(i)-isub

c clear the higher node number

if(ki.ne.i)nstor(i)=0

mp=1

if(nodold(i).lt.0)mp=-1

x(ki)=x(i)

y(ki)=y(i)

z(ki)=z(i)

ki=ki+1

40 continue

c change element connection for new node numbering

do 100 i=1,limit

if(icon(i,1).eq.0)go to 100

do 60 j=2,8

if(icon(i,j).eq.0)go to 60

do 70 k=1,ico

if(icon(i,j).ne.iold(k))go to 70

icon(i,j)=new(k)

70 continue

60 continue

100 continue

write(6,*)' coincident nodes deleted',itotn

c write to tape 7

500 write(7,*)' nodal point coordinates and element connections'

itot=0

do 510 i=1,limit

if(nstor(i).le.0)go to 510

itot=itot+1

510 continue

c write total number of nodes

write(7,511)itot

511 format(5x,i5)