triang.f

Go to the documentation of this file.

C> @file grid2D/triang.f
C> Triangulation of 2D grid
C GRIDMAN, grid managment library. Author: Vladislav Kotov, v.kotov@fz-juelich.de

!    Copyright (c) 2017  Forschungszentrum Juelich GmbH
!                        Vladislav Kotov 
!
!    This file is part of GRIDMAN.
!
!    GRIDMAN is free software: you can redistribute it and/or modify
!    it under the terms of the GNU General Public License as published by
!    the Free Software Foundation, either version 3 of the License, or
!    (at your option) any later version.
!
!    GRIDMAN is distributed in the hope that it will be useful,
!    but WITHOUT ANY WARRANTY; without even the implied warranty of
!    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!    GNU General Public License for more details.
!
!    You should have received a copy of the GNU General Public License
!    along with GRIDMAN.  If not, see <http://www.gnu.org/licenses/>.

C**********************************************************************************************
C> Triangulation of 2D grid
C>
C> WARNING: if the object GRID already exist it will be overwritten!
C> Same is true for mapping tables.
C>
C> The algorithm works iteratively. On each iteration one new edge 
C> is added to the cells with number of vertices > 3.  
C> Iterations stop when all cells have only 3 vertices.
C> 
C> The edge is added in the following way. For each vertex
C> a triangle is built by connecting its both neighbour vertices. 
C> Triangles are rejected for which:
C>  1. New added edge intersect other edges of the cells
C>  2. Basis vertex is concave 
C>  3. Triangle area is larger than 0.51*(total cell area) 
C>     
C> For each triangle which is not rejected the smallest angle 
C> is calculated. Triangle with the largest smallest angle is 
C> finally accepted.
C>
C> Condition 3 is required to avoid the following situation:
C>
C>            B         Here the algorithm may choose triangle ABC 
C>           / \        which leads to "flat" triangle ACD on the other side. 
C>          /   \       Condition S<0.51S_tot forces the algorithm to chose
C>         /     \      triangle ABD or BDC which leaves "non flat" rest 
C>        /___ ___\     of the cell
C>       A    D    C
C>
C**********************************************************************************************
      SUBROUTINE gridman_grid2d_triang(TRIA,GRID,IERR)
      USE gridman,ONLY:gridman_grid,gridman_sp,gridman_unit,
     u                 gridman_dbg,gridman_check,gridman_indlist
      USE gridman_lib,ONLY:gridman_grid_check,gridman_grid2d_chains,
     u                     gridman_grid_allocate,gridman_grid_take,
     u                     gridman_indlist_allocate,
     u                     gridman_indlist_deallocate,
     u                     gridman_index_allocate,gridman_index_copy
      IMPLICIT NONE
      INTRINSIC trim,abs
C> Resulting triangular grid     
      TYPE(gridman_grid) :: TRIA
C> Initial grid which is divided into triangles
      TYPE(gridman_grid) :: GRID  
C> Error code
      INTEGER,INTENT(OUT) :: IERR   
      
      INTEGER :: RES,IERR0,II,ST
      INTEGER(GRIDMAN_SP)  :: ICELL,NP,NCELLS0,NEDGES0,ICELL0,
     i                        dnc,dne,ncells,nedges,n,ie,i,nelements
      INTEGER(GRIDMAN_SP),ALLOCATABLE ::  IOLDCELL(:)
      TYPE(gridman_indlist) :: CHAINS,INEWCELL       

      IF(gridman_dbg) 
     f   WRITE(gridman_unit,*) "Starting GRIDMAN_GRID2D_TRIANG"

      ierr=0

      IF(grid%TYPE.NE.2.OR.grid%PDIM.NE.2.OR.grid%EDIM.NE.2) THEN
       ierr=100
       WRITE(gridman_unit,*) "GRIDMAN_GRID2D_TRIANG: ",
     w                       "the grid is not of 2D type"
       WRITE(gridman_unit,*) " TYPE, PDIM, EDIM ",
     w                         grid%TYPE,grid%PDIM,grid%EDIM
       RETURN 
      END IF 

      IF(grid%NCELLS.LT.1) GOTO 100

      IF(gridman_check) THEN
       CALL gridman_grid2d_check(grid,res,ierr)
       IF(res.NE.0.OR.ierr.GT.0) THEN
        ierr=100
        WRITE(gridman_unit,*) "ERROR in GRIDMAN_GRID2D_TRIANG: ",
     w                        "incorrect input grid"
        RETURN
       END IF
      END IF      

