contri.f

Programs in the irregular grid design package described in this manual are used to carry out five ma

          X1=X(V1)
          Y1=Y(V1)
          X2=X(V2)
          Y2=Y(V2)
          IF((X3-X1)*(Y4-Y1).LT.(X4-X1)*(Y3-Y1))THEN
            IF((X3-X2)*(Y4-Y2).GT.(X4-X2)*(Y3-Y2))THEN
*
*             Quad is convex so swap diagonal arcs
*             Update vertex and adjacency lists for triangle L
*
              IF(ELR.EQ.1)THEN
                V(2,L)=V4
                C     =T(2,L)
                T(1,L)=A
                T(2,L)=R
              ELSEIF(ELR.EQ.2)THEN
                V(3,L)=V4
                C     =T(3,L)
                T(2,L)=A
                T(3,L)=R             
              ELSE
                V(1,L)=V4
                C     =T(1,L)
                T(3,L)=A
                T(1,L)=R             
              ENDIF                 
*
*             Update vertex and adjacency lists for triangle R
*
              IF(ERL.EQ.1)THEN
                V(2,R)=V3
                T(1,R)=C
                T(2,R)=L
              ELSEIF(ERL.EQ.2)THEN
                V(3,R)=V3
                T(2,R)=C
                T(3,R)=L             
              ELSE
                V(1,R)=V3
                T(3,R)=C
                T(1,R)=L             
              ENDIF
*
*             Update adjacency lists for triangles A and C
*
              IF(T(1,C).EQ.L)THEN
                T(1,C)=R
              ELSEIF(T(2,C).EQ.L)THEN
                T(2,C)=R
              ELSE
                T(3,C)=R
              ENDIF
              IF(T(1,A).EQ.R)THEN
                T(1,A)=L
              ELSEIF(T(2,A).EQ.R)THEN
                T(2,A)=L
              ELSE
                T(3,A)=L
              ENDIF
*
*             Update vertex-triangle list
*
              TN(V1)=L
              TN(V2)=R
*
*             Test if new diagonal arc crosses VI-VJ and store it if it
*             does
*
              IF(((XI-X3)*(YJ-Y3)-(XJ-X3)*(YI-Y3))*
     +          ((XI-X4)*(YJ-Y4)-(XJ-X4)*(YI-Y4)).LT.C00000)THEN
                W(1,FIRST)=V4
                W(2,FIRST)=V3
                FIRST=FIRST+1
              ELSE
                W(1,FIRST)=W(1,LAST)
                W(1,LAST) =V3
                W(2,FIRST)=W(2,LAST)
                W(2,LAST) =V4        
                LAST=LAST-1                    
              ENDIF
              GOTO 40
            ENDIF
          ENDIF
*
*         Arc cannot be swapped, so move to next intersecting arc
*
          FIRST=FIRST+1
          GOTO 40
        ENDIF
        GOTO 35
      ENDIF
*----------------------------------------------------------------------
*     Optimise all new arcs (except VI-VJ)
*        
      SWAP=.TRUE.
   50 IF(SWAP)THEN
        SWAP=.FALSE.
        DO 70 I=2,NC
*
*         Find triangle L which is left of V1-V2 
*         Find triangle R which is right of V1-V2
*
          V1=W(1,I)
          V2=W(2,I)
*
*         Exchange V1 and V2 if V1 is a supertriangle vertex
*
          IF(V1.GT.NPTS)THEN
            IF(V2.GT.NPTS)THEN
              WRITE(66,'(//,''***ERROR IN EDGE***'',
     +                   /,''ARC BETWEEN VERTICES'',i8,'' AND'',i8,
     +                   /,''CANNOT BE OPTIMISED SINCE IT IS A '',
     +                     ''SUPERTRIANGLE BOUNDARY'')')V1,V2
               RETURN !STOP
            ENDIF
            W(1,I)=V2
            W(2,I)=V1
            V2=V1
            V1=W(1,I)
          ENDIF
          L=TN(V1)
   60     IF(V(1,L).EQ.V1)THEN
            IF(V(2,L).EQ.V2)THEN
              V3 =V(3,L)
              ELR=1
              R  =T(1,L)
            ELSE
              L=T(3,L)
              GOTO 60
            ENDIF
          ELSEIF(V(2,L).EQ.V1)THEN
            IF(V(3,L).EQ.V2)THEN
              V3 =V(1,L)
              ELR=2
              R  =T(2,L)
            ELSE
              L=T(1,L)
              GOTO 60
            ENDIF
          ELSEIF(V(1,L).EQ.V2)THEN
              V3 =V(2,L)
              ELR=3
              R  =T(3,L)
          ELSE
              L=T(2,L)
              GOTO 60
          ENDIF
