read_mesh_info.f

Go to the documentation of this file.

!                   *************************
                    SUBROUTINE read_mesh_info
!                   *************************
!
     &(fformat,nfic,title,nvar,npoin,typ_elem,nelem,nptfr,nptir,
     & ndp,nplan,x_orig,y_orig,typ_bnd_elem,nelebd)
!
!***********************************************************************
! HERMES   V6P3                                   21/08/2010
!***********************************************************************
!
!BRIEF    READS OR COMPUTES THE VALUES OF NPOIN, NELEM, NPTFR,
!+                MXPTVS, MXELVS IN THE GEOMETRY FILE (CHANNEL NGEO).
!
!
!history  Y AUDOUIN (LNHE)
!+        21/05/2015
!+        V7P0
!+   First version.
!
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| FFORMAT      |-->| FORMAT OF THE FILE
!| TITLE        |<--| TITLE OF THE MESH
!| NVAR         |<--| NUMBER OF VARIABLES
!| NPOIN        |<--| NUMBER OF MESH NODES
!| TYP_ELEM     |<--| TYPE OF ELEMENT
!| NELEM        |<--| NUMBER OF ELEMENTS
!| NDP          |<--| NUMBER OF ELEMENT FACES
!| NPLAN        |<--| NUMBER OF PLAN
!| NPTFR        |<--| NUMBER OF BOUNDARY NODES
!| NPTIR        |<--| NUMBER OF INTERFACES
!| NFIC         |<--| FILE TO READ
!| X_ORIG       |<--| Off set of the X coordinates
!| Y_ORIG       |<--| Off set of the Y coordinates
!| TYP_BND_ELEM |<--| TYPE OF BOUNDARY ELEMENTS
!| NELEBD       |<--| NUMBER OF BOUNDARY ELEMENTS
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE interface_hermes
      USE bief, ONLY: ncsize
!
      USE declarations_special
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      CHARACTER(LEN=8), INTENT(IN)   :: FFORMAT ! FILE FORMAT
      CHARACTER(LEN=80), INTENT(OUT) :: TITLE
      INTEGER, INTENT(OUT)           :: NVAR   ! NUMBER OF VARIABLES
      INTEGER, INTENT(OUT)           :: NPOIN  ! NUMBER OF MESH NODES
      INTEGER, INTENT(OUT)           :: TYP_ELEM  ! TYPE OF ELEMENT
      INTEGER, INTENT(OUT)           :: NELEM  ! NUMBER OF ELEMENTS
      INTEGER, INTENT(OUT)           :: NDP    ! NUMBER OF ELEMENT FACES
      INTEGER, INTENT(OUT)           :: NPLAN  ! NUMBER OF PLAN
      INTEGER, INTENT(OUT)           :: NPTFR  ! NUMBER OF BOUNDARY NODES
      INTEGER, INTENT(OUT)           :: NPTIR  ! NUMBER OF INTERFACES
      INTEGER, INTENT(OUT)           :: X_ORIG  ! X ORIGIN
      INTEGER, INTENT(OUT)           :: Y_ORIG  ! Y_ORIGIN
      INTEGER, INTENT(IN)            :: NFIC   ! FILE TO READ
      INTEGER,OPTIONAL,INTENT(OUT)   :: TYP_BND_ELEM ! TYPE OF BND ELEMENTS
      INTEGER,OPTIONAL,INTENT(OUT)   :: NELEBD ! NUMBER OF BOUNDARY ELEMENTS
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER IERR,NDIM
      INTEGER NNELEBD, TYP, TYP_BND
      INTEGER NNELEM
!
!-----------------------------------------------------------------------
!
      CALL get_mesh_title(fformat,nfic,title,ierr)
      CALL check_call(ierr,'READ_MESH_INFO:GET_MESH_TITLE')
!
      CALL get_data_nvar(fformat,nfic,nvar,ierr)
      CALL check_call(ierr,'READ_MESH_INFO:GET_DATA_NVAR')
!
!     IDENTIFY THE MESH ELEMENT TYPE AND READ THE NUMBER OF ELEMENT
!
      CALL get_mesh_dimension(fformat,nfic,ndim,ierr)
      CALL check_call(ierr,'READ_MESH_INFO:GET_MESH_DIMENSION')
      IF(ndim.EQ.3) THEN
        ! TEST FOR TETRAEDRON
        typ = tetrahedron_elt_type
        CALL get_mesh_nelem(fformat,nfic,typ,nnelem,ierr)
        CALL check_call(ierr,
     &                  'READ_MESH_INFO:GET_MESH_NELEM:TETRAHEDRON')
        IF(nnelem.EQ.0) THEN
          ! TEST FOR PRISM
          typ = prism_elt_type
          CALL get_mesh_nelem(fformat,nfic,typ,nnelem,ierr)
          CALL check_call(ierr,
     &                    'READ_MESH_INFO:GET_MESH_NELEM:PRISM')
          IF(nnelem.EQ.0) THEN
            WRITE(lu,*) 'NO 3D ELEMENTS IN A 3D MESH'
            CALL plante(1)
          ENDIF
        ENDIF
      ELSE
        ! TEST FOR TRIANGLE
        typ = triangle_elt_type
        CALL get_mesh_nelem(fformat,nfic,typ,nnelem,ierr)
        CALL check_call(ierr,'READ_MESH_INFO:GET_MESH_NELEM:TRIANGLE')
        IF(nnelem.EQ.0) THEN
          ! TEST FOR QUADRANGLE
          typ = quadrangle_elt_type
          CALL get_mesh_nelem(fformat,nfic,typ,nnelem,ierr)
          CALL check_call(ierr,
     &                    'READ_MESH_INFO:GET_MESH_NELEM:QUADRANGLE')
          IF(nnelem.EQ.0) THEN
            WRITE(lu,*) 'NO 2D ELEMENTS IN A 2D MESH'
            CALL plante(1)
          ENDIF
        ENDIF
      ENDIF
      typ_elem = typ
      nelem = nnelem
