223031012
/
DAMASK_EICMD
mirror of https://github.com/achalhp/DAMASK_EICMD.git
1
0
Fork 0
Code Issues Projects Releases Wiki Activity

cleaning

This commit is contained in:
Martin Diehl 2019-05-31 20:48:16 +02:00
parent 95d93ce58c
commit c045bd6355
1 changed files with 309 additions and 340 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

283
src/mesh_grid.f90
View File

@ -10,10 +10,19 @@ module mesh
use prec use prec
use geometry_plastic_nonlocal use geometry_plastic_nonlocal
use mesh_base use mesh_base
#include <petsc/finclude/petscsys.h>
use PETScsys
use DAMASK_interface
use IO
use debug
use numerics
use FEsolving
implicit none implicit none
private private
include 'fftw3-mpi.f03'
integer, public, protected :: & integer, public, protected :: &
mesh_Nnodes mesh_Nnodes
@ -51,7 +60,7 @@ module mesh
logical, dimension(3), public, parameter :: mesh_periodicSurface = .true. !< flag indicating periodic outer surfaces (used for fluxes) logical, dimension(3), public, parameter :: mesh_periodicSurface = .true. !< flag indicating periodic outer surfaces (used for fluxes)
! grid specific ! grid specific
integer, dimension(3), public, protected :: & integer, dimension(3), public, protected :: &
grid !< (global) grid grid !< (global) grid
integer, public, protected :: & integer, public, protected :: &
@ -103,7 +112,7 @@ subroutine tMesh_grid_init(self,nodes)
class(tMesh_grid) :: self class(tMesh_grid) :: self
real(pReal), dimension(:,:), intent(in) :: nodes real(pReal), dimension(:,:), intent(in) :: nodes
call self%tMesh%init('grid',10_pInt,nodes) call self%tMesh%init('grid',10,nodes)
end subroutine tMesh_grid_init end subroutine tMesh_grid_init
@ -113,25 +122,6 @@ end subroutine tMesh_grid_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine mesh_init(ip,el) subroutine mesh_init(ip,el)
#include <petsc/finclude/petscsys.h>
use PETScsys
use DAMASK_interface
use IO, only: &
IO_error
use debug, only: &
debug_e, &
debug_i, &
debug_level, &
debug_mesh, &
debug_levelBasic
use numerics, only: &
numerics_unitlength
use FEsolving, only: &
FEsolving_execElem, &
FEsolving_execIP
include 'fftw3-mpi.f03'
integer(C_INTPTR_T) :: devNull, local_K, local_K_offset integer(C_INTPTR_T) :: devNull, local_K, local_K_offset
integer :: ierr, worldsize, j integer :: ierr, worldsize, j
integer, intent(in), optional :: el, ip integer, intent(in), optional :: el, ip
@ -141,15 +131,15 @@ subroutine mesh_init(ip,el)
mesh_unitlength = numerics_unitlength ! set physical extent of a length unit in mesh mesh_unitlength = numerics_unitlength ! set physical extent of a length unit in mesh
myDebug = (iand(debug_level(debug_mesh),debug_levelBasic) /= 0_pInt) myDebug = (iand(debug_level(debug_mesh),debug_levelBasic) /= 0)
call fftw_mpi_init() call fftw_mpi_init()
call mesh_spectral_read_grid() call mesh_spectral_read_grid()
call MPI_comm_size(PETSC_COMM_WORLD, worldsize, ierr) call MPI_comm_size(PETSC_COMM_WORLD, worldsize, ierr)
if(ierr /=0_pInt) call IO_error(894_pInt, ext_msg='MPI_comm_size') if(ierr /=0) call IO_error(894, ext_msg='MPI_comm_size')
if(worldsize>grid(3)) call IO_error(894_pInt, ext_msg='number of processes exceeds grid(3)') if(worldsize>grid(3)) call IO_error(894, ext_msg='number of processes exceeds grid(3)')
devNull = fftw_mpi_local_size_3d(int(grid(3),C_INTPTR_T), & devNull = fftw_mpi_local_size_3d(int(grid(3),C_INTPTR_T), &
@ -165,7 +155,7 @@ subroutine mesh_init(ip,el)
mesh_NcpElemsGlobal = product(grid) mesh_NcpElemsGlobal = product(grid)
mesh_Nnodes = product(grid(1:2) + 1_pInt)*(grid3 + 1_pInt) mesh_Nnodes = product(grid(1:2) + 1)*(grid3 + 1)
mesh_node0 = mesh_spectral_build_nodes() mesh_node0 = mesh_spectral_build_nodes()
mesh_node = mesh_node0 mesh_node = mesh_node0
@ -196,13 +186,13 @@ subroutine mesh_init(ip,el)
if (myDebug) write(6,'(a)') ' Built IP neighborhood'; flush(6) if (myDebug) write(6,'(a)') ' Built IP neighborhood'; flush(6)
if (debug_e < 1 .or. debug_e > theMesh%nElems) & if (debug_e < 1 .or. debug_e > theMesh%nElems) &
call IO_error(602_pInt,ext_msg='element') ! selected element does not exist call IO_error(602,ext_msg='element') ! selected element does not exist
if (debug_i < 1 .or. debug_i > theMesh%elem%nIPs) & if (debug_i < 1 .or. debug_i > theMesh%elem%nIPs) &
call IO_error(602_pInt,ext_msg='IP') ! selected element does not have requested IP call IO_error(602,ext_msg='IP') ! selected element does not have requested IP
FEsolving_execElem = [ 1_pInt,theMesh%nElems ] ! parallel loop bounds set to comprise all DAMASK elements FEsolving_execElem = [ 1,theMesh%nElems ] ! parallel loop bounds set to comprise all DAMASK elements
allocate(FEsolving_execIP(2_pInt,theMesh%nElems), source=1_pInt) ! parallel loop bounds set to comprise from first IP... allocate(FEsolving_execIP(2,theMesh%nElems), source=1) ! parallel loop bounds set to comprise from first IP...
forall (j = 1_pInt:theMesh%nElems) FEsolving_execIP(2,j) = theMesh%elem%nIPs ! ...up to own IP count for each element forall (j = 1:theMesh%nElems) FEsolving_execIP(2,j) = theMesh%elem%nIPs ! ...up to own IP count for each element
!!!! COMPATIBILITY HACK !!!! !!!! COMPATIBILITY HACK !!!!
@ -219,22 +209,13 @@ end subroutine mesh_init
! supposed to be called only once! ! supposed to be called only once!
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine mesh_spectral_read_grid() subroutine mesh_spectral_read_grid()
use IO, only: &
IO_stringPos, &
IO_lc, &
IO_stringValue, &
IO_intValue, &
IO_floatValue, &
IO_error
use DAMASK_interface, only: &
geometryFile
character(len=:), allocatable :: rawData character(len=:), allocatable :: rawData
character(len=65536) :: line character(len=65536) :: line
integer, allocatable, dimension(:) :: chunkPos integer, allocatable, dimension(:) :: chunkPos
integer :: h =- 1_pInt integer :: h =- 1
integer :: & integer :: &
headerLength = -1_pInt, & !< length of header (in lines) headerLength = -1, & !< length of header (in lines)
fileLength, & !< length of the geom file (in characters) fileLength, & !< length of the geom file (in characters)
fileUnit, & fileUnit, &
startPos, endPos, & startPos, endPos, &
@ -245,7 +226,7 @@ subroutine mesh_spectral_read_grid()
e, & !< "element", i.e. spectral collocation point e, & !< "element", i.e. spectral collocation point
i, j i, j
grid = -1_pInt grid = -1
geomSize = -1.0_pReal geomSize = -1.0_pReal
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -253,7 +234,7 @@ subroutine mesh_spectral_read_grid()
inquire(file = trim(geometryFile), size=fileLength) inquire(file = trim(geometryFile), size=fileLength)
open(newunit=fileUnit, file=trim(geometryFile), access='stream',& open(newunit=fileUnit, file=trim(geometryFile), access='stream',&
status='old', position='rewind', action='read',iostat=myStat) status='old', position='rewind', action='read',iostat=myStat)
if(myStat /= 0_pInt) call IO_error(100_pInt,ext_msg=trim(geometryFile)) if(myStat /= 0) call IO_error(100,ext_msg=trim(geometryFile))
allocate(character(len=fileLength)::rawData) allocate(character(len=fileLength)::rawData)
read(fileUnit) rawData read(fileUnit) rawData
close(fileUnit) close(fileUnit)
@ -263,106 +244,106 @@ subroutine mesh_spectral_read_grid()
endPos = index(rawData,new_line('')) endPos = index(rawData,new_line(''))
if(endPos <= index(rawData,'head')) then if(endPos <= index(rawData,'head')) then
startPos = len(rawData) startPos = len(rawData)
call IO_error(error_ID=841_pInt, ext_msg='mesh_spectral_read_grid') call IO_error(error_ID=841, ext_msg='mesh_spectral_read_grid')
else else
chunkPos = IO_stringPos(rawData(1:endPos)) chunkPos = IO_stringPos(rawData(1:endPos))
if (chunkPos(1) < 2_pInt) call IO_error(error_ID=841_pInt, ext_msg='mesh_spectral_read_grid') if (chunkPos(1) < 2) call IO_error(error_ID=841, ext_msg='mesh_spectral_read_grid')
headerLength = IO_intValue(rawData(1:endPos),chunkPos,1_pInt) headerLength = IO_intValue(rawData(1:endPos),chunkPos,1)
startPos = endPos + 1_pInt startPos = endPos + 1
endif endif
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! read and interprete header ! read and interprete header
l = 0 l = 0
do while (l < headerLength .and. startPos < len(rawData)) do while (l < headerLength .and. startPos < len(rawData))
endPos = startPos + index(rawData(startPos:),new_line('')) - 1_pInt endPos = startPos + index(rawData(startPos:),new_line('')) - 1
if (endPos < startPos) endPos = len(rawData) ! end of file without new line if (endPos < startPos) endPos = len(rawData) ! end of file without new line
line = rawData(startPos:endPos) line = rawData(startPos:endPos)
startPos = endPos + 1_pInt startPos = endPos + 1
l = l + 1_pInt l = l + 1
chunkPos = IO_stringPos(trim(line)) chunkPos = IO_stringPos(trim(line))
if (chunkPos(1) < 2) cycle ! need at least one keyword value pair if (chunkPos(1) < 2) cycle ! need at least one keyword value pair
select case ( IO_lc(IO_StringValue(trim(line),chunkPos,1_pInt,.true.)) ) select case ( IO_lc(IO_StringValue(trim(line),chunkPos,1,.true.)) )
case ('grid') case ('grid')
if (chunkPos(1) > 6) then if (chunkPos(1) > 6) then
do j = 2_pInt,6_pInt,2_pInt do j = 2,6,2
select case (IO_lc(IO_stringValue(line,chunkPos,j))) select case (IO_lc(IO_stringValue(line,chunkPos,j)))
case('a') case('a')
grid(1) = IO_intValue(line,chunkPos,j+1_pInt) grid(1) = IO_intValue(line,chunkPos,j+1)
case('b') case('b')
grid(2) = IO_intValue(line,chunkPos,j+1_pInt) grid(2) = IO_intValue(line,chunkPos,j+1)
case('c') case('c')
grid(3) = IO_intValue(line,chunkPos,j+1_pInt) grid(3) = IO_intValue(line,chunkPos,j+1)
end select end select
enddo enddo
endif endif
case ('size') case ('size')
if (chunkPos(1) > 6) then if (chunkPos(1) > 6) then
do j = 2_pInt,6_pInt,2_pInt do j = 2,6,2
select case (IO_lc(IO_stringValue(line,chunkPos,j))) select case (IO_lc(IO_stringValue(line,chunkPos,j)))
case('x') case('x')
geomSize(1) = IO_floatValue(line,chunkPos,j+1_pInt) geomSize(1) = IO_floatValue(line,chunkPos,j+1)
case('y') case('y')
geomSize(2) = IO_floatValue(line,chunkPos,j+1_pInt) geomSize(2) = IO_floatValue(line,chunkPos,j+1)
case('z') case('z')
geomSize(3) = IO_floatValue(line,chunkPos,j+1_pInt) geomSize(3) = IO_floatValue(line,chunkPos,j+1)
end select end select
enddo enddo
endif endif
case ('homogenization') case ('homogenization')
if (chunkPos(1) > 1) h = IO_intValue(line,chunkPos,2_pInt) if (chunkPos(1) > 1) h = IO_intValue(line,chunkPos,2)
end select end select
enddo enddo
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! sanity checks ! sanity checks
if(h < 1_pInt) & if(h < 1) &
call IO_error(error_ID = 842_pInt, ext_msg='homogenization (mesh_spectral_read_grid)') call IO_error(error_ID = 842, ext_msg='homogenization (mesh_spectral_read_grid)')
if(any(grid < 1_pInt)) & if(any(grid < 1)) &
call IO_error(error_ID = 842_pInt, ext_msg='grid (mesh_spectral_read_grid)') call IO_error(error_ID = 842, ext_msg='grid (mesh_spectral_read_grid)')
if(any(geomSize < 0.0_pReal)) & if(any(geomSize < 0.0_pReal)) &
call IO_error(error_ID = 842_pInt, ext_msg='size (mesh_spectral_read_grid)') call IO_error(error_ID = 842, ext_msg='size (mesh_spectral_read_grid)')
allocate(microGlobal(product(grid)), source = -1_pInt) allocate(microGlobal(product(grid)), source = -1)
allocate(mesh_homogenizationAt(product(grid)), source = h) ! too large in case of MPI (shrink later, not very elegant) allocate(mesh_homogenizationAt(product(grid)), source = h) ! too large in case of MPI (shrink later, not very elegant)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! read and interpret content ! read and interpret content
e = 1_pInt e = 1
do while (startPos < len(rawData)) do while (startPos < len(rawData))
endPos = startPos + index(rawData(startPos:),new_line('')) - 1_pInt endPos = startPos + index(rawData(startPos:),new_line('')) - 1
if (endPos < startPos) endPos = len(rawData) ! end of file without new line if (endPos < startPos) endPos = len(rawData) ! end of file without new line
line = rawData(startPos:endPos) line = rawData(startPos:endPos)
startPos = endPos + 1_pInt startPos = endPos + 1
l = l + 1_pInt l = l + 1
chunkPos = IO_stringPos(trim(line)) chunkPos = IO_stringPos(trim(line))
noCompression: if (chunkPos(1) /= 3) then noCompression: if (chunkPos(1) /= 3) then
c = chunkPos(1) c = chunkPos(1)
microGlobal(e:e+c-1_pInt) = [(IO_intValue(line,chunkPos,i+1_pInt), i=0_pInt, c-1_pInt)] microGlobal(e:e+c-1) = [(IO_intValue(line,chunkPos,i+1), i=0, c-1)]
else noCompression else noCompression
compression: if (IO_lc(IO_stringValue(line,chunkPos,2)) == 'of') then compression: if (IO_lc(IO_stringValue(line,chunkPos,2)) == 'of') then
c = IO_intValue(line,chunkPos,1) c = IO_intValue(line,chunkPos,1)
microGlobal(e:e+c-1_pInt) = [(IO_intValue(line,chunkPos,3),i = 1_pInt,IO_intValue(line,chunkPos,1))] microGlobal(e:e+c-1) = [(IO_intValue(line,chunkPos,3),i = 1,IO_intValue(line,chunkPos,1))]
else if (IO_lc(IO_stringValue(line,chunkPos,2)) == 'to') then compression else if (IO_lc(IO_stringValue(line,chunkPos,2)) == 'to') then compression
c = abs(IO_intValue(line,chunkPos,3) - IO_intValue(line,chunkPos,1)) + 1_pInt c = abs(IO_intValue(line,chunkPos,3) - IO_intValue(line,chunkPos,1)) + 1
o = merge(+1_pInt, -1_pInt, IO_intValue(line,chunkPos,3) > IO_intValue(line,chunkPos,1)) o = merge(+1, -1, IO_intValue(line,chunkPos,3) > IO_intValue(line,chunkPos,1))
microGlobal(e:e+c-1_pInt) = [(i, i = IO_intValue(line,chunkPos,1),IO_intValue(line,chunkPos,3),o)] microGlobal(e:e+c-1) = [(i, i = IO_intValue(line,chunkPos,1),IO_intValue(line,chunkPos,3),o)]
else compression else compression
c = chunkPos(1) c = chunkPos(1)
microGlobal(e:e+c-1_pInt) = [(IO_intValue(line,chunkPos,i+1_pInt), i=0_pInt, c-1_pInt)] microGlobal(e:e+c-1) = [(IO_intValue(line,chunkPos,i+1), i=0, c-1)]
endif compression endif compression
endif noCompression endif noCompression
e = e+c e = e+c
end do end do
if (e-1 /= product(grid)) call IO_error(error_ID = 843_pInt, el=e) if (e-1 /= product(grid)) call IO_error(error_ID = 843, el=e)
end subroutine mesh_spectral_read_grid end subroutine mesh_spectral_read_grid
@ -413,29 +394,29 @@ end function mesh_build_ipCoordinates
!> @brief Store FEid, type, material, texture, and node list per element. !> @brief Store FEid, type, material, texture, and node list per element.
!! Allocates global array 'mesh_element' !! Allocates global array 'mesh_element'
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine mesh_spectral_build_elements() subroutine mesh_spectral_build_elements
integer :: & integer :: &
e, & e, &
elemOffset elemOffset
allocate(mesh_element (4+8,theMesh%nElems), source = 0)
allocate(mesh_element (4_pInt+8_pInt,theMesh%nElems), source = 0_pInt)
elemOffset = product(grid(1:2))*grid3Offset elemOffset = product(grid(1:2))*grid3Offset
do e=1, theMesh%nElems do e=1, theMesh%nElems
mesh_element( 1,e) = -1_pInt ! DEPRECATED mesh_element( 1,e) = -1 ! DEPRECATED
mesh_element( 2,e) = -1_pInt ! DEPRECATED mesh_element( 2,e) = -1 ! DEPRECATED
mesh_element( 3,e) = mesh_homogenizationAt(e) mesh_element( 3,e) = mesh_homogenizationAt(e)
mesh_element( 4,e) = microGlobal(e+elemOffset) ! microstructure mesh_element( 4,e) = microGlobal(e+elemOffset) ! microstructure
mesh_element( 5,e) = e + (e-1_pInt)/grid(1) + & mesh_element( 5,e) = e + (e-1)/grid(1) + &
((e-1_pInt)/(grid(1)*grid(2)))*(grid(1)+1_pInt) ! base node ((e-1)/(grid(1)*grid(2)))*(grid(1)+1) ! base node
mesh_element( 6,e) = mesh_element(5,e) + 1_pInt mesh_element( 6,e) = mesh_element(5,e) + 1
mesh_element( 7,e) = mesh_element(5,e) + grid(1) + 2_pInt mesh_element( 7,e) = mesh_element(5,e) + grid(1) + 2
mesh_element( 8,e) = mesh_element(5,e) + grid(1) + 1_pInt mesh_element( 8,e) = mesh_element(5,e) + grid(1) + 1
mesh_element( 9,e) = mesh_element(5,e) +(grid(1) + 1_pInt) * (grid(2) + 1_pInt) ! second floor base node mesh_element( 9,e) = mesh_element(5,e) +(grid(1) + 1) * (grid(2) + 1) ! second floor base node
mesh_element(10,e) = mesh_element(9,e) + 1_pInt mesh_element(10,e) = mesh_element(9,e) + 1
mesh_element(11,e) = mesh_element(9,e) + grid(1) + 2_pInt mesh_element(11,e) = mesh_element(9,e) + grid(1) + 2
mesh_element(12,e) = mesh_element(9,e) + grid(1) + 1_pInt mesh_element(12,e) = mesh_element(9,e) + grid(1) + 1
enddo enddo
end subroutine mesh_spectral_build_elements end subroutine mesh_spectral_build_elements
@ -450,44 +431,44 @@ subroutine mesh_spectral_build_ipNeighborhood
integer :: & integer :: &
x,y,z, & x,y,z, &
e e
allocate(mesh_ipNeighborhood(3,6,1,theMesh%nElems),source=0_pInt) allocate(mesh_ipNeighborhood(3,6,1,theMesh%nElems),source=0)
e = 0_pInt e = 0
do z = 0_pInt,grid3-1_pInt do z = 0,grid3-1
do y = 0_pInt,grid(2)-1_pInt do y = 0,grid(2)-1
do x = 0_pInt,grid(1)-1_pInt do x = 0,grid(1)-1
e = e + 1_pInt e = e + 1
mesh_ipNeighborhood(1,1,1,e) = z * grid(1) * grid(2) & mesh_ipNeighborhood(1,1,1,e) = z * grid(1) * grid(2) &
+ y * grid(1) & + y * grid(1) &
+ modulo(x+1_pInt,grid(1)) & + modulo(x+1,grid(1)) &
+ 1_pInt + 1
mesh_ipNeighborhood(1,2,1,e) = z * grid(1) * grid(2) & mesh_ipNeighborhood(1,2,1,e) = z * grid(1) * grid(2) &
+ y * grid(1) & + y * grid(1) &
+ modulo(x-1_pInt,grid(1)) & + modulo(x-1,grid(1)) &
+ 1_pInt + 1
mesh_ipNeighborhood(1,3,1,e) = z * grid(1) * grid(2) & mesh_ipNeighborhood(1,3,1,e) = z * grid(1) * grid(2) &
+ modulo(y+1_pInt,grid(2)) * grid(1) & + modulo(y+1,grid(2)) * grid(1) &
+ x & + x &
+ 1_pInt + 1
mesh_ipNeighborhood(1,4,1,e) = z * grid(1) * grid(2) & mesh_ipNeighborhood(1,4,1,e) = z * grid(1) * grid(2) &
+ modulo(y-1_pInt,grid(2)) * grid(1) & + modulo(y-1,grid(2)) * grid(1) &
+ x & + x &
+ 1_pInt + 1
mesh_ipNeighborhood(1,5,1,e) = modulo(z+1_pInt,grid3) * grid(1) * grid(2) & mesh_ipNeighborhood(1,5,1,e) = modulo(z+1,grid3) * grid(1) * grid(2) &
+ y * grid(1) & + y * grid(1) &
+ x & + x &
+ 1_pInt + 1
mesh_ipNeighborhood(1,6,1,e) = modulo(z-1_pInt,grid3) * grid(1) * grid(2) & mesh_ipNeighborhood(1,6,1,e) = modulo(z-1,grid3) * grid(1) * grid(2) &
+ y * grid(1) & + y * grid(1) &
+ x & + x &
+ 1_pInt + 1
mesh_ipNeighborhood(2,1:6,1,e) = 1_pInt mesh_ipNeighborhood(2,1:6,1,e) = 1
mesh_ipNeighborhood(3,1,1,e) = 2_pInt mesh_ipNeighborhood(3,1,1,e) = 2
mesh_ipNeighborhood(3,2,1,e) = 1_pInt mesh_ipNeighborhood(3,2,1,e) = 1
mesh_ipNeighborhood(3,3,1,e) = 4_pInt mesh_ipNeighborhood(3,3,1,e) = 4
mesh_ipNeighborhood(3,4,1,e) = 3_pInt mesh_ipNeighborhood(3,4,1,e) = 3
mesh_ipNeighborhood(3,5,1,e) = 6_pInt mesh_ipNeighborhood(3,5,1,e) = 6
mesh_ipNeighborhood(3,6,1,e) = 5_pInt mesh_ipNeighborhood(3,6,1,e) = 5
enddo enddo
enddo enddo
enddo enddo
@ -499,10 +480,6 @@ end subroutine mesh_spectral_build_ipNeighborhood
!> @brief builds mesh of (distorted) cubes for given coordinates (= center of the cubes) !> @brief builds mesh of (distorted) cubes for given coordinates (= center of the cubes)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function mesh_nodesAroundCentres(gDim,Favg,centres) result(nodes) function mesh_nodesAroundCentres(gDim,Favg,centres) result(nodes)
use debug, only: &
debug_mesh, &
debug_level, &
debug_levelBasic
real(pReal), intent(in), dimension(:,:,:,:) :: & real(pReal), intent(in), dimension(:,:,:,:) :: &
centres centres
@ -518,22 +495,22 @@ function mesh_nodesAroundCentres(gDim,Favg,centres) result(nodes)
integer :: & integer :: &
i,j,k,n i,j,k,n
integer, dimension(3), parameter :: & integer, dimension(3), parameter :: &
diag = 1_pInt diag = 1
integer, dimension(3) :: & integer, dimension(3) :: &
shift = 0_pInt, & shift = 0, &
lookup = 0_pInt, & lookup = 0, &
me = 0_pInt, & me = 0, &
iRes = 0_pInt iRes = 0
integer, dimension(3,8) :: & integer, dimension(3,8) :: &
neighbor = reshape([ & neighbor = reshape([ &
0_pInt, 0_pInt, 0_pInt, & 0, 0, 0, &
1_pInt, 0_pInt, 0_pInt, & 1, 0, 0, &
1_pInt, 1_pInt, 0_pInt, & 1, 1, 0, &
0_pInt, 1_pInt, 0_pInt, & 0, 1, 0, &
0_pInt, 0_pInt, 1_pInt, & 0, 0, 1, &
1_pInt, 0_pInt, 1_pInt, & 1, 0, 1, &
1_pInt, 1_pInt, 1_pInt, & 1, 1, 1, &
0_pInt, 1_pInt, 1_pInt ], [3,8]) 0, 1, 1 ], [3,8])
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! initializing variables ! initializing variables
@ -541,40 +518,32 @@ function mesh_nodesAroundCentres(gDim,Favg,centres) result(nodes)
nodes = 0.0_pReal nodes = 0.0_pReal
wrappedCentres = 0.0_pReal wrappedCentres = 0.0_pReal
!--------------------------------------------------------------------------------------------------
! report
if (iand(debug_level(debug_mesh),debug_levelBasic) /= 0_pInt) then
write(6,'(a)') ' Meshing cubes around centroids'
write(6,'(a,3(e12.5))') ' Dimension: ', gDim
write(6,'(a,3(i5))') ' Resolution:', iRes
endif
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! building wrappedCentres = centroids + ghosts ! building wrappedCentres = centroids + ghosts
wrappedCentres(1:3,2_pInt:iRes(1)+1_pInt,2_pInt:iRes(2)+1_pInt,2_pInt:iRes(3)+1_pInt) = centres wrappedCentres(1:3,2:iRes(1)+1,2:iRes(2)+1,2:iRes(3)+1) = centres
do k = 0_pInt,iRes(3)+1_pInt do k = 0,iRes(3)+1
do j = 0_pInt,iRes(2)+1_pInt do j = 0,iRes(2)+1
do i = 0_pInt,iRes(1)+1_pInt do i = 0,iRes(1)+1
if (k==0_pInt .or. k==iRes(3)+1_pInt .or. & ! z skin if (k==0 .or. k==iRes(3)+1 .or. & ! z skin
j==0_pInt .or. j==iRes(2)+1_pInt .or. & ! y skin j==0 .or. j==iRes(2)+1 .or. & ! y skin
i==0_pInt .or. i==iRes(1)+1_pInt ) then ! x skin i==0 .or. i==iRes(1)+1 ) then ! x skin
me = [i,j,k] ! me on skin me = [i,j,k] ! me on skin
shift = sign(abs(iRes+diag-2_pInt*me)/(iRes+diag),iRes+diag-2_pInt*me) shift = sign(abs(iRes+diag-2*me)/(iRes+diag),iRes+diag-2*me)
lookup = me-diag+shift*iRes lookup = me-diag+shift*iRes
wrappedCentres(1:3,i+1_pInt, j+1_pInt, k+1_pInt) = & wrappedCentres(1:3,i+1, j+1, k+1) = &
centres(1:3,lookup(1)+1_pInt,lookup(2)+1_pInt,lookup(3)+1_pInt) & centres(1:3,lookup(1)+1,lookup(2)+1,lookup(3)+1) &
- matmul(Favg, real(shift,pReal)*gDim) - matmul(Favg, real(shift,pReal)*gDim)
endif endif
enddo; enddo; enddo enddo; enddo; enddo
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! averaging ! averaging
do k = 0_pInt,iRes(3); do j = 0_pInt,iRes(2); do i = 0_pInt,iRes(1) do k = 0,iRes(3); do j = 0,iRes(2); do i = 0,iRes(1)
do n = 1_pInt,8_pInt do n = 1,8
nodes(1:3,i+1_pInt,j+1_pInt,k+1_pInt) = & nodes(1:3,i+1,j+1,k+1) = &
nodes(1:3,i+1_pInt,j+1_pInt,k+1_pInt) + wrappedCentres(1:3,i+1_pInt+neighbor(1,n), & nodes(1:3,i+1,j+1,k+1) + wrappedCentres(1:3,i+1+neighbor(1,n), &
j+1_pInt+neighbor(2,n), & j+1+neighbor(2,n), &
k+1_pInt+neighbor(3,n) ) k+1+neighbor(3,n) )
enddo enddo
enddo; enddo; enddo enddo; enddo; enddo
nodes = nodes/8.0_pReal nodes = nodes/8.0_pReal