identify_liq_bnd.f

Go to the documentation of this file.

!                   ************************
                    PROGRAM identify_liq_bnd
!                   ************************
!
!
!***********************************************************************
! IDENTIFY_LIQ_BND
!***********************************************************************
!
!brief    Returns the liquid boundary informations
!
!history  Y. AUDOUIN (EDF)
!+        01/2018
!+        V8P0
!+        FIRST  VERSION
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE declarations_partel
      USE mod_numbering_open_boundaries
      USE bief, ONLY : nbmaxnshare, nptir,
     &                 read_mesh_info, front2, ncsize
      USE interface_hermes
      USE declarations_special
      IMPLICIT NONE
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      CHARACTER(LEN=PATH_LEN) :: NAMEGEO
      CHARACTER(LEN=PATH_LEN) :: NAMECLI
      CHARACTER(LEN=8) :: FFORMAT
!
      CHARACTER(LEN=80)  :: TITLE
      LOGICAL :: IS
      INTEGER NVAR, NPLAN, NPTFR, NELEBD
      INTEGER NELEM, NPOIN, NDP, NELEM2, NPOIN2, NDP_BND
      INTEGER TYP_ELEM, TYP_BND_ELEM
      INTEGER DIM_MESH
!
      INTEGER, ALLOCATABLE :: IKLES(:)
      INTEGER, ALLOCATABLE :: IKLES3D(:)
      INTEGER, ALLOCATABLE :: LIHBOR(:),LIUBOR(:)
      INTEGER, ALLOCATABLE :: LIVBOR(:),LITBOR(:),COLOR(:)
      DOUBLE PRECISION, ALLOCATABLE :: HBOR(:),UBOR(:),VBOR(:)
      DOUBLE PRECISION, ALLOCATABLE :: CHBORD(:)
      DOUBLE PRECISION, ALLOCATABLE :: TBOR(:),ATBOR(:),BTBOR(:)
      INTEGER, ALLOCATABLE :: NBOR(:),IKLE_BND(:)
      INTEGER, ALLOCATABLE :: NUMLIQ(:)
      INTEGER, ALLOCATABLE :: IFABOR(:,:), NELBOR(:)
      INTEGER, ALLOCATABLE :: IKLE(:,:)
      INTEGER, ALLOCATABLE :: KP1BOR(:,:)
!
!     FOR DOUBLE PRECISION SERAFIN FORMAT
!
      DOUBLE PRECISION, ALLOCATABLE    :: F(:,:)
!
      INTEGER :: NINP
      INTEGER :: NOUT
      INTEGER TIME(3),DATE(3), DATE_TMP(6)
      INTEGER :: I,J,K,IERR
      INTEGER :: FIRST, LAST, NFRLIQ
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
!
      lng=2 ! JE NE PARLE PAS FRANCAIS, JE SUIS BARBARIEN
      lu=6  ! FORTRAN STANDARD OUPUT CHANNEL
      li=5  ! FORTRAN STANDARD INPUT CHANNEL
!
!----------------------------------------------------------------------
! INTRODUCE YOURSELF
!
!----------------------------------------------------------------------
! NAMES OF THE INPUT FILES:
!
!
      WRITE(lu,*) '--INPUT GEOMETRY FILE NAME <INPUT_NAME>: '
      READ(li,*) namegeo
!
      IF (namegeo.EQ.' ') THEN
        WRITE (lu,*) ' NO FILENAME '
        CALL plante(1)
        stop
      ELSE
        WRITE(lu,*) 'INPUT: ',trim(namegeo)
      ENDIF

      INQUIRE (file=namegeo,exist=is)
      IF (.NOT.is) THEN
        WRITE (lu,*)' FILE DOES NOT EXIST: ',trim(namegeo)
        CALL plante(1)
        stop
      ENDIF
!
      WRITE(lu,*)
     & '--INPUT FILE FORMAT <FFORMAT> [MED,SERAFIN,SERAFIND]: '
      READ(li,*) fformat
      IF ( (fformat .NE. 'MED     ') .AND.
     &     (fformat(1:7) .NE. 'SERAFIN') ) THEN
        WRITE(lu,*)
     &  ' FILE FORMAT MUST BE "MED" OR "SERAFIN" OR "SERAFIND" '
        CALL plante(1)
        stop
      ELSE
        WRITE(lu,*) ' INPUT: ', fformat
      ENDIF
!
      WRITE(lu, *) '--BOUNDARY CONDITIONS FILE NAME: '
      READ(li,*) namecli
      IF (namecli.EQ.' ') THEN
        WRITE (lu,*) ' NO FILENAME '
        CALL plante(1)
        stop
      ELSE
        WRITE(lu,*) 'INPUT: ',trim(namecli)
      ENDIF
!
      INQUIRE (file=namecli,exist=is)
      IF (.NOT.is) THEN
        WRITE (lu,*) 'FILE DOES NOT EXIST: ',trim(namecli)
        CALL plante(1)
        stop
      ENDIF