C     CALCULATE THE NUMBER OF CELLS AND EDGES IN THE TRIANGULATED GRID 
      CALL gridman_grid2d_chains(grid,chains,ierr)
      IF(ierr.GT.0) GOTO 150
      ncells=grid%NCELLS
      nedges=grid%NEDGES
      DO icell=1,chains%N
       np=chains%ILAST(icell)-chains%IFIRST(icell) !NUMBER OF POINTS 
       IF(np.LT.3) GOTO 200
       ncells=ncells+np-3
       nedges=nedges+np-3
      END DO
C     ALLOCATE GRID
      CALL gridman_grid_allocate(tria,grid%TYPE,nedges,grid%NPOINTS, 
     s                           ncells,ierr,grid%NEDGEINDEX)
      IF(ierr.GT.0) GOTO 1000
C     TAKE POINTS AND EDGES OF THE INITIAL GRID INTO NEW GRID 
      CALL gridman_grid_take(tria,grid,ierr)    
      IF(ierr.GT.0) GOTO 1000
C     EDGE INDEXES - OLD EDGES ARE NOT CHANGED
      DO ii=1,grid%NEDGEINDEX
       CALL gridman_index_copy(tria%EDGEINDEX(ii),
     c                         grid%EDGEINDEX(ii),ierr) 
       IF(ierr.NE.0) GOTO 600
      END DO
C     RESET CURRENT NUMBER OF CELLS AND EDGES IN THE NEW GRID
      ncells0=grid%NCELLS
      nedges0=grid%NEDGES
      tria%NCELLS=grid%NCELLS
      tria%NEDGES=grid%NEDGES
      IF(grid%NCELLINDEX.GT.0) THEN
C      ALLOCATE TEMPORARY ARRAY
       ALLOCATE(ioldcell(ncells),stat=st) 
       IF(st.NE.0) GOTO 300
       DO icell=1,ncells
        ioldcell(icell)=icell
       END DO
      END IF

 77   CONTINUE
C     ONE STEP OF TRIANGULATION: EACH (NON-TRIANGLE) CELL
C     IS DIVIDED INTO A TRIANGLE AND THE REST             
      IF(ALLOCATED(ioldcell)) THEN
       CALL gridman_grid2d_triang_step(tria,ncells,nedges,chains,ierr,
     c                                 ioldcell)
      ELSE
       CALL gridman_grid2d_triang_step(tria,ncells,nedges,chains,ierr)
      END IF 
      IF(ierr.GT.0) GOTO 1000
C     DIFFERENCE BETWEEN OLD AND NEW NUMBER OF CELLS (AND EDGES)
      dnc=tria%NCELLS-ncells0
      dne=tria%NEDGES-nedges0
      IF(dnc.LT.0.OR.dne.LT.0.OR.dnc.NE.dne) GOTO 400
C     IF NON-ZERO INCREMENT - THEN GO FOR THE NEXT STEP
      IF(dnc.GT.0) THEN
       ncells0=tria%NCELLS
       nedges0=tria%NEDGES 
       CALL gridman_grid2d_chains(tria,chains,ierr)
       IF(ierr.GT.0) GOTO 1000
       GOTO 77
      END IF

      CALL gridman_indlist_deallocate(chains,ierr0) 

C     CHECK 
      IF(tria%NCELLS.NE.ncells.OR.tria%NEDGES.NE.nedges) GOTO 500

