triang.f

Go to the documentation of this file.

C> @file apps/triang.f
C> Combine plasma grid (fort.30) with triangular grid in the "void" region.
C> Replacement for TRIAGEOM. 
C>
C> Based on GRIDMAN library
C 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/>.

      PROGRAM triang
       USE gridman,ONLY:gridman_grid,gridman_dp,
     u                  gridman_unit,gridman_dbg,gridman_check
       USE gridman_lib,ONLY:gridman_carre_read30,gridman_grid2d_triang,
     u                      gridman_tria_read,
     u                      gridman_grid2d_merge,gridman_tria_write,
     u                      gridman_vtk_grid2d_write,
     u                      gridman_grid_write,gridman_template_write,
     u                      gridman_index_deallocate
       IMPLICIT NONE
       INTRINSIC len_trim,get_command_argument,adjustl,trim    

       INTEGER,PARAMETER :: file_length=1024
C      INPUT PARAMETERS
       CHARACTER(LEN=FILE_LENGTH) :: 
     c           fort30_in,fort33_in,fort34_in,fort35_in  
       REAL(GRIDMAN_DP) :: tol
       CHARACTER(LEN=FILE_LENGTH) :: 
     c           fort33_out,fort34_out,fort35_out, 
     c           vtk_out,grid_out,template_out 
       LOGICAL :: dbgmod,lcheck 

C      LOCAL VARIABLES
       TYPE(gridman_grid) :: plasma,combined,tria,tria_out
       INTEGER :: ierr,res,lt
       CHARACTER(LEN=8) :: ftmp

       CALL get_command_argument(1,ftmp)
       IF(len_trim(ftmp).GT.0) THEN
        IF(adjustl(trim(ftmp)).EQ.'help') THEN
         CALL print_description 
         stop
        END IF
       ELSE
C      0. READ INPUT FROM NAMELIST
        CALL read_input
       END IF

C      1. READ FORT.30
       WRITE(gridman_unit,*) "TRIANG: read ",trim(fort30_in)
       CALL gridman_carre_read30(plasma,fort30_in,ierr,.true.)
       IF(ierr.NE.0) CALL triang_error("failed to read fort.30")  

C      2. READ INPUT TRIANGULAR GRID
       WRITE(gridman_unit,*) "TRIANG: read "//trim(fort33_in)
       WRITE(gridman_unit,*)  "      "//trim(fort34_in)//" "
     w                               //trim(fort35_in)
       CALL gridman_tria_read(tria,
     c                        (/fort33_in,fort34_in,fort35_in/),ierr)
       IF(ierr.NE.0) CALL triang_error("failed to read triangular grid")  

C      3. PREPARE INDEXES: cell indexes will be taken from the plasma grid
       WRITE(gridman_unit,*) "TRIANG: preparing indexes" 
       CALL gridman_index_deallocate(tria%CELLINDEX(1),ierr)
       IF(ierr.NE.0) CALL triang_error("")  
       DEALLOCATE(tria%CELLINDEX)
       tria%NCELLINDEX=0

C      4. MERGE GRIDS
       WRITE(gridman_unit,*) "TRIANG: merging grids"
       CALL gridman_grid2d_merge(combined,tria,plasma,tol,ierr)
       IF(ierr.NE.0) 
     f    CALL triang_error("merging failed. Try to change TOL")  

C      5. TRIANGULATE
       WRITE(gridman_unit,*) "TRIANG: triangulate plasma grid"
       CALL gridman_grid2d_triang(tria_out,combined,ierr)
       IF(ierr.NE.0) CALL triang_error("triangulation failed")  
       tria_out%DESCRIPTION=trim(tria_out%DESCRIPTION)//
     /                      ' - Created by TRIANG '

C      6. WRITE DOWN COMBINED GRID
       WRITE(gridman_unit,*) "TRIANG: writing combined grid into"
       WRITE(gridman_unit,*) "  ",trim(fort33_out)
       WRITE(gridman_unit,*) "  ",trim(fort34_out)
       WRITE(gridman_unit,*) "  ",trim(fort35_out)
       
       CALL gridman_tria_write(tria_out,
     c                (/fort33_out,fort34_out,fort35_out/),ierr)   
       IF(ierr.NE.0) CALL triang_error("writing combined grid failed")  

C      7. CHECK THE CREATED GRID
       WRITE(gridman_unit,*) "TRIANG: checking the created grid"
       CALL gridman_tria_read(tria,
     c            (/fort33_out,fort34_out,fort35_out/),ierr)
       IF(ierr.NE.0) 
     f    CALL triang_error("the created grid cannot be read")  
       CALL gridman_grid2d_check(tria,res,ierr)
       IF(res.NE.0.OR.ierr.NE.0) 
     c    CALL triang_error("the created grid does not pass checks")  

