 SUBROUTINE read_gmsh(fichmsh,calc)
!=========================== DOCUMENTATION ============================= 
!DOCDEBh*    UTI/PRE/read_gmsh  *
!
!    NOM
! read_gmsh(ficmsh,msh)
!    DESCRIPTION  
!*
! Lecture d un fichier de maillage au format Gmsh et affectation dans les blocs au format CERF
!
! Read a mesh file in Gmsh format and assign it to blocks in CERF format
!*
! 
!    ENTREES / INPUT
! fichmsh: nom de fichier de maillage au format Gmsh v2.2 / Mesh file name in Gmsh v2.2 format
!
!    SORTIES / OUTPUT
! calc: objet calcul modifie / modified calculation object
!
!    SOURCE
! Ce logiciel est regi par la licence [CeCILL-B]
! This software is governed by the [CeCILL-B] license
!=========================== DEBUT DES DECLARATIONS ==================== 
!.1-----  Implicit, Use
  use num
  use zorder
  use msh
  implicit none
!.2-----  Declaration
  type(calcul), intent(inout) :: calc !! Objet calcul / Calcul object
  character( *), intent(in)   :: fichmsh !! Nom du fichier de maillage / Mesh file name
  
  integer, dimension(:,:), allocatable :: iconnec
  integer, dimension(15)               :: nbn=(/2,3,4,4,8,1,1,1,1,1,1,1,1,1,1/) ! nbre de noeud en fct du type d elt
  real(kind=kind(0.d0))                :: x,y,z,d
  integer                              :: in,mortonx,mortony,mortonz,ie,ier,nb_elem,jn,itemp,iphys,ityp,j,nbelem,nbf,ii
  integer                              :: if1,if2,ic1,ic2,k
  integer(kind=16)                     :: nmorton
  integer(kind=16) ,allocatable,dimension(:)     :: lm
  character*80                         :: ligne
  type(chaine_MF), pointer             :: lcf,dblechf,totof
  type(vertex)                         :: v
  type(mface)                          :: f
  type(noe), dimension(:), allocatable :: lst_noe
  integer                              :: nb_noe
  type(noe),dimension(8)               :: ln,lnf
  type(noe)                            :: noe_null
  type(bloc)                           :: b_null
!=========================== DEBUT DU CODE EXECUTABLE ==================
!
!-------  Initialisations
!

noe_null%coord%x= 0.d0 ; noe_null%coord%y= 0.d0 ; noe_null%coord%z= 0.d0
noe_null%morton= 0 ; noe_null%num= 0
b_null%numbl= 0 ; b_null%idom= 0 ;b_null%morton= 0 ;b_null%nrb= 0 ;b_null%nx= 1 ;b_null%ny= 1 ;b_null%nz= 1
b_null%s(:)%x= 0.d0 ; b_null%s(:)%y= 0.d0 ; b_null%s(:)%z= 0.d0 ; b_null%volc= 0.d0
b_null%cl(:)= 0 ; b_null%nrf(:)= 0 ; b_null%dom(:)= 0
b_null%first_cell= 0 ; b_null%nb_cell= 0
b_null%first_face= 0 ; b_null%nb_face= 0
b_null%first_vertex= 0 ; b_null%nb_vertex= 0 ; b_null%isol= .false.
    allocate(iconnec(1,1))
    iconnec = 0
    allocate(lst_noe(1))
    lst_noe(:)=noe_null


open(11,file=fichmsh,status='old')
10    continue 
read(11,'(a)', end=99) ligne 
! lecture des noeuds
if (ligne(1:6) == '$Nodes') then
    read(11,*) nb_noe
    write(*,*)'read_gmsh-INFO Reading ',nb_noe,' vertexes'
    deallocate(lst_noe)
    allocate(lst_noe(nb_noe))
    lst_noe(:)=noe_null
    x_min_morton=1.d20
    y_min_morton=1.d20
    z_min_morton=1.d20
    do in = 1,nb_noe
        read(11,*) jn,x,y,z
        lst_noe(in)%coord%x = x
        lst_noe(in)%coord%y = y
        lst_noe(in)%coord%z = z
        lst_noe(in)%num = in
        if (x < x_min_morton) x_min_morton = x
        if (y < y_min_morton) y_min_morton = y
        if (z < z_min_morton) z_min_morton = z
    enddo
