tria.f

Go to the documentation of this file.

C> @file convert/tria.f
C> Converters of GRIDMAN_GRID object to and from EIRENE triangular 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> Convert EIRENE triangular grid into GRIDMAN_GRID grid object (type=GRID2D)
C**********************************************************************************************
C>
C> WARNING: if the object GRID already exists it will be overwritten!
C>
C> Correct XYTRIAN, NECKE, NCHBAR, NSEITE, ITRI, are produced by
C> GRIDMAN_TRIA_READ
C>
C> WARNING: if different ITRI(3:5) are defined for the same edge in two different
C> triangles, then one of the two values (the last to appear in the list) is taken.
C> Only non-zero edge tags are stored.
C
      SUBROUTINE gridman_tria2grid(GRID,NRKNOT,NTRII,XYTRIAN,NECKE,
     s                             nchbar,nseite,itri,ierr)
      USE gridman,ONLY:gridman_grid,gridman_dp,gridman_sp,
     u                 gridman_unit,gridman_dbg,gridman_check,
     u                 gridman_indlist
      USE gridman_lib,ONLY:gridman_grid_allocate,
     u                     gridman_grid_deallocate,
     u                     gridman_index_allocate,
     u                     gridman_index_deallocate,
     u                     gridman_grid2d_check
      IMPLICIT NONE
      INTRINSIC mod
C> Resulting grid object
      TYPE(gridman_grid) :: GRID           
C> Number of nodes 
      INTEGER(GRIDMAN_SP),INTENT(IN) :: NRKNOT
C> Number of triangles 
      INTEGER(GRIDMAN_SP),INTENT(IN) :: NTRII  
C> Coordinates of the nodes
C>
C> XYTRIAN(1,:) is the X coordinate, XYTRIAN(2,:) is the Y coordinate
      REAL(GRIDMAN_DP),INTENT(IN) :: XYTRIAN(2,nrknot) 
C> Indices of nodes for each triangle              
      INTEGER(GRIDMAN_SP),INTENT(IN) :: NECKE(3,ntrii)  
C> Indices of neigbour cell for each triangle
      INTEGER(GRIDMAN_SP),INTENT(IN) :: NCHBAR(3,ntrii)  
C> Indices of sides of neighbours for each triangle    
      INTEGER(GRIDMAN_SP),INTENT(IN) :: NSEITE(3,ntrii) 
C> Tags from fort.35
C>
C> ITRI(1,:) and ITRI(2,:) are IXTRI and IYTRI.
C> ITRI(3:5,:) are tags of each side of the triange.
      INTEGER(GRIDMAN_SP),INTENT(IN) :: ITRI(5,ntrii)    
C> Error code
      INTEGER,INTENT(OUT) :: IERR 

      INTEGER(GRIDMAN_SP) :: IND(ntrii*3)  
      INTEGER(GRIDMAN_SP) :: ICELL,ICELL1,IS,IS1,IEDGE,IEDGE1,IP,IPOINT,
     i                       ie,nedges,nelements
      INTEGER :: IERR0,RES
      INTEGER(GRIDMAN_SP),ALLOCATABLE :: IEDGES(:)
      LOGICAL :: LWITRI
     
      IF(gridman_dbg) 
     f   WRITE(gridman_unit,*) "Starting GRIDMAN_TRIA2GRID"

      ierr=0

C     CHECKING INPUT
      IF(nrknot.LE.0) THEN
       ierr=100
       WRITE(gridman_unit,*) "ERROR in GRIDMAN_TRIA2GRID: ",
     w                       .LE."NRKNOT0"
       RETURN 
      END IF
      IF(ntrii.LE.0) THEN
       ierr=100
       WRITE(gridman_unit,*) "ERROR in GRIDMAN_TRIA2GRID: ",
     w                       .LE."NTRII0"
       RETURN 
      END IF

C     CALCULATE THE NUMBER OF EDGES AND AUXILIARY TABLE 
C     TO CONVERT TRIANGLE-CENTERED EDGE INDICES = (ICELL-1)*3+IS
C     TO THE EDGE INDICES IN THE RESULTING GRID 
      nedges=0  
      ind=0
      DO icell=1,ntrii
       DO is=1,3
        iedge=(icell-1)*3+is
        icell1=nchbar(is,icell)  
        is1=nseite(is,icell)
        IF(icell1.GT.0) THEN