*
*         Find vertices V3 and V4 where:
*         triangle L is defined by V3-V1-V2 
*         triangle R is defined by V4-V2-V1
*
          IF(T(1,R).EQ.L)THEN
            V4=V(3,R)
            A =T(2,R)
            ERL=1
          ELSEIF(T(2,R).EQ.L)THEN
            V4=V(1,R)
            A =T(3,R)
            ERL=2
          ELSE
            V4=V(2,R)
            A =T(1,R)
            ERL=3
          ENDIF
          X13=X(V1)-X(V3)
          Y13=Y(V1)-Y(V3)
          X14=X(V1)-X(V4)
          Y14=Y(V1)-Y(V4)
          X23=X(V2)-X(V3)
          Y23=Y(V2)-Y(V3)
          X24=X(V2)-X(V4)
          Y24=Y(V2)-Y(V4)
          COSA=X13*X23+Y23*Y13
          COSB=X24*X14+Y24*Y14
          IF(COSA.LT.C00000.OR.COSB.LT.C00000)THEN
            IF((COSA.LT.C00000.AND.COSB.LT.C00000).OR.
     +        ((X13*Y23-X23*Y13)*COSB-(X14*Y24-X24*Y14)*COSA.LT.
     +         -TOL*SQRT((X13*X13+Y13*Y13)*(X23*X23+Y23*Y23)*
     +                   (X24*X24+Y24*Y24)*(X14*X14+Y14*Y14))))THEN
*
*             V4 is inside circumcircle for triangle L
*             Swap diagonal for convex quad formed by V3-V1-V4-V2
*             Update vertex and adjacency lists for triangle L
*
              SWAP=.TRUE.
              IF(ELR.EQ.1)THEN
                V(2,L)=V4
                C     =T(2,L)
                T(1,L)=A
                T(2,L)=R
              ELSEIF(ELR.EQ.2)THEN
                V(3,L)=V4
                C     =T(3,L)
                T(2,L)=A
                T(3,L)=R             
              ELSE
                V(1,L)=V4
                C     =T(1,L)
                T(3,L)=A
                T(1,L)=R             
              ENDIF                 
*
*             Update vertex and adjacency lists for triangle R
*
              IF(ERL.EQ.1)THEN
                V(2,R)=V3
                T(1,R)=C
                T(2,R)=L
              ELSEIF(ERL.EQ.2)THEN
                V(3,R)=V3
                T(2,R)=C
                T(3,R)=L             
              ELSE
                V(1,R)=V3
                T(3,R)=C
                T(1,R)=L             
              ENDIF
*
*             Update adjacency lists for triangles A and C
*
              IF(T(1,C).EQ.L)THEN
                T(1,C)=R
              ELSEIF(T(2,C).EQ.L)THEN
                T(2,C)=R
              ELSE
                T(3,C)=R
              ENDIF
              IF(T(1,A).EQ.R)THEN
                T(1,A)=L
              ELSEIF(T(2,A).EQ.R)THEN
                T(2,A)=L
              ELSE
                T(3,A)=L
              ENDIF
*
*             Update vertex-triangle list and arc list
*
              TN(V1)=L
              TN(V2)=R
              W(1,I)=V3
              W(2,I)=V4
            ENDIF
          ENDIF
   70   CONTINUE
        GOTO 50
      ENDIF
	EDGE = .TRUE.
*---------------------------------------------------------------------
 1000 FORMAT(//,'***ERROR IN EDGE***',
     +        /,'   VERTEX',i8,' NOT IN ANY TRIANGLE')
      END
c----------
      FUNCTION TCHECK(NPTS,N,X,Y,LIST,V,T,NTRI,NEF,NB,NCE,NCB,ELIST,W,
     +    xmin,ymin,dmax)
	LOGICAL TCHECK
