标题： 一个简单的三角形网格生成程序，供初学者参考
作者:xurui [3月28日,8:32]
这是一个简单的采用advanced front方法编制的三角形网格生成程序，供初学者参考。我看周华的软件下载区一直空着，就抛砖引玉一下，希望有更多的人能贡献程序共享。若周华看到后，希望能把这个小小的程序放到下载区，因为在这里很快会被淹没。
下面是源码
program bili
open(10,file='aa.map')
open(3,file='aa')
call coor0
close(10)
end
subroutine coor0
implicit real*4 (a-h,o-z)
common/nxy/nx,ny,numnod,nume,numb,m
common /enod/ nod1(80000),nod2(80000),nod3(80000)
common /coor/ x(80000),y(80000)
integer n,nod1,nod2,nod3
c write(*,*) 'input ≠&nx;&ny; '
open(1,file='aa.dat')
read(1,*) ne
do i=1,ne
read(1,*)x(i),y(i)
enddo
ns = 1
c ne = (nx+ny)*2
nume = 0
nums = 0
call tricb(nume,ns,ne)
numnod = nume
call trimesh(nums,numnod,nume,5,x,y)
call trinod(nume,nt,x,y,nod1,nod2,nod3)
mtij=numnod
mtik= nt
write(10,*) mtij,mtik
write(10,'(2f8.2)')(x(mtii),y(mtii),mtii=1,mtij)
write(10,'(3i8)')(nod1(mtii),nod2(mtii),nod3(mtii),mtii=1,mtik)
return
end
c..................................................................
c the following subroutines for generating triangular mesh
c according to circle boundary nodal points
c trimesh generates nodal points of mesh
c trinod generates nodal No. of each triangle
c tricb generates boundary edges according to
c continuous nodal No. of circle boundary
c tricbn generates boundary edges according to
c discontinuous nodal No. of circle boundary
c..................................................................
c trimesh(nums,numnod,numedge,mc,x,y)
c trinod(numedge,nt,x,y,nod1,nod2,nod3)
c tricb(numedge,ns,ne)
c tricbn(numedge,numnod,node)
c..................................................................
c nums ....... starting edge No. - 1
c numnod ..... number of nodes
c numedge .... number of edges
c mc ......... number of iteration for modifying nodal points
c nt ......... number of triangles
c ns ......... first nodal point
c ne ......... last nodal point
c..................................................................
c x,y ............... coordinate arrays of nodal points
c nod1,nod2,nod3 .... nodal No. arrays for each triangle
c node .............. circle boundary nodal No. array
c..................................................................
subroutine trimesh(nums,numnod,nume,mc,x,y)
c .... generate new nodal points for triagle elements
c .... store all nodal No. of edges in arrays nod1 & nod2
implicit real*4 (a-h,o-z)
common /edge/ nod1(80000),nod2(80000)
dimension x(*),y(*),nv(100)
COMMON /TMDATA/ index(80000)
integer edge1,edge2
numb = numnod
c write(*,*) 'numb =',numb,' x,y ='
c write(*,'(1x,6f12.4)') (x(n),n=1,numnod)
c write(*,'(1x,6f12.4)') (y(n),n=1,nu
mnod)
c write(*,*) 'nume =',nume,' nod1,nod2 ='
c write(*,'(1x,15i5)') (nod1(n),n=1,nume)
c write(*,'(1x,15i5)') (nod2(n),n=1,nume)
dmin = 1.e20
dmax = 0.0
do 50 n=1,nume
index(n) = 0
i = nod1(n)
j = nod2(n)
d = (x(i)-x(j))**2+(y(i)-y(j))**2
if (d.lt.dmin) dmin = d
if (d.gt.dmax) dmax = d
50 continue
dmin = sqrt(dmin)
dmax = sqrt(dmax)
c write(*,*) 'dmin,dmax =',dmin,dmax
c write(*,*) 'index ='
write(*,'(1x,15i5)') (index(n),n=1,nume)
c .... m is the middle point of n1 & n2
c .... m_h is the normal vector of n1_n2
c .... h is the relative high with the bottom n1_n2
cccc nums = 0
h = sqrt(3.0)
cccc h = 1.d0
100 continue
nums = nums+1
write(*,*) 'nums, index =',nums,index(nums)
if (index(nums).gt.0) goto 2400
write(*,*) 'nums ,nume =',nums,nume
n1 = nod1(nums)
n2 = nod2(nums)
write(*,*) 'n1,n2 =',n1,n2
xx1 = x(n1)
yy1 = y(n1)
xx2 = x(n2)
yy2 = y(n2)
xm = (xx1+xx2)/2
ym = (yy1+yy2)/2
a = xm-xx1
b = ym-yy1
dd = (a*a+b*b)*4
d = sqrt(dd)
hh = h
if (d.lt.dmin) hh = h*dmin/d
if (d.gt.dmax) hh = h*dmax/d
xh = xm - b*hh
yh = ym + a*hh
write(*,*) 'xh,yh =',xh,yh
e = d*hh*1.0e-4
ee = e*1.e-30
nn = 0
smin = dmax*dmax
do 140 i=1,numnod
cccc do 140 n=nums+1,nume
cccc i = nod2(n)
xx = x(i)
yy = y(i)
s = (xx-xh)**2+(yy-yh)**2
if (s.lt.smin) then
j = i
smin = s
endif
140 continue
smin = sqrt(smin)
write(*,*) 'j,smin =',j,smin
cccc if (smin.lt.dmin*0.5) then
if (smin.lt.dmin*0.5 .or. smin.lt.d*0.5) then
nn = j
write(*,*) 'nn,smin,d =',nn,smin,d
xh = x(nn)
yh = y(nn)
write(*,*) 'xh,yh =',xh,yh
index(nums) = 2
endif
cccc do 1000 n=nums+1,nume
do 1000 n=1,nume
c write(*,*) 'n =',n
if (n.eq.nums) goto 1000
k1 = nod1(n)
k2 = nod2(n)
if (k1.eq.nn .or. k2.eq.nn) goto 1000
c write(*,*) 'k1,k2 =',k1,k2
x1 = x(k1)
x2 = x(k2)
y1 = y(k1)
y2 = y(k2)
dm = det(x1,x2,xm,y1,y2,ym)
dh = det(x1,x2,xh,y1,y2,yh)
if (dm*dh.gt.-ee) goto 1000
d1 = det(xm,xh,x1,ym,yh,y1)
d2 = det(xm,xh,x2,ym,yh,y2)
if (d1*d2.gt.-ee) goto 1000
write(*,*) 'dm,dh =',dm,dh
write(*,*) 'd1,d2 =',d1,d2
write(*,*) 'n,k1,k2 =',n,k1,k2
write(*,*) 'x1,x2,y1,y2 =',x1,x2,y1,y2
if (k1.eq.n1 .or. k1.eq.n2) then
nn = k2
goto 500
endif
if (k2.eq.n1 .or. k2.eq.n2) then
nn = k1
goto 500
endif
if (d1*d1.lt.d2*d2) then
nn = k1
else
nn = k2
endif
500 continue
xh = x(nn)
yh = y(nn)
index(nums) = 4
1000 continue
do 400 n=1,numnod
if (n.eq.n1 .or. n.eq.n2 .or. n.eq.nn) goto 400
xn = x(n)
yn = y(n)
if (det(xn,xx1,xx2,yn,yy1,yy2).lt.-e) goto 400
if (det(xn,xx2,xh,yn,yy2,yh).lt.-e) goto 400
if (det(xn,xh,xx1,yn,yh,yy1).lt.-e) goto 400
nn = n
xh = xn
yh = yn
index(nums) = 3
write(*,*) 'n1,nn,n2,e =',n1,nn,n2,e
write(*,*) 'nn,xh,yh ============',nn,xh,yh
400 continue
write(*,*) 'nn,nums,index =',nn,nums,index(nums)
if (nn.ne.0) goto 2100
index(nums) = 1
numnod = nu
mnod+1
x(numnod) = xh
y(numnod) = yh
write(*,*) 'numnod,x,y =',numnod,x(numnod),y(numnod)
nume = nume+1
index(nume) = 0
nod1(nume) = nod1(nums)
nod2(nume) = numnod
nod1(nume+1) = numnod
nod2(nume+1) = nod2(nums)
write(*,*) 'numnod,nums,nume,nod1,nod2 =',numnod,nums,nume,
1 nod1(nume),nod2(nume),nod1(nume+1),nod2(nume+1)
write(3,*) 'numnod,nums,nume,nod1,nod2 =',numnod,nums,nume,
1 nod1(nume),nod2(nume),nod1(nume+1),nod2(nume+1)
nume = nume+1
index(nume) = 0
cccc nums = nums+1
goto 2400
2100 continue
edge1 = 1
edge2 = 1
write(*,*) 'n1,nn,n2 =',n1,nn,n2
cccc do 2200 n=nums+1,nume
do 2200 n=1,nume
write(*,*) 'n,nod1,nod2 =',n,nod1(n),nod2(n)
if ((n1.eq.nod1(n) .and. nn.eq.nod2(n)) .or.
& (n1.eq.nod2(n) .and. nn.eq.nod1(n))) then
index(n) = 5
write(*,*) 'n ..........',n
edge1 = 0
endif
if ((n2.eq.nod1(n) .and. nn.eq.nod2(n)) .or.
& (n2.eq.nod2(n) .and. nn.eq.nod1(n))) then
index(n) = 5
write(*,*) 'n ..........',n
edge2 = 0
endif
write(*,*) 'edge1,edge2 =',edge1,edge2
2200 continue
write(*,*) 'edge1,edge2 =',edge1,edge2
if (edge1.eq.1) then
nume = nume+1
index(nume) = 0
nod1(nume) = n1
nod2(nume) = nn
endif
if (edge2.eq.1) then
nume = nume+1
index(nume) = 0
nod1(nume) = nn
nod2(nume) = n2
endif
2400 continue
if (nume.gt.nums) goto 100
call delnod(numb,numnod,nume,nod1,nod2,x,y,nv)
c mc = 20
do 2800 k=1,mc
do 2700 n=numb+1,numnod
m = 0
do 2600 i=1,nume
if (n.eq.nod1(i)) then
m = m+1
nv(m) = nod2(i)
goto 2600
endif
if (n.eq.nod2(i)) then
m = m+1
nv(m) = nod1(i)
endif
2600 continue
xc = 0.0
yc = 0.0
do 2650 i=1,m
j = nv(i)
xc = xc+x(j)
yc = yc+y(j)
2650 continue
x(n) = xc/m
y(n) = yc/m
2700 continue
2800 continue
write(*,*) 'numnod =',numnod
write(*,'(1x,6f12.4)') (x(n),n=1,numnod)
write(*,'(1x,6f12.4)') (y(n),n=1,numnod)
write(*,*) 'nume =',nume,' nod1,nod2 ='
write(*,'(1x,15i5)') (nod1(n),n=1,nume)
write(*,'(1x,15i5)') (nod2(n),n=1,nume)
write(*,*) ' n nod1 nod2 index '
do 3000 n=1,nume
write(*,'(1x,4i5)') n,nod1(n),nod2(n),index(n)
3000 continue
return
end
subroutine delnod(numb,numnod,nume,nod1,nod2,x,y,nv)
implicit real*4 (a-h,o-z)
dimension nod1(nume),nod2(nume),x(numnod),y(numnod),nv(*)
250 continue
do 270 n=numb+1,numnod
m = 0
do 260 i=1,nume
if (n.eq.nod1(i)) then
m = m+1
nv(m) = i
goto 260
endif
if (n.eq.nod2(i)) then
m = m+1
nv(m) = i
endif
if (m.gt.3) goto 270
260 continue
if (m.le.3) goto 300
270 continue
return
300 continue
do 310 i=n+1,numnod
x(i) = x(i+1)
y(i) = y(i+1)
310 continue
nn = 0
do 320 i=1,nume
do 315 j=1,m
if (i.eq.nv(j)) goto 320
315 continue
nn = nn+1
nod1(nn) = nod1(i)
nod2(nn) = nod2(i)
if (nod1(nn).gt.n) nod1(nn) = nod1(nn)-1
if (nod2(nn).gt.n) nod2(nn) = nod2(nn)-1
320 continue
numnod = numnod-1
nume = nn
goto 250
end
subroutine trinod(nume,
nt,x,y,nod1,nod2,nod3)
c .... generate nodal No. of each triangle
c .... according to nodal No. of edges
implicit real*4 (a-h,o-z)
common /edge/ n1(80000),n2(80000)
dimension nod1(*),nod2(*),nod3(*)
dimension x(*),y(*)
write(*,*) 'nume =',nume,' n1,n2 ='
write(*,'(1x,15i5)') (n1(n),n=1,nume)
write(*,'(1x,15i5)') (n2(n),n=1,nume)
nt = 0
do 500 n=1,nume
i = n1(n)
j = n2(n)
do 300 m=n+1,nume
m1 = n1(m)
m2 = n2(m)
if (m1.eq.i .or. m2.eq.i) then
if (m1.eq.i) then
m12 = m2
else
m12 = m1
endif
do 100 k=n+1,nume
k1 = n1(k)
k2 = n2(k)
if ((k1.eq.j .and. k2.eq.m12) .or.
& (k2.eq.j .and. k1.eq.m12)) then
x1 = x(i)
x2 = x(j)
x3 = x(m12)
y1 = y(i)
y2 = y(j)
y3 = y(m12)
d = det(x1,x2,x3,y1,y2,y3)
nt = nt+1
if (d.gt.0.0d0) then
nod1(nt) = i
nod2(nt) = j
else
nod1(nt) = j
nod2(nt) = i
endif
nod3(nt) = m12
endif
100 continue
endif
300 continue
500 continue
write(*,*) 'nt =',nt,' nod1,nod2,nod3 ='
do 1000 n=1,nt
write(*,'(1x,3i5)') nod1(n),nod2(n),nod3(n)
1000 continue
return
end
real*4 function det(x1,x2,x3,y1,y2,y3)
implicit real*4 (a-h,o-z)
det = x1*y2+x2*y3+x3*y1
& -y1*x2-y2*x3-y3*x1
c write(*,*) 'x1,x2,x3,y1,y2,y3 ='
c write(*,'(1x,3f12.4)') x1,x2,x3
c write(*,'(1x,3f12.4)') y1,y2,y3
c write(*,*) 'det =',det
return
end
subroutine tricb(n,ns,ne)
c .... generate circle boundary edges
c .... according to circle boundary nodes
c .... ns first nodal No.
c .... ne last nodal No.
c .... n starting edge No. - 1
common /edge/ nod1(80000),nod2(80000)
do 100 i=ns,ne
n = n+1
nod1(n) = i
nod2(n) = i+1
100 continue
nod2(n) = ns
return
end
以下是一个简单的输入文件
16
0.000000 0.000000
2.50000 0.000000
5.00000 0.000000
7.50000 0.000000
10.0000 0.000000
10.0000 2.50000
10.0000 5.00000
10.0000 7.50000
10.0000 10.0000
7.50000 10.0000
5.00000 10.0000
2.50000 10.0000
0.000000 10.0000
0.000000 7.50000
0.000000 5.00000
0.000000 2.50000