utils_vtk.F

Go to the documentation of this file.

!                 ************
                  MODULE utils_vtk
!                 ************
!
!***********************************************************************
! HERMES  V7P0                                                 2015
!***********************************************************************
!
!brief    INTERFACES OF VTK PUBLIC SUBROUTINES
!
!history YOANN AUDOUIN
!+       13/11/2016
!+       V7PR2
!+       First version
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE declarations_special
      USE hash_table
      USE lib_vtk_io
!
      IMPLICIT NONE
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!     NEW OBJECT RELATED TO MESHES IN A VTK FILE
      TYPE vtk_info
        INTEGER :: file_id
        CHARACTER(LEN=MAXLENFILE) :: file_name
        CHARACTER(LEN=80) :: title
        INTEGER, ALLOCATABLE :: cell_type(:)
        INTEGER, ALLOCATABLE :: connect(:)
        INTEGER :: nelem, npoin
        INTEGER :: ivar, nvar
        INTEGER :: mesh_dim
        DOUBLE PRECISION, ALLOCATABLE :: x(:)
        DOUBLE PRECISION, ALLOCATABLE :: y(:)
        DOUBLE PRECISION, ALLOCATABLE :: z(:)
      END TYPE vtk_info
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      INTEGER :: nfiles = 0
!
!     HASH TABLE FOR FILES
      INTEGER :: hash(max_file) = 0
!
      TYPE(vtk_info) :: vtk_obj_tab(max_file)
!
!-----------------------------------------------------------------------
!
      CONTAINS
!
!-----------------------------------------------------------------------
!
!
      ! CONVERT AN INT INTO A STRING WITH ZEROS
      SUBROUTINE number2string(RES,VAL)
!
!brief    Fill a 00000000 string with val
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| RES         |-->| OUTPUT STRING
!| VAL         |<--| INTEGER TO CONVERT
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      USE declarations_special
!
      IMPLICIT NONE
!
      INTEGER, INTENT(IN) :: VAL
!
      CHARACTER(LEN=8), INTENT(OUT) :: RES
!
      IF(val.LT.10) THEN
        WRITE(res,'(A7,I1)') '0000000',val
      ELSE IF(val.LT.100) THEN
        WRITE(res,'(A6,I2)') '000000',val
      ELSE IF(val.LT.1000) THEN
        WRITE(res,'(A5,I3)') '00000',val
      ELSE IF(val.LT.10000) THEN
        WRITE(res,'(A4,I4)') '0000',val
      ELSE IF(val.LT.100000) THEN
        WRITE(res,'(A3,I5)') '000',val
      ELSE IF(val.LT.1000000) THEN
        WRITE(res,'(A2,I6)') '00',val
      ELSE IF(val.LT.10000000) THEN
        WRITE(res,'(A1,I7)') '0',val
      ELSE
        WRITE(res,'(I8)') val
      ENDIF
      END SUBROUTINE number2string
!
!-----------------------------------------------------------------------
!
      SUBROUTINE open_mesh_vtk(FILE_NAME,FILE_ID,OPEN_MODE,IERR)
!
!BRIEF    OPENS A MESH FILE
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| FILE_NAME      |<--| NAME OF THE VTK FILE
!| FILE_ID        |<--| VTK FILE DESCRIPTOR
!| OPEN_MODE      |<--| OPENING MODE (READ/READWRITE/WRITE)
!| IERR           |-->| ERROR TAG
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      IMPLICIT NONE
!
      CHARACTER(LEN=*), INTENT(IN)  :: FILE_NAME
      INTEGER,          INTENT(OUT)  :: FILE_ID
      CHARACTER(LEN=9), INTENT(IN)  :: OPEN_MODE
      INTEGER,          INTENT(OUT) :: IERR
!
#if !defined (DISABLE_VTK)
      INTEGER :: VTK_ID
!
!-----------------------------------------------------------------------
!
      SELECT CASE(open_mode)
        CASE('READ     ')
          ierr = hermes_write_only_err
        CASE('READWRITE')
        CASE('WRITE    ')
        CASE DEFAULT
