Commit 938d9dad authored by Roy Fabrice's avatar Roy Fabrice
Browse files

add selection when halo larger than cell cube

parent 19dc3046
......@@ -66,7 +66,6 @@ contains
#endif
#ifdef DEBUG
......@@ -170,19 +169,23 @@ contains
use char_utils_m
integer, dimension(*), intent(inout) :: cells_cubes_list
integer, intent(in) :: cells_cubes_number
integer, dimension(0:cells_cubes_number), intent(inout) :: cells_cubes_list
real(kind=8), dimension(3), intent(in) :: position
real(kind=8), intent(in) :: radius
integer, intent(in) :: cells_cubes_number
real(kind=8), dimension(2,3,*), intent(in) :: cells_cubes_boundaries
real(kind=8), dimension(2,3,cells_cubes_number), intent(in) :: cells_cubes_boundaries
integer :: alloc_stat
integer :: cells_cubes_dim
real(kind=8) :: cells_cubes_size
integer(kind=4), dimension(3) :: coords
integer(kind=4), dimension(3) :: coords_center
integer(kind=4), allocatable, dimension(:,:) :: coords_loop
integer(kind=4), dimension(3) :: coords_minus
integer(kind=4), dimension(3) :: coords_plus
integer, allocatable, dimension(:,:,:), save :: coords_to_cubes_id
character(len=ERR_MSG_LEN) :: error_message
integer :: icube
integer :: icube, jcube, kcube
integer :: idim, jdim, kdim
real(kind=8) :: sqrt3m1
......@@ -204,6 +207,7 @@ contains
coords(:) = int( 0.5d0*(cells_cubes_boundaries(2,:,icube)+cells_cubes_boundaries(1,:,icube)) / cells_cubes_size, kind=4) + 1
coords_to_cubes_id(coords(1),coords(2),coords(3)) = icube
end do
do kdim = 1, cells_cubes_dim
do jdim = 1, cells_cubes_dim
do idim = 1, cells_cubes_dim
......@@ -218,41 +222,89 @@ contains
end if
! check center
coords(:) = int( position(:) / cells_cubes_size, kind=4) + 1
cells_cubes_list(coords_to_cubes_id(coords(1),coords(2),coords(3))) = 1
! right and left
coords(:) = int( modulo(position(:)+[radius,0.0d0,0.d0],1.0d0) / cells_cubes_size, kind=4) + 1
cells_cubes_list(coords_to_cubes_id(coords(1),coords(2),coords(3))) = 1
coords(:) = int( modulo(position(:)-[radius-1.0d0,0.0d0,0.d0],1.0d0) / cells_cubes_size, kind=4) + 1
cells_cubes_list(coords_to_cubes_id(coords(1),coords(2),coords(3))) = 1
! front and back
coords(:) = int( modulo(position(:)+[0.0d0,radius,0.d0],1.0d0) / cells_cubes_size, kind=4) + 1
cells_cubes_list(coords_to_cubes_id(coords(1),coords(2),coords(3))) = 1
coords(:) = int( modulo(position(:)-[0.0d0,radius-1.0d0,0.d0],1.0d0) / cells_cubes_size, kind=4) + 1
cells_cubes_list(coords_to_cubes_id(coords(1),coords(2),coords(3))) = 1
! up and down
coords(:) = int( modulo(position(:)+[0.0d0,0.d0,radius],1.0d0) / cells_cubes_size, kind=4) + 1
cells_cubes_list(coords_to_cubes_id(coords(1),coords(2),coords(3))) = 1
coords(:) = int( modulo(position(:)-[0.0d0,0.d0,radius-1.0d0],1.0d0) / cells_cubes_size, kind=4) + 1
cells_cubes_list(coords_to_cubes_id(coords(1),coords(2),coords(3))) = 1
! corners
sqrt3m1 = radius/sqrt(3.d0)
coords(:) = int( modulo(position(:)+[sqrt3m1,sqrt3m1,sqrt3m1],1.d0) / cells_cubes_size, kind=4) + 1
cells_cubes_list(coords_to_cubes_id(coords(1),coords(2),coords(3))) = 1
coords(:) = int( modulo(position(:)+[sqrt3m1,sqrt3m1,1.d0-sqrt3m1],1.d0) / cells_cubes_size, kind=4) + 1
cells_cubes_list(coords_to_cubes_id(coords(1),coords(2),coords(3))) = 1
coords(:) = int( modulo(position(:)+[sqrt3m1,1.d0-sqrt3m1,sqrt3m1],1.d0) / cells_cubes_size, kind=4) + 1
cells_cubes_list(coords_to_cubes_id(coords(1),coords(2),coords(3))) = 1
coords(:) = int( modulo(position(:)+[sqrt3m1,1.d0-sqrt3m1,1.d0-sqrt3m1],1.d0) / cells_cubes_size, kind=4) + 1
cells_cubes_list(coords_to_cubes_id(coords(1),coords(2),coords(3))) = 1
coords(:) = int( modulo(position(:)+[1.d0-sqrt3m1,sqrt3m1,sqrt3m1],1.d0) / cells_cubes_size, kind=4) + 1
cells_cubes_list(coords_to_cubes_id(coords(1),coords(2),coords(3))) = 1
coords(:) = int( modulo(position(:)+[1.d0-sqrt3m1,sqrt3m1,1.d0-sqrt3m1],1.d0) / cells_cubes_size, kind=4) + 1
cells_cubes_list(coords_to_cubes_id(coords(1),coords(2),coords(3))) = 1
coords(:) = int( modulo(position(:)+[1.d0-sqrt3m1,1.d0-sqrt3m1,sqrt3m1],1.d0) / cells_cubes_size, kind=4) + 1
cells_cubes_list(coords_to_cubes_id(coords(1),coords(2),coords(3))) = 1
coords(:) = int( modulo(position(:)+[1.d0-sqrt3m1,1.d0-sqrt3m1,1.d0-sqrt3m1],1.d0) / cells_cubes_size, kind=4) + 1
coords_center(:) = int( position(:) / cells_cubes_size, kind=4) + 1
cells_cubes_list(coords_to_cubes_id(coords(1),coords(2),coords(3))) = 1
if(radius < cells_cubes_size) then
! right and left
coords(:) = int( modulo(position(:)+[radius,0.0d0,0.d0],1.0d0) / cells_cubes_size, kind=4) + 1
cells_cubes_list(coords_to_cubes_id(coords(1),coords(2),coords(3))) = 1
coords(:) = int( modulo(position(:)-[radius-1.0d0,0.0d0,0.d0],1.0d0) / cells_cubes_size, kind=4) + 1
cells_cubes_list(coords_to_cubes_id(coords(1),coords(2),coords(3))) = 1
! front and back
coords(:) = int( modulo(position(:)+[0.0d0,radius,0.d0],1.0d0) / cells_cubes_size, kind=4) + 1
cells_cubes_list(coords_to_cubes_id(coords(1),coords(2),coords(3))) = 1
coords(:) = int( modulo(position(:)-[0.0d0,radius-1.0d0,0.d0],1.0d0) / cells_cubes_size, kind=4) + 1
cells_cubes_list(coords_to_cubes_id(coords(1),coords(2),coords(3))) = 1
! up and down
coords(:) = int( modulo(position(:)+[0.0d0,0.d0,radius],1.0d0) / cells_cubes_size, kind=4) + 1
cells_cubes_list(coords_to_cubes_id(coords(1),coords(2),coords(3))) = 1
coords(:) = int( modulo(position(:)-[0.0d0,0.d0,radius-1.0d0],1.0d0) / cells_cubes_size, kind=4) + 1
cells_cubes_list(coords_to_cubes_id(coords(1),coords(2),coords(3))) = 1
! corners
sqrt3m1 = radius/sqrt(3.d0)
coords(:) = int( modulo(position(:)+[sqrt3m1,sqrt3m1,sqrt3m1],1.d0) / cells_cubes_size, kind=4) + 1
cells_cubes_list(coords_to_cubes_id(coords(1),coords(2),coords(3))) = 1
coords(:) = int( modulo(position(:)+[sqrt3m1,sqrt3m1,1.d0-sqrt3m1],1.d0) / cells_cubes_size, kind=4) + 1
cells_cubes_list(coords_to_cubes_id(coords(1),coords(2),coords(3))) = 1
coords(:) = int( modulo(position(:)+[sqrt3m1,1.d0-sqrt3m1,sqrt3m1],1.d0) / cells_cubes_size, kind=4) + 1
cells_cubes_list(coords_to_cubes_id(coords(1),coords(2),coords(3))) = 1
coords(:) = int( modulo(position(:)+[sqrt3m1,1.d0-sqrt3m1,1.d0-sqrt3m1],1.d0) / cells_cubes_size, kind=4) + 1
cells_cubes_list(coords_to_cubes_id(coords(1),coords(2),coords(3))) = 1
coords(:) = int( modulo(position(:)+[1.d0-sqrt3m1,sqrt3m1,sqrt3m1],1.d0) / cells_cubes_size, kind=4) + 1
cells_cubes_list(coords_to_cubes_id(coords(1),coords(2),coords(3))) = 1
coords(:) = int( modulo(position(:)+[1.d0-sqrt3m1,sqrt3m1,1.d0-sqrt3m1],1.d0) / cells_cubes_size, kind=4) + 1
cells_cubes_list(coords_to_cubes_id(coords(1),coords(2),coords(3))) = 1
coords(:) = int( modulo(position(:)+[1.d0-sqrt3m1,1.d0-sqrt3m1,sqrt3m1],1.d0) / cells_cubes_size, kind=4) + 1
cells_cubes_list(coords_to_cubes_id(coords(1),coords(2),coords(3))) = 1
coords(:) = int( modulo(position(:)+[1.d0-sqrt3m1,1.d0-sqrt3m1,1.d0-sqrt3m1],1.d0) / cells_cubes_size, kind=4) + 1
cells_cubes_list(coords_to_cubes_id(coords(1),coords(2),coords(3))) = 1
else
! if the radius of the halo is larger than the size of a cube of cells, we select every cubes from "center - radius"
! to "center + radius" in every direction
coords_minus(:) = int( modulo(position(:)-[radius-1.0d0,radius-1.0d0,radius-1.0d0],1.0d0) / cells_cubes_size, kind=4) + 1
coords_plus(:) = int( modulo(position(:)+[radius,radius,radius],1.0d0) / cells_cubes_size, kind=4) + 1
! coords_loop = loop indices from 1 to cells_cube_dim in each direction
! needed to loop over cubes only for selected cubes
! coords_loop(3,1:cells_cubes_dim) = 1 if cube is selected, = 0 if not selected
allocate(coords_loop(3,cells_cubes_dim),stat=alloc_stat, errmsg = error_message)
if(alloc_stat /= 0) then
call Allocate_error('coords_loop','Select_cells_cubes', error_message, alloc_stat, mpi_process%rank)
end if
coords_loop = 0
do idim = 1, 3
! the halo is in the middle of the simulation domain
if(coords_minus(idim) <= coords_plus(idim)) then
do icube = coords_minus(idim), coords_plus(idim)
coords_loop(idim,icube) = 1
end do
! the halo is on an edge (perdiodic conditions)
else
do icube = coords_minus(idim), cells_cubes_dim
coords_loop(idim,icube) = 1
end do
do icube = 1, coords_plus(idim)
coords_loop(idim,icube) = 1
end do
end if
end do
! loop over all cubes
! if the indices correspond to a selected cube then the cube is added to the list
do icube = 1, cells_cubes_dim
do jcube = 1, cells_cubes_dim
do kcube = 1, cells_cubes_dim
if(coords_loop(1,icube) == 1 .and. coords_loop(2,jcube)==1 .and. coords_loop(3,kcube)==1) then
cells_cubes_list(coords_to_cubes_id(icube,jcube,kcube)) = 1
end if
end do
end do
end do
end if
! count number of cubes to read
do icube = 1, cells_cubes_number
if(cells_cubes_list(icube)/= 0) cells_cubes_list(0) = cells_cubes_list(0)+1
end do
#ifdef DEBUG
write(ERROR_UNIT,'(a,i0)') 'Select_cells_cubes ends on process', mpi_process%rank
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment