contri.f

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

        VJ=ELIST(2,I)
        IF(.NOT. EDGE(VI,VJ,NPTS,NTRI,NEF,JW,X,Y,V,T,W,W(1,2),
     +        XMIN,YMIN,DMAX))THEN
          call PigPutMessage('Triangulation error. See gridit.log.')
          call PigUWait(3.0)
		  RETURN
	  END IF
  140 CONTINUE
	CONTRI = .TRUE.
*---------------------------------------------------------------------
*     Clip triangulation and check it
*
      IF(.NOT. BCLIP(NPTS,NCB,ELIST,W,V,T,NTRI,NB,
     +           nerrors, maxerrors, err_elist) )THEN
        call PigPutMessage('Triangulation error. See gridit.log.')
C        call PigUWait(3.0)
		CONTRI = .FALSE.
	END IF

      if (contri)then
        IF(.NOT.
     +  	TCHECK(NPTS,N,X,Y,LIST,V,T,NTRI,NEF,NB,NCE,NCB,ELIST,W,
     +        xmin,ymin,dmax)
     +	  ) THEN
          call PigPutMessage('Triangulation error. See gridit.log.')
C          call PigUWait(3.0)
		  CONTRI = .FALSE.
	  END IF
      end if
*---------------------------------------------------------------------
*     Reset x-y coords to original values
*
      DO 150 I=1,N
        P=LIST(I)
        X(P)=X(P)*DMAX+XMIN
        Y(P)=Y(P)*DMAX+YMIN
  150 CONTINUE
c loop below added by a g dolling 17 October 1997
c to return supertriangle in original coordinates
      DO P=NPTS+1,NPTS+3
        X(P)=X(P)*DMAX+XMIN
        Y(P)=Y(P)*DMAX+YMIN
      end do
*---------------------------------------------------------------------
 1000 FORMAT(//,'***ERROR IN CONTRI***',
     +        /,'ILLEGAL VALUE IN LIST',
     +        /,'LIST(',i8,' )=',i8)
 2000 FORMAT(//,'***ERROR IN CONTRI***',
     +        /,'ILLEGAL VALUE IN ELIST',
     +        /,'ELIST(',i8,',',i8,' )=',i8)
 3000 FORMAT(//,'***ERROR IN CONTRI***',
     +        /,'ILLEGAL VALUE IN ELIST',
     +        /,'ELIST(',i8,',',i8,' )=',i8,
     +        /,'THIS POINT IS NOT IN LIST')
      END 
c----------
      FUNCTION BCLIP(NPTS,NCB,BLIST,TN,V,T,NTRI,NB,
     +           nerrors, maxerrors, err_elist)
	LOGICAL BCLIP
***********************************************************************
*
*     PURPOSE:
*     --------
*
*     Clip constrained Delaunay triangulation to boundaries in BLIST
*     If BLIST is empty, then the triangulation is clipped to a convex
*     hull by removing all triangles that are formed with the 
*     supertriangle vertices
*     The triangulation MUST initially include the supertriangle 
*     vertices
*
*     INPUT:
*     ------
*
*     NPTS   - Total number of points in data set (NPTS ge N)
*     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
*     BLIST  - List of edges which must be present in triangulation
*              and define boundaries
*            - Edge I defined by vertices in BLIST(1,I) and BLIST(2,I)
*              where I ranges from 1,...,NCB
*            - Edges which define an external boundary must be listed
*              anticlockwise but may be presented in any order
*            - Edges which define an internal boundary (hole) must be 
*              listed clockwise but may be presented in any order
*     TN     - List of triangle numbers such that vertex I can be 
*              found in triangle TN(I)
*     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
*              and J=1,2,...,NTRI
*            - Adjacent triangles listed in anticlockwise sequence
*            - Zero denotes no adjacent triangle
*     NTRI   - Number of triangles, including those formed with
*              vertices of supertriangle
*     NB     - Not defined
*
*     OUTPUT:
*     -------
*
*     NPTS   - Unchanged
*     NCB    - Unchanged
*     BLIST  - Unchanged
*     TN     - Not defined 
*     V      - Updated 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      - Updated adjacency array for triangulation
*            - Triangles adjacent to J are found in T(I,J) for I=1,2,3
*              and J=1,2,...,NTRI
*            - Adjacent triangles listed in anticlockwise sequence
*            - Zero denotes no adjacent triangle
*     NTRI   - Updated number of triangles in final 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
*
*     PROGRAMMER:           Scott Sloan
*     -----------
*
*     LAST MODIFIED:        3 march 1991        Scott Sloan
*     --------------
*                           2 october 1997      Adrian Dolling
*                           Added NTRI_OUT to allow debugging / release
*                           with array bounds checking turned on under
*                           MS Fortran PowerStation 4.0
*                           V is dimensioned V(3, NTRI), then NTRI is set
*                           to zero which triggers the bounds checking
*                           error detection
*
*     COPYRIGHT 1990:       Scott Sloan
*     ---------------       Department of Civil Engineering
*                           University of Newcastle
*                           NSW 2308
*
***********************************************************************
      INTEGER A,C,E,I,J,L,R,S
      INTEGER NB,TS,V1,V2
      INTEGER NCB,NEV
      INTEGER NPTS,NTRI
      INTEGER NNTRI,NPTS1,NTRI1
      INTEGER T(3,NTRI),V(3,NTRI)
      INTEGER TN(NPTS+3)
      INTEGER BLIST(2,*)
	INTEGER NTRI_OUT

	integer nerrors, maxerrors
	integer err_elist(maxerrors)