!         OPEN MODE IS INVALID
          ierr = hermes_invalid_open_mode_err
      END SELECT
!
      CALL get_free_id(file_id)
!     CREATE A NEW OBJECT FOR THE NEW VTK FILE
      CALL add_obj(hash,file_id,vtk_id,ierr)
      IF(ierr.NE.0) THEN
        error_message = 'ERROR IN '//
     &       trim(file_name)//': '//
     &       'OPEN_MESH_MED:ADD_OBJ'
        RETURN
      ENDIF
!
      vtk_obj_tab(vtk_id)%FILE_NAME = file_name
!     Opening bogus file to secure the id
      OPEN(unit=file_id,file='bogus.txt',
     &     action='WRITE',iostat=ierr)
      IF(ierr.NE.0) THEN
        error_message = 'ERROR IN '//
     &       trim(vtk_obj_tab(vtk_id)%FILE_NAME)//': '//
     &       'OPEN_MESH_VTK:BOGUS'
        RETURN
      ENDIF
!
#else
!
!     VTK LIBRARY NOT LOADED
      ierr = hermes_vtk_not_loaded_err
!
#endif
!
      RETURN
      END SUBROUTINE
!
!-----------------------------------------------------------------------
!
      SUBROUTINE close_mesh_vtk (FILE_ID, IERR)
!
!BRIEF    CLOSES A MESH FILE
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| FILE_ID        |<--| VTK FILE DESCRIPTOR
!| IERR           |-->| ERROR TAG
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      IMPLICIT NONE
!
      INTEGER, INTENT(IN)  :: FILE_ID
      INTEGER, INTENT(OUT) :: IERR
!
      INTEGER VTK_ID
!
!-----------------------------------------------------------------------
!
#if !defined (DISABLE_VTK)
!
      CALL get_obj(hash,file_id,vtk_id,ierr)
      IF(ierr.NE.0) THEN
        error_message = 'ERROR WITH ID '//i2char(file_id)//': '//
     &       'CLOSE_MESH_VTK:GET_OBJ'
        RETURN
      ENDIF
!
      ! Clearing id in the hash table
      hash(vtk_id) = 0
      nfiles = nfiles - 1
!
      ! CLEAN UP OF STRUCTURE
      DEALLOCATE(vtk_obj_tab(vtk_id)%CONNECT)
      DEALLOCATE(vtk_obj_tab(vtk_id)%CELL_TYPE)
      DEALLOCATE(vtk_obj_tab(vtk_id)%X)
      DEALLOCATE(vtk_obj_tab(vtk_id)%Y)
      CLOSE(file_id)
!
#else
!
!     VTK LIBRARY NOT LOADED
      ierr = hermes_vtk_not_loaded_err
!
#endif
!
      RETURN
      END SUBROUTINE
!
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
      SUBROUTINE set_header_vtk(FILE_ID,TITLE,NVAR,IERR)
!
!BRIEF    WRITES THE TITLE OF THE MESH
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| FILE_ID        |<--| VTK FILE DESCRIPTOR
!| TITLE          |<--| TITLE OF THE MESH
!| NVAR           |<--| Number of variable
!| IERR           |-->| ERROR TAG
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      IMPLICIT NONE
!
      INTEGER, INTENT(IN)  :: FILE_ID
      INTEGER, INTENT(IN)  :: NVAR
      CHARACTER(LEN=80), INTENT(IN)  :: TITLE
      INTEGER, INTENT(OUT) :: IERR
!
#if !defined (DISABLE_VTK)
      INTEGER :: VTK_ID
!
!-----------------------------------------------------------------------
!
      CALL get_obj(hash,file_id,vtk_id,ierr)
      IF(ierr.NE.0) THEN
        error_message = 'ERROR WITH ID '//i2char(file_id)//': '//
     &       'SET_HEADER_VTK:GET_OBJ'
        RETURN
      ENDIF
!
      vtk_obj_tab(vtk_id)%TITLE = title
      vtk_obj_tab(vtk_id)%NVAR = nvar
!
#else
!
!     VTK LIBRARY NOT LOADED
      ierr = hermes_vtk_not_loaded_err
