57 INTEGER,
PARAMETER :: MNODE = 100000, &
68 INTEGER :: MXBPTS,MXCAV,MXCELL,MXEDGE,MXFACE,MXNODE,MXOCTR,MXRING,MXTEST
106 INTEGER,
INTENT(INOUT) :: nbface,nbpts,ncell,nnode
107 INTEGER,
INTENT(INOUT) :: nfcei(3,nbface),ndci(4,ncell)
108 DOUBLE PRECISION,
INTENT(INOUT) :: xi(3,nnode),xbndyi(3,nbpts)
109 LOGICAL,
INTENT(IN) :: modflag
113 INTEGER :: ioctr,
n,nedge,nfail
114 INTEGER :: idgp(mnode),idone(mnode),iedkp(4,medkp),iflag(mnode), &
115 ikeep(mtest),ipoint(mnode),iprot(mcell),iring(mring), &
116 ishk(mring),hk(mring), &
117 ityp(mnode),kshk(mring),ksrch(mtest),ldel(mtest), &
118 listf(mcell),lnbr(mring),lnkdn(mcell),lnkup(mcell), &
119 mnbr(mring),nacpt(mcell),nbh(4,mcell),nbhkp(3,mring), &
120 ncav(4,mcav),ncavfc(3,mtest),ndc(4,mcell),ndg(2,medge), &
121 ndgp(medge),nedgrm(mtest),newc(mring),newcel(mtest), &
122 newnbh(4,mtest),nfad(3,mring),nfce(3,mbface), &
123 nfill(mtest),nflag(mcell),nlink(mnode),noctr(2,moctr), &
124 nold(mtest),npoint(mcell),npp(mring),nprop(mbface), &
125 nptet(mnode),nref(mnode),nshake(mtest),nsrch(mtest), &
126 ntetkp(mring),ntri(3,mtest),nvcnt(mnode)
127 DOUBLE PRECISION :: rcmx,tolv,volmin
128 DOUBLE PRECISION :: count(mnode),dens(mnode),ds(mring),
dx(mring), &
129 dy(mring),
dz(mring),fac(mnode),rad(mtest), &
130 rat(mcell),rc(mcell),rcrin(mtest), &
131 resid(mnode),sig1(mcell),sig2(mcell), &
132 sig3(mcell),sv(medge),
v(mtest),vlt(medkp), &
133 vol(mcell),xc(mtest),xcen(mcell),xch(3,mnode), &
134 xfar(2),xhold(2,mtest),xkeep(2), &
135 xoctr(2,mtest),yc(mtest),ycen(mcell),yfar(2), &
136 yhold(2,mtest),ykeep(2),yoctr(2,mtest), &
137 zc(mtest),zcen(mcell),zfar(2),zhold(2,mtest), &
138 zoctr(2,mtest),zkeep(2),
x(3,mnode), &
139 xbndy(3,mbpts),xnewbn(3,mnode)
162 nfce(1,
n) = nfcei(1,
n)
163 nfce(2,
n) = nfcei(2,
n)
164 nfce(3,
n) = nfcei(3,
n)
181 xbndy(1,
n) = xbndyi(1,
n)
182 xbndy(2,
n) = xbndyi(2,
n)
183 xbndy(3,
n) = xbndyi(3,
n)
196 5 CALL
input(
x,nnode,ndc,ncell,nfce,nbpts,nbface, &
197 ityp,nprop,xbndy,xfar,yfar,zfar)
201 CALL
struct(
x,nnode,ndc,nbh,iprot,ncell,nfce,nbface, &
202 nedge,ndg,idgp,ndgp,ipoint,npoint,nptet, &
203 xcen,ycen,zcen,vol,rc,rat, &
205 ioctr,nlink,noctr,idone,nref,volmin,rcmx)
209 CALL
radrat(
x,nnode,ndc,nbh,iprot,ncell,nfce,nbface, &
210 ityp,ipoint,vol,rc,rat)
214 CALL
tetmv(
x,xnewbn,nnode,ndc,nbh,ncell,nfce,nbface,nfail, &
215 ityp,xcen,ycen,zcen,vol,rc,rat, &
216 nedge,ndg,idgp,ndgp,ipoint,count,xch,resid,sv, &
217 sig1,sig2,sig3,xfar,yfar,zfar,tolv)
225 IF ( modflag .EQV. .true. )
THEN
226 CALL
tetmod(
x,ityp,nnode,ndc,nbh,iprot,ncell, &
227 ndg,idgp,ndgp,nedge,nfce,nbface, &
228 vol,xcen,ycen,zcen,rc,rat,dens,nptet,nacpt, &
229 sig1,sig2,sig3,nvcnt,resid,count,fac, &
230 idone,nref,nlink,noctr,ioctr,xfar,yfar,zfar, &
231 xoctr,yoctr,zoctr,xhold,yhold,zhold, &
232 xkeep,ykeep,zkeep,ksrch,nsrch, &
233 ipoint,npoint,iflag,nflag, &
234 dx,
dy,
dz,ds,vlt,iring,ntetkp,nfad,newc, &
235 nbhkp,iedkp,lnbr,ishk,mnbr,kshk,npp, &
237 ncav,nshake,newnbh,nold,ncavfc,ikeep,ldel, &
238 nedgrm,xc,yc,zc,
v,rad,rcrin,lnkup,lnkdn, &
239 listf,volmin,rcmx,tolv)
243 CALL
radrat(
x,nnode,ndc,nbh,iprot,ncell,nfce,nbface, &
244 ityp,ipoint,vol,rc,rat)
254 610
FORMAT(//5
x,
'FAILURE IN MESH MOVEMENT PHASE')
255 620
FORMAT(/5
x,
'****************************************************'/ &
256 5
x,
'****************************************************'/ &
258 5
x,
'** CYCLE ',i3,
' OF MESH MODIFICATION IS COMPLETE **'/ &
259 5
x,
'** PARAMETRIC TIME = ',f8.4,
' **'/ &
261 5
x,
'****************************************************'/ &
262 5
x,
'****************************************************')
263 630
FORMAT(/5
x,
'MESH MOVEMENT SCHEME HAS FAILED IN STEP ',i4,
' OF A ', &
264 i4,
' STAGE CYCLE.'/ &
265 5
x,
'DOUBLE THE NUMBER OF STAGES AND REPEAT PROCEDURE.'/)
275 SUBROUTINE input (X,NNODE,NDC,NCELL,NFCE,NBPTS,NBFACE, &
276 ityp,nprop,xbndy,xfar,yfar,zfar)
292 INTEGER :: nbface,nbpts,ncell,nnode
294 INTEGER :: ityp(*),nprop(*)
296 DOUBLE PRECISION ::
x(3,*)
297 DOUBLE PRECISION :: xbndy(3,*)
298 DOUBLE PRECISION :: xfar(2),yfar(2),zfar(2)
300 INTEGER ::
j,jamin,jmax,jmin,l,
n,n1,n2,n3
301 DOUBLE PRECISION :: angl,angl1,angl2,angl3,angmax,angmin,farea, &
302 fmax,fmin,
q,qmax,qmin
315 IF ( nbpts > mxbpts )
THEN
319 IF ( nnode > mxnode )
THEN
323 IF ( nbface > mxface )
THEN
327 IF ( ncell > mxcell )
THEN
357 WRITE (6,700) nbface,nnode,ncell
373 CALL
fangle(
j,
x,nfce,angl1,angl2,angl3,
q)
375 angl =
min(angl1,angl2,angl3)
376 IF (angl.GT.angmin) go to 60
379 60 angmax =
max(angl1,angl2,angl3,angmax)
383 IF (
j.EQ.1) fmax = farea
384 IF (
j.EQ.1) fmin = farea
385 IF (farea.LT.fmax) go to 70
388 70
IF (farea.GT.fmin) go to 75
392 WRITE (6,705) angmin,angmax,qmin,qmax
396 IF (angmin.LT.1.)
WRITE (6,706) jamin,n1,
x(1,n1),
x(2,n1),
x(3,n1), &
397 n2,
x(1,n2),
x(2,n2),
x(3,n2), &
398 n3,
x(1,n3),
x(2,n3),
x(3,n3)
402 WRITE (6,707) fmax,fmin,jmin,n1,
x(1,n1),
x(2,n1),
x(3,n1), &
403 n2,
x(1,n2),
x(2,n2),
x(3,n2), &
404 n3,
x(1,n3),
x(2,n3),
x(3,n3)
415 xfar(1) =
min(xfar(1),
x(1,
n))
416 xfar(2) =
max(xfar(2),
x(1,
n))
417 yfar(1) =
min(yfar(1),
x(2,
n))
418 yfar(2) =
max(yfar(2),
x(2,
n))
419 zfar(1) =
min(zfar(1),
x(3,
n))
420 zfar(2) =
max(zfar(2),
x(3,
n))
422 WRITE (6,920) xfar(1),xfar(2),yfar(1),yfar(2),zfar(1),zfar(2)
423 920
FORMAT(
'XFAR ',f10.2,1
x,f10.2,
' YFAR ',f10.2,1
x,f10.2, &
424 ' ZFAR ',f10.2,1
x,f10.2)
440 300
WRITE (6,600) nbpts
442 305
WRITE (6,605) nnode
444 310
WRITE (6,610) nbface
446 315
WRITE (6,615) ncell
453 515
FORMAT(3e13.5,i10)
456 600
FORMAT(///5
x,
'NUMBER OF SURFACE POINTS ',i7/ &
457 5
x,
'EXCEEDS MAXIMUM ALLOWED. INCREASE SIZE OF MBPTS'/ &
458 5
x,
'PROGRAM STOPPED IN ROUTINE INPUT')
459 605
FORMAT(///5
x,
'TOTAL NUMBER OF MESH POINTS ',i7/ &
460 5
x,
'EXCEEDS MAXIMUM ALLOWED. INCREASE SIZE OF MNODE'/ &
461 5
x,
'PROGRAM STOPPED IN ROUTINE INPUT')
462 610
FORMAT(///5
x,
'NUMBER OF SURFACE FACES ',i7/ &
463 5
x,
'EXCEEDS MAXIMUM ALLOWED. INCREASE SIZE OF MBPTS'/ &
464 5
x,
'PROGRAM STOPPED IN ROUTINE INPUT')
465 615
FORMAT(///5
x,
'TOTAL NUMBER OF MESH CELLS ',i7/ &
466 5
x,
'EXCEEDS MAXIMUM ALLOWED. INCREASE SIZE OF MCELL'/ &
467 5
x,
'PROGRAM STOPPED IN ROUTINE INPUT')
468 700
FORMAT(/5
x,
'SURFACE AND VOLUME MESH READ'// &
469 5
x,
'NUMBER OF BOUNDARY SURFACE TRIANGLES = ',i7/ &
470 5
x,
'TOTAL NUMBER OF MESH POINTS = ',i7/ &
471 5
x,
'TOTAL NUMBER OF MESH CELLS = ',i7)
472 705
FORMAT(/5
x,
'MINIMUM BOUNDARY FACE ANGLE = ',f6.2/ &
473 5
x,
'MAXIMUM BOUNDARY FACE ANGLE = ',f6.2/ &
474 5
x,
'MINIMUM BOUNDARY RADIUS RATIO = ',f6.2/ &
475 5
x,
'MAXIMUM BOUNDARY RADIUS RATIO = ',f6.2)
476 706
FORMAT(/5
x,
'WARNING ! MINIMUM FACE ANGLE IS LESS THAN 1 DEGREE'/ &
477 5
x,
'FACE ADDRESS ',i6/ &
478 5
x,
'VERTEX ',i6,
' X = ',f12.4,
' Y = ',f12.4,
' Z = ',f12.4/ &
479 5
x,
'VERTEX ',i6,
' X = ',f12.4,
' Y = ',f12.4,
' Z = ',f12.4/ &
480 5
x,
'VERTEX ',i6,
' X = ',f12.4,
' Y = ',f12.4,
' Z = ',f12.4)
481 707
FORMAT(/5
x,
'MAX BOUNDARY FACE AREA = ',e13.5/ &
482 5
x,
'MIN BOUNDARY FACE AREA = ',e13.5/ &
483 5
x,
'ADDRESS OF BOUNDARY FACE WITH MINIMUM AREA',i6/ &
484 5
x,
'VERTEX ',i6,
' X = ',f12.4,
' Y = ',f12.4,
' Z = ',f12.4/ &
485 5
x,
'VERTEX ',i6,
' X = ',f12.4,
' Y = ',f12.4,
' Z = ',f12.4/ &
486 5
x,
'VERTEX ',i6,
' X = ',f12.4,
' Y = ',f12.4,
' Z = ',f12.4/)
496 SUBROUTINE struct (X,NNODE,NDC,NBH,IPROT,NCELL,NFCE,NBFACE, &
497 nedge,ndg,idgp,ndgp,ipoint,npoint,nptet, &
498 xcen,ycen,zcen,vol,rc,rat, &
500 ioctr,nlink,noctr,idone,nref,volmin,rcmx)
516 INTEGER :: ioctr,nbface,ncell,nedge,nnode
517 INTEGER :: ndc(4,*),ndg(2,*),idgp(*),ndgp(*)
518 INTEGER :: nbh(4,*),iprot(*),nptet(*),npoint(*),ipoint(*)
520 INTEGER :: idone(*),nref(*),nlink(*),noctr(2,*)
521 DOUBLE PRECISION :: volmin,rcmx
522 DOUBLE PRECISION ::
x(3,*)
523 DOUBLE PRECISION :: vol(*),xcen(*),ycen(*),zcen(*),rc(*),rat(*)
524 DOUBLE PRECISION :: xfar(2),yfar(2),zfar(2)
526 INTEGER ::
i,iedg,
j,
k,kk,l,lk,lm,ln,l1,l2,mm,m1,m2,m3,m4,
n,ncelm, &
527 nchg,ncnt,nend,ninc,nlook,nn,nold,n1,n2,n3,n4
528 DOUBLE PRECISION :: area,c5,c6,c9,fac,h1,h2,h3,h4,rnx,rny,rnz,tol, &
529 vcell,vmax,vmin,vrat,vsep,xshf,yshf,zshf
530 DOUBLE PRECISION ::
v(4),ar(4)
553 rnx =
cofact(
x(2,n1),
x(2,n2),
x(2,n3),
x(3,n1),
x(3,n2),
x(3,n3))
554 rny =
cofact(
x(3,n1),
x(3,n2),
x(3,n3),
x(1,n1),
x(1,n2),
x(1,n3))
555 rnz =
cofact(
x(1,n1),
x(1,n2),
x(1,n3),
x(2,n1),
x(2,n2),
x(2,n3))
559 v(ncnt) = rnx*(
x(1,n4) -
x(1,n1)) +rny*(
x(2,n4) -
x(2,n1)) &
560 +rnz*(
x(3,n4) -
x(3,n1))
561 ar(ncnt) =
sqrt(rnx*rnx +rny*rny +rnz*rnz)
562 IF (n1.EQ.nend) go to 20
569 20 vol(l) = .25*(abs(
v(1)) +abs(
v(2)) +abs(
v(3)) +abs(
v(4)))
570 area = ar(1) +ar(2) +ar(3) +ar(4)
576 IF (vol(l).LT.tol) go to 30
578 xshf = .5*(
x(1,n1) +
x(1,n2) +
x(1,n3) +
x(1,n4))
579 yshf = .5*(
x(2,n1) +
x(2,n2) +
x(2,n3) +
x(2,n4))
580 zshf = .5*(
x(3,n1) +
x(3,n2) +
x(3,n3) +
x(3,n4))
581 h1 = (
x(1,n4) -
x(1,n1))*(
x(1,n4) -xshf +
x(1,n1)) &
582 +(
x(2,n4) -
x(2,n1))*(
x(2,n4) -yshf +
x(2,n1)) &
583 +(
x(3,n4) -
x(3,n1))*(
x(3,n4) -zshf +
x(3,n1))
584 h2 = (
x(1,n3) -
x(1,n1))*(
x(1,n3) -xshf +
x(1,n1)) &
585 +(
x(2,n3) -
x(2,n1))*(
x(2,n3) -yshf +
x(2,n1)) &
586 +(
x(3,n3) -
x(3,n1))*(
x(3,n3) -zshf +
x(3,n1))
587 h3 = (
x(1,n2) -
x(1,n1))*(
x(1,n2) -xshf +
x(1,n1)) &
588 +(
x(2,n2) -
x(2,n1))*(
x(2,n2) -yshf +
x(2,n1)) &
589 +(
x(3,n2) -
x(3,n1))*(
x(3,n2) -zshf +
x(3,n1))
590 c5 = h2*(
x(3,n2) -
x(3,n1)) -h3*(
x(3,n3) -
x(3,n1))
591 c6 = h2*(
x(2,n2) -
x(2,n1)) -h3*(
x(2,n3) -
x(2,n1))
592 c9 = h3*(
x(1,n3) -
x(1,n1)) -h2*(
x(1,n2) -
x(1,n1))
593 xcen(l) = .5*xshf +(h1*rnx +(
x(2,n4) -
x(2,n1))*c5 &
594 -(
x(3,n4) -
x(3,n1))*c6)*fac
595 ycen(l) = .5*yshf +(-(
x(1,n4) -
x(1,n1))*c5 +h1*rny &
596 -(
x(3,n4) -
x(3,n1))*c9)*fac
597 zcen(l) = .5*zshf +((
x(1,n4) -
x(1,n1))*c6 &
598 +(
x(2,n4) -
x(2,n1))*c9 +h1*rnz)*fac
599 rc(l) =
sqrt((
x(1,n1) -xcen(l))**2 +(
x(2,n1) -ycen(l))**2 &
600 +(
x(3,n1) -zcen(l))**2)
601 rat(l) = rc(l)*area/vol(l)
602 vmin =
min(abs(
v(1)),abs(
v(2)),abs(
v(3)),abs(
v(4)))
603 vmax =
max(abs(
v(1)),abs(
v(2)),abs(
v(3)),abs(
v(4)))
607 IF (l.EQ.1) volmin = vol(1)
608 IF (l.EQ.1) rcmx = rc(1)
609 volmin =
min(volmin,vol(l))
610 rcmx =
max(rcmx,rc(l))
621 ndc(1,ncell) = nfce(1,
j)
622 ndc(2,ncell) = nfce(2,
j)
623 ndc(3,ncell) = nfce(3,
j)
635 IF (nptet(
n).GT.0) go to 55
644 IF (npoint(l).EQ.4) go to 120
649 IF (n1.LT.0) go to 100
651 IF (
j.EQ.0.OR.
j.EQ.l) go to 100
656 70
IF (m1.EQ.n1) go to 75
664 IF (n2.EQ.m2.OR.n2.EQ.m3.OR.n2.EQ.m4) ninc = ninc +1
665 IF (n3.EQ.m2.OR.n3.EQ.m3.OR.n3.EQ.m4) ninc = ninc +1
666 IF (n4.EQ.m2.OR.n4.EQ.m3.OR.n4.EQ.m4) ninc = ninc +1
667 IF (ninc.LT.2) go to 100
668 80
IF (npoint(l).EQ.0) go to 90
671 IF (ln.EQ.
j) go to 100
673 90 npoint(l) = npoint(l) +1
675 npoint(
j) = npoint(
j) +1
677 100
IF (n1.EQ.ndc(4,l)) go to 120
683 IF (n1.LT.0) go to 100
692 IF (npoint(l).EQ.4) go to 125
696 IF (nchg.EQ.0) go to 140
702 IF (
n.GT.0) nptet(
n) = 0
708 IF (
n.LT.0) go to 135
709 IF (nptet(
n).EQ.0) nptet(
n) = l
718 IF (npoint(l).NE.4) go to 320
721 IF (l.GT.ncelm) go to 155
725 IF (
n.GT.0) nptet(
n) = l
740 CALL
octfil(
n,
x,noctr,ioctr,nlink,nref,xfar,yfar,zfar)
753 IF (ndc(4,
j).EQ.-1) go to 200
756 IF (
i.EQ.5) go to 200
760 IF (
k.EQ.
i) go to 180
765 190
IF (nold.EQ.0) go to 195
766 IF (l2.EQ.
max(ndg(1,nold),ndg(2,nold))) go to 185
774 IF (idgp(l1).NE.0) ndgp(iedg) = nedge
775 IF (idgp(l1).EQ.0) idgp(l1) = nedge
779 320
WRITE (6,620) l,(ndc(kk,l),kk=1,4),npoint(l)
781 600
FORMAT(//5
x,
'TETRAHEDRON WITH AN EXTREMELY SMALL VOLUME FOUND'// &
782 ' IN ROUTINE STRUCT'//5
x,
'VOLUME = ',e13.5/)
783 610
FORMAT(5
x,
'IMPRECISE ESTIMATE OF TETRAHEDRON VOLUME', &
784 5
x,
'VSEP = ',e13.5,
' VOLUME = ',e13.5, &
785 ' VSEP/VOLUME = ',e13.5)
786 620
FORMAT(//5
x,
'NPOINT IS NOT EQUAL TO FOUR FOR TETRAHEDRON ',i7, &
787 5
x,
'TETRAHEDRON VERTICES ARE ',4i6,
' AND NPOINT = ',i2, &
788 5
x,
'PROGRAM STOPPED IN ROUTINE STRUCT')
796 SUBROUTINE tetmv (X,XNEWBN,NNODE,NDC,NBH,NCELL,NFCE,NBFACE,NFAIL, &
797 ityp,xcen,ycen,zcen,vol,rc,rat, &
798 nedge,ndg,idgp,ndgp,ipoint,fcount,xc,resid,sv, &
799 sig1,sig2,sig3,xfar,yfar,zfar,tolv)
827 INTEGER :: nbface,ncell,nedge,nfail,nnode
828 INTEGER :: ndc(4,*),nbh(4,*),ityp(*),nfce(3,*)
829 INTEGER :: ndg(2,*),idgp(*),ndgp(*),ipoint(*)
830 DOUBLE PRECISION :: tolv
831 DOUBLE PRECISION ::
x(3,*),xnewbn(3,*)
832 DOUBLE PRECISION :: xcen(*),ycen(*),zcen(*),vol(*),rc(*),rat(*)
833 DOUBLE PRECISION :: fcount(*),xc(3,*),resid(*),sv(*), &
834 sig1(*),sig2(*),sig3(*)
835 DOUBLE PRECISION :: rmax(3),itmax(3)
836 DOUBLE PRECISION :: xfar(2),yfar(2),zfar(2)
837 DOUBLE PRECISION :: xcn,ycn,zcn,vl,radc
838 DOUBLE PRECISION :: vmin,vmax
840 INTEGER ::
i,ibig,iedg,ismall,it,
j,
k,kp1,kp2,kp3,l,loop,lsub, &
841 l1,l2,
n,nneg,nold,npass,n1,n2,n3,n4
842 DOUBLE PRECISION :: area,cndmax,cndmin,difx,eps,fac,rbig,resm, &
843 rmax1,rnx,rny,rnz,sigmx,sigmn,strmax,strmin, &
844 strtet,third,vtet,vtet1,vtet2,xbar,ybar,zbar
869 IF (nbh(1,
j).EQ.0) go to 22
870 IF (ndc(4,
j).EQ.-1) go to 22
880 xbar = third*(
x(1,n1) +
x(1,n2) +
x(1,n3))
881 ybar = third*(
x(2,n1) +
x(2,n2) +
x(2,n3))
882 zbar = third*(
x(3,n1) +
x(3,n2) +
x(3,n3))
883 rnx =
cofact(
x(2,n1),
x(2,n2),
x(2,n3),
x(3,n1),
x(3,n2),
x(3,n3))
884 rny =
cofact(
x(3,n1),
x(3,n2),
x(3,n3),
x(1,n1),
x(1,n2),
x(1,n3))
885 rnz =
cofact(
x(1,n1),
x(1,n2),
x(1,n3),
x(2,n1),
x(2,n2),
x(2,n3))
886 vtet = (
x(1,n4) -xbar)*rnx +(
x(2,n4) -ybar)*rny &
888 vol(
j) = vol(
j) +abs(vtet)
893 IF (nbh(1,
j).EQ.0) go to 25
894 IF (ndc(4,
j).EQ.-1) go to 25
900 vmin =
min(vmin,vol(
j))
901 vmax =
max(vmax,vol(
j))
905 IF (nbh(1,
j).EQ.0) go to 30
906 IF (ndc(4,
j).EQ.-1) go to 30
907 vol(
j) = 1. +(vmax -vmin)/vol(
j)
909 WRITE (6,940) vmin,vmax
910 940
FORMAT(/
'VMIN = ',e13.5,
' VMAX = ',e13.5/)
916 IF (nbh(1,
j).EQ.0) go to 70
917 IF (ndc(4,
j).EQ.-1) go to 70
929 60
IF (nold.EQ.0) go to 65
930 IF (l2.EQ.
max(ndg(1,nold),ndg(2,nold)))
THEN
931 sv(nold) = sv(nold) +vol(
j)
941 IF (idgp(l1).NE.0) ndgp(iedg) = nedge
942 IF (idgp(l1).EQ.0) idgp(l1) = nedge
957 IF (ityp(
n).LT.0) go to 85
967 IF (ityp(
n).LT.0) go to 90
968 IF (ipoint(
n).GT.0) go to 90
970 xc(loop,
n) = (xnewbn(loop,
n) -
x(loop,
n))*l*fac
972 IF (loop.GT.1) go to 105
979 IF (n1.LE.0.OR.n2.LE.0) go to 95
980 IF (ityp(n1).LT.0.OR.ityp(n2).LT.0) go to 95
981 fcount(n1) = fcount(n1) +sv(
i)
982 fcount(n2) = fcount(n2) +sv(
i)
985 IF (ityp(
n).LT.0) go to 100
986 fcount(
n) = 1./fcount(
n)
1005 IF (n1.LE.0.OR.n2.LE.0) go to 115
1006 IF (ityp(n1).LT.0.OR.ityp(n2).LT.0) go to 115
1007 difx = (xc(loop,n1) -xc(loop,n2))*sv(
i)
1008 resid(n1) = resid(n1) -difx
1009 resid(n2) = resid(n2) +difx
1016 IF (ityp(
n).LT.0) go to 120
1017 IF (ipoint(
n).EQ.0)
THEN
1018 xc(loop,
n) = xc(loop,
n) +eps*fcount(
n)*resid(
n)
1019 resm = abs(resid(
n))
1020 rmax(loop) =
max(rmax(loop),fcount(
n)*resm)
1027 IF (it.EQ.1) rmax1 = rmax(loop)
1028 IF (rmax(loop).GT.
max(0.0001*rmax1,1.0d0-9).AND.it.LT.2000) &
1031 IF (loop.LT.3) go to 80
1032 rbig =
max(rmax(1),rmax(2),rmax(3))
1033 ibig =
max(itmax(1),itmax(2),itmax(3))
1034 WRITE (6,600) rbig,ibig
1043 IF (nbh(1,
j).EQ.0) go to 130
1044 IF (ndc(4,
j).EQ.-1) go to 130
1049 xbar = third*(
x(1,n1) +
x(1,n2) +
x(1,n3))
1050 ybar = third*(
x(2,n1) +
x(2,n2) +
x(2,n3))
1051 zbar = third*(
x(3,n1) +
x(3,n2) +
x(3,n3))
1052 rnx =
cofact(
x(2,n1),
x(2,n2),
x(2,n3),
x(3,n1),
x(3,n2),
x(3,n3))
1053 rny =
cofact(
x(3,n1),
x(3,n2),
x(3,n3),
x(1,n1),
x(1,n2),
x(1,n3))
1054 rnz =
cofact(
x(1,n1),
x(1,n2),
x(1,n3),
x(2,n1),
x(2,n2),
x(2,n3))
1055 vtet1 = rnx*(
x(1,n4) -xbar) +rny*(
x(2,n4) -ybar) &
1056 +rnz*(
x(3,n4) -zbar)
1057 xbar = xbar +third*(xc(1,n1) +xc(1,n2) +xc(1,n3))
1058 ybar = ybar +third*(xc(2,n1) +xc(2,n2) +xc(2,n3))
1059 zbar = zbar +third*(xc(3,n1) +xc(3,n2) +xc(3,n3))
1060 rnx =
cofact(
x(2,n1)+xc(2,n1),
x(2,n2)+xc(2,n2),
x(2,n3)+xc(2,n3), &
1061 x(3,n1)+xc(3,n1),
x(3,n2)+xc(3,n2),
x(3,n3)+xc(3,n3))
1062 rny =
cofact(
x(3,n1)+xc(3,n1),
x(3,n2)+xc(3,n2),
x(3,n3)+xc(3,n3), &
1063 x(1,n1)+xc(1,n1),
x(1,n2)+xc(1,n2),
x(1,n3)+xc(1,n3))
1064 rnz =
cofact(
x(1,n1)+xc(1,n1),
x(1,n2)+xc(1,n2),
x(1,n3)+xc(1,n3), &
1065 x(2,n1)+xc(2,n1),
x(2,n2)+xc(2,n2),
x(2,n3)+xc(2,n3))
1066 vtet2 = rnx*(
x(1,n4) +xc(1,n4) -xbar) &
1067 +rny*(
x(2,n4) +xc(2,n4) -ybar) &
1068 +rnz*(
x(3,n4) +xc(3,n4) -zbar)
1069 strtet = abs(vtet2/vtet1)
1070 strmax =
max(strmax,strtet)
1071 strmin =
min(strmin,strtet)
1072 IF (vtet1*vtet2.LE.0.) nneg = nneg +1
1075 CALL
snglar(nnode,
x,xc,ncell,ndc,nbh,sigmn,sigmx, &
1076 sig1,sig2,sig3,cndmin,cndmax)
1078 WRITE (6,700) strmax,strmin,sigmx,sigmn,cndmax,cndmin
1079 IF (nneg.GT.0)
WRITE (6,705) nneg
1089 IF (ityp(
n).LT.0) go to 140
1091 x(1,
n) =
x(1,
n) +xc(1,
n)
1092 x(2,
n) =
x(2,
n) +xc(2,
n)
1093 x(3,
n) =
x(3,
n) +xc(3,
n)
1096 920
FORMAT(/
'SUB-ITERATION ',i3,
' IS COMPLETE'/)
1103 IF (ityp(
n).LT.0) go to 155
1104 IF (npass.EQ.0)
THEN
1113 xfar(1) =
min(xfar(1),
x(1,
n))
1114 xfar(2) =
max(xfar(2),
x(1,
n))
1115 yfar(1) =
min(yfar(1),
x(2,
n))
1116 yfar(2) =
max(yfar(2),
x(2,
n))
1117 zfar(1) =
min(zfar(1),
x(3,
n))
1118 zfar(2) =
max(zfar(2),
x(3,
n))
1121 WRITE (6,925) xfar(1),xfar(2),yfar(1),yfar(2),zfar(1),zfar(2)
1122 925
FORMAT(
'XFAR ',f10.2,1
x,f10.2,
' YFAR ',f10.2,1
x,f10.2, &
1123 ' ZFAR ',f10.2,1
x,f10.2)
1129 IF (nbh(1,
j).EQ.0) go to 160
1134 IF (n4.EQ.-1) go to 160
1136 CALL
circum(
x,n1,n2,n3,n4,xcn,ycn,zcn,vl,radc,ismall,tolv)
1138 IF (ismall.EQ.1) go to 310
1145 rat(
j) = rc(
j)*area/vol(
j)
1152 600
FORMAT(/5
x,
'MAXIMUM RESIDUAL IS ',f10.6, &
1153 ' MAXIMUM ITERATIONS = ',i5)
1154 610
FORMAT(///5
x,
'AT LEAST ONE NEW TETRAHEDRON HAS TOO SMALL A VOLUME' &
1155 /5
x,
'PROGRAM STOPPED IN TETMV')
1156 700
FORMAT(/5
x,
'****************************************************'/ &
1158 5
x,
'** MAXIMUM CELL STRETCHING IS ',f10.4,
' **'/ &
1159 5
x,
'** MAXIMUM CELL COMPRESSION IS ',f10.4,
' **'/ &
1161 5
x,
'** MAXIMUM SINGULAR VALUE IS ',f10.4,
' **'/ &
1162 5
x,
'** MINIMUM SINGULAR VALUE IS ',f10.4,
' **'/ &
1164 5
x,
'** MAXIMUM CONDITION NUMBER ',f10.4,
' **'/ &
1165 5
x,
'** MINIMUM CONDITION NUMBER ',f10.4,
' **'/ &
1167 5
x,
'****************************************************'/)
1168 705
FORMAT(/5
x,
'THERE ARE ',i6,
' CELLS THAT HAVE BEEN INVERTED.'/ &
1169 5
x,
'THE NEW MESH IS THEREFORE INVALID.')
1171 END SUBROUTINE tetmv
1180 SUBROUTINE tetmod (X,ITYP,NNODE,NDC,NBH,IPROT,NCELL, &
1181 ndg,idgp,ndgp,nedge,nfce,nbface, &
1182 vol,xcen,ycen,zcen,rc,rat,dens,nptet,nacpt, &
1183 sig1,sig2,sig3,nvcnt,resid,count,fac, &
1184 idone,nref,nlink,noctr,ioctr,xfar,yfar,zfar, &
1185 xoctr,yoctr,zoctr,xhold,yhold,zhold, &
1186 xkeep,ykeep,zkeep,ksrch,nsrch, &
1187 ipoint,npoint,iflag,nflag, &
1188 dx,
dy,
dz,ds,vlt,iring,ntetkp,nfad,newc, &
1189 nbhkp,iedkp,lnbr,ishk,mnbr,kshk,npp, &
1190 nfill,newcel,ntri, &
1191 ncav,nshake,newnbh,nold,ncavfc,ikeep,ldel, &
1192 nedgrm,xc,yc,zc,
v,rad,rcrin,lnkup,lnkdn, &
1193 listf,volmin,rcmx,tolv)
1209 INTEGER :: ioctr,nbface,ncell,nedge,nnode
1210 INTEGER :: lnkup(*),lnkdn(*)
1211 INTEGER :: nsrch(*),ksrch(*)
1213 INTEGER :: idone(*),nref(*),nlink(*),noctr(2,*)
1214 INTEGER :: ityp(*),nptet(*),listf(*),nacpt(*)
1215 INTEGER :: npoint(*),ipoint(*)
1216 INTEGER :: iflag(*),nflag(*)
1217 INTEGER :: ndc(4,*),nbh(4,*),iprot(*),ndg(2,*),idgp(*),ndgp(*)
1218 INTEGER :: nfce(3,*)
1219 INTEGER :: ntri(3,*),nfill(*),newnbh(4,*),nold(*),newcel(*), &
1220 nshake(*),ncav(4,*),nedgrm(*)
1221 INTEGER :: ldel(*),ncavfc(3,*),ikeep(*)
1222 INTEGER :: iring(*),ntetkp(*),nfad(3,*),newc(*), &
1223 nbhkp(3,*),iedkp(4,*),npp(*), &
1224 lnbr(*),ishk(*),mnbr(*),kshk(*)
1225 DOUBLE PRECISION :: volmin,tolv,rcmx
1226 DOUBLE PRECISION ::
x(3,*),dens(*)
1227 DOUBLE PRECISION :: vol(*),xcen(*),ycen(*),zcen(*),rc(*),rat(*)
1228 DOUBLE PRECISION :: sig1(*),sig2(*),sig3(*)
1229 DOUBLE PRECISION :: resid(*),count(*),fac(*)
1230 DOUBLE PRECISION :: xoctr(2,*),yoctr(2,*),zoctr(2,*), &
1231 xhold(2,*),yhold(2,*),zhold(2,*), &
1232 xfar(2),yfar(2),zfar(2),xkeep(2), &
1234 DOUBLE PRECISION :: xc(*),yc(*),zc(*),
v(*),rad(*),rcrin(*)
1235 DOUBLE PRECISION ::
dx(*),
dy(*),
dz(*),ds(*),vlt(*)
1237 INTEGER :: jfirst,jlast,lc,nbpts,ntrack
1243 CALL
densfn(
x,nnode,nfce,nbface,nedge,ndg,dens,resid,count)
1247 CALL
edglen(
x,nnode,ityp,nedge,ndg,dens,nvcnt)
1253 CALL
coarsn(lc,
x,nnode,ndc,nbh,iprot,ncell, &
1254 ityp,xcen,ycen,zcen,vol,rc,rat, &
1255 nvcnt,dens,iflag,nflag,nptet, &
1256 nedge,ndg,idgp,ndgp, &
1257 noctr,ioctr,nlink,xfar,yfar,zfar,idone,nref, &
1258 fac,sig1,sig2,sig3, &
1259 ksrch,nsrch,iring,ntetkp, &
1260 lnbr,ishk,mnbr,kshk,tolv)
1264 CALL
smooth(
x,ityp,nnode,ndc,nbh,iprot,ncell, &
1265 ndg,idgp,ndgp,nedge, &
1266 vol,xcen,ycen,zcen,rc,rat,dens,nptet,nacpt, &
1267 idone,nref,nlink,noctr,ioctr,xfar,yfar,zfar, &
1268 xoctr,yoctr,zoctr,xhold,yhold,zhold, &
1269 xkeep,ykeep,zkeep,ksrch,nsrch, &
1270 ipoint,npoint,iflag,nflag, &
1271 dx,
dy,
dz,ds,vlt,iring,ntetkp,nfad,newc, &
1272 nbhkp,iedkp,lnbr,ishk,mnbr,kshk,npp, &
1273 nfill,newcel,ntri, &
1274 ncav,nshake,newnbh,nold,ncavfc,ikeep,ldel, &
1275 nedgrm,xc,yc,zc,
v,rad,rcrin,lnkup,lnkdn, &
1276 listf,volmin,rcmx,tolv)
1280 CALL
edglen(
x,nnode,ityp,nedge,ndg,dens,nvcnt)
1286 CALL
densfn(
x,nnode,nfce,nbface,nedge,ndg,dens,resid,count)
1292 CALL
volput(
x,ityp,nbpts,nnode,ndc,nbh,iprot,ncell, &
1293 ndg,idgp,ndgp,nedge, &
1294 vol,xcen,ycen,zcen,rc,rat,dens,nptet,nacpt, &
1295 idone,nref,nlink,noctr,ioctr,xfar,yfar,zfar, &
1296 xoctr,yoctr,zoctr,xhold,yhold,zhold, &
1297 xkeep,ykeep,zkeep,ksrch,nsrch, &
1298 ipoint,npoint,iflag,nflag,nfill,newcel,ntri, &
1299 ncav,nshake,newnbh,nold,ncavfc,ikeep,ldel, &
1300 nedgrm,xc,yc,zc,
v,rad,rcrin,lnkup,lnkdn, &
1301 jlast,jfirst,ntrack,volmin,rcmx,tolv)
1305 CALL
smooth(
x,ityp,nnode,ndc,nbh,iprot,ncell, &
1306 ndg,idgp,ndgp,nedge, &
1307 vol,xcen,ycen,zcen,rc,rat,dens,nptet,nacpt, &
1308 idone,nref,nlink,noctr,ioctr,xfar,yfar,zfar, &
1309 xoctr,yoctr,zoctr,xhold,yhold,zhold, &
1310 xkeep,ykeep,zkeep,ksrch,nsrch, &
1311 ipoint,npoint,iflag,nflag, &
1312 dx,
dy,
dz,ds,vlt,iring,ntetkp,nfad,newc, &
1313 nbhkp,iedkp,lnbr,ishk,mnbr,kshk,npp, &
1314 nfill,newcel,ntri, &
1315 ncav,nshake,newnbh,nold,ncavfc,ikeep,ldel, &
1316 nedgrm,xc,yc,zc,
v,rad,rcrin,lnkup,lnkdn, &
1317 listf,volmin,rcmx,tolv)
1321 CALL
edglen(
x,nnode,ityp,nedge,ndg,dens,nvcnt)
1333 SUBROUTINE octfil (N,X,NOCTR,IOCTR,NLINK,NREF,XFAR,YFAR,ZFAR)
1350 INTEGER :: nref(*),nlink(*),noctr(2,*)
1351 DOUBLE PRECISION ::
x(3,*)
1352 DOUBLE PRECISION :: xfar(2),yfar(2),zfar(2)
1353 DOUBLE PRECISION :: xroot(2),yroot(2),zroot(2)
1355 INTEGER ::
i,iroot,
j,
k,l,ltest,
next,ntest,nxsgn,nysgn,nzsgn
1356 INTEGER :: jj(8),nstore(8)
1357 DOUBLE PRECISION :: tol,tolpt, &
1358 xhalf,xhigh,xlow,xsgn,xshift,xsize, &
1359 yhalf,yhigh,ylow,ysgn,yshift,ysize, &
1360 zhalf,zhigh,zlow,zsgn,zshift,zsize
1366 IF (
x(1,
n).LT.xfar(1)-tolpt.OR.
x(1,
n).GT.xfar(2)+tolpt) go to 200
1367 IF (
x(2,
n).LT.yfar(1)-tolpt.OR.
x(2,
n).GT.yfar(2)+tolpt) go to 200
1368 IF (
x(3,
n).LT.zfar(1)-tolpt.OR.
x(3,
n).GT.zfar(2)+tolpt) go to 200
1376 IF (noctr(1,
i).LT.0) go to 80
1381 IF (
next.NE.0) go to 6
1389 IF (l.EQ.8) go to 30
1393 IF (l.GT.1) go to 15
1396 15
next = noctr(1,
i)
1397 20 ntest = nlink(
next)
1398 IF (ntest.EQ.0) go to 25
1409 IF (ioctr.GT.mxoctr) go to 210
1411 35 noctr(1,ioctr) = 0
1422 xshift =
x(1,
j) -.5*(xlow +xhigh)
1423 xsize =
max(1.0d0-9,abs(xshift))
1424 nxsgn = (int(
REAL(tol*xshift/xsize)) +1)/2*2
1426 yshift =
x(2,
j) -.5*(ylow +yhigh)
1427 ysize =
max(1.0d0-9,abs(yshift))
1428 nysgn = (int(
REAL(tol*yshift/ysize)) +1)/2*2
1430 zshift =
x(3,
j) -.5*(zlow +zhigh)
1431 zsize =
max(1.0d0-9,abs(zshift))
1432 nzsgn = (int(
REAL(tol*zshift/zsize)) +1)/2*2
1434 l = 1 +nxsgn/2 +nysgn +2*nzsgn
1437 nref(
j) = ioctr -8 +l
1438 IF (jj(l).GT.1) go to 41
1439 noctr(1,ioctr-8+l) =
j
1441 41
next = noctr(1,ioctr-8+l)
1442 42 ntest = nlink(
next)
1443 IF (ntest.EQ.0) go to 43
1451 noctr(1,
i) = 7 -ioctr
1454 IF (jj(l).EQ.8) ltest = l
1456 IF (ltest.EQ.0) go to 70
1458 xhalf = .5*(xlow +xhigh)
1460 xlow = xsgn*xlow +(1. -xsgn)*xhalf
1461 xhigh = xsgn*xhalf +(1. -xsgn)*xhigh
1462 yhalf = .5*(ylow +yhigh)
1463 ysgn = isign(1,2*mod(ltest-1,4)-3)
1464 ylow = .5*((1. -ysgn)*ylow +(1. +ysgn)*yhalf)
1465 yhigh = .5*((1. -ysgn)*yhalf +(1. +ysgn)*yhigh)
1466 zhalf = .5*(zlow +zhigh)
1467 zsgn = isign(1,2*ltest-9)
1468 zlow = .5*((1. -zsgn)*zlow +(1. +zsgn)*zhalf)
1469 zhigh = .5*((1. -zsgn)*zhalf +(1. +zsgn)*zhigh)
1481 80 xhalf = .5*(xlow +xhigh)
1482 yhalf = .5*(ylow +yhigh)
1483 zhalf = .5*(zlow +zhigh)
1484 xshift =
x(1,
n) -xhalf
1485 xsize =
max(1.0d0-9,abs(xshift))
1486 nxsgn = (int(
REAL(tol*xshift/xsize)) +1)/2*2
1488 yshift =
x(2,
n) -yhalf
1489 ysize =
max(1.0d0-9,abs(yshift))
1490 nysgn = (int(
REAL(tol*yshift/ysize)) +1)/2*2
1492 zshift =
x(3,
n) -zhalf
1493 zsize =
max(1.0d0-9,abs(zshift))
1494 nzsgn = (int(
REAL(tol*zshift/zsize)) +1)/2*2
1496 l = 1 +nxsgn/2 +nysgn +2*nzsgn
1497 i = -noctr(1,
i) +l -1
1499 xlow = xsgn*xlow +(1. -xsgn)*xhalf
1500 xhigh = xsgn*xhalf +(1. -xsgn)*xhigh
1501 ysgn = isign(1,2*mod(l-1,4)-3)
1502 ylow = .5*((1. -ysgn)*ylow +(1. +ysgn)*yhalf)
1503 yhigh = .5*((1. -ysgn)*yhalf +(1. +ysgn)*yhigh)
1504 zsgn = isign(1,2*l-9)
1505 zlow = .5*((1. -zsgn)*zlow +(1. +zsgn)*zhalf)
1506 zhigh = .5*((1. -zsgn)*zhalf +(1. +zsgn)*zhigh)
1507 IF (noctr(1,
i).LT.0) go to 80
1509 IF (noctr(1,
i).EQ.0) go to 10
1513 IF (
next.NE.0) go to 85
1517 WRITE (6,900)
n,
x(1,
n),
x(2,
n),
x(3,
n), &
1518 xfar(1),xfar(2),yfar(1),yfar(2),zfar(1),zfar(2)
1519 900
FORMAT(
'N ',i4,
' X ',f8.3,
' Y ',f8.3,
' Z ',f8.3/ &
1520 'XFAR',2(1
x,f8.3),
' YFAR',2(1
x,f8.3),
' ZFAR',2(1
x,f8.3))
1525 600
FORMAT(5
x,
'POINT N= ',i6,
' LIES OUTSIDE CONVEX HULL')
1526 610
FORMAT(5
x,
'DIMENSION OF NOCTR ARRAY EXCEEDED'/ &
1527 5
x,
'INCREASE SIZE OF MOCTR')
1537 SUBROUTINE octfnd (NCODE,NCLOSE,XPT,YPT,ZPT,N1,N2,N3,DISMIN, &
1538 x,noctr,nlink,xfar,yfar,zfar,idone, &
1539 xoctr,yoctr,zoctr,xhold,yhold,zhold, &
1540 xkeep,ykeep,zkeep,ksrch,nsrch)
1557 INTEGER :: ncode,nclose,n1,n2,n3
1558 INTEGER :: idone(*),nlink(*),noctr(2,*)
1559 INTEGER :: nsrch(*),ksrch(*)
1560 DOUBLE PRECISION :: dismin,xpt,ypt,zpt
1561 DOUBLE PRECISION ::
x(3,*)
1562 DOUBLE PRECISION :: xoctr(2,*),yoctr(2,*),zoctr(2,*), &
1563 xhold(2,*),yhold(2,*),zhold(2,*), &
1564 xkeep(2),ykeep(2),zkeep(2),xfar(2), &
1567 INTEGER ::
i,
ic,iflag,ii,itry,
j,jj,
k,
kc,l,lflag,
next,nxsgn, &
1569 DOUBLE PRECISION :: c1,c2,c3,det,
dist,dmin,tol,tolpt, &
1570 xhalf,xhigh,xl,xlow,xshift,xsize,xsgn,xu, &
1571 yhalf,yhigh,yl,ylow,yshift,ysize,ysgn,yu, &
1572 zhalf,zhigh,zl,zlow,zshift,zsize,zsgn,zu
1578 IF (xpt.LT.xfar(1)-tolpt.OR.xpt.GT.xfar(2)+tolpt) go to 200
1579 IF (ypt.LT.yfar(1)-tolpt.OR.ypt.GT.yfar(2)+tolpt) go to 200
1580 IF (zpt.LT.zfar(1)-tolpt.OR.zpt.GT.zfar(2)+tolpt) go to 200
1581 IF (n1.EQ.0) go to 5
1585 c1 =
cofact(
x(2,n1),
x(2,n2),
x(2,n3),
x(3,n1),
x(3,n2),
x(3,n3))
1586 c2 =
cofact(
x(1,n1),
x(1,n2),
x(1,n3),
x(3,n1),
x(3,n2),
x(3,n3))
1587 c3 =
cofact(
x(1,n1),
x(1,n2),
x(1,n3),
x(2,n1),
x(2,n2),
x(2,n3))
1598 IF (noctr(1,
i).LT.0) go to 20
1601 IF (
next.NE.0) go to 10
1606 20 xhalf = .5*(xlow +xhigh)
1607 yhalf = .5*(ylow +yhigh)
1608 zhalf = .5*(zlow +zhigh)
1610 xsize =
max(1.0d0-9,abs(xshift))
1611 nxsgn = (int(
REAL(tol*xshift/xsize)) +1)/2*2
1614 ysize =
max(1.0d0-9,abs(yshift))
1615 nysgn = (int(
REAL(tol*yshift/ysize)) +1)/2*2
1618 zsize =
max(1.0d0-9,abs(zshift))
1619 nzsgn = (int(
REAL(tol*zshift/zsize)) +1)/2*2
1621 l = 1 +nxsgn/2 +nysgn +2*nzsgn
1622 i = -noctr(1,
i) +l -1
1624 xlow = xsgn*xlow +(1. -xsgn)*xhalf
1625 xhigh = xsgn*xhalf +(1. -xsgn)*xhigh
1626 ysgn = isign(1,2*mod(l-1,4)-3)
1627 ylow = .5*((1. -ysgn)*ylow +(1. +ysgn)*yhalf)
1628 yhigh = .5*((1. -ysgn)*yhalf +(1. +ysgn)*yhigh)
1629 zsgn = isign(1,2*l-9)
1630 zlow = .5*((1. -zsgn)*zlow +(1. +zsgn)*zhalf)
1631 zhigh = .5*((1. -zsgn)*zhalf +(1. +zsgn)*zhigh)
1632 IF (noctr(1,
i).LT.0) go to 20
1633 IF (noctr(1,
i).EQ.0) go to 50
1636 IF (
next.NE.0) go to 25
1641 dismin = (xfar(2) -xfar(1))**2 +(yfar(2) -yfar(1))**2 &
1642 +(zfar(2) -zfar(1))**2
1643 IF (noctr(1,
i).EQ.0) go to 65
1645 58
IF (idone(
j).EQ.0.AND.ncode.EQ.0) go to 60
1646 IF (idone(
j).EQ.2) go to 60
1647 IF (n1.EQ.0) go to 59
1648 det = c1*(
x(1,
j) -xpt) -c2*(
x(2,
j) -ypt) &
1650 IF (det.LT.tolpt) go to 60
1651 59
dist = (xpt -
x(1,
j))**2 +(ypt -
x(2,
j))**2 &
1653 IF (
dist.GE.dismin) go to 60
1657 IF (
j.NE.0) go to 58
1658 65 dmin =
sqrt(dismin)
1665 70
IF (xl.GT.xlow.AND.xu.LT.xhigh.AND. &
1666 yl.GT.ylow.AND.yu.LT.yhigh.AND. &
1667 zl.GT.zlow.AND.zu.LT.zhigh)
RETURN
1671 lflag = iflag +noctr(1,
i) +1
1673 xoctr(1,1) = (2. -xsgn)*xlow -(1. -xsgn)*xhigh
1674 xoctr(2,1) = (1. +xsgn)*xhigh -xsgn*xlow
1675 ysgn = isign(1,2*mod(lflag-1,4)-3)
1676 yoctr(1,1) = .5*((3. +ysgn)*ylow -(1. +ysgn)*yhigh)
1677 yoctr(2,1) = .5*((3. -ysgn)*yhigh -(1. -ysgn)*ylow)
1678 zsgn = isign(1,2*lflag-9)
1679 zoctr(1,1) = .5*((3. +zsgn)*zlow -(1. +zsgn)*zhigh)
1680 zoctr(2,1) = .5*((3. -zsgn)*zhigh -(1. -zsgn)*zlow)
1681 xkeep(1) = xoctr(1,1)
1682 xkeep(2) = xoctr(2,1)
1683 ykeep(1) = yoctr(1,1)
1684 ykeep(2) = yoctr(2,1)
1685 zkeep(1) = zoctr(1,1)
1686 zkeep(2) = zoctr(2,1)
1696 itry = -noctr(1,ii) +
k -1
1697 IF (itry.EQ.iflag) go to 90
1698 xhalf = .5*(xoctr(1,
j) +xoctr(2,
j))
1700 xlow = xsgn*xoctr(1,
j) +(1. -xsgn)*xhalf
1701 xhigh = xsgn*xhalf +(1. -xsgn)*xoctr(2,
j)
1702 yhalf = .5*(yoctr(1,
j) +yoctr(2,
j))
1703 ysgn = isign(1,2*mod(
k-1,4)-3)
1704 ylow = .5*((1. -ysgn)*yoctr(1,
j) +(1. +ysgn)*yhalf)
1705 yhigh = .5*((1. -ysgn)*yhalf +(1. +ysgn)*yoctr(2,
j))
1706 zhalf = .5*(zoctr(1,
j) +zoctr(2,
j))
1707 zsgn = isign(1,2*
k-9)
1708 zlow = .5*((1. -zsgn)*zoctr(1,
j) +(1. +zsgn)*zhalf)
1709 zhigh = .5*((1. -zsgn)*zhalf +(1. +zsgn)*zoctr(2,
j))
1710 IF (xl.GT.xhigh.OR.xu.LT.xlow.OR. &
1711 yl.GT.yhigh.OR.yu.LT.ylow.OR. &
1712 zl.GT.zhigh.OR.zu.LT.zlow) go to 90
1713 IF (noctr(1,itry).GE.0) go to 80
1715 IF (
ic.GT.mxtest) go to 210
1724 80
IF (noctr(1,itry).EQ.0) go to 90
1726 82
IF (idone(jj).EQ.0.AND.ncode.EQ.0) go to 85
1727 IF (idone(jj).EQ.2) go to 85
1728 IF (n1.EQ.0) go to 83
1729 det = c1*(
x(1,jj) -xpt) -c2*(
x(2,jj) -ypt) &
1731 IF (det.LT.tolpt) go to 85
1732 83
dist = (xpt -
x(1,jj))**2 +(ypt -
x(2,jj))**2 &
1734 IF (
dist.GE.dismin) go to 85
1738 IF (jj.NE.0) go to 82
1747 IF (
ic.GT.0) go to 92
1758 xoctr(1,
j) = xhold(1,
j)
1759 xoctr(2,
j) = xhold(2,
j)
1760 yoctr(1,
j) = yhold(1,
j)
1761 yoctr(2,
j) = yhold(2,
j)
1762 zoctr(1,
j) = zhold(1,
j)
1763 zoctr(2,
j) = zhold(2,
j)
1766 200
WRITE (6,600) xpt,ypt,zpt
1770 600
FORMAT(5
x,
'POINT XPT = ',f12.4,
' YPT = ',f12.4,
' ZPT = ',f12.4/ &
1771 5
x,
'LIES OUTSIDE CONVEX HULL')
1772 610
FORMAT(5
x,
'DIMENSION OF NSRCH AND KSRCH ARRAYS EXCEEDED'/ &
1773 5
x,
'IN ROUTINE OCTFND. INCREASE SIZE OF MTEST')
1782 SUBROUTINE octrmv (N,X,NOCTR,NLINK,IDONE,NREF)
1798 INTEGER :: idone(*),nref(*),nlink(*),noctr(2,*)
1799 DOUBLE PRECISION ::
x(3,*)
1801 INTEGER ::
i,
next,npre
1807 IF (
i.EQ.0) go to 200
1809 IF (
next.NE.
n) go to 10
1815 10
IF (
next.EQ.0) go to 210
1818 IF (
next.NE.
n) go to 10
1819 nlink(npre) = nlink(
n)
1823 200
WRITE (6,600)
n,
x(1,
n),
x(2,
n),
x(3,
n)
1825 210
WRITE (6,610)
n,
x(1,
n),
x(2,
n),
x(3,
n)
1827 600
FORMAT(5
x,
'POINT ',i6,
' HAS AN NREF VALUE OF ZERO'/ &
1828 5
x,
'X = ',f10.4,
' Y = ',f10.4,
' Z = ',f10.4)
1829 610
FORMAT(5
x,
'UNABLE TO FIND POINT WITH ADDRESS ',i6,
' IN OCTREE'/ &
1830 5
x,
'X = ',f10.4,
' Y = ',f10.4,
' Z = ',f10.4)
1855 INTEGER :: ndg(2,*),idgp(*),ndgp(*)
1857 INTEGER ::
i,iedg,n1
1861 IF (nedg.EQ.0) go to 200
1862 n1 =
min(ndg(1,nedg),ndg(2,nedg))
1866 IF (
i.NE.nedg) go to 10
1874 IF (
i.EQ.0) go to 210
1875 IF (
i.NE.nedg) go to 10
1885 600
FORMAT(//5
x,
'EDGE REQUESTED HAS ADDRESS ZERO'/ &
1886 5
x,
'PROGRAM STOPPED IN EDGERM')
1887 610
FORMAT(//5
x,
'UNABLE TO FIND EDGE IN LINKED LIST'/ &
1888 5
x,
'PROGRAM STOPPED IN EDGERM')
1897 SUBROUTINE tree (PROP,NLEFT,NRIGHT,NBACK,LISTF,NACPT,NTOT, &
1914 INTEGER :: ncell,ntot
1915 INTEGER :: nbh(4,*),iprot(*)
1916 INTEGER :: nleft(*),nright(*),nback(*),listf(*),nacpt(*)
1917 DOUBLE PRECISION :: prop(*)
1919 INTEGER ::
j,jback,jj,jnext,jprobe,jstart,l,ncnt,
nj
1920 DOUBLE PRECISION :: tol
1931 IF (
j.GT.ncell) go to 90
1932 IF (iprot(
j).EQ.1) go to 10
1933 IF (nbh(1,
j).EQ.0) go to 10
1940 IF (
j.GT.ncell) go to 90
1941 IF (iprot(
j).EQ.1) go to 20
1942 IF (nbh(1,
j).EQ.0) go to 20
1944 30 jnext = nleft(jj)
1945 IF (prop(
j).GT.prop(jj)) jnext = nright(jj)
1946 IF (jnext.EQ.0) go to 40
1949 40
IF (prop(
j).GT.prop(jj)) go to 50
1962 90
IF (ntot.EQ.0)
RETURN
1965 100
j = nright(jprobe)
1966 IF (
j.EQ.0) go to 110
1977 120 jprobe = nleft(
j)
1978 IF (jprobe.GT.0) go to 100
1979 125 jback = nback(
j)
1980 IF (jback.EQ.0) go to 140
1983 IF (
nj.GT.0) go to 115
1985 140
DO 150 l=1,ncell
1998 SUBROUTINE snglar (NNODE,X,XD,NCELL,NDC,NBH,SIGMN,SIGMX, &
1999 sig1,sig2,sig3,cndmin,cndmax)
2055 INTEGER :: ncell,nnode
2056 INTEGER ndc(4,*),nbh(4,*)
2057 DOUBLE PRECISION :: cndmax,cndmin,sigmn,sigmx
2058 DOUBLE PRECISION ::
x(3,*),xd(3,*),sig1(*),sig2(*),sig3(*)
2060 INTEGER :: l,lmax,lmin,npass,n1,n2,n3,n4
2061 DOUBLE PRECISION :: a0,a1,a11,a12,a13,a2,a21,a22,a23,a31,a32,a33, &
2062 b1,b2,b3,b4,b4sq,b5,b5sq,b6,b6sq,cond,ct, &
2063 d,delx1,delx2,delx3,dely1,dely2,dely3, &
2064 delz1,delz2,delz3,det,dsq,dx1,dx2,dx3, &
2065 dy1,dy2,dy3,dz1,dz2,dz3,fac,
p1,p2,
q,qa, &
2066 qroot,
r,
s,siglow,siglrg,slide,slidcb,slidsq, &
2067 st,s1x,s1y,s1z,s2x,s2y,s2z,s3x,s3y,s3z, &
2076 IF (nbh(1,l).EQ.0) go to 100
2077 IF (ndc(4,l).LE.0) go to 100
2086 dx1 =
x(1,n2) -
x(1,n1)
2087 dy1 =
x(2,n2) -
x(2,n1)
2088 dz1 =
x(3,n2) -
x(3,n1)
2089 dx2 =
x(1,n3) -
x(1,n1)
2090 dy2 =
x(2,n3) -
x(2,n1)
2091 dz2 =
x(3,n3) -
x(3,n1)
2092 dx3 =
x(1,n4) -
x(1,n1)
2093 dy3 =
x(2,n4) -
x(2,n1)
2094 dz3 =
x(3,n4) -
x(3,n1)
2095 s1x = dy2*dz3 -dy3*dz2
2096 s1y = dz2*dx3 -dz3*dx2
2097 s1z = dx2*dy3 -dx3*dy2
2098 s2x = dy3*dz1 -dy1*dz3
2099 s2y = dz3*dx1 -dz1*dx3
2100 s2z = dx3*dy1 -dx1*dy3
2101 s3x = dy1*dz2 -dy2*dz1
2102 s3y = dz1*dx2 -dz2*dx1
2103 s3z = dx1*dy2 -dx2*dy1
2104 det = dx1*s1x +dy1*s1y +dz1*s1z
2105 IF (abs(det).LT.tol) go to 300
2111 delx1 = xd(1,n2) -xd(1,n1)
2112 dely1 = xd(2,n2) -xd(2,n1)
2113 delz1 = xd(3,n2) -xd(3,n1)
2114 delx2 = xd(1,n3) -xd(1,n1)
2115 dely2 = xd(2,n3) -xd(2,n1)
2116 delz2 = xd(3,n3) -xd(3,n1)
2117 delx3 = xd(1,n4) -xd(1,n1)
2118 dely3 = xd(2,n4) -xd(2,n1)
2119 delz3 = xd(3,n4) -xd(3,n1)
2120 a11 = (delx1*s1x +delx2*s2x +delx3*s3x)*fac
2121 a12 = (delx1*s1y +delx2*s2y +delx3*s3y)*fac
2122 a13 = (delx1*s1z +delx2*s2z +delx3*s3z)*fac
2123 a21 = (dely1*s1x +dely2*s2x +dely3*s3x)*fac
2124 a22 = (dely1*s1y +dely2*s2y +dely3*s3y)*fac
2125 a23 = (dely1*s1z +dely2*s2z +dely3*s3z)*fac
2126 a31 = (delz1*s1x +delz2*s2x +delz3*s3x)*fac
2127 a32 = (delz1*s1y +delz2*s2y +delz3*s3y)*fac
2128 a33 = (delz1*s1z +delz2*s2z +delz3*s3z)*fac
2132 b1 = (1. +a11)**2 +a21*a21 +a31*a31
2133 b2 = a12*a12 +(1. +a22)**2 +a32*a32
2134 b3 = a13*a13 +a23*a23 +(1. +a33)**2
2135 b4 = a12*(1. +a11) +a21*(1. +a22) +a31*a32
2136 b5 = a13*(1. +a11) +a21*a23 +a31*(1. +a33)
2137 b6 = a12*a13 +a23*(1. +a22) +a32*(1. +a33)
2145 a0 = b1*b6sq +b2*b5sq +b3*b4sq -b1*b2*b3 -2.*b4*b5*b6
2146 a1 = b1*b2 +b1*b3 +b2*b3 -b4sq -b5sq -b6sq
2155 slidsq = slide*slide
2156 slidcb = slidsq*slide
2157 q = a1*third -slidsq
2158 r = .5*third*(a1*a2 -3.*a0) -slidcb
2160 IF (dsq+tol.LT.0.)
WRITE (6,900) a0,a1,a2,
r,
q,dsq
2162 IF (dsq.LT.tol)
THEN
2168 sig1(l) = 2.*
s -slide
2174 ct =
cos(theta*third)
2175 st =
sin(theta*third)
2179 p2 = qroot*st*
sqrt(3.0d0)
2180 sig1(l) = 2.*
p1 -slide
2181 sig2(l) = -
p1 -p2 -slide
2182 sig3(l) = -
p1 +p2 -slide
2184 IF (
min(sig1(l),sig2(l),sig3(l)).LT.0.)
THEN
2185 IF (
min(sig1(l),sig2(l),sig3(l))+tol.LT.0.)
THEN
2188 sig1(l) =
max(0.0d0,sig1(l))
2189 sig2(l) =
max(0.0d0,sig2(l))
2190 sig3(l) =
max(0.0d0,sig3(l))
2193 sig1(l) =
max(tol,sig1(l))
2194 sig2(l) =
max(tol,sig2(l))
2195 sig3(l) =
max(tol,sig3(l))
2196 sig1(l) =
sqrt(sig1(l))
2197 sig2(l) =
sqrt(sig2(l))
2198 sig3(l) =
sqrt(sig3(l))
2199 siglrg =
max(sig1(l),sig2(l),sig3(l))
2200 siglow =
min(sig1(l),sig2(l),sig3(l))
2201 IF (siglow.GT.1.e-8)
THEN
2202 cond = siglrg/siglow
2206 IF (npass.EQ.0)
THEN
2215 IF (siglow.LT.sigmn) lmin = l
2216 IF (siglrg.GT.sigmx) lmax = l
2217 sigmn =
min(sigmn,siglow)
2218 sigmx =
max(sigmx,siglrg)
2219 cndmin =
min(cndmin,cond)
2220 cndmax =
max(cndmax,cond)
2224 300
WRITE (6,600) n1,n2,n3,n4
2227 WRITE (6,900) a0,a1,a2,
r,
q,dsq
2229 320
WRITE (6,620) l,sig1(l),sig2(l),sig3(l)
2231 600
FORMAT(//5
x,
'A ZERO VOLUME TETRAHEDRON HAS BEEN FOUND IN THE '// &
2232 'UNDEFORMED MESH'//5
x,
'THE VERTICES OF THE ELEMENT '// &
2233 'ARE ',i7,
' , ',i7,
' , ',i7,
' AND ',i7// &
2234 5
x,
'PROGRAM STOPPED IN ROUTINE SNGLAR.')
2235 610
FORMAT(//5
x,
'THE DISCRIMINANT ASSOCIATED WITH THE SOLUTION OF', &
2236 5
x,
'CUBIC EQUATION FOR THE SINGULAR VALUES OF THE', &
2237 5
x,
'DEFORMATION MATRIX IS POSITIVE, THUS INDICATING', &
2238 5
x,
'ONLY ONE REAL ROOT. THERE APPEARS TO BE AN ERROR IN', &
2239 5
x,
'THE CODING FOR THE COMPUTATION OF SINGULAR VALUES.', &
2240 5
x,
'PROGRAM STOPPED IN ROUTINE SNGLAR.')
2241 620
FORMAT(//5
x,
'AT LEAST ONE OF THE SINGULAR VALUES SQUARED IS', &
2242 5
x,
'NEGATIVE. THERE THUS APPEARS TO BE AN ERROR IN', &
2243 5
x,
'THE CODING FOR THE COMPUTATION OF SINGULAR VALUES.', &
2244 5
x,
'PROGRAM STOPPED IN ROUTINE SNGLAR.'//5
x,
'CELL ',i7, &
2245 ' SIG1 = ',e13.5,
' SIG2 = ',e13.5,
' SIG3 = ',e13.5)
2246 900
FORMAT(/5
x,
'A POSITIVE DISCRIMINANT HAS BEEN COMPUTED WHEN', &
2247 5
x,
'SOLVING THE CUBIC EQUATION FOR THE SINGULAR VALUES', &
2248 5
x,
'OF THE DEFORMATION MATRIX, THUS INDICATING ONLY', &
2249 5
x,
'ONE REAL ROOT.'// 5
x,
'COEFFICIENTS ARE ', &
2250 'A0 = ',f10.4,
' A1 = ',f10.4,
' A2 = ',f10.4, &
2251 5
x,
'R = ',e13.5,
' Q = ',e13.5,
' DSQ = ',e13.5)
2260 SUBROUTINE densfn (X,NNODE,NFCE,NBFACE,NEDGE,NDG,DENS,RESID,FVCNT)
2275 INTEGER :: nbface,nedge,nnode
2276 INTEGER :: ndg(2,*),nfce(3,*)
2277 DOUBLE PRECISION ::
x(3,*),dens(*)
2278 DOUBLE PRECISION :: fvcnt(*),resid(*)
2280 INTEGER ::
i,it,
j,
k,
n,n1,n2
2281 DOUBLE PRECISION :: difden,
dist,eps,rmax,rmax1
2309 IF (n1.LE.0.OR.n2.LE.0) go to 20
2310 dist =
sqrt((
x(1,n1) -
x(1,n2))**2 +(
x(2,n1) -
x(2,n2))**2 &
2311 +(
x(3,n1) -
x(3,n2))**2)
2312 IF (fvcnt(n1).GE.0.)
THEN
2313 fvcnt(n1) = fvcnt(n1) +1.
2314 dens(n1) = dens(n1) +
dist
2316 IF (fvcnt(n2).LT.0.)
THEN
2317 fvcnt(n1) = fvcnt(n1) -1.
2318 dens(n1) = dens(n1) +
dist
2319 fvcnt(n2) = fvcnt(n2) -1.
2320 dens(n2) = dens(n2) +
dist
2323 IF (fvcnt(n2).GE.0.)
THEN
2324 fvcnt(n2) = fvcnt(n2) +1.
2325 dens(n2) = dens(n2) +
dist
2332 IF (fvcnt(
n).LT.0.) fvcnt(
n) = fvcnt(
n) +1.
2333 IF (abs(fvcnt(
n)).LT.1.) go to 30
2334 fvcnt(
n) = 1./fvcnt(
n)
2335 dens(
n) = abs(fvcnt(
n))*dens(
n)
2355 IF (n1.LE.0.OR.n2.LE.0) go to 75
2356 difden = dens(n1) -dens(n2)
2357 resid(n1) = resid(n1) -difden
2358 resid(n2) = resid(n2) +difden
2365 IF (fvcnt(
n).LT.0.) go to 80
2366 dens(
n) = dens(
n) +eps*fvcnt(
n)*resid(
n)
2367 rmax =
max(rmax,abs(resid(
n)))
2372 IF (it.EQ.1) rmax1 = rmax
2373 IF (rmax.GT.
max(0.0001*rmax1,1.0d0-9).AND.it.LT.1000) go to 65
2374 WRITE (6,900) it,rmax
2375 900
FORMAT(
'IN DENSFN..., IT = ',i4,
' RMAX = ',f10.6)
2392 SUBROUTINE edglen (X,NNODE,ITYP,NEDGE,NDG,DENS,NVCNT)
2407 INTEGER :: nedge,nnode
2408 INTEGER :: ndg(2,*),ityp(*),nvcnt(*)
2409 DOUBLE PRECISION ::
x(3,*),dens(*)
2411 INTEGER ::
i,
n,nhigh,nlow,nmid,npass,ntot,nvmax,nvmin,n1,n2
2412 DOUBLE PRECISION :: densav,dsqav
2428 IF (n1.LE.0.OR.n2.LE.0) go to 20
2429 nvcnt(n1) = nvcnt(n1) +1
2430 nvcnt(n2) = nvcnt(n2) +1
2431 densav = .5*(dens(n1) +dens(n2))
2432 dsqav =
sqrt((
x(1,n1) -
x(1,n2))**2 +(
x(2,n1) -
x(2,n2))**2 &
2433 +(
x(3,n1) -
x(3,n2))**2)
2434 IF (dsqav.LT.0.5*densav) nlow = nlow +1
2435 IF (dsqav.GT.1.5*densav) nhigh = nhigh +1
2436 IF (dsqav.GE.0.5*densav.AND.dsqav.LE.1.5*densav) nmid = nmid +1
2438 ntot = nlow +nmid +nhigh
2439 WRITE (6,600) nlow,nmid,nhigh,ntot
2444 IF (ityp(
n).GE.0)
THEN
2445 IF (npass.EQ.0)
THEN
2450 nvmax =
max(nvmax,nvcnt(
n))
2451 nvmin =
min(nvmin,nvcnt(
n))
2458 WRITE (6,610) nvmin,nvmax
2460 600
FORMAT(/5
x,
'NLOW = ',i6,
' NMID = ',i6,
' NHIGH = ',i6, &
2462 610
FORMAT(/5
x,
'MIN EDGE VALENCE = ',i3,
' MAX EDGE VALENCE = ',i3)
2472 SUBROUTINE coarsn (LC,X,NNODE,NDC,NBH,IPROT,NCELL, &
2473 ityp,xcen,ycen,zcen,vol,rc,rat, &
2474 nvcnt,dens,iflag,nflag,nptet, &
2475 nedge,ndg,idgp,ndgp, &
2476 noctr,ioctr,nlink,xfar,yfar,zfar,idone,nref, &
2477 fac,sig1,sig2,sig3, &
2478 ksrch,nsrch,iring,ntetkp, &
2479 lnbr,ishk,mnbr,kshk,tolv)
2495 INTEGER :: ioctr,lc,ncell,nedge,nnode
2496 INTEGER :: ityp(*),ndc(4,*),nbh(4,*),iprot(*),ndg(2,*)
2497 INTEGER :: idgp(*),ndgp(*),nvcnt(*),nptet(*)
2498 INTEGER :: idone(*),nref(*),nlink(*),noctr(2,*)
2499 INTEGER :: iflag(*),nflag(*),nsrch(*),ksrch(*)
2500 INTEGER :: iring(*),ntetkp(*),lnbr(*),ishk(*),mnbr(*),kshk(*)
2501 DOUBLE PRECISION :: tolv
2502 DOUBLE PRECISION ::
x(3,*),dens(*),fac(*)
2503 DOUBLE PRECISION :: vol(*),xcen(*),ycen(*),zcen(*),rc(*),rat(*)
2504 DOUBLE PRECISION :: sig1(*),sig2(*),sig3(*)
2505 DOUBLE PRECISION :: xfar(2),yfar(2),zfar(2)
2507 INTEGER ::
i,
j,l,
n,
ncol,nfail,npts,ntest,ntet,n1,n2,n3,n4
2508 DOUBLE PRECISION :: densav,
dist,sigmax,sigmin,excess,expand
2522 IF (n1.LE.0.OR.n2.LE.0) go to 20
2523 nvcnt(n1) = nvcnt(n1) +1
2524 nvcnt(n2) = nvcnt(n2) +1
2531 IF (nbh(1,
j).LE.0) go to 30
2532 IF (ndc(4,
j).EQ.-1) go to 30
2534 sigmax =
max(sig1(
j),sig2(
j),sig3(
j))
2535 sigmin =
min(sig1(
j),sig2(
j),sig3(
j))
2536 expand =
min(1.0d0,sigmax/sigmin -1.)
2538 excess =
max(0.0d0,rat(
j)-20.)
2539 expand =
min(1.0d0,excess/100.)
2545 fac(n1) = fac(n1) +expand
2546 fac(n2) = fac(n2) +expand
2547 fac(n3) = fac(n3) +expand
2548 fac(n4) = fac(n4) +expand
2551 IF (ityp(
n).GE.0)
THEN
2552 dens(
n) = dens(
n)*(1. +fac(
n)/float(nvcnt(
n)))
2562 IF (n1.LE.0.OR.n2.LE.0) go to 40
2563 densav = .5*(dens(n1) +dens(n2))
2564 dist =
sqrt((
x(1,n1) -
x(1,n2))**2 +(
x(2,n1) -
x(2,n2))**2 &
2565 +(
x(3,n1) -
x(3,n2))**2)
2566 IF (
dist.GE.0.4*densav) go to 40
2570 ityp,xcen,ycen,zcen,vol,rc,rat, &
2571 nvcnt,iflag,nflag,nptet, &
2573 noctr,ioctr,nlink,xfar,yfar,zfar,idone,nref, &
2574 ksrch,nsrch,iring,ntetkp, &
2575 lnbr,ishk,mnbr,kshk,tolv)
2583 IF (ityp(
n).LT.0) go to 50
2588 IF (nbh(1,l).EQ.0) go to 60
2589 IF (ndc(4,l).EQ.-1) go to 60
2592 WRITE (6,730) npts,
ncol,ntet
2594 730
FORMAT(/
'MESH COARSENING COMPLETE'/ &
2595 5
x,i7,
' TOTAL MESH POINTS'/ &
2596 5
x,i7,
' FIELD POINTS REMOVED'/ &
2597 5
x,i7,
' MESH CELLS'/)
2609 SUBROUTINE radrat (X,NNODE,NDC,NBH,IPROT,NCELL,NFCE,NBFACE, &
2610 ityp,ipoint,vol,rc,rat)
2626 INTEGER :: nbface,ncell,nnode
2627 INTEGER :: nbh(4,*),iprot(*),ityp(*),ipoint(*)
2628 INTEGER :: nfce(3,*),ndc(4,*)
2629 DOUBLE PRECISION ::
x(3,*)
2630 DOUBLE PRECISION :: vol(*),rc(*),rat(*)
2632 INTEGER ::
j,jamin,jmax,jmin,
k,l,
n,nbfc,ncnt,nsurpt,nvolpt, &
2634 DOUBLE PRECISION ::
v(4),ar(4)
2635 DOUBLE PRECISION :: angl,angl1,angl2,angl3,angmax,angmin,dismin, &
2636 farea,fmax,fmin,
q,qmax,qmin,vmax,vmin, &
2637 ratavr,ratmax,ratmin,rcmax,rcmin,sigma,sigsq
2644 IF (iprot(l).EQ.1) go to 20
2645 IF (nbh(1,l).EQ.0) go to 20
2646 IF (ndc(4,l).LE.0) go to 20
2648 ratavr = ratavr +rat(l)
2649 IF (ncnt.GT.1) go to 10
2657 10 vmin =
min(vol(l),vmin)
2658 vmax =
max(vol(l),vmax)
2659 rcmin =
min(rcmin,rc(l))
2660 rcmax =
max(rcmax,rc(l))
2661 ratmin =
min(ratmin,rat(l))
2662 ratmax =
max(ratmax,rat(l))
2664 ratavr = ratavr/float(ncnt)
2673 IF (iprot(l).EQ.1) go to 30
2674 IF (nbh(1,l).EQ.0) go to 30
2675 IF (ndc(4,l).LE.0) go to 30
2676 sigsq = sigsq +(rat(l) -ratavr)**2
2678 sigma =
sqrt(sigsq)/float(ncnt)
2679 WRITE (6,600) vmin,vmax,rcmin,rcmax,ratmin,ratmax,ratavr,sigma
2681 CALL
tetang(
x,ndc,nbh,iprot,ncell)
2699 IF (n1.LT.0) go to 75
2702 CALL
fangle(
j,
x,nfce,angl1,angl2,angl3,
q)
2704 angl =
min(angl1,angl2,angl3)
2705 IF (angl.GT.angmin) go to 60
2708 60 angmax =
max(angl1,angl2,angl3,angmax)
2712 IF (fmax.LT.0.) fmax = farea
2713 IF (fmin.LT.0.) fmin = farea
2714 IF (farea.LT.fmax) go to 70
2717 70
IF (farea.GT.fmin) go to 75
2721 WRITE (6,610) angmin,angmax,qmin,qmax
2729 IF (ityp(
n).GE.0) nvolpt = nvolpt +1
2734 IF (ipoint(
n).EQ.0)
THEN
2742 WRITE (6,620) nsurpt,nbfc,nvolpt,ncnt
2744 600
FORMAT(//5
x,
'*************************************************'/ &
2746 5
x,
'* STATISTICS OF MESH TETRAHEDRA *'/ &
2748 5
x,
'* MINIMUM VOLUME = ',e13.5,
' *'/ &
2749 5
x,
'* MAXIMUM VOLUME = ',e13.5,
' *'/ &
2751 5
x,
'* MINIMUM CIRCUMRADIUS = ',e13.5,
' *'/ &
2752 5
x,
'* MAXIMUM CIRCUMRADIUS = ',e13.5,
' *'/ &
2754 5
x,
'* MINIMUM RADIUS RATIO = ',e13.5,
' *'/ &
2755 5
x,
'* MAXIMUM RADIUS RATIO = ',e13.5,
' *'/ &
2756 5
x,
'* AVERAGE RADIUS RATIO = ',e13.5,
' *'/ &
2757 5
x,
'* STANDARD DEVIATION = ',e13.5,
' *'/ &
2759 5
x,
'*************************************************')
2760 610
FORMAT( 5
x,
'*************************************************'/ &
2762 5
x,
'* STATISTICS OF BOUNDARY SURFACE TRIANGLES *'/ &
2764 5
x,
'* MINIMUM BOUNDARY FACE ANGLE = ',f6.2,
' *'/ &
2765 5
x,
'* MAXIMUM BOUNDARY FACE ANGLE = ',f6.2,
' *'/ &
2766 5
x,
'* MINIMUM BOUNDARY RADIUS RATIO = ',f6.2,
' *'/ &
2767 5
x,
'* MAXIMUM BOUNDARY RADIUS RATIO = ',f6.2,
' *'/ &
2769 5
x,
'*************************************************')
2770 620
FORMAT( 5
x,
'*************************************************'/ &
2772 5
x,
'* TOTAL NUMBER OF SURFACE POINTS = ',i6,
' *'/ &
2773 5
x,
'* TOTAL NUMBER OF SURFACE FACES = ',i6,
' *'/ &
2774 5
x,
'* TOTAL NUMBER OF MESH POINTS = ',i6,
' *'/ &
2775 5
x,
'* TOTAL NUMBER OF TETRAHEDRA = ',i6,
' *'/ &
2777 5
x,
'*************************************************')
2788 SUBROUTINE volput (X,ITYP,NBPTS,NNODE,NDC,NBH,IPROT,NCELL, &
2789 ndg,idgp,ndgp,nedge, &
2790 vol,xcen,ycen,zcen,rc,rat,dens,nptet,nacpt, &
2791 idone,nref,nlink,noctr,ioctr,xfar,yfar,zfar, &
2792 xoctr,yoctr,zoctr,xhold,yhold,zhold, &
2793 xkeep,ykeep,zkeep,ksrch,nsrch, &
2794 ipoint,npoint,iflag,nflag,nfill,newcel,ntri, &
2795 ncav,nshake,newnbh,nold,ncavfc,ikeep,ldel, &
2796 nedgrm,xc,yc,zc,
v,rad,rcrin,lnkup,lnkdn, &
2797 jlast,jfirst,ntrack,volmin,rcmx,tolv)
2812 INTEGER :: ioctr,jfirst,jlast,nbpts,ncell,nedge,nnode,ntrack
2813 INTEGER :: ndc(4,*),ndg(2,*),nbh(4,*),iprot(*),idgp(*),ndgp(*)
2814 INTEGER :: ityp(*),nptet(*),nacpt(*)
2815 INTEGER :: idone(*),nref(*),nlink(*),noctr(2,*)
2816 INTEGER :: npoint(*),ipoint(*)
2817 INTEGER :: iflag(*),nflag(*)
2818 INTEGER :: lnkup(*),lnkdn(*)
2819 INTEGER :: nsrch(*),ksrch(*)
2820 INTEGER :: ntri(3,*),nfill(*),newnbh(4,*),nold(*), &
2821 newcel(*),nshake(*),ncav(4,*)
2822 INTEGER :: nedgrm(*)
2823 INTEGER :: ldel(*),ncavfc(3,*),ikeep(*)
2824 DOUBLE PRECISION :: volmin,rcmx,tolv
2825 DOUBLE PRECISION ::
x(3,*),dens(*)
2826 DOUBLE PRECISION :: vol(*),xcen(*),ycen(*),zcen(*),rc(*),rat(*)
2827 DOUBLE PRECISION :: xoctr(2,*),yoctr(2,*),zoctr(2,*), &
2828 xhold(2,*),yhold(2,*),zhold(2,*), &
2829 xkeep(2),ykeep(2),zkeep(2),xfar(2), &
2831 DOUBLE PRECISION :: xc(*),yc(*),zc(*),
v(*),rad(*),rcrin(*)
2833 INTEGER ::
j,jpast,jpre,
k,k1,k2,l,lbrk,lnbh,
n,nclose,ncyc,ndiff, &
2834 nfail,ninp,npass,npts,nsrfpt,nstart,ntet,n1,n2,n3,n4,n5
2835 DOUBLE PRECISION :: dent,dismin,rnx,rny,rnz,scldis,vtet1,vtet2, &
2858 IF (
j.GT.ncell) go to 15
2859 IF (iprot(
j).EQ.1) go to 10
2860 IF (nbh(1,
j).LE.0) go to 10
2866 dent = dens(n1) +dens(n2) +dens(n3) +dens(n4)
2867 IF (rc(
j).LT.0.22*dent) go to 10
2870 IF (jlast.NE.0) lnkdn(jlast) =
j
2874 IF (jfirst.NE.0) go to 10
2877 15
WRITE (6,910) ntrack
2878 910
FORMAT(
'IN VOLPUT, NTRACK = ',i6)
2879 IF (ntrack.EQ.0) go to 70
2885 IF (nacpt(
j).EQ.0) go to 215
2886 IF (iprot(
j).EQ.1) go to 200
2890 IF (xpt.LT.xfar(1).OR.xpt.GT.xfar(2)) go to 45
2891 IF (ypt.LT.yfar(1).OR.ypt.GT.yfar(2)) go to 45
2892 IF (zpt.LT.zfar(1).OR.zpt.GT.zfar(2)) go to 45
2895 IF (iprot(lnbh).EQ.0) go to 32
2897 CALL
lock(
j,lnbh,n1,n2,n3,n4,n5,k1,k2,ndc,nbh)
2899 rnx =
cofact(
x(2,n1),
x(2,n2),
x(2,n3),
x(3,n1),
x(3,n2),
x(3,n3))
2900 rny =
cofact(
x(3,n1),
x(3,n2),
x(3,n3),
x(1,n1),
x(1,n2),
x(1,n3))
2901 rnz =
cofact(
x(1,n1),
x(1,n2),
x(1,n3),
x(2,n1),
x(2,n2),
x(2,n3))
2902 vtet1 = rnx*(
x(1,n4) -
x(1,n1)) +rny*(
x(2,n4) -
x(2,n1)) &
2903 +rnz*(
x(3,n4) -
x(3,n1))
2904 vtet2 = rnx*(xpt -
x(1,n1)) +rny*(ypt -
x(2,n1)) &
2906 IF (vtet1*vtet2.LT.0.) go to 45
2909 IF (nnode.GT.mxnode) go to 230
2919 CALL
octfnd(0,nclose,
x(1,
n),
x(2,
n),
x(3,
n),0,0,0,dismin, &
2920 x,noctr,nlink,xfar,yfar,zfar,idone, &
2921 xoctr,yoctr,zoctr,xhold,yhold,zhold, &
2922 xkeep,ykeep,zkeep,ksrch,nsrch)
2924 IF (dismin.LT.1.e-15) go to 40
2925 IF (nclose.LE.nsrfpt) go to 35
2926 IF (nclose.LE.nstart) go to 35
2927 IF (
sqrt(dismin).LT.scldis*dens(nclose)) go to 40
2929 lbrk = nptet(nclose)
2931 CALL
insert(
n,nclose,lbrk,
j,nfail, &
2932 x,ndc,nbh,iprot,ncell,ndg,idgp,ndgp,nedge, &
2933 vol,xcen,ycen,zcen,rc,rat,dens,nptet,nacpt, &
2934 ipoint,npoint,iflag,nflag,nfill,newcel,ntri, &
2935 ncav,nshake,newnbh,nold,ncavfc,ikeep,ldel, &
2936 nedgrm,xc,yc,zc,
v,rad,rcrin,lnkup,lnkdn, &
2937 jlast,jfirst,ntrack,volmin,rcmx,tolv)
2939 IF (nfail.EQ.5.AND.npass.EQ.0)
THEN
2944 xpt = .25*(
x(1,n1) +
x(1,n2) +
x(1,n3) +
x(1,n4))
2945 ypt = .25*(
x(2,n1) +
x(2,n2) +
x(2,n3) +
x(2,n4))
2946 zpt = .25*(
x(3,n1) +
x(3,n2) +
x(3,n3) +
x(3,n4))
2950 IF (nfail.GE.1) go to 40
2960 CALL
octfil(
n,
x,noctr,ioctr,nlink,nref,xfar,yfar,zfar)
2970 IF (jpre.NE.0) lnkdn(jpre) = jpast
2971 IF (jpast.NE.0) lnkup(jpast) = jpre
2972 IF (jpre.EQ.0) jfirst = jpast
2973 IF (jpast.EQ.0) jlast = jpre
2975 50
IF (mod(ncyc,10000).EQ.0)
WRITE (6,710) ncyc,ninp,ntrack
2976 IF (jfirst.NE.0) go to 30
2980 70
WRITE (6,710) ncyc,ninp,ntrack
2984 IF (ityp(
n).GE.0) npts = npts +1
2987 ndiff = nnode -nstart
2990 IF (iprot(l).EQ.1) go to 80
2991 IF (nbh(1,l).EQ.0) go to 80
2994 WRITE (6,720) npts,ndiff,ntet
2998 210
WRITE (6,610) nfail,
n,
x(1,
n),
x(2,
n),
x(3,
n)
3004 710
FORMAT(/
' NCYC = ',i7,
' POINTS INSERTED = ',i7,
' NTRACK = ',i7)
3005 720
FORMAT(/
'ADAPTIVE REFINEMENT COMPLETE'/ &
3006 5
x,i7,
' TOTAL MESH POINTS'/ &
3007 5
x,i7,
' FIELD POINTS INSERTED'/ &
3008 5
x,i7,
' MESH CELLS'/)
3009 600
FORMAT(//
'A PROTECTED TETRAHEDRON IS AMONG THOSE TO BE REFINED')
3010 610
FORMAT(//5
x,
'NFAIL = ',i6,
' N ',i6,
' X,Y,Z ',3(2
x,f6.2)/ &
3011 5
x,
'PROGRAM STOPPED IN VOLPUT')
3012 615
FORMAT(//5
x,
'LABEL NACPT IS ZERO FOR CELL FROM ACTIVE LIST'/ &
3013 5
x,
'PROGRAM STOPPED IN VOLPUT')
3014 630
FORMAT(//
'NUMBER OF POINTS INSERTED EXCEEDS DIMENSION OF ARRAY X'/ &
3015 'INCREASE SIZE OF MNODE.'/ &
3016 'PROGRAM STOPPED IN ROUTINE VOLPUT.')
3027 x,ityp,nnode,ndc,nbh,iprot,ncell, &
3028 ndg,idgp,ndgp,nedge, &
3029 vol,xcen,ycen,zcen,rc,rat,dens,nptet,nacpt, &
3030 idone,nref,nlink,noctr,ioctr,xfar,yfar,zfar, &
3031 xoctr,yoctr,zoctr,xhold,yhold,zhold, &
3032 xkeep,ykeep,zkeep,ksrch,nsrch, &
3033 ipoint,npoint,iflag,nflag,nfill,newcel,ntri, &
3034 ncav,nshake,newnbh,nold,ncavfc,ikeep,ldel, &
3035 nedgrm,xc,yc,zc,
v,rad,rcrin,lnkup,lnkdn, &
3052 INTEGER :: ioctr,
j,ncell,nedge,nfail,nnode
3053 INTEGER :: ndc(4,*),nbh(4,*),iprot(*),ndg(2,*),idgp(*),ndgp(*)
3054 INTEGER :: ntri(3,*),nfill(*),newnbh(4,*),nold(*),newcel(*), &
3056 INTEGER :: nedgrm(*),ldel(*),ncavfc(3,*),ikeep(*)
3057 INTEGER :: ityp(*),nptet(*),nacpt(*)
3058 INTEGER :: idone(*),nref(*),nlink(*),noctr(2,*)
3059 INTEGER :: npoint(*),ipoint(*)
3060 INTEGER :: iflag(*),nflag(*)
3061 INTEGER :: lnkup(*),lnkdn(*)
3062 INTEGER :: nsrch(*),ksrch(*)
3063 DOUBLE PRECISION :: volmin,rcmx,tolv
3064 DOUBLE PRECISION ::
x(3,*),dens(*)
3065 DOUBLE PRECISION :: vol(*),xcen(*),ycen(*),zcen(*),rc(*),rat(*)
3066 DOUBLE PRECISION :: xoctr(2,*),yoctr(2,*),zoctr(2,*), &
3067 xhold(2,*),yhold(2,*),zhold(2,*), &
3068 xfar(2),yfar(2),zfar(2),xkeep(2), &
3070 DOUBLE PRECISION :: xc(*),yc(*),zc(*),
v(*),rad(*),rcrin(*)
3072 INTEGER :: jfirst,jlast,
k,k1,k2,lbrk,lnbh,
n,nclose,npass,ntrack, &
3074 DOUBLE PRECISION :: dismin,rnx,rny,rnz,vtet1,vtet2,xpt,ypt,zpt
3082 IF (xpt.LT.xfar(1).OR.xpt.GT.xfar(2)) go to 50
3083 IF (ypt.LT.yfar(1).OR.ypt.GT.yfar(2)) go to 50
3084 IF (zpt.LT.zfar(1).OR.zpt.GT.zfar(2)) go to 50
3087 IF (iprot(lnbh).EQ.0) go to 32
3089 CALL
lock(
j,lnbh,n1,n2,n3,n4,n5,k1,k2,ndc,nbh)
3091 rnx =
cofact(
x(2,n1),
x(2,n2),
x(2,n3),
x(3,n1),
x(3,n2),
x(3,n3))
3092 rny =
cofact(
x(3,n1),
x(3,n2),
x(3,n3),
x(1,n1),
x(1,n2),
x(1,n3))
3093 rnz =
cofact(
x(1,n1),
x(1,n2),
x(1,n3),
x(2,n1),
x(2,n2),
x(2,n3))
3094 vtet1 = rnx*(
x(1,n4) -
x(1,n1)) +rny*(
x(2,n4) -
x(2,n1)) &
3095 +rnz*(
x(3,n4) -
x(3,n1))
3096 vtet2 = rnx*(xpt -
x(1,n1)) +rny*(ypt -
x(2,n1)) &
3098 IF (vtet1*vtet2.LT.0.) go to 50
3101 IF (nnode.GT.mxnode) go to 200
3111 CALL
octfnd(0,nclose,
x(1,
n),
x(2,
n),
x(3,
n),0,0,0,dismin, &
3112 x,noctr,nlink,xfar,yfar,zfar,idone, &
3113 xoctr,yoctr,zoctr,xhold,yhold,zhold, &
3114 xkeep,ykeep,zkeep,ksrch,nsrch)
3116 IF (dismin.LT.1.e-15) go to 40
3118 lbrk = nptet(nclose)
3123 CALL
insert(
n,nclose,lbrk,
j,nfail, &
3124 x,ndc,nbh,iprot,ncell,ndg,idgp,ndgp,nedge, &
3125 vol,xcen,ycen,zcen,rc,rat,dens,nptet,nacpt, &
3126 ipoint,npoint,iflag,nflag,nfill,newcel,ntri, &
3127 ncav,nshake,newnbh,nold,ncavfc,ikeep,ldel, &
3128 nedgrm,xc,yc,zc,
v,rad,rcrin,lnkup,lnkdn, &
3129 jlast,jfirst,ntrack,volmin,rcmx,tolv)
3131 IF (nfail.EQ.5.AND.npass.EQ.0)
THEN
3136 xpt = .25*(
x(1,n1) +
x(1,n2) +
x(1,n3) +
x(1,n4))
3137 ypt = .25*(
x(2,n1) +
x(2,n2) +
x(2,n3) +
x(2,n4))
3138 zpt = .25*(
x(3,n1) +
x(3,n2) +
x(3,n3) +
x(3,n4))
3142 IF (nfail.GE.1)
THEN
3149 CALL
octfil(
n,
x,noctr,ioctr,nlink,nref,xfar,yfar,zfar)
3160 600
FORMAT(//
'NUMBER OF POINTS INSERTED EXCEEDS DIMENSION OF ARRAY X'/ &
3161 'INCREASE SIZE OF MNODE.'/ &
3162 'PROGRAM STOPPED IN ROUTINE PUTPNT.')
3173 SUBROUTINE insert (NP,NCLOSE,LBRK,JCONT,NFAIL, &
3174 x,ndc,nbh,iprot,ncell,ndg,idgp,ndgp,nedge, &
3175 vol,xcen,ycen,zcen,rc,rat,dens,nptet,nacpt, &
3176 ipoint,npoint,iflag,nflag,nfill,newcel,ntri, &
3177 ncav,nshake,newnbh,nold,ncavfc,ikeep,ldel, &
3178 nedgrm,xc,yc,zc,
v,rad,rcrin,lnkup,lnkdn, &
3179 jlast,jfirst,ntrack,volmin,rcmx,tolv)
3197 INTEGER :: jcont,jfirst,jlast,lbrk,ncell,nclose,nedge,nfail,np, &
3199 INTEGER :: ndc(4,*),nbh(4,*),iprot(*),ndg(2,*),idgp(*),ndgp(*)
3200 INTEGER :: nptet(*),nacpt(*)
3201 INTEGER :: npoint(*),ipoint(*),iflag(*),nflag(*)
3202 INTEGER :: ntri(3,*),nfill(*),newnbh(4,*),nold(*), &
3203 newcel(*),nshake(*),ncav(4,*)
3204 INTEGER :: nedgrm(*)
3205 INTEGER :: lnkup(*),lnkdn(*)
3206 INTEGER :: ldel(*),ncavfc(3,*),ikeep(*)
3207 DOUBLE PRECISION :: volmin,rcmx,tolv
3208 DOUBLE PRECISION ::
x(3,*),dens(*)
3209 DOUBLE PRECISION :: vol(*),xcen(*),ycen(*),zcen(*),rc(*),rat(*)
3210 DOUBLE PRECISION :: xc(*),yc(*),zc(*),
v(*),rad(*),rcrin(*)
3212 INTEGER ::
i,iedgrm,ismall,
j,jpast,jpre,
k,l,lcont,
m,mm,nchk, &
3213 ncnt,ncpnt,ndel,nedg,ntot1,ntot2,n1,n2,n3,n4
3214 DOUBLE PRECISION :: area,a1,a2,a3,b1,b2,b3,dent, &
3215 radmax,ratmax,sum1,sum2,vdiff,xpt,ypt,zpt
3226 CALL
tetloc(np,nclose,lbrk,lcont,nfail, &
3227 x,ndc,nbh,iprot,dens,vol, &
3228 iflag,nflag,nfill,newcel,tolv)
3230 IF (nfail.NE.0)
RETURN
3235 CALL
cavity(np,lcont,ndel,ldel, &
3236 x,ndc,nbh,iprot,iflag,nflag, &
3237 xcen,ycen,zcen,rc,nfill,newcel,tolv)
3242 CALL
cavbnd(np,lcont,ndel,ldel,ncnt, &
3243 x,ndc,nbh,vol,iflag,nflag, &
3244 nfill,newcel,ntri,newnbh,nold,tolv)
3252 radmax =
max(radmax,rc(ldel(
k)))
3253 ratmax =
max(ratmax,rat(ldel(
k)))
3254 IF (ldel(
k).EQ.jcont) nchk = 1
3256 IF (nchk.EQ.0) go to 295
3268 CALL
circum(
x,n1,n2,n3,n4,xc(
k),yc(
k),zc(
k),
v(
k),rad(
k), &
3271 IF (ismall.EQ.1) go to 300
3272 IF (
v(
k).LT.tolv) go to 310
3273 IF (rad(
k).GT.radmax) go to 305
3275 rcrin(
k) = rad(
k)*area/
v(
k)
3276 IF (rcrin(
k).GT.ratmax) go to 305
3307 sum2 = sum2 +vol(ldel(
k))
3309 IF (sum1.GT.sum2*(1.+tolv)) go to 350
3313 CALL
recon(ldel,ndel,ncnt, &
3314 ndc,nbh,iprot,ndg,idgp,ndgp,nflag, &
3315 ntri,ncavfc,ikeep,nedgrm,iedgrm)
3321 IF (
k.GT.ndel) go to 35
3325 IF (ncell.GT.mxcell) go to 370
3337 IF (nflag(ncpnt).EQ.1) go to 44
3339 iflag(ntot1) = ncpnt
3349 IF (nflag(ncpnt).EQ.1) go to 48
3351 iflag(ntot2) = ncpnt
3357 IF (ntot1.NE.ntot2) go to 390
3366 IF (nbh(
m,l).EQ.nold(
k)) go to 60
3373 CALL
cavedg(ncnt,nedg,ipoint,npoint, &
3374 ntri,ncav,newnbh,newcel)
3378 CALL
datsrf(np,ncnt,nedg,ndg,idgp,ndgp,nedge, &
3379 ipoint,ncav,nedgrm,iedgrm)
3386 70 nbh(
i,ldel(
j)) = 0
3398 nbh(nshake(
k),newnbh(1,
k)) = newcel(
k)
3402 nbh(
i,newcel(
k)) = newnbh(
i,
k)
3404 ndc(4,newcel(
k)) = np
3405 ndc(1,newcel(
k)) =
min(ntri(1,
k),ntri(2,
k),ntri(3,
k))
3406 ndc(2,newcel(
k)) =
max(ntri(1,
k),ntri(2,
k),ntri(3,
k))
3407 ndc(3,newcel(
k)) = ntri(1,
k) +ntri(2,
k) +ntri(3,
k) &
3408 -ndc(1,newcel(
k)) -ndc(2,newcel(
k))
3409 nptet(ntri(1,
k)) = newcel(
k)
3410 nptet(ntri(2,
k)) = newcel(
k)
3411 nptet(ntri(3,
k)) = newcel(
k)
3412 xcen(newcel(
k)) = xc(
k)
3413 ycen(newcel(
k)) = yc(
k)
3414 zcen(newcel(
k)) = zc(
k)
3415 vol(newcel(
k)) =
v(
k)
3416 rc(newcel(
k)) = rad(
k)
3417 rat(newcel(
k)) = rcrin(
k)
3418 volmin =
min(volmin,
v(
k))
3419 rcmx =
max(rcmx,rad(
k))
3421 nptet(np) = newcel(1)
3426 IF (nacpt(ldel(
j)).EQ.0) go to 140
3427 jpre = lnkup(ldel(
j))
3428 jpast = lnkdn(ldel(
j))
3429 IF (jpre.NE.0) lnkdn(jpre) = jpast
3430 IF (jpast.NE.0) lnkup(jpast) = jpre
3431 IF (jpre.EQ.0) jfirst = jpast
3432 IF (jpast.EQ.0) jlast = jpre
3437 n1 = ndc(1,newcel(
k))
3438 n2 = ndc(2,newcel(
k))
3439 n3 = ndc(3,newcel(
k))
3440 n4 = ndc(4,newcel(
k))
3441 nacpt(newcel(
k)) = 0
3442 dent = dens(n1) +dens(n2) +dens(n3) +dens(n4)
3443 IF (rc(newcel(
k)).LT.0.22*dent) go to 150
3444 nacpt(newcel(
k)) = 1
3446 IF (jlast.NE.0) lnkdn(jlast) = newcel(
k)
3447 lnkdn(newcel(
k)) = 0
3448 lnkup(newcel(
k)) = jlast
3450 IF (jfirst.EQ.0) jfirst = newcel(
k)
3462 WRITE (6,888) n1,
x(1,n1),
x(2,n1),
x(3,n1), &
3463 n2,
x(1,n2),
x(2,n2),
x(3,n2), &
3464 n3,
x(1,n3),
x(2,n3),
x(3,n3), &
3465 n4,
x(1,n4),
x(2,n4),
x(3,n4),
v(
k),tolv
3466 888
FORMAT(
'N1 ',i5,
' X = ',f8.4,
' Y = ',f8.4,
' Z = ',f8.4/ &
3467 'N2 ',i5,
' X = ',f8.4,
' Y = ',f8.4,
' Z = ',f8.4/ &
3468 'N3 ',i5,
' X = ',f8.4,
' Y = ',f8.4,
' Z = ',f8.4/ &
3469 'N4 ',i5,
' X = ',f8.4,
' Y = ',f8.4,
' Z = ',f8.4/ &
3470 'VOLUME = ',e13.5,
' TOLV = ',e13.5)
3476 WRITE (6,610)
v(
k),tolv
3507 WRITE (6,690) ntot1,ntot2
3509 595
FORMAT(5
x,
'CAVITY DOES NOT CONTAIN ORIGINATING TETRAHEDRON')
3510 600
FORMAT(//5
x,
'AT LEAST ONE NEW TETRAHEDRON HAS TOO SMALL A VOLUME'/ &
3511 5
x,
'ADDRESS OF INSERTED POINT IS ',i6)
3512 605
FORMAT(/5
x,
'RADMAX ',e13.5,
' RC ',e13.5, &
3513 ' RATMAX ',e13.5,
' RAT ',e13.5)
3514 610
FORMAT(5
x,
'VOLUME OF A NEW TETRAHEDRON IS LESS THAN TOLV'/ &
3515 5
x,
'VOLUME = ',e13.5,
' TOLV = ',e13.5/)
3516 620
FORMAT(5
x,
'NEW POINT CREATES A TETRAHEDRON TOO CLOSE TO BOUNDARY')
3517 650
FORMAT(5
x,
'VOLUME VISIBILITY CHECK FAILED, VDIFF = ',e13.5)
3518 670
FORMAT(//5
x,
'DIMENSION OF NDC EXCEEDED IN ROUTINE INSERT'/ &
3519 5
x,
'INCREASE SIZE OF MCELL.')
3520 660
FORMAT(//5
x,
'UNABLE TO FIND EDGE ADDRESS FOR A NEW TETRAHEDRON'/ &
3521 5
x,
'PROGRAM STOPPED IN ROUTINE INSERT')
3522 680
FORMAT(//5
x,
'UNABLE TO FIND A CONTIGUITY BETWEEN A NEW'/ &
3523 5
x,
'TETRAHEDRON AND A NON-CAVITY TETRAHEDRON'/ &
3524 5
x,
'PROGRAM STOPPED IN ROUTINE INSERT')
3525 690
FORMAT(//5
x,
'NUMBER OF CAVITY POINTS = ',i6,
' IS DIFFERENT FROM'/ &
3526 5
x,
'THE NUMBER OF NEW CELL BOUNDARY POINTS = ',i6// &
3527 5
x,
'PROBABLE CAUSE IS A FAILURE IN THE TOLERANCE FOR'/ &
3528 5
x,
'THE DELAUNAY SPHERE TEST IN ROUTINE CAVITY.'/ &
3529 5
x,
'PROGRAM STOPPED IN ROUTINE INSERT.')
3540 SUBROUTINE tetloc (NP,NCLOSE,LBRK,LCONT,NFAIL, &
3541 x,ndc,nbh,iprot,dens,vol, &
3542 iflag,nflag,nfill,newcel,tolv)
3558 INTEGER :: lbrk,lcont,nclose,nfail,np
3559 INTEGER :: ndc(4,*),nbh(4,*),iprot(*)
3560 INTEGER :: iflag(*),nflag(*)
3561 INTEGER :: nfill(*),newcel(*)
3562 DOUBLE PRECISION :: tolv
3563 DOUBLE PRECISION ::
x(3,*),dens(*),vol(*)
3565 INTEGER ::
k,k1,k2,l,lbest,lchk,lcnt,lnext,lseek,l1,l2,l3,l4, &
3566 l5,nchk,ncont,ncnt,nmon,nprox,n1,n2,n3,n4
3567 DOUBLE PRECISION :: denb,rtest,tols,tolz,vbest,vdiff,vtet, &
3568 vthres,v1,v2,v3,v4,v5,xpt,ypt,zpt
3581 CALL
volcom(lseek,np,vdiff,ncont,
x,ndc)
3583 vthres =
max(tolv,tolz*vol(lseek))
3584 IF (vdiff.LT.vthres) go to 50
3585 IF (ncont.EQ.0) go to 50
3586 IF (ncont.EQ.-1) go to 360
3606 IF (nflag(lseek).EQ.1) go to 20
3607 IF (ndc(4,lseek).EQ.-1.AND.nprox.EQ.0) go to 20
3608 IF (ndc(4,lseek).EQ.-1.AND.nprox.GT.0) go to 15
3609 IF (nprox.EQ.0.AND.n1.NE.nclose.AND.n2.NE.nclose &
3610 .AND.n3.NE.nclose.AND.n4.NE.nclose) go to 20
3613 CALL
volcom(lseek,np,vdiff,ncont,
x,ndc)
3615 vthres =
max(tolv,tolz*vol(lseek))
3616 IF (vdiff.LT.vthres) go to 50
3617 IF (ncont.EQ.0) go to 50
3618 IF (ncont.EQ.-1) go to 360
3619 IF (vdiff.GT.vbest) go to 15
3623 IF (ncnt.GT.mxtest) go to 310
3624 newcel(ncnt) = lseek
3626 IF (nmon.GT.mxnode) go to 315
3627 IF (nmon.GT.mxnode) go to 320
3631 IF (ncnt.GT.0.AND.nprox.EQ.0) go to 25
3632 IF (nprox.EQ.0) go to 45
3633 IF (nprox.EQ.15) go to 320
3635 IF (ncnt.GT.0) go to 25
3637 rtest = vbest/vol(lbest)
3638 IF (rtest.GT.tolz)
WRITE (6,920) nprox,vbest,vol(lbest),rtest
3639 920
FORMAT(
'DIFFICULTY, NPROX = ',i4,
' MIN VDIFF = ',e13.5, &
3640 ' VOL ',e13.5,
' RTEST ',e13.5)
3642 IF (rtest.GT.tolz) go to 320
3647 nfill(
k) = newcel(
k)
3663 IF (iprot(lcont).EQ.1) go to 330
3669 CALL
tetcof(l1,l2,l3,l4,xpt,ypt,zpt,v1,v2,v3,v4,vtet,
x)
3671 dens(np) = (v1*dens(l1) +v2*dens(l2) &
3672 +v3*dens(l3) +v4*dens(l4))/vtet
3673 IF (v1.GE.tolv.AND.v2.GE.tolv.AND.v3.GE.tolv &
3674 .AND.v4.GE.tolv) go to 100
3675 IF (v2.LT.tols.AND.v3.LT.tols.AND.v4.LT.tols) go to 340
3676 IF (v1.LT.tols.AND.v3.LT.tols.AND.v4.LT.tols) go to 340
3677 IF (v1.LT.tols.AND.v2.LT.tols.AND.v4.LT.tols) go to 340
3678 IF (v1.LT.tols.AND.v2.LT.tols.AND.v3.LT.tols) go to 340
3679 IF (v1.LT.tolv) go to 80
3680 IF (v2.LT.tolv) go to 85
3681 IF (v3.LT.tolv) go to 90
3683 CALL
neighb(lcont,l1,l2,l3,lnext,k1,k2,ndc,nbh)
3685 l5 = ndc(1,lnext) +ndc(2,lnext) +ndc(3,lnext) &
3686 +ndc(4,lnext) -l1 -l2 -l3
3688 CALL
tetcof(l1,l2,l3,l5,xpt,ypt,zpt,v1,v2,v3,v5,vtet,
x)
3690 denb = (v1*dens(l1) +v2*dens(l2) +v3*dens(l3) &
3692 dens(np) = .5*(dens(np) +denb)
3695 80 CALL
neighb(lcont,l2,l3,l4,lnext,k1,k2,ndc,nbh)
3697 l5 = ndc(1,lnext) +ndc(2,lnext) +ndc(3,lnext) &
3698 +ndc(4,lnext) -l2 -l3 -l4
3700 CALL
tetcof(l2,l3,l4,l5,xpt,ypt,zpt,v2,v3,v4,v5,vtet,
x)
3702 denb = (v2*dens(l2) +v3*dens(l3) +v4*dens(l4) &
3704 dens(np) = .5*(dens(np) +denb)
3707 85 CALL
neighb(lcont,l3,l4,l1,lnext,k1,k2,ndc,nbh)
3709 l5 = ndc(1,lnext) +ndc(2,lnext) +ndc(3,lnext) &
3710 +ndc(4,lnext) -l3 -l4 -l1
3712 CALL
tetcof(l3,l4,l1,l5,xpt,ypt,zpt,v3,v4,v1,v5,vtet,
x)
3714 denb = (v3*dens(l3) +v4*dens(l4) +v1*dens(l1) &
3716 dens(np) = .5*(dens(np) +denb)
3719 90 CALL
neighb(lcont,l4,l1,l2,lnext,k1,k2,ndc,nbh)
3721 l5 = ndc(1,lnext) +ndc(2,lnext) +ndc(3,lnext) &
3722 +ndc(4,lnext) -l4 -l1 -l2
3724 CALL
tetcof(l4,l1,l2,l5,xpt,ypt,zpt,v4,v1,v2,v5,vtet,
x)
3726 denb = (v4*dens(l4) +v1*dens(l1) +v2*dens(l2) &
3728 dens(np) = .5*(dens(np) +denb)
3747 WRITE (6,640) np,
x(1,np),
x(2,np),
x(3,np), &
3748 l1,l2,l3,l4,vtet,v1,v2,v3,v4
3756 WRITE (6,910) np,
x(1,np),
x(2,np),
x(3,np),ncont
3757 910
FORMAT(
'NP = ',i6,
' X = ',f8.3,
' Y = ',f8.3,
' Z = ',f8.3, &
3761 610
FORMAT(//5
x,
'DIMENSION OF NEWCEL EXCEEDED IN ROUTINE TETLOC'/ &
3762 5
x,
'INCREASE SIZE OF MTEST')
3763 615
FORMAT(//5
x,
'DIMENSION OF IFLAG EXCEEDED IN ROUTINE TETLOC'/ &
3764 5
x,
'INCREASE SIZE OF MBPTS')
3765 620
FORMAT(5
x,
'UNABLE TO FIND TETRAHEDRON THAT CONTAINS NEW POINT', &
3766 ' IN ROUTINE TETLOC')
3767 630
FORMAT(5
x,
'POINT LIES OUTSIDE DOMAIN TO BE MESHED')
3768 640
FORMAT(//5
x,
'NEW POINT, NP = ',i6,
' X = ',f6.2,
' Y = ',f6.2, &
3769 ' Z = ',f6.2/5
x,
'APPEARS TO BE COINCIDENT WITH AN', &
3770 ' EXISTING POINT IN THE MESH'/ &
3771 5
x,
'VERTEX ADDRESSES OF CONTAINING TETRAHEDRON'/ &
3772 5
x,
'L1 ',i6,
' L2 ',i6,
' L3 ',i6,
' L4 ',i6,
' VOL = ',e13.5/ &
3773 5
x,
'V1 = ',e13.5,
' V2 = ',e13.5,
' V3 = ',e13.5, &
3774 ' V4 = ',e13.5/5
x,
'PROGRAM STOPPED IN ROUTINE TETLOC')
3775 660
FORMAT(//5
x,
'A GEOMETRIC INCONSISTENCY HAS BEEN FOUND IN ROUTINE', &
3776 ' TETLOC.'/5
x,
'THE INSERTED POINT HAS A NEGATIVE', &
3777 ' ORIENTATION WITH RESPECT TO EVERY TETRAHEDRON FACE')
3789 x,ndc,nbh,iprot,iflag,nflag, &
3790 xcen,ycen,zcen,rc,nfill,newcel,tolv)
3807 INTEGER :: lcont,ndel,np
3808 INTEGER :: iflag(*),nflag(*)
3809 INTEGER :: nfill(*),newcel(*)
3811 INTEGER :: ndc(4,*),nbh(4,*),iprot(*)
3812 DOUBLE PRECISION :: tolv
3813 DOUBLE PRECISION ::
x(3,*)
3814 DOUBLE PRECISION :: xcen(*),ycen(*),zcen(*),rc(*)
3816 INTEGER :: l,lseek,l1,l2,l3,
k,k1,k2,lchk,lcnt,lnext,ncnt,nmon,
m, &
3818 DOUBLE PRECISION :: a1,a2,a3,b1,b2,b3,dcen,rnx,rny,rnz, &
3819 thres,tolp,vtet1,vtet2,xpt,ypt,zpt
3846 IF (nflag(lseek).EQ.1) go to 20
3851 IF (ndc(4,lseek).EQ.-1) go to 20
3853 IF (nmon.GT.mxnode) go to 315
3860 dcen =
sqrt((xpt -xcen(lseek))**2 +(ypt -ycen(lseek))**2 &
3861 +(zpt -zcen(lseek))**2)
3862 IF (dcen.GE.rc(lseek)*(1.+tolp))go to 20
3866 IF (iprot(lseek).EQ.1) go to 20
3873 lnext = nbh(
m,lseek)
3874 IF (iprot(lnext).EQ.0) go to 15
3876 CALL
lock(lseek,lnext,n1,n2,n3,n4,n5,k1,k2,ndc,nbh)
3880 l2 = n1 +n2 +n3 -l1 -l3
3881 rnx =
cofact(
x(2,l1),
x(2,l2),
x(2,l3),
x(3,l1),
x(3,l2),
x(3,l3))
3882 rny =
cofact(
x(3,l1),
x(3,l2),
x(3,l3),
x(1,l1),
x(1,l2),
x(1,l3))
3883 rnz =
cofact(
x(1,l1),
x(1,l2),
x(1,l3),
x(2,l1),
x(2,l2),
x(2,l3))
3884 vtet1 = rnx*(xpt -
x(1,l1)) +rny*(ypt -
x(2,l1)) &
3886 vtet2 = rnx*(
x(1,n4) -
x(1,l1)) +rny*(
x(2,n4) -
x(2,l1)) &
3887 +rnz*(
x(3,n4) -
x(3,l1))
3888 thres =
max(1.0d0,abs(vtet2))
3889 IF (abs(vtet1).LT.tolv*thres) go to 20
3890 IF (vtet1*vtet2.LT.0.) go to 20
3896 IF (ncnt.GT.mxtest) go to 310
3897 newcel(ncnt) = lseek
3899 IF (ndel.GT.mxtest) go to 320
3902 IF (ncnt.EQ.0) go to 40
3905 nfill(
k) = newcel(
k)
3924 610
FORMAT(//5
x,
'DIMENSION OF NEWCEL EXCEEDED IN ROUTINE CAVITY'/ &
3925 5
x,
'INCREASE SIZE OF MTEST')
3926 615
FORMAT(//5
x,
'DIMENSION OF IFLAG EXCEEDED IN ROUTINE CAVITY'/ &
3927 5
x,
'INCREASE SIZE OF MNODE')
3928 620
FORMAT(//5
x,
'DIMENSION OF LDEL EXCEEDED IN ROUTINE CAVITY'/ &
3929 5
x,
'INCREASE SIZE OF MTEST')
3938 SUBROUTINE cavbnd (NP,LCONT,NDEL,LDEL,NCNT, &
3939 x,ndc,nbh,vol,iflag,nflag, &
3940 nfill,newcel,ntri,newnbh,nold,tolv)
3956 INTEGER :: lcont,ncnt,ndel,np
3957 INTEGER :: ndc(4,*),nbh(4,*)
3958 INTEGER :: iflag(*),nflag(*)
3959 INTEGER :: ntri(3,*),nfill(*),newnbh(4,*),nold(*),newcel(*)
3961 DOUBLE PRECISION :: tolv
3962 DOUBLE PRECISION ::
x(3,*),vol(*)
3964 INTEGER ::
j,
k,l,lchk,lcnt,ll,ltry,
m,mend,mmax,mmin,msum, &
3965 m1,m2,m3,m4,
n,ncont,nend,njoin,nmax,nmon,nmin, &
3966 nsum,nwatch,n1,n2,n3,n4
3967 DOUBLE PRECISION :: a1,a2,a3,b1,b2,b3,fac,rnx,rny,rnz,vdiff, &
3968 vtet,vtet1,vtet2,v1,v2,v3,v4,xpt,ypt,zpt
3985 IF (nflag(l).GT.0) go to 65
3992 nmin =
min(n1,n2,n3)
3993 nmax =
max(n1,n2,n3)
3995 30 mmin =
min(m1,m2,m3)
3996 mmax =
max(m1,m2,m3)
3998 IF (mmin.EQ.nmin.AND.mmax.EQ.nmax.AND.msum.EQ.nsum) go to 40
3999 IF (m1.EQ.mend) go to 35
4006 35
IF (n1.EQ.nend) go to 300
4017 40 rnx =
cofact(
x(2,n1),
x(2,n2),
x(2,n3),
x(3,n1),
x(3,n2),
x(3,n3))
4018 rny =
cofact(
x(3,n1),
x(3,n2),
x(3,n3),
x(1,n1),
x(1,n2),
x(1,n3))
4019 rnz =
cofact(
x(1,n1),
x(1,n2),
x(1,n3),
x(2,n1),
x(2,n2),
x(2,n3))
4020 fac = 1./
sqrt(rnx*rnx +rny*rny +rnz*rnz)
4021 vtet1 = fac*(rnx*(xpt -
x(1,n1)) +rny*(ypt -
x(2,n1)) &
4022 +rnz*(zpt -
x(3,n1)))
4023 vtet2 = rnx*(
x(1,n4) -
x(1,n1)) +rny*(
x(2,n4) -
x(2,n1)) &
4024 +rnz*(
x(3,n4) -
x(3,n1))
4027 IF (ldel(
k).NE.lcont.AND.abs(vtet1).LT.tolv) go to 45
4028 IF (vtet1*vtet2.GE.0.) go to 50
4029 45 nflag(ldel(
k)) = 2
4037 IF (ncnt.GT.mxtest) go to 310
4043 nold(ncnt) = ldel(
k)
4045 IF (njoin.EQ.4.AND.ndel.GT.1) go to 320
4047 IF (nwatch.EQ.0)
RETURN
4048 IF (nflag(lcont).EQ.2) go to 330
4063 IF (nflag(ltry).NE.1) go to 85
4068 IF (nmon.GT.mxnode) go to 335
4071 IF (ncnt.EQ.0) go to 95
4074 nfill(
k) = newcel(
k)
4101 CALL
volcom(lcont,np,vdiff,ncont,
x,ndc)
4108 CALL
tetcof(n1,n2,n3,n4,xpt,ypt,zpt,v1,v2,v3,v4,vtet,
x)
4110 WRITE (6,998) np,ndel,lcont,n1,n2,n3,n4, &
4111 v1,v2,v3,v4,vtet,vdiff,ncont
4112 998
FORMAT(
'NP ',i6,
' NDEL ',i4,
' LDEL ',i7,
' N1,N2,N3,N4 ',4i6/ &
4113 'V1 ',e13.5,
' V2 ',e13.5,
' V3 ',e13.5,
' V4 ',e13.5/ &
4114 ' VTET ',e13.5,
' VDIFF ',e13.5,
' NCONT ',i2)
4121 996
WRITE (6,997) ll,nflag(ll),n1,n2,n3,n4,vol(ll)
4122 997
FORMAT(
'CELL ',i7,
' NFLAG ',i2,
' N1,N2,N3,N4 ',4i6,
' VOL ',e13.5)
4127 600
FORMAT(//5
x,
'UNABLE TO FIND A CONTIGUITY BETWEEN A CAVITY'/ &
4128 5
x,
'TETRAHEDRON AND A NON-CAVITY TETRAHEDRON'/ &
4129 5
x,
'PROGRAM STOPPED IN ROUTINE CAVBND')
4130 610
FORMAT(//5
x,
'DIMENSION OF NTRI ARRAY EXCEEDED IN ROUTINE CAVBND'/ &
4131 5
x,
'INCREASE SIZE OF MTEST')
4132 620
FORMAT(//5
x,
'FOUR CONTIGUITIES HAVE BEEN FOUND FOR A NEW' &
4134 5
x,
'NEW TETRAHEDRON IS THEREFORE ISOLATED'/ &
4135 5
x,
'PROGRAM STOPPED IN ROUTINE CAVBND')
4136 630
FORMAT(//5
x,
'TETRAHEDRON CONTAINING POINT HAS FAILED'/ &
4137 5
x,
'VISIBILITY TEST'/ &
4138 5
x,
'PROGRAM STOPPED IN ROUTINE CAVBND')
4139 635
FORMAT(//5
x,
'DIMENSION OF IFLAG EXCEEDED IN ROUTINE CAVBND.'/ &
4140 5
x,
'INCREASE SIZE OF MBPTS')
4152 ndc,nbh,iprot,ndg,idgp,ndgp,nflag, &
4153 ntri,ncavfc,ikeep,nedgrm,iedgrm)
4168 INTEGER :: iedgrm,ncnt,ndel
4169 INTEGER :: ndc(4,*),nbh(4,*),iprot(*),ndg(2,*),nflag(*)
4170 INTEGER :: ntri(3,*),nedgrm(*)
4171 INTEGER :: idgp(*),ndgp(*)
4172 INTEGER :: ldel(*),ncavfc(3,*),ikeep(*)
4174 INTEGER ::
j,
k,kcnt,l,l1,l2,
m,mend,mmax,mmin,msum,m1,m2,m3,m4, &
4175 n,nedg,nend,nf,nmax,nmin,nsum,n1,n2,n3,n4
4184 IF (iprot(ldel(
k)).EQ.1) go to 370
4188 IF (nflag(l).EQ.0) go to 65
4199 nmin =
min(n1,n2,n3)
4200 nmax =
max(n1,n2,n3)
4202 45 mmin =
min(m1,m2,m3)
4203 mmax =
max(m1,m2,m3)
4205 IF (mmin.EQ.nmin.AND.mmax.EQ.nmax.AND.msum.EQ.nsum) go to 60
4206 IF (m1.EQ.mend) go to 50
4213 50
IF (n1.EQ.nend) go to 300
4221 IF (nf.GT.mxtest) go to 310
4239 95
IF (nedg.EQ.0) go to 320
4240 IF (l2.EQ.ndg(2,nedg)) go to 97
4243 97
IF (nflag(nedg).EQ.1) go to 98
4246 IF (kcnt.GT.mxtest) go to 330
4248 98
IF (n1.EQ.ntri(3,l)) go to 100
4265 120
IF (nedg.EQ.0) go to 320
4266 IF (l2.EQ.ndg(2,nedg)) go to 130
4269 130
IF (nflag(nedg).NE.0) go to 135
4272 IF (kcnt.GT.mxtest) go to 330
4274 135
IF (n1.EQ.ncavfc(3,
j)) go to 140
4281 IF (kcnt.EQ.0)
RETURN
4284 IF (nflag(nedg).NE.2) go to 170
4286 CALL
edgerm(nedg,ndg,idgp,ndgp)
4288 IF (iedgrm.EQ.mxtest) go to 170
4290 nedgrm(iedgrm) = nedg
4309 600
FORMAT(//5
x,
'UNABLE TO FIND A CONTIGUITY BETWEEN CAVITY'/ &
4311 5
x,
'PROGRAM STOPPED IN ROUTINE RECON')
4312 610
FORMAT(//5
x,
'DIMENSION OF NTRI ARRAY EXCEEDED IN ROUTINE RECON'/ &
4313 5
x,
'INCREASE SIZE OF MTEST')
4314 620
FORMAT(5
x,
'UNABLE TO FIND EDGE IN CAVITY AMONG NDG ARRAY',/ &
4315 'PROGRAM STOPPED IN ROUTINE RECON')
4316 630
FORMAT(//5
x,
'DIMENSION OF IKEEP ARRAY EXCEEDED IN ROUTINE RECON'/ &
4317 5
x,
'INCREASE SIZE OF MTEST')
4318 650
FORMAT(5
x,
'AT LEAST ONE INTERNAL CAVITY FACE IS A PROTECTED'/ &
4319 5
x,
'BOUNDARY FACE. PROGRAM STOPPED IN ROUTINE RECON')
4320 670
FORMAT(/5
x,
'AT LEAST ONE CAVITY TETRAHEDRON IS PROTECTED.'/ &
4321 5
x,
'PROGRAM STOPPED IN ROUTINE RECON.')
4322 680
FORMAT(//5
x,
'UNABLE TO FIND FACE IN LINKED LIST'/ &
4323 5
x,
'PROGRAM STOPPED IN RECON')
4324 END SUBROUTINE recon
4334 SUBROUTINE cavedg (NCNT,NEDG,IPOINT,NPOINT, &
4335 ntri,ncav,newnbh,newcel)
4351 INTEGER :: ncnt,nedg
4352 INTEGER :: npoint(*),ipoint(*)
4353 INTEGER :: ntri(3,*),newnbh(4,*),newcel(*),ncav(4,*)
4355 INTEGER ::
i,j1,j2,
k,l,l1,l2,l3,
next,n1,n2
4367 IF (
i.EQ.0) go to 30
4368 20
IF (n2.EQ.ncav(2,
i)) go to 40
4370 IF (
next.EQ.0) go to 30
4374 IF (nedg.GT.mxcav) go to 300
4380 IF (
i.NE.0) npoint(
i) = nedg
4381 IF (
i.EQ.0) ipoint(n1) = nedg
4384 50
IF (l1.EQ.ntri(3,
k)) go to 60
4400 npoint(j1) = npoint(j1) +1
4401 npoint(j2) = npoint(j2) +1
4402 newnbh(npoint(j1),j1) = newcel(j2)
4403 newnbh(npoint(j2),j2) = newcel(j1)
4406 IF (npoint(
k).NE.4) go to 310
4412 IF (ntri(1,
k).GT.0) ipoint(ntri(1,
k)) = 0
4413 IF (ntri(2,
k).GT.0) ipoint(ntri(2,
k)) = 0
4414 IF (ntri(3,
k).GT.0) ipoint(ntri(3,
k)) = 0
4423 600
FORMAT(//5
x,
'DIMENSION OF NCAV EXCEEDED IN ROUTINE CAVEDG'/ &
4424 5
x,
'INCREASE SIZE OF MCAV')
4425 610
FORMAT(//5
x,
'UNABLE TO FIND ALL CONTIGUITIES BETWEEN NEW', &
4426 ' TETRAHEDRA'/5
x,
'PROGRAM STOPPED IN ROUTINE CAVEDG')
4436 SUBROUTINE datsrf (NP,NCNT,NEDG,NDG,IDGP,NDGP,NEDGE, &
4437 ipoint,ncav,nedgrm,iedgrm)
4452 INTEGER :: iedgrm,ncnt,nedg,nedge,np
4453 INTEGER :: ndg(2,*),idgp(*),ndgp(*)
4454 INTEGER :: ipoint(*)
4455 INTEGER :: ncav(4,*),nedgrm(*)
4457 INTEGER ::
i,iedg,j1,j2,l1,l2,lp,nedcnt,newed,n1,n2
4470 IF (ipoint(n1).GT.0) go to 65
4474 IF (iedg.EQ.0) go to 50
4475 40
IF (ndgp(iedg).EQ.0) go to 50
4478 50 nedcnt = nedcnt +1
4479 IF (nedcnt.GT.iedgrm) go to 55
4480 newed = nedgrm(nedcnt)
4483 IF (nedge.GT.mxedge) go to 310
4485 60 ndg(1,newed) = l1
4488 IF (iedg.EQ.0) idgp(l1) = newed
4489 IF (iedg.GT.0) ndgp(iedg) = newed
4491 65
IF (ipoint(n2).GT.0) go to 80
4495 IF (iedg.EQ.0) go to 68
4496 66
IF (ndgp(iedg).EQ.0) go to 68
4499 68 nedcnt = nedcnt +1
4500 IF (nedcnt.GT.iedgrm) go to 70
4501 newed = nedgrm(nedcnt)
4504 IF (nedge.GT.mxedge) go to 310
4506 75 ndg(1,newed) = l2
4509 IF (iedg.EQ.0) idgp(l2) = newed
4510 IF (iedg.GT.0) ndgp(iedg) = newed
4516 610
FORMAT(//5
x,
'DIMENSION OF NDG ARRAY EXCEEDED IN ROUTINE DATSRF.'/ &
4517 5
x,
'INCREASE SIZE OF MEDGE.')
4529 ityp,xcen,ycen,zcen,vol,rc,rat, &
4530 nvcnt,iflag,nflag,nptet, &
4532 noctr,ioctr,nlink,xfar,yfar,zfar,idone,nref, &
4533 ksrch,nsrch,iring,ntetkp, &
4534 lnbr,ishk,mnbr,kshk,tolv)
4551 INTEGER :: ityp(*),ndc(4,*),nbh(4,*),iprot(*),ndg(2,*)
4552 INTEGER :: idgp(*),ndgp(*)
4553 INTEGER :: idone(*),nref(*),nlink(*),noctr(2,*)
4554 INTEGER :: nvcnt(*),nptet(*)
4555 INTEGER :: iflag(*),nflag(*)
4556 INTEGER :: nsrch(*),ksrch(*)
4557 INTEGER :: iring(*),ntetkp(*),lnbr(*),ishk(*),mnbr(*),kshk(*)
4558 DOUBLE PRECISION :: tolv
4559 DOUBLE PRECISION ::
x(3,*)
4560 DOUBLE PRECISION :: vol(*),xcen(*),ycen(*),zcen(*),rc(*),rat(*)
4561 DOUBLE PRECISION :: xfar(2),yfar(2),zfar(2)
4563 INTEGER ::
i,iedg,imon,inad,ip1,ismall,itypt,jnext,jpre,
k,kcnt, &
4564 kk,k1,k2,l,lc,ln,l1,l2,
m,m1,m2,m3,m4,
n,nabdy,nbbdy, &
4565 ncnt,nexch,nkeep,np,npass,npnext,npre,nrbdy,nrngp1, &
4567 DOUBLE PRECISION :: area,a1,a2,a3,b1,b2,b3,radc,radmax,ratmax, &
4568 rcrin,rnx,rny,rnz,tols,vl,vtet1,vtet2,xcn, &
4576 IF (na.LT.0.OR.nb.LT.0)
RETURN
4577 IF (ityp(na).LT.0.OR.ityp(nb).LT.0)
RETURN
4581 ksrch(1) = nptet(na)
4587 IF (lc.LE.0) go to 300
4588 IF (nflag(lc).EQ.1) go to 22
4593 IF (m1.NE.nb.AND.m2.NE.nb.AND.m3.NE.nb.AND.m4.NE.nb) go to 15
4597 IF (imon.GT.mxnode) go to 310
4602 IF (nflag(ln).EQ.1) go to 20
4607 IF (m1.NE.na.AND.m2.NE.na.AND.m3.NE.na.AND.m4.NE.na) go to 20
4609 IF (ncnt.GT.mxtest) go to 320
4613 IF (ncnt.EQ.0) go to 200
4622 30
IF (imon.EQ.0) go to 40
4631 IF (
m.NE.na.AND.
m.NE.nb) go to 50
4634 n2 = ndc(1,lc) +ndc(2,lc) +ndc(3,lc) +ndc(4,lc) &
4637 IF (nvcnt(n1).EQ.4)
THEN
4649 CALL
neighb(lc,n1,na,nb,jnext,k1,k2,ndc,nbh)
4651 npnext = ndc(1,jnext) +ndc(2,jnext) +ndc(3,jnext) &
4652 +ndc(4,jnext) -n1 -na -nb
4657 IF (iprot(jnext).EQ.1) nrbdy = 1
4661 CALL
neighb(jpre,npre,na,nb,jnext,k1,k2,ndc,nbh)
4663 npnext = ndc(1,jnext) +ndc(2,jnext) +ndc(3,jnext) &
4664 +ndc(4,jnext) -npre -na -nb
4666 IF (nvcnt(npre).EQ.4)
THEN
4675 IF (nring.GT.mxring) go to 340
4677 ntetkp(nring) = jnext
4679 iflag(nring) = jnext
4680 IF (iprot(jnext).EQ.1) nrbdy = 1
4681 IF (npnext.EQ.n1) go to 350
4682 IF (npnext.NE.n2) go to 60
4684 IF (nvcnt(n2).EQ.4)
THEN
4693 IF (nring.GT.mxring) go to 340
4705 65 ksrch(1) = nptet(np)
4712 IF (lc.LE.0) go to 300
4713 IF (nflag(lc).GT.0) go to 70
4714 IF (n1.NE.np.AND.n2.NE.np.AND.n3.NE.np.AND.n4.NE.np) go to 70
4717 IF (imon.GT.mxnode) go to 310
4719 IF (iprot(lc).EQ.1.AND.np.EQ.na) nabdy = 1
4720 IF (iprot(lc).EQ.1.AND.np.EQ.nb) nbbdy = 1
4726 IF (lc.LE.0) go to 300
4727 IF (nflag(lc).GT.0) go to 75
4732 IF (n1.NE.np.AND.n2.NE.np.AND.n3.NE.np.AND.n4.NE.np) go to 75
4734 IF (ncnt.GT.mxtest) go to 320
4738 IF (imon.GT.mxnode) go to 310
4740 IF (iprot(lc).EQ.1.AND.np.EQ.na) nabdy = 1
4741 IF (iprot(lc).EQ.1.AND.np.EQ.nb) nbbdy = 1
4743 IF (ncnt.EQ.0) go to 85
4749 85
IF (np.EQ.nb) go to 90
4758 radmax =
max(radmax,rc(ntetkp(
i)))
4759 ratmax =
max(ratmax,rat(ntetkp(
i)))
4760 ip1 = mod(
i,nring) +1
4762 CALL
neighb(ntetkp(
i),iring(
i),iring(ip1),na,lnbr(
i),kk,ishk(
i), &
4765 CALL
neighb(ntetkp(
i),iring(
i),iring(ip1),nb,mnbr(
i),kk,kshk(
i), &
4774 IF (nabdy.EQ.1.AND.nbbdy.EQ.1) go to 230
4777 xpt = .5*(
x(1,na) +
x(1,nb))
4778 ypt = .5*(
x(2,na) +
x(2,nb))
4779 zpt = .5*(
x(3,na) +
x(3,nb))
4780 itypt =
max(ityp(na),ityp(nb))
4781 IF (nrbdy.EQ.nabdy.AND.nrbdy.EQ.nbbdy) go to 118
4786 IF (nrbdy.EQ.0.AND.nabdy.EQ.1.AND.nbbdy.EQ.0) go to 118
4794 IF (nrbdy.EQ.0.AND.nabdy.EQ.0.AND.nbbdy.EQ.1)
THEN
4804 118 nrngp1 = nring +1
4805 DO 125
i=nrngp1,imon
4807 IF (iprot(l).EQ.1) go to 125
4812 120
IF (n4.EQ.na.OR.n4.EQ.nb) go to 123
4819 123 rnx =
cofact(
x(2,n1),
x(2,n2),
x(2,n3),
x(3,n1),
x(3,n2),
x(3,n3))
4820 rny =
cofact(
x(3,n1),
x(3,n2),
x(3,n3),
x(1,n1),
x(1,n2),
x(1,n3))
4821 rnz =
cofact(
x(1,n1),
x(1,n2),
x(1,n3),
x(2,n1),
x(2,n2),
x(2,n3))
4822 vtet1 = rnx*(
x(1,n4) -
x(1,n1)) +rny*(
x(2,n4) -
x(2,n1)) &
4823 +rnz*(
x(3,n4) -
x(3,n1))
4824 vtet2 = rnx*(xpt -
x(1,n1)) +rny*(ypt -
x(2,n1)) &
4826 IF (abs(vtet2).LT.tols.OR.vtet1*vtet2.LT.0.0)
THEN
4827 IF (npass.EQ.1)
THEN
4831 IF (n4.EQ.na.AND.nabdy.EQ.0.AND.nbbdy.EQ.1) go to 220
4832 IF (n4.EQ.nb.AND.nabdy.EQ.1.AND.nbbdy.EQ.0) go to 220
4841 CALL
circum2(
x,n1,n2,n3,xpt,ypt,zpt, &
4842 xcn,ycn,zcn,vl,radc,ismall,tolv)
4844 IF (radc.GT.50.*radmax) go to 220
4845 area =
tetar3(n1,n2,n3,xpt,ypt,zpt,
x)
4846 rcrin = radc*area/vl
4847 IF (rcrin.GT.50.*ratmax) go to 220
4857 CALL
octrmv(na,
x,noctr,nlink,idone,nref)
4859 CALL
octrmv(nb,
x,noctr,nlink,idone,nref)
4867 CALL
octfil(na,
x,noctr,ioctr,nlink,nref,xfar,yfar,zfar)
4875 nbh(1,ntetkp(
i)) = 0
4876 nbh(2,ntetkp(
i)) = 0
4877 nbh(3,ntetkp(
i)) = 0
4878 nbh(4,ntetkp(
i)) = 0
4879 nbh(ishk(
i),lnbr(
i)) = mnbr(
i)
4880 nbh(kshk(
i),mnbr(
i)) = lnbr(
i)
4881 IF (iring(
i).GT.0) nvcnt(iring(
i)) = nvcnt(iring(
i)) -1
4883 nvcnt(na) = nvcnt(na) +nvcnt(nb) -nring -2
4891 135
IF (iedg.EQ.0) go to 330
4892 IF (l2.EQ.ndg(2,iedg)) go to 140
4896 140 CALL
edgerm(iedg,ndg,idgp,ndgp)
4899 IF (iring(
i).LT.0) go to 152
4900 l1 =
min(nb,iring(
i))
4901 l2 =
max(nb,iring(
i))
4903 145
IF (iedg.EQ.0) go to 330
4904 IF (l2.EQ.ndg(2,iedg)) go to 150
4908 150 CALL
edgerm(iedg,ndg,idgp,ndgp)
4914 DO 162
i=nrngp1,imon
4916 IF (nflag(l).EQ.na) go to 162
4921 IF (n1.EQ.nb) ndc(1,l) = na
4922 IF (n2.EQ.nb) ndc(2,l) = na
4923 IF (n3.EQ.nb) ndc(3,l) = na
4924 IF (n4.EQ.nb) ndc(4,l) = na
4929 DO 177
i=nrngp1,imon
4931 IF (nflag(l).EQ.na) go to 177
4934 IF (
n.LT.0.OR.
n.EQ.na) go to 175
4938 165
IF (iedg.EQ.0) go to 175
4939 IF (l2.EQ.ndg(2,iedg)) go to 170
4943 170 CALL
edgerm(iedg,ndg,idgp,ndgp)
4948 IF (inad.EQ.0) go to 174
4949 172
IF (ndgp(inad).EQ.0) go to 174
4952 174 ndg(1,iedg) = m1
4955 IF (inad.EQ.0) idgp(m1) = iedg
4956 IF (inad.GT.0) ndgp(inad) = iedg
4962 DO 195
k=nrngp1,imon
4964 IF (iprot(l).EQ.1) go to 195
4970 CALL
circum(
x,n1,n2,n3,n4,xcn,ycn,zcn,vl,radc,ismall,tolv)
4972 IF (ismall.EQ.1) go to 360
4979 rat(l) = rc(l)*area/vol(l)
4991 700
FORMAT(
'EDGE WITH VERTICES ',i6,
' , ',i6,
' DOES NOT EXIST.'/ &
4992 'COLLAPSE OF THESE TWO POINTS IS NOT POSSIBLE.')
4996 720
FORMAT(
'ATTEMPTED COLLAPSE OF EDGE WITH VERTICES ',i6,
' , ',i6/ &
4997 'CREATES A NON-CONVEX TETRAHEDRON')
5001 725
FORMAT(
'ATTEMPTED COLLAPSE OF AN EDGE THAT DOES NOT LIE IN THE'/ &
5002 'BOUNDARY SURFACE BUT WHOSE END-POINTS ',i6,
' AND ',i6/ &
5003 'BOTH LIE IN THE BOUNDARY SURFACE')
5011 300
WRITE (6,600) na,nb,lc
5025 600
FORMAT(///5
x,
'INVALID TETRAHEDRON ADDRESS FOUND WHILE SEARCHING'/ &
5026 5
x,
'FOR A TETRAHEDRON INCIDENT TO EDGE WITH VERTICES', &
5027 i6,
' , ',i6/5
x,
'TETRAHEDRON ADDRESS IS ',i7/ &
5028 5
x,
'PROGRAM STOPPED IN ROUTINE COLAPS')
5029 610
FORMAT(///5
x,
'DIMENSION OF ARRAY IFLAG EXCEEDED. INCREASE SIZE'/ &
5030 5
x,
'OF MNODE. PROGRAM STOPPED IN COLAPS')
5031 620
FORMAT(///5
x,
'DIMENSION OF ARRAY NSRCH EXCEEDED. INCREASE SIZE'/ &
5032 5
x,
'OF MTEST. PROGRAM STOPPED IN COLAPS')
5033 650
FORMAT(///5
x,
'SEARCH FOR TETRAHEDRAL RING HAS RETURNED TO THE'/ &
5034 5
x,
'STARTING FACE. THIS INDICATES AN INCONSISTENCY IN'/ &
5035 5
x,
'THE TETRAHEDRAL ENSEMBLE.'/ &
5036 5
x,
'PROGRAM STOPPED IN COLAPS')
5037 630
FORMAT(///5
x,
'UNABLE TO FIND EDGE ADDRESS FOR A NEW TETRAHEDRON'/ &
5038 5
x,
'PROGRAM STOPPED IN ROUTINE COLAPS')
5039 640
FORMAT(///5
x,
'DIMENSION OF ARRAY IRING EXCEEDED. INCREASE SIZE'/ &
5040 5
x,
'OF MRING. PROGRAM STOPPED IN COLAPS')
5041 660
FORMAT(///5
x,
'AT LEAST ONE NEW TETRAHEDRON HAS TOO SMALL A VOLUME' &
5042 /5
x,
'PROGRAM STOPPED IN COLAPS')
5053 SUBROUTINE smooth (X,ITYP,NNODE,NDC,NBH,IPROT,NCELL, &
5054 ndg,idgp,ndgp,nedge, &
5055 vol,xcen,ycen,zcen,rc,rat,dens,nptet,nacpt, &
5056 idone,nref,nlink,noctr,ioctr,xfar,yfar,zfar, &
5057 xoctr,yoctr,zoctr,xhold,yhold,zhold, &
5058 xkeep,ykeep,zkeep,ksrch,nsrch, &
5059 ipoint,npoint,iflag,nflag, &
5060 dx,
dy,
dz,ds,vlt,iring,ntetkp,nfad,newc, &
5061 nbhkp,iedkp,lnbr,ishk,mnbr,kshk,npp, &
5062 nfill,newcel,ntri, &
5063 ncav,nshake,newnbh,nold,ncavfc,ikeep,ldel, &
5064 nedgrm,xc,yc,zc,
v,rad,rcrin,lnkup,lnkdn, &
5065 listf,volmin,rcmx,tolv)
5081 INTEGER :: ioctr,ncell,nedge,nnode
5082 INTEGER :: ndc(4,*),nbh(4,*),iprot(*),ndg(2,*),idgp(*),ndgp(*)
5083 INTEGER :: idone(*),nref(*),nlink(*),noctr(2,*)
5084 INTEGER :: ityp(*),nptet(*),nacpt(*),listf(*)
5085 INTEGER :: npoint(*),ipoint(*)
5086 INTEGER :: iflag(*),nflag(*)
5087 INTEGER :: lnkup(*),lnkdn(*)
5088 INTEGER :: nsrch(*),ksrch(*)
5089 INTEGER :: iring(*),ntetkp(*),nfad(3,*),newc(*), &
5090 nbhkp(3,*),iedkp(4,*),lnbr(*),ishk(*), &
5091 mnbr(*),kshk(*),npp(*)
5092 INTEGER :: ntri(3,*),nfill(*),newnbh(4,*),nold(*), &
5093 newcel(*),nshake(*),ncav(4,*)
5094 INTEGER :: nedgrm(*)
5095 INTEGER :: ldel(*),ncavfc(3,*),ikeep(*)
5096 DOUBLE PRECISION :: volmin,rcmx,tolv
5097 DOUBLE PRECISION ::
x(3,*),dens(*)
5098 DOUBLE PRECISION :: vol(*),xcen(*),ycen(*),zcen(*),rc(*),rat(*)
5099 DOUBLE PRECISION :: xoctr(2,*),yoctr(2,*),zoctr(2,*), &
5100 xhold(2,*),yhold(2,*),zhold(2,*), &
5101 xfar(2),yfar(2),zfar(2),xkeep(2), &
5103 DOUBLE PRECISION ::
dx(*),
dy(*),
dz(*),ds(*),vlt(*)
5104 DOUBLE PRECISION :: xc(*),yc(*),zc(*),
v(*),rad(*),rcrin(*)
5105 DOUBLE PRECISION :: ang(6)
5107 INTEGER ::
i,i1,i2,i3,
j,
k,l,loop,m1,m2,
n,na,nb,nchk,nfail, &
5108 nmax,nmin,nn,np,npass,npt,nring,nsum,ntet, &
5110 DOUBLE PRECISION :: angmax,angmx1,angmx2,angmx3
5122 CALL
tree(rat,lnkdn,lnkup,nflag,listf,nacpt,nwait, &
5125 30
WRITE (6,910) nwait,loop
5126 910
FORMAT(//
'IN SMOOTH, NWAIT = ',i7,
' ITERATION ',i2)
5134 IF (nbh(1,
j).LE.0) go to 130
5135 IF (nacpt(
j).EQ.1) go to 130
5141 ang(1) =
dihed(n1,n2,n3,n4,
x)
5142 IF (ang(1).GT.180.) ang(1) = 360. -ang(1)
5143 ang(2) =
dihed(n2,n3,n4,n1,
x)
5144 IF (ang(2).GT.180.) ang(2) = 360. -ang(2)
5145 ang(3) =
dihed(n3,n4,n1,n2,
x)
5146 IF (ang(3).GT.180.) ang(3) = 360. -ang(3)
5147 ang(4) =
dihed(n4,n1,n2,n3,
x)
5148 IF (ang(4).GT.180.) ang(4) = 360. -ang(4)
5149 ang(5) =
dihed(n1,n3,n2,n4,
x)
5150 IF (ang(5).GT.180.) ang(5) = 360. -ang(5)
5151 ang(6) =
dihed(n2,n4,n1,n3,
x)
5152 IF (ang(6).GT.180.) ang(6) = 360. -ang(6)
5159 IF (ang(
k).LT.angmx1) go to 50
5165 IF (
k.EQ.i1) go to 60
5166 IF (i2.EQ.0) go to 55
5167 IF (ang(
k).LT.angmx2) go to 60
5173 IF (
k.EQ.i1.OR.
k.EQ.i2) go to 70
5174 IF (i3.EQ.0) go to 65
5175 IF (ang(
k).LT.angmx3) go to 70
5183 nmin =
min(i1,i2,i3)
5184 nmax =
max(i1,i2,i3)
5186 IF (nmin.EQ.1.AND.nmax.EQ.5.AND.nsum.EQ.10) np = n1
5187 IF (nmin.EQ.1.AND.nmax.EQ.6.AND.nsum.EQ.9) np = n2
5188 IF (nmin.EQ.2.AND.nmax.EQ.5.AND.nsum.EQ.10) np = n3
5189 IF (nmin.EQ.3.AND.nmax.EQ.6.AND.nsum.EQ.13) np = n4
5190 IF (np.EQ.0.OR.angmx3.LT.120.) go to 110
5194 75
IF (np.EQ.n4) go to 80
5202 80 CALL
triswp(n4,n1,n2,n3,
j,nfail, &
5203 x,ndc,nbh,iprot,ncell,ndg,idgp,ndgp,nedge, &
5204 vol,xcen,ycen,zcen,rc,rat,nptet,nacpt,tolv)
5206 IF (nfail.EQ.0)
THEN
5214 110
IF (
max(angmx1,angmx2).LT.120.) go to 135
5218 115
IF (angmax.LT.120.) go to 130
5222 IF (
i.EQ.1) go to 120
5226 IF (
i.EQ.2) go to 120
5230 IF (
i.EQ.3) go to 120
5234 IF (
i.EQ.4) go to 120
5238 IF (
i.EQ.5) go to 120
5242 120 m2 = n1 +n2 +n3 +n4 -na -nb -m1
5244 CALL
replace(m1,m2,na,nb,
j,nring,npass,nfail, &
5245 x,dens,ityp,nnode,ndc,nbh,iprot,ncell, &
5246 ndg,idgp,ndgp,nedge,ipoint,npoint, &
5247 vol,xcen,ycen,zcen,rc,rat,nptet,nacpt, &
5248 noctr,ioctr,nlink,nref,xfar,yfar,zfar, &
5249 iring,ntetkp,lnbr,ishk,mnbr,kshk, &
5250 nfad,nbhkp,iedkp,newc,
dx,
dy,
dz,ds,npp,vlt,tolv)
5252 IF (nfail.EQ.0) go to 125
5253 IF (npass.EQ.1.OR.nfail.EQ.2) go to 128
5260 128
IF (loop.LE.4)
THEN
5263 x,ityp,nnode,ndc,nbh,iprot,ncell, &
5264 ndg,idgp,ndgp,nedge, &
5265 vol,xcen,ycen,zcen,rc,rat,dens,nptet,nacpt, &
5266 idone,nref,nlink,noctr,ioctr,xfar,yfar,zfar, &
5267 xoctr,yoctr,zoctr,xhold,yhold,zhold, &
5268 xkeep,ykeep,zkeep,ksrch,nsrch, &
5269 ipoint,npoint,iflag,nflag,nfill,newcel,ntri, &
5270 ncav,nshake,newnbh,nold,ncavfc,ikeep,ldel, &
5271 nedgrm,xc,yc,zc,
v,rad,rcrin,lnkup,lnkdn, &
5274 IF (nfail.EQ.0) go to 125
5277 IF (l.GT.nwait) go to 135
5281 950
FORMAT(i6,
' SLIVERS REMOVED')
5286 IF (nchk.EQ.0) go to 145
5290 IF (iprot(l).EQ.1) go to 160
5291 IF (nbh(1,l).EQ.0) go to 160
5296 IF (ityp(
n).LT.0) go to 165
5299 WRITE (6,600) npt,ntet
5301 600
FORMAT(/
'MESH OPTIMIZATION COMPLETE'/ &
5302 5
x,i7,
' MESH POINTS'/ &
5303 5
x,i7,
' MESH CELLS'/)
5304 610
FORMAT(/
'BEGIN MESH OPTIMIZATION')
5314 SUBROUTINE replace (N1,N2,NA,NB,J,NRING,NPASS,NFAIL, &
5315 x,dens,ityp,nnode,ndc,nbh,iprot,ncell, &
5316 ndg,idgp,ndgp,nedge,ipoint,npoint, &
5317 vol,xcen,ycen,zcen,rc,rat,nptet,nacpt, &
5318 noctr,ioctr,nlink,nref,xfar,yfar,zfar, &
5319 iring,ntetkp,lnbr,ishk,mnbr,kshk, &
5320 nfad,nbhkp,iedkp,newc,
dx,
dy,
dz,ds,np,vlt,tolv)
5337 INTEGER :: ioctr,
j,na,nb,ncell,nedge,nfail,nnode,npass,nring,n1,n2
5338 INTEGER :: ityp(*),npoint(*),ipoint(*),nptet(*),nacpt(*)
5339 INTEGER :: ndc(4,*),nbh(4,*),iprot(*),ndg(2,*),idgp(*),ndgp(*)
5340 INTEGER :: nref(*),nlink(*),noctr(2,*)
5341 INTEGER :: iring(*),ntetkp(*),nfad(3,*),newc(*), &
5342 nbhkp(3,*),iedkp(4,*),lnbr(*),ishk(*), &
5343 mnbr(*),kshk(*),np(*)
5344 DOUBLE PRECISION :: tolv
5345 DOUBLE PRECISION ::
x(3,*),dens(*)
5346 DOUBLE PRECISION :: vol(*),xcen(*),ycen(*),zcen(*),rc(*),rat(*)
5347 DOUBLE PRECISION :: xfar(2),yfar(2),zfar(2)
5348 DOUBLE PRECISION ::
dx(*),
dy(*),
dz(*),ds(*),vlt(*)
5350 INTEGER ::
i,iedg,ip1,ismall,jnext,jpre,j1,j2,
k,kcnt,kend,khalf, &
5351 kk,kmax,kmin,km1,kplus,kplusm,kplusp,kp1,ksum,k1,k2, &
5352 k3,k4,l,loop,l1,l2,
m,mend,mmax,mmin,msum,m1,m2,m3,m4, &
5353 m5,
n,nedk,npnext,npre,nrngp1
5354 DOUBLE PRECISION :: area,a1,a2,a3,b1,b2,b3,radc,radmax,ratmax, &
5355 rcrin,rnx,rny,rnz,vl,vmin,vorig,vtet1,vtet2, &
5368 CALL
neighb(
j,n1,na,nb,jnext,k1,k2,ndc,nbh)
5370 npnext = ndc(1,jnext) +ndc(2,jnext) +ndc(3,jnext) &
5371 +ndc(4,jnext) -n1 -na -nb
5377 CALL
neighb(jpre,npre,na,nb,jnext,k1,k2,ndc,nbh)
5379 npnext = ndc(1,jnext) +ndc(2,jnext) +ndc(3,jnext) &
5380 +ndc(4,jnext) -npre -na -nb
5382 IF (nring.GT.mxring) go to 310
5384 ntetkp(nring) = jnext
5385 IF (npnext.EQ.n1) go to 300
5386 IF (npnext.NE.n2) go to 10
5396 IF (iprot(ntetkp(
i)).EQ.1) go to 200
5397 radmax =
max(radmax,rc(ntetkp(
i)))
5398 ratmax =
max(ratmax,rat(ntetkp(
i)))
5399 ip1 = mod(
i,nring) +1
5401 CALL
neighb(ntetkp(
i),iring(
i),iring(ip1),na,lnbr(
i),kk,ishk(
i), &
5404 CALL
neighb(ntetkp(
i),iring(
i),iring(ip1),nb,mnbr(
i),kk,kshk(
i), &
5410 CALL
tessel(iring,nring,nfad,nbhkp,iedkp,nedk,0,0, &
5411 x,ipoint,npoint,
dx,
dy,
dz,ds,np)
5420 CALL
circum(
x,m1,m2,m3,na,xcn,ycn,zcn,vl,radc,ismall,tolv)
5422 IF (radc.GT.radmax) go to 200
5424 rcrin = radc*area/vl
5425 IF (rcrin.GT.ratmax) go to 200
5427 CALL
circum(
x,m1,m2,m3,nb,xcn,ycn,zcn,vl,radc,ismall,tolv)
5429 IF (radc.GT.radmax) go to 200
5431 rcrin = radc*area/vl
5432 IF (rcrin.GT.ratmax) go to 200
5433 rnx =
cofact(
x(2,m1),
x(2,m2),
x(2,m3),
x(3,m1),
x(3,m2),
x(3,m3))
5434 rny =
cofact(
x(3,m1),
x(3,m2),
x(3,m3),
x(1,m1),
x(1,m2),
x(1,m3))
5435 rnz =
cofact(
x(1,m1),
x(1,m2),
x(1,m3),
x(2,m1),
x(2,m2),
x(2,m3))
5436 vtet1 = rnx*(
x(1,na) -
x(1,m1)) +rny*(
x(2,na) -
x(2,m1)) &
5437 +rnz*(
x(3,na) -
x(3,m1))
5438 vtet2 = rnx*(
x(1,nb) -
x(1,m1)) +rny*(
x(2,nb) -
x(2,m1)) &
5439 +rnz*(
x(3,nb) -
x(3,m1))
5440 IF (abs(vtet1).LT.vorig.OR.abs(vtet2).LT.vorig) go to 100
5441 IF (vtet1*vtet2.GE.0.0) go to 100
5447 IF (l.LT.0) go to 38
5449 IF (loop.EQ.2) j2 = mnbr(l)
5450 IF (iprot(j2).EQ.1) go to 38
5455 IF (loop.EQ.2) k4 = nb
5462 kmin =
min(k1,k2,k3)
5463 kmax =
max(k1,k2,k3)
5465 32 mmin =
min(m1,m2,m3)
5466 mmax =
max(m1,m2,m3)
5468 IF (mmin.EQ.kmin.AND.mmax.EQ.kmax.AND.msum.EQ.ksum) go to 34
5469 IF (m1.EQ.mend) go to 33
5476 33
IF (k1.EQ.kend) go to 340
5484 rnx =
cofact(
x(2,m1),
x(2,m2),
x(2,m3),
x(3,m1),
x(3,m2),
x(3,m3))
5485 rny =
cofact(
x(3,m1),
x(3,m2),
x(3,m3),
x(1,m1),
x(1,m2),
x(1,m3))
5486 rnz =
cofact(
x(1,m1),
x(1,m2),
x(1,m3),
x(2,m1),
x(2,m2),
x(2,m3))
5487 vtet1 = rnx*(
x(1,m4) -
x(1,m1)) +rny*(
x(2,m4) -
x(2,m1)) &
5488 +rnz*(
x(3,m4) -
x(3,m1))
5489 vtet2 = rnx*(
x(1,m5) -
x(1,m1)) +rny*(
x(2,m5) -
x(2,m1)) &
5490 +rnz*(
x(3,m5) -
x(3,m1))
5491 IF (vtet1*vtet2.GE.0.) go to 100
5499 IF (
k.GT.nring) go to 35
5511 ndc(1,newc(
k)) = nfad(1,
k)
5512 ndc(2,newc(
k)) = nfad(2,
k)
5513 ndc(3,newc(
k)) = nfad(3,
k)
5516 ndc(1,newc(kplus)) = nfad(1,
k)
5517 ndc(2,newc(kplus)) = nfad(2,
k)
5518 ndc(3,newc(kplus)) = nfad(3,
k)
5519 ndc(4,newc(kplus)) = nb
5525 nbh(1,ntetkp(
k)) = 0
5531 IF (l.GT.0) go to 80
5533 nbh(
i,newc(
k)) = newc(l)
5534 nbh(
i,newc(kplus)) = newc(l+khalf)
5536 80 nbh(
i,newc(
k)) = lnbr(l)
5537 nbh(
i,newc(kplus)) = mnbr(l)
5538 nbh(ishk(l),lnbr(l)) = newc(
k)
5539 nbh(kshk(l),mnbr(l)) = newc(kplus)
5541 nbh(4,newc(
k)) = newc(kplus)
5542 nbh(4,newc(kplus)) = newc(
k)
5547 IF (nring.EQ.nedk) go to 145
5555 IF (iedg.EQ.0) go to 94
5556 92
IF (ndgp(iedg).EQ.0) go to 94
5560 IF (nedge.GT.mxedge) go to 370
5564 IF (iedg.EQ.0) idgp(l1) = nedge
5565 IF (iedg.GT.0) ndgp(iedg) = nedge
5572 100 CALL
neighb(
j,n1,n2,na,j1,k1,k2,ndc,nbh)
5574 CALL
neighb(
j,n1,n2,nb,j2,k1,k2,ndc,nbh)
5576 IF (iprot(j1).EQ.0.AND.iprot(j2).EQ.0.AND.npass.EQ.0) go to 210
5578 CALL
sliver(iring,nring,na,nb,nfail,vorig, &
5579 x,dens,ityp,ipoint,nnode, &
5580 noctr,ioctr,nlink,nref,xfar,yfar,zfar,vlt)
5582 IF (nfail.NE.0)
RETURN
5588 IF (
k.GT.nring) go to 105
5592 105 ncell = ncell +1
5600 kp1 = mod(
k,nring) +1
5605 ndc(3,newc(
k)) = nnode
5609 ndc(1,newc(kplus)) = m1
5610 ndc(2,newc(kplus)) = m2
5611 ndc(3,newc(kplus)) = nnode
5612 ndc(4,newc(kplus)) = nb
5613 iprot(newc(kplus)) = 0
5619 nbh(1,ntetkp(
k)) = 0
5622 km1 = mod(nring-2+
k,nring) +1
5623 kp1 = mod(
k,nring) +1
5627 nbh(1,newc(
k)) = newc(km1)
5628 nbh(2,newc(
k)) = newc(kp1)
5629 nbh(3,newc(
k)) = lnbr(
k)
5630 nbh(4,newc(
k)) = newc(kplus)
5631 nbh(ishk(
k),lnbr(
k)) = newc(
k)
5632 nbh(1,newc(kplus)) = newc(kplusm)
5633 nbh(2,newc(kplus)) = newc(kplusp)
5634 nbh(3,newc(kplus)) = mnbr(
k)
5635 nbh(4,newc(kplus)) = newc(
k)
5636 nbh(kshk(
k),mnbr(
k)) = newc(kplus)
5642 iedg = idgp(iring(
k))
5643 IF (iedg.EQ.0) go to 134
5644 132
IF (ndgp(iedg).EQ.0) go to 134
5647 134 nedge = nedge +1
5648 IF (nedge.GT.mxedge) go to 370
5649 ndg(1,nedge) = iring(
k)
5650 ndg(2,nedge) = nnode
5652 IF (iedg.EQ.0) idgp(iring(
k)) = nedge
5653 IF (iedg.GT.0) ndgp(iedg) = nedge
5656 IF (iedg.EQ.0) go to 140
5657 138
IF (ndgp(iedg).EQ.0) go to 140
5660 140 nedge = nedge +1
5661 IF (nedge.GT.mxedge) go to 370
5663 ndg(2,nedge) = nnode
5665 IF (iedg.EQ.0) idgp(na) = nedge
5666 IF (iedg.GT.0) ndgp(iedg) = nedge
5668 IF (iedg.EQ.0) go to 144
5669 142
IF (ndgp(iedg).EQ.0) go to 144
5672 144 nedge = nedge +1
5673 IF (nedge.GT.mxedge) go to 370
5675 ndg(2,nedge) = nnode
5677 IF (iedg.EQ.0) idgp(nb) = nedge
5678 IF (iedg.GT.0) ndgp(iedg) = nedge
5685 146
IF (l2.EQ.ndg(2,iedg)) go to 148
5689 148 CALL
edgerm(iedg,ndg,idgp,ndgp)
5699 CALL
circum(
x,m1,m2,m3,m4,xcn,ycn,zcn,vl,radc,ismall,tolv)
5701 IF (
k.EQ.1) vmin = vl
5702 IF (
k.GT.1) vmin =
min(vmin,vl)
5704 IF (ismall.EQ.1)
WRITE (6,960)
k,vl,radc, &
5705 m1,
x(1,m1),
x(2,m1),
x(3,m1), &
5706 m2,
x(1,m2),
x(2,m2),
x(3,m2), &
5707 m3,
x(1,m3),
x(2,m3),
x(3,m3), &
5708 m4,
x(1,m4),
x(2,m4),
x(3,m4)
5709 960
FORMAT(
'RING ',i2,
' VOLUME = ',e13.5,
' RC = ',e13.5/ &
5710 'M1 ',i6,
' X = ',f6.3,
' Y = ',f6.3,
' Z = ',f6.3/ &
5711 'M2 ',i6,
' X = ',f6.3,
' Y = ',f6.3,
' Z = ',f6.3/ &
5712 'M3 ',i6,
' X = ',f6.3,
' Y = ',f6.3,
' Z = ',f6.3/ &
5713 'M4 ',i6,
' X = ',f6.3,
' Y = ',f6.3,
' Z = ',f6.3)
5715 IF (ismall.EQ.1) go to 320
5722 rat(newc(
k)) = rc(newc(
k))*area/vol(newc(
k))
5759 600
FORMAT(///5
x,
'SEARCH FOR TETRAHEDRAL RING HAS RETURNED TO THE'/ &
5760 5
x,
'STARTING FACE. THIS INDICATES AN INCONSISTENCY IN'/ &
5761 5
x,
'THE TETRAHEDRAL ENSEMBLE.'/ &
5762 5
x,
'PROGRAM STOPPED IN REPLACE')
5763 610
FORMAT(///5
x,
'DIMENSION OF ARRAY IRING EXCEEDED. INCREASE SIZE'/ &
5764 5
x,
'OF MRING. PROGRAM STOPPED IN REPLACE')
5765 620
FORMAT(///5
x,
'AT LEAST ONE NEW TETRAHEDRON HAS TOO SMALL A VOLUME' &
5766 /5
x,
'PROGRAM STOPPED IN REPLACE')
5767 640
FORMAT(///5
x,
'UNABLE TO FIND COMMON FACE BETWEEN ADJACENT CELLS' &
5768 /5
x,
'PROGRAM STOPPED IN REPLACE')
5769 660
FORMAT(//5
x,
'UNABLE TO FIND EDGE ADDRESS FOR A NEW TETRAHEDRON'/ &
5770 5
x,
'PROGRAM STOPPED IN ROUTINE REPLACE')
5771 670
FORMAT(//5
x,
'DIMENSION OF NDG ARRAY EXCEEDED IN ROUTINE REPLACE.'/ &
5772 5
x,
'INCREASE SIZE OF MBPTS.')
5784 x,ndc,nbh,iprot,ncell,ndg,idgp,ndgp,nedge, &
5785 vol,xcen,ycen,zcen,rc,rat,nptet,nacpt,tolv)
5801 INTEGER ::
j,nedge,nfail,npt,n1,n2,n3
5802 INTEGER :: nptet(*),nacpt(*),ndc(4,*),nbh(4,*),iprot(*),ndg(2,*), &
5804 DOUBLE PRECISION :: tolv
5805 DOUBLE PRECISION ::
x(3,*),vol(*),xcen(*),ycen(*),zcen(*), &
5808 INTEGER ::
i,iedg,im1,ip1,ismall,
k,kk,km1,kp1,k1,k2,la,l1,l2, &
5809 m1,m2,m3,m4,na,ncell,nold
5810 INTEGER :: iring(3),newc(3),lnbr(3),ishk(3),mnbr(3),kshk(3)
5811 DOUBLE PRECISION :: area,a1,a2,a3,b1,b2,b3,radmax,ratmax,rcrin, &
5812 rnx,rny,rnz,xcn,ycn,zcn,vl,vtet1,vtet2,radc
5820 CALL
neighb(
j,n1,n2,n3,la,k1,k2,ndc,nbh)
5822 IF (iprot(la).EQ.1)
RETURN
5823 na = ndc(1,la) +ndc(2,la) +ndc(3,la) +ndc(4,la) &
5831 radmax =
max(rc(
j),rc(la))
5832 ratmax =
max(rat(
j),rat(la))
5840 CALL
circum(
x,m1,m2,npt,na,xcn,ycn,zcn,vl,radc,ismall,tolv)
5842 IF (radc.GT.radmax) go to 300
5844 rcrin = radc*area/vl
5845 IF (rcrin.GT.ratmax) go to 300
5846 rnx =
cofact(
x(2,m1),
x(2,m2),
x(2,na),
x(3,m1),
x(3,m2),
x(3,na))
5847 rny =
cofact(
x(3,m1),
x(3,m2),
x(3,na),
x(1,m1),
x(1,m2),
x(1,na))
5848 rnz =
cofact(
x(1,m1),
x(1,m2),
x(1,na),
x(2,m1),
x(2,m2),
x(2,na))
5849 vtet1 = rnx*(
x(1,m3) -
x(1,m1)) +rny*(
x(2,m3) -
x(2,m1)) &
5850 +rnz*(
x(3,m3) -
x(3,m1))
5851 vtet2 = rnx*(
x(1,npt) -
x(1,m1)) +rny*(
x(2,npt) -
x(2,m1)) &
5852 +rnz*(
x(3,npt) -
x(3,m1))
5853 IF (vtet1*vtet2.LT.0) go to 310
5854 IF (abs(vtet1).LT.10.*tolv) go to 310
5855 IF (abs(vtet2).LT.10.*tolv) go to 310
5863 CALL
neighb(
j,iring(
k),iring(kp1),npt,lnbr(
k),kk,ishk(
k),ndc,nbh)
5865 CALL
neighb(la,iring(
k),iring(kp1),na,mnbr(
k),kk,kshk(
k),ndc,nbh)
5874 IF (ncell.GT.mxcell) go to 360
5881 ndc(1,newc(
k)) = npt
5882 ndc(2,newc(
k)) = iring(
k)
5883 ndc(3,newc(
k)) = iring(kp1)
5897 CALL
circum(
x,m1,m2,m3,m4,xcn,ycn,zcn,vl,radc,ismall,tolv)
5900 IF (ismall.EQ.1)
WRITE (6,910) npt,na,(iring(kk),kk=1,3)
5901 910
FORMAT(
'NPT ',i6,
' NA ',i6,
' IRING ',4i6)
5902 IF (ismall.EQ.1)
WRITE (6,960)
k,vl, &
5903 m1,
x(1,m1),
x(2,m1),
x(3,m1), &
5904 m2,
x(1,m2),
x(2,m2),
x(3,m2), &
5905 m3,
x(1,m3),
x(2,m3),
x(3,m3), &
5906 m4,
x(1,m4),
x(2,m4),
x(3,m4)
5907 960
FORMAT(
'RING ',i2,
' VOLUME = ',e13.5/ &
5908 'M1 ',i6,
' X = ',f10.4,
' Y = ',f10.4,
' Z = ',f10.4/ &
5909 'M2 ',i6,
' X = ',f10.4,
' Y = ',f10.4,
' Z = ',f10.4/ &
5910 'M3 ',i6,
' X = ',f10.4,
' Y = ',f10.4,
' Z = ',f10.4/ &
5911 'M4 ',i6,
' X = ',f10.4,
' Y = ',f10.4,
' Z = ',f10.4)
5913 IF (ismall.EQ.1) go to 350
5920 rat(newc(
k)) = rc(newc(
k))*area/vol(newc(
k))
5928 nbh(ishk(
k),lnbr(
k)) = newc(
k)
5929 nbh(kshk(
k),mnbr(
k)) = newc(
k)
5930 nbh(1,newc(
k)) = newc(km1)
5931 nbh(2,newc(
k)) = newc(kp1)
5932 nbh(3,newc(
k)) = lnbr(
k)
5933 nbh(4,newc(
k)) = mnbr(
k)
5941 IF (iedg.EQ.0) go to 60
5942 55
IF (ndgp(iedg).EQ.0) go to 60
5946 IF (nedge.GT.mxedge) go to 330
5950 IF (iedg.EQ.0) idgp(l1) = nedge
5951 IF (iedg.GT.0) ndgp(iedg) = nedge
5955 nptet(npt) = newc(1)
5958 nptet(iring(
k)) = newc(
k)
5967 320
WRITE (6,620) nold,l1,l2,npt,na
5975 600
FORMAT(5
x,
'A NEW TETRAHEDRON IS WORSE THAN AN EXISTING ONE')
5976 610
FORMAT(5
x,
'TETRAHEDRAL RING IS NON-CONVEX')
5977 620
FORMAT(///5
x,
'EDGE TO BE INSERTED ALREADY EXISTS'/ &
5978 5
x,
'EDGE ',i6,
' VERTICES ',2i6,
' NPT ',i6,
' NA ',i6/ &
5979 //5
x,
'PROGRAM STOPPED IN TRISWP')
5980 630
FORMAT(///5
x,
'DIMENSION OF ARRAY NDG EXCEEDED IN ROUTINE TRISWP.'/ &
5981 //5
x,
'INCREASE SIZE OF MNODE.'/ &
5982 //5
x,
'PROGRAM STOPPED IN TRISWP')
5983 650
FORMAT(///5
x,
'AT LEAST ONE OF THE NEW TETRAHEDRA HAS TOO'/ &
5984 5
x,
'SMALL A VOLUME.'/ &
5985 //5
x,
'PROGRAM STOPPED IN TRISWP')
5986 660
FORMAT(///5
x,
'DIMENSION OF ARRAY NDC EXCEEDS ALLOWED VALUE.'/ &
5987 //5
x,
'INCREASE SIZE OF MCELL.'/ &
5988 //5
x,
'PROGRAM STOPPED IN TRISWP')
5999 SUBROUTINE tessel (IRING,NRING,NFAD,NBHKP,IEDKP,NEDK,IMEET,ITOUCH, &
6000 x,ipoint,npoint,
dx,
dy,
dz,ds,np)
6015 INTEGER :: imeet,itouch,nedk,nring
6016 INTEGER :: npoint(*),ipoint(*)
6017 INTEGER :: iring(*),nfad(3,*),nbhkp(3,*),iedkp(4,*),np(*)
6018 DOUBLE PRECISION ::
x(3,*)
6019 DOUBLE PRECISION ::
dx(*),
dy(*),
dz(*),ds(*)
6021 INTEGER ::
i,ia,ib,imin,imin1,im1,ip1,
k,kpre,l,la,lb,ncnt,nedtot, &
6033 ipoint(iring(
i)) = 0
6037 ip1 = mod(
i,nring) +1
6038 ia =
min(iring(
i),iring(ip1))
6039 ib =
max(iring(
i),iring(ip1))
6045 6
IF (
k.EQ.0) go to 8
6049 8
IF (ipoint(ia).NE.0) npoint(kpre) =
i
6050 IF (ipoint(ia).EQ.0) ipoint(ia) =
i
6059 20
IF (num.EQ.3) go to 60
6063 im1 = mod(num+imin-2,num) +1
6064 ip1 = mod(imin,num) +1
6066 nfad(1,nswit) = iring(imin)
6067 nfad(2,nswit) = iring(ip1)
6068 nfad(3,nswit) = iring(im1)
6070 ia =
min(iring(im1),iring(ip1))
6071 ib =
max(iring(im1),iring(ip1))
6074 iedkp(3,nedk) = -nswit
6077 22
IF (
k.EQ.0) go to 24
6081 24
IF (ipoint(ia).NE.0) npoint(kpre) = nedk
6082 IF (ipoint(ia).EQ.0) ipoint(ia) = nedk
6083 ia =
min(iring(im1),iring(imin))
6084 ib =
max(iring(im1),iring(imin))
6086 25
IF (ib.EQ.iedkp(2,
k)) go to 30
6088 IF (
k.EQ.0) go to 330
6090 30 iedkp(4,
k) = -nswit
6091 ia =
min(iring(ip1),iring(imin))
6092 ib =
max(iring(ip1),iring(imin))
6094 35
IF (ib.EQ.iedkp(2,
k)) go to 40
6096 IF (
k.EQ.0) go to 340
6098 40 iedkp(4,
k) = -nswit
6101 IF (
i.EQ.imin) go to 50
6103 iring(ncnt) = iring(
i)
6111 IF (nswit.NE.nring-2) go to 320
6112 nfad(1,nswit) = iring(1)
6113 nfad(2,nswit) = iring(2)
6114 nfad(3,nswit) = iring(3)
6115 ia =
min(iring(1),iring(2))
6116 ib =
max(iring(1),iring(2))
6118 65
IF (ib.EQ.iedkp(2,
k)) go to 70
6121 70 iedkp(4,
k) = -nswit
6122 ia =
min(iring(2),iring(3))
6123 ib =
max(iring(2),iring(3))
6125 75
IF (ib.EQ.iedkp(2,
k)) go to 80
6128 80 iedkp(4,
k) = -nswit
6129 ia =
min(iring(3),iring(1))
6130 ib =
max(iring(3),iring(1))
6132 85
IF (ib.EQ.iedkp(2,
k)) go to 90
6135 90 iedkp(4,
k) = -nswit
6147 la = iabs(iedkp(3,
k))
6148 lb = iabs(iedkp(4,
k))
6149 IF (la.EQ.0.OR.lb.EQ.0) go to 300
6150 npoint(lb) = npoint(lb) +1
6151 nbhkp(npoint(lb),lb) = iedkp(3,
k)
6152 IF (iedkp(3,
k).GT.0) go to 120
6153 npoint(la) = npoint(la) +1
6154 nbhkp(npoint(la),la) = iedkp(4,
k)
6157 IF (npoint(l).NE.3) go to 310
6161 300
WRITE (6,600) la,lb
6165 320
WRITE (6,620) nring,nswit
6167 330
WRITE (6,630) ia,ib
6169 340
WRITE (6,640) ia,ib
6171 600
FORMAT(/5
x,
'A ZERO ADDRESS HAS BEEN FOUND FOR A TRIANGLE', &
6172 5
x,
'ADDRESSES ARE ',i6,
' AND ',i6, &
6173 5
x,
'PROGRAM STOPPED IN ROUTINE TESSEL')
6174 610
FORMAT(/5
x,
'INCORRECT ASSEMBLY OF TRIANGLES', &
6175 5
x,
'PROGRAM STOPPED IN ROUTINE TESSEL')
6176 620
FORMAT(/5
x,
'NRING IS NOT EQUAL TO NSWIT PLUS TWO', &
6177 5
x,
'NRING = ',i3,
' NSWIT = ',i3, &
6178 5
x,
'PROGRAM STOPPED IN ROUTINE TESSEL')
6179 630
FORMAT(/5
x,
'EDGE WITH ZERO ADDRESS FOUND AFTER LABEL 25', &
6180 5
x,
'IA = ',i6,
' IB = ',i6, &
6181 5
x,
'PROGRAM STOPPED IN ROUTINE TESSEL')
6182 640
FORMAT(/5
x,
'EDGE WITH ZERO ADDRESS FOUND AFTER LABEL 35', &
6183 5
x,
'IA = ',i6,
' IB = ',i6, &
6184 5
x,
'PROGRAM STOPPED IN ROUTINE TESSEL')
6197 SUBROUTINE angfnd (NUM,IRING,IMIN,IMEET,ITOUCH,X,DX,DY,DZ,DS)
6213 INTEGER :: imeet,imin,itouch,num
6215 DOUBLE PRECISION ::
x(3,*)
6216 DOUBLE PRECISION ::
dx(*),
dy(*),
dz(*),ds(*)
6218 INTEGER ::
i,im1,ip1
6219 DOUBLE PRECISION :: ang,angmin,beta,cang,det,
div,ds1sq,ds2sq,fac, &
6220 fact1,fact2,gamma,
pi,prod,rnx,rpx,rny,rpy, &
6229 im1 = mod(num+
i-2,num) +1
6230 dx(
i) =
x(1,iring(im1)) -
x(1,iring(
i))
6231 dy(
i) =
x(2,iring(im1)) -
x(2,iring(
i))
6232 dz(
i) =
x(3,iring(im1)) -
x(3,iring(
i))
6241 rpx = rpx -fac*(
dy(ip1)*
dz(
i) -
dy(
i)*
dz(ip1))
6242 rpy = rpy -fac*(
dz(ip1)*
dx(
i) -
dz(
i)*
dx(ip1))
6243 rpz = rpz -fac*(
dx(ip1)*
dy(
i) -
dx(
i)*
dy(ip1))
6245 fac = 1./
sqrt(rpx*rpx +rpy*rpy +rpz*rpz)
6251 im1 = mod(num+
i-2,num) +1
6253 fact1 = ds(
i)*ds(ip1) -prod
6254 fact2 = ds(
i)*ds(ip1) +prod
6255 det = 1./(fact1*fact2)
6257 ds2sq = ds(ip1)*ds(ip1)
6258 beta = ds2sq*(rpx*
dx(
i) +rpy*
dy(
i) +rpz*
dz(
i)) &
6259 +prod*(rpx*
dx(ip1) +rpy*
dy(ip1) +rpz*
dz(ip1))*det
6260 gamma = ds1sq*(rpx*
dx(ip1) +rpy*
dy(ip1) +rpz*
dz(ip1)) &
6261 +prod*(rpx*
dx(
i) +rpy*
dy(
i) +rpz*
dz(
i))*det
6262 rnx = rpx -beta*
dx(
i) -gamma*
dx(ip1)
6263 rny = rpy -beta*
dy(
i) -gamma*
dy(ip1)
6264 rnz = rpz -beta*
dz(
i) -gamma*
dz(ip1)
6265 fac = 1./
sqrt(rnx*rnx +rny*rny +rnz*rnz)
6274 ang =
atan2(sang,cang)
6275 IF (ang.LT.0.) ang = ang +2.*
pi
6276 IF (imin.EQ.0) go to 25
6277 IF (ang.GT.angmin) go to 30
6278 IF (iring(
i).EQ.imeet.AND.iring(im1).EQ.itouch) go to 25
6279 IF (iring(
i).EQ.imeet.AND.iring(ip1).EQ.itouch) go to 25
6280 IF (iring(
i).EQ.itouch.AND.iring(im1).EQ.imeet) go to 25
6281 IF (iring(
i).EQ.itouch.AND.iring(ip1).EQ.imeet) go to 25
6282 IF (iring(
i).EQ.imeet.OR.iring(
i).EQ.itouch) go to 30
6297 SUBROUTINE sliver (IRING,NRING,NA,NB,NFAIL,VORIG, &
6298 x,dens,ityp,ipoint,nnode, &
6299 noctr,ioctr,nlink,nref,xfar,yfar,zfar,vlt)
6315 INTEGER :: ioctr,na,nb,nfail,nnode,nring
6316 INTEGER :: ipoint(*),iring(*),ityp(*),nref(*),nlink(*),noctr(2,*)
6317 DOUBLE PRECISION :: vorig
6318 DOUBLE PRECISION ::
x(3,*),dens(*)
6319 DOUBLE PRECISION :: xfar(2),yfar(2),zfar(2),vlt(*)
6321 INTEGER ::
i,im1,m1,m2,
n,ncnt,ncyc,nrm1
6322 DOUBLE PRECISION :: a1,a2,a3,b1,b2,b3,densum,fac,rnx,rny,rnz,tol, &
6323 vtet1,vtet2,xpt,ypt,zpt
6334 xpt = xpt +
x(1,iring(
i))
6335 ypt = ypt +
x(2,iring(
i))
6336 zpt = zpt +
x(3,iring(
i))
6337 densum = densum +dens(iring(
i))
6339 fac = 1./float(nring-2)
6340 xpt = .5*fac*xpt +.25*(
x(1,na) +
x(1,nb))
6341 ypt = .5*fac*ypt +.25*(
x(2,na) +
x(2,nb))
6342 zpt = .5*fac*zpt +.25*(
x(3,na) +
x(3,nb))
6343 densum = .5*fac*densum +.25*(dens(na) +dens(nb))
6347 im1 = mod(nring-2+
i,nring) +1
6350 rnx =
cofact(
x(2,m1),
x(2,m2),
x(2,na),
x(3,m1),
x(3,m2),
x(3,na))
6351 rny =
cofact(
x(3,m1),
x(3,m2),
x(3,na),
x(1,m1),
x(1,m2),
x(1,na))
6352 rnz =
cofact(
x(1,m1),
x(1,m2),
x(1,na),
x(2,m1),
x(2,m2),
x(2,na))
6353 vtet2 = rnx*(xpt -
x(1,m1)) +rny*(ypt -
x(2,m1)) &
6357 IF (abs(vtet2).LT.vorig) go to 210
6358 IF (
i.EQ.1) go to 15
6359 IF (vtet1*vtet2.LE.tol) go to 100
6364 im1 = mod(nring-2+
i,nring) +1
6367 rnx =
cofact(
x(2,m1),
x(2,m2),
x(2,nb),
x(3,m1),
x(3,m2),
x(3,nb))
6368 rny =
cofact(
x(3,m1),
x(3,m2),
x(3,nb),
x(1,m1),
x(1,m2),
x(1,nb))
6369 rnz =
cofact(
x(1,m1),
x(1,m2),
x(1,nb),
x(2,m1),
x(2,m2),
x(2,nb))
6370 vtet2 = rnx*(xpt -
x(1,m1)) +rny*(ypt -
x(2,m1)) &
6374 IF (abs(vtet2).LT.vorig) go to 210
6375 IF (vtet1*vtet2.GE.tol) go to 100
6378 IF (nnode.GT.mxnode) go to 230
6387 CALL
octfil(
n,
x,noctr,ioctr,nlink,nref,xfar,yfar,zfar)
6390 100
IF (ncyc.EQ.2) go to 200
6392 xpt = .5*xpt +.25*(
x(1,na) +
x(1,nb))
6393 ypt = .5*ypt +.25*(
x(2,na) +
x(2,nb))
6394 zpt = .5*zpt +.25*(
x(3,na) +
x(3,nb))
6395 densum = .5*densum +.25*(dens(na) +dens(nb))
6405 600
FORMAT(//
'POINT TO BE INSERTED LIES OUTSIDE TETRAHEDRAL RING')
6406 610
FORMAT(//
'A NEW TETRAHEDRON IS TOO SMALL'/ &
6407 'VORIG = ',e13.5,
' VTET = ',e13.5)
6408 630
FORMAT(//5
x,
'NNODE EXCEEDS DIMENSION OF ARRAY X.'/ &
6409 5
x,
'INCREASE SIZE OF MNODE.'/ &
6410 5
x,
'PROGRAM STOPPED IN ROUTINE SLIVER')
6421 SUBROUTINE circum (X,N1,N2,N3,N4,XCN,YCN,ZCN,VTET,RCIR,IFLAG,TOL)
6439 INTEGER :: iflag,n1,n2,n3,n4
6440 DOUBLE PRECISION :: xcn,ycn,zcn,vtet,rcir,tol
6441 DOUBLE PRECISION ::
x(3,*)
6443 DOUBLE PRECISION :: a1,a2,a3,b1,b2,b3,c5,c6,c9,det,fac,h1,h2,h3, &
6444 rnx,rny,rnz,vcell,xshf,yshf,zshf
6452 rnx =
cofact(
x(2,n1),
x(2,n2),
x(2,n3),
x(3,n1),
x(3,n2),
x(3,n3))
6453 rny =
cofact(
x(3,n1),
x(3,n2),
x(3,n3),
x(1,n1),
x(1,n2),
x(1,n3))
6454 rnz =
cofact(
x(1,n1),
x(1,n2),
x(1,n3),
x(2,n1),
x(2,n2),
x(2,n3))
6455 det = rnx*(
x(1,n4) -
x(1,n1)) +rny*(
x(2,n4) -
x(2,n1)) &
6456 +rnz*(
x(3,n4) -
x(3,n1))
6462 IF (vtet.LT.tol) go to 30
6464 xshf = .5*(
x(1,n1) +
x(1,n2) +
x(1,n3) +
x(1,n4))
6465 yshf = .5*(
x(2,n1) +
x(2,n2) +
x(2,n3) +
x(2,n4))
6466 zshf = .5*(
x(3,n1) +
x(3,n2) +
x(3,n3) +
x(3,n4))
6467 h1 = (
x(1,n4) -
x(1,n1))*(
x(1,n4) -xshf +
x(1,n1)) &
6468 +(
x(2,n4) -
x(2,n1))*(
x(2,n4) -yshf +
x(2,n1)) &
6469 +(
x(3,n4) -
x(3,n1))*(
x(3,n4) -zshf +
x(3,n1))
6470 h2 = (
x(1,n3) -
x(1,n1))*(
x(1,n3) -xshf +
x(1,n1)) &
6471 +(
x(2,n3) -
x(2,n1))*(
x(2,n3) -yshf +
x(2,n1)) &
6472 +(
x(3,n3) -
x(3,n1))*(
x(3,n3) -zshf +
x(3,n1))
6473 h3 = (
x(1,n2) -
x(1,n1))*(
x(1,n2) -xshf +
x(1,n1)) &
6474 +(
x(2,n2) -
x(2,n1))*(
x(2,n2) -yshf +
x(2,n1)) &
6475 +(
x(3,n2) -
x(3,n1))*(
x(3,n2) -zshf +
x(3,n1))
6476 c5 = h2*(
x(3,n2) -
x(3,n1)) -h3*(
x(3,n3) -
x(3,n1))
6477 c6 = h2*(
x(2,n2) -
x(2,n1)) -h3*(
x(2,n3) -
x(2,n1))
6478 c9 = h3*(
x(1,n3) -
x(1,n1)) -h2*(
x(1,n2) -
x(1,n1))
6479 xcn = .5*xshf +(h1*rnx +(
x(2,n4) -
x(2,n1))*c5 &
6480 -(
x(3,n4) -
x(3,n1))*c6)*fac
6481 ycn = .5*yshf +(-(
x(1,n4) -
x(1,n1))*c5 +h1*rny &
6482 -(
x(3,n4) -
x(3,n1))*c9)*fac
6483 zcn = .5*zshf +((
x(1,n4) -
x(1,n1))*c6 &
6484 +(
x(2,n4) -
x(2,n1))*c9 +h1*rnz)*fac
6488 rcir =
sqrt((
x(1,n1) -xcn)**2 +(
x(2,n1) -ycn)**2 &
6495 600
FORMAT(//5
x,
'TETRAHEDRON WITH AN EXTREMELY SMALL VOLUME FOUND'// &
6496 ' IN ROUTINE CIRCUM'//5
x,
'VOLUME = ',e13.5/)
6508 xcn,ycn,zcn,vtet,rcir,iflag,tol)
6528 DOUBLE PRECISION :: xcn,xpt,ycn,ypt,zcn,zpt,vtet,rcir,tol
6529 DOUBLE PRECISION ::
x(3,*)
6532 DOUBLE PRECISION :: a1,a2,a3,b1,b2,b3,c5,c6,c9,det,fac,h1,h2,h3, &
6533 rnx,rny,rnz,vcell,xshf,yshf,zshf
6541 rnx =
cofact(
x(2,n1),
x(2,n2),
x(2,n3),
x(3,n1),
x(3,n2),
x(3,n3))
6542 rny =
cofact(
x(3,n1),
x(3,n2),
x(3,n3),
x(1,n1),
x(1,n2),
x(1,n3))
6543 rnz =
cofact(
x(1,n1),
x(1,n2),
x(1,n3),
x(2,n1),
x(2,n2),
x(2,n3))
6544 det = rnx*(xpt -
x(1,n1)) +rny*(ypt -
x(2,n1)) &
6551 IF (vtet.LT.tol) go to 30
6553 xshf = .5*(
x(1,n1) +
x(1,n2) +
x(1,n3) +xpt)
6554 yshf = .5*(
x(2,n1) +
x(2,n2) +
x(2,n3) +ypt)
6555 zshf = .5*(
x(3,n1) +
x(3,n2) +
x(3,n3) +zpt)
6556 h1 = (xpt -
x(1,n1))*(xpt -xshf +
x(1,n1)) &
6557 +(ypt -
x(2,n1))*(ypt -yshf +
x(2,n1)) &
6558 +(zpt -
x(3,n1))*(zpt -zshf +
x(3,n1))
6559 h2 = (
x(1,n3) -
x(1,n1))*(
x(1,n3) -xshf +
x(1,n1)) &
6560 +(
x(2,n3) -
x(2,n1))*(
x(2,n3) -yshf +
x(2,n1)) &
6561 +(
x(3,n3) -
x(3,n1))*(
x(3,n3) -zshf +
x(3,n1))
6562 h3 = (
x(1,n2) -
x(1,n1))*(
x(1,n2) -xshf +
x(1,n1)) &
6563 +(
x(2,n2) -
x(2,n1))*(
x(2,n2) -yshf +
x(2,n1)) &
6564 +(
x(3,n2) -
x(3,n1))*(
x(3,n2) -zshf +
x(3,n1))
6565 c5 = h2*(
x(3,n2) -
x(3,n1)) -h3*(
x(3,n3) -
x(3,n1))
6566 c6 = h2*(
x(2,n2) -
x(2,n1)) -h3*(
x(2,n3) -
x(2,n1))
6567 c9 = h3*(
x(1,n3) -
x(1,n1)) -h2*(
x(1,n2) -
x(1,n1))
6568 xcn = .5*xshf +(h1*rnx +(ypt -
x(2,n1))*c5 &
6569 -(zpt -
x(3,n1))*c6)*fac
6570 ycn = .5*yshf +(-(xpt -
x(1,n1))*c5 +h1*rny &
6571 -(zpt -
x(3,n1))*c9)*fac
6572 zcn = .5*zshf +((xpt -
x(1,n1))*c6 &
6573 +(ypt -
x(2,n1))*c9 +h1*rnz)*fac
6577 rcir =
sqrt((xpt -xcn)**2 +(ypt -ycn)**2 &
6584 600
FORMAT(//5
x,
'TETRAHEDRON WITH AN EXTREMELY SMALL VOLUME FOUND'// &
6585 ' IN ROUTINE CIRCUM2'//5
x,
'VOLUME = ',e13.5/)
6597 SUBROUTINE neighb (L1,N1,N2,N3,L2,K1,K2,NDC,NBH)
6615 INTEGER :: k1,k2,l1,l2,n1,n2,n3
6616 INTEGER :: ndc(4,*),nbh(4,*)
6618 INTEGER ::
k,kk,
m,mend,mmax,mmin,msum,m1,m2,m3,m4,nmax,nmin,nsum
6622 nmin =
min(n1,n2,n3)
6623 nmax =
max(n1,n2,n3)
6633 10 mmin =
min(m1,m2,m3)
6634 mmax =
max(m1,m2,m3)
6636 IF (mmin.EQ.nmin.AND.mmax.EQ.nmax.AND.msum.EQ.nsum) go to 30
6637 IF (m1.EQ.mend) go to 20
6651 IF (nbh(k2,l2).EQ.l1)
RETURN
6657 WRITE (6,900) n1,n2,n3, &
6658 l1,(ndc(kk,l1),kk=1,4),l2,(ndc(kk,l2),kk=1,4)
6659 900
FORMAT(
'N1,N2,N3 ',3i6,
' L1 ',i6,
' VERTS ',4i6/ &
6660 27
x,
' L2 ',i6,
' VERTS ',4i6)
6663 600
FORMAT(//5
x,
'UNABLE TO FIND NBH INDEX OF NEIGHBORING TETRAHEDRON', &
6664 5
x,
'TO ORIGINAL TETRAHEDRON. PROGRAM STOPPED IN NEIGHB')
6665 610
FORMAT(//5
x,
'UNABLE TO FIND A NEIGHBORING TETRAHEDRON THAT HAS', &
6666 5
x,
'THE REQUIRED FACE. PROGRAM STOPPED IN NEIGHB')
6679 SUBROUTINE lock (L1,L2,N1,N2,N3,N4,N5,K1,K2,NDC,NBH)
6698 INTEGER :: k1,k2,l1,l2,n1,n2,n3,n4,n5
6699 INTEGER :: ndc(4,*),nbh(4,*)
6701 INTEGER ::
m,mend,mmax,mmin,msum,m1,m2,m3,m4,
n,nend,nmax,nmin,nsum
6715 nmin =
min(n1,n2,n3)
6716 nmax =
max(n1,n2,n3)
6718 20 mmin =
min(m1,m2,m3)
6719 mmax =
max(m1,m2,m3)
6721 IF (mmin.EQ.nmin.AND.mmax.EQ.nmax.AND.msum.EQ.nsum) go to 40
6722 IF (m1.EQ.mend) go to 30
6729 30
IF (n1.EQ.nend) go to 100
6741 IF (nbh(1,l1).EQ.l2) k1 = 1
6742 IF (nbh(2,l1).EQ.l2) k1 = 2
6743 IF (nbh(3,l1).EQ.l2) k1 = 3
6744 IF (nbh(4,l1).EQ.l2) k1 = 4
6746 IF (nbh(1,l2).EQ.l1) k2 = 1
6747 IF (nbh(2,l2).EQ.l1) k2 = 2
6748 IF (nbh(3,l2).EQ.l1) k2 = 3
6749 IF (nbh(4,l2).EQ.l1) k2 = 4
6753 600
FORMAT(//5
x,
'UNABLE TO FIND THE COMMON FACE IN ROUTINE LOCK')
6775 DOUBLE PRECISION ::
facear
6778 DOUBLE PRECISION :: ax,ay,az
6779 DOUBLE PRECISION ::
x(3,*)
6783 ax = (
x(2,n2) -
x(2,n1))*(
x(3,n3) -
x(3,n1)) &
6784 -(
x(3,n2) -
x(3,n1))*(
x(2,n3) -
x(2,n1))
6785 ay = (
x(3,n2) -
x(3,n1))*(
x(1,n3) -
x(1,n1)) &
6786 -(
x(1,n2) -
x(1,n1))*(
x(3,n3) -
x(3,n1))
6787 az = (
x(1,n2) -
x(1,n1))*(
x(2,n3) -
x(2,n1)) &
6788 -(
x(2,n2) -
x(2,n1))*(
x(1,n3) -
x(1,n1))
6812 DOUBLE PRECISION ::
dihed
6815 DOUBLE PRECISION ::
x(3,*)
6817 DOUBLE PRECISION :: asiz,ax,ay,az,a1,a2,a3, &
6818 bsiz,bx,
by,bz,b1,b2,b3, &
6819 cosang,fac,
pi,px,py,pz,rad, &
6824 pi = 4.0d0*atan(1.0d0)
6830 ax =
cofact(
x(2,n1),
x(2,n2),
x(2,n3),
x(3,n1),
x(3,n2),
x(3,n3))
6831 ay =
cofact(
x(3,n1),
x(3,n2),
x(3,n3),
x(1,n1),
x(1,n2),
x(1,n3))
6832 az =
cofact(
x(1,n1),
x(1,n2),
x(1,n3),
x(2,n1),
x(2,n2),
x(2,n3))
6833 asiz =
sqrt(ax*ax +ay*ay +az*az)
6834 IF (asiz.LT.1.0
d-6)
RETURN
6842 bx =
cofact(
x(2,n1),
x(2,n4),
x(2,n2),
x(3,n1),
x(3,n4),
x(3,n2))
6843 by =
cofact(
x(3,n1),
x(3,n4),
x(3,n2),
x(1,n1),
x(1,n4),
x(1,n2))
6844 bz =
cofact(
x(1,n1),
x(1,n4),
x(1,n2),
x(2,n1),
x(2,n4),
x(2,n2))
6846 IF (bsiz.LT.1.0
d-6)
RETURN
6854 cosang = ax*bx +ay*
by +az*bz
6858 px =
x(1,n2) -
x(1,n1)
6859 py =
x(2,n2) -
x(2,n1)
6860 pz =
x(3,n2) -
x(3,n1)
6861 fac = 1.0d0/
sqrt(px*px +py*py +pz*pz)
6872 sinang = sx*px +sy*py +sz*pz
6896 DOUBLE PRECISION ::
tetar
6900 DOUBLE PRECISION ::
x(3,*)
6902 INTEGER ::
n,nend,n1,n2,n3,n4
6903 DOUBLE PRECISION :: area,a1,a2,a3,b1,b2,b3,rnx,rny,rnz
6913 10 rnx =
cofact(
x(2,n1),
x(2,n2),
x(2,n3),
x(3,n1),
x(3,n2),
x(3,n3))
6914 rny =
cofact(
x(3,n1),
x(3,n2),
x(3,n3),
x(1,n1),
x(1,n2),
x(1,n3))
6915 rnz =
cofact(
x(1,n1),
x(1,n2),
x(1,n3),
x(2,n1),
x(2,n2),
x(2,n3))
6916 area = area +
sqrt(rnx*rnx +rny*rny +rnz*rnz)
6917 IF (n1.EQ.nend) go to 20
6947 DOUBLE PRECISION ::
tetar2
6949 INTEGER :: m1,m2,m3,m4
6950 DOUBLE PRECISION ::
x(3,*)
6952 INTEGER ::
n,nend,n1,n2,n3,n4
6953 DOUBLE PRECISION :: area,a1,a2,a3,b1,b2,b3,rnx,rny,rnz
6963 10 rnx =
cofact(
x(2,n1),
x(2,n2),
x(2,n3),
x(3,n1),
x(3,n2),
x(3,n3))
6964 rny =
cofact(
x(3,n1),
x(3,n2),
x(3,n3),
x(1,n1),
x(1,n2),
x(1,n3))
6965 rnz =
cofact(
x(1,n1),
x(1,n2),
x(1,n3),
x(2,n1),
x(2,n2),
x(2,n3))
6966 area = area +
sqrt(rnx*rnx +rny*rny +rnz*rnz)
6967 IF (n1.EQ.nend) go to 20
6998 DOUBLE PRECISION ::
tetar3
7001 DOUBLE PRECISION :: xpt,ypt,zpt
7002 DOUBLE PRECISION ::
x(3,*)
7004 INTEGER ::
n,nend,n1,n2,n3
7005 DOUBLE PRECISION :: area,a1,a2,a3,b1,b2,b3,rnx,rny,rnz
7014 10 rnx =
cofact(
x(2,n1),
x(2,n2),ypt,
x(3,n1),
x(3,n2),zpt)
7015 rny =
cofact(
x(3,n1),
x(3,n2),zpt,
x(1,n1),
x(1,n2),xpt)
7016 rnz =
cofact(
x(1,n1),
x(1,n2),xpt,
x(2,n1),
x(2,n2),ypt)
7017 area = area +
sqrt(rnx*rnx +rny*rny +rnz*rnz)
7018 IF (n1.EQ.nend) go to 20
7024 20 rnx =
cofact(
x(2,n1),
x(2,n2),
x(2,n3),
x(3,n1),
x(3,n2),
x(3,n3))
7025 rny =
cofact(
x(3,n1),
x(3,n2),
x(3,n3),
x(1,n1),
x(1,n2),
x(1,n3))
7026 rnz =
cofact(
x(1,n1),
x(1,n2),
x(1,n3),
x(2,n1),
x(2,n2),
x(2,n3))
7027 tetar3 = area +
sqrt(rnx*rnx +rny*rny +rnz*rnz)
7039 SUBROUTINE fangle (J,X,NFCE,ANGL1,ANGL2,ANGL3,Q)
7055 INTEGER :: nfce(3,*)
7056 DOUBLE PRECISION :: angl1,angl2,angl3,
q
7057 DOUBLE PRECISION ::
x(3,*)
7060 DOUBLE PRECISION :: chalf,c1,c2,c3,da,db,dc,prod1,prod2,prod3, &
7061 rad,s1,s2,s3,t1,t2,t3
7065 rad = 45.0d0/atan(1.0d0)
7069 prod1 = (
x(1,n2) -
x(1,n1))*(
x(1,n3) -
x(1,n1)) &
7070 +(
x(2,n2) -
x(2,n1))*(
x(2,n3) -
x(2,n1)) &
7071 +(
x(3,n2) -
x(3,n1))*(
x(3,n3) -
x(3,n1))
7072 da =
sqrt((
x(1,n2) -
x(1,n1))**2 +(
x(2,n2) -
x(2,n1))**2 &
7073 +(
x(3,n2) -
x(3,n1))**2)
7074 db =
sqrt((
x(1,n3) -
x(1,n1))**2 +(
x(2,n3) -
x(2,n1))**2 &
7075 +(
x(3,n3) -
x(3,n1))**2)
7078 IF (chalf.LT.1.e-6) go to 200
7079 t1 =
sqrt(1./chalf -1)
7083 prod2 = (
x(1,n3) -
x(1,n2))*(
x(1,n1) -
x(1,n2)) &
7084 +(
x(2,n3) -
x(2,n2))*(
x(2,n1) -
x(2,n2)) &
7085 +(
x(3,n3) -
x(3,n2))*(
x(3,n1) -
x(3,n2))
7086 da =
sqrt((
x(1,n3) -
x(1,n2))**2 +(
x(2,n3) -
x(2,n2))**2 &
7087 +(
x(3,n3) -
x(3,n2))**2)
7088 db =
sqrt((
x(1,n1) -
x(1,n2))**2 +(
x(2,n1) -
x(2,n2))**2 &
7089 +(
x(3,n1) -
x(3,n2))**2)
7092 IF (chalf.LT.1.e-6) go to 200
7093 t2 =
sqrt(1./chalf -1)
7097 prod3 = (
x(1,n1) -
x(1,n3))*(
x(1,n2) -
x(1,n3)) &
7098 +(
x(2,n1) -
x(2,n3))*(
x(2,n2) -
x(2,n3)) &
7099 +(
x(3,n1) -
x(3,n3))*(
x(3,n2) -
x(3,n3))
7100 da =
sqrt((
x(1,n1) -
x(1,n3))**2 +(
x(2,n1) -
x(2,n3))**2 &
7101 +(
x(3,n1) -
x(3,n3))**2)
7102 db =
sqrt((
x(1,n2) -
x(1,n3))**2 +(
x(2,n2) -
x(2,n3))**2 &
7103 +(
x(3,n2) -
x(3,n3))**2)
7106 IF (chalf.LT.1.e-6) go to 200
7107 t3 =
sqrt(1./chalf -1)
7111 q = .5*(s1 +s2 +s3)/(s1*s2*s3)
7113 200
WRITE (6,600)
j,n1,
x(1,n1),
x(2,n1),
x(3,n1), &
7114 n2,
x(1,n2),
x(2,n2),
x(3,n2), &
7115 n3,
x(1,n3),
x(2,n3),
x(3,n3)
7117 600
FORMAT(//5
x,
'FACE ',i6,
' HAS AN ANGLE OF 180 DEGREES'// &
7118 5
x,
'VERTEX ',i6,
' X = ',f10.4,
' Y = ',f10.4,
' Z = ',f10.4/ &
7119 5
x,
'VERTEX ',i6,
' X = ',f10.4,
' Y = ',f10.4,
' Z = ',f10.4/ &
7120 5
x,
'VERTEX ',i6,
' X = ',f10.4,
' Y = ',f10.4,
' Z = ',f10.4// &
7121 5
x,
'PROGRAM STOPPED IN FANGLE')
7150 INTEGER :: ndc(4,*),nbh(4,*),iprot(*)
7151 DOUBLE PRECISION ::
x(3,*)
7153 INTEGER ::
j,npass,n1,n2,n3,n4
7154 DOUBLE PRECISION :: angmax,angmin,ang1,ang2,ang3,ang4,ang5,ang6, &
7155 dsmax,dsmin,ds1,ds2,ds3,ds4,ds5,ds6
7164 IF (iprot(
j).EQ.1) go to 20
7165 IF (ndc(1,
j).EQ.0) go to 20
7166 IF (nbh(1,
j).EQ.0) go to 20
7171 IF (n4.EQ.-1) go to 20
7172 ang1 =
dihed(n1,n2,n3,n4,
x)
7173 IF (ang1.GT.180.) ang1 = 360. -ang1
7174 ang2 =
dihed(n2,n3,n1,n4,
x)
7175 IF (ang2.GT.180.) ang2 = 360. -ang2
7176 ang3 =
dihed(n3,n1,n2,n4,
x)
7177 IF (ang3.GT.180.) ang3 = 360. -ang3
7178 ang4 =
dihed(n1,n4,n2,n3,
x)
7179 IF (ang4.GT.180.) ang4 = 360. -ang4
7180 ang5 =
dihed(n2,n4,n1,n3,
x)
7181 IF (ang5.GT.180.) ang5 = 360. -ang5
7182 ang6 =
dihed(n3,n4,n1,n2,
x)
7183 IF (ang6.GT.180.) ang6 = 360. -ang6
7184 angmin =
min(angmin,ang1,ang2,ang3,ang4,ang5,ang6)
7185 angmax =
max(angmax,ang1,ang2,ang3,ang4,ang5,ang6)
7186 ds1 = (
x(1,n1) -
x(1,n2))**2 +(
x(2,n1) -
x(2,n2))**2 &
7187 +(
x(3,n1) -
x(3,n2))**2
7188 ds2 = (
x(1,n1) -
x(1,n3))**2 +(
x(2,n1) -
x(2,n3))**2 &
7189 +(
x(3,n1) -
x(3,n3))**2
7190 ds3 = (
x(1,n1) -
x(1,n4))**2 +(
x(2,n1) -
x(2,n4))**2 &
7191 +(
x(3,n1) -
x(3,n4))**2
7192 ds4 = (
x(1,n2) -
x(1,n3))**2 +(
x(2,n2) -
x(2,n3))**2 &
7193 +(
x(3,n2) -
x(3,n3))**2
7194 ds5 = (
x(1,n2) -
x(1,n4))**2 +(
x(2,n2) -
x(2,n4))**2 &
7195 +(
x(3,n2) -
x(3,n4))**2
7196 ds6 = (
x(1,n3) -
x(1,n4))**2 +(
x(2,n3) -
x(2,n4))**2 &
7197 +(
x(3,n3) -
x(3,n4))**2
7198 IF (npass.EQ.1) go to 10
7200 dsmin =
min(ds1,ds2,ds3,ds4,ds5,ds6)
7201 dsmax =
max(ds1,ds2,ds3,ds4,ds5,ds6)
7203 10 dsmin =
min(dsmin,ds1,ds2,ds3,ds4,ds5,ds6)
7204 dsmax =
max(dsmax,ds1,ds2,ds3,ds4,ds5,ds6)
7208 WRITE (6,600) angmin,angmax,dsmin,dsmax
7210 600
FORMAT( 5
x,
'*************************************************'/ &
7212 5
x,
'* MINIMUM DIHEDRAL ANGLE IS ',f6.2,
' DEGREES *'/ &
7213 5
x,
'* MAXIMUM DIHEDRAL ANGLE IS ',f6.2,
' DEGREES *'/ &
7215 5
x,
'* MINIMUM EDGE DISTANCE IS ',f12.4,
' *'/ &
7216 5
x,
'* MAXIMUM EDGE DISTANCE IS ',f12.4,
' *'/ &
7218 5
x,
'*************************************************')
7229 SUBROUTINE tetcof (N1,N2,N3,N4,XPT,YPT,ZPT,V1,V2,V3,V4,VTET,X)
7246 DOUBLE PRECISION :: vtet,v1,v2,v3,v4,xpt,ypt,zpt
7247 DOUBLE PRECISION ::
x(3,*)
7249 INTEGER :: n1,n2,n3,n4
7250 DOUBLE PRECISION :: a1,a2,a3,b1,b2,b3,rnx,rny,rnz
7254 rnx =
cofact(
x(2,n1),
x(2,n2),
x(2,n3),
x(3,n1),
x(3,n2),
x(3,n3))
7255 rny =
cofact(
x(3,n1),
x(3,n2),
x(3,n3),
x(1,n1),
x(1,n2),
x(1,n3))
7256 rnz =
cofact(
x(1,n1),
x(1,n2),
x(1,n3),
x(2,n1),
x(2,n2),
x(2,n3))
7257 v4 = abs(rnx*(xpt -
x(1,n1)) +rny*(ypt -
x(2,n1)) &
7258 +rnz*(zpt -
x(3,n1)))
7259 vtet = abs(rnx*(
x(1,n4) -
x(1,n1)) +rny*(
x(2,n4) -
x(2,n1)) &
7260 +rnz*(
x(3,n4) -
x(3,n1)))
7261 rnx =
cofact(
x(2,n2),
x(2,n3),
x(2,n4),
x(3,n2),
x(3,n3),
x(3,n4))
7262 rny =
cofact(
x(3,n2),
x(3,n3),
x(3,n4),
x(1,n2),
x(1,n3),
x(1,n4))
7263 rnz =
cofact(
x(1,n2),
x(1,n3),
x(1,n4),
x(2,n2),
x(2,n3),
x(2,n4))
7264 v1 = abs(rnx*(xpt -
x(1,n2)) +rny*(ypt -
x(2,n2)) &
7265 +rnz*(zpt -
x(3,n2)))
7266 rnx =
cofact(
x(2,n3),
x(2,n4),
x(2,n1),
x(3,n3),
x(3,n4),
x(3,n1))
7267 rny =
cofact(
x(3,n3),
x(3,n4),
x(3,n1),
x(1,n3),
x(1,n4),
x(1,n1))
7268 rnz =
cofact(
x(1,n3),
x(1,n4),
x(1,n1),
x(2,n3),
x(2,n4),
x(2,n1))
7269 v2 = abs(rnx*(xpt -
x(1,n3)) +rny*(ypt -
x(2,n3)) &
7270 +rnz*(zpt -
x(3,n3)))
7271 rnx =
cofact(
x(2,n4),
x(2,n1),
x(2,n2),
x(3,n4),
x(3,n1),
x(3,n2))
7272 rny =
cofact(
x(3,n4),
x(3,n1),
x(3,n2),
x(1,n4),
x(1,n1),
x(1,n2))
7273 rnz =
cofact(
x(1,n4),
x(1,n1),
x(1,n2),
x(2,n4),
x(2,n1),
x(2,n2))
7274 v3 = abs(rnx*(xpt -
x(1,n4)) +rny*(ypt -
x(2,n4)) &
7275 +rnz*(zpt -
x(3,n4)))
7305 DOUBLE PRECISION ::
x(3,*)
7306 INTEGER :: l,
n,ncnt,ncont,nend,npt,n1,n2,n3,n4
7309 DOUBLE PRECISION :: a1,a2,a3,b1,b2,b3,ptot,rnx,rny,rnz,sc,sp,vdiff
7310 DOUBLE PRECISION ::
v(4),vcell(4),ptest(4)
7324 rnx =
cofact(
x(2,n1),
x(2,n2),
x(2,n3),
x(3,n1),
x(3,n2),
x(3,n3))
7325 rny =
cofact(
x(3,n1),
x(3,n2),
x(3,n3),
x(1,n1),
x(1,n2),
x(1,n3))
7326 rnz =
cofact(
x(1,n1),
x(1,n2),
x(1,n3),
x(2,n1),
x(2,n2),
x(2,n3))
7327 v(ncnt) = rnx*(
x(1,npt) -
x(1,n1)) +rny*(
x(2,npt) -
x(2,n1)) &
7328 +rnz*(
x(3,npt) -
x(3,n1))
7329 vcell(ncnt) = rnx*(
x(1,n4) -
x(1,n1)) +rny*(
x(2,n4) -
x(2,n1)) &
7330 +rnz*(
x(3,n4) -
x(3,n1))
7331 sp =
sign(1.0d0,
v(ncnt))
7332 sc =
sign(1.0d0,vcell(ncnt))
7333 v(ncnt) = sp*
v(ncnt)
7334 vcell(ncnt) = sc*vcell(ncnt)
7336 IF (n1.EQ.nend) go to 20
7347 20 vdiff =
v(1)+
v(2)+
v(3)+
v(4) &
7348 -0.25d0*(vcell(1)+vcell(2)+vcell(3)+vcell(4))
7349 ptot = ptest(1)+ptest(2)+ptest(3)+ptest(4)
7351 IF (ptot.GT.3.0d0)
RETURN
7353 IF (ptot.GT.-3.0d0)
RETURN
7356 WRITE (6,900) n1,
x(1,n1),
x(2,n1),
x(3,n1), &
7357 n2,
x(1,n2),
x(2,n2),
x(3,n2), &
7358 n3,
x(1,n3),
x(2,n3),
x(3,n3), &
7359 n4,
x(1,n4),
x(2,n4),
x(3,n4), &
7360 npt,
x(1,npt),
x(2,npt),
x(3,npt), &
7361 v(1),
v(2),
v(3),
v(4), &
7362 vcell(1),vcell(2),vcell(3),vcell(4), &
7363 ptest(1),ptest(2),ptest(3),ptest(4)
7364 900
FORMAT(
'N1 ',i5,
' X = ',f10.4,
' Y = ',f10.4,
'Z = ',f10.4/ &
7365 'N2 ',i5,
' X = ',f10.4,
' Y = ',f10.4,
'Z = ',f10.4/ &
7366 'N3 ',i5,
' X = ',f10.4,
' Y = ',f10.4,
'Z = ',f10.4/ &
7367 'N4 ',i5,
' X = ',f10.4,
' Y = ',f10.4,
'Z = ',f10.4/ &
7368 'NPT ',i5,
' X = ',f10.4,
' Y = ',f10.4,
'Z = ',f10.4/ &
7369 'V1 = ',e13.5,
' V2 = ',e13.5,
' V3 = ',e13.5,
' V4 = ',e13.5/ &
7370 'VCELL = ',4(2
x,e13.5)/ &
7371 'P1 = ',e13.5,
' P2 = ',e13.5,
' P3 = ',e13.5,
' P4 = ',e13.5)
7385 DOUBLE PRECISION ::
cofact
7387 DOUBLE PRECISION,
INTENT(IN) :: a1,a2,a3,b1,b2,b3
7389 cofact = (a2 -a1)*(b3 -b1) -(a3 -a1)*(b2 -b1)
subroutine tessel(IRING, NRING, NFAD, NBHKP, IEDKP, NEDK, IMEET, ITOUCH, X, IPOINT, NPOINT, DX, DY, DZ, DS, NP)
*********************************************************************Illinois Open Source License ****University of Illinois NCSA **Open Source License University of Illinois All rights reserved ****Developed by
double precision function dihed(N1, N2, N3, N4, X)
**********************************************************************Rocstar Simulation Suite Illinois Rocstar LLC All rights reserved ****Illinois Rocstar LLC IL **www illinoisrocstar com **sales illinoisrocstar com WITHOUT WARRANTY OF ANY **EXPRESS OR INCLUDING BUT NOT LIMITED TO THE WARRANTIES **OF FITNESS FOR A PARTICULAR PURPOSE AND **NONINFRINGEMENT IN NO EVENT SHALL THE CONTRIBUTORS OR **COPYRIGHT HOLDERS BE LIABLE FOR ANY DAMAGES OR OTHER WHETHER IN AN ACTION OF TORT OR **Arising OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE **USE OR OTHER DEALINGS WITH THE SOFTWARE **********************************************************************INTERFACE SUBROUTINE kc
Tfloat sum() const
Return the sum of all the pixel values in an image.
double precision function tetar2(M1, M2, M3, M4, X)
subroutine rflu_repair3d(NBPTS, NBFACE, NNODE, NCELL, XI, NFCEI, NDCI, XBNDYI, modFlag)
static SURF_BEGIN_NAMESPACE double sign(double x)
subroutine densfn(X, NNODE, NFCE, NBFACE, NEDGE, NDG, DENS, RESID, FVCNT)
subroutine volput(X, ITYP, NBPTS, NNODE, NDC, NBH, IPROT, NCELL, NDG, IDGP, NDGP, NEDGE, VOL, XCEN, YCEN, ZCEN, RC, RAT, DENS, NPTET, NACPT, IDONE, NREF, NLINK, NOCTR, IOCTR, XFAR, YFAR, ZFAR, XOCTR, YOCTR, ZOCTR, XHOLD, YHOLD, ZHOLD, XKEEP, YKEEP, ZKEEP, KSRCH, NSRCH, IPOINT, NPOINT, IFLAG, NFLAG, NFILL, NEWCEL, NTRI, NCAV, NSHAKE, NEWNBH, NOLD, NCAVFC, IKEEP, LDEL, NEDGRM, XC, YC, ZC, V, RAD, RCRIN, LNKUP, LNKDN, JLAST, JFIRST, NTRACK, VOLMIN, RCMX, TOLV)
Vector_n max(const Array_n_const &v1, const Array_n_const &v2)
**********************************************************************Rocstar Simulation Suite Illinois Rocstar LLC All rights reserved ****Illinois Rocstar LLC IL **www illinoisrocstar com **sales illinoisrocstar com WITHOUT WARRANTY OF ANY **EXPRESS OR INCLUDING BUT NOT LIMITED TO THE WARRANTIES **OF FITNESS FOR A PARTICULAR PURPOSE AND **NONINFRINGEMENT IN NO EVENT SHALL THE CONTRIBUTORS OR **COPYRIGHT HOLDERS BE LIABLE FOR ANY DAMAGES OR OTHER WHETHER IN AN ACTION OF TORT OR **Arising OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE **USE OR OTHER DEALINGS WITH THE SOFTWARE **********************************************************************INTERFACE SUBROUTINE ic
subroutine sliver(IRING, NRING, NA, NB, NFAIL, VORIG, X, DENS, ITYP, IPOINT, NNODE, NOCTR, IOCTR, NLINK, NREF, XFAR, YFAR, ZFAR, VLT)
subroutine tetmv(X, XNEWBN, NNODE, NDC, NBH, NCELL, NFCE, NBFACE, NFAIL, ITYP, XCEN, YCEN, ZCEN, VOL, RC, RAT, NEDGE, NDG, IDGP, NDGP, IPOINT, FCOUNT, XC, RESID, SV, SIG1, SIG2, SIG3, XFAR, YFAR, ZFAR, TOLV)
subroutine colaps(NA, NB, NCOL, NVERT, NFAIL, X, NDC, NBH, IPROT, ITYP, XCEN, YCEN, ZCEN, VOL, RC, RAT, NVCNT, IFLAG, NFLAG, NPTET, NDG, IDGP, NDGP, NOCTR, IOCTR, NLINK, XFAR, YFAR, ZFAR, IDONE, NREF, KSRCH, NSRCH, IRING, NTETKP, LNBR, ISHK, MNBR, KSHK, TOLV)
double precision function facear(X, N1, N2, N3)
subroutine smooth(X, ITYP, NNODE, NDC, NBH, IPROT, NCELL, NDG, IDGP, NDGP, NEDGE, VOL, XCEN, YCEN, ZCEN, RC, RAT, DENS, NPTET, NACPT, IDONE, NREF, NLINK, NOCTR, IOCTR, XFAR, YFAR, ZFAR, XOCTR, YOCTR, ZOCTR, XHOLD, YHOLD, ZHOLD, XKEEP, YKEEP, ZKEEP, KSRCH, NSRCH, IPOINT, NPOINT, IFLAG, NFLAG, DX, DY, DZ, DS, VLT, IRING, NTETKP, NFAD, NEWC, NBHKP, IEDKP, LNBR, ISHK, MNBR, KSHK, NPP, NFILL, NEWCEL, NTRI, NCAV, NSHAKE, NEWNBH, NOLD, NCAVFC, IKEEP, LDEL, NEDGRM, XC, YC, ZC, V, RAD, RCRIN, LNKUP, LNKDN, LISTF, VOLMIN, RCMX, TOLV)
subroutine angfnd(NUM, IRING, IMIN, IMEET, ITOUCH, X, DX, DY, DZ, DS)
*********************************************************************Illinois Open Source License ****University of Illinois NCSA **Open Source License University of Illinois All rights reserved ****Developed free of to any person **obtaining a copy of this software and associated documentation to deal with the Software without including without limitation the rights to and or **sell copies of the and to permit persons to whom the **Software is furnished to do subject to the following this list of conditions and the following disclaimers ****Redistributions in binary form must reproduce the above **copyright this list of conditions and the following **disclaimers in the documentation and or other materials **provided with the distribution ****Neither the names of the Center for Simulation of Advanced the University of nor the names of its **contributors may be used to endorse or promote products derived **from this Software without specific prior written permission ****THE SOFTWARE IS PROVIDED AS WITHOUT WARRANTY OF ANY **EXPRESS OR INCLUDING BUT NOT LIMITED TO THE WARRANTIES **OF FITNESS FOR A PARTICULAR PURPOSE AND **NONINFRINGEMENT IN NO EVENT SHALL THE CONTRIBUTORS OR **COPYRIGHT HOLDERS BE LIABLE FOR ANY DAMAGES OR OTHER WHETHER IN AN ACTION OF TORT OR **ARISING OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE **USE OR OTHER DEALINGS WITH THE SOFTWARE v
subroutine cavbnd(NP, LCONT, NDEL, LDEL, NCNT, X, NDC, NBH, VOL, IFLAG, NFLAG, NFILL, NEWCEL, NTRI, NEWNBH, NOLD, TOLV)
subroutine fangle(J, X, NFCE, ANGL1, ANGL2, ANGL3, Q)
subroutine octfnd(NCODE, NCLOSE, XPT, YPT, ZPT, N1, N2, N3, DISMIN, X, NOCTR, NLINK, XFAR, YFAR, ZFAR, IDONE, XOCTR, YOCTR, ZOCTR, XHOLD, YHOLD, ZHOLD, XKEEP, YKEEP, ZKEEP, KSRCH, NSRCH)
subroutine tetang(X, NDC, NBH, IPROT, NCELL)
subroutine tetcof(N1, N2, N3, N4, XPT, YPT, ZPT, V1, V2, V3, V4, VTET, X)
subroutine snglar(NNODE, X, XD, NCELL, NDC, NBH, SIGMN, SIGMX, SIG1, SIG2, SIG3, CNDMIN, CNDMAX)
subroutine cavity(NP, LCONT, NDEL, LDEL, X, NDC, NBH, IPROT, IFLAG, NFLAG, XCEN, YCEN, ZCEN, RC, NFILL, NEWCEL, TOLV)
subroutine radrat(X, NNODE, NDC, NBH, IPROT, NCELL, NFCE, NBFACE, ITYP, IPOINT, VOL, RC, RAT)
subroutine input(X, NNODE, NDC, NCELL, NFCE, NBPTS, NBFACE, ITYP, NPROP, XBNDY, XFAR, YFAR, ZFAR)
subroutine putpnt(J, NFAIL, X, ITYP, NNODE, NDC, NBH, IPROT, NCELL, NDG, IDGP, NDGP, NEDGE, VOL, XCEN, YCEN, ZCEN, RC, RAT, DENS, NPTET, NACPT, IDONE, NREF, NLINK, NOCTR, IOCTR, XFAR, YFAR, ZFAR, XOCTR, YOCTR, ZOCTR, XHOLD, YHOLD, ZHOLD, XKEEP, YKEEP, ZKEEP, KSRCH, NSRCH, IPOINT, NPOINT, IFLAG, NFLAG, NFILL, NEWCEL, NTRI, NCAV, NSHAKE, NEWNBH, NOLD, NCAVFC, IKEEP, LDEL, NEDGRM, XC, YC, ZC, V, RAD, RCRIN, LNKUP, LNKDN, VOLMIN, RCMX, TOLV)
subroutine edglen(X, NNODE, ITYP, NEDGE, NDG, DENS, NVCNT)
subroutine tetmod(X, ITYP, NNODE, NDC, NBH, IPROT, NCELL, NDG, IDGP, NDGP, NEDGE, NFCE, NBFACE, VOL, XCEN, YCEN, ZCEN, RC, RAT, DENS, NPTET, NACPT, SIG1, SIG2, SIG3, NVCNT, RESID, COUNT, FAC, IDONE, NREF, NLINK, NOCTR, IOCTR, XFAR, YFAR, ZFAR, XOCTR, YOCTR, ZOCTR, XHOLD, YHOLD, ZHOLD, XKEEP, YKEEP, ZKEEP, KSRCH, NSRCH, IPOINT, NPOINT, IFLAG, NFLAG, DX, DY, DZ, DS, VLT, IRING, NTETKP, NFAD, NEWC, NBHKP, IEDKP, LNBR, ISHK, MNBR, KSHK, NPP, NFILL, NEWCEL, NTRI, NCAV, NSHAKE, NEWNBH, NOLD, NCAVFC, IKEEP, LDEL, NEDGRM, XC, YC, ZC, V, RAD, RCRIN, LNKUP, LNKDN, LISTF, VOLMIN, RCMX, TOLV)
subroutine neighb(L1, N1, N2, N3, L2, K1, K2, NDC, NBH)
subroutine circum(X, N1, N2, N3, N4, XCN, YCN, ZCN, VTET, RCIR, IFLAG, TOL)
subroutine lock(L1, L2, N1, N2, N3, N4, N5, K1, K2, NDC, NBH)
subroutine cavedg(NCNT, NEDG, IPOINT, NPOINT, NTRI, NCAV, NEWNBH, NEWCEL)
static void div(const Attribute *x, const Attribute *y, Attribute *z)
Operation wrapper for division.
subroutine octrmv(N, X, NOCTR, NLINK, IDONE, NREF)
subroutine struct(X, NNODE, NDC, NBH, IPROT, NCELL, NFCE, NBFACE, NEDGE, NDG, IDGP, NDGP, IPOINT, NPOINT, NPTET, XCEN, YCEN, ZCEN, VOL, RC, RAT, XFAR, YFAR, ZFAR, IOCTR, NLINK, NOCTR, IDONE, NREF, VOLMIN, RCMX)
subroutine tetloc(NP, NCLOSE, LBRK, LCONT, NFAIL, X, NDC, NBH, IPROT, DENS, VOL, IFLAG, NFLAG, NFILL, NEWCEL, TOLV)
subroutine tree(PROP, NLEFT, NRIGHT, NBACK, LISTF, NACPT, NTOT, NBH, IPROT, NCELL)
subroutine volcom(L, NPT, VDIFF, NCONT, X, NDC)
Vector_n min(const Array_n_const &v1, const Array_n_const &v2)
double precision function cofact(A1, A2, A3, B1, B2, B3)
subroutine replace(N1, N2, NA, NB, J, NRING, NPASS, NFAIL, X, DENS, ITYP, NNODE, NDC, NBH, IPROT, NCELL, NDG, IDGP, NDGP, NEDGE, IPOINT, NPOINT, VOL, XCEN, YCEN, ZCEN, RC, RAT, NPTET, NACPT, NOCTR, IOCTR, NLINK, NREF, XFAR, YFAR, ZFAR, IRING, NTETKP, LNBR, ISHK, MNBR, KSHK, NFAD, NBHKP, IEDKP, NEWC, DX, DY, DZ, DS, NP, VLT, TOLV)
double precision function tetar(L, X, NDC)
subroutine circum2(X, N1, N2, N3, XPT, YPT, ZPT, XCN, YCN, ZCN, VTET, RCIR, IFLAG, TOL)
subroutine octfil(N, X, NOCTR, IOCTR, NLINK, NREF, XFAR, YFAR, ZFAR)
void next()
Go to the next element within the connectivity tables of a pane.
long double dist(long double *coord1, long double *coord2, int size)
subroutine datsrf(NP, NCNT, NEDG, NDG, IDGP, NDGP, NEDGE, IPOINT, NCAV, NEDGRM, IEDGRM)
void insert(Attribute *attr)
Insert an attribute onto the pane.
subroutine coarsn(LC, X, NNODE, NDC, NBH, IPROT, NCELL, ITYP, XCEN, YCEN, ZCEN, VOL, RC, RAT, NVCNT, DENS, IFLAG, NFLAG, NPTET, NEDGE, NDG, IDGP, NDGP, NOCTR, IOCTR, NLINK, XFAR, YFAR, ZFAR, IDONE, NREF, FAC, SIG1, SIG2, SIG3, KSRCH, NSRCH, IRING, NTETKP, LNBR, ISHK, MNBR, KSHK, TOLV)
subroutine triswp(NPT, N1, N2, N3, J, NFAIL, X, NDC, NBH, IPROT, NCELL, NDG, IDGP, NDGP, NEDGE, VOL, XCEN, YCEN, ZCEN, RC, RAT, NPTET, NACPT, TOLV)
CImg< T > & atan2(const CImg< t > &img)
Compute the arc-tangent of each pixel.
subroutine recon(LDEL, NDEL, NCNT, NDC, NBH, IPROT, NDG, IDGP, NDGP, NFLAG, NTRI, NCAVFC, IKEEP, NEDGRM, IEDGRM)
subroutine edgerm(NEDG, NDG, IDGP, NDGP)
double precision function tetar3(M1, M2, M3, XPT, YPT, ZPT, X)