C     CELL INDEXES 
      IF(ALLOCATED(ioldcell)) THEN
       CALL gridman_indlist_allocate(inewcell,ncells0,ncells,ierr)
       IF(ierr.NE.0) GOTO 700
       n=0
       DO icell0=1,ncells0
        inewcell%IFIRST(icell0)=n+1
        DO icell=1,ncells
         IF(icell0.EQ.ioldcell(icell)) THEN
          n=n+1
          IF(n.GT.ncells) GOTO 710
          inewcell%IND(n)=icell
         END IF 
        END DO 
        inewcell%ILAST(icell0)=n
       END DO
       IF(n.NE.ncells) GOTO 720 
       DEALLOCATE(ioldcell)
       ALLOCATE(tria%CELLINDEX(grid%NCELLINDEX),stat=st)  
       IF(st.NE.0) GOTO 725     
       tria%NCELLINDEX=grid%NCELLINDEX
       DO ii=1,grid%NCELLINDEX
        nelements=0
        DO ie=1,grid%CELLINDEX(ii)%NELEMENTS
         icell0=grid%CELLINDEX(ii)%INDEXES(0,ie) 
         nelements=nelements+inewcell%ILAST(icell0)
     -                      -inewcell%IFIRST(icell0)+1
        END DO
        CALL gridman_index_allocate(tria%CELLINDEX(ii),
     c                              grid%CELLINDEX(ii)%NINDEX,
     c                              nelements,ierr) 
        IF(ierr.GT.0) GOTO 1000
        n=0
        DO ie=1,grid%CELLINDEX(ii)%NELEMENTS
         icell0=grid%CELLINDEX(ii)%INDEXES(0,ie) 
         DO i=inewcell%IFIRST(icell0),inewcell%ILAST(icell0)
          icell=inewcell%IND(i)
          n=n+1  
          IF(n.GT.tria%CELLINDEX(ii)%NELEMENTS) GOTO 730
          tria%CELLINDEX(ii)%INDEXES(0,n)=icell
          tria%CELLINDEX(ii)%INDEXES(1:,n)=
     =                       grid%CELLINDEX(ii)%INDEXES(1:,ie) 
         END DO
        END DO
        IF(n.NE.tria%CELLINDEX(ii)%NELEMENTS) GOTO 730
        tria%CELLINDEX(ii)%DESCRIPTION=grid%CELLINDEX(ii)%DESCRIPTION
        tria%CELLINDEX(ii)%COLUMNS=grid%CELLINDEX(ii)%COLUMNS
       END DO
       CALL gridman_indlist_deallocate(inewcell,ierr0) 
      END IF !IF(ALLOCATED(IOLDCELL))

C     META DATA  
      tria%DESCRIPTION=trim(grid%DESCRIPTION)//
     /                 " >> triangulated by GRIDMAN_GRID2D_TRIANG"