endif

! lecture de la connectivite des elements
if (ligne(1:9) == '$Elements') then
    read(11,*) nbelem
    write(*,*)'read_gmsh-INFO Reading ',nbelem,' elements'
    deallocate(iconnec)
    allocate(iconnec(nbelem,11))
    iconnec = 0
    do ie = 1,nbelem
        read(11,*) jn,ityp,in,iphys,(itemp,ii=1,in-1),(iconnec(ie,ii+2),ii=1,nbn(ityp))
        iconnec(ie,1) = ityp !ityp=1,2,3,4,5,6 -> l2,t3,q4,th4,h8,point
        iconnec(ie,2) = iphys  !propriete physique
    enddo
endif

goto 10 

99    continue 
close(11)

! Estimation de petit_morton et du nombre d elements
petit_morton=1.d20
nb_elem = 0
do ie=1,nbelem
    ! nbn(iconnec(ie,1)) : ityp
    if (iconnec(ie,1) == 2 ) cycle ! t3
    if (iconnec(ie,1) == 3 ) cycle ! q4
    if (iconnec(ie,1) == 4 ) cycle ! th4
    if (iconnec(ie,1) == 5)then
        nb_elem = nb_elem + 1 ! h8
        iconnec(ie,11)=nb_elem
    endif
    if (iconnec(ie,1) >= 6) cycle ! on boucle
    do j=1,nbn(iconnec(ie,1))
        do k=1,nbn(iconnec(ie,1))
            if (k /= j) then
                d=sqrt((lst_noe(iconnec(ie,k+2))%coord%x-lst_noe(iconnec(ie,j+2))%coord%x)**2 + &
                       (lst_noe(iconnec(ie,k+2))%coord%y-lst_noe(iconnec(ie,j+2))%coord%y)**2 + &
                       (lst_noe(iconnec(ie,k+2))%coord%z-lst_noe(iconnec(ie,j+2))%coord%z)**2)
                if (d < petit_morton) petit_morton = d
            endif
        enddo
    enddo
enddo
petit_morton = petit_morton/20.d0
write(*,*)'read_gmsh-INFO Small_Morton:',petit_morton


 ! Determination du code Morton des noeuds
do in=1,nb_noe
    mortonx = coord_real2int(lst_noe(in)%coord%x,x_min_morton,petit_morton)
    mortony = coord_real2int(lst_noe(in)%coord%y,y_min_morton,petit_morton)
    !if (calc%msh%ndim == 2) then
    !    call zordercurve2d(mortonx,mortony,nmorton)
    !else
        mortonz = coord_real2int(lst_noe(in)%coord%z,z_min_morton,petit_morton)
        call zordercurve3d(mortonx,mortony,mortonz,nmorton)
    !endif
    lst_noe(in)%morton = nmorton
enddo

! Initialisations de la liste chainee des faces
allocate(lcf)
nullify(lcf%suiv); nullify(lcf%prec)
lcf%f%morton = 0
lcf%f%num = 0