C        EDGE BETWEEN TWO TRIANGLES 
         IF(is1.LT.1) THEN
          ierr=100
          WRITE(gridman_unit,*) "ERROR in GRIDMAN_TRIA2GRID: ",
     w                          "no open side, but ISIDE<=0"
          WRITE(gridman_unit,*) " ITRIA, ISIDE ",icell,is
          RETURN         
         END IF
         iedge1=(icell1-1)*3+is1
         IF(ind(iedge1).EQ.0) THEN
C         ADD CLOSED SIDE ONLY IF IT WAS NOT ALREADY ADDED
          nedges=nedges+1
          ind(iedge)=nedges
         ELSE
C         OTHERWISE DUBLICATE THE INDEX
          ind(iedge)=ind(iedge1)
         END IF
        ELSE
C        ADD OPEN SIDE OF TRIANGLE
         nedges=nedges+1
         ind(iedge)=nedges
        END IF
       END DO !DO IS=1,3
      END DO !DO ICELL=1,NTRII  
  
C     ALLOCATE GRID OBJECT
      CALL gridman_grid_allocate(grid,2,nedges,nrknot,ntrii,ierr,1,1)  
      IF(ierr.GT.0) GOTO 100
      CALL gridman_index_allocate(grid%CELLINDEX(1),2,ntrii,ierr)
      IF(ierr.NE.0) THEN
       WRITE(gridman_unit,*) "ERROR in GRIDMAN_TRIA2GRID: ",
     w                       "cannot allocate cell index"
       GOTO 200
      END IF
C      ... edge index is allocated later


C     DEFINE ARRAY GRID%CELLS
      grid%CELLS=0
      DO icell=1,ntrii
       DO is=1,3      
        iedge=ind((icell-1)*3+is)
        IF(iedge.LT.1) THEN
         WRITE(gridman_unit,*) "ERROR in GRIDMAN_TRIA2GRID: ",
     w                         "internal error IND<0"
         WRITE(gridman_unit,*) " ITRIA, ISIDE, IEDGE ",icell,is,iedge
         ierr=400
         GOTO 200
        END IF
        IF(grid%CELLS(1,iedge).EQ.0) THEN
         grid%CELLS(1,iedge)=icell
        ELSE IF(grid%CELLS(2,iedge).EQ.0) THEN  
         grid%CELLS(2,iedge)=icell
        ELSE
         WRITE(gridman_unit,*) "ERROR in GRIDMAN_TRIA2GRID: ",
     w                         "wrong triangle data"
         WRITE(gridman_unit,*) "Edge belongs to more that 2 triangles"
         WRITE(gridman_unit,*) " ITRIA, ISIDE, IEDGE ",icell,is,iedge
         WRITE(gridman_unit,*) " CELLS ",grid%CELLS(:,iedge),icell
         ierr=100
         GOTO 200
        END IF
       END DO
      END DO 
C     DEFINE ARRAY GRID%POINTS
      grid%POINTS=0
      ALLOCATE(iedges(nedges),stat=is)
      IF(is.NE.0) THEN 
       WRITE(gridman_unit,*) "ERROR in GRIDMAN_TRIA2GRID: ",
     w                       "cannot allocate temporary array"
       WRITE(gridman_unit,*) "NEDGES ",nedges
       GOTO 200
      END IF
      iedges=0
      nelements=0
      DO icell=1,ntrii
       DO is=1,3
        iedge=ind((icell-1)*3+is)
        is1=is
        DO ip=1,2
         ipoint=necke(is1,icell)  
         IF(grid%POINTS(ip,iedge).EQ.0) THEN
          grid%POINTS(ip,iedge)=ipoint  
         ELSE
          IF(ipoint.NE.grid%POINTS(1,iedge).AND.
     f       ipoint.NE.grid%POINTS(2,iedge)) THEN
           WRITE(gridman_unit,*) "ERROR in GRIDMAN_TRIA2GRID: ",
     w          "point indices do not match previous definition"
           WRITE(gridman_unit,*) " ITRIA, ISIDE, IEDGE ",icell,is1,iedge
           WRITE(gridman_unit,*) " NECKE, POINTS ",
     w                             ipoint,grid%POINTS(:,iedge)                      
           ierr=100
           DEALLOCATE(iedges)
           GOTO 200
          END IF 
         END IF        
         is1=is1+1
         IF(is1.GT.3) is1=1
        END DO
        IF(itri(2+is,icell).NE.0.AND.iedges(iedge).EQ.0) THEN
         nelements=nelements+1
         iedges(iedge)=nelements
        END IF    
       END DO !DO IS=1,3
       grid%CELLINDEX(1)%INDEXES(0,icell)=icell
       grid%CELLINDEX(1)%INDEXES(1,icell)=itri(1,icell) !IXTRI
       grid%CELLINDEX(1)%INDEXES(2,icell)=itri(2,icell) !IYTRI
      END DO !DO ICELL=1,NTRII

