Create a fully triply-periodic structured hexahedral mesh and store it in SELF's unstructured mesh format. Element geometry and ordering are identical to UniformStructuredMesh; the only difference is connectivity. The faces on the six domain boundaries are wired as interior faces whose neighbor is the element on the opposite side of the domain. This realises the triply periodic box T^3 = [0,Lx] x [0,Ly] x [0,Lz] required by, e.g., the Arnold-Beltrami-Childress (ABC) flow benchmark.
Input - this : Fresh/empty mesh3d_t object - nxPerTile,nyPerTile,nzPerTile : Elements per tile in each direction - nTileX,nTileY,nTileZ : Number of tiles in each direction - dx,dy,dz : Element widths in each direction
Total elements: nX = nxPerTilenTileX, etc. Domain lengths: Lx = dxnX, etc.
Local side numbering is 1=bottom, 2=south, 3=east, 4=north, 5=west, 6=top. Because the mesh is uniform and axis-aligned, each periodic face pair has the same orientation as the corresponding interior face pair generated by UniformStructuredMesh, so the side "flip" index is 0 throughout (the identity face mapping in SideExchange). This matches the convention used by UniformStructuredMesh, which assigns flip 0 to every interior face and never calls RecalculateFlip.
Global side IDs are assigned analytically so that the two faces of every periodic pair share a single ID. With nXYZ = nXnYnZ: x-faces (east/west) : 1 .. nXYZ y-faces (south/north) : nXYZ+1 .. 2nXYZ z-faces (bottom/top) : 2nXYZ+1 .. 3nXYZ The east face of the element at global position (gx,gy,gz) and the west face of its +x neighbor both reference xface(gx,gy,gz); the periodic wrap identifies gx=nX (east of the last element) with the west face of the first. There are therefore exactly 3nElem unique sides and no domain boundaries.
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| class(Mesh3D_t), | intent(out) | :: | this | |||
| integer, | intent(in) | :: | nxPerTile | |||
| integer, | intent(in) | :: | nyPerTile | |||
| integer, | intent(in) | :: | nzPerTile | |||
| integer, | intent(in) | :: | nTileX | |||
| integer, | intent(in) | :: | nTileY | |||
| integer, | intent(in) | :: | nTileZ | |||
| real(kind=prec), | intent(in) | :: | dx | |||
| real(kind=prec), | intent(in) | :: | dy | |||
| real(kind=prec), | intent(in) | :: | dz |
subroutine UniformPeriodicMesh_Mesh3D_t(this,nxPerTile,nyPerTile,nzPerTile, &
nTileX,nTileY,nTileZ,dx,dy,dz)
!!
!! Create a fully triply-periodic structured hexahedral mesh and store it in
!! SELF's unstructured mesh format. Element geometry and ordering are identical
!! to UniformStructuredMesh; the only difference is connectivity. The faces on
!! the six domain boundaries are wired as interior faces whose neighbor is the
!! element on the opposite side of the domain. This realises the triply
!! periodic box T^3 = [0,Lx] x [0,Ly] x [0,Lz] required by, e.g., the
!! Arnold-Beltrami-Childress (ABC) flow benchmark.
!!
!! Input
!! - this : Fresh/empty mesh3d_t object
!! - nxPerTile,nyPerTile,nzPerTile : Elements per tile in each direction
!! - nTileX,nTileY,nTileZ : Number of tiles in each direction
!! - dx,dy,dz : Element widths in each direction
!!
!! Total elements: nX = nxPerTile*nTileX, etc.
!! Domain lengths: Lx = dx*nX, etc.
!!
!! Connectivity notes
!! ------------------
!! Local side numbering is 1=bottom, 2=south, 3=east, 4=north, 5=west, 6=top.
!! Because the mesh is uniform and axis-aligned, each periodic face pair has
!! the same orientation as the corresponding interior face pair generated by
!! UniformStructuredMesh, so the side "flip" index is 0 throughout (the
!! identity face mapping in SideExchange). This matches the convention used by
!! UniformStructuredMesh, which assigns flip 0 to every interior face and never
!! calls RecalculateFlip.
!!
!! Global side IDs are assigned analytically so that the two faces of every
!! periodic pair share a single ID. With nXYZ = nX*nY*nZ:
!! x-faces (east/west) : 1 .. nXYZ
!! y-faces (south/north) : nXYZ+1 .. 2*nXYZ
!! z-faces (bottom/top) : 2*nXYZ+1 .. 3*nXYZ
!! The east face of the element at global position (gx,gy,gz) and the west face
!! of its +x neighbor both reference xface(gx,gy,gz); the periodic wrap
!! identifies gx=nX (east of the last element) with the west face of the first.
!! There are therefore exactly 3*nElem unique sides and no domain boundaries.
!!
implicit none
class(Mesh3D_t),intent(out) :: this
integer,intent(in) :: nxPerTile
integer,intent(in) :: nyPerTile
integer,intent(in) :: nzPerTile
integer,intent(in) :: nTileX
integer,intent(in) :: nTileY
integer,intent(in) :: nTileZ
real(prec),intent(in) :: dx
real(prec),intent(in) :: dy
real(prec),intent(in) :: dz
! Local
integer :: nX,nY,nZ,nGeo,nBCs
integer :: nGlobalElem
integer :: nUniqueSides
integer :: nUniqueNodes
integer :: nLocalElems
integer :: nLocalSides
integer :: nLocalNodes
integer :: nXYZ
real(prec),allocatable :: nodeCoords(:,:,:,:,:)
integer,allocatable :: globalNodeIDs(:,:,:,:)
integer,allocatable :: sideInfo(:,:,:)
integer :: i,j,k,ti,tj,tk
integer :: ix,iy,iz,iel
integer :: ni,nj,nk
integer :: e1,e2
integer :: gx,gy,gz,gxm,gxp,gym,gyp,gzm,gzp
call this%decomp%init()
nX = nTileX*nxPerTile
nY = nTileY*nyPerTile
nZ = nTileZ*nzPerTile
nGeo = 1 ! Force the geometry to be linear
nBCs = 1 ! No domain boundaries; allocate a single (unused) slot
nXYZ = nX*nY*nZ
nGlobalElem = nX*nY*nZ
nUniqueSides = 3*nGlobalElem ! periodic: nXYZ x-faces + nXYZ y-faces + nXYZ z-faces
nUniqueNodes = (nX+1)*(nY+1)*(nZ+1)
allocate(nodeCoords(1:3,1:nGeo+1,1:nGeo+1,1:nGeo+1,1:nGlobalElem))
allocate(globalNodeIDs(1:nGeo+1,1:nGeo+1,1:nGeo+1,1:nGlobalElem))
allocate(sideInfo(1:5,1:6,1:nGlobalElem))
! Node coordinates and global node IDs (identical to UniformStructuredMesh)
do tk = 1,nTileZ
do tj = 1,nTileY
do ti = 1,nTileX
do k = 1,nzPerTile
iz = k+nzPerTile*(tk-1)
do j = 1,nyPerTile
iy = j+nyPerTile*(tj-1)
do i = 1,nxPerTile
iel = elementid(i,j,k,ti,tj,tk, &
nxpertile,nypertile,nzpertile, &
ntilex,ntiley,ntilez)
ix = i+nxPerTile*(ti-1)
do nk = 1,nGeo+1
do nj = 1,nGeo+1
do ni = 1,nGeo+1
nodeCoords(1,ni,nj,nk,iel) = real(ni-1+ix-1,prec)*dx
nodeCoords(2,ni,nj,nk,iel) = real(nj-1+iy-1,prec)*dy
nodeCoords(3,ni,nj,nk,iel) = real(nk-1+iz-1,prec)*dz
globalNodeIDs(ni,nj,nk,iel) = ni-1+i+(nxPerTile+1)*( &
nj-1+j-1+(nyPerTile+1)*( &
nk-1+k-1+(nzPerTile+1)*( &
(ti-1+nTileX*( &
tj-1+nTileY*(tk-1))))))
enddo
enddo
enddo
enddo
enddo
enddo
enddo
enddo
enddo
! Fill in face information with full triple periodicity.
! sideInfo(1:5,iSide,iEl)
! 1 - Side Type (currently unused in SELF)
! 2 - Global Side ID (shared by both faces of a periodic pair)
! 3 - Neighbor Element ID (always > 0; periodic mesh has no boundaries)
! 4 - 10*( neighbor local side ) + flip (flip is 0 for axis-aligned wrap)
! 5 - Boundary Condition ID (always 0; interior face)
do gz = 1,nZ
do gy = 1,nY
do gx = 1,nX
iel = gid2eid(gx,gy,gz,nxPerTile,nyPerTile,nzPerTile, &
nTileX,nTileY,nTileZ)
! +/- neighbor positions with periodic wrap (1..n)
gxp = modulo(gx,nX)+1; gxm = modulo(gx-2,nX)+1
gyp = modulo(gy,nY)+1; gym = modulo(gy-2,nY)+1
gzp = modulo(gz,nZ)+1; gzm = modulo(gz-2,nZ)+1
! bottom (s1=1) <-> top (s2=6) of the -z neighbor
sideInfo(1,1,iel) = 0
sideInfo(2,1,iel) = 2*nXYZ+gzm+nZ*((gx-1)+nX*(gy-1)) ! zface(gx,gy,gzm)
sideInfo(3,1,iel) = gid2eid(gx,gy,gzm,nxPerTile,nyPerTile,nzPerTile, &
nTileX,nTileY,nTileZ)
sideInfo(4,1,iel) = 10*6
sideInfo(5,1,iel) = 0
! south (s1=2) <-> north (s2=4) of the -y neighbor
sideInfo(1,2,iel) = 0
sideInfo(2,2,iel) = nXYZ+gym+nY*((gx-1)+nX*(gz-1)) ! yface(gx,gym,gz)
sideInfo(3,2,iel) = gid2eid(gx,gym,gz,nxPerTile,nyPerTile,nzPerTile, &
nTileX,nTileY,nTileZ)
sideInfo(4,2,iel) = 10*4
sideInfo(5,2,iel) = 0
! east (s1=3) <-> west (s2=5) of the +x neighbor
sideInfo(1,3,iel) = 0
sideInfo(2,3,iel) = gx+nX*((gy-1)+nY*(gz-1)) ! xface(gx,gy,gz)
sideInfo(3,3,iel) = gid2eid(gxp,gy,gz,nxPerTile,nyPerTile,nzPerTile, &
nTileX,nTileY,nTileZ)
sideInfo(4,3,iel) = 10*5
sideInfo(5,3,iel) = 0
! north (s1=4) <-> south (s2=2) of the +y neighbor
sideInfo(1,4,iel) = 0
sideInfo(2,4,iel) = nXYZ+gy+nY*((gx-1)+nX*(gz-1)) ! yface(gx,gy,gz)
sideInfo(3,4,iel) = gid2eid(gx,gyp,gz,nxPerTile,nyPerTile,nzPerTile, &
nTileX,nTileY,nTileZ)
sideInfo(4,4,iel) = 10*2
sideInfo(5,4,iel) = 0
! west (s1=5) <-> east (s2=3) of the -x neighbor
sideInfo(1,5,iel) = 0
sideInfo(2,5,iel) = gxm+nX*((gy-1)+nY*(gz-1)) ! xface(gxm,gy,gz)
sideInfo(3,5,iel) = gid2eid(gxm,gy,gz,nxPerTile,nyPerTile,nzPerTile, &
nTileX,nTileY,nTileZ)
sideInfo(4,5,iel) = 10*3
sideInfo(5,5,iel) = 0
! top (s1=6) <-> bottom (s2=1) of the +z neighbor
sideInfo(1,6,iel) = 0
sideInfo(2,6,iel) = 2*nXYZ+gz+nZ*((gx-1)+nX*(gy-1)) ! zface(gx,gy,gz)
sideInfo(3,6,iel) = gid2eid(gx,gy,gzp,nxPerTile,nyPerTile,nzPerTile, &
nTileX,nTileY,nTileZ)
sideInfo(4,6,iel) = 10*1
sideInfo(5,6,iel) = 0
enddo
enddo
enddo
call this%decomp%GenerateDecomposition(nGlobalElem,nUniqueSides)
e1 = this%decomp%offsetElem(this%decomp%rankId+1)+1
e2 = this%decomp%offsetElem(this%decomp%rankId+2)
nLocalElems = e2-e1+1
nLocalSides = nLocalElems*6
nLocalNodes = nLocalElems*8
call this%Init(nGeo,nLocalElems,nLocalSides,nLocalNodes,nBCs)
this%nUniqueSides = nUniqueSides
this%quadrature = UNIFORM
this%nodeCoords(1:3,1:nGeo+1,1:nGeo+1,1:nGeo+1,1:nLocalElems) = nodeCoords(1:3,1:nGeo+1,1:nGeo+1,1:nGeo+1,e1:e2)
this%globalNodeIDs(1:nGeo+1,1:nGeo+1,1:nGeo+1,1:nLocalElems) = globalNodeIDs(1:nGeo+1,1:nGeo+1,1:nGeo+1,e1:e2)
this%sideInfo(1:5,1:6,1:nLocalElems) = sideInfo(1:5,1:6,e1:e2)
deallocate(nodeCoords)
deallocate(globalNodeIDs)
deallocate(sideInfo)
call this%UpdateDevice()
endsubroutine UniformPeriodicMesh_Mesh3D_t