! On range les faces  dans la liste chainee  et on range les blocs
calc%nb_bloc=nb_elem
allocate(calc%list_bloc(calc%nb_bloc))
allocate(lm(calc%nb_bloc))
nb_elem = 0
do ie = 1,nbelem
    select case (iconnec(ie,1))
        case(1,15)
        !~case(1) ! l2 pour conserver le numero de propriete des frontieres  --------------------
        !~    if (calc%msh%ndim == 2) then
        !~        call ini_mface(f)
        !~        lnf(1)=lst_noe(iconnec(ie,3))
        !~        lnf(2)=lst_noe(iconnec(ie,4)) 
        !~        call mortonface(f,lnf,2,calc%msh%ndim)
        !~        f%be2=-iconnec(ie,2)
        !~        call range_mface(f,lcf,ier)
        !~    endif

        case(3) ! q4 ---------------------------------------------------------------------
            !~if (calc%msh%ndim == 2) then ! alors c est un element
                !~nb_elem=nb_elem+1
                !~calc%list_bloc(nb_elem)=b_null
                !~ln(:)=noe_null
                !~do j=1,nbn(iconnec(ie,1))
                !~    ln(j)=lst_noe(iconnec(ie,j+2))
                !~enddo
                !~call permnoe(ln,nbn(iconnec(ie,1))) ! on permute les noe pour demarrer du plus petit morton
                !~call noe2bloc(ln,calc%list_bloc(nb_elem),nbn(iconnec(ie,1)),2)
                !~do j=1,nbn(iconnec(ie,1))
                !~    lnf(:)=noe_null
                !~    lnf(2)=ln(j)
                !~    if (j /= nbn(iconnec(ie,1))) then
                !~        lnf(1)=ln(j+1)
                !~    else
                !~        lnf(1)=ln(1)
                !~    endif
                !~    call ini_mface(f)
                !~    call mortonface(f,lnf,2,2)
                !~    f%be1=nb_elem
                !~    f%fac_loc(1)=j
                !~    call perm4(f%mv,nbn(iconnec(ie,1)))
                !~    call range_mface(f,lcf,ier)
                !~enddo
                !~calc%list_bloc(nb_elem)%numbl   = nb_elem
            !~else ! sinon c est une face
                v%x = 0.d0; v%y = 0.d0; v%z = 0.d0
                call ini_mface(f) 
                lnf(:) = noe_null
                do j = 1,nbn(iconnec(ie,1))
                    lnf(j) = lst_noe(iconnec(ie,j+2))
                enddo
                call permnoe(lnf,nbn(iconnec(ie,1))) ! on permute les noe pour demarrer du plus petit morton
                call mortonface(f,lnf,nbn(iconnec(ie,1)),3)
                f%be2 = -iconnec(ie,2)
                !call affiche_mface(f)
                call range_mface(f,lcf,ier)
            !~endif
            
        case(5) ! h8 -------------------------------------------------------------------------------
            nb_elem=nb_elem+1
            ! bloc
            calc%list_bloc(nb_elem)=b_null
            ln(:)=noe_null
            do j=1,nbn(iconnec(ie,1))
                ln(j)=lst_noe(iconnec(ie,j+2))
            enddo
            !call permnoe(ln,8) ! on permute les noe pour demarrer du plus petit morton
            call noe2bloc(ln,calc%list_bloc(nb_elem),nbn(iconnec(ie,1)),3,lm(nb_elem))
            calc%list_bloc(nb_elem)%numbl   = nb_elem

            ! face 1 5 8 4
            call ini_mface(f)
            lnf(:)=noe_null ; lnf(1)=ln(1) ; lnf(2)=ln(5) ; lnf(3)=ln(8) ; lnf(4)=ln(4)
            call mortonface(f,lnf,4,3)
            f%be1=nb_elem
            f%fac_loc(1)=1
            call perm4(f%mv,4)
            call range_mface(f,lcf,ier)
            ! face 2 3 7 6
            call ini_mface(f)
            lnf(:)=noe_null ; lnf(1)=ln(2) ; lnf(2)=ln(3) ; lnf(3)=ln(7) ; lnf(4)=ln(6)
            call mortonface(f,lnf,4,3)
            f%be1=nb_elem
            f%fac_loc(1)=2
            call perm4(f%mv,4)
            call range_mface(f,lcf,ier)
            ! face 1 2 6 5
            call ini_mface(f)
            lnf(:)=noe_null ; lnf(1)=ln(1) ; lnf(2)=ln(2) ; lnf(3)=ln(6) ; lnf(4)=ln(5)
            call mortonface(f,lnf,4,3)
            f%be1=nb_elem
            f%fac_loc(1)=3
            call perm4(f%mv,4)
            call range_mface(f,lcf,ier)
            ! face 3 4 8 7
            call ini_mface(f)
            lnf(:)=noe_null ; lnf(1)=ln(3) ; lnf(2)=ln(4) ; lnf(3)=ln(8) ; lnf(4)=ln(7)
            call mortonface(f,lnf,4,3)
            f%be1=nb_elem
            f%fac_loc(1)=4
            call perm4(f%mv,4)
            call range_mface(f,lcf,ier)
            ! face 1 4 3 2
            call ini_mface(f)
            lnf(:)=noe_null ; lnf(1)=ln(1) ; lnf(2)=ln(4) ; lnf(3)=ln(3) ; lnf(4)=ln(2)
            call mortonface(f,lnf,4,3)
            f%be1=nb_elem
            f%fac_loc(1)=5
            call perm4(f%mv,4)
            call range_mface(f,lcf,ier)
            ! face 5 6 7 8
            call ini_mface(f)
            lnf(:)=noe_null ; lnf(1)=ln(5) ; lnf(2)=ln(6) ; lnf(3)=ln(7) ; lnf(4)=ln(8)
            call mortonface(f,lnf,4,3)
            f%be1=nb_elem
            f%fac_loc(1)=6
            call perm4(f%mv,4)
            call range_mface(f,lcf,ier)
            calc%list_bloc(nb_elem)%volc=volbloc(calc%list_bloc(nb_elem))
            
        case default
            write(*,*)' INFO-read_gmsh, type element',iconnec(ie,1),' non prevu!'
    end select