!
#endif
      RETURN
      END SUBROUTINE
!
!-----------------------------------------------------------------------
!
      SUBROUTINE set_mesh_vtk(FILE_ID,NB_DIM_MESH,TYPE_ELEM,
     &  NDP,NELEM,NPOIN,IKLE,X,Y,IERR)
!
!BRIEF    WRITES A MESH IN A VTK FILE
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| FILE_ID        |<--| VTK FILE DESCRIPTOR
!| NB_DIM_PB      |<--| NUMBER OF DIMENSION OF THE PHYSICAL DOMAIN
!| NB_DIM_MESH    |<--| NUMBER OF DIMENSION OF THE MESH
!| TYPE_ELEM      |<--| TYPE OF THE ELEMENT IN SLF FORMAT
!| NDP            |<--| NUMBER OF POINTS PER ELEMENT
!| NPTFR          |<--| NUMBER OF BOUNDARY NODES
!| NPTIR          |<--| NUMBER OF INTERFACE NODES FOR THE SUB-DOMAIN
!| NELEM          |<--| TOTAL NUMBER OF ELEMENTS
!| NPOIN          |<--| TOTAL NUMBER OF NODES
!| IKLE           |<--| CONNECTIVITY TABLE
!| IPOBO          |<--| LOGICAL TABLE FOR BOUNDARY NODES
!| KNOLG          |<--| NODES INDEX TABLE FROM LOCAL TO GLOBAL
!| X              |<--| X COORDINATES OF THE NODES
!| Y              |<--| Y COORDINATES OF THE NODES
!| NPLAN          |<--| NUMBER OF PLANES
!| IERR           |-->| ERROR TAG
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      IMPLICIT NONE
!
      INTEGER, INTENT(IN) :: FILE_ID, NB_DIM_MESH
      INTEGER, INTENT(IN) :: NELEM, NPOIN, TYPE_ELEM, NDP
      INTEGER, INTENT(IN) :: IKLE(ndp*nelem)
      DOUBLE PRECISION, INTENT(IN) :: X(npoin)
      DOUBLE PRECISION, INTENT(IN) :: Y(npoin)
      INTEGER, INTENT(OUT) :: IERR
!
#if !defined (DISABLE_VTK)
!
      INTEGER :: VTK_ID
      INTEGER :: I,J
!
!-----------------------------------------------------------------------
!
      CALL get_obj(hash,file_id,vtk_id,ierr)
      IF(ierr.NE.0) THEN
        error_message = 'ERROR WITH ID '//i2char(file_id)//': '//
     &       'SET_MESH_VTK:GET_OBJ'
        RETURN
      ENDIF