*---------------------------------------------------------------------
*     Skip to triangle deletion step if triangulation has no
*     boundary constraints
*
	BCLIP = .FALSE.
      IF(NCB.EQ.0)THEN
        NB=1
        GOTO 55
      ENDIF
*---------------------------------------------------------------------
*     Mark all edges which define the boundaries
*     S=triangle in which search starts
*     R=triangle to right of edge V1-V2
*     
      DO 20 I=1,NCB
        V1=BLIST(1,I)
        V2=BLIST(2,I)
        S=TN(V1)
        R=S
*
*       Circle anticlockwise round V1 until edge V1-V2 is found
*
   10   IF(V(1,R).EQ.V1)THEN
          IF(V(3,R).EQ.V2)THEN
            T(3,R)=-T(3,R)
            GOTO 20
          ELSE
            R=ABS(T(3,R))
          ENDIF
        ELSEIF(V(2,R).EQ.V1)THEN
          IF(V(1,R).EQ.V2)THEN
            T(1,R)=-T(1,R)
            GOTO 20
          ELSE
            R=ABS(T(1,R))
          ENDIF
        ELSEIF(V(2,R).EQ.V2)THEN
            T(2,R)=-T(2,R)
            GOTO 20
        ELSE
            R=ABS(T(2,R))
        ENDIF
        IF(R.EQ.S)THEN
          WRITE(66,'(//,''***ERROR IN BCLIP***'',
     +               /,''CONSTRAINED BOUNDARY EDGE NOT FOUND'',
     +               /,''V1='',i8,'' V2='',i8,
     +               /,''CHECK THAT THIS EDGE IS NOT CROSSED'',
     +               /,''BY ANOTHER BOUNDARY EDGE'')')V1,V2
	    if(nerrors.lt.maxerrors)then
	      nerrors = nerrors + 1
	      err_elist(nerrors) = i
	    end if
          go to 20
c          RETURN !STOP
        ENDIF
        GOTO 10
   20 CONTINUE
      if(nerrors.ne.0)return
*--------------------------------------------------------------------
*     Mark all triangles which are to right of boundaries by
*     circling anticlockwise around the outside of each boundary
*     Loop while some boundary edges have not been visited
*     S = triangle from which search starts
*     NEV = number of edges visited
*
      S  =0
      NB =0
      NEV=0
      NTRI1=NTRI+1
      NPTS1=NPTS+1
   30 IF(NEV.LT.NCB)THEN
*
*       Find an edge on a boundary
*
   40   S=S+1
        IF(T(1,S).LT.0)THEN
          E=1
        ELSEIF(T(2,S).LT.0)THEN
          E=2
        ELSEIF(T(3,S).LT.0)THEN
          E=3
        ELSE
          GOTO 40
        ENDIF
*
*       Store and mark starting edge
*       Mark starting triangle for deletion 
*       Increment count of edges visited
*       C = current triangle
*
        TS =T(E,S)
        T(E,S)=NTRI1
        V(1,S)=NPTS1
        NEV=NEV+1
        C=S
*
*       Increment to next edge
*
        E=MOD(E+1,3)+1
*
*       Loop until trace of current boundary is complete
*
   50   IF(T(E,C).NE.NTRI1)THEN
          IF(T(E,C).LT.0)THEN
