contri.f

作者：Scott Sloan 相关文献参考： Adv Eng Software 1987 V9No1,34-

            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
*---------------------------------------------------------------------
 1000 FORMAT('Vertex',i5,' not in any triangle in subroutine EDGE')
      END
************************************************************************
      SUBROUTINE TCHECK(NPTS,N,X,Y,LIST,V,T,NTRI,NEF,NB,NCE,NCB,ELIST,W,
     &                  maxelist)
      use mod_output
      implicit double precision (a-h,o-z)
************************************************************************
*
*     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 maxelist
      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,maxelist)
      character(len=80) expl(19)
*
      DOUBLE PRECISION X1,X2,X3,X4,Y1,Y2,Y3,Y4
      DOUBLE PRECISION DET,DIS,R21,R31,RAD,TOL,X21,X31,XCC,Y21,Y31,YCC
      DOUBLE PRECISION C00000,CP5000
      DOUBLE PRECISION X(NPTS+3),Y(NPTS+3)
*
      PARAMETER (TOL=1.d-5)
      PARAMETER (C00000=0.d0, CP5000=0.5d0)

*-----------------------------------------------------------------------
*     Load error explanation in case its needed
*
      expl(1)='Check the following conditions:'
      expl(2)='* Edges which define boundaries must come first in'
      expl(3)='  ELIST and thus occupy the first NCB columns'
      expl(4)='* Edges which define an external boundary must be'
      expl(5)='  listed anticlockwise but may be presented in any order'
      expl(6)='* Edges which define an internal boundary (hole) must'
      expl(7)='  be listed clockwise but may be presented in any order'
      expl(8)='* An internal boundary (hole) cannot be specified'
      expl(9)='  unless an external boundary is also specified'
      expl(10)='* All boundaries must form closed loops'
      expl(11)='* An edge may not appear more than once in ELIST'
      expl(12)='* An external or internal boundary may not cross'
      expl(13)='  itself and may not share a common edge with any'
      expl(14)='  other boundary'
      expl(15)='* Internal edges, which are not meant to define'
      expl(16)='  boundaries but must be present in the final'
      expl(17)='  triangulation, occupy columns NCB+1,... ,NCE of ELIST'
      expl(18)='* No point in the list vector may lie outside any'
      expl(19)='  external boundary or inside any internal boundary'
*-----------------------------------------------------------------------
*     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,L)
              V2=V(1,L)
              V3=V(2,L)
            ENDIF
*
*           Triangle L is left of V1-V2 and is defined by V3-V1-V2 
*           Triangle R is right of V1-V2 and is defined by V4-V2-V1
*
            IF(T(1,R).EQ.L)THEN
              V4=V(3,R)
            ELSEIF(T(2,R).EQ.L)THEN
              V4=V(1,R)
            ELSE
              V4=V(2,R)
            ENDIF
*
*           Find circumcentre of triangle L and its circumcircle radius
*
            X1 =X(V1)
            Y1 =Y(V1)
            X21=X(V2)-X1
            Y21=Y(V2)-Y1
            X31=X(V3)-X1
            Y31=Y(V3)-Y1
            DET=X21*Y31-X31*Y21
            IF(DET.LE.C00000)THEN
              write(string,
     &        '(''* WARNING: Zero or negative area for triangle'',i5)')L
              call output
              write(string,'(11x,a)')'in subroutine TCHECK'
              call output
            ELSE
              DET=CP5000/DET
              R21=X21*X21+Y21*Y21
              R31=X31*X31+Y31*Y31
              XCC=DET*(R21*Y31-R31*Y21)
              YCC=DET*(X21*R31-X31*R21)
              RAD=SQRT(XCC*XCC+YCC*YCC)
              XCC=XCC+X1
              YCC=YCC+Y1
*
*             Check if V4 is inside circumcircle for triangle L
*
              DIS=SQRT((XCC-X(V4))**2.d0+(YCC-Y(V4))**2.d0)
              IF(RAD-DIS.GT.TOL*RAD)THEN
                NOPT=NOPT+1
              ENDIF
            ENDIF
          ENDIF
   10   CONTINUE
   20 CONTINUE
*-----------------------------------------------------------------------
*     Check triangulation is valid
*
      IF(NBOV.LT.3)THEN
        nerror = 2
        error(1) =
     &    'Number of boundary vertices less than 3 in subroutine TCHECK'
        write(error(2),'(''NBOV='',i5)')NBOV
        call erroroutput
        STOP
      ENDIF
      IF(NTRI.NE.2*(N+NB)-NBOV-4)THEN
        nerror = 21
        error(1) = 'Invalid triangulation in subroutine TCHECK'
        error(2) = 'NTRI is not equal to 2*(N+NB)-NBOV-4'
        do i=1,19
           error(2+i)=expl(i)
        enddo
        call erroroutput
        STOP
      ENDIF
      IF(NEDG.NE.N+NTRI+NB-2)THEN
        nerror = 21
        error(1) = 'Invalid triangulation in subroutine TCHECK'
        error(2) = 'NEDG is not equal to N+NTRI+NB-2'
        do i=1,19
           error(2+i)=expl(i)
        enddo
        call erroroutput
        STOP
      ENDIF
      IF(NOPT.GT.NEF)THEN
        nerror = 21
        error(1) = 'Invalid triangulation in subroutine TCHECK'
        error(2) = 'Too many non-optimal edges'
        do i=1,19
           error(2+i)=expl(i)
        enddo
        call erroroutput
        STOP    
      ENDIF
      IF(NCB.GT.0)THEN
        IF(NCB.NE.NBOV)THEN
          nerror = 21
          error(1) = 'Invalid triangulation in subroutine TCHECK'
          error(2) = 'Number of boundary vertices not equal to NBC'
          do i=1,19
             error(2+i)=expl(i)
          enddo
          call erroroutput
          STOP    
        ENDIF
      ENDIF
*----------------------------------------------------------------------
*     Check that each node appears in at least one triangle
*
      DO 25 I=1,NPTS
        W(I)=0
   25 CONTINUE
      DO 35 J=1,NTRI