C      8. WRITE IN OTHER FORMATS
       lt=len_trim(vtk_out)
       IF(lt.GT.0) THEN
        WRITE(gridman_unit,*) "TRIANG: writing combined grid"
        WRITE(gridman_unit,*) " into VTK file ",vtk_out(1:lt)
        CALL gridman_vtk_grid2d_write(
     c       vtk_out,tria_out%DESCRIPTION,tria_out,
     c       tria_out%NCELLS,0,0,0,0,ierr)
        IF(ierr.NE.0) CALL triang_error("could not write VTK file")  
       END IF
       lt=len_trim(template_out)
       IF(lt.GT.0) THEN
        WRITE(gridman_unit,*) "TRIANG: writing combined grid"
        WRITE(gridman_unit,*) 
     w      " into DG template file ",template_out(1:lt)
        CALL gridman_template_write(tria_out,template_out,ierr)
        IF(ierr.NE.0) CALL triang_error("could not write template file")  
       END IF
       lt=len_trim(grid_out)
       IF(lt.GT.0) THEN
        WRITE(gridman_unit,*) "TRIANG: writing combined grid"
        WRITE(gridman_unit,*) " into GRIDMAN-grid file ",grid_out(1:lt)
        CALL gridman_grid_write(tria_out,grid_out,ierr)
        IF(ierr.NE.0) CALL triang_error("could not write grid file")  
       END IF
 
C      DEALLOCATE MEMORY
       CALL clean_all 
       WRITE(gridman_unit,*) "TRIANG COMPLETED"

      CONTAINS

C
      SUBROUTINE read_input
       INTEGER :: io
       namelist /triang/ fort30_in,fort33_in,fort34_in,fort35_in,
     n                   tol,fort33_out,fort34_out,fort35_out, 
     n                   vtk_out,grid_out,template_out,dbgmod,lcheck 

C      DEFAULTS
       fort30_in=''
       fort33_in=''
       fort34_in=''
       fort35_in=''
       fort33_out=''
       fort34_out=''
       fort35_out=''
       vtk_out=''
       grid_out=''   
       template_out=''
       tol=1e-5
       dbgmod=.false.
       lcheck =.false.
   
       rewind(5)
       READ(5,nml=triang,iostat=io)
       IF(io.NE.0) CALL triang_error("can not read namelist TRIANG")

       gridman_dbg=dbgmod
       gridman_check=lcheck       

      END SUBROUTINE read_input
C 
      SUBROUTINE clean_all
       CALL gridman_grid_deallocate(plasma,ierr)
       CALL gridman_grid_deallocate(combined,ierr)
       CALL gridman_grid_deallocate(tria,ierr)
       CALL gridman_grid_deallocate(tria_out,ierr)
      END SUBROUTINE clean_all

C     Print error message and exit
      SUBROUTINE triang_error(STR)           
        USE gridman
        IMPLICIT NONE
        CHARACTER(*),INTENT(IN) :: str
        INTEGER :: lt
        INTRINSIC len_trim
        lt=len_trim(str)
        IF(lt.GT.0) 
     w     WRITE(gridman_unit,*) "ERROR in TRIANG: "//str(1:lt)
        stop "ERROR in TRIANG - see log output. "
       END SUBROUTINE triang_error

C
      SUBROUTINE print_description
       INTRINSIC index,len,trim
       CHARACTER(LEN=FILE_LENGTH) :: path
       CHARACTER(LEN=128) :: str
       INTEGER :: i,io
       CALL get_command_argument(0,path)
       i=index(path,'/',.true.) 
       IF(i.GT.len(path)) THEN 
        path=''
       ELSE 
        path=path(1:i)
       END IF

       IF(gridman_dbg)  WRITE(gridman_unit,*) " PATH:",trim(path)

       OPEN(unit=3,file=trim(path)//'triang.parameters.description',
     o      status='OLD',iostat=io) 
       IF(io.NE.0) GOTO 100

       DO
        READ(3,'(A)',iostat=io,end=200) str
        IF(io.NE.0) GOTO 200
        WRITE(*,*) trim(str)
       END DO
 200   RETURN

 100   WRITE(gridman_unit,*)
     w      "Could not find triang.parameters.description",
     w      trim(path)
       WRITE(gridman_unit,*)
     w  "Use 'which triang' to invoke via the full path"

       END SUBROUTINE print_description

      END PROGRAM triang