*
*           Have found next boundary edge
*           Increment count of boundary edges visited, unmark the edge
*           and move to next edge 
*
            NEV=NEV+1
            T(E,C)=-T(E,C)
            E=MOD(E+1,3)+1
          ELSE
*
*           Non-boundary edge, circle anticlockwise around boundary 
*           vertex to move to next triangle, and mark next
*           triangle for deletion
*
            L=C
            C=T(E,L)
            IF(T(1,C).EQ.L)THEN
              E=3
            ELSEIF(T(2,C).EQ.L)THEN
              E=1
            ELSE
              E=2
            ENDIF
            V(1,C)=NPTS1
          ENDIF
          GOTO 50
        ENDIF
*
*       Trace of current boundary is complete
*       Reset marked edge to original value and check for any more
*       boundaries
*
        T(E,C)=-TS
        NB=NB+1
        GOTO 30
      ENDIF
*---------------------------------------------------------------------
*     Remove all triangles which have been marked for deletion
*     Any triangle with a vertex number greater than NPTS is deleted
*     
   55 CONTINUE
      NNTRI=NTRI
      NTRI_OUT =0
      DO 80 J=1,NNTRI
        IF(MAX(V(1,J),V(2,J),V(3,J)).GT.NPTS)THEN
*
*         Triangle J is to be deleted
*         Update adjacency lists for triangles adjacent to J
* 
          DO 60 I=1,3
            A=T(I,J)
            IF(A.NE.0)THEN
              IF(T(1,A).EQ.J)THEN
                T(1,A)=0
              ELSEIF(T(2,A).EQ.J)THEN
                T(2,A)=0
              ELSE
                T(3,A)=0
              ENDIF
            ENDIF
   60     CONTINUE
        ELSE
*
*         Triangle J is not to be deleted
*         Update count of triangles
*
          NTRI_OUT=NTRI_OUT+1
          IF(NTRI_OUT.LT.J)THEN
*
*           At least one triangle has already been deleted
*           Relabel triangle J as triangle NTRI
*
            DO 70 I=1,3
              A=T(I,J)
              T(I,NTRI_OUT)=A
              V(I,NTRI_OUT)=V(I,J)
              IF(A.NE.0)THEN
                IF(T(1,A).EQ.J)THEN
                  T(1,A)=NTRI_OUT
                ELSEIF(T(2,A).EQ.J)THEN
                  T(2,A)=NTRI_OUT
                ELSE
                  T(3,A)=NTRI_OUT
                ENDIF
              ENDIF
  70        CONTINUE
          ENDIF
        ENDIF
  80  CONTINUE 
      NTRI = NTRI_OUT
	BCLIP = .TRUE.
      END
c----------
      FUNCTION BSORT(N,NPTS,X,Y,XMIN,XMAX,YMIN,YMAX,DMAX,BIN,LIST,W)
	LOGICAL BSORT
************************************************************************
*
*     PURPOSE:
*     --------
*
*     Sort points such that consecutive points are close to one another 
*     in the x-y plane using a bin sort
*
*     INPUT:
*     ------
*
*     N      - Total number of points to be triangulated (N le NPTS)
*     NPTS   - Total number of points in data set
*     X      - X-coords of all points in data set
*            - If point is in list,the coordinate must be normalised
*              according to X=(X-XMIN)/DMAX
*            - 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
*            - If point is in list,the coordinate must be normalised
*              according to Y=(Y-YMIN)/DMAX
*            - Y-coord of point I given by Y(I)
*            - Last three locations are used to store y-coords of
*              supertriangle vertices in subroutine delaun
*     XMIN   - Min x-coord of points in LIST
*     XMAX   - Max x-coord of points in LIST
*     YMIN   - Min y-coord of points in LIST
*     YMAX   - Max y-coord of points in LIST
*     DMAX   - DMAX=MAX(XMAX-XMIN,YMAX-YMIN)
*     BIN    - Not defined
*     LIST   - List of points to be triangulated
*     W      - Undefined workspace
*
*     OUTPUT:
*     -------
*
*     N      - Unchanged
*     NPTS   - Unchanged
*     X      - Unchanged
*     Y      - Unchanged
*     XMIN   - Unchanged
*     XMAX   - Unchanged
*     YMIN   - Unchanged
*     YMAX   - Unchanged
*     DMAX   - Unchanged
*     BIN    - Not used
*     LIST   - List of points to be triangulated
*            - Points ordered such that consecutive points are close