!
!     Reading mesh informations
!
      ! SET NCSIZE TO 1 TO USE VOISIN AND READ_MESH_INFO IN SERIAL MODE
      ncsize=1
!
      CALL open_mesh(fformat, namegeo, ninp, 'READ     ', ierr)
      CALL check_call(ierr, 'PARTEL:OPENMESH:INP')
!
      CALL open_bnd(fformat,namecli,ninp,'READ     ',ierr)
      CALL check_call(ierr,'PARTEL:OPEN_BND:NCLI')
!
!----------------------------------------------------------------------
!
! START READING THE GEOMETRY OR RESULT FILE
!
      CALL read_mesh_info(fformat,ninp,title,nvar,npoin,typ_elem,nelem,
     &                    nptfr,nptir,ndp,nplan,x_orig,y_orig,
     &                    typ_bnd_elem,nelebd)
!
! READ THE REST OF THE SELAFIN FILE
! 10 INTEGERS, THE FIRST IS THE NUMBER OF RECORDS (TIMESTEPS)
!
      CALL get_mesh_date(fformat,ninp,date_tmp,ierr)
      CALL check_call(ierr,'PARTEL:GET_MESH_DATE:INP')
      DO i=1,3
        date(i) = date_tmp(i)
        time(i) = date_tmp(i+3)
      ENDDO
!
      IF(nplan.GT.1) THEN
        npoin2 = npoin/nplan
        nelem2 = nelem/(nplan-1)
        dim_mesh = 3
      ELSE
        npoin2 = npoin
        nelem2 = nelem
        dim_mesh = 2
      ENDIF
!
! NOW LET US ALLOCATE
!
      ALLOCATE (ikles(nelem2*3),stat=ierr)
      CALL check_allocate(ierr, 'IKLES')
      IF(nplan.GT.1) THEN
        ALLOCATE (ikles3d(nelem*ndp),stat=ierr)
        CALL check_allocate(ierr, 'IKLES3D')
      ENDIF
!
!     SIZE 3: FIRST TWO FUNCTIONS ARE X AND Y, 3 IS ALL OTHER
!             VARIABLES (THEY WILL BE COPIED AND WRITTEN
!             ONE AFTER THE OTHER...)
!     NPOIN IS 3D HERE IN 3D
!
      ALLOCATE (f(npoin,2),stat=ierr)
      CALL check_allocate(ierr, 'F')
!
! CONNECTIVITY TABLE:
!
      IF(nplan.LE.1) THEN
        CALL get_mesh_connectivity(fformat,ninp,typ_elem,ikles,
     &                             nelem,ndp,ierr)
        CALL check_call(ierr,'PARTEL:GET_MESH_CONNECTIVITY:2D')
      ELSE
        CALL get_mesh_connectivity(fformat,ninp,typ_elem,ikles3d,
     &                             nelem,ndp,ierr)
        CALL check_call(ierr,'PARTEL:GET_MESH_CONNECTIVITY:3D')
!       BUILDING IKLES
        DO j=1,3
          DO k=1,nelem2
            ikles((k-1)*3+j)=ikles3d((k-1)*6+j)
          ENDDO
        ENDDO
      ENDIF
!
! X-, Y-COORDINATES
!
      CALL get_mesh_coord(fformat,ninp,1,2,npoin,f(:,1),ierr)
      CALL check_call(ierr,'PARTEL:GET_MESH_COORD:X')
      CALL get_mesh_coord(fformat,ninp,2,2,npoin,f(:,2),ierr)
      CALL check_call(ierr,'PARTEL:GET_MESH_COORD:Y')
