SUBROUTINE SUPERT(XPT,YPT,NVERT) USE BLKMAP ! INCLUDE 'BLK1.COM' REAL*8 XPT(*),YPT(*) REAL*8 XMINM,YMINM,X45 DATA VDX9/-9.E9/ ! Find minimum x and y xminm=1.e20 yminm=1.e20 x45=-1.e20 DO J=1,NVERT IF(XPT(J) .GT. VDX9) THEN if(xminm .GT. XPT(j) ) then xminm=XPT(j) end if IF(yminm .GT. YPT(j)) then yminm=YPT(j) endif ENDIF ENDDO ! Find max at 45 degrees DO J=1,NVERT IF(XPT(J) .GT. VDX9) THEN X45T=((XPT(J)-XMINM)+(YPT(J)-YMINM))/1.414 IF(x45 .LT. X45T) THEN X45=X45T ENDIF ENDIF END DO XPT(NVERT+1)=XMINM-5 YPT(NVERT+1)=YMINM-5. XPT(NVERT+2)=XMINM+1.414*X45+10. YPT(NVERT+2)=YMINM-5. XPT(NVERT+3)=XMINM-5. YPT(NVERT+3)=YMINM+1.414*X45+10. NELT=1 NOPEL(1,1)=NVERT+1 NOPEL(1,2)=NVERT+2 NOPEL(1,3)=NVERT+3 NVERT=NVERT+3 CALL CCENTRE(XPT(NOPEL(1,1)),XPT(NOPEL(1,2)),XPT(NOPEL(1,3)) & &,YPT(NOPEL(1,1)),YPT(NOPEL(1,2)),YPT(NOPEL(1,3)) & &,XCEN(1),YCEN(1),RADS(1)) RETURN END SUBROUTINE SUBROUTINE INSIDCIRC(XPT,YPT,J,N,ISWT) ! Test for point inside circumcircle USE BLKMAP ! INCLUDE 'BLK1.COM' REAL*8 XPT(*),YPT(*) REAL*8 DISQ ! Get the distance for this element DISQ=(XCEN(J)-XPT(N))**2+(YCEN(J)-YPT(N))**2 ! Test against the radius IF(DISQ .GT. RADS(J)*RADS(J)) THEN ISWT=0 ELSE ISWT=1 ENDIF RETURN END SUBROUTINE SUBROUTINE PROCESS(J,NEDGE,NGAP) ! Drop triangle and form edge buffers USE BLKMAP USE BLK1MOD ! INCLUDE 'BLK1.COM' NEDGE=NEDGE+3 IEDGE(NEDGE-2,1)=NOPEL(J,1) IEDGE(NEDGE-1,1)=NOPEL(J,2) IEDGE(NEDGE,1) =NOPEL(J,3) IEDGE(NEDGE-2,2)=NOPEL(J,2) IEDGE(NEDGE-1,2)=NOPEL(J,3) IEDGE(NEDGE,2) =NOPEL(J,1) NOPEL(J,1)=0 NOPEL(J,2)=0 NOPEL(J,3)=0 NGAP=NGAP+1 IGAP(NGAP)=J RETURN END SUBROUTINE SUBROUTINE FORMT(XPT,YPT,J,N,NGAP,K,WD) ! Form the triangle USE BLKMAP REAL*8 XPT(*),YPT(*) REAL WD(*) ! INCLUDE 'BLK1.COM' IF(NGAP .GT. 0) THEN K=IGAP(NGAP) NGAP=NGAP-1 ELSE NELTS=NELTS+1 K=NELTS ENDIF NOPEL(K,1)=IEDGE(J,1) NOPEL(K,2)=IEDGE(J,2) NOPEL(K,3)=N CALL TESTANG(XPT,YPT,K,WD) ! Now get circumcircle data CALL CCENTRE(XPT(NOPEL(K,1)),XPT(NOPEL(K,2)),XPT(NOPEL(K,3)) & &,YPT(NOPEL(K,1)),YPT(NOPEL(K,2)),YPT(NOPEL(K,3)) & &,XCEN(K),YCEN(K),RADS(K)) RETURN END SUBROUTINE SUBROUTINE CCENTRE(X1,X2,X3,Y1,Y2,Y3,XC,YC,RC) ! get circumcentre and radius REAL*8 X1,Y1,X2,Y2,X3,Y3,A,B,C,D,AF,R1,R2,RC,XC,YC A=X2-X1 B=Y2-Y1 C=X3-X1 D=Y3-Y1 AF=2.*(B*C-A*D) R1=(-D*(A**2+B**2) + B*(C**2+D**2))/AF R2=( C*(A**2+B**2) - A*(C**2+D**2))/AF RC=SQRT(R1**2+R2**2) XC=X1+R1 YC=Y1+R2 RETURN END SUBROUTINE SUBROUTINE RIDPOINT(NVERT) USE BLKMAP NCOUNT=0 DO N=1,NELTS DO K=1,3 IF(NOPEL(N,K) .GT. NVERT-3) THEN DO L=1,3 NOPEL(N,L)=0 ENDDO GO TO 500 ENDIF ENDDO NCOUNT=NCOUNT+1 DO K=1,3 NOPEL(NCOUNT,K)=NOPEL(N,K) ENDDO XCEN(NCOUNT)=XCEN(N) YCEN(NCOUNT)=YCEN(N) RADS(NCOUNT)=RADS(N) 500 CONTINUE ENDDO NELTS=NCOUNT RETURN END SUBROUTINE SORTDB(A,NKEY,N) !*********************************** .....SORT..... !- !......SORT IS A SIMPLE SHELL SORT ROUTINE IN DOUBLE PRECISION !- ! SHELL SORT SAVE ! !IPK JAN94 INTEGER*2 NKEY REAL*8 A(*) INTEGER NKEY(*) IF(N.LT.2) RETURN DO 90 J=1,N NKEY(J)=J 90 END DO ID = N 100 ID = ID / 2 110 IB = 1 120 GO TO 200 130 IB = IB + 1 IF( IB .LE. ID ) GO TO 200 IF( ID .GT. 1 ) GO TO 100 RETURN 200 I = IB 210 K = I + ID 220 IF( A(NKEY(I)) .LE. A(NKEY(K)) ) GO TO 250 NKT = NKEY(K) NKEY(K) = NKEY(I) J = I 230 K = J - ID IF( K .LT. 1 ) GO TO 240 IF( A(NKT) .GT. A(NKEY(K)) ) GO TO 240 NKEY(J) = NKEY(K) J = K GO TO 230 240 NKEY(J) = NKT 250 I = I + ID IF( I + ID .LE. N ) GO TO 210 GO TO 130 END SUBROUTINE SETEDG(NEDGE) USE BLKMAP ! Setup to form new triangles DO J=1,NEDGE IF(J .LT. NEDGE) THEN DO K=J+1,NEDGE IF(IEDGE(K,1) .EQ. IEDGE(J,1)) THEN IF(IEDGE(K,2) .EQ. IEDGE(J,2)) THEN IEDGE(J,1)=0 IEDGE(J,2)=0 IEDGE(K,1)=0 IEDGE(K,2)=0 ENDIF ELSEIF(IEDGE(K,1) .EQ. IEDGE(J,2)) THEN IF(IEDGE(K,2) .EQ. IEDGE(J,1)) THEN IEDGE(J,1)=0 IEDGE(J,2)=0 IEDGE(K,1)=0 IEDGE(K,2)=0 ENDIF ENDIF ENDDO ENDIF ENDDO RETURN END SUBROUTINE TESTANG(XPT,YPT,K,WD) USE BLKMAP REAL*8 XPT(*),YPT(*) REAL WD(*) DATA PI/3.14159/ ! IF(WD(NOPEL(K,1)) .EQ. WD(NOPEL(K,2)) .and. wd(nopel(k,1)) .gt. -9990. ) THEN ! RETURN ! ENDIF ! IFD=0 DO N=1,NELTS IF(N .NE. K) THEN DO J=1,3 IF(NOPEL(K,1) .EQ. NOPEL(N,J)) THEN IF(J .GT. 1) THEN IF(NOPEL(K,2) .EQ. NOPEL(N,J-1)) THEN IFD=N ISIDE=J GO TO 400 ENDIF ELSE IF(NOPEL(K,2) .EQ. NOPEL(N,3)) THEN IFD=N ISIDE=J GO TO 400 ENDIF ENDIF ENDIF ENDDO ENDIF ENDDO RETURN 400 CONTINUE J1=ISIDE+1 IF(J1 .GT. 3) J1=1 ! ! IF(WD(NOPEL(K,3)) .EQ. WD(NOPEL(IFD,J1)) .and. wd(nopel(k,1)) .gt. -9990. ) THEN ! WRITE(148,'(12I8)') K,J,IFD,NOPEL(K,1),NOPEL(K,2),NOPEL(K,3),NOPEL(IFD,1),NOPEL(IFD,2),NOPEL(IFD,3) ! IF(NELTS .GT. 20) WRITE(148,'(12I8)') NOPEL(21,1),NOPEL(21,2),NOPEL(21,3) ! NOPEL(IFD,1)=NOPEL(K,3) ! NOPEL(IFD,2)=NOPEL(K,1) ! NOPEL(IFD,3)=NOPEL(IFD,J1) ! NOPEL(K,1)=NOPEL(IFD,3) ! ! WRITE(148,'(12I8)') K,J,IFD,NOPEL(K,1),NOPEL(K,2),NOPEL(K,3),NOPEL(IFD,1),NOPEL(IFD,2),NOPEL(IFD,3) ! IF(NELTS .GT. 20) WRITE(148,'(12I8)') NOPEL(21,1),NOPEL(21,2),NOPEL(21,3) ! ! CALL CCENTRE(XPT(NOPEL(IFD,1)),XPT(NOPEL(IFD,2)),XPT(NOPEL(IFD,3)) & !& ,YPT(NOPEL(IFD,1)),YPT(NOPEL(IFD,2)),YPT(NOPEL(IFD,3)) & !& ,XCEN(IFD),YCEN(IFD),RADS(K)) ! CALL CCENTRE(XPT(NOPEL(K,1)),XPT(NOPEL(K,2)),XPT(NOPEL(K,3)) & !& ,YPT(NOPEL(K,1)),YPT(NOPEL(K,2)),YPT(NOPEL(K,3)) & !& ,XCEN(K),YCEN(K),RADS(K)) ! ! RETURN ! ENDIF A1=ATAN2(YPT(NOPEL(K,1))-YPT(NOPEL(K,3)),XPT(NOPEL(K,1))-XPT(NOPEL(K,3))) A2=ATAN2(YPT(NOPEL(K,2))-YPT(NOPEL(K,3)),XPT(NOPEL(K,2))-XPT(NOPEL(K,3))) IF(A1 .LT. 0) A1=A1+2.*PI IF(A2 .LT. 0) A2=A2+2.*PI DIFFA=A2-A1 ! WRITE(148,*) 'DIFFA',K,DIFFA,NOPEL(K,1),NOPEL(K,2),NOPEL(K,3) IF(DIFFA .LT. 0) DIFFA=DIFFA+PI*2. IF(DIFFA .LT. 2./3.*PI) RETURN ! WRITE(148,'(9I8)') K,J,IFD,NOPEL(K,1),NOPEL(K,2),NOPEL(K,3),NOPEL(IFD,1),NOPEL(IFD,2),NOPEL(IFD,3) B1=ATAN2(YPT(NOPEL(IFD,J1))-YPT(NOPEL(K,2)),XPT(NOPEL(IFD,J1))-XPT(NOPEL(K,2))) B2=ATAN2(YPT(NOPEL(K, 3))-YPT(NOPEL(K,2)),XPT(NOPEL(K, 3))-XPT(NOPEL(K,2))) IF(B1 .LT. 0) B1=B1+2.*PI IF(B2 .LT. 0) B2=B2+2.*PI DIFFB=B2-B1 ! WRITE(148,*) 'DIFFB',DIFFB,B2,B1 IF(DIFFB .LT. 0) DIFFB=DIFFB+2.*PI IF(DIFFB .GT. DIFFA) RETURN C1=ATAN2(YPT(NOPEL(K, 3))-YPT(NOPEL(K,1)),XPT(NOPEL(K, 3))-XPT(NOPEL(K,1))) C2=ATAN2(YPT(NOPEL(IFD,J1))-YPT(NOPEL(K,1)),XPT(NOPEL(IFD,J1))-XPT(NOPEL(K,1))) IF(C1 .LT. 0) C1=C1+2.*PI IF(C2 .LT. 0) C2=C2+2.*PI DIFFC=C2-C1 ! WRITE(148,*) 'DIFFC',DIFFC,C2,C1 IF(DIFFC .LT. 0) DIFFC=DIFFC+2.*PI IF(DIFFC .GT. DIFFA) RETURN NOPEL(IFD,1)=NOPEL(K,3) NOPEL(IFD,2)=NOPEL(K,1) NOPEL(IFD,3)=NOPEL(IFD,J1) NOPEL(K,1)=NOPEL(IFD,3) ! WRITE(148,'(9I8)') K,J,IFD,NOPEL(K,1),NOPEL(K,2),NOPEL(K,3),NOPEL(IFD,1),NOPEL(IFD,2),NOPEL(IFD,3) CALL CCENTRE(XPT(NOPEL(IFD,1)),XPT(NOPEL(IFD,2)),XPT(NOPEL(IFD,3)) & &,YPT(NOPEL(IFD,1)),YPT(NOPEL(IFD,2)),YPT(NOPEL(IFD,3)) & &,XCEN(IFD),YCEN(IFD),RADS(IFD)) RETURN END SUBROUTINE TESTTR(XPT,YPT,K,WD) USE BLKMAP REAL WD(*) REAL*8 XPT(*),YPT(*) DATA PI/3.14159/ IF(WD(NOPEL(K,1)) .EQ. WD(NOPEL(K,2)) .and. wd(nopel(k,1)) .gt. -9990. ) THEN RETURN ENDIF IFD=0 DO N=1,NELTS IF(N .NE. K) THEN DO J=1,3 IF(NOPEL(K,2) .EQ. NOPEL(N,J)) THEN IF(J .GT. 1) THEN IF(NOPEL(K,3) .EQ. NOPEL(N,J-1)) THEN IFD=N ISIDE=J GO TO 400 ENDIF ELSE IF(NOPEL(K,3) .EQ. NOPEL(N,3)) THEN IFD=N ISIDE=J GO TO 400 ENDIF ENDIF ENDIF ENDDO ENDIF ENDDO RETURN 400 CONTINUE WRITE(148,'(9I8)') K,J,IFD,NOPEL(K,1),NOPEL(K,2),NOPEL(K,3),NOPEL(IFD,1),NOPEL(IFD,2),NOPEL(IFD,3) write(148,'(9x,6f8.0)')wd(NOPEL(K,1)),wd(NOPEL(K,2)),wd(NOPEL(K,3)),wd(NOPEL(IFD,1)),wd(NOPEL(IFD,2)),wd(NOPEL(IFD,3)) J1=ISIDE+1 IF(J1 .GT. 3) J1=1 WRITE(148,*) J1 IF(WD(NOPEL(K,1)) .EQ. WD(NOPEL(IFD,J1)) .and. wd(nopel(k,1)) .gt. -9990. ) THEN ITEMP=NOPEL(IFD,J1) NOPEL(IFD,1)=NOPEL(K,3) NOPEL(IFD,2)=NOPEL(K,1) NOPEL(IFD,3)=ITEMP NOPEL(K,1)=NOPEL(IFD,3) WRITE(148,'(9I8)') K,J,IFD,NOPEL(K,1),NOPEL(K,2),NOPEL(K,3),NOPEL(IFD,1),NOPEL(IFD,2),NOPEL(IFD,3) CALL CCENTRE(XPT(NOPEL(IFD,1)),XPT(NOPEL(IFD,2)),XPT(NOPEL(IFD,3)) & & ,YPT(NOPEL(IFD,1)),YPT(NOPEL(IFD,2)),YPT(NOPEL(IFD,3)) & & ,XCEN(IFD),YCEN(IFD),RADS(IFD)) RETURN ENDIF A1=ATAN2(YPT(NOPEL(K,2))-YPT(NOPEL(K,1)),XPT(NOPEL(K,2))-XPT(NOPEL(K,1))) A2=ATAN2(YPT(NOPEL(K,3))-YPT(NOPEL(K,1)),XPT(NOPEL(K,3))-XPT(NOPEL(K,1))) IF(A1 .LT. 0) A1=A1+2.*PI IF(A2 .LT. 0) A2=A2+2.*PI DIFFA=A2-A1 ! WRITE(148,*) 'DIFFA',K,DIFFA,NOPEL(K,1),NOPEL(K,2),NOPEL(K,3) IF(DIFFA .LT. 0) DIFFA=DIFFA+PI*2. IF(DIFFA .LT. 2./3.*PI) RETURN B1=ATAN2(YPT(NOPEL(IFD,J1))-YPT(NOPEL(K,3)),XPT(NOPEL(IFD,J1))-XPT(NOPEL(K,3))) B2=ATAN2(YPT(NOPEL(K, 1))-YPT(NOPEL(K,3)),XPT(NOPEL(K, 1))-XPT(NOPEL(K,3))) IF(B1 .LT. 0) B1=B1+2.*PI IF(B2 .LT. 0) B2=B2+2.*PI DIFFB=B2-B1 ! WRITE(148,*) 'DIFFB',DIFFB,B2,B1 IF(DIFFB .LT. 0) DIFFB=DIFFB+2.*PI IF(DIFFB .GT. DIFFA) RETURN C1=ATAN2(YPT(NOPEL(K, 1))-YPT(NOPEL(K,2)),XPT(NOPEL(K, 1))-XPT(NOPEL(K,2))) C2=ATAN2(YPT(NOPEL(IFD,J1))-YPT(NOPEL(K,2)),XPT(NOPEL(IFD,J1))-XPT(NOPEL(K,2))) IF(C1 .LT. 0) C1=C1+2.*PI IF(C2 .LT. 0) C2=C2+2.*PI DIFFC=C2-C1 ! WRITE(148,*) 'DIFFC',DIFFC,C2,C1 IF(DIFFC .LT. 0) DIFFC=DIFFC+2.*PI IF(DIFFC .GT. DIFFA) RETURN ITEMP=NOPEL(IFD,J1) NOPEL(IFD,1)=NOPEL(K,1) NOPEL(IFD,2)=NOPEL(K,2) NOPEL(IFD,3)=ITEMP NOPEL(K,2)=NOPEL(IFD,3) ! WRITE(148,'(9I8)') K,J,IFD,NOPEL(K,1),NOPEL(K,2),NOPEL(K,3),NOPEL(IFD,1),NOPEL(IFD,2),NOPEL(IFD,3) CALL CCENTRE(XPT(NOPEL(IFD,1)),XPT(NOPEL(IFD,2)),XPT(NOPEL(IFD,3)) & &,YPT(NOPEL(IFD,1)),YPT(NOPEL(IFD,2)),YPT(NOPEL(IFD,3)) & &,XCEN(IFD),YCEN(IFD),RADS(IFD)) RETURN END