template.f

Go to the documentation of this file.

C> @file convert/template.f
C> Converters of GRIDMAN_GRID object to and from simple template (DG template)
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 data read from a file in simple template format (DG template) into GRIDMAN_GRID object
C**********************************************************************************************
C>
C> WARNING: object GRID will be overwritten if already exists 
C>
C> The code tries to figure out if the polygone is closed or not. For
C> closed polygones last point of the chain is moved to the first point.
C> WARNING: single points are not taken into the grid object
C>
C> Edge index EDGEINDEX(1) contains index of polygons for each edge
      SUBROUTINE gridman_template2grid(GRID,M,N,L,X,Y,IERR)
      USE gridman,ONLY:gridman_grid,gridman_sp,gridman_dp,
     u                 gridman_dbg,gridman_unit,gridman_check,
     u                 gridman_tol
      USE gridman_lib,ONLY:gridman_grid_allocate,gridman_index_allocate,
     u                     gridman_grid_deallocate,gridman_grid2d_check
      IMPLICIT NONE
      INTRINSIC abs
C> Resulting grid object 
      TYPE(gridman_grid) :: GRID
C> Number of polygons 
      INTEGER(GRIDMAN_SP),INTENT(IN) :: M
C> Number of vertices in each polygon
      INTEGER(GRIDMAN_SP),INTENT(IN) :: N(m)
C> Total number of points, must be L=SUM(N)
      INTEGER(GRIDMAN_SP),INTENT(IN) :: L
C> X coordinates of the vertices
      REAL(GRIDMAN_DP),INTENT(IN) :: X(l)
C> Y coordinates of the vertices 
      REAL(GRIDMAN_DP),INTENT(IN) :: Y(l)
C> Error code
      INTEGER,INTENT(OUT) :: IERR

      INTEGER(GRIDMAN_SP) :: L1,IP,IV,IL,IPOINT,IPOINT0,
     i                       iedge,nedges,npoints
      INTEGER :: IERR0,RES
      REAL(GRIDMAN_DP) :: X1,Y1,X2,Y2
       
      IF(gridman_dbg) 
     w   WRITE(gridman_unit,*) "Starting GRIDMAN_TEMPLATE2GRID"

      ierr=0 

C     CHECK INPUT
      IF(m.LT.1) THEN
       WRITE(gridman_unit,*) "ERROR in GRIDMAN_TEMPLATE2GRID: ",
     w                       "incorrect number of polygons, M ",m
       ierr=100
       RETURN
      END IF

      l1=sum(n)
      IF(l1.NE.l) THEN
       WRITE(gridman_unit,*) "ERROR in GRIDMAN_TEMPLATE2GRID: ",
     w                       "incorrect number of points - L=/=sum(N)"
       WRITE(gridman_unit,*) "L, SUM(N) ",l,l1
       ierr=100
       RETURN
      END IF

C     FIND, WHICH POLYGONS ARE CLOSED AND
C     CALCULATE NUMBER OF POINTS AND EDGES
      il=0
      npoints=0
      nedges=0
      DO ip=1,m
       IF(n(ip).GT.1) THEN !SINGLE POINTS ARE SKIPPED
        x1=x(il+1)
        y1=y(il+1)
        x2=x(il+n(ip))
        y2=y(il+n(ip))
        IF(abs(x1-x2).GT.gridman_tol*(abs(x1)+abs(x2)).OR.
     f     abs(y1-y2).GT.gridman_tol*(abs(y1)+abs(y2))) THEN
C        OPEN POLYGON
         nedges=nedges+n(ip)-1 
         npoints=npoints+n(ip) 
        ELSE
C        CLOSED POLYGON  
         IF(n(ip).LT.4) THEN
          WRITE(gridman_unit,*) "ERROR in GRIDMAN_TEMPLATE2GRID: ",
     w                 "close polygon must have at least 4 points"
          WRITE(gridman_unit,*) " IP, N ",ip,n(ip)
          ierr=100
          RETURN
         END IF
         nedges=nedges+n(ip)-1 
         npoints=npoints+n(ip)-1       
        END IF
       END IF !IF(N(IP).GT.1)
       il=il+n(ip)
      END DO      