!
      ! DEFINING THE CELLS TYPE
      ALLOCATE(vtk_obj_tab(vtk_id)%CELL_TYPE(nelem),stat=ierr)
      IF(ierr.NE.0) THEN
        error_message = 'ERROR IN '//
     &       trim(vtk_obj_tab(vtk_id)%FILE_NAME)//': '//
     &       'ALLOCATING CELL_TYPE'
        RETURN
      ENDIF
      ! VALUE CHOSEN IN PDF VTK FILE FORMAT (GOOGLE)
      SELECT CASE(type_elem)
      CASE(triangle_elt_type)
        vtk_obj_tab(vtk_id)%CELL_TYPE(:) = 5
      CASE(quadrangle_elt_type)
        vtk_obj_tab(vtk_id)%CELL_TYPE(:) = 9
      CASE(tetrahedron_elt_type)
        vtk_obj_tab(vtk_id)%CELL_TYPE(:) = 10
      CASE(prism_elt_type)
        vtk_obj_tab(vtk_id)%CELL_TYPE(:) = 13
      END SELECT
      ! BUILDING THE CONNECTION TABLE FOLLOWING THE VTK STANDARD
      ! FOR EACH CELL :
      !   FIRST THE NUMBER OF NODES COMPOSING THE CELL
      !   THEN THE NUMBER OF EACH NODE
      ! WE NEED TO RENUMBER TO NODES 0 TO NELEM-1
      ALLOCATE(vtk_obj_tab(vtk_id)%CONNECT(nelem*(ndp+1)),stat=ierr)
      IF(ierr.NE.0) THEN
        error_message = 'ERROR IN '//
     &       trim(vtk_obj_tab(vtk_id)%FILE_NAME)//': '//
     &       'ALLOCATING CONNECT'
        RETURN
      ENDIF
      DO i=1,nelem
        vtk_obj_tab(vtk_id)%CONNECT((i-1)*(ndp+1)+1) = ndp
        DO j=1,ndp
          vtk_obj_tab(vtk_id)%CONNECT((i-1)*(ndp+1)+j+1) =
     &            ikle((j-1)*nelem+i)-1
        ENDDO
      ENDDO
      vtk_obj_tab(vtk_id)%NELEM = nelem
      vtk_obj_tab(vtk_id)%NPOIN = npoin

      ALLOCATE(vtk_obj_tab(vtk_id)%X(npoin),stat=ierr)
      IF(ierr.NE.0) THEN
        error_message = 'ERROR IN '//
     &       trim(vtk_obj_tab(vtk_id)%FILE_NAME)//': '//
     &       'ALLOCATING X'
        RETURN
      ENDIF
      vtk_obj_tab(vtk_id)%X = x
      ALLOCATE(vtk_obj_tab(vtk_id)%Y(npoin),stat=ierr)
      IF(ierr.NE.0) THEN
        error_message = 'ERROR IN '//
     &       trim(vtk_obj_tab(vtk_id)%FILE_NAME)//': '//
     &       'ALLOCATING Y'
        RETURN
      ENDIF
      vtk_obj_tab(vtk_id)%Y = y

      vtk_obj_tab(vtk_id)%MESH_DIM = nb_dim_mesh
!
#endif
!
      RETURN
      END SUBROUTINE
!
!-----------------------------------------------------------------------
!
      SUBROUTINE add_data_vtk (FILE_ID, VAR_NAME, RECORD,
     &  FIRST_VAR, VAR_VALUE, N, IERR)
!
!BRIEF    WRITE DATA VALUES FOR A GIVEN VARIABLE
!         ONLY FOR VALUES WITH DOUBLE PRECISION, ONE COMPONENT AND VALUE
!         ON THE NODES
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| FILE_ID        |<--| VTK FILE DESCRIPTOR
!| VAR_NAME       |<--| NAME OF THE DATA
!| RECORD         |<--| ITERATION NUMBER
!| FIRST_VAR      |<--| TRUE IF FIRST VARIABLE
!| VAR_VALUE      |<--| TABLE OF VARIABLE VALUES
!| N              |<--| NUMBER OF ELEMENTS
!| IERR           |-->| ERROR TAG
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      IMPLICIT NONE
!
      INTEGER, INTENT(IN) :: FILE_ID, N
      INTEGER, INTENT(IN) :: RECORD
      CHARACTER(LEN=32), INTENT(IN) :: VAR_NAME
      DOUBLE PRECISION, INTENT(IN) :: VAR_VALUE(n)
      LOGICAL, INTENT(IN) :: FIRST_VAR
      INTEGER, INTENT(OUT) :: IERR
!
#if !defined (DISABLE_VTK)
      INTEGER :: VTK_ID
      DOUBLE PRECISION, ALLOCATABLE :: Z(:)
      CHARACTER(LEN=MAXLENFILE) :: FILENAME
      CHARACTER(LEN=8) :: TIMESTR
      CHARACTER(LEN=32) :: VAR_NAME_VTK
!
!-----------------------------------------------------------------------
!
      CALL get_obj(hash,file_id,vtk_id,ierr)
      IF(ierr.NE.0) THEN
        error_message = 'ERROR WITH ID '//i2char(file_id)//': '//
     &       'ADD_DATA:GET_OBJ'
        RETURN
      ENDIF
