TABLE OF CONTENTS
UTI/PRE/block2mesh [ 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