************************************************************************
*
*     PURPOSE:
*     --------
*
*     Check Delaunay triangulation which may be constrained
*
*     INPUT:
*     ------
*
*     NPTS   - Total number of points in data set
*     N      - Total number of points to be triangulated (N le NPTS)
*     X      - X-coords of all points in data set
*            - X-coord of point I given by X(I)
*            - Last three locations are used to store x-coords of
*              supertriangle vertices in subroutine delaun
*     Y      - Y-coords of all points in data set
*            - Y-coord of point I given by Y(I)
*            - Last three locations are used to store y-coords of
*              supertriangle vertices in subroutine delaun          
*     LIST   - List of points to be triangulated
*     V      - Vertex array for triangulation
*            - Vertices listed in anticlockwise sequence
*            - Vertices for triangle J are found in V(I,J) for I=1,2,3
*              and J=1,2,...,NTRI
*            - First vertex is at point of contact of first and third
*              adjacent triangles
*     T      - Adjacency array for triangulation
*            - Triangles adjacent to J are found in T(I,J) for I=1,2,3
*              J=1,2,...,NTRI
*            - Adjacent triangles listed in anticlockwise sequence
*            - Zero denotes no adjacent triangle
*     NTRI   - Number of triangles in triangulation
*     NEF    - Number of forced edges in triangulation
*     NB     - Number of boundaries defining the triangulation
*            - NB=1 for a simply connected domain with no holes 
*            - NB=H+1 for a domain with H holes
*     NCE    - Total number of constrained edges which must be 
*              present, including those which define boundaries
*            - NCE=0 indicates triangulation is unconstrained so that
*              ELIST is empty and the code will produce a 
*              triangulation of the convex hull
*     NCB    - Total number of constrained edges which define one
*              external boundary and any internal boundaries (holes)
*            - NCB=0 indicates there are no boundary edge constraints
*              and the code will produce a triangulation of the convex
*              hull
*     ELIST  - List of edges which must be present in triangulation
*            - These may be internal edges or edges which define 
*              boundaries
*     W      - Undefined, vector used as workspace
*
*     OUTPUT:
*     ------- 
*
*     NPTS   - Unchanged
*     N      - Unchanged
*     X      - Unchanged
*     Y      - Unchanged
*     LIST   - Unchanged
*     V      - Unchanged
*     T      - Unchanged
*     NTRI   - Unchanged
*     NEF    - Unchanged
*     NB     - Unchanged
*     NCE    - Unchanged
*     NCB    - Unchanged
*     ELIST  - Unchanged
*     W      - Not used
*
*     NOTES:
*     ------
*     
*     - This routine performs a number of elementary checks to test
*       the integrity of the triangulation
*     - NTRI=2*(N+NB)-NBOV-4 for a valid triangulation, where NBOV is
*       the number of boundary vertices
*     - NEDG=N+NTRI+NB-2 for a valid triangulation, where NEDG is the
*       number of edges in the triangulation
*     - NOPT le NEF for a valid triangulation, where NOPT is the number
*       of non-optimal edges in the triangulation
*     - The triangulation is tested to ensure that each non-boundary
*       constrained edge (if there are any) is present
*
*     PROGRAMMER:             Scott Sloan
*     -----------
*
*     LAST MODIFIED:          3 march 1991          Scott Sloan
*     --------------
*
*     COPYRIGHT 1990:         Scott Sloan
*     ---------------         Department of Civil Engineering
*                             University of Newcastle
*                             NSW 2308
*
************************************************************************
      INTEGER I,J,L,N,R,S
      INTEGER NB,V1,V2,V3,V4
      INTEGER NBOV,NEDG,NOPT,NPTS,NTRI
      INTEGER NCB,NCE,NEF
      INTEGER T(3,NTRI),V(3,NTRI),W(NPTS)
      INTEGER LIST(N)
      INTEGER ELIST(2,NCE)
c      INTEGER ELIST(2,*)
      real xmin,ymin,dmax
*
      REAL X1,X2,X3,X4,Y1,Y2,Y3,Y4
      REAL DET,DIS,R21,R31,RAD,TOL,X21,X31,XCC,Y21,Y31,YCC
      REAL C00000,CP5000
      REAL X(NPTS+3),Y(NPTS+3)
*
c      PARAMETER (TOL=1.E-5)  !32 bit
      PARAMETER (TOL=1.E-16)   !64 bit
      PARAMETER (C00000=0.0, CP5000=0.5)

	TCHECK = .FALSE.
	TCHECK = .TRUE.
*-----------------------------------------------------------------------
*     Loop over each triangle and count the following
*     NOPT=number of edges which are not optimal
*     NEDG=number of edges in triangulation
*     NBOV=number of boundary vertices
*
      NOPT=0
      NEDG=0
      NBOV=0
      DO 20 L=1,NTRI
*
*       Loop over each side of triangle
*
        DO 10 I=1,3
          R=T(I,L)
          IF(R.EQ.0)THEN
            NEDG=NEDG+1
            NBOV=NBOV+1
          ELSEIF(R.LT.L)THEN
            NEDG=NEDG+1                                   
            IF(I.EQ.1)THEN
              V1=V(1,L)
              V2=V(2,L)
              V3=V(3,L)
            ELSEIF(I.EQ.2)THEN
              V1=V(2,L)
              V2=V(3,L)
              V3=V(1,L)
            ELSE
              V1=V(3,