!
!     IDENTIFY THE BOUNDARY ELEMENT TYPE AND READ THE NUMBER OF ELEMENT
!
      IF(ndim.EQ.3) THEN
        ! TEST FOR TRIANGLE
        typ_bnd = triangle_bnd_elt_type
        CALL get_bnd_nelem(fformat,nfic,typ_bnd,nnelebd,ierr)
        CALL check_call(ierr,
     &                  'READ_MESH_INFO:GET_MESH_NELEM:TRIANGLE_BND')
        IF(nnelebd.EQ.0) THEN
          WRITE(lu,*) 'NO BND ELEMENTS IN A 3D MESH'
        ENDIF
      ELSE
        ! TEST FOR EDGE
        typ_bnd = edge_bnd_elt_type
        IF(ncsize.GT.1) typ_bnd = point_bnd_elt_type
        CALL get_bnd_nelem(fformat,nfic,typ_bnd,nnelebd,ierr)
        CALL check_call(ierr,
     &                  'READ_MESH_INFO:GET_MESH_NELEM:EDGE')
        IF(nnelebd.EQ.0) THEN
          ! TEST FOR POINTS
          typ_bnd = point_bnd_elt_type
          CALL get_bnd_nelem(fformat,nfic,typ_bnd,nnelebd,ierr)
          CALL check_call(ierr,
     &                    'READ_MESH_INFO:GET_MESH_NELEM:POINT')
          IF(nnelebd.EQ.0) THEN
            WRITE(lu,*) 'NO BND ELEMENTS IN A 2D MESH'
            ! Only stop if in serial cause the
            typ_bnd = type_null
          ENDIF
        ENDIF
      ENDIF
      IF(PRESENT(nelebd)) nelebd = nnelebd
      IF(PRESENT(typ_bnd_elem)) typ_bnd_elem = typ_bnd
!
      CALL get_mesh_npoin(fformat,nfic,typ,npoin,ierr)
      CALL check_call(ierr,'READ_MESH_INFO:GET_MESH_NPOIN')
!
      CALL get_mesh_npoin_per_element(fformat,nfic,typ,ndp,ierr)
      CALL check_call(ierr,'READ_MESH_INFO:GET_MESH_NPOIN_PER_ELEMENT')
!
      CALL get_mesh_nplan(fformat,nfic,nplan,ierr)
      CALL check_call(ierr,'READ_MESH_INFO:GET_MESH_NPLAN')
!
      ! Getting the number of boundary points
      ! If we are in parallel partel will have compute boundary on node
      IF(ncsize.GT.1) THEN
        typ_bnd = point_bnd_elt_type
      ENDIF
      CALL get_bnd_npoin(fformat,nfic,typ_bnd,nptfr,ierr)
      CALL check_call(ierr,'READ_MESH_INFO:GET_MESH_NPOIN:NPTFR')
!
      ! TODO: See If nptir is actually used before leclim same for nptfr
      nptir = 0
      IF(ncsize.GT.1) THEN
        CALL get_mesh_nptir(fformat,nfic,nptir,ierr)
        CALL check_call(ierr,'READ_MESH_INFO:GET_MESH_NPTIR')
      ENDIF

      CALL get_mesh_orig(fformat,nfic,x_orig,y_orig,ierr)
      CALL check_call(ierr,'READ_MESH_INFO:GET_MESH_ORIG')
!
!     IF(PRESENT(NELEBD).AND.) THEN
!       CALL GET_MESH_NELEM(FFORMAT,NFIC,TYP_BND,NELEBD,IERR)
!       CALL CHECK_CALL(IERR,'READ_MESH_INFO:GET_MESH_NELEM:BND_ELEM')
!     ENDIF
!
!  PRINTOUT FORMATS:
!
      WRITE(lu,301) title(1:72)
      WRITE(lu,501) nelem,npoin
      IF(title(73:80).EQ.'SERAFIN ') THEN
        WRITE(lu,*) '           SINGLE PRECISION FORMAT (R4)'
      ELSEIF(title(73:80).EQ.'SERAFIND') THEN
        WRITE(lu,*) '           DOUBLE PRECISION FORMAT (R8)'
      ELSE
        WRITE(lu,*) '           FORMAT NOT INDICATED IN TITLE'
      ENDIF
      WRITE(lu,*) ' '
!
      IF(npoin.LT.3) THEN
        WRITE(lu,24) npoin
        CALL plante(1)
        stop
      ENDIF
!
24    FORMAT(1x,'READ_MESH_INFO : NUMBER OF POINTS IN THE MESH: ',
     &       1i9,/,1x,
     &          '           NUMBER OF BOUNDARY POINTS: ',1i8,/,1x,
     &          '           WRONG DATA, PROGRAMME STOPPED')
301   FORMAT(1x,//,1x,'READ_MESH_INFO: TITLE= ',a72)
501   FORMAT(12x,'NUMBER OF ELEMENTS:',1i9,/,
     &       12x,'NUMBER OF POINTS:',1i9,/)
!
!-----------------------------------------------------------------------
!
      RETURN
      END