read_gmsh.f90

### TABLE OF CONTENTS

- 1. UTI/PRE/read_gmsh

## UTI/PRE/read_gmsh [ Modules ]

[ Top ] [ Modules ]

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