C     (UNITS were written in GRIDMAN_GRID_TAKE)

      IF(gridman_check) THEN
       CALL gridman_grid2d_check(tria,res,ierr)
       IF(res.NE.0.OR.ierr.GT.0) THEN
        ierr=400
        WRITE(gridman_unit,*) "ERROR in GRIDMAN_GRID2D_TRIANG: ",
     w                        "wrong resulting grid"
       END IF
      END IF 

      IF(gridman_dbg) 
     f  WRITE(gridman_unit,*) "GRIDMAN_GRID2D_TRIANG finished"

      RETURN

  100 WRITE(gridman_unit,*) "ERROR in GRIDMAN_GRID2D_TRIANG: ",
     w                      "grid has no cells"
      ierr=100
      RETURN

  150 WRITE(gridman_unit,*) "ERROR in GRIDMAN_GRID2D_TRIANG: ",
     w            "cannon create temporar index list CHAINS"
      RETURN

  200 WRITE(gridman_unit,*) "ERROR in GRIDMAN_GRID2D_TRIANG: ",
     w                      "internal error"
      WRITE(gridman_unit,*) 
     w      " Cell has less than 3 points, ICELL ",icell 
      ierr=400
      CALL gridman_indlist_deallocate(chains,ierr0) 
      RETURN

  300 WRITE(gridman_unit,*) "ERROR in GRIDMAN_GRID2D_TRIANG: ",
     w                      "cannot allocate temporary array"
      ierr=200
      RETURN   
    
  400 ierr=400
      WRITE(gridman_unit,*) "ERROR in GRIDMAN_GRID2D_TRIANG: ",
     w                      "internal error"
      WRITE(gridman_unit,*) 
     w    " Wrong increment of the number of edges or cells"   
      WRITE(gridman_unit,*) " NCELLS, NCELLS0 ",tria%NCELLS,ncells0
      WRITE(gridman_unit,*) " NEDGES, NEDGES0 ",tria%NEDGES,nedges0
      CALL gridman_indlist_deallocate(chains,ierr0) 
      IF(ALLOCATED(ioldcell)) DEALLOCATE(ioldcell)
      RETURN

  500 ierr=400
      WRITE(gridman_unit,*) "ERROR in GRIDMAN_GRID2D_TRIANG: ",
     w          "mismatch between expected and actual grid size"
      WRITE(gridman_unit,*) "Expected: NCELLS, NEDGES ",ncells,nedges
      WRITE(gridman_unit,*) "Actual: NCELLS, NEDGES ",
     w                                tria%NCELLS,tria%NEDGES
      IF(ALLOCATED(ioldcell)) DEALLOCATE(ioldcell)

  600 WRITE(gridman_unit,*) "ERROR in GRIDMAN_GRID2D_TRIANG: ",
     w                      "cannon create edge index II ",ii
      IF(ALLOCATED(ioldcell)) DEALLOCATE(ioldcell)
      CALL gridman_indlist_deallocate(chains,ierr0) 
      RETURN

  700 WRITE(gridman_unit,*) "ERROR in GRIDMAN_GRID2D_TRIANG: ",
     w            "cannon create temporar index list INEWCELL"
      RETURN

  710 WRITE(gridman_unit,*) "ERROR in GRIDMAN_GRID2D_TRIANG: ",
     w            "internal error - wrong IOLDCELL"
      WRITE(gridman_unit,*) " ICELL0, ICELL, N, NCELLS ",
     w                        icell0,icell,n,ncells
      ierr=410 
      CALL gridman_indlist_deallocate(inewcell,ierr0) 
      IF(ALLOCATED(ioldcell)) DEALLOCATE(ioldcell)
      RETURN
  720 WRITE(gridman_unit,*) "ERROR in GRIDMAN_GRID2D_TRIANG: ",
     w                      "internal error - wrong IOLDCELL"
      WRITE(gridman_unit,*) " It must be N==NCELLS, N, NCELLS ",n,ncells
      ierr=420 
      CALL gridman_indlist_deallocate(inewcell,ierr0) 
      IF(ALLOCATED(ioldcell)) DEALLOCATE(ioldcell)
      RETURN
  725 WRITE(gridman_unit,*) "ERROR in GRIDMAN_GRID2D_TRIANG: ",
     w                      "cannot allocate TRIA%CELLINDEX"
      WRITE(gridman_unit,*) " NCELLINDEX ",grid%NCELLINDEX
      ierr=200  
      RETURN 
  730 WRITE(gridman_unit,*) "ERROR in GRIDMAN_GRID2D_TRIANG: ",
     w                      "internal error - unexpected cell index"
      WRITE(gridman_unit,*) " II, N, NELEMENTS ",ii,n,nelements 
      ierr=430 
      CALL gridman_indlist_deallocate(inewcell,ierr0) 
      IF(ALLOCATED(ioldcell)) DEALLOCATE(ioldcell)
      RETURN

 1000 CALL gridman_indlist_deallocate(chains,ierr0) 
      IF(ALLOCATED(ioldcell)) DEALLOCATE(ioldcell)
      WRITE(gridman_unit,*) "GRIDMAN_GRID2D_TRIANG is terminated"

      CONTAINS
C
C     ONE STEP OF TRIANGULATION 
C     EACH (NON-TRUANGLE) CELL IS DIVIDED INTO A TRIANGLE AND THE REST    
C     
C
      SUBROUTINE gridman_grid2d_triang_step(GRID,NCELLS_MAX,NEDGES_MAX,
     c                                      chains,ierr,ioldcell) 
      USE gridman,ONLY:gridman_dp
      USE gridman_lib,ONLY:gridman_grid_cells
      INTRINSIC tiny       
      TYPE(gridman_grid) :: GRID
      INTEGER(GRIDMAN_SP),INTENT(IN) :: NCELLS_MAX,NEDGES_MAX
      TYPE(gridman_indlist) :: CHAINS
      INTEGER,INTENT(OUT) :: IERR   
      INTEGER(GRIDMAN_SP),OPTIONAL :: IOLDCELL(ncells_max)

      INTEGER(GRIDMAN_SP) :: ICELL,IPOINT(3),IP1,IP,IP2,I,N,ID,IE,IEDGE,
     i                       ncells0,ncells,nedges,npmax,np,st
      INTEGER :: IERR0
      REAL(GRIDMAN_DP) :: VTVS,VTVT,STOT,MMCOS 
      REAL(GRIDMAN_DP) :: X(3),Y(3)
      TYPE(gridman_indlist) :: CELLS