C     EDGE INDEX
      CALL gridman_index_allocate(grid%EDGEINDEX(1),1,nelements,ierr)
      IF(ierr.NE.0) THEN
       WRITE(gridman_unit,*) "ERROR in GRIDMAN_TRIA2GRID: ",
     w                       "cannot allocate edge index"
       GOTO 200
      END IF

      lwitri=.false.
      iedges=0
      ie=0
      DO icell=1,ntrii
       DO is=1,3
        iedge=ind((icell-1)*3+is)
        IF(itri(2+is,icell).NE.0) THEN
         IF(iedges(iedge).EQ.0) THEN
          ie=ie+1
          grid%EDGEINDEX(1)%INDEXES(0,ie)=iedge
          grid%EDGEINDEX(1)%INDEXES(1,ie)=itri(2+is,icell)
          iedges(iedge)=ie
         ELSE        
          IF(grid%EDGEINDEX(1)%INDEXES(1,iedges(iedge)).NE. 
     f       itri(2+is,icell)) lwitri=.true.
         END IF
        END IF
       END DO 
      END DO 

      DEALLOCATE(iedges,stat=is)
      IF(is.NE.0) THEN
       WRITE(gridman_unit,*) "WARNING from GRIDMAN_TRIA2GRID: ",
     w               "cannot deallocate temporary array IEDGES"
       ierr=-200
      END IF
     
      IF(lwitri) THEN
       WRITE(gridman_unit,*) "WARNING from GRIDMAN_TRIA2GRID: "
       WRITE(gridman_unit,*) 
     w      "Some edge tags defined for the same side ",
     w      "in different triangles do not match"
       ierr=-100
      END IF

C     CORDINATES
      grid%X=xytrian 

C     UNITS
      grid%UNIT2SI=1e-2_gridman_dp        !CM TO M
      grid%UNITS='CENTIMETER'
C     DEFAULT DESCRIPTION
      grid%DESCRIPTION='EIRENE triangular grid'      
      grid%CELLINDEX(1)%DESCRIPTION='IX, IY of the plasma grid'      
      grid%CELLINDEX(1)%COLUMNS(1)='TRIA_IX'
      grid%CELLINDEX(1)%COLUMNS(2)='TRIA_IY'
      grid%EDGEINDEX(1)%DESCRIPTION='Edge tags of TRIA grid'      
      grid%EDGEINDEX(1)%COLUMNS(1)='TRIA_ISRF'

      IF(gridman_check) THEN
       CALL gridman_grid2d_check(grid,res,ierr0)
       IF(res.NE.0.OR.ierr0.GT.0) THEN  
        ierr=100
        WRITE(gridman_unit,*) "ERROR in GRIDMAN_TRIA2GRID: ", 
     w                        "the resulting grid is incorrect" 
        GOTO 200
       END IF
      END IF  

      IF(gridman_dbg) 
     f  WRITE(gridman_unit,*) "GRIDMAN_TRIA2GRID finished"
      RETURN

 100  WRITE(gridman_unit,*) "GRIDMAN_TRIA2GRID terminated"
 200  CALL gridman_grid_deallocate(grid,ierr0)
      RETURN

      END SUBROUTINE gridman_tria2grid

C**********************************************************************************************
C> Convert GRIDMAN_GRID grid object into EIRENE triangular grid
C>
C> All grid cells must be triangles - this subroutine does not 
C> do triangulation. For triangulation see \ref GRIDMAN_GRID2D_TRIANG
C>
C> XYTRIAN, NECKE, NCHBAR, NSEITE, ITRI
C> are produced in exactly same format as by GRIDMAN_TRIA_READ.
C> Arrays have to be allocated in advance.
C>  
C> ITRI is taken from the 1st cell and edge indices.
C> ITRI is not changed if this information is missing. 
C**********************************************************************************************
      SUBROUTINE gridman_grid2tria(GRID,XYTRIAN,NECKE,
     s                             nchbar,nseite,itri,ierr)
      USE gridman,ONLY:gridman_grid,gridman_indlist,
     u                 gridman_dp,gridman_sp,
     u                 gridman_unit,gridman_dbg,gridman_check
      USE gridman_lib,ONLY:gridman_grid2d_check,gridman_grid_cells,
     u            gridman_grid2d_chains,gridman_indlist_deallocate
      IMPLICIT NONE
      INTRINSIC mod,min