!
!----------------------------------------------------------------------
!
!     READ THE BOUNDARY CONDITIONS FILE
!
      !
      CALL get_nodes_per_element(typ_bnd_elem,ndp_bnd)
      ! GET THE NUMBER OF BOUNDARY POINTS AND ELEMENTS
      CALL get_bnd_nelem(fformat,ninp,typ_bnd_elem,nelebd,ierr)
      CALL check_call(ierr,'PARTEL:GET_BND_NELEBD:NCLI')
      CALL get_bnd_npoin(fformat,ninp,typ_bnd_elem,nptfr,ierr)
      CALL check_call(ierr,'PARTEL:GET_BND_NPOIN:NCLI')
      !
      ALLOCATE(ikle_bnd(nelebd*ndp_bnd),stat=ierr)
      CALL check_allocate(ierr,'PARTEL:IKLE_BND')
      ALLOCATE(nbor(nptfr),stat=ierr)
      CALL check_allocate(ierr,'PARTEL:NBOR')
      ! ALLOCATING ARRAY FOR THE BOUNDARY CONDITIONS
      ALLOCATE(lihbor(nptfr),stat=ierr)
      CALL check_allocate(ierr,'PARTEL:LIHBOR')
      ALLOCATE(liubor(nptfr),stat=ierr)
      CALL check_allocate(ierr,'PARTEL:LIUBOR')
      ALLOCATE(livbor(nptfr),stat=ierr)
      CALL check_allocate(ierr,'PARTEL:LIVBOR')
      ALLOCATE(hbor(nptfr),stat=ierr)
      CALL check_allocate(ierr,'PARTEL:HBOR')
      ALLOCATE(ubor(nptfr),stat=ierr)
      CALL check_allocate(ierr,'PARTEL:UBOR')
      ALLOCATE(vbor(nptfr),stat=ierr)
      CALL check_allocate(ierr,'PARTEL:VBOR')
      ALLOCATE(chbord(nptfr),stat=ierr)
      CALL check_allocate(ierr,'PARTEL:CHBORD')
      ALLOCATE(litbor(nptfr),stat=ierr)
      CALL check_allocate(ierr,'PARTEL:LITBOR')
      ALLOCATE(tbor(nptfr),stat=ierr)
      CALL check_allocate(ierr,'PARTEL:TBOR')
      ALLOCATE(atbor(nptfr),stat=ierr)
      CALL check_allocate(ierr,'PARTEL:ATBOR')
      ALLOCATE(btbor(nptfr),stat=ierr)
      CALL check_allocate(ierr,'PARTEL:BTBOR')
      ALLOCATE (numliq(nptfr),stat=ierr)
      CALL check_allocate(ierr, 'NUMLIQ')
      ALLOCATE (color(nptfr),stat=ierr)
      CALL check_allocate(ierr, 'COLOR')
      ! Get the connectivity table for the boundary elements
      CALL get_bnd_connectivity(fformat,ninp,typ_bnd_elem,nelebd,
     &                    ndp_bnd,ikle_bnd,ierr)
      CALL check_call(ierr,'PARTEL:GET_BND_CONNECTIVITY:NCLI')
      ! FILL NBOR
      CALL get_bnd_numbering(fformat,ninp,typ_bnd_elem,nptfr,nbor,ierr)
      ! GET THE VALUE OF EACH BOUNDARY
      CALL get_bnd_value(fformat,ninp,typ_bnd_elem,nelebd,
     &                   lihbor,liubor,
     &                   livbor,hbor,ubor,vbor,chbord,.true.,
     &                   litbor,tbor,atbor,btbor,nptfr,ierr)
      CALL check_call(ierr,'PARTEL:GET_BND_VALUE:NCLI')
      CALL get_bnd_color(fformat,ninp,typ_bnd_elem,nelebd,
     &                   color,ierr)
      CALL check_call(ierr,'PARTEL:GET_BND_COLOR:NCLI')
!
!     Computing liquid bounndaries info
!
      CALL numbering_open_boundaries
     & (namegeo, ikle, ikles,
     &  kp1bor, numliq, dim_mesh, npoin2, nptfr, npoin, nelem2,
     &  nelbor, liubor, lihbor, nbor, ifabor, f, .true.)

      ! Closing the files
      CALL close_bnd(fformat,ninp,ierr)
      CALL check_call(ierr,'PARTEL:CLOSE_BND:NCLI')
      CALL close_mesh(fformat,ninp,ierr)
      CALL check_call(ierr,'PARTEL:CLOSE_MESH:NINP')
!
!     Output informations in file
!
      CALL get_free_id(nout)
      OPEN(file="liq_bnd.txt", action='WRITE', unit=nout,
     &     form='FORMATTED', iostat=ierr)
      CALL check_call(ierr, 'opening liq_bnd.txt')
      ! Finding the first liquid boundary
      nfrliq = maxval(numliq)
      WRITE(nout,*) nfrliq
      DO i=1,nfrliq
!
!  MARKS THE NUMBERS OF THE LIQUID BOUNDARIES
!
        ! LOKKING FOR FIRST POINT
        DO j=1,nptfr
          IF(numliq(j).NE.i .AND. numliq(mod(j,nptfr)+1).EQ.i) THEN
            first = mod(j,nptfr)+1
            EXIT
          ENDIF
        ENDDO
        ! LOOKING FOR LAST POINT
        DO j=1,nptfr
          IF(numliq(j).EQ.i .AND. numliq(mod(j,nptfr)+1).NE.i) THEN
            last = j
            EXIT
          ENDIF
        ENDDO
        WRITE(nout,*) i, first, nbor(first),
     &                f(nbor(first),1), f(nbor(first),2),
     &                last, nbor(last),
     &                f(nbor(last),1), f(nbor(last),2)
      ENDDO

      CLOSE(nout)
!
!     Deallocation
!
      DEALLOCATE(ikles)
      IF(nplan.GT.1) THEN
        DEALLOCATE(ikles3d)
      ENDIF
      DEALLOCATE(f)
      DEALLOCATE(ikle_bnd)
      DEALLOCATE(nbor)
      DEALLOCATE(lihbor)
      DEALLOCATE(liubor)
      DEALLOCATE(livbor)
      DEALLOCATE(hbor)
      DEALLOCATE(ubor)
      DEALLOCATE(vbor)
      DEALLOCATE(chbord)
      DEALLOCATE(litbor)
      DEALLOCATE(tbor)
      DEALLOCATE(atbor)
      DEALLOCATE(btbor)
      DEALLOCATE (numliq)
      DEALLOCATE (color)
      stop 0
      END PROGRAM