program sild3 * ------------------------------------------------------------ * * * * Generate the Laguere Delaunay diagram in the * * three-dimensional space * * * * Copyright in 1997 by Kokichi Sugihara * * (Revised on 2003-05-29) * * * * ------------------------------------------------------------ * common /ch40/ kpara(30) common /ch41a/ kx(1024),ky(1024),kz(1024),ku(1024) common /ch41b/ kvt(1024),kname(1024),mol(1024) common /ch41c/ kvact(1024) common /ch42/ kesv(20480),keev(20480),kecol(20480) common /ch43a/ kfv(3,20480),kft(2,20480),kfnt(20480) common /ch43b/ idface(20480) common /ch44/ ktv(4,10240),ktf(4,10240),ktcol(10240) common /ch4f1/ fltonv(333) common /ch4f2/ fltstk(10240),fltdel(10240) common /ch4f3/ flvert(1024) common /ch4w1/ istack(333),itlist(333),idlist(333) common /ch4w2/ islist(333),idflst(333) common /ch4w3/ ivlist(333) common /lague/ radius(1024) logical fltonv,fltstk,fltdel,flvert logical kvact * * --- set the basic parameters for the data structure --- * kvmax=1024 kpara(1)=kvmax kpara(2)=20*kvmax kpara(3)=20*kvmax kpara(4)=10*kvmax kpara(5)=333 kpara(6)=2**28-1 * * -------------------------------------------------------------- * kpara(1) maximum number of vertices * kpara(2) maximum number of edges * kpara(3) maximum number of faces * kpara(4) maximum number of tetrahedra * kpara(5) size of the stacks * kpara(6) maximum absolute value of the coordinate * kpara(7) number of the input vertices * -------------------------------------------------------------- * * --- read the coordinates of the spheres --- * umax=3.001*float(kpara(6)) open(2,file='spdata.dat') read(2,*) nofp kpara(7)=nofp do 10 i=1,nofp read(2,*) x,y,z,radius(i) kx(i)=ifix(x) ky(i)=ifix(y) kz(i)=ifix(z) ku(i)=ifix((x**2+y**2+z**2-radius(i)**2)/umax) kname(i)=i 10 continue * open(7,file='monitor') * call chull4 * call monit * open(1,file='draw.ps') call drlag3 call drdel3 close(1) * c call monit call statis * stop end * * subroutine drdel3 * ------------------------------------------------------------ * * * * draw the three-dimensional Laguere triangulation * * * * ------------------------------------------------------------ * common /ch40/ kpara(30) common /ch41a/ kx(1024),ky(1024),kz(1024),ku(1024) common /ch41b/ kvt(1024),kname(1024),mol(1024) common /ch42/ kesv(20480),keev(20480),kecol(20480) common /ch43a/ kfv(3,20480),kft(2,20480),kfnt(20480) common /ch43b/ idface(20480) common /ch44/ ktv(4,10240),ktf(4,10240),ktcol(10240) common /ch4f1/ fltonv(333) common /ch4f2/ fltstk(10240),fltdel(10240) common /ch4f3/ flvert(1024) common /ch4w1/ istack(333),itlist(333),idlist(333) common /ch4w2/ islist(333),idflst(333) common /ch4w3/ ivlist(333) logical fltonv,fltstk,fltdel,flvert * DATA WIDTH/0.5/,RADIUS/0.5/ DATA EYEX/1.0/,EYEY/0.0/,EYEZ/8.0/ DATA THETA/1.0/ c write(7,611) c 611 format('enter drdel3') ixymax=kpara(6) nofp=kpara(7) XYMAX=FLOAT(IXYMAX)*0.9 EX=XYMAX*EYEX EY=XYMAX*EYEY EZ=XYMAX*EYEZ * c open(1,file='draw.ps') * * --- classify the edges --- * call etcolo * * --- draw the front side in the above half and * the rear side in the below half --- * do 1 ifront=1,2 iup=ifront if(ifront.eq.1) then inhibi=2 else inhibi=1 endif * * --- draw the left diagram --- * CALL DRSTLE(-IXYMAX,IXYMAX,-IXYMAX,IXYMAX,iup) CALL LWIDTH(WIDTH) * * --- plot the input points --- * DO 10 II=1,NOFP X1=FLOAT(kX(II)) Y1=FLOAT(kY(II)) Z1=FLOAT(kZ(II)) CALL CENTPR(EX,EY,EZ,X1,Y1,Z1,XX,YY) IXX=IFIX(XX) IYY=IFIX(YY) CALL DRDOT(IXX,IYY,RADIUS) 10 CONTINUE * * --- draw the edges --- * DO 20 II=1,KPARA(22) * --- skip unnecessary edges --- if(kecol(ii).eq.inhibi) goto 20 * --- compute the start vertex --- IVS=KeSV(II) X1=kX(IVS) Y1=kY(IVS) Z1=kZ(IVS) CALL CENTPR(EX,EY,EZ,X1,Y1,Z1,XX,YY) IX1=IFIX(XX) IY1=IFIX(YY) * --- compute the end vertex --- IVE=KeEV(II) X2=kX(IVE) Y2=kY(IVE) Z2=kZ(IVE) CALL CENTPR(EX,EY,EZ,X2,Y2,Z2,XX,YY) IX2=IFIX(XX) IY2=IFIX(YY) CALL DRLINE(IX1,IY1,IX2,IY2) 20 CONTINUE * * --- draw the right diagram --- * CALL DRSTRI * * --- plot the input points --- * DO 30 II=1,NOFP X1=FLOAT(kX(II)) Y1=FLOAT(kY(II)) Z1=FLOAT(kZ(II)) CALL CENTPR(-EX,EY,EZ,X1,Y1,Z1,XX,YY) IXX=IFIX(XX) IYY=IFIX(YY) CALL DRDOT(IXX,IYY,RADIUS) 30 CONTINUE * * --- draw the edges --- * DO 40 II=1,KPARA(22) * --- skip unnecessary edges --- if(kecol(ii).eq.inhibi) goto 40 * --- compute the start vertex --- IVS=KeSV(II) X1=kX(IVS) Y1=kY(IVS) Z1=kZ(IVS) CALL CENTPR(-EX,EY,EZ,X1,Y1,Z1,XX,YY) IX1=IFIX(XX) IY1=IFIX(YY) * --- compute the end vertex --- IVE=KeEV(II) X2=kX(IVE) Y2=kY(IVE) Z2=kZ(IVE) CALL CENTPR(-EX,EY,EZ,X2,Y2,Z2,XX,YY) IX2=IFIX(XX) IY2=IFIX(YY) CALL DRLINE(IX1,IY1,IX2,IY2) 40 CONTINUE 1 continue CALL DRCOMP * return end * * subroutine drlag3 * ------------------------------------------------------------ * * * * draw the three-dimensional Laguere Delaunay * * tetrahedrization * * * * ------------------------------------------------------------ * common /ch40/ kpara(30) common /ch41a/ kx(1024),ky(1024),kz(1024),ku(1024) common /ch41b/ kvt(1024),kname(1024),mol(1024) common /ch42/ kesv(20480),keev(20480),kecol(20480) common /ch43a/ kfv(3,20480),kft(2,20480),kfnt(20480) common /ch43b/ idface(20480) common /ch44/ ktv(4,10240),ktf(4,10240),ktcol(10240) common /ch4f1/ fltonv(333) common /ch4f2/ fltstk(10240),fltdel(10240) common /ch4f3/ flvert(1024) common /ch4w1/ istack(333),itlist(333),idlist(333) common /ch4w2/ islist(333),idflst(333) common /ch4w3/ ivlist(333) common /lague/ radius(1024) logical fltonv,fltstk,fltdel,flvert * DATA WIDTH/0.5/ DATA EYEX/1.0/,EYEY/0.0/,EYEZ/8.0/ data x0/0.0/,y0/0.0/,z0/0.0/ DATA THETA/1.0/ ixymax=kpara(6) nofp=kpara(7) XYMAX=FLOAT(IXYMAX)*0.9 EX=XYMAX*EYEX EY=XYMAX*EYEY EZ=XYMAX*EYEZ * * * --- classify the edges --- * call etcolo * * --- draw the front side in the above half and * the rear side in the below half --- * do 1 ifront=1,2 iup=ifront if(ifront.eq.1) then inhibi=2 else inhibi=1 endif * * --- draw the left diagram --- * CALL DRSTLE(-IXYMAX,IXYMAX,-IXYMAX,IXYMAX,iup) CALL LWIDTH(WIDTH) * * --- plot the input spheres --- * c if(ifront.eq.2) goto 11 DO 10 II=1,NOFP X1=FLOAT(kX(II)) Y1=FLOAT(kY(II)) Z1=FLOAT(kZ(II)) CALL CENTPR(EX,EY,EZ,X1,Y1,Z1,XX,YY) IXX=IFIX(XX) IYY=IFIX(YY) r1=radius(ii) CALL CENTPR(EX,EY,EZ,r1,y0,z0,rr,YY) irr=ifix(rr) CALL drcirc(IXX,IYY,irr) 10 CONTINUE 11 continue * * --- draw the edges --- * DO 20 II=1,KPARA(22) * --- skip unnecessary edges --- if(kecol(ii).eq.inhibi) goto 20 * --- fetch the spheres --- IVS=KeSV(II) X1=kX(IVS) Y1=kY(IVS) Z1=kZ(IVS) r1=radius(ivs) IVE=KeEV(II) X2=kX(IVE) Y2=kY(IVE) Z2=kZ(IVE) c r2=radius(ive) c d=sqrt((x1-x2)**2+(y1-y2)**2+(z1-z2)**2) c if(r1+r2.lt.d) goto 20 * --- compute the start vertex --- CALL CENTPR(EX,EY,EZ,X1,Y1,Z1,XX,YY) IX1=IFIX(XX) IY1=IFIX(YY) * --- compute the end vertex --- CALL CENTPR(EX,EY,EZ,X2,Y2,Z2,XX,YY) IX2=IFIX(XX) IY2=IFIX(YY) CALL DRLINE(IX1,IY1,IX2,IY2) 20 CONTINUE * * --- draw the right diagram --- * CALL DRSTRI * * --- plot the input points --- * c if(ifront.eq.2) goto 31 DO 30 II=1,NOFP X1=FLOAT(kX(II)) Y1=FLOAT(kY(II)) Z1=FLOAT(kZ(II)) CALL CENTPR(-EX,EY,EZ,X1,Y1,Z1,XX,YY) IXX=IFIX(XX) IYY=IFIX(YY) r1=radius(ii) CALL CENTPR(EX,EY,EZ,r1,y0,z0,rr,YY) irr=ifix(rr) CALL drcirc(IXX,IYY,irr) 30 CONTINUE 31 continue * * --- draw the edges --- * DO 40 II=1,KPARA(22) * --- skip unnecessary edges --- if(kecol(ii).eq.inhibi) goto 40 * --- fetch the spheres --- IVS=KeSV(II) X1=kX(IVS) Y1=kY(IVS) Z1=kZ(IVS) r1=radius(ivs) IVE=KeEV(II) X2=kX(IVE) Y2=kY(IVE) Z2=kZ(IVE) c r2=radius(ive) c d=sqrt((x1-x2)**2+(y1-y2)**2+(z1-z2)**2) c if(r1+r2.lt.d) goto 40 * --- compute the start vertex --- CALL CENTPR(-EX,EY,EZ,X1,Y1,Z1,XX,YY) IX1=IFIX(XX) IY1=IFIX(YY) * --- compute the end vertex --- CALL CENTPR(-EX,EY,EZ,X2,Y2,Z2,XX,YY) IX2=IFIX(XX) IY2=IFIX(YY) CALL DRLINE(IX1,IY1,IX2,IY2) 40 CONTINUE 1 continue CALL DRCOMP * return end subroutine chull4 * ------------------------------------------------------------ * * * * Generate the convex hull of points in the * * four-dimensional space * * * * Copyright in 1997 by Kokichi Sugihara * * * * ------------------------------------------------------------ * common /ch40/ kpara(30) common /ch41a/ kx(1024),ky(1024),kz(1024),ku(1024) common /ch41b/ kvt(1024),kname(1024),mol(1024) common /ch41c/ kvact(1024) common /ch42/ kesv(20480),keev(20480),kecol(20480) common /ch43a/ kfv(3,20480),kft(2,20480),kfnt(20480) common /ch43b/ idface(20480) common /ch44/ ktv(4,10240),ktf(4,10240),ktcol(10240) common /ch4f1/ fltonv(333) common /ch4f2/ fltstk(10240),fltdel(10240) common /ch4f3/ flvert(1024) common /ch4w1/ istack(333),itlist(333),idlist(333) common /ch4w2/ islist(333),idflst(333) common /ch4w3/ ivlist(333) logical fltonv,fltstk,fltdel,flvert logical kvact * * --- check the number of points --- * nofp=kpara(7) IF(NOFP.lt.5) GOTO 700 IF(KPARA(7)+1.GT.KPARA(1)) GOTO 701 * * --- check the range of the coordinates --- * ixymax=kpara(6) DO 4 II=1,NOFP IF(IABS(kx(II)).GT.IXYMAX) THEN NP=II IXII=kx(II) GOTO 702 ENDIF IF(IABS(ky(II)).GT.IXYMAX) THEN NP=II IYII=ky(II) GOTO 702 ENDIF IF(IABS(kz(II)).GT.IXYMAX) THEN NP=II IZII=kz(II) GOTO 702 ENDIF IF(IABS(ku(II)).GT.IXYMAX) THEN NP=II IuII=ku(II) GOTO 702 ENDIF 4 CONTINUE * * --- check whether kname is consistent --- * CALL HPSORT(NOFP,kname,kname,MOL) DO 5 II=1,NOFP-1 IF(kname(MOL(II)).EQ.kname(MOL(II+1))) GOTO 6 5 CONTINUE GOTO 10 * * --- assign kname automatically --- * 6 CONTINUE DO 7 II=1,NOFP kname(II)=II 7 CONTINUE * * --- sort the input points in the x-coordinates --- * 10 continue call hpsort(nofp,kx,kname,mol) * * --- set the initial state of the data structure --- * call kinit * * --- construct the simplex for the first five points --- * call make5(mol) c write(6,678) c 678 format('exit from make5') * call monit * * --- add the other points one by one --- * if(nofp.eq.5) goto 30 ivlast=mol(5) do 20 i=6,nofp ivnew=mol(i) call vadd(ivlast,ivnew) if(mod(i,10).eq.0) write(6,660) i 660 format('n =',i8) c write(7,780) ivnew c 780 format(/'new point is added; ivnew=',i4) c call monit ivlast=ivnew 20 continue * 30 continue * --- discriminate inactive vertices --- call empchk c write(6,679) c 679 format('exit from empchk') * * --- create edges --- call makeed c write(6,680) c 680 format('exit from makeed') * 100 continue return * 700 continue write(6,710) nofp 710 format('chull4: too few points. nofp=',i3) goto 777 701 continue write(6,711) nofp 711 format('chull4: too many points. nofp=',i10) goto 777 702 continue write(6,712) np,ixii,iyii,izii,iuii 712 format('chull4: coordinate is out of range.'/ + 'ii,ix,iy,iz,iu=',i5,4i15) 777 continue stop end * * subroutine kinit * ------------------------------------------------------------ * * * * generate the initial state for the data structure * * * * ------------------------------------------------------------ * common /ch40/ kpara(30) common /ch41a/ kx(1024),ky(1024),kz(1024),ku(1024) common /ch41b/ kvt(1024),kname(1024),mol(1024) common /ch42/ kesv(20480),keev(20480),kecol(20480) common /ch43a/ kfv(3,20480),kft(2,20480),kfnt(20480) common /ch43b/ idface(20480) common /ch44/ ktv(4,10240),ktf(4,10240),ktcol(10240) common /ch4f1/ fltonv(333) common /ch4f2/ fltstk(10240),fltdel(10240) common /ch4f3/ flvert(1024) common /ch4w1/ istack(333),itlist(333),idlist(333) common /ch4w2/ islist(333),idflst(333) logical fltonv,fltstk,fltdel,flvert * -------------------------------------------------------------- * kpara(1)=kvmax maxmum number of the vertices * kpara(2)=kemax maxmum number of the edges * kpara(3)=kfmax maxmum number of the faces * kpara(4)=ktmax maxmum number of the tetrahedra * kpara(5)=kpmax maxmum number of the items * kpara(12) pointer to the empty edge list * kpara(13) pointer to the empty face list * kpara(14) pointer to the empty tetrahedron list * kpara(22) maximum index of edges ever used * kpara(23) maximum index of faces ever used * kpara(24) maximum index of tetrahedra ever used * kvt vertex-tetra pointer * kname distinct number assigned to vertices * mol list of vertices sorted in the x coordinates * kesv start vertex of an edge * keev end vertex of an edge * kecol 1: Delaunay edge, 2: farthest-point Delaunay edge, * 3: edge belong to the both * kfv three vertices on a face * kft two tetras incident to a face * kfnt new tetra created in one side of the face * idface type of a face (0: unchanged, 1: boundary face, * 3: face to be deleted) * ktv four vertices on a tetra * ktf four faces on a tetra * ktcol color of the tetrahedra (1:visible, 2:invisible) * fltonv flag used in subroutine tonv * fltstk flag for the tetrahedra that are put in istack * (used in vadd) * fltdel flag for the tetrahedra to be deleted * (used in vadd) * flvert flag for vertices that are checked * (used in makeed) * istack stack for tetrahedra to be checked * itlist list of tetrahedra incident to a vertex * idlist list of tetrahedra whose fltdel is changed * islist list of tetrahedra whose fltstk is changed * idflst list of faces whose idface is changed * ivlist list of vertices whose flvert is changed * (used in makeed) * -------------------------------------------------------------- kemax=kpara(2) kfmax=kpara(3) ktmax=kpara(4) kpmax=kpara(5) kpara(12)=1 kpara(13)=1 kpara(14)=1 kpara(22)=0 kpara(23)=0 kpara(24)=0 * do 10 i=1,kemax-1 kesv(i)=i+1 10 continue kesv(kemax)=0 * do 20 i=1,kfmax-1 kfv(1,i)=i+1 20 continue kfv(1,kfmax)=0 * do 30 i=1,ktmax-1 ktv(1,i)=i+1 30 continue ktv(1,ktmax)=0 * do 40 i=1,kpmax fltonv(i)=.FALSE. 40 continue * return end * * function knewe(idummy) * ------------------------------------------------------------ * * * * generate a new edge * * * * ------------------------------------------------------------ * common /ch40/ kpara(30) common /ch41a/ kx(1024),ky(1024),kz(1024),ku(1024) common /ch41b/ kvt(1024),kname(1024),mol(1024) common /ch42/ kesv(20480),keev(20480),kecol(20480) common /ch43a/ kfv(3,20480),kft(2,20480),kfnt(20480) common /ch43b/ idface(20480) common /ch44/ ktv(4,10240),ktf(4,10240),ktcol(10240) common /ch4f1/ fltonv(333) common /ch4f2/ fltstk(10240),fltdel(10240) common /ch4f3/ flvert(1024) common /ch4w1/ istack(333),itlist(333),idlist(333) common /ch4w2/ islist(333),idflst(333) logical fltonv,fltstk,fltdel,flvert * -------------------------------------------------------------- ip=idummy ip=kpara(12) if(ip.ne.0) then knewe=ip kpara(12)=kesv(ip) if(ip.gt.kpara(22)) kpara(22)=ip return else write(6,700) write(7,700) 700 format('knewe: edges overflow') stop endif end * * function knewf(idummy) * ------------------------------------------------------------ * * * * generate a new face * * * * ------------------------------------------------------------ * common /ch40/ kpara(30) common /ch41a/ kx(1024),ky(1024),kz(1024),ku(1024) common /ch41b/ kvt(1024),kname(1024),mol(1024) common /ch42/ kesv(20480),keev(20480),kecol(20480) common /ch43a/ kfv(3,20480),kft(2,20480),kfnt(20480) common /ch43b/ idface(20480) common /ch44/ ktv(4,10240),ktf(4,10240),ktcol(10240) common /ch4f1/ fltonv(333) common /ch4f2/ fltstk(10240),fltdel(10240) common /ch4f3/ flvert(1024) common /ch4w1/ istack(333),itlist(333),idlist(333) common /ch4w2/ islist(333),idflst(333) logical fltonv,fltstk,fltdel,flvert * -------------------------------------------------------------- ip=idummy ip=kpara(13) if(ip.eq.0) then write(6,700) write(7,700) 700 format('knewf: faces overflow') stop endif knewf=ip kpara(13)=kfv(1,ip) if(ip.gt.kpara(23)) kpara(23)=ip * do 10 i=1,3 kfv(i,ip)=0 10 continue do 11 i=1,2 kft(i,ip)=0 11 continue kfnt(ip)=0 idface(ip)=0 * return end * * function knewt(idummy) * ------------------------------------------------------------ * * * * generate a new tetrahedron * * * * ------------------------------------------------------------ * common /ch40/ kpara(30) common /ch41a/ kx(1024),ky(1024),kz(1024),ku(1024) common /ch41b/ kvt(1024),kname(1024),mol(1024) common /ch42/ kesv(20480),keev(20480),kecol(20480) common /ch43a/ kfv(3,20480),kft(2,20480),kfnt(20480) common /ch43b/ idface(20480) common /ch44/ ktv(4,10240),ktf(4,10240),ktcol(10240) common /ch4f1/ fltonv(333) common /ch4f2/ fltstk(10240),fltdel(10240) common /ch4f3/ flvert(1024) common /ch4w1/ istack(333),itlist(333),idlist(333) common /ch4w2/ islist(333),idflst(333) logical fltonv,fltstk,fltdel,flvert * -------------------------------------------------------------- ip=idummy ip=kpara(14) if(ip.eq.0) then write(6,700) write(7,700) 700 format('knewt: tetrahedra overflow') stop endif knewt=ip kpara(14)=ktv(1,ip) if(ip.gt.kpara(24)) kpara(24)=ip * do 10 i=1,4 ktv(i,ip)=0 ktf(i,ip)=0 10 continue fltstk(ip)=.false. fltdel(ip)=.false. * return end * * subroutine kbacke(iedge) * ------------------------------------------------------------ * * * * return the edge to the empty-edge list * * * * ------------------------------------------------------------ * common /ch40/ kpara(30) common /ch41a/ kx(1024),ky(1024),kz(1024),ku(1024) common /ch41b/ kvt(1024),kname(1024),mol(1024) common /ch42/ kesv(20480),keev(20480),kecol(20480) common /ch43a/ kfv(3,20480),kft(2,20480),kfnt(20480) common /ch43b/ idface(20480) common /ch44/ ktv(4,10240),ktf(4,10240),ktcol(10240) common /ch4f1/ fltonv(333) common /ch4f2/ fltstk(10240),fltdel(10240) common /ch4f3/ flvert(1024) common /ch4w1/ istack(333),itlist(333),idlist(333) common /ch4w2/ islist(333),idflst(333) logical fltonv,fltstk,fltdel,flvert * -------------------------------------------------------------- kesv(iedge)=kpara(12) kpara(12)=iedge return end * * subroutine kbackf(iface) * ------------------------------------------------------------ * * * * return the face to the empty-face list * * * * ------------------------------------------------------------ * common /ch40/ kpara(30) common /ch41a/ kx(1024),ky(1024),kz(1024),ku(1024) common /ch41b/ kvt(1024),kname(1024),mol(1024) common /ch42/ kesv(20480),keev(20480),kecol(20480) common /ch43a/ kfv(3,20480),kft(2,20480),kfnt(20480) common /ch43b/ idface(20480) common /ch44/ ktv(4,10240),ktf(4,10240),ktcol(10240) common /ch4f1/ fltonv(333) common /ch4f2/ fltstk(10240),fltdel(10240) common /ch4f3/ flvert(1024) common /ch4w1/ istack(333),itlist(333),idlist(333) common /ch4w2/ islist(333),idflst(333) logical fltonv,fltstk,fltdel,flvert * -------------------------------------------------------------- kfv(1,iface)=kpara(13) kpara(13)=iface return end * * subroutine kbackt(itetra) * ------------------------------------------------------------ * * * * return the tetrahedron to the empty-tetrahedron list * * * * ------------------------------------------------------------ * common /ch40/ kpara(30) common /ch41a/ kx(1024),ky(1024),kz(1024),ku(1024) common /ch41b/ kvt(1024),kname(1024),mol(1024) common /ch42/ kesv(20480),keev(20480),kecol(20480) common /ch43a/ kfv(3,20480),kft(2,20480),kfnt(20480) common /ch43b/ idface(20480) common /ch44/ ktv(4,10240),ktf(4,10240),ktcol(10240) common /ch4f1/ fltonv(333) common /ch4f2/ fltstk(10240),fltdel(10240) common /ch4f3/ flvert(1024) common /ch4w1/ istack(333),itlist(333),idlist(333) common /ch4w2/ islist(333),idflst(333) logical fltonv,fltstk,fltdel,flvert * -------------------------------------------------------------- ktv(1,itetra)=kpara(14) kpara(14)=itetra return end * * subroutine make5(iv5) * ------------------------------------------------------------ * * * * generate the four-dimensional simplex * * * * ------------------------------------------------------------ * common /ch40/ kpara(30) common /ch41a/ kx(1024),ky(1024),kz(1024),ku(1024) common /ch41b/ kvt(1024),kname(1024),mol(1024) common /ch42/ kesv(20480),keev(20480),kecol(20480) common /ch43a/ kfv(3,20480),kft(2,20480),kfnt(20480) common /ch43b/ idface(20480) common /ch44/ ktv(4,10240),ktf(4,10240),ktcol(10240) common /ch4f1/ fltonv(333) common /ch4f2/ fltstk(10240),fltdel(10240) common /ch4f3/ flvert(1024) common /ch4w1/ istack(333),itlist(333),idlist(333) common /ch4w2/ islist(333),idflst(333) logical fltonv,fltstk,fltdel,flvert * -------------------------------------------------------------- dimension iv5(1) * * --- check the orientation --- * iori=jor5(iv5) if(iori.eq.1) then jv3=iv5(4) jv4=iv5(3) else jv3=iv5(3) jv4=iv5(4) endif jv1=iv5(1) jv2=iv5(2) jv5=iv5(5) * * --- new faces --- * if1=knewf(idummy) if2=knewf(idummy) if3=knewf(idummy) if4=knewf(idummy) if5=knewf(idummy) if6=knewf(idummy) if7=knewf(idummy) if8=knewf(idummy) if9=knewf(idummy) if10=knewf(idummy) * * --- new tetrahedra --- * it1=knewt(idummy) it2=knewt(idummy) it3=knewt(idummy) it4=knewt(idummy) it5=knewt(idummy) * ktv(1,it1)=jv1 ktv(2,it1)=jv2 ktv(3,it1)=jv3 ktv(4,it1)=jv4 ktv(1,it2)=jv2 ktv(2,it2)=jv3 ktv(3,it2)=jv4 ktv(4,it2)=jv5 ktv(1,it3)=jv3 ktv(2,it3)=jv4 ktv(3,it3)=jv5 ktv(4,it3)=jv1 ktv(1,it4)=jv4 ktv(2,it4)=jv5 ktv(3,it4)=jv1 ktv(4,it4)=jv2 ktv(1,it5)=jv5 ktv(2,it5)=jv1 ktv(3,it5)=jv2 ktv(4,it5)=jv3 * --- vertices --- kvt(jv1)=it1 kvt(jv2)=it1 kvt(jv3)=it1 kvt(jv4)=it1 kvt(jv5)=it2 * --- faces --- kfv(1,if1)=jv1 kfv(2,if1)=jv2 kfv(3,if1)=jv3 kfv(1,if2)=jv1 kfv(2,if2)=jv2 kfv(3,if2)=jv4 kfv(1,if3)=jv2 kfv(2,if3)=jv3 kfv(3,if3)=jv4 kfv(1,if4)=jv1 kfv(2,if4)=jv3 kfv(3,if4)=jv4 kfv(1,if5)=jv1 kfv(2,if5)=jv2 kfv(3,if5)=jv5 kfv(1,if6)=jv2 kfv(2,if6)=jv3 kfv(3,if6)=jv5 kfv(1,if7)=jv1 kfv(2,if7)=jv3 kfv(3,if7)=jv5 kfv(1,if8)=jv1 kfv(2,if8)=jv4 kfv(3,if8)=jv5 kfv(1,if9)=jv2 kfv(2,if9)=jv4 kfv(3,if9)=jv5 kfv(1,if10)=jv3 kfv(2,if10)=jv4 kfv(3,if10)=jv5 * kft(1,if1)=it1 kft(2,if1)=it5 kft(1,if2)=it1 kft(2,if2)=it4 kft(1,if3)=it1 kft(2,if3)=it2 kft(1,if4)=it1 kft(2,if4)=it3 kft(1,if5)=it4 kft(2,if5)=it5 kft(1,if6)=it2 kft(2,if6)=it5 kft(1,if7)=it3 kft(2,if7)=it5 kft(1,if8)=it3 kft(2,if8)=it4 kft(1,if9)=it2 kft(2,if9)=it4 kft(1,if10)=it2 kft(2,if10)=it3 * --- tetrahedra --- ktf(1,it1)=if1 ktf(2,it1)=if2 ktf(3,it1)=if3 ktf(4,it1)=if4 ktf(1,it2)=if3 ktf(2,it2)=if6 ktf(3,it2)=if9 ktf(4,it2)=if10 ktf(1,it3)=if4 ktf(2,it3)=if7 ktf(3,it3)=if8 ktf(4,it3)=if10 ktf(1,it4)=if2 ktf(2,it4)=if5 ktf(3,it4)=if8 ktf(4,it4)=if9 ktf(1,it5)=if1 ktf(2,it5)=if5 ktf(3,it5)=if6 ktf(4,it5)=if7 * return end * * subroutine vadd(ivlast,ivnew) * ------------------------------------------------------------ * * * * add one point and augment the convex hull * * * * ------------------------------------------------------------ * common /ch40/ kpara(30) common /ch41a/ kx(1024),ky(1024),kz(1024),ku(1024) common /ch41b/ kvt(1024),kname(1024),mol(1024) common /ch42/ kesv(20480),keev(20480),kecol(20480) common /ch43a/ kfv(3,20480),kft(2,20480),kfnt(20480) * kfnt --- new tetra incident to a face common /ch43b/ idface(20480) common /ch44/ ktv(4,10240),ktf(4,10240),ktcol(10240) common /ch4f1/ fltonv(333) common /ch4f2/ fltstk(10240),fltdel(10240) common /ch4f3/ flvert(1024) common /ch4w1/ istack(333),itlist(333),idlist(333) common /ch4w2/ islist(333),idflst(333) logical fltonv,fltstk,fltdel,flvert * -------------------------------------------------------------- * dimension iv5(5) * * --- push the tetrahedra incident to ivlast --- * kpmax=kpara(5) call tonv(ivlast,ntetra) do 1 ii=1,ntetra it=itlist(ii) fltstk(it)=.true. * tetra it is visited istack(ii)=it islist(ii)=it * memo for clearing the flag fltstk 1 continue ipoint=ntetra * pointer for istack nofs=ntetra * pointer for islist nofd=0 * pointer for idlist (tetras to be deleted) nofdf=0 * pointer for idflst (faces to be deleted) * * --- pop a tetrahedron and check --- * 10 continue if(ipoint.eq.0) goto 50 itet=istack(ipoint) ipoint=ipoint-1 call cansee(ivnew,itet,iyes) c write(7,789) ivnew,itet,iyes c 789 format('ivnew,itet,iyes=',3i4) if(iyes.eq.0) goto 20 fltdel(itet)=.true. * tetra to be deleted nofd=nofd+1 if(nofd.gt.kpmax) goto 704 idlist(nofd)=itet do 11 ii=1,4 iface=ktf(ii,itet) jtet=kft(1,iface) if(jtet.eq.itet) jtet=kft(2,iface) if(fltstk(jtet).eqv..false.) then fltstk(jtet)=.true. ipoint=ipoint+1 if(ipoint.gt.kpmax) goto 700 istack(ipoint)=jtet nofs=nofs+1 if(nofs.gt.kpmax) goto 701 islist(nofs)=jtet endif 11 continue 20 continue goto 10 * * --- find relevant faces --- * 50 continue do 51 ii=1,nofd it=idlist(ii) do 52 jj=1,4 iface=ktf(jj,it) idface(iface)=idface(iface)+1 if(idface(iface).eq.1) then nofdf=nofdf+1 if(nofdf.gt.kpmax) goto 702 idflst(nofdf)=iface * idflst = 2 --- to be deleted * = 1 --- one side tetra should be deleted endif 52 continue 51 continue c do 787 iface=1,kpara(23) c write(7,788) iface,idface(iface) c 788 format('iface,idface=',2i4) c 787 continue * * --- create new tetrahedra --- * nofnt=0 do 60 ii=1,nofdf iface=idflst(ii) if(idface(iface).eq.2) goto 60 newt=knewt(idummy) nofnt=nofnt+1 if(nofnt.gt.kpmax) goto 705 itlist(nofnt)=newt kfnt(iface)=newt iv1=kfv(1,iface) iv2=kfv(2,iface) iv3=kfv(3,iface) ktv(3,newt)=iv3 ktv(4,newt)=ivnew c write(7,785) ii,iface,iv1,iv2,iv3,newt c 785 format('ii,iface,vvv,newt=',6i4) * * --- decide the orientation of the new tetrahedron --- * itoppo=kft(1,iface) if(fltdel(itoppo).eqv..true.) itoppo=kft(2,iface) do 61 mm=1,4 jjv=ktv(mm,itoppo) if(jjv.eq.iv1) goto 61 if(jjv.eq.iv2) goto 61 if(jjv.eq.iv3) goto 61 ivoppo=jjv goto 62 61 continue goto 703 * 62 continue iv5(1)=iv1 iv5(2)=iv2 iv5(3)=iv3 iv5(4)=ivnew iv5(5)=ivoppo iori=jor5(iv5) if(iori.eq.1) then ktv(1,newt)=iv2 ktv(2,newt)=iv1 else ktv(1,newt)=iv1 ktv(2,newt)=iv2 endif * * orientation of tetra = 1 when seen from outside * of the 4-d convex hull * = -1 when seen from inside * ktf(4,newt)=iface * --- change the vertex-tetrahedron pointer --- do 63 mm=1,3 mv=kfv(mm,iface) kvt(mv)=newt 63 continue 60 continue kvt(ivnew)=newt * * --- create new faces --- * do 70 ii=1,nofdf iface=idflst(ii) if(idface(iface).eq.2) goto 70 itet=kfnt(iface) do 71 m=1,3 c write(7,786) iface,itet,m,ktf(m,itet) c 786 format('iface,itet,m,ktf=',4i4) if(ktf(m,itet).ne.0) goto 71 iv1=ktv(m,itet) if(m.lt.3) then iv2=ktv(m+1,itet) else iv2=ktv(1,itet) endif call opface(iface,iv1,iv2,ifoppo) c write(7,783) iface,iv1,iv2,ifoppo c 783 format('iface,iv1,iv2,ifoppo=',4i4) jtet=kfnt(ifoppo) newf=knewf(idummy) kft(1,newf)=itet kft(2,newf)=jtet kfv(1,newf)=iv1 kfv(2,newf)=iv2 kfv(3,newf)=ivnew ktf(m,itet)=newf if(iv1.eq.ktv(1,jtet)) then if(iv2.eq.ktv(2,jtet)) then mitem=1 else mitem=3 endif elseif(iv1.eq.ktv(2,jtet)) then if(iv2.eq.ktv(3,jtet)) then mitem=2 else mitem=1 endif else if(iv2.eq.ktv(1,jtet)) then mitem=3 else mitem=2 endif endif ktf(mitem,jtet)=newf c write(7,784) jtet,mitem,newf c 784 format('vadd:jtet,mitem,newf=',3i4) 71 continue 70 continue * * --- delete old tetrahedra and old faces --- * do 80 jj=1,nofd it=idlist(jj) call kbackt(it) 80 continue do 81 ii=1,nofdf iface=idflst(ii) if(idface(iface).eq.2) call kbackf(iface) 81 continue * * --- clear f-t pointers of the face on the boundary --- * do 82 ii=1,nofdf iface=idflst(ii) if(idface(iface).eq.2) goto 82 itet=kft(1,iface) if(fltdel(itet).eqv..true.) then mitem=1 else mitem=2 endif kft(mitem,iface)=kfnt(iface) 82 continue * * --- clear flags --- * do 91 ii=1,nofd it=idlist(ii) fltdel(it)=.false. 91 continue do 92 ii=1,nofs it=islist(ii) fltstk(it)=.false. 92 continue do 93 ii=1,nofdf iface=idflst(ii) idface(iface)=0 93 continue * return 700 continue write(6,710) 710 format('vadd: istack overflows') goto 777 701 continue write(6,711) 711 format('vadd: islist overflows') goto 777 702 continue write(6,712) 712 format('vadd: idflst overflows') goto 777 703 continue write(6,713) 713 format('vadd: tetra itoppo does not contain jjv') goto 777 704 continue write(6,714) nofd 714 format('vadd: idlist overflows, nofd =',i10) goto 777 705 continue write(6,715) nofnt 715 format('vadd: itlist overflows, nofnt =',i10) 777 continue stop end * * subroutine tonv(iv,ntetra) * ------------------------------------------------------------ * * * * list up the tetrahedra incident to the vertex iv * * * * ntetra = number of tetrahedra * * The result of enumeration is stored in itlist * * * * ------------------------------------------------------------ * common /ch40/ kpara(30) common /ch41a/ kx(1024),ky(1024),kz(1024),ku(1024) common /ch41b/ kvt(1024),kname(1024),mol(1024) common /ch42/ kesv(20480),keev(20480),kecol(20480) common /ch43a/ kfv(3,20480),kft(2,20480),kfnt(20480) common /ch43b/ idface(20480) common /ch44/ ktv(4,10240),ktf(4,10240),ktcol(10240) common /ch4f1/ fltonv(333) common /ch4f2/ fltstk(10240),fltdel(10240) common /ch4f3/ flvert(1024) common /ch4w1/ istack(333),itlist(333),idlist(333) common /ch4w2/ islist(333),idflst(333) logical fltonv,fltstk,fltdel,flvert * kpmax=kpara(5) * its=kvt(iv) ntetra=1 itlist(ntetra)=its fltonv(its)=.TRUE. ipoint=1 istack(ipoint)=its * * --- pop a tetrahedron from istack --- * 1 continue if(ipoint.eq.0) goto 90 it=istack(ipoint) ipoint=ipoint-1 do 10 ii=1,4 iface=ktf(ii,it) * * --- check if iface has vertex iv --- * do 11 jj=1,3 jv=kfv(jj,iface) if(jv.eq.iv) goto 12 11 continue * --- iface does not have vertex iv --- goto 10 * --- iface has vertex iv --- 12 continue jt=kft(1,iface) if(jt.eq.it) jt=kft(2,iface) if(fltonv(jt).eqv..false.) then fltonv(jt)=.true. ntetra=ntetra+1 if(ntetra.gt.kpmax) goto 700 itlist(ntetra)=jt ipoint=ipoint+1 if(ipoint.gt.kpmax) goto 701 istack(ipoint)=jt endif 10 continue goto 1 * * --- clear the fltonv --- * 90 continue do 91 ii=1,ntetra fltonv(itlist(ii))=.false. 91 continue return * 700 continue write(6,710) ntetra 710 format('tonv: ntetra overflows. ntetra=',i6) goto 777 701 continue write(6,711) 711 format('tonv: istack overflows') 777 continue stop end * * subroutine opface(ifaces,iv1,iv2,ifoppo) * ------------------------------------------------------------ * * * * find the face that is the other side of the edge * * (iv1,iv2) on the boundary of the vacant space * * * * ------------------------------------------------------------ * common /ch40/ kpara(30) common /ch41a/ kx(1024),ky(1024),kz(1024),ku(1024) common /ch41b/ kvt(1024),kname(1024),mol(1024) common /ch42/ kesv(20480),keev(20480),kecol(20480) common /ch43a/ kfv(3,20480),kft(2,20480),kfnt(20480) common /ch43b/ idface(20480) common /ch44/ ktv(4,10240),ktf(4,10240),ktcol(10240) common /ch4f1/ fltonv(333) common /ch4f2/ fltstk(10240),fltdel(10240) common /ch4f3/ flvert(1024) common /ch4w1/ istack(333),itlist(333),idlist(333) common /ch4w2/ islist(333),idflst(333) logical fltonv,fltstk,fltdel,flvert * dimension ivaxis(2) c write(7,780) ifaces,iv1,iv2 c 780 format('enter opface -- ifaces,iv1,iv2=',3i4) ivaxis(1)=iv1 ivaxis(2)=iv2 it=kft(1,ifaces) iface=ifaces 1 continue do 10 ii=1,4 jface=ktf(ii,it) if(jface.eq.iface) goto 10 do 11 jj=1,2 do 13 jjj=1,3 if(ivaxis(jj).eq.kfv(jjj,jface)) goto 11 13 continue goto 10 11 continue goto 20 10 continue goto 700 20 continue c write(7,781) jface,idface(jface) c 781 format('jface,idf=',2i4) if(idface(jface).eq.1) then * --- the opposite face is found --- ifoppo=jface return endif * * --- next face --- * jt=kft(1,jface) if(jt.eq.it) jt=kft(2,jface) it=jt if(jface.eq.ifaces) goto 701 iface=jface goto 1 * 700 continue write(6,710) jface 710 format('opface: next face cannot be found, jface=',i7) goto 777 701 continue write(6,711) jface 711 format('opface: return to the start face, jface=',i8) 777 continue stop end * * function jor5(iv5) * ------------------------------------------------------------ * * * * compute the orientation using the symbolic perturbation * * * * ------------------------------------------------------------ * common /ch40/ kpara(30) common /ch41a/ kx(1024),ky(1024),kz(1024),ku(1024) common /ch41b/ kvt(1024),kname(1024),mol(1024) common /ch42/ kesv(20480),keev(20480),kecol(20480) common /ch43a/ kfv(3,20480),kft(2,20480),kfnt(20480) common /ch43b/ idface(20480) common /ch44/ ktv(4,10240),ktf(4,10240),ktcol(10240) common /ch4f1/ fltonv(333) common /ch4f2/ fltstk(10240),fltdel(10240) common /ch4f3/ flvert(1024) common /ch4w1/ istack(333),itlist(333),idlist(333) common /ch4w2/ islist(333),idflst(333) common /ch4w3/ ivlist(333) logical fltonv,fltstk,fltdel,flvert * dimension iv5(5),jv5(5) dimension l211(6),l212(6),l213(6),l214(6) dimension l311(6),l312(6),l313(6),l314(6) dimension l411(6),l412(6),l413(6),l414(6) dimension l511(6),l512(6),l513(6),l514(6) dimension l321(6),l322(6),l323(6),l324(6) dimension l421(6),l422(6),l423(6),l424(6) dimension l521(6),l522(6),l523(6),l524(6) dimension l431(6),l432(6),l433(6),l434(6) dimension l531(6),l532(6),l533(6),l534(6) dimension l541(6),l542(6),l543(6),l544(6) dimension lm1(6),lm2(6),lm3(6),ldeta(6),ldetb(6) data length/6/ * * --- reorder the vertices --- * call sortbn(5,iv5,kname,jv5,ichang) * * --- floating-point acceleration --- * isign=jacce5(jv5) c write(6,677) isign,ichang c 677 format('jor5: float comp. isign,ichang =',2i4) isign=0 if(isign.ne.0) goto 100 * * --- set the long integers --- * ival=kx(jv5(2))-kx(jv5(1)) call ltrans(ival,l211,length) ival=kx(jv5(3))-kx(jv5(1)) call ltrans(ival,l311,length) ival=kx(jv5(4))-kx(jv5(1)) call ltrans(ival,l411,length) ival=kx(jv5(5))-kx(jv5(1)) call ltrans(ival,l511,length) * ival=ky(jv5(2))-ky(jv5(1)) call ltrans(ival,l212,length) ival=ky(jv5(3))-ky(jv5(1)) call ltrans(ival,l312,length) ival=ky(jv5(4))-ky(jv5(1)) call ltrans(ival,l412,length) ival=ky(jv5(5))-ky(jv5(1)) call ltrans(ival,l512,length) * ival=kz(jv5(2))-kz(jv5(1)) call ltrans(ival,l213,length) ival=kz(jv5(3))-kz(jv5(1)) call ltrans(ival,l313,length) ival=kz(jv5(4))-kz(jv5(1)) call ltrans(ival,l413,length) ival=kz(jv5(5))-kz(jv5(1)) call ltrans(ival,l513,length) * ival=ku(jv5(2))-ku(jv5(1)) call ltrans(ival,l214,length) ival=ku(jv5(3))-ku(jv5(1)) call ltrans(ival,l314,length) ival=ku(jv5(4))-ku(jv5(1)) call ltrans(ival,l414,length) ival=ku(jv5(5))-ku(jv5(1)) call ltrans(ival,l514,length) * * --- general case --- * call exmult(l211,l312,lm1) call exmult(l413,l514,lm2) call exmult(lm1,lm2,ldeta) * call exmult(l211,l312,lm1) call exmult(l513,l414,lm2) call exmult(lm1,lm2,lm3) call exsub(ldeta,lm3,ldetb) * call exmult(l211,l412,lm1) call exmult(l313,l514,lm2) call exmult(lm1,lm2,lm3) call exsub(ldetb,lm3,ldeta) * call exmult(l211,l412,lm1) call exmult(l513,l314,lm2) call exmult(lm1,lm2,lm3) call exadd(ldeta,lm3,ldetb) * call exmult(l211,l512,lm1) call exmult(l313,l414,lm2) call exmult(lm1,lm2,lm3) call exadd(ldetb,lm3,ldeta) * call exmult(l211,l512,lm1) call exmult(l413,l314,lm2) call exmult(lm1,lm2,lm3) call exsub(ldeta,lm3,ldetb) * call exmult(l311,l212,lm1) call exmult(l413,l514,lm2) call exmult(lm1,lm2,lm3) call exsub(ldetb,lm3,ldeta) * call exmult(l311,l212,lm1) call exmult(l513,l414,lm2) call exmult(lm1,lm2,lm3) call exadd(ldeta,lm3,ldetb) * call exmult(l311,l412,lm1) call exmult(l213,l514,lm2) call exmult(lm1,lm2,lm3) call exadd(ldetb,lm3,ldeta) * call exmult(l311,l412,lm1) call exmult(l513,l214,lm2) call exmult(lm1,lm2,lm3) call exsub(ldeta,lm3,ldetb) * call exmult(l311,l512,lm1) call exmult(l213,l414,lm2) call exmult(lm1,lm2,lm3) call exsub(ldetb,lm3,ldeta) * call exmult(l311,l512,lm1) call exmult(l413,l214,lm2) call exmult(lm1,lm2,lm3) call exadd(ldeta,lm3,ldetb) * call exmult(l411,l212,lm1) call exmult(l313,l514,lm2) call exmult(lm1,lm2,lm3) call exadd(ldetb,lm3,ldeta) * call exmult(l411,l212,lm1) call exmult(l513,l314,lm2) call exmult(lm1,lm2,lm3) call exsub(ldeta,lm3,ldetb) * call exmult(l411,l312,lm1) call exmult(l213,l514,lm2) call exmult(lm1,lm2,lm3) call exsub(ldetb,lm3,ldeta) * call exmult(l411,l312,lm1) call exmult(l513,l214,lm2) call exmult(lm1,lm2,lm3) call exadd(ldeta,lm3,ldetb) * call exmult(l411,l512,lm1) call exmult(l213,l314,lm2) call exmult(lm1,lm2,lm3) call exadd(ldetb,lm3,ldeta) * call exmult(l411,l512,lm1) call exmult(l313,l214,lm2) call exmult(lm1,lm2,lm3) call exsub(ldeta,lm3,ldetb) * call exmult(l511,l212,lm1) call exmult(l313,l414,lm2) call exmult(lm1,lm2,lm3) call exsub(ldetb,lm3,ldeta) * call exmult(l511,l212,lm1) call exmult(l413,l314,lm2) call exmult(lm1,lm2,lm3) call exadd(ldeta,lm3,ldetb) * call exmult(l511,l312,lm1) call exmult(l213,l414,lm2) call exmult(lm1,lm2,lm3) call exadd(ldetb,lm3,ldeta) * call exmult(l511,l312,lm1) call exmult(l413,l214,lm2) call exmult(lm1,lm2,lm3) call exsub(ldeta,lm3,ldetb) * call exmult(l511,l412,lm1) call exmult(l213,l314,lm2) call exmult(lm1,lm2,lm3) call exsub(ldetb,lm3,ldeta) * call exmult(l511,l412,lm1) call exmult(l313,l214,lm2) call exmult(lm1,lm2,lm3) call exadd(ldeta,lm3,ldetb) * isign=ldetb(2) if(isign.eq.1.or.isign.eq.-1) goto 100 * * --- degenerate case 1 --- * ival=kx(jv5(3))-kx(jv5(2)) call ltrans(ival,l321,length) ival=kx(jv5(4))-kx(jv5(2)) call ltrans(ival,l421,length) ival=kx(jv5(5))-kx(jv5(2)) call ltrans(ival,l521,length) * ival=ky(jv5(3))-ky(jv5(2)) call ltrans(ival,l322,length) ival=ky(jv5(4))-ky(jv5(2)) call ltrans(ival,l422,length) ival=ky(jv5(5))-ky(jv5(2)) call ltrans(ival,l522,length) * ival=kz(jv5(3))-kz(jv5(2)) call ltrans(ival,l323,length) ival=kz(jv5(4))-kz(jv5(2)) call ltrans(ival,l423,length) ival=kz(jv5(5))-kz(jv5(2)) call ltrans(ival,l523,length) * ival=ku(jv5(3))-ku(jv5(2)) call ltrans(ival,l324,length) ival=ku(jv5(4))-ku(jv5(2)) call ltrans(ival,l424,length) ival=ku(jv5(5))-ku(jv5(2)) call ltrans(ival,l524,length) * isign=-jor3(l322,l323,l324,l422,l423,l424, + l522,l523,l524) if(isign.eq.1.or.isign.eq.-1) goto 100 * * --- degenerate case 2 --- * isign=jor3(l312,l313,l314,l412,l413,l414, + l512,l513,l514) if(isign.eq.1.or.isign.eq.-1) goto 100 * * --- degenerate case 3 --- * isign=-jor3(l212,l213,l214,l412,l413,l414, + l512,l513,l514) if(isign.eq.1.or.isign.eq.-1) goto 100 * * --- degenerate case 4 --- * isign=jor3(l212,l213,l214,l312,l313,l314, + l512,l513,l514) if(isign.eq.1.or.isign.eq.-1) goto 100 * * --- degenerate case 5 --- * isign=jor3(l321,l323,l324,l421,l423,l424, + l521,l523,l524) if(isign.eq.1.or.isign.eq.-1) goto 100 * * --- degenerate case 6 --- * ival=kx(jv5(4))-kx(jv5(3)) call ltrans(ival,l431,length) ival=kx(jv5(5))-kx(jv5(3)) call ltrans(ival,l531,length) * ival=ky(jv5(4))-ky(jv5(3)) call ltrans(ival,l432,length) ival=ky(jv5(5))-ky(jv5(3)) call ltrans(ival,l532,length) * ival=kz(jv5(4))-kz(jv5(3)) call ltrans(ival,l433,length) ival=kz(jv5(5))-kz(jv5(3)) call ltrans(ival,l533,length) * ival=ku(jv5(4))-ku(jv5(3)) call ltrans(ival,l434,length) ival=ku(jv5(5))-ku(jv5(3)) call ltrans(ival,l534,length) * isign=-jor2(l433,l434,l533,l534) if(isign.eq.1.or.isign.eq.-1) goto 100 * * --- degenerate case 7 --- * isign=jor2(l423,l424,l523,l524) if(isign.eq.1.or.isign.eq.-1) goto 100 * * --- degenerate case 8 --- * isign=-jor2(l323,l324,l523,l524) if(isign.eq.1.or.isign.eq.-1) goto 100 * * --- degenerate case 9 --- * isign=-jor3(l311,l313,l314,l411,l413,l414, + l511,l513,l514) if(isign.eq.1.or.isign.eq.-1) goto 100 * * --- degenerate case 10 --- * isign=-jor2(l413,l414,l513,l514) if(isign.eq.1.or.isign.eq.-1) goto 100 * * --- degenerate case 11 --- * isign=jor2(l313,l314,l513,l514) if(isign.eq.1.or.isign.eq.-1) goto 100 * * --- degenerate case 12 --- * isign=jor3(l211,l213,l214,l411,l413,l414, + l511,l513,l514) if(isign.eq.1.or.isign.eq.-1) goto 100 * * --- degenerate case 13 --- * isign=-jor2(l213,l214,l513,l514) if(isign.eq.1.or.isign.eq.-1) goto 100 * * --- degenerate case 14 --- * isign=-jor3(l211,l213,l214,l311,l313,l314, + l511,l513,l514) if(isign.eq.1.or.isign.eq.-1) goto 100 * * --- degenerate case 15 --- * isign=-jor3(l321,l322,l324,l421,l422,l424, + l521,l522,l524) if(isign.eq.1.or.isign.eq.-1) goto 100 * * --- degenerate case 16 --- * isign=jor2(l432,l434,l532,l534) if(isign.eq.1.or.isign.eq.-1) goto 100 * * --- degenerate case 17 --- * isign=-jor2(l422,l424,l522,l524) if(isign.eq.1.or.isign.eq.-1) goto 100 * * --- degenerate case 18 --- * isign=jor2(l322,l324,l522,l524) if(isign.eq.1.or.isign.eq.-1) goto 100 * * --- degenerate case 19 --- * isign=-jor2(l431,l434,l531,l534) if(isign.eq.1.or.isign.eq.-1) goto 100 * * --- degenerate case 20 --- * ival=kx(jv5(5))-kx(jv5(4)) call ltrans(ival,l541,length) * ival=ky(jv5(5))-ky(jv5(4)) call ltrans(ival,l542,length) * ival=kz(jv5(5))-kz(jv5(4)) call ltrans(ival,l543,length) * ival=ku(jv5(5))-ku(jv5(4)) call ltrans(ival,l544,length) * isign=l544(2) if(isign.eq.1.or.isign.eq.-1) goto 100 * * --- degenerate case 21 --- * isign=-l534(2) if(isign.eq.1.or.isign.eq.-1) goto 100 * * --- degenerate case 22 --- * isign=jor2(l421,l424,l521,l524) if(isign.eq.1.or.isign.eq.-1) goto 100 * * --- degenerate case 23 --- * isign=l524(2) if(isign.eq.1.or.isign.eq.-1) goto 100 * * --- degenerate case 24 --- * isign=-jor2(l321,l324,l521,l524) if(isign.eq.1.or.isign.eq.-1) goto 100 * * --- degenerate case 25 --- * isign=jor3(l311,l312,l314,l411,l412,l414, + l511,l512,l514) if(isign.eq.1.or.isign.eq.-1) goto 100 * * --- degenerate case 26 --- * isign=jor2(l412,l414,l512,l514) if(isign.eq.1.or.isign.eq.-1) goto 100 * * --- degenerate case 27 --- * isign=-jor2(l312,l314,l512,l514) if(isign.eq.1.or.isign.eq.-1) goto 100 * * --- degenerate case 28 --- * isign=-jor2(l411,l414,l511,l514) if(isign.eq.1.or.isign.eq.-1) goto 100 * * --- degenerate case 29 --- * isign=-l514(2) if(isign.eq.1.or.isign.eq.-1) goto 100 * * --- degenerate case 30 --- * isign=jor2(l311,l314,l511,l514) if(isign.eq.1.or.isign.eq.-1) goto 100 * * --- degenerate case 31 --- * isign=-jor3(l211,l212,l214,l411,l412,l414, + l511,l512,l514) if(isign.eq.1.or.isign.eq.-1) goto 100 * * --- degenerate case 32 --- * isign=jor2(l212,l214,l512,l514) if(isign.eq.1.or.isign.eq.-1) goto 100 * * --- degenerate case 33 --- * isign=-jor2(l211,l214,l511,l514) if(isign.eq.1.or.isign.eq.-1) goto 100 * * --- degenerate case 34 --- * isign=jor3(l211,l212,l214,l311,l312,l314, + l511,l512,l514) if(isign.eq.1.or.isign.eq.-1) goto 100 * * --- degenerate case 35 --- * isign=jor3(l321,l322,l323,l421,l422,l423, + l521,l522,l523) if(isign.eq.1.or.isign.eq.-1) goto 100 * * --- degenerate case 36 --- * isign=-jor2(l432,l433,l532,l533) if(isign.eq.1.or.isign.eq.-1) goto 100 * * --- degenerate case 37 --- * isign=jor2(l422,l423,l522,l523) if(isign.eq.1.or.isign.eq.-1) goto 100 * * --- degenerate case 38 --- * isign=-jor2(l322,l323,l522,l523) if(isign.eq.1.or.isign.eq.-1) goto 100 * * --- degenerate case 39 --- * isign=jor2(l431,l433,l531,l533) if(isign.eq.1.or.isign.eq.-1) goto 100 * * --- degenerate case 40 --- * isign=-l543(2) if(isign.eq.1.or.isign.eq.-1) goto 100 * * --- degenerate case 41 --- * isign=l533(2) if(isign.eq.1.or.isign.eq.-1) goto 100 * * --- degenerate case 42 --- * isign=-jor2(l421,l423,l521,l523) if(isign.eq.1.or.isign.eq.-1) goto 100 * * --- degenerate case 43 --- * isign=-l523(2) if(isign.eq.1.or.isign.eq.-1) goto 100 * * --- degenerate case 44 --- * isign=jor2(l321,l323,l521,l523) if(isign.eq.1.or.isign.eq.-1) goto 100 * * --- degenerate case 45 --- * isign=-jor2(l431,l432,l531,l532) if(isign.eq.1.or.isign.eq.-1) goto 100 * * --- degenerate case 46 --- * isign=l542(2) if(isign.eq.1.or.isign.eq.-1) goto 100 * * --- degenerate case 47 --- * isign=-l532(2) if(isign.eq.1.or.isign.eq.-1) goto 100 * * --- degenerate case 48 --- * isign=-l541(2) if(isign.eq.1.or.isign.eq.-1) goto 100 * * --- degenerate case 49 --- * c write(6,600) isign c 600 format('jor5: sign cannot be determined yet; isign=',i2) isign=1 * 100 continue jor5=isign*ichang return end * * function jor3(l11,l12,l13,l21,l22,l23,l31,l32,l33) * ------------------------------------------------------------ * * * * compute the sign of 3 by 3 matrix with long integers * * * * ------------------------------------------------------------ * dimension l11(6),l12(6),l13(6) dimension l21(6),l22(6),l23(6) dimension l31(6),l32(6),l33(6) dimension la(6),lb(6),lc(6),ldeta(6),ldetb(6) * call exmult(l22,l33,la) call exmult(l32,l23,lb) call exsub(la,lb,lc) call exmult(l11,lc,ldeta) * call exmult(l12,l33,la) call exmult(l32,l13,lb) call exsub(la,lb,lc) call exmult(l21,lc,la) call exsub(ldeta,la,ldetb) * call exmult(l12,l23,la) call exmult(l22,l13,lb) call exsub(la,lb,lc) call exmult(l31,lc,la) call exadd(ldetb,la,ldeta) * jor3=ldeta(2) return end * * function jor2(l11,l12,l21,l22) * ------------------------------------------------------------ * * * * compute the sign of 2 by 2 matrix with long integers * * * * ------------------------------------------------------------ * dimension l11(6),l12(6),l21(6),l22(6) dimension la(6),lb(6),ldeta(6) * call exmult(l11,l22,la) call exmult(l21,l12,lb) call exsub(la,lb,ldeta) jor2=ldeta(2) return end * * SUBROUTINE SORTBN(NOFP,ILIST,NAME,JLIST,ICHANG) * ------------------------------------------------------------ * * * * Sort ILIST with the key NAME, and store the result into * * JLIST. Also set * * * * ICHANG = 1 if JLIST is an even permutation of ILIST * * = -1 otherwise * * * * ------------------------------------------------------------ * DIMENSION ILIST(1),JLIST(1),NAME(1) DO 10 I=1,NOFP JLIST(I)=ILIST(I) 10 CONTINUE ICHANG=1 DO 20 II=NOFP-1,1,-1 DO 21 JJ=1,II IF(NAME(JLIST(JJ)).GT.NAME(JLIST(JJ+1))) THEN JTEMP=JLIST(JJ) JLIST(JJ)=JLIST(JJ+1) JLIST(JJ+1)=JTEMP ICHANG=(-1)*ICHANG ENDIF 21 CONTINUE 20 CONTINUE RETURN END * * function jacce5(jv5) * ------------------------------------------------------------ * * * * compute the orientation in floating-point arithmetic * * * * ------------------------------------------------------------ * common /ch40/ kpara(30) common /ch41a/ kx(1024),ky(1024),kz(1024),ku(1024) common /ch41b/ kvt(1024),kname(1024),mol(1024) common /ch42/ kesv(20480),keev(20480),kecol(20480) common /ch43a/ kfv(3,20480),kft(2,20480),kfnt(20480) common /ch43b/ idface(20480) common /ch44/ ktv(4,10240),ktf(4,10240),ktcol(10240) common /ch4f1/ fltonv(333) common /ch4f2/ fltstk(10240),fltdel(10240) common /ch4f3/ flvert(1024) common /ch4w1/ istack(333),itlist(333),idlist(333) common /ch4w2/ islist(333),idflst(333) common /ch4w3/ ivlist(333) logical fltonv,fltstk,fltdel,flvert * -------------------------------------------------------------- * * zmax maximum absolute value of the terms * amax maximum absolute value of the pivots * * -------------------------------------------------------------- dimension jv5(5) double precision a(4,4),w,ww,copy data delta0/0.00000000001/ do 10 i=1,4 a(i,1)=float(kx(jv5(i+1))-kx(jv5(1))) a(i,2)=float(ky(jv5(i+1))-ky(jv5(1))) a(i,3)=float(kz(jv5(i+1))-kz(jv5(1))) a(i,4)=float(ku(jv5(i+1))-ku(jv5(1))) 10 continue delta=delta0 jsign=1 zmax=1.0 * --- execute the Gauss elimination --- do 20 i=1,3 c write(6,730) i c write(7,730) i c 730 format('i=',i5) * --- find the maximum absolute value in the i-th column --- amax=abs(a(i,i)) mmax=i do 21 m=i+1,4 if(abs(a(m,i)).gt.amax) then amax=abs(a(m,i)) mmax=m endif 21 continue if(mmax.eq.i) goto 23 * --- change the pivot --- jsign=-jsign do 22 j=1,4 copy=a(i,j) a(i,j)=a(mmax,j) a(mmax,j)=copy 22 continue 23 continue * if(abs(a(i,i)).gt.zmax) zmax=abs(a(i,i)) if(abs(a(i,i)).lt.zmax*delta) then jacce5=0 goto 100 elseif(a(i,i).lt.0.0) then jsign=-jsign endif * do 25 j=i+1,4 w=a(j,i)/a(i,i) do 25 k=i+1,4 ww=w*a(i,k) if(abs(ww).gt.zmax) zmax=abs(ww) if(abs(a(j,k)).gt.zmax) zmax=abs(a(j,k)) a(j,k)=a(j,k)-ww delta=4.0*delta 25 continue 20 continue if(abs(a(4,4)).gt.zmax) zmax=abs(a(4,4)) if(abs(a(4,4)).lt.delta*zmax) then jacce5=0 goto 100 elseif(a(4,4).lt.0.0) then jsign=-jsign endif jacce5=jsign 100 continue return end * * subroutine cansee(ivnew,itet,iyes) * ------------------------------------------------------------ * * * * check whether itet is visible from ivnew * * * * ------------------------------------------------------------ * common /ch40/ kpara(30) common /ch41a/ kx(1024),ky(1024),kz(1024),ku(1024) common /ch41b/ kvt(1024),kname(1024),mol(1024) common /ch42/ kesv(20480),keev(20480),kecol(20480) common /ch43a/ kfv(3,20480),kft(2,20480),kfnt(20480) common /ch43b/ idface(20480) common /ch44/ ktv(4,10240),ktf(4,10240),ktcol(10240) common /ch4f1/ fltonv(333) common /ch4f2/ fltstk(10240),fltdel(10240) common /ch4f3/ flvert(1024) common /ch4w1/ istack(333),itlist(333),idlist(333) common /ch4w2/ islist(333),idflst(333) logical fltonv,fltstk,fltdel,flvert * dimension iv5(5) do 10 i=1,4 iv5(i)=ktv(i,itet) 10 continue iv5(5)=ivnew iori=jor5(iv5) if(iori.eq.1) then iyes=1 else iyes=0 endif return end * * subroutine empchk * ------------------------------------------------------------ * * * * discriminate between active elements and inactive ones * * * * ------------------------------------------------------------ * common /ch40/ kpara(30) common /ch41a/ kx(1024),ky(1024),kz(1024),ku(1024) common /ch41b/ kvt(1024),kname(1024),mol(1024) common /ch41c/ kvact(1024) common /ch42/ kesv(20480),keev(20480),kecol(20480) common /ch43a/ kfv(3,20480),kft(2,20480),kfnt(20480) common /ch43b/ idface(20480) common /ch44/ ktv(4,10240),ktf(4,10240),ktcol(10240) common /ch4f1/ fltonv(333) common /ch4f2/ fltstk(10240),fltdel(10240) common /ch4f3/ flvert(1024) common /ch4w1/ istack(333),itlist(333),idlist(333) common /ch4w2/ islist(333),idflst(333) common /ch4w3/ ivlist(333) logical fltonv,fltstk,fltdel,flvert logical kvact * -------------------------------------------------------------- * kvact(i) = 0 if the i-th vertex is not on the boundary * of the four-dimensional convex hull * = 1 if it is on the boundary * kfv(2,j) = 0 if the j-th face is not used * = 1 if the j-th face is used * ktv(2,k) = 0 if the k-th tetrahedron is not used * = 1 if the k-th tetrahedron is used * -------------------------------------------------------------- * * --- classify tetra into active and inactive --- * ip=kpara(14) 10 continue if(ip.eq.0) goto 11 if(ip.gt.kpara(24)) goto 11 ktv(2,ip)=0 ip=ktv(1,ip) goto 10 * * --- classify faces into active and inactive --- * 11 continue ip=kpara(13) 15 continue if(ip.eq.0) goto 20 if(ip.gt.kpara(23)) goto 20 kfv(2,ip)=0 ip=kfv(1,ip) goto 15 * * --- classify vertices into active and inactive --- * 20 continue nofv=kpara(7) do 21 iv=1,nofv kvact(iv)=.true. 21 continue if(nofv.eq.kpara(1)) goto 30 do 22 iv=nofv+1,kpara(1) kvact(iv)=.false. 22 continue * 30 continue do 31 iv=1,nofv it=kvt(iv) if(ktv(2,it).eq.0) then kvact(iv)=.false. goto 31 endif do 32 jj=1,4 if(ktv(jj,it).eq.iv) goto 31 32 continue kvact(iv)=.false. 31 continue * return end * * subroutine makeed * ------------------------------------------------------------ * * * * make the edges * * * * ------------------------------------------------------------ * common /ch40/ kpara(30) common /ch41a/ kx(1024),ky(1024),kz(1024),ku(1024) common /ch41b/ kvt(1024),kname(1024),mol(1024) common /ch41c/ kvact(1024) common /ch42/ kesv(20480),keev(20480),kecol(20480) common /ch43a/ kfv(3,20480),kft(2,20480),kfnt(20480) common /ch43b/ idface(20480) common /ch44/ ktv(4,10240),ktf(4,10240),ktcol(10240) common /ch4f1/ fltonv(333) common /ch4f2/ fltstk(10240),fltdel(10240) common /ch4f3/ flvert(1024) common /ch4w1/ istack(333),itlist(333),idlist(333) common /ch4w2/ islist(333),idflst(333) common /ch4w3/ ivlist(333) logical fltonv,fltstk,fltdel,flvert logical kvact * kpmax=kpara(5) nofv=kpara(7) do 10 i=1,nofv flvert(i)=.false. 10 continue do 20 iv=1,nofv * * --- check if iv is active --- if(kvact(iv).eqv..false.) goto 20 * ipoint=0 call tonv(iv,ntetra) c write(7,781) iv,ntetra,(itlist(ijk),ijk=1,ntetra) c 781 format('iv,ntetra,itl=',10i3) do 30 jj=1,ntetra it=itlist(jj) do 31 kk=1,4 jv=ktv(kk,it) if(jv.gt.iv.and.(flvert(jv).eqv..false.)) then if(ipoint.ge.kpmax) then call compr(ipoint,ivlist) endif ipoint=ipoint+1 if(ipoint.gt.kpmax) goto 700 ivlist(ipoint)=jv flvert(jv)=.true. c write(7,790) iv,ipoint,jv c 790 format('makeed; iv,ipoint,jv=',3i6) endif 31 continue 30 continue call compr(ipoint,ivlist) c write(7,780) ipoint,(ivlist(ijk),ijk=1,ipoint) c 780 format('ip,ivl=',10i3) if(ipoint.eq.0) goto 20 * --- generate edges from iv to its neighbors --- do 32 jj=1,ipoint iedge=knewe(idummy) kesv(iedge)=iv keev(iedge)=ivlist(jj) * --- clear the flag flvert --- flvert(ivlist(jj))=.false. 32 continue 20 continue return * 700 continue write(6,710) ipoint 710 format('makeed: ipoint overflows, ipoint=',i6) stop end * * subroutine compr(ipoint,ivlist) * ------------------------------------------------------------ * * * * delete duplicated vertices * * * * ------------------------------------------------------------ * dimension ivlist(1) * if(ipoint.eq.0) goto 100 * * --- construct the heap structure in ivlist --- * DO 10 I=1,ipoint * * --- add a new element to the heap --- * J1=I 11 CONTINUE J2=J1/2 IF(J2.EQ.0) GOTO 10 IF(ivlist(J1).le.ivlist(J2)) GOTO 10 MV=ivlist(J1) ivlist(J1)=ivlist(J2) ivlist(J2)=MV J1=J2 GOTO 11 10 CONTINUE * * --- construct the sorted list in ivlist --- * if(ipoint.eq.1) goto 31 DO 20 I=ipoint,2,-1 * * --- extract the largest element from the heap --- * MV=ivlist(I) ivlist(I)=ivlist(1) ivlist(1)=MV * * --- reconstruct the heap for the remaining elements --- * J1=1 21 CONTINUE J2=2*J1 j3=J2+1 IF(J2.GE.I) GOTO 20 IF(J3.GE.I) THEN JLARGE=J2 ELSE IF(ivlist(J2).gt.ivlist(J3)) THEN JLARGE=J2 ELSE JLARGE=J3 ENDIF ENDIF * IF(ivlist(J1).gt.ivlist(JLARGE)) GOTO 20 MV=ivlist(J1) ivlist(J1)=ivlist(JLARGE) ivlist(JLARGE)=MV J1=JLARGE GOTO 21 * 20 CONTINUE * * --- delete duplicated elements --- * 31 continue jpoint=1 ivold=ivlist(1) do 30 i=2,ipoint if(ivlist(i).ne.ivold) then jpoint=jpoint+1 ivlist(jpoint)=ivlist(i) ivold=ivlist(jpoint) endif 30 continue ipoint=jpoint * 100 continue RETURN END * * subroutine tone(ie,ntetra) * ------------------------------------------------------------ * * * * list up the tetrahedra incident to edge ie * * * * ------------------------------------------------------------ * common /ch40/ kpara(30) common /ch41a/ kx(1024),ky(1024),kz(1024),ku(1024) common /ch41b/ kvt(1024),kname(1024),mol(1024) common /ch42/ kesv(20480),keev(20480),kecol(20480) common /ch43a/ kfv(3,20480),kft(2,20480),kfnt(20480) common /ch43b/ idface(20480) common /ch44/ ktv(4,10240),ktf(4,10240),ktcol(10240) common /ch4f1/ fltonv(333) common /ch4f2/ fltstk(10240),fltdel(10240) common /ch4f3/ flvert(1024) common /ch4w1/ istack(333),itlist(333),idlist(333) common /ch4w2/ islist(333),idflst(333) common /ch4w3/ ivlist(333) logical fltonv,fltstk,fltdel,flvert * c write(7,610) c 610 format('enter tone') kpmax=kpara(5) * * --- find an initial tetrahedron --- * ivs=kesv(ie) ive=keev(ie) call tonv(ivs,ntetra) c write(7,689) ie,ivs,ive,ntetra c 689 format('ie,ivs,ive,ntetra=',4i5) c do 888 i=1,ntetra c write(7,688) i,itlist(i) c 688 format('tonv,i,itet=',2i5) c 888 continue do 10 ii=1,ntetra itet=itlist(ii) do 11 jj=1,4 if(ktv(jj,itet).eq.ivs) goto 12 11 continue goto 10 12 continue do 13 jj=1,4 if(ktv(jj,itet).eq.ive) goto 20 13 continue 10 continue goto 700 * * --- initial tetrahedron is found --- * 20 continue itets=itet ntetra=1 itlist(ntetra)=itets c write(7,782) ntetra,itets c 782 format('tone: ntetra,jtet=',2i4) * * --- find an initial face --- * do 21 ii=1,4 iface=ktf(ii,itets) do 22 jj=1,3 if(kfv(jj,iface).eq.ivs) goto 23 22 continue goto 21 23 continue do 24 jj=1,3 if(kfv(jj,iface).eq.ive) goto 30 24 continue 21 continue goto 701 * * --- initail face is found --- * 30 continue jtet=kft(1,iface) if(jtet.eq.itet) jtet=kft(2,iface) if(jtet.eq.itets) goto 100 * * --- next tetra is found --- * ntetra=ntetra+1 if(ntetra.gt.kpmax) goto 702 itlist(ntetra)=jtet c write(7,781) ntetra,jtet,iface c 781 format('tone: ntetra,jtet,iface=',3i4) * * --- next face --- * itet=jtet do 31 jj=1,4 jface=ktf(jj,itet) if(jface.eq.iface) goto 31 do 32 jk=1,3 if(kfv(jk,jface).eq.ivs) goto 33 32 continue goto 31 33 continue do 34 jk=1,3 if(kfv(jk,jface).eq.ive) goto 40 34 continue 31 continue goto 701 40 continue iface=jface goto 30 * 100 continue return * 700 continue write(6,710) ie,ivs,ive 710 format('tone: initial tetra cannot be found'/ +' ie,ivs,ive=',3i10) goto 777 701 continue write(6,711) 711 format('tone: next tetra cannot be found') goto 777 702 continue write(6,712) ntetra 712 format('tone: list area overflows, ntetra=',i10) 777 continue stop end * * subroutine etcolo * ------------------------------------------------------------ * * * * classify the tetrahedra and edges * * * * kecol 1: edge of the Delaunay tetrahedraizaion * * 2: edge of the farthest-point version * * 3: edge on both the tetrahedrizaions * * ktcol 1: Delaunay tetrahedron * * 2: farthest-point Delaunay tetrahedron * * 0: unused cell * * ------------------------------------------------------------ * common /ch40/ kpara(30) common /ch41a/ kx(1024),ky(1024),kz(1024),ku(1024) common /ch41b/ kvt(1024),kname(1024),mol(1024) common /ch42/ kesv(20480),keev(20480),kecol(20480) common /ch43a/ kfv(3,20480),kft(2,20480),kfnt(20480) common /ch43b/ idface(20480) common /ch44/ ktv(4,10240),ktf(4,10240),ktcol(10240) common /ch4f1/ fltonv(333) common /ch4f2/ fltstk(10240),fltdel(10240) common /ch4f3/ flvert(1024) common /ch4w1/ istack(333),itlist(333),idlist(333) common /ch4w2/ islist(333),idflst(333) common /ch4w3/ ivlist(333) logical fltonv,fltstk,fltdel,flvert * dimension iv5(5) ivmax=kpara(1) ixymax=kpara(6) * * --- select active tetrahedra --- * ktcol=1: active, =0: inactive * do 10 itet=1,kpara(24) ktcol(itet)=1 10 continue ip=kpara(14) 11 continue if(ip.eq.0) goto 20 if(ip.gt.kpara(24)) goto 20 ktcol(ip)=0 ip=ktv(1,ip) goto 11 * * --- classify the tetrahedra --- * * This is dangerous because address ivmax is used for * the test of the orientation; so essential ivamx should * ivmax-1. * 20 continue minf=-(ixymax+1) iv5(5)=ivmax do 21 itet=1,kpara(24) if(ktcol(itet).eq.0) goto 21 do 22 ii=1,4 iv5(ii)=ktv(ii,itet) 22 continue c write(6,676) itet,(iv5(jjj),jjj=1,5) c 676 format('itet,iv5=',6i7) do 25 ippp=1,4 kx(ivmax)=kx(iv5(ippp)) ky(ivmax)=ky(iv5(ippp)) kz(ivmax)=kz(iv5(ippp)) ku(ivmax)=minf iori=jor5(iv5) c do 26 jjj=1,5 c write(6,677) kx(iv5(jjj)),ky(iv5(jjj)),kz(iv5(jjj)), c + ku(iv5(jjj)) c 677 format(4i16) c 26 continue c write(6,678) iori c 678 format('iori = ',i5) 25 continue if(iori.eq.1) then ktcol(itet)=1 else ktcol(itet)=2 endif c write(7,780) itet,ktcol(itet) c 780 format('etcolo: itet,tcol=',2i4) 21 continue * * --- classify the edges --- * do 30 ie=1,kpara(22) call tone(ie,ntetra) icolor=ktcol(itlist(1)) do 31 ii=2,ntetra itet=itlist(ii) if(ktcol(itet).ne.icolor) then kecol(ie)=3 goto 30 endif 31 continue if(icolor.eq.1) then kecol(ie)=1 else kecol(ie)=2 endif 30 continue * return end * * subroutine monit * ------------------------------------------------------------ * * * * dump the data structure * * * * ------------------------------------------------------------ * common /ch40/ kpara(30) common /ch41a/ kx(1024),ky(1024),kz(1024),ku(1024) common /ch41b/ kvt(1024),kname(1024),mol(1024) common /ch41c/ kvact(1024) common /ch42/ kesv(20480),keev(20480),kecol(20480) common /ch43a/ kfv(3,20480),kft(2,20480),kfnt(20480) common /ch43b/ idface(20480) common /ch44/ ktv(4,10240),ktf(4,10240),ktcol(10240) common /ch4f1/ fltonv(333) common /ch4f2/ fltstk(10240),fltdel(10240) common /ch4f3/ flvert(1024) common /ch4w1/ istack(333),itlist(333),idlist(333) common /ch4w2/ islist(333),idflst(333) common /ch4w3/ ivlist(333) logical fltonv,fltstk,fltdel,flvert logical kvact * * --- set the basic parameters for the data structure --- * write(7,700)(kpara(ii),ii=1,8),(kpara(ii),ii=12,14), + (kpara(ii),ii=22,24) 700 format('monitor the data structure'/ + 'kpara(1,2,3,4)=',4i12/ + 'kpara(5,6,7,8)=',4i12/ + 'kpara(12,13,14)=',3i12/ + 'kpara(22,23,24)=',3i12) * write(7,701) 701 format('--- vertices ---') do 10 ii=1,kpara(7) write(7,702) ii,kvact(ii),kvt(ii), + kx(ii),ky(ii),kz(ii),ku(ii) 702 format('i,kvact,kvt,x=',i4,i2,i4,4i12) 10 continue write(7,703) 703 format('--- edges ---') do 20 ii=1,kpara(22) write(7,704) ii,kesv(ii),keev(ii),kecol(ii) 704 format('i,sv,ev,col=',4i5) 20 continue write(7,705) 705 format('--- faces ---') do 30 i=1,kpara(23) write(7,706) i,(kfv(j,i),j=1,3),(kft(k,i),k=1,2) 706 format('i,vvv,tt=',6i5) 30 continue write(7,707) 707 format('--- tetrahedra ---') do 40 i=1,kpara(24) write(7,708) i,(ktv(j,i),j=1,4),(ktf(j,i),j=1,4), + ktcol(i) 708 format('i,vvvv,ffff,col=',10i4) 40 continue return end * * subroutine statis * ------------------------------------------------------------ * * * * print the statistics data * * * * ------------------------------------------------------------ * common /ch40/ kpara(30) common /ch41a/ kx(1024),ky(1024),kz(1024),ku(1024) common /ch41b/ kvt(1024),kname(1024),mol(1024) common /ch42/ kesv(20480),keev(20480),kecol(20480) common /ch43a/ kfv(3,20480),kft(2,20480),kfnt(20480) common /ch43b/ idface(20480) common /ch44/ ktv(4,10240),ktf(4,10240),ktcol(10240) common /ch4f1/ fltonv(333) common /ch4f2/ fltstk(10240),fltdel(10240) common /ch4f3/ flvert(1024) common /ch4w1/ istack(333),itlist(333),idlist(333) common /ch4w2/ islist(333),idflst(333) common /ch4w3/ ivlist(333) logical fltonv,fltstk,fltdel,flvert * * --- set the basic parameters for the data structure --- * write(7,700) 700 format(/'--- statistics ---') write(7,701)(kpara(ii),ii=1,8),(kpara(ii),ii=12,14), + (kpara(ii),ii=22,24) 701 format('monitor the data structure'/ + 'kpara(1,2,3,4)=',4i12/ + 'kpara(5,6,7,8)=',4i12/ + 'kpara(12,13,14)=',3i12/ + 'kpara(22,23,24)=',3i12) write(7,702) kpara(22) 702 format('no. of edges =',i10) * --- count the number of active faces --- noffac=kpara(23) ip=kpara(13) 10 continue if(ip.eq.0) goto 11 if(ip.le.kpara(23)) then noffac=noffac-1 ip=kfv(1,ip) goto 10 endif 11 continue write(7,703) noffac 703 format('no. of faces =',i10) * --- count the number of active tetrahedra --- noftet=kpara(24) ip=kpara(14) 20 continue if(ip.eq.0) goto 21 if(ip.le.kpara(24)) then noftet=noftet-1 ip=ktv(1,ip) goto 20 endif 21 continue write(7,704) noftet 704 format('no. of tetrahedra =',i10) return end SUBROUTINE EXADD(LONGA,LONGB,LONGC) * ------------------------------------------------------------ * * * * Exact addition of variable-length integers * * * * LONGC = LONGA + LONGB * * * * LONG(1) length of the array LONG for a long integer * * LONG(2) sign of the integer represented in LONG * * LONG(K) K=3,4,...,LONG(1); absolute values divided * * into thirty-bit integers * * * * ------------------------------------------------------------ * DIMENSION LONGA(1),LONGB(1),LONGC(1) DATA LBASE/1073741824/ * -------------------------------------------------------------- * * LBASE=2**30 * * -------------------------------------------------------------- * N=LONGA(1) IF(LONGB(1).NE.N) GOTO 700 LONGC(1)=N * * --- check whether LONGA or LONGB is 0 --- * IF(LONGA(2).EQ.0) THEN DO 10 I=2,N LONGC(I)=LONGB(I) 10 CONTINUE GOTO 100 ELSEIF(LONGB(2).EQ.0) THEN DO 20 I=2,N LONGC(I)=LONGA(I) 20 CONTINUE GOTO 100 ENDIF * * --- check whether LONGA and LONGB are of the same sign --- * IF(LONGA(2)*LONGB(2).EQ.-1) GOTO 50 * * --- two numbers with the same sign --- * LONGC(2)=LONGA(2) ICARRY=0 DO 30 II=3,N IVAL=LONGA(II)+LONGB(II)+ICARRY ICARRY=IVAL/LBASE LONGC(II)=MOD(IVAL,LBASE) 30 CONTINUE GOTO 100 * * --- two numbers with different signs --- * 50 CONTINUE * * --- compare the absolute values --- * DO 51 II=N,3,-1 IF(LONGA(II).GT.LONGB(II)) GOTO 60 IF(LONGB(II).GT.LONGA(II)) GOTO 70 51 CONTINUE * * --- The absolute values are the same --- * LONGC(2)=0 DO 52 II=3,N LONGC(II)=0 52 CONTINUE GOTO 100 * * --- The absolute value of LONGA is larger --- * --- than the absolute value of LONGB --- * 60 CONTINUE LONGC(2)=LONGA(2) ICARRY=0 DO 61 II=3,N KVAL=LONGA(II)-LONGB(II)-ICARRY IF(KVAL.LT.0) THEN ICARRY=1 KVAL=KVAL+LBASE ELSE ICARRY=0 ENDIF LONGC(II)=KVAL 61 CONTINUE GOTO 100 * * --- The absolute value of LONGB is larger --- * --- than the absolute value of LONGA --- * 70 CONTINUE LONGC(2)=LONGB(2) ICARRY=0 DO 71 II=3,N KVAL=LONGB(II)-LONGA(II)-ICARRY IF(KVAL.LT.0) THEN ICARRY=1 KVAL=KVAL+LBASE ELSE ICARRY=0 ENDIF LONGC(II)=KVAL 71 CONTINUE GOTO 100 * 100 CONTINUE RETURN * * --- error message --- * 700 CONTINUE WRITE(6,701) LONGA(1),LONGB(1) 701 FORMAT(' Error in EXADD: The integer sizes do not match.'/ + ' LONGA: ',I3,' words, LONGB: ',I3,' words') STOP * END * * SUBROUTINE EXSUB(LONGA,LONGB,LONGC) * ------------------------------------------------------------ * * * * Exact subtraction of variable-length integers * * * * LONGC = LONGA - LONGB * * * * -------------------------------------------------------------* DIMENSION LONGA(1),LONGB(1),LONGC(1) * IF(LONGB(2).EQ.0) THEN DO 10 I=1,LONGA(1) LONGC(I)=LONGA(I) 10 CONTINUE ELSE LONGB(2)=(-1)*LONGB(2) CALL EXADD(LONGA,LONGB,LONGC) LONGB(2)=(-1)*LONGB(2) ENDIF RETURN END * * SUBROUTINE EXMULT(LONGA,LONGB,LONGC) * ------------------------------------------------------------ * * * * Exact multiplication of variable-length integers * * * * LONGC = LONGA * LONGB * * * * ------------------------------------------------------------ * DIMENSION LONGA(1),LONGB(1),LONGC(1) DATA LBASE/1073741824/, LHBASE/32768/ * -------------------------------------------------------------- * * LBASE=2**30, LHBASE=2**15 * * -------------------------------------------------------------- * * --- set the initial values --- * N=LONGA(1) IF(LONGB(1).NE.N) GOTO 700 LONGC(1)=N M=(N-2)*2 LONGC(2)=LONGA(2)*LONGB(2) * * check whether LONGC is 0 --- * IF(LONGC(2).EQ.0) THEN DO 10 II=3,N LONGC(II)=0 10 CONTINUE GOTO 100 ENDIF * ICARRY=0 * * --- larger loop (multiplication) --- * DO 20 II=1,M * KVAL=ICARRY JCARRY=0 * * --- smaller loop (compute the II-th half-word integer) --- * DO 30 JJ=1,II * KK=II-JJ+1 KVAL=KVAL+LHALF(LONGA,JJ)*LHALF(LONGB,KK) IF(KVAL.GE.LBASE) THEN JCARRY=JCARRY+1 KVAL=KVAL-LBASE ENDIF * 30 CONTINUE * ICARRY=JCARRY*LHBASE+KVAL/LHBASE KVAL=MOD(KVAL,LHBASE) IF(MOD(II,2).EQ.1) THEN LONGC((II+1)/2+2)=KVAL ELSE LONGC(II/2+2)=LONGC(II/2+2)+KVAL*LHBASE ENDIF * 20 CONTINUE * 100 CONTINUE RETURN * * --- error message --- * 700 CONTINUE WRITE(6,701) LONGA(1),LONGB(1) 701 FORMAT(' Error in EXMULT:', + ' The integer sizes do not match.'/ + ' LONGA: ',I3,' words, LONGB: ',I3,' words') STOP * END * * FUNCTION LHALF(LONG,I) * -------------------------------------------------------------- * * * * Extract the I-th half-word integer of the * * variable-length integer LONG * * * * -------------------------------------------------------------- * DIMENSION LONG(1) DATA LHBASE/32768/ * M=(I+1)/2 IF(M.LE.0.OR.M.GT.LONG(1)-2) GOTO 700 IF(MOD(I,2).EQ.1) THEN LHALF=MOD(LONG(M+2),LHBASE) ELSE LHALF=LONG(M+2)/LHBASE ENDIF * 100 CONTINUE RETURN * * --- error message --- * 700 CONTINUE WRITE(6,701) I 701 FORMAT(' Error in LHALF: Integer I is out of range.', + ' I = ',I5) STOP END * * FUNCTION FLOATL(LONG) * ------------------------------------------------------------ * * * * convert a variable-length integer into a floating- * * point number * * * * ------------------------------------------------------------ * DIMENSION LONG(1) DATA BASE/1073741824.0/ * VAL=0.0 DO 10 I=3,LONG(1) VAL=VAL+FLOAT(LONG(I))*(BASE**(I-3)) 10 CONTINUE FLOATL=VAL*LONG(2) RETURN END * * SUBROUTINE LTRANS(IA,LONG,N) * ------------------------------------------------------------ * * * * Convert a single-precision integer IA into * * variable-length integer LONG, which uses N * * words to represent an integer * * * * ------------------------------------------------------------ * DIMENSION LONG(1) DATA LBASE/1073741824/ * -------------------------------------------------------------- * * LBASE = 2**30 * * -------------------------------------------------------------- * * --- store the length of the array LONG --- * IF(N.LE.3) GOTO 700 LONG(1)=N * * --- store the sign of the integer --- * IF(IA.GT.0) THEN LONG(2)=1 ELSEIF(IA.EQ.0) THEN LONG(2)=0 ELSE LONG(2)=-1 ENDIF * * --- store the absolute value of the integer --- * IB=IABS(IA) IC=IB/LBASE IF(IC.EQ.1) IB=IB-LBASE LONG(3)=IB LONG(4)=IC IF(N.LE.4) GOTO 100 DO 10 I=5,N LONG(I)=0 10 CONTINUE * 100 CONTINUE RETURN 700 CONTINUE WRITE(6,701) N 701 FORMAT('Error in LTRANS: N should be greater than 3.'/ + 'Present N is',I3) GOTO 100 END * * SUBROUTINE LTRANF(A,LONG,N) * ------------------------------------------------------------ * * * * translate a floating-point number into a varialbe- * * length integer * * * * ------------------------------------------------------------ * DIMENSION LONG(1) DATA BASE/1073741824.0/ * * --- set the length of the array --- * LONG(1)=N * * --- set the sign LONG(2) and the absolute value B of A --- * IF(A.GT.0.0) THEN LONG(2)=1 B=A ELSEIF(A.LT.0.0) THEN LONG(2)=-1 B=-A ELSE LONG(2)=0 DO 10 I=3,N LONG(I)=0 10 CONTINUE GOTO 100 ENDIF * * --- compute the minimum size M of the array LONG to * represent the integer part of the value B --- * C=B M=3 20 CONTINUE IF(C.LT.BASE) GOTO 30 C=C/BASE M=M+1 GOTO 20 * * --- put the absolute value B into the array LONG --- * 30 CONTINUE IF(M.GT.N) GOTO 700 LONG(M)=IFIX(C) IF(M.EQ.3) GOTO 40 DO 31 I=M-1,3,-1 B=B-FLOAT(IFIX(C))*(BASE**(I-2)) C=B/(BASE**(I-3)) LONG(I)=IFIX(C) 31 CONTINUE * 40 CONTINUE IF(M.EQ.N) GOTO 100 DO 41 I=M+1,N LONG(I)=0 41 CONTINUE * 100 CONTINUE RETURN * 700 CONTINUE WRITE(6,710) A,N 710 FORMAT('Error in LTRANF: A is too large.'/ + ' A, N =',E30.15,I3) STOP END SUBROUTINE HPSORT(NOFP,IX,NAME,MOL) * ---------------------------------------------------------------- * * * * Heap Sort * * * * ---------------------------------------------------------------- * DIMENSION IX(1),NAME(1),MOL(1) * ------------------------------------------------------------------ * * NOFP number of points to be sorted * IX key for the sorting * NAME ordinal number (name) assigned to the points * MOL sorted list (this area is also used tentatively * to represent the heap structure) * * ------------------------------------------------------------------ * * --- construct the heap structure in MOL --- * DO 10 I=1,NOFP * * --- add a new element to the heap --- * MOL(I)=I J1=I 11 CONTINUE J2=J1/2 IF(J2.EQ.0) GOTO 10 IF(LARGER(MOL(J1),MOL(J2),IX,NAME).EQ.0) GOTO 10 MV=MOL(J1) MOL(J1)=MOL(J2) MOL(J2)=MV J1=J2 GOTO 11 10 CONTINUE * * --- construct the sorted list in MOL --- * DO 20 I=NOFP,2,-1 * * --- extract the largest element from the heap --- * MV=MOL(I) MOL(I)=MOL(1) MOL(1)=MV * * --- reconstruct the heap for the remaining elements --- * J1=1 21 CONTINUE J2=2*J1 j3=J2+1 IF(J2.GE.I) GOTO 20 IF(J3.GE.I) THEN JLARGE=J2 ELSE IF(LARGER(MOL(J2),MOL(J3),IX,NAME).EQ.1) THEN JLARGE=J2 ELSE JLARGE=J3 ENDIF ENDIF * IF(LARGER(MOL(J1),MOL(JLARGE),IX,NAME).EQ.1) GOTO 20 MV=MOL(J1) MOL(J1)=MOL(JLARGE) MOL(JLARGE)=MV J1=JLARGE GOTO 21 * 20 CONTINUE * RETURN END * * FUNCTION LARGER(IV,JV,IX,NAME) * ------------------------------------------------------------ * * * * Check whether the x coordinate of vertex IV is larger * * than that of JV in the sense of symbolic perturbation * * * * ------------------------------------------------------------ * DIMENSION IX(1),NAME(1) * IF(IX(IV).GT.IX(JV)) THEN LARGER=1 ELSEIF(IX(IV).LT.IX(JV)) THEN LARGER=0 ELSE IF(NAME(IV).GT.NAME(JV)) THEN LARGER=0 ELSE LARGER=1 ENDIF ENDIF * RETURN END SUBROUTINE DRAWC2(NOFP,IX,IY,IXYMAX,NOFV,NVLIST) * ---------------------------------------------------------------- * * * * Draw the two-dimensional convex hull * * (construct the ps file 'draw.ps') * * * * ---------------------------------------------------------------- * DIMENSION IX(1),IY(1),NVLIST(1) DATA RADIUS/1.0/ OPEN(1,FILE='draw.ps') * CALL DRSTAR(-IXYMAX,IXYMAX,-IXYMAX,IXYMAX) * * --- plot the input points --- * DO 10 II=1,NOFP CALL DRDOT(IX(II),IY(II),RADIUS) 10 CONTINUE * * --- draw the boundary of the convex hull --- * IV=NVLIST(1) IX1=IX(IV) IY1=IY(IV) DO 20 II=2,NOFV IV=NVLIST(II) IX2=IX(IV) IY2=IY(IV) CALL DRLINE(IX1,IY1,IX2,IY2) IX1=IX2 IY1=IY2 20 CONTINUE IV=NVLIST(1) IX2=IX(IV) IY2=IY(IV) CALL DRLINE(IX1,IY1,IX2,IY2) * CALL DRCOMP * CLOSE(1) RETURN END * * SUBROUTINE DRSTAR(IXMIN,IXMAX,IYMIN,IYMAX) * ---------------------------------------------------------------- * * * * Set the initial state of the ps file * * * * ---------------------------------------------------------------- * COMMON /DRAWPA/SCALE DATA IWIDTH/450/,XSHIFT/50.0/,YSHIFT/150.0/ DATA LWIDTH/0.4/ IXS=IXMAX-IXMIN IYS=IYMAX-IYMIN IF(IXS.GE.IYS) THEN MS=IXS ELSE MS=IYS ENDIF SCALE=FLOAT(IWIDTH)/FLOAT(MS) XMIN=FLOAT(IXMIN)*SCALE YMIN=FLOAT(IYMIN)*SCALE WRITE(1,100) -XMIN+XSHIFT,-YMIN+YSHIFT 100 FORMAT(2F10.3,' translate') WRITE(1,101) LWIDTH 101 FORMAT(I4,' setlinewidth') RETURN END * * SUBROUTINE LWIDTH(WIDTH) * ---------------------------------------------------------------- * * * * set the line width * * * * ---------------------------------------------------------------- * WRITE(1,101) WIDTH 101 FORMAT(F8.4,' setlinewidth') RETURN END * * SUBROUTINE DRLINE(JX1,JY1,JX2,JY2) * ---------------------------------------------------------------- * * * * Draw a line from (JX1,JY1) to (JX2,JY2) * * * * ---------------------------------------------------------------- * COMMON /DRAWPA/SCALE X=FLOAT(JX1)*SCALE Y=FLOAT(JY1)*SCALE WRITE(1,100) X,Y 100 FORMAT('newpath ',2F10.3,' moveto') X=FLOAT(JX2)*SCALE Y=FLOAT(JY2)*SCALE WRITE(1,101) X,Y 101 FORMAT(' ',2F10.3,' lineto stroke') RETURN END * * SUBROUTINE DRDOT(JX,JY,R) * ---------------------------------------------------------------- * * * * Draw a dot with radius R centered at (JX,JY) * * * * ---------------------------------------------------------------- * COMMON /DRAWPA/SCALE X=FLOAT(JX)*SCALE Y=FLOAT(JY)*SCALE WRITE(1,100) X,Y,R 100 FORMAT('newpath ',3F10.3,' 0 360 arc closepath fill') RETURN END * * SUBROUTINE DRCIRC(JX,JY,JR) * ---------------------------------------------------------------- * * * * Draw a circle with radius JR centered at (JX,JY) * * * * ---------------------------------------------------------------- * COMMON /DRAWPA/SCALE X=FLOAT(JX)*SCALE Y=FLOAT(JY)*SCALE R=FLOAT(JR)*SCALE WRITE(1,100) X,Y,R 100 FORMAT('newpath ',3F10.3,' 0 360 arc closepath stroke') RETURN END * * SUBROUTINE DRARC(JX0,JY0,JR,A1,A2) * ---------------------------------------------------------------- * * * * Draw an arc with the center (JX0,JY0) and radius JR * * A1 and A2 are start and end angles (degrees) * * * * ---------------------------------------------------------------- * COMMON /DRAWPA/SCALE X=FLOAT(JX0)*SCALE Y=FLOAT(JY0)*SCALE R=FLOAT(JR)*SCALE WRITE(1,100) X,Y,R,A1,A2 100 FORMAT('newpath ',5F10.3,' arc stroke') RETURN END * * SUBROUTINE DRELLI(X0,Y0,RATIO,THETA,RLONG) * ---------------------------------------------------------------- * * * * Draw an ellipse * * * * ---------------------------------------------------------------- * DEG=3.1416/180.0 ST=SIN(THETA) CT=COS(THETA) A11I=RATIO*(CT**2)+ST**2 A22I=RATIO*(ST**2)+CT**2 A12I=(RATIO-1.0)*CT*ST IX1=A11I*RLONG+X0 IY1=A12I*RLONG+Y0 DO 10 II=1,360 X=RLONG*COS(DEG*II) Y=RLONG*SIN(DEG*II) IX2=A11I*X+A12I*Y+X0 IY2=A12I*X+A22I*Y+Y0 CALL DRLINE(IX1,IY1,IX2,IY2) IX1=IX2 IY1=IY2 10 CONTINUE RETURN END * * SUBROUTINE DRCOMP * ---------------------------------------------------------------- * * * * Complete the ps file * * * * ---------------------------------------------------------------- * WRITE(1,100) 100 FORMAT('showpage') RETURN END * * **************************************************************** **************************************************************** * Subroutines for making ps files * **************************************************************** **************************************************************** SUBROUTINE DRSTLE(IXMIN,IXMAX,IYMIN,IYMAX,iup) * ------------------------------------------------------------ * * * * Set the environment for drawing the left diagram * * * * ------------------------------------------------------------ * COMMON /DRAWPA/SCALE DATA IWIDTH/200/,XSHIFT/50.0/,YSHIFT/400.0/ data ydown/-250.0/,xback/-220.0/ IXS=IXMAX-IXMIN IYS=IYMAX-IYMIN IF(IXS.GE.IYS) THEN MS=IXS ELSE MS=IYS ENDIF SCALE=FLOAT(IWIDTH)/FLOAT(MS) XMIN=FLOAT(IXMIN)*SCALE YMIN=FLOAT(IYMIN)*SCALE if(iup.eq.1) then WRITE(1,100) -XMIN+XSHIFT,-YMIN+YSHIFT else WRITE(1,100) xback,ydown endif 100 FORMAT(2F10.3,' translate') RETURN END * * SUBROUTINE DRSTRI * ------------------------------------------------------------ * * * * Set the environment for drawing the left diagram * * * * ------------------------------------------------------------ * COMMON /DRAWPA/SCALE DATA IWIDTH/200/,XSHIFT/50.0/,YSHIFT/150.0/ DATA GAP/20.0/ WRITE(1,100) FLOAT(IWIDTH)+GAP,0.0 100 FORMAT(2F10.3,' translate') RETURN END * * SUBROUTINE CENTPR(EX,EY,EZ,X,Y,Z,XX,YY) * ------------------------------------------------------------ * * * * Central projection of point (x,y,z) with respect to * * the eye at (ex,ey,ez) onto the xy plane * * * * ------------------------------------------------------------ * XX=X-Z*(EX-X)/(EZ-Z) YY=Y-Z*(EY-Y)/(EZ-Z) RETURN END * * SUBROUTINE DRDELA(NOFP,IX,IY,IXYMAX,RADIUS, + KPARA,KVE,KFE,KSV,KEV,KRF,KLF,KSCE,KSCCE,KECE,KECCE, + KVCOLO,KECOLO) * ------------------------------------------------------------ * * * * Draw the Delaunay triangulation * * * * ------------------------------------------------------------ * DIMENSION IX(1),IY(1) DIMENSION KPARA(1),KVE(1),KFE(1),KSV(1),KEV(1) DIMENSION KSCE(1),KSCCE(1),KECE(1),KECCE(1) DIMENSION KVCOLO(1),KECOLO(1) * * --- plot the input points --- * DO 10 II=1,NOFP CALL DRDOT(IX(II),IY(II),RADIUS) 10 CONTINUE * * --- draw the edges --- * DO 20 II=1,KPARA(3) IF(KECOLO(II).EQ.0.OR.KECOLO(II).EQ.2) GOTO 20 IVS=KSV(II) IX1=IX(IVS) IY1=IY(IVS) IVE=KEV(II) IX2=IX(IVE) IY2=IY(IVE) CALL DRLINE(IX1,IY1,IX2,IY2) 20 CONTINUE * RETURN END * * SUBROUTINE DRFDEL(NOFP,IX,IY,IXYMAX,RADIUS, + KPARA,KVE,KFE,KSV,KEV,KRF,KLF,KSCE,KSCCE,KECE,KECCE, + KVCOLO,KECOLO) * ------------------------------------------------------------ * * * * Draw the farthest Delaunay triangulation * * * * ------------------------------------------------------------ * DIMENSION IX(1),IY(1) DIMENSION KPARA(1),KVE(1),KFE(1),KSV(1),KEV(1) DIMENSION KSCE(1),KSCCE(1),KECE(1),KECCE(1) DIMENSION KVCOLO(1),KECOLO(1) * * --- plot the input points --- * DO 10 II=1,NOFP CALL DRDOT(IX(II),IY(II),RADIUS) 10 CONTINUE * * --- draw the edges --- * DO 20 II=1,KPARA(3) IF(KECOLO(II).EQ.0.OR.KECOLO(II).EQ.1) GOTO 20 IVS=KSV(II) IX1=IX(IVS) IY1=IY(IVS) IVE=KEV(II) IX2=IX(IVE) IY2=IY(IVE) CALL DRLINE(IX1,IY1,IX2,IY2) 20 CONTINUE * RETURN END * * SUBROUTINE DRVORO(NOFP,IX,IY,IXYMAX,RADIUS, + KPARA,KVE,KFE,KSV,KEV,KRF,KLF,KSCE,KSCCE,KECE,KECCE, + XVORO,YVORO,KVCOLO,KECOLO,XEDIRE,YEDIRE,KFCOLO,KFCENT) * ------------------------------------------------------------ * * * * Draw the ordinary Voronoi diagram * * * * ------------------------------------------------------------ * DIMENSION IX(1),IY(1) DIMENSION KPARA(1),KVE(1),KFE(1) DIMENSION KSV(1),KEV(1),KRF(1),KLF(1) DIMENSION KSCE(1),KSCCE(1),KECE(1),KECCE(1) DIMENSION XVORO(1),YVORO(1) DIMENSION KVCOLO(1),KECOLO(1),KFCOLO(1),KFCENT(1) DIMENSION XEDIRE(1),YEDIRE(1) DATA SCA/2.0/ * * --- plot the input points --- * DO 10 II=1,NOFP CALL DRDOT(IX(II),IY(II),RADIUS) 10 CONTINUE * * --- draw the edges --- * DO 20 IEDGE=1,KPARA(3) IF(KECOLO(IEDGE).EQ.0.OR.KECOLO(IEDGE).EQ.2) GOTO 20 IFS=KRF(IEDGE) IFE=KLF(IEDGE) IF(KFCENT(IFS).EQ.9.AND.KFCENT(IFE).EQ.9) GOTO 20 IF(KECOLO(IEDGE).EQ.3) THEN IF(KFCOLO(IFS).EQ.1.AND.KFCENT(IFS).EQ.9) GOTO 20 IF(KFCOLO(IFE).EQ.1.AND.KFCENT(IFE).EQ.9) GOTO 20 ENDIF IF(KFCENT(IFS).EQ.9.OR.KFCENT(IFE).EQ.9) THEN * * --- almost infinite edge --- * DX=XEDIRE(IEDGE) DY=YEDIRE(IEDGE) IF(KFCENT(IFS).NE.9) THEN * * --- IFS is a finite vertex --- * JX1=IFIX(XVORO(IFS)) JY1=IFIX(YVORO(IFS)) ELSE * * --- IFE is a finite vertex --- * JX1=IFIX(XVORO(IFE)) JY1=IFIX(YVORO(IFE)) DX=-DX DY=-DY ENDIF JX2=IFIX(SCA*DX*IXYMAX)+JX1 JY2=IFIX(SCA*DY*IXYMAX)+JY1 * ELSEIF(KECOLO(IEDGE).EQ.1) THEN * * --- finite edge --- * JX1=IFIX(XVORO(IFS)) JY1=IFIX(YVORO(IFS)) JX2=IFIX(XVORO(IFE)) JY2=IFIX(YVORO(IFE)) * ELSE * * --- infinite edge --- * DX=XEDIRE(IEDGE) DY=YEDIRE(IEDGE) IF(KFCOLO(IFS).EQ.1) THEN * * --- IFS is a finite vertex --- * JX1=IFIX(XVORO(IFS)) JY1=IFIX(YVORO(IFS)) ELSE * * --- IFE is a finite vertex --- * JX1=IFIX(XVORO(IFE)) JY1=IFIX(YVORO(IFE)) DX=-DX DY=-DY ENDIF JX2=IFIX(SCA*DX*IXYMAX)+JX1 JY2=IFIX(SCA*DY*IXYMAX)+JY1 ENDIF CALL DRLINE(JX1,JY1,JX2,JY2) 20 CONTINUE * RETURN END * * SUBROUTINE DRFVOR(NOFP,IX,IY,IXYMAX,RADIUS, + KPARA,KVE,KFE,KSV,KEV,KRF,KLF,KSCE,KSCCE,KECE,KECCE, + XVORO,YVORO,KVCOLO,KECOLO,XEDIRE,YEDIRE,KFCOLO,KFCENT) * ------------------------------------------------------------ * * * * Draw the farthest-point Voronoi diagram * * * * ------------------------------------------------------------ * DIMENSION IX(1),IY(1) DIMENSION KPARA(1),KVE(1),KFE(1) DIMENSION KSV(1),KEV(1),KRF(1),KLF(1) DIMENSION KSCE(1),KSCCE(1),KECE(1),KECCE(1) DIMENSION XVORO(1),YVORO(1) DIMENSION KVCOLO(1),KECOLO(1),KFCOLO(1),KFCENT(1) DIMENSION XEDIRE(1),YEDIRE(1) DATA SCA/2.0/ * * --- plot the input points --- * DO 10 II=1,NOFP CALL DRDOT(IX(II),IY(II),RADIUS) 10 CONTINUE * * --- draw the edges --- * DO 20 IEDGE=1,KPARA(3) IF(KECOLO(IEDGE).EQ.0.OR.KECOLO(IEDGE).EQ.1) GOTO 20 IFS=KRF(IEDGE) IFE=KLF(IEDGE) IF(KFCENT(IFS).EQ.9.AND.KFCENT(IFE).EQ.9) GOTO 20 IF(KECOLO(IEDGE).EQ.3) THEN IF(KFCOLO(IFS).EQ.-1.AND.KFCENT(IFS).EQ.9) GOTO 20 IF(KFCOLO(IFE).EQ.-1.AND.KFCENT(IFE).EQ.9) GOTO 20 ENDIF * IF(KECOLO(IEDGE).EQ.3) THEN * * --- infinite edge --- * DX=XEDIRE(IEDGE) DY=YEDIRE(IEDGE) IF(KFCOLO(IFS).EQ.-1) THEN * * --- IFS is a finite vertex --- * JX1=IFIX(XVORO(IFS)) JY1=IFIX(YVORO(IFS)) ELSE * * --- IFE is a finite vertex --- * JX1=IFIX(XVORO(IFE)) JY1=IFIX(YVORO(IFE)) DX=-DX DY=-DY ENDIF JX2=IFIX(SCA*DX*IXYMAX)+JX1 JY2=IFIX(SCA*DY*IXYMAX)+JY1 * ELSEIF(KFCENT(IFS).EQ.9.OR.KFCENT(IFE).EQ.9) THEN * * --- almost infinite edge --- * DX=XEDIRE(IEDGE) DY=YEDIRE(IEDGE) IF(KFCENT(IFS).NE.9) THEN * * --- IFS is a finite vertex --- * JX1=IFIX(XVORO(IFS)) JY1=IFIX(YVORO(IFS)) ELSE * * --- IFE is a finite vertex --- * JX1=IFIX(XVORO(IFE)) JY1=IFIX(YVORO(IFE)) DX=-DX DY=-DY ENDIF JX2=IFIX(SCA*DX*IXYMAX)+JX1 JY2=IFIX(SCA*DY*IXYMAX)+JY1 * ELSE * * --- finite edge --- * JX1=IFIX(XVORO(IFS)) JY1=IFIX(YVORO(IFS)) JX2=IFIX(XVORO(IFE)) JY2=IFIX(YVORO(IFE)) ENDIF CALL DRLINE(JX1,JY1,JX2,JY2) 20 CONTINUE * RETURN END