C> 2D grid object to be converted into triangular grid 
      TYPE(gridman_grid) :: GRID                 
C> Coordinates of the nodes
      REAL(GRIDMAN_DP),INTENT(OUT) :: XYTRIAN(2,grid%npoints) 
C> Indices of nodes for each triangle                  
      INTEGER(GRIDMAN_SP),INTENT(OUT) :: NECKE(3,grid%ncells) 
C> Indices of neigbour cell for each triangle
      INTEGER(GRIDMAN_SP),INTENT(OUT) :: NCHBAR(3,grid%ncells) 
C> Indices of sides of neighbours for each triangle    
      INTEGER(GRIDMAN_SP),INTENT(OUT) :: NSEITE(3,grid%ncells)
C> Tags from fort.35
      INTEGER(GRIDMAN_SP),INTENT(OUT) :: ITRI(5,grid%ncells)  
C> Error code
      INTEGER,INTENT(OUT) :: IERR 

      TYPE(gridman_indlist) :: CELLS,POINTS  
      INTEGER(GRIDMAN_SP) :: ICELL,N,IPOINT,I,IE,IEDGE,ICELL2,
     i                       ipoint1,ipoint2,iside1,iside2,itmp
      INTEGER(GRIDMAN_SP) :: INDCELL(grid%nedges),INDSIDE(grid%nedges)
      INTEGER :: RES,IERR0,L
      REAL(GRIDMAN_DP) :: S,X(3),Y(3)

      IF(gridman_dbg) WRITE(gridman_unit,*) "Starting GRIDMAN_GRID2TRIA"

      ierr=0