C     ARRAYS WHICH STORE CANDIDADE NEW EDGES
      INTEGER(GRIDMAN_DP),ALLOCATABLE :: IEP(:,:)
      REAL(GRIDMAN_DP),ALLOCATABLE :: VTV(:),S(:),MAXCOS(:)

      IF(gridman_dbg) 
     f   WRITE(gridman_unit,*) "Starting GRIDMAN_GRID2D_TRIANG_STEP"

      ierr=0

C     DETERMINE THE MAXIMAL NUMBER OF POINTS (EDGES) OF THE CELLS 
      npmax=0
      DO icell=1,chains%N
       np=chains%ILAST(icell)-chains%IFIRST(icell)
       IF(np.GT.npmax) npmax=np 
      END DO

C     ALLOCATE TEMPORARY ARRAYS
      ALLOCATE(iep(3,npmax),vtv(npmax),s(npmax),maxcos(npmax),stat=st)
      IF(st.NE.0) THEN
       ierr=200
       WRITE(gridman_unit,*) "Erroor in GRIDMAN_GRID2D_TRIANG_STEP: ",  
     w                       "can not allocate temporary arrays"
       GOTO 500
      END IF

C     ARRAY OF EDGES FOR EACH CELL
      CALL gridman_grid_cells(cells,grid,ierr)   
      IF(ierr.GT.0) GOTO 500
      
      ncells=grid%NCELLS
      nedges=grid%NEDGES
      ncells0=ncells
      DO icell=1,ncells0 !CYCLE OVER CELLS
       n=chains%ILAST(icell)-chains%IFIRST(icell) !NUMBER OF EDGES
       IF(n.LT.4) cycle !SKIP TRIANGLES
       vtvs=0.
       stot=0.
       n=0
       ip1=chains%ILAST(icell)-1   
       DO ip=chains%IFIRST(icell),chains%ILAST(icell)-1
C       EDGE WHICH CONNECTS IP1 AND IP+1
C       CHECK FOR INTERSECTIONS WITH OTHER EDGES OF THIS CELL    
        ipoint(1)=chains%IND(ip1)
        ipoint(2)=chains%IND(ip)
        ipoint(3)=chains%IND(ip+1)  
        x(1:3)=grid%X(1,ipoint(1:3))
        y(1:3)=grid%X(2,ipoint(1:3))
        vtvt=(x(2)-x(1))*(y(3)-y(2))-(x(3)-x(2))*(y(2)-y(1))
        vtvs=vtvs+vtvt
        stot=stot+(x(3)+x(2))*(y(3)-y(2))
        IF(self_intersections(ipoint(1),ipoint(3),icell)) THEN
         ip1=ip
         cycle
        END IF            
C       ADD EDGE (IPOINT1,IPOINT3) INTO TEMPORARY ARRAY
        n=n+1  
        iep(1:3,n)=ipoint(1:3)
C       PERPENDICULAR DOT PRODUCT: TO FIND OUT IF THE VERTEX IS CONVEX
        vtv(n)=vtvt
C       AREAS OF THE RESULTING TRIANGLES
        s(n)=(x(3)+x(2))*(y(3)-y(2))+
     +       (x(2)+x(1))*(y(2)-y(1))+(x(1)+x(3))*(y(1)-y(3))
C       MAX COSINUS (SMALLEST ANGLE)
        maxcos(n)=maximum_cosine(x,y)
        IF(maxcos(n).LT.-1.0) THEN
         ierr=100
         GOTO 500
        END IF 
C
        ip1=ip
       END DO !DO IP=CHAINS%IFIRST(ICELL),CHAINS%ILAST(ICELL)-1
