block2mesh.f90

### TABLE OF CONTENTS

- 1. UTI/PRE/block2mesh

## UTI/PRE/block2mesh [ Modules ]

[ Top ] [ Modules ]

NOM

     block2mesh(cb,nzone,zone,nzonei,zonei)

DESCRIPTION

     A partir du maillage maitre (les blocs) on construit le maillage de chaque domaine.
     From the master mesh (the blocks) we construct the mesh for each domain.

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 uti, only: inter,init_def_bloc,ini_w,splitdom_Z ! ,bathy_sv1,initsu_sv1
      IMPLICIT NONE
    !.2-----  Declaration
      type(calcul)                                      :: cb !! Objet calcul contenant le maillage maitre / Calcul object containing the master mesh
      integer                                           :: nzonei
    !! Parametres de definition des maillages et conditions aux limites / Mesh definition parameters and boundary conditions
      real(kind=kind(0.d0)),allocatable,dimension(:,:)  :: zone,zonei
    !! Parametres de definition des maillages et conditions aux limites / Mesh definition parameters and boundary conditions
      integer,allocatable,dimension(:,:)                :: nzone
    !! Parametres de definition des maillages et conditions aux limites / Mesh definition parameters and boundary conditions

      type(bloc)                                        :: blv(6),bl_null
      type(calcul)                                      :: calc
      integer                                           :: imac(6),imav(6),imaf(6)
      integer                                           :: idom,nb_dom,nx,ny,nz,i,nb_bloc,ib
      integer                                           :: icb,ibv
      integer                                           :: k,nvi
      integer                                           :: l
      real(kind=kind(0.d0))                             :: hmin,dd,x,y
      type(list_int), allocatable,dimension(:)          :: interdom
      integer, allocatable,dimension(:)                 :: nbl

    !=========================== DEBUT DU CODE EXECUTABLE ==================
    bl_null%numbl =0 ;bl_null%idom =0 ;bl_null%morton  =0 ;bl_null%nrb =0 ;bl_null%nx =0 ;bl_null%ny =0 ;bl_null%nz  =0 ;
    bl_null%volc  =0.d0 ;bl_null%cl(:)  =0 ;bl_null%nrf(:)  =0 ;bl_null%dom(:)  =0 ;bl_null%first_cell =0 ;
    bl_null%nb_cell   =0 ;   bl_null%first_face =0 ;bl_null%nb_face    =0 ;  bl_null%first_vertex =0 ;bl_null%nb_vertex  =0

      write(*,*)
      write(*,*)achar(27)//'[34m================================================'
      write(*,*)'Creation of the domain mesh from the master mesh'
      write(*,*)'================================================'//achar(27)//'[0m'
      write(*,*)
      nb_bloc=cb%nb_bloc
      nb_dom=cb%nbds
      write(*,*)'BLOCK2MESH-INFO Problem dimension :',cb%msh%ndim
    ! estimation de hmin
        hmin=1.d20
        DO ib=1,nb_bloc
            dd=sqrt((cb%list_bloc(ib)%s(1)%x-cb%list_bloc(ib)%s(2)%x)**2+ &
                    (cb%list_bloc(ib)%s(1)%y-cb%list_bloc(ib)%s(2)%y)**2+(cb%list_bloc(ib)%s(1)%z-cb%list_bloc(ib)%s(2)%z)**2)
            if(dd.lt.hmin)hmin=dd
            dd=sqrt((cb%list_bloc(ib)%s(1)%x-cb%list_bloc(ib)%s(4)%x)**2+ &
                    (cb%list_bloc(ib)%s(1)%y-cb%list_bloc(ib)%s(4)%y)**2+(cb%list_bloc(ib)%s(1)%z-cb%list_bloc(ib)%s(4)%z)**2)
            if(dd.lt.hmin)hmin=dd
            dd=sqrt((cb%list_bloc(ib)%s(1)%x-cb%list_bloc(ib)%s(5)%x)**2+ &
                    (cb%list_bloc(ib)%s(1)%y-cb%list_bloc(ib)%s(5)%y)**2+(cb%list_bloc(ib)%s(1)%z-cb%list_bloc(ib)%s(5)%z)**2)
            if(dd.lt.hmin)hmin=dd
        END DO

    !----------------------------------------------------------------------------------
    ! on definit la discretisation des blocs et l appartenance a un domaine
    !
      CALL init_def_bloc(cb,nzone,zone)

    !----------------------------------------------------------------------------------
    ! on verifie la condition sur le niveau de discretisation de 2 blocs voisins
      DO l=1,cb%nrma
        DO ib=1,nb_bloc
            DO i=1,6
                IF (cb%list_bloc(ib)%cl(i).gt.0) THEN
                   ibv=cb%list_bloc(ib)%cl(i)
                   IF(cb%list_bloc(ib)%nrb.eq.(cb%nrma+1-l).and.cb%list_bloc(ibv)%nrb.lt.(cb%nrma-l)) cb%list_bloc(ibv)%nrb=(cb%nrma-l)
                END IF
            END DO
        END DO
      END DO
    !----------------------------------------------------------------------------------
    ! on decompose en sous-domaines par code de Morton

      write( *,1011)nb_bloc,nb_dom
    1011 format('Aggregation of ',i6,' blocks into',i3,' sub-domains')

      CALL splitdom_Z(cb%list_bloc,nb_bloc,nb_dom, cb%msh%ndim)

    ! liste des voisins: -1/ le voisin est + grossier
    !                     0/ le voisin est identique
    !                     1/ le voisin est + fin
      DO ib=1,nb_bloc
          cb%list_bloc(ib)%nrf=0
          DO i=1,6
              ibv=cb%list_bloc(ib)%cl(i)
              IF(ibv.lt.-100000)ibv=-ibv-100000
              IF (ibv.gt.0) THEN
                  IF (cb%list_bloc(ibv)%nrb.gt.cb%list_bloc(ib)%nrb) cb%list_bloc(ib)%nrf(i)=1
                  IF (cb%list_bloc(ibv)%nrb.lt.cb%list_bloc(ib)%nrb) cb%list_bloc(ib)%nrf(i)=-1
              END IF
          END DO
      END DO
    !
    !-------  Determination de la taille des tableaux du bloc
      DO ib=1,nb_bloc
          cb%list_bloc(ib)%dom=0
          DO i=1,6
              ibv=cb%list_bloc(ib)%cl(i)
              IF(ibv.lt.-100000)ibv=-ibv-100000
              IF (ibv.gt.0) cb%list_bloc(ib)%dom(i)=cb%list_bloc(ibv)%idom
          END DO
          imac=1
          imav=1
          imaf=1
          DO i=1,6 ! ima: indice de maillage de l interface entre 2 blocs voisins
              ibv=cb%list_bloc(ib)%cl(i)
              IF(ibv.lt.-100000)ibv=-ibv-100000
              IF (cb%list_bloc(ib)%idom.eq.cb%list_bloc(ib)%dom(i).and.ib.lt.ibv) imaf(i)=0
              IF (cb%list_bloc(ib)%idom.eq.cb%list_bloc(ib)%dom(i).and.ib.lt.ibv) imav(i)=0
              IF (cb%list_bloc(ib)%idom.eq.cb%list_bloc(ib)%dom(i).or.ibv.le.0)   imac(i)=0
          END DO

          SELECT CASE (cb%msh%ndim)
            CASE (1) ! Cas 1D
              nx=cb%list_bloc(ib)%nx*2**(cb%list_bloc(ib)%nrb)
              ny=1
              nz=1
              cb%list_bloc(ib)%nb_cell   = nx+imac(1)+imac(2)
              cb%list_bloc(ib)%nb_vertex = (nx+1)*(ny+1)*2
              cb%list_bloc(ib)%nb_face   = (nx-1)+nx*2*2 + imaf(1)+imaf(2)

             CASE (2) ! Cas 2D
              nx=cb%list_bloc(ib)%nx*2**(cb%list_bloc(ib)%nrb)
              ny=cb%list_bloc(ib)%ny*2**(cb%list_bloc(ib)%nrb)
              nz=1

              cb%list_bloc(ib)%nb_cell   = nx*ny + int(ny*2.**cb%list_bloc(ib)%nrf(1)*imac(1)) &
                                       + int(ny*2.**cb%list_bloc(ib)%nrf(2)*imac(2)) &
                                       + int(nx*2.**cb%list_bloc(ib)%nrf(3)*imac(3)) &
                                       + int(nx*2.**cb%list_bloc(ib)%nrf(4)*imac(4))
              cb%list_bloc(ib)%nb_vertex = (nx+1)*(ny+1)*2 &
                                       + cb%list_bloc(ib)%nrf(1)*(cb%list_bloc(ib)%nrf(1)+1)*ny*imav(1) &
                                       + cb%list_bloc(ib)%nrf(2)*(cb%list_bloc(ib)%nrf(2)+1)*ny*imav(2) &
                                       + cb%list_bloc(ib)%nrf(3)*(cb%list_bloc(ib)%nrf(3)+1)*nx*imav(3) &
                                       + cb%list_bloc(ib)%nrf(4)*(cb%list_bloc(ib)%nrf(4)+1)*nx*imav(4)
              cb%list_bloc(ib)%nb_face   = (nx-1)*ny+nx*(ny-1)+nx*ny*2 &
                                       + ny*imaf(1)+ny*imaf(2)+nx*imaf(3)+nx*imaf(4)  &
                                       + cb%list_bloc(ib)%nrf(1)*(cb%list_bloc(ib)%nrf(1)+1)*ny/2*imaf(1) &
                                       + cb%list_bloc(ib)%nrf(2)*(cb%list_bloc(ib)%nrf(2)+1)*ny/2*imaf(2) &
                                       + cb%list_bloc(ib)%nrf(3)*(cb%list_bloc(ib)%nrf(3)+1)*nx/2*imaf(3) &
                                       + cb%list_bloc(ib)%nrf(4)*(cb%list_bloc(ib)%nrf(4)+1)*nx/2*imaf(4)
             CASE (3) ! Cas 3D
              nx=cb%list_bloc(ib)%nx*2**(cb%list_bloc(ib)%nrb)
              ny=cb%list_bloc(ib)%ny*2**(cb%list_bloc(ib)%nrb)
              nz=cb%list_bloc(ib)%nz*2**(cb%list_bloc(ib)%nrb)
              cb%list_bloc(ib)%nb_cell   = nx*ny*nz &
                                       + int(ny*nz*4.**cb%list_bloc(ib)%nrf(1)*imac(1)) &
                                       + int(ny*nz*4.**cb%list_bloc(ib)%nrf(2)*imac(2)) &
                                       + int(nx*nz*4.**cb%list_bloc(ib)%nrf(3)*imac(3)) &
                                       + int(nx*nz*4.**cb%list_bloc(ib)%nrf(4)*imac(4)) &
                                       + int(nx*ny*4.**cb%list_bloc(ib)%nrf(5)*imac(5)) &
                                       + int(nx*ny*4.**cb%list_bloc(ib)%nrf(6)*imac(6))
              cb%list_bloc(ib)%nb_vertex = (nx+1)*(ny+1)*(nz+1) &
                                       + cb%list_bloc(ib)%nrf(1)*(cb%list_bloc(ib)%nrf(1)+1)/2*(3*ny*nz+ny+nz)*imav(1) &
                                       + cb%list_bloc(ib)%nrf(2)*(cb%list_bloc(ib)%nrf(2)+1)/2*(3*ny*nz+ny+nz)*imav(2) &
                                       + cb%list_bloc(ib)%nrf(3)*(cb%list_bloc(ib)%nrf(3)+1)/2*(3*nx*nz+nx+nz)*imav(3) &
                                       + cb%list_bloc(ib)%nrf(4)*(cb%list_bloc(ib)%nrf(4)+1)/2*(3*nx*nz+nx+nz)*imav(4) &
                                       + cb%list_bloc(ib)%nrf(5)*(cb%list_bloc(ib)%nrf(5)+1)/2*(3*nx*ny+nx+ny)*imav(5) &
                                       + cb%list_bloc(ib)%nrf(6)*(cb%list_bloc(ib)%nrf(6)+1)/2*(3*nx*ny+nx+ny)*imav(6)
              cb%list_bloc(ib)%nb_face   = (nx-1)*ny*nz+nx*(ny-1)*nz+nx*ny*(nz-1) &
                                       + ny*nz*imaf(1)+ny*nz*imaf(2)+nx*nz*imaf(3)+nx*nz*imaf(4) &
                                       + nx*ny*imaf(5)+nx*ny*imaf(6) &
                                       + cb%list_bloc(ib)%nrf(1)*(cb%list_bloc(ib)%nrf(1)+1)/2*3*ny*nz*imaf(1) &
                                       + cb%list_bloc(ib)%nrf(2)*(cb%list_bloc(ib)%nrf(2)+1)/2*3*ny*nz*imaf(2) &
                                       + cb%list_bloc(ib)%nrf(3)*(cb%list_bloc(ib)%nrf(3)+1)/2*3*nx*nz*imaf(3) &
                                       + cb%list_bloc(ib)%nrf(4)*(cb%list_bloc(ib)%nrf(4)+1)/2*3*nx*nz*imaf(4) &
                                       + cb%list_bloc(ib)%nrf(5)*(cb%list_bloc(ib)%nrf(5)+1)/2*3*nx*ny*imaf(5) &
                                       + cb%list_bloc(ib)%nrf(6)*(cb%list_bloc(ib)%nrf(6)+1)/2*3*nx*ny*imaf(6)

          END SELECT
          !call display(cb%list_bloc(ib)) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      END DO

    ! Creation de la liste des interfaces
      allocate(interdom(int(nb_dom*(nb_dom-1)/2)))
      allocate(nbl(int(nb_dom*(nb_dom-1)/2)))
      CALL inter(cb%list_bloc,nb_bloc,nb_dom,interdom,nbl,cb%msh%ndim)

    ! Determination de petit
      petit=1.d-3
      DO ib=1,nb_bloc
          IF (hmin/(cb%list_bloc(ib)%nx*2**cb%list_bloc(ib)%nrb).lt.petit) &
              petit=hmin/(cb%list_bloc(ib)%nx*2**cb%list_bloc(ib)%nrb)
          IF (hmin/(cb%list_bloc(ib)%ny*2**cb%list_bloc(ib)%nrb).lt.petit.and.cb%msh%ndim.ne.1) &
              petit=hmin/(cb%list_bloc(ib)%ny*2**cb%list_bloc(ib)%nrb)
          IF (hmin/(cb%list_bloc(ib)%nz*2**cb%list_bloc(ib)%nrb).lt.petit.and.cb%msh%ndim.eq.3) &
              petit=hmin/(cb%list_bloc(ib)%nz*2**cb%list_bloc(ib)%nrb)
      END DO
      petit=petit*1.d-5
      write( *,1212)petit
     1212  format('Small dimension',e13.5)

     !----------- Cas particulier de Saint Venant -----
    ! On stocke eventuellement la bathy sur chaque bloc
    if(imodel_phy.eq.1.and.ibathy.eq.2)then
    open(14,file=fich_bathy,status='old')
       do ib=1,nb_bloc
           do i=1,4
              read(14,*)x,y,cb%list_bloc(ib)%vb(i), cb%list_bloc(ib)%vb(i+4)
           enddo
       enddo
       close(14)
    endif

    !-----------------------------------------------------------------------
    ! Boucle sur les domaines pour creer calc
    !-----------------------------------------------------------------------
      DO idom=1,nb_dom
    ! Initialisations
          calc%msh%ndim  =cb%msh%ndim
          calc%msh%nb_dom=nb_dom
          nullify(calc%msh%list_vertex)
          nullify(calc%msh%list_face)
          nullify(calc%msh%list_cell)
          nullify(calc%msh%list_send)
          nullify(calc%msh%list_recv)
          CALL empty_mesh(calc%msh,calc%msh%nb_dom,idom)
          IF (associated(calc%list_bloc)) THEN
              deallocate(calc%list_bloc)
              nullify(calc%list_bloc)
          END IF
          nullify(calc%list_bloc)

    ! Determination de la taille des tableaux du domaine
          calc%msh%nb_cell   = 0
          calc%msh%nb_vertex = 0
          calc%msh%nb_face   = 0
          calc%nb_bloc       = 0

          DO ib=1,nb_bloc
              IF (cb%list_bloc(ib)%idom.eq.idom) THEN
                  calc%nb_bloc=calc%nb_bloc+1
                  cb%list_bloc(ib)%first_cell   = calc%msh%nb_cell
                  calc%msh%nb_cell    = calc%msh%nb_cell + cb%list_bloc(ib)%nb_cell
                  cb%list_bloc(ib)%first_face   = calc%msh%nb_face
                  calc%msh%nb_face    = calc%msh%nb_face + cb%list_bloc(ib)%nb_face
                  cb%list_bloc(ib)%first_vertex = calc%msh%nb_vertex
                  calc%msh%nb_vertex  = calc%msh%nb_vertex + cb%list_bloc(ib)%nb_vertex
              END IF
          END DO

          allocate(calc%list_bloc(calc%nb_bloc))
          allocate(calc%msh%list_cell(calc%msh%nb_cell))
          allocate(calc%msh%list_vertex(calc%msh%nb_vertex))
          allocate(calc%msh%list_face(calc%msh%nb_face))
          allocate(calc%msh%list_send(calc%msh%nb_dom))
          allocate(calc%msh%list_recv(calc%msh%nb_dom))

          DO i=1,idom-1
              k=idom-i+(i-1)*(2*nb_dom-i)/2
              IF (nbl(k).ne.0) allocate(calc%msh%list_recv(i)%L(interdom(k)%L(4*nbl(k)-1)))
              IF (nbl(k).ne.0) allocate(calc%msh%list_send(i)%L(interdom(k)%L(4*nbl(k)-0)))
          END DO
          nullify(calc%msh%list_send(idom)%L)
          nullify(calc%msh%list_recv(idom)%L)
          DO i=idom+1,nb_dom
              k=i-idom+(idom-1)*(2*nb_dom-idom)/2
              IF (nbl(k).ne.0) allocate(calc%msh%list_recv(i)%L(interdom(k)%L(4*nbl(k)-0)))
              IF (nbl(k).ne.0) allocate(calc%msh%list_send(i)%L(interdom(k)%L(4*nbl(k)-1)))
          END DO
          CALL initialise(calc%msh)

    ! Maillage des blocs du domaine
          icb=0
          DO ib=1,nb_bloc
              IF (cb%list_bloc(ib)%idom.eq.idom) THEN
                  icb=icb+1
                  calc%list_bloc(icb)=cb%list_bloc(ib)
                  DO i=1,6
                      blv(i)=bl_null
                      ibv=cb%list_bloc(ib)%cl(i)
                      if(ibv.lt.-100000)ibv=-ibv-100000
                      IF (ibv.gt.0) blv(i)=cb%list_bloc(ibv)
                  END DO

                  SELECT CASE (cb%msh%ndim)
                     CASE (1)
                       CALL mesh1d(cb%list_bloc(ib),calc%msh,ib,nb_dom,blv,interdom)

                     CASE (2)
                       CALL MESH2D(cb%list_bloc(ib),calc%msh,ib,nb_dom,blv,interdom)

                     CASE (3)
                       CALL MESH3D(cb%list_bloc(ib),calc%msh,ib,nb_dom,blv,interdom)

                END SELECT
              END IF
          END DO ! Fin boucle sur les blocs

          CALL init_mesh(calc%msh)

    ! On elimine les sommets surabondants
          nvi=calc%msh%nb_vertex
          CALL trivertex(calc%msh)
          write( *,1111) idom,nvi-calc%msh%nb_vertex,nvi,calc%msh%nb_vertex
    1111  format('Domain',i3,' elimination of ',i7,' vertexes out of',i7,', i.e.',i7,' vertexes.')
          write( *,2222) idom,calc%msh%nb_vertex
          write( *,3333) idom,calc%msh%nb_face
          write( *,4444) idom,calc%msh%nb_cell
    2222  format('Domain',i3,' Number of vertexes :',i7)
    3333  format('Domain',i3,' Number of faces    :',i7)
    4444  format('Domain',i3,' Number of cells    :',i7)

    ! Recuperation des donnÃ©es de cb
          calc%tmin           = cb%tmin
          calc%tmax           = cb%tmax
          calc%dt             = cb%dt
          calc%pdtp           = cb%pdtp
          calc%cfl            = cb%cfl
          calc%vmax           = cb%vmax
          calc%isauv          = cb%isauv
          calc%istart         = cb%istart
          calc%ipas           = cb%ipas
          calc%iordre_e       = cb%iordre_e
          calc%iordre_t       = cb%iordre_t
          calc%imeth_int_temps= cb%imeth_int_temps
          calc%point_sondes   = cb%point_sondes
          calc%nb_sondes      = cb%nb_sondes
          calc%num_point_min  = cb%num_point_min
          calc%ityp_flux      = cb%ityp_flux
          calc%ilimiteur      = cb%ilimiteur
          calc%nbds           = cb%nbds
          calc%nrma           = cb%nrma
          calc%vcde           = cb%vcde
          calc%vcra           = cb%vcra
          if(cb%nb_sol.ne.0)call copy_sol(cb,calc)

    ! Initialisation sur tout le domaine

          calc%msh%numdom=idom
          call ini_w(calc,nzonei,zonei)
           !call display(calc%msh)  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
          CALL save_calc(calc)
          CALL empty_mesh(calc%msh,0,0)  ! liberation de la memoire

      END DO ! Fin boucle sur les domaines

      !deallocate(bl)

    !===========================   FIN DE LA ROUTINE    ====================

    END SUBROUTINE block2mesh