C     ALLOCATE GRID OBJECT
      CALL gridman_grid_allocate(grid,2,nedges,npoints,
     c                           0_gridman_sp,ierr,1)
      IF(ierr.NE.0) THEN
       WRITE(gridman_unit,*) "GRIDMAN_TEMPLATE2GRID terminated"
       RETURN
      END IF
      CALL gridman_index_allocate(grid%EDGEINDEX(1),1,nedges,ierr)
      IF(ierr.NE.0) THEN
       WRITE(gridman_unit,*) "ERROR in GRIDMAN_TEMPLATE2GRID: ",
     w                       "cannot allocate index"
       CALL gridman_grid_deallocate(grid,ierr0)
       RETURN
      END IF

C     DEFINE GRID OBJECT
      il=0
      iedge=0
      ipoint=0
      DO ip=1,m
       IF(n(ip).GT.1) THEN !SINGLE POINTS ARE SKIPPED
        x1=x(il+1)
        y1=y(il+1)
        x2=x(il+n(ip))
        y2=y(il+n(ip))
        IF(abs(x1-x2).GT.gridman_tol*(abs(x1)+abs(x2)).OR.
     f     abs(y1-y2).GT.gridman_tol*(abs(y1)+abs(y2))) THEN
C        OPEN POLYGON
         DO iv=1,n(ip)-1
          iedge=iedge+1 
          IF(iedge.GT.nedges) GOTO 200
          ipoint=ipoint+1    
          IF(ipoint.GT.npoints) GOTO 300
          grid%X(1,ipoint)=x(il+iv)
          grid%X(2,ipoint)=y(il+iv)
          grid%POINTS(1,iedge)=ipoint
          grid%POINTS(2,iedge)=ipoint+1
          grid%EDGEINDEX(1)%INDEXES(0,iedge)=iedge
          grid%EDGEINDEX(1)%INDEXES(1,iedge)=ip
         END DO
         ipoint=ipoint+1    
         IF(ipoint.GT.npoints) GOTO 300
         grid%X(1,ipoint)=x(il+n(ip))
         grid%X(2,ipoint)=y(il+n(ip))
        ELSE
C        CLOSED POLYGON     
         ipoint0=ipoint+1  
         DO iv=1,n(ip)-1
          iedge=iedge+1
          IF(iedge.GT.nedges) GOTO 200
          ipoint=ipoint+1
          IF(ipoint.GT.npoints) GOTO 300
          grid%X(1,ipoint)=x(il+iv)
          grid%X(2,ipoint)=y(il+iv)
          grid%POINTS(1,iedge)=ipoint
          grid%POINTS(2,iedge)=ipoint+1
          grid%EDGEINDEX(1)%INDEXES(0,iedge)=iedge
          grid%EDGEINDEX(1)%INDEXES(1,iedge)=ip
         END DO  
         grid%POINTS(2,iedge)=ipoint0
        END IF
       END IF !IF(N(IP).GT.1)
       il=il+n(ip)
      END DO     
      grid%CELLS=0

C     UNITS
      grid%UNIT2SI=1e-3_gridman_dp
      grid%UNITS='MILLIMETER'
C     DEFAULT DESCRIPTION
      grid%DESCRIPTION='Converted from DG template'
      grid%EDGEINDEX(1)%DESCRIPTION='Index of polygon'
      grid%EDGEINDEX(1)%COLUMNS(1)='DG_IPOL'

C     CHECK
      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_TEMPLATE2GRID: ", 
     w                        "the resulting grid is incorrect" 
        GOTO 100
       END IF
      END IF  

      IF(gridman_dbg) 
     w   WRITE(gridman_unit,*) "GRIDMAN_TEMPLATE2GRID finished"

      RETURN
 100  WRITE(gridman_unit,*) "GRIDMAN_TEMPLATE2GRID terminated"
      CALL gridman_grid_deallocate(grid,ierr0)
      RETURN
 200  WRITE(gridman_unit,*) "ERROR in GRIDMAN_TEMPLATE2GRID: ", 
     w                      "internal error"  
      WRITE(gridman_unit,*)
     w     " Edge index is out of range, IEDGE, NEDGES ",iedge,nedges
      ierr=400
      CALL gridman_grid_deallocate(grid,ierr0)
      RETURN
 300  WRITE(gridman_unit,*) "ERROR in GRIDMAN_TEMPLATE2GRID: ", 
     w                      "internal error"  
      WRITE(gridman_unit,*)
     w     " Point index is out of range, IPOINT, NPOINTS ",
     w                                    ipoint, npoints
      ierr=400
      CALL gridman_grid_deallocate(grid,ierr0)
      RETURN

      END SUBROUTINE gridman_template2grid