C          
C      TEMPORARY ARRAY CONTAINS THE CANDIDADE EDGES 
C      THE ONE WILL BE CHOSEN FOR WHICH
C      1. THE VORTEX IS CONCAVE
C      2. THE TRIANGLE AREA S IS SMALLER THAT 0.51*STOT
C
C      AMONG THOSE THE ONE WITH THE LARGEST SMALLEST ANGLE IS TAKEN
C
       mmcos=2.0
       stot=0.51*abs(stot)
       id=0
       DO i=1,n
        IF(( vtvs.LT.0.AND.vtv(i).GT.-tiny(vtvs)).OR. !CLOCKWISE DIRECTION
     f     ( vtvs.GT.0.AND.vtv(i).LT.tiny(vtvs)).OR.  !COUNTER-CLOCKWISE
     f       abs(s(i)).GT.stot) cycle !SKIP CONVEX CORNERS AND TOO LARGE CELLS
        IF(maxcos(i).LT.mmcos) THEN
         ipoint(1:3)=iep(1:3,i)
         id=id+1
         mmcos=maxcos(i)
        END IF   
       END DO              
       IF(id.EQ.0) THEN
        ierr=400
        WRITE(gridman_unit,*) "ERROR in GRIDMAN_GRID2_TRIANG_STEP: ",
     w                        "internal error - no new edges found"
        WRITE(gridman_unit,*) " ICELL, N ",icell,n
        WRITE(gridman_unit,*) " MAXCOS ",maxcos(1:n)
        WRITE(gridman_unit,*) " STOT, S ",stot,abs(s(1:n))
        GOTO 500
       END IF
C      ADD NEW EDGE AND NEW CELL
       ncells=ncells+1
       nedges=nedges+1
       IF(ncells.GT.ncells_max.OR.nedges.GT.nedges_max) THEN
        ierr=400
        WRITE(gridman_unit,*) "ERROR in GRIDMAN_GRID2_TRIANG_STEP: ",
     w                    "internal error - index exceeds the maximum"
        WRITE(gridman_unit,*) " ICELL ",icell
        WRITE(gridman_unit,*) " NCELLS, NCELLS_MAX ",ncells,ncells_max
        WRITE(gridman_unit,*) " NEDGES, NEDGES_MAX ",nedges,nedges_max
        GOTO 500
       END IF
       grid%CELLS(1,nedges)=icell  
       grid%CELLS(2,nedges)=ncells
       grid%POINTS(1,nedges)=ipoint(1)
       grid%POINTS(2,nedges)=ipoint(3)
C      CORRECT CELL INDEX FOR TWO OLD EDGES OF THE NEW TRIANGLE
       DO ie=cells%IFIRST(icell),cells%ILAST(icell)
        iedge=cells%IND(ie)
        IF(iedge.LT.1.OR.iedge.GT.nedges_max) THEN
         ierr=400
         WRITE(gridman_unit,*) "ERROR in GRIDMAN_GRID2_TRIANG_STEP: ",
     w                       "internal erorr - edge index out of range"
         WRITE(gridman_unit,*) " ICELL, IEDGE, NEDGES_MAX ",
     w                           icell, iedge, nedges_max
        END IF
        ip1=grid%POINTS(1,iedge)
        ip2=grid%POINTS(2,iedge)
        DO i=1,2
         IF((ip1.EQ.ipoint(i).AND.ip2.EQ.ipoint(i+1)).OR.
     f      (ip2.EQ.ipoint(i).AND.ip1.EQ.ipoint(i+1))) THEN  
          IF(grid%CELLS(1,iedge).EQ.icell ) THEN
             grid%CELLS(1,iedge)=ncells
             cycle
          ELSE IF(grid%CELLS(2,iedge).EQ.icell ) THEN
                  grid%CELLS(2,iedge)=ncells
                  cycle
          ELSE
           ierr=400
           WRITE(gridman_unit,*) 
     w               "ERROR in GRIDMAN_GRID2_TRIANG_STEP: ",
     w               "internal erorr - cell index mismatch"
           WRITE(gridman_unit,*) " ICELL, IEDGE ",icell,iedge
           WRITE(gridman_unit,*) " IPOINT1, IPOINT2 ",
     w                             ipoint(i),ipoint(i+1)
           WRITE(gridman_unit,*) " IP1, IP2 ",ip2,ip2
           WRITE(gridman_unit,*) " CELLS ",grid%CELLS(:,iedge)
           GOTO 500
          END IF        
         END IF
        END DO
       END DO 
       grid%NCELLS=ncells
       grid%NEDGES=nedges
       IF(PRESENT(ioldcell)) ioldcell(ncells)=ioldcell(icell)
      END DO !DO ICELL=1,NCELLS

      DEALLOCATE(iep,vtv,s,maxcos,stat=st)
      IF(st.NE.0) THEN 
       ierr=-200
       WRITE(gridman_unit,*) 
     w      "Warning from GRIDMAN_GRID2D_TRIANG_STEP ",
     w      "could not deallocate temporary arrays"
      END IF
      CALL gridman_indlist_deallocate(cells,ierr)

      IF(gridman_dbg) 
     f   WRITE(gridman_unit,*) "GRIDMAN_GRID2_TRIANG_STEP finished"

      RETURN 

 500  DEALLOCATE(iep,vtv,s,maxcos,stat=st)
      CALL gridman_indlist_deallocate(cells,ierr0)
      WRITE(gridman_unit,*) "GRIDMAN_GRID2_TRIANG_STEP terminated"
      RETURN

      END SUBROUTINE gridman_grid2d_triang_step