C     THE GRID MUST BE OF 2D TYPE
      IF(grid%TYPE.NE.2.OR.grid%PDIM.NE.2.OR.grid%EDIM.NE.2) THEN
       ierr=100
       WRITE(gridman_unit,*) "GRIDMAN_GRID2TRIA: ",
     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(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_GRID2TRIA: ", 
     w                        "the grid is incorrect" 
        RETURN
       END IF
      END IF  
 
      CALL gridman_grid_cells(cells,grid,ierr)
      IF(ierr.NE.0) GOTO 100

C     ALL CELLS MUST CONSIST OF EXACTLY 3 EDGES   
      DO icell=1,grid%NCELLS
       n=cells%ILAST(icell)-cells%IFIRST(icell)
       IF(n.NE.2) THEN
        WRITE(gridman_unit,*) "ERROR in GRIDMAN_GRID2TRIA: ", 
     w                    "cell is not triangle "
        WRITE(gridman_unit,*) " ICELL, N" ,icell,n+1 
        ierr=100
        GOTO 200   
       END IF
      END DO

C     COORDINATES OF THE CORNERS
      xytrian=grid%X*grid%UNIT2SI*100_gridman_dp !CONVERT UNITS TO CM

C     EACH TRIANGLE IS DEFINED IN FORT.34 BY ITS 3 POINTS:
C         POINT1 POINT2 POINT3
C     EIRENE ASSUMES THAT SIDE 1 OF THE TRIANGLE LIES BETWEEN POINT1 AND POINT2
C     SIDE 2 LIES BETWEEN POINT2 AND POINT3 AND    
C     SIDE 2 LIES BETWEEN POINT3 AND POINT1
C     THOSE SIDE INDECES ARE ASSUMED IN THE LIST OF NEIGHBOURS IN FORT.35
C     THEN REFERING TO THE SIDE OF THE NEIGHBOURING TRIANGLES
C     ALSO, EIRENE ASSUMES THAT POINT1 POINT2 POINT3 ARE ARRANGED IN
C     THE COUNTERCLOCKWISE ORDER     
      CALL gridman_grid2d_chains(grid,points,ierr)
      IF(ierr.NE.0) GOTO 100

C     CORNERS OF EACH TRIANGLE
      DO icell=1,grid%NCELLS
       i=0
       DO ie=points%IFIRST(icell),points%ILAST(icell)-1
        ipoint=points%IND(ie)
        i=i+1
        IF(i.GT.3) THEN
         WRITE(gridman_unit,*) "ERROR in GRIDMAN_GRID2TRIA: ", 
     w        "internal erorr - cell has more than 3 points "
         WRITE(gridman_unit,*) " ICELL" ,icell 
         ierr=400
         GOTO 200   
        END IF
        necke(i,icell)=ipoint       
       END DO
      END DO

C     CHANGE ORDER OF CORNERS IF THEY ARE NOT IN COUNTERCLOCKWISE ORDER
      DO icell=1,grid%NCELLS       
       DO i=1,3
        ipoint=necke(i,icell)
        IF(ipoint.LT.1.OR.ipoint.GT.grid%NPOINTS) THEN
          WRITE(gridman_unit,*) "ERROR in GRIDMAN_GRID2TRIA: ", 
     w    "internal error - triangle corner index is out of range"
          WRITE(gridman_unit,*) " ICELL, ISIDE, IPOINT, NPOINTS " ,
     w                            icell,i,ipoint,grid%NPOINTS
          ierr=400
          GOTO 200  
        END IF
        x(i)=grid%X(1,ipoint)
        y(i)=grid%X(2,ipoint)
       END DO    
       s=(x(1)*y(2)-x(2)*y(1))+(x(2)*y(3)-x(3)*y(2))+
     +   (x(3)*y(1)-x(1)*y(3))
       IF(s.LT.0.) THEN
C       1 2 3 -> 1 3 2, EXCHANGE 2 and 3
        itmp=necke(3,icell)
        necke(3,icell)=necke(2,icell)
        necke(2,icell)=itmp
       END IF 
      END DO

C     TABLE OF NEIGHBOURS
      indcell=0
      indside=0
      DO icell=1,grid%NCELLS
       DO ie=cells%IFIRST(icell),cells%ILAST(icell)
        iedge=cells%IND(ie)
        ipoint1=grid%POINTS(1,iedge)
        ipoint2=grid%POINTS(2,iedge)
        iside1=find_side(icell,ipoint1,ipoint2)
        IF(iside1.LT.1) THEN
         WRITE(gridman_unit,*) "ERROR in GRIDMAN_GRID2TRIA: ", 
     w           "internal error - cannot find triangle side"
         WRITE(gridman_unit,*) " ICELL, IEDGE " ,icell,iedge 
         WRITE(gridman_unit,*) " IPOINT1, IPOINT2 ",ipoint1,ipoint2 
         ierr=400
         GOTO 200   
        END IF
        icell2=grid%CELLS(1,iedge)
        IF(icell2.EQ.icell) icell2=grid%CELLS(2,iedge)
        IF(icell2.LT.1) THEN
         icell2=0  
         iside2=0
        ELSE
         iside2=find_side(icell2,ipoint1,ipoint2)
         IF(iside2.LT.1) THEN
          WRITE(gridman_unit,*) "ERROR in GRIDMAN_GRID2TRIA: ", 
     w            "internal error - cannot find triangle side"
          WRITE(gridman_unit,*) " ICELL2, IEDGE " ,icell2,iedge 
          WRITE(gridman_unit,*) " IPOINT1, IPOINT2 ",ipoint1,ipoint2 
          ierr=400
          GOTO 200           
         END IF
        END IF
        nchbar(iside1,icell)=icell2
        nseite(iside1,icell)=iside2
        indcell(iedge)=icell
        indside(iedge)=iside1
       END DO !DO IE=CELLS%IFIRST(ICELL),CELLS%ILAST(ICELL)
      END DO

      IF(grid%NCELLINDEX.GT.0) THEN
       l=min(grid%CELLINDEX(1)%NINDEX,2)
       DO ie=1,grid%CELLINDEX(1)%NELEMENTS
        icell=grid%CELLINDEX(1)%INDEXES(0,ie)
        IF(icell.GT.0.AND.icell.LE.grid%NCELLS) THEN
         itri(1:l,icell)=grid%CELLINDEX(1)%INDEXES(1:l,ie) !IX,IY
        END IF
       END DO
      END IF
      IF(grid%NEDGEINDEX.GT.0) THEN
       DO ie=1,grid%EDGEINDEX(1)%NELEMENTS
        iedge=grid%EDGEINDEX(1)%INDEXES(0,ie)   
        IF(iedge.GT.0.AND.iedge.LE.grid%NEDGES) THEN
         icell=indcell(iedge)
         iside1=indside(iedge)
         IF(icell.GT.0.AND.iside1.GT.0) THEN
          itri(2+iside1,icell)=grid%EDGEINDEX(1)%INDEXES(1,ie)   
          icell2=nchbar(iside1,icell)
          iside2=nseite(iside1,icell)
          IF(icell2.GT.0.AND.iside2.GT.0)
     f      itri(2+iside2,icell2)=itri(2+iside1,icell)
         END IF      
        END IF
       END DO
      END IF
 
      CALL gridman_indlist_deallocate(cells,ierr)
      CALL gridman_indlist_deallocate(points,ierr)

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

      RETURN

 100  WRITE(gridman_unit,*) "GRIDMAN_GRID2TRIA terminated" 
 200  CALL gridman_indlist_deallocate(cells,ierr0)
      CALL gridman_indlist_deallocate(points,ierr0)

      CONTAINS

C     FIND WHICH SIDE OF TRIANGLE ICELL 
C     CORRESPONDS TO POINTS IPOINT1, IPOINT2
      FUNCTION find_side(ICELL,IPOINT1,IPOINT2)
       INTEGER(GRIDMAN_SP),INTENT(IN) :: ICELL,IPOINT1,IPOINT2
       INTEGER(GRIDMAN_SP) :: FIND_SIDE
       INTEGER(GRIDMAN_SP) :: IP1,IP2,IP3

C     EIRENE ASSUMES THAT SIDE 1 OF THE TRIANGLE LIES BETWEEN IP1 AND IP2
C     SIDE 2 LIES BETWEEN IP2 AND IP3 AND SIDE 3 LIES BETWEEN IP3 AND IP1

       ip1=necke(1,icell)
       ip2=necke(2,icell)
       ip3=necke(3,icell)
       IF((ipoint1.EQ.ip1.AND.ipoint2.EQ.ip2).OR.
     f    (ipoint1.EQ.ip2.AND.ipoint2.EQ.ip1)) THEN
        find_side=1
       ELSE IF((ipoint1.EQ.ip2.AND.ipoint2.EQ.ip3).OR.
     f         (ipoint1.EQ.ip3.AND.ipoint2.EQ.ip2)) THEN
        find_side=2   
       ELSE IF((ipoint1.EQ.ip3.AND.ipoint2.EQ.ip1).OR.
     f         (ipoint1.EQ.ip1.AND.ipoint2.EQ.ip3)) THEN
        find_side=3
       ELSE
        find_side=0
       END IF
      END FUNCTION find_side

      END SUBROUTINE gridman_grid2tria


C**********************************************************************************************
C> Read EIRENE triangular grid from fort.33-35, returns GRID object
C>
C> WARNING: if GRID object elready exists it will be overwritten,
C> array ITRI will be re-allocated as well
C>
C> This subroutine is a shell for GRIDMAN_TRIA_READ_ARRAY
C>
C> GRIDMAN_TRIA_READ_ARRAY and GRIDMAN_TRIA_READ_GRID
C> are combined in the interface GRIDMAN_TRIA_READ
C>
C**********************************************************************************************
      SUBROUTINE gridman_tria_read_grid(GRID,FNAMES,IERR)
       USE gridman,ONLY:gridman_grid,gridman_dp,gridman_sp,
     u                  gridman_dbg,gridman_unit
       USE gridman_lib,ONLY:gridman_tria_read_array,gridman_tria2grid
       IMPLICIT NONE
C> Resulting grid   
       TYPE(gridman_grid) :: GRID  
C> Strungs with names of three input files
C>
C> If FNAMES has zero length than standard names are used: fort.33, fort.34, fort.35
       CHARACTER(*) :: FNAMES(3)     
C> Error code
      INTEGER,INTENT(OUT) :: IERR 

      INTEGER(GRIDMAN_SP) :: NRKNOT,NTRII
      INTEGER :: IS
      INTEGER(GRIDMAN_SP),ALLOCATABLE :: NECKE(:,:),
     i            nchbar(:,:),nseite(:,:),itri(:,:)
      REAL(GRIDMAN_DP),ALLOCATABLE :: XYTRIAN(:,:) 

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

      ierr=0

      CALL gridman_tria_read_array(fnames,nrknot,ntrii,xytrian,
     s                             necke,nchbar,nseite,itri,ierr)
      IF(ierr.NE.0) GOTO 100

      CALL gridman_tria2grid(grid,nrknot,ntrii,xytrian,necke,
     s                       nchbar,nseite,itri,ierr)

      IF(ALLOCATED(xytrian)) DEALLOCATE(xytrian,stat=is)
      IF(ALLOCATED(necke)) DEALLOCATE(necke,stat=is)
      IF(ALLOCATED(nchbar)) DEALLOCATE(nchbar,stat=is)
      IF(ALLOCATED(nseite)) DEALLOCATE(nseite,stat=is)
      IF(ALLOCATED(itri)) DEALLOCATE(itri,stat=is)

      IF(gridman_dbg) 
     f   WRITE(gridman_unit,*) "GRIDMAN_TRIA_READ_GRID finished"
      RETURN 

 100  WRITE(gridman_unit,*) "GRIDMAN_TRIA_READ_GRID terminated"
      RETURN

      END SUBROUTINE gridman_tria_read_grid

C**********************************************************************************************
C> Write EIRENE triangular grid - defined as GRID object - into fort.33-35 
C>
C> Output of GRIDMAN_TRIA_READ is input of GRIDMAN_TRIA_WRITE
C>
C> GRIDMAN_TRIA_WRITE_ARRAY and GRIDMAN_TRIA_WRITE_GRID
C> are combined in the interface GRIDMAN_TRIA_WRITE
C**********************************************************************************************
      SUBROUTINE gridman_tria_write_grid(GRID,FNAMES,IERR)
       USE gridman,ONLY:gridman_grid,gridman_dp,gridman_sp,
     u                  gridman_dbg,gridman_unit
       USE gridman_lib,ONLY:gridman_tria_write_array,gridman_grid2tria
       IMPLICIT NONE
C> Grid object to be stored
       TYPE(gridman_grid) :: GRID     
C> Strungs with names of three input files
C>
C> If FNAMES has zero length than standard names are used: fort.33, fort.34, fort.35
       CHARACTER(*) :: FNAMES(3)  
C> Error code
       INTEGER,INTENT(OUT) :: IERR 

       INTEGER(GRIDMAN_SP) :: NRKNOT,NTRII
       INTEGER :: IS
       INTEGER(GRIDMAN_SP),ALLOCATABLE :: NECKE(:,:),NCHBAR(:,:),
     i                                    nseite(:,:),itri(:,:)
       REAL(GRIDMAN_DP),ALLOCATABLE :: XYTRIAN(:,:)

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

       ntrii=grid%NCELLS
       nrknot=grid%NPOINTS
       ALLOCATE(xytrian(2,nrknot),necke(3,ntrii), 
     a          nchbar(3,ntrii),nseite(3,ntrii),
     a          itri(5,ntrii),stat=is)
       IF(is.NE.0) THEN
        WRITE(gridman_unit,*) "ERROR in GRIDMAN_TRIA_WRITE_GRID",
     w                        "can not allocate temporary arrays"
        WRITE(gridman_unit,*) "  NTRII, NRKNOT ", ntrii,nrknot
        ierr=200
        CALL deallocate_all()
        RETURN
       END IF
       itri=0
       CALL gridman_grid2tria(grid,xytrian,necke,
     s                        nchbar,nseite,itri,ierr)
       IF(ierr.NE.0) THEN
        CALL deallocate_all()
        WRITE(gridman_unit,*) "GRIDMAN_TRIA_WRITE_GRID terminated"
        RETURN
       END IF
       CALL gridman_tria_write_array(fnames,nrknot,ntrii,xytrian,
     s                               necke,nchbar,nseite,itri,ierr)
       CALL deallocate_all()

       IF(ierr.NE.0) THEN
        WRITE(gridman_unit,*) "GRIDMAN_TRIA_WRITE_GRID terminated"
        RETURN
       END IF

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

       CONTAINS 
C
       SUBROUTINE deallocate_all
        IF(ALLOCATED(necke)) DEALLOCATE(necke)
        IF(ALLOCATED(nchbar)) DEALLOCATE(nchbar)
        IF(ALLOCATED(nseite)) DEALLOCATE(nseite)
        IF(ALLOCATED(xytrian)) DEALLOCATE(xytrian)
        IF(ALLOCATED(itri)) DEALLOCATE(itri)
       END SUBROUTINE deallocate_all

      END SUBROUTINE gridman_tria_write_grid