C**********************************************************************************************
C> Convert GRIDMAN_GRID object into arrays which can be stored in simple template format
C> (DG template)
C**********************************************************************************************
C>
C> WARNING: arrays N, X, Y must NOT be allcoated - 
C> they are allocated inside the subroutine. 
C> Arrays which already exist will be re-allocated
C>
C> Each edge is stored as a separate polygon. This subroutine 
C> does not try to build chains of points
C
      SUBROUTINE gridman_grid2template(GRID,M,N,L,X,Y,IERR)
      USE gridman,ONLY:gridman_grid,gridman_sp,gridman_dp,
     u                 gridman_dbg,gridman_unit,gridman_check
      USE gridman_lib,ONLY:gridman_grid2d_check
      IMPLICIT NONE
C> Grid object to be converted
      TYPE(gridman_grid) :: GRID     
C> Number of polygons 
      INTEGER(GRIDMAN_SP),INTENT(OUT) :: M
C> Number of vertices in each polygon, N(M) 
      INTEGER(GRIDMAN_SP),ALLOCATABLE :: N(:)
C> Total number of points, L=SUM(N)
      INTEGER(GRIDMAN_SP),INTENT(OUT) :: L
C> X coordinates of the vertices, X(L) 
      REAL(GRIDMAN_DP),ALLOCATABLE :: X(:)
C> Y coordinates of the vertices, Y(L) 
      REAL(GRIDMAN_DP),ALLOCATABLE :: Y(:)
C> Error code
      INTEGER,INTENT(OUT) :: IERR

      INTEGER(GRIDMAN_SP) :: IEDGE,IP1,IP2,IL
      INTEGER :: RES,ST
      REAL(GRIDMAN_DP) :: FC

      IF(gridman_dbg) 
     w   WRITE(gridman_unit,*) "Starting GRIDMAN_GRID2TEMPLATE"

      ierr=0
      m=0
      l=0

C     CHECK INPUT
      IF(grid%TYPE.NE.2.OR.grid%PDIM.NE.2.OR.grid%EDIM.NE.2) THEN
       ierr=100
       WRITE(gridman_unit,*) "ERROR in GRIDMAN_GRID2TEMPLATE: ", 
     w                       "grid is not of 2D type" 
       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_GRID2TEMPLATE: ", 
     w                        "the input grid is incorrect" 
        RETURN
       END IF
      END IF  

C     ALLOCATE ARRAYS
      m=grid%NEDGES
      l=2*grid%NEDGES
      IF(ALLOCATED(n)) DEALLOCATE(n,stat=st)
      IF(ALLOCATED(x)) DEALLOCATE(x,stat=st)
      IF(ALLOCATED(y)) DEALLOCATE(y,stat=st)
      ALLOCATE(n(m),x(l),y(l),stat=st)
      IF(st.NE.0) THEN
       ierr=200
       WRITE(gridman_unit,*) "ERROR in GRIDMAN_GRID2TEMPLATE: ", 
     w                       "cannon perform allocation M, L ",m,l 
       m=0
       l=0 
       RETURN      
      END IF 
      
C     ASSIGN POINTS
      fc=grid%UNIT2SI*1e3_gridman_dp !TO TRANSFER INTO MM
      il=0
      DO iedge=1,m 
       n(iedge)=2      
       ip1=grid%POINTS(1,iedge)
       ip2=grid%POINTS(2,iedge)
       il=il+1
       x(il)=grid%X(1,ip1)*fc
       y(il)=grid%X(2,ip1)*fc
       il=il+1
       x(il)=grid%X(1,ip2)*fc
       y(il)=grid%X(2,ip2)*fc
      END DO

      IF(gridman_dbg) 
     w   WRITE(gridman_unit,*) "GRIDMAN_GRID2TEMPLATE finished"

      END SUBROUTINE gridman_grid2template

C**********************************************************************************************
C> Read from a simple template file (DG template) into GRIDMAN_GRID object
C>
C> WARNING: the GRID object  will be overwritten if elready exists
C>
C> This subroutine is a shell for GRIDMAN_TEMPLATE_READ_ARRAY
C>
C> GRIDMAN_TEMPLATE_READ_ARRAY AND GRIDMAN_TEMPLATE_READ_GRID
C> are combined in the interface GRIDMAN_TEMPLATE_READ
C>
C**********************************************************************************************
      SUBROUTINE gridman_template_read_grid(GRID,FNAME,IERR)
       USE gridman,ONLY:gridman_grid,gridman_sp,gridman_dp,
     u                  gridman_dbg,gridman_unit
       USE gridman_lib,ONLY:gridman_template_read_array,
     u                      gridman_template2grid
       IMPLICIT NONE