C
C     CALCULATE LARGEST COSINE OF THE TRIANGLES ANGLES:
C     THAT IS, COSINE OF THE SMALLEST ANGLE
C
      FUNCTION maximum_cosine(X,Y)
      USE gridman,ONLY:gridman_dp
      REAL(GRIDMAN_DP),INTENT(IN) :: X(3),Y(3)
      REAL(GRIDMAN_DP) :: MAXIMUM_COSINE
      REAL(GRIDMAN_DP) :: VX(3),VY(3),V(3),CS
      INTEGER :: I,I1 

      INTRINSIC sqrt

      maximum_cosine=-2.0
      DO i=1,3
       i1=i+1
       IF(i1.GT.3) i1=1
       vx(i)=x(i1)-x(i)
       vy(i)=y(i1)-y(i)
       v(i)=sqrt(vx(i)**2+vy(i)**2)
       IF(.NOT.v(i).GT.tiny(cs)) THEN
        WRITE(gridman_unit,*) "ERROR in GRIDMAN_GRID2_TRIANG: ",
     w                        "edge has zero length"   
        WRITE(gridman_unit,*) " X ",x

        WRITE(gridman_unit,*) " Y ",y

        RETURN
       END IF
      END DO
      DO i=1,3
       i1=i+1
       IF(i1.GT.3) i1=1
       cs=vx(i)*vx(i1)+vy(i)*vy(i1)
       cs=-cs/(v(i)*v(i1))  
       IF(cs.GT.maximum_cosine) maximum_cosine=cs
      END DO
      
      END FUNCTION maximum_cosine

C
      FUNCTION self_intersections(IPOINT10,IPOINT20,ICELL)
      USE gridman,ONLY:gridman_dp
      USE gridman_lib,ONLY:gridman_intersect2d
      INTEGER(GRIDMAN_SP),INTENT(IN) :: IPOINT10,IPOINT20,ICELL
      LOGICAL :: SELF_INTERSECTIONS
      REAL(GRIDMAN_DP) :: X01,X02,Y01,Y02,X1,X2,Y1,Y2,XI,YI
      INTEGER(GRIDMAN_SP) :: IP,IPOINT1,IPOINT2
      
      x01=grid%X(1,ipoint10)
      y01=grid%X(2,ipoint10)
      x02=grid%X(1,ipoint20)
      y02=grid%X(2,ipoint20)
      DO ip=chains%IFIRST(icell),chains%ILAST(icell)-1           
       ipoint1=chains%IND(ip)
       ipoint2=chains%IND(ip+1)
       IF(ipoint1.EQ.ipoint10.OR.ipoint1.EQ.ipoint20.OR.
     f    ipoint2.EQ.ipoint10.OR.ipoint2.EQ.ipoint20) cycle
        x1=grid%X(1,ipoint1)
        y1=grid%X(2,ipoint1)
        x2=grid%X(1,ipoint2)
        y2=grid%X(2,ipoint2)
        IF(gridman_intersect2d(x01,y01,x02,y02,
     f                         x1,y1,x2,y2,xi,yi)) THEN
         self_intersections=.true.
         RETURN
        END IF
      END DO
      self_intersections=.false.

      END FUNCTION self_intersections
 
      END SUBROUTINE gridman_grid2d_triang