!
      IF(first_var) THEN
        filename = repeat(' ',maxlenfile)
        CALL number2string(timestr,record)
        filename = trim(vtk_obj_tab(vtk_id)%FILE_NAME)//timestr
        ! INITIALISING THE VTK FILE
        ierr=vtk_ini('ASCII',filename,vtk_obj_tab(vtk_id)%TITLE,
     &               'UNSTRUCTURED_GRID')
        IF(ierr.NE.0) THEN
          error_message = 'ERROR IN '//
     &         trim(vtk_obj_tab(vtk_id)%FILE_NAME)//': '//
     &         'ADD_DATA:VTK_INI'
          RETURN
        ENDIF
        ! WRITTING THE NODES COORDINATES
        ALLOCATE(z(vtk_obj_tab(vtk_id)%NPOIN),stat=ierr)
        IF(ierr.NE.0) THEN
          error_message = 'ERROR IN '//
     &         trim(vtk_obj_tab(vtk_id)%FILE_NAME)//': '//
     &         'ALLOCATING Z'
          RETURN
        ENDIF
        IF(vtk_obj_tab(vtk_id)%MESH_DIM.EQ.3) THEN
          z = var_value
        ELSE
          z = 0.d0
        ENDIF
        ierr=vtk_geo(vtk_obj_tab(vtk_id)%NPOIN,
     &               vtk_obj_tab(vtk_id)%X,
     &               vtk_obj_tab(vtk_id)%Y,
     &               z)
        DEALLOCATE(z)
        IF(ierr.NE.0) THEN
          error_message = 'ERROR IN '//
     &         trim(vtk_obj_tab(vtk_id)%FILE_NAME)//': '//
     &         'ADD_DATA:VTK_GEO'
          RETURN
        ENDIF
        ierr=vtk_con(vtk_obj_tab(vtk_id)%NELEM,
     &               vtk_obj_tab(vtk_id)%CONNECT,
     &               vtk_obj_tab(vtk_id)%CELL_TYPE)
        IF(ierr.NE.0) THEN
          error_message = 'ERROR IN '//
     &         trim(vtk_obj_tab(vtk_id)%FILE_NAME)//': '//
     &         'ADD_DATA:VTK_CON'
          RETURN
        ENDIF
        ! ADDING THE NUMBER OF NODES ON WHICH THERE IS GONNA BE A
        ! VARAIBLE HERE EVERY ONE
        ierr=vtk_dat(vtk_obj_tab(vtk_id)%NPOIN,'node')
        IF(ierr.NE.0) THEN
          error_message = 'ERROR IN '//
     &         trim(vtk_obj_tab(vtk_id)%FILE_NAME)//': '//
     &         'ADD_DATA:VTK_DAT'
          RETURN
        ENDIF
        ! Initializing variable counter
        vtk_obj_tab(vtk_id)%IVAR = 0
      ENDIF

      ! REMOVE WHITE SPACES IN VAR NAMES
      ! TO BE READABLE IN PARAVIEW 5.0
      vtk_obj_tab(vtk_id)%IVAR = vtk_obj_tab(vtk_id)%IVAR + 1
      var_name_vtk = var_name
      CALL blanc2uscore(var_name_vtk, len(var_name_vtk))
      ierr=vtk_var(vtk_obj_tab(vtk_id)%NPOIN, var_name_vtk, var_value)
      IF(ierr.NE.0) THEN
        error_message = 'ERROR IN '//
     &       trim(vtk_obj_tab(vtk_id)%FILE_NAME)//': '//
     &       'ADD_DATA:VTK_VAR '
        RETURN
      ENDIF

      ! Closing the file after the last varialbe was written
      IF(vtk_obj_tab(vtk_id)%IVAR.EQ.vtk_obj_tab(vtk_id)%NVAR) THEN
!       CLOSE VTK FILE
        ! CLOSING THE FILE
        ierr=vtk_end()
        IF(ierr.NE.0) THEN
          error_message = 'ERROR IN '//
     &         trim(vtk_obj_tab(vtk_id)%FILE_NAME)//': '//
     &         'ADD_DATA:VTK_END'
          RETURN
        ENDIF
      ENDIF
!
#else
!
!     VTK LIBRARY NOT LOADED
      ierr = hermes_vtk_not_loaded_err
!
#endif
      RETURN
      END SUBROUTINE
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
      END MODULE utils_vtk