C> Resulting grid object 
       TYPE(gridman_grid) :: GRID
C> String containing the name of input file 
       CHARACTER(*),INTENT(IN) :: FNAME  

C> Error code
       INTEGER,INTENT(OUT) :: IERR

       INTEGER(GRIDMAN_SP) :: M,L
       INTEGER :: ST 
       INTEGER(GRIDMAN_SP),ALLOCATABLE :: N(:)
       REAL(GRIDMAN_DP),ALLOCATABLE :: X(:),Y(:)

       IF(gridman_dbg) 
     w   WRITE(gridman_unit,*) "Starting GRIDMAN_TEMPLATE_READ_GRID"

       ierr=0

       CALL gridman_template_read_array(fname,m,n,l,x,y,ierr) 
       IF(ierr.NE.0) GOTO 100

       CALL gridman_template2grid(grid,m,n,l,x,y,ierr)
       IF(ierr.NE.0) GOTO 200

       DEALLOCATE(n,stat=st)
       DEALLOCATE(x,stat=st)
       DEALLOCATE(y,stat=st)

       IF(gridman_dbg) 
     w   WRITE(gridman_unit,*) "GRIDMAN_TEMPLATE_READ_GRID finished"        
       RETURN

 200   IF(ALLOCATED(n)) DEALLOCATE(n,stat=st)
       IF(ALLOCATED(x)) DEALLOCATE(x,stat=st)
       IF(ALLOCATED(y)) DEALLOCATE(y,stat=st)
 100   WRITE(gridman_unit,*) "GRIDMAN_TEMPLATE_READ_GRID terminated"

      END SUBROUTINE gridman_template_read_grid

C**********************************************************************************************
C> Write GRIDMAN_GRID object in a simple template format (DG template)
C>
C> This subroutine is a shell for GRIDMAN_TEMPLATE_WRITE_ARRAY,
C> it uses GRIDMAN_GRID2TEMPLATE to convert grid into template
C>
C> GRIDMAN_TEMPLATE_WRITE_ARRAY and GRIDMAN_TEMPLATE_WRITE_GRID
C> are combined in the interface GRIDMAN_TEMPLATE_WRITE
C**********************************************************************************************
      SUBROUTINE gridman_template_write_grid(GRID,FNAME,IERR)
       USE gridman,ONLY:gridman_grid,gridman_sp,gridman_dp,
     u                  gridman_dbg,gridman_unit
       USE gridman_lib,ONLY:gridman_grid2template,
     u                      gridman_template_write_array
       IMPLICIT NONE
C> Grid object to be written
       TYPE(gridman_grid) :: GRID
C> String containing the name of input file 
       CHARACTER(*),INTENT(IN) :: FNAME  
C> Error code
       INTEGER,INTENT(OUT) :: IERR

       INTEGER(GRIDMAN_SP) :: M,L
       INTEGER :: ST
       INTEGER(GRIDMAN_SP),ALLOCATABLE :: N(:)
       REAL(GRIDMAN_DP),ALLOCATABLE :: X(:),Y(:)

       IF(gridman_dbg) 
     w   WRITE(gridman_unit,*) "Starting GRIDMAN_TEMPLATE_WRITE_GRID" 
    
       ierr=0

       CALL gridman_grid2template(grid,m,n,l,x,y,ierr)
       IF(ierr.NE.0) GOTO 100

       CALL gridman_template_write_array(fname,m,n,l,x,y,ierr) 
       IF(ierr.NE.0) GOTO 200

       DEALLOCATE(n,stat=st)
       DEALLOCATE(x,stat=st)
       DEALLOCATE(y,stat=st)

       IF(gridman_dbg) 
     w   WRITE(gridman_unit,*) "GRIDMAN_TEMPLATE_WRITE_GRID finished"  
       RETURN

 200   IF(ALLOCATED(n)) DEALLOCATE(n,stat=st)
       IF(ALLOCATED(x)) DEALLOCATE(x,stat=st)
       IF(ALLOCATED(y)) DEALLOCATE(y,stat=st)
 100   WRITE(gridman_unit,*) "GRIDMAN_TEMPLATE_WRITE_GRID terminated"

      END SUBROUTINE gridman_template_write_grid