enddo

!---------------------------------------------------------
! On parcourt les faces pour determiner les voisins des blocs
!
nbf = 0
dblechf=>lcf
do while (associated(dblechf%prec))
    dblechf=>dblechf%prec
enddo
dblechf=>dblechf%suiv  ! on saute le premier qui est nul
do while (associated(dblechf))
    nbf=nbf+1
    dblechf%f%num=nbf
    dblechf=>dblechf%suiv
enddo
 
nbf= 0
dblechf=>lcf
do while (associated(dblechf%prec))
    dblechf=>dblechf%prec
enddo
dblechf=>dblechf%suiv  ! on saute le premier qui est nul
do while (associated(dblechf))
    if (dblechf%f%num > 0) then
        nbf=nbf+1   
        ic1 = dblechf%f%be1
        ic2 = dblechf%f%be2
        if1 = dblechf%f%fac_loc(1)
        if2 = dblechf%f%fac_loc(2)
        calc%list_bloc(ic1)%cl(if1)  =ic2
        if (ic2 > 0) then
            calc%list_bloc(ic2)%cl(if2)  =ic1
        endif
    endif
    dblechf=>dblechf%suiv
enddo
    write(*,*)'read_gmsh-INFO creation of ',nbf,' faces'

! on ordonne les numeros (locaux) des blocs selon le code de morton (lm)
k=0
do ii=1,calc%nb_bloc
    j=int(minloc(lm,1,mask=lm.ne.0))
    k=k+1
    calc%list_bloc(ii)%morton    = k
    lm(j)=0
enddo
deallocate(lm)
    write(*,*)'read_gmsh-INFO creation of ',calc%nb_bloc,' master elements'

!----------------------------------------------------------------------------
! LIBERATION DE LA MEMOIRE
!----------------------------------------------------------------------------
dblechf=>lcf
do while (associated(dblechf%prec))
    dblechf => dblechf%prec
enddo
do while (associated(dblechf))
    totof => dblechf
    dblechf => dblechf%suiv
    deallocate(totof)
enddo
!
!===========================   FIN DE LA ROUTINE    ====================
END SUBROUTINE read_gmsh
!DOCFIN 

