! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ! ! ! Maintainers : support@fluidnumerics.com ! Official Repository : https://github.com/FluidNumerics/self/ ! ! Copyright © 2024 Fluid Numerics LLC ! ! Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: ! ! 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. ! ! 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in ! the documentation and/or other materials provided with the distribution. ! ! 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from ! this software without specific prior written permission. ! ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS “AS IS” AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT ! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT ! HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT ! LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY ! THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF ! THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. ! ! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ! module SELF_Mesh_3D_t use SELF_Constants use SELF_Lagrange use SELF_Quadrature use SELF_SupportRoutines use SELF_HDF5 use SELF_Mesh use SELF_Mesh_2D_t,only:Mesh2D_t use SELF_DomainDecomposition ! External Libs ! use HDF5 use iso_c_binding implicit none ! ========================================================================= ! ! Node, Edge, Face, Element and Connectivity Standard ! ========================================================================= ! ! ! To define the element corner nodes, the side order and side connectivity, ! we follow the standard from CGNS SIDS (CFD General Notation System, ! Standard Interface Data Structures, http: //cgns.sourceforge.net/ ). ! ! Computational coordinate directions are defined as follows ! ! xi1 direction points from "West" (xi1=-1) to "East" (xi1=1) ! xi2 direction points from "South" (xi2=-1) to "North" (xi2=1) ! xi3 direction points from "Bottom" (xi3=-1) to "Top" (xi3=1) ! ! 3-D Hexahedreal Element sides are defined as ! ! Side 1 = Bottom (xi3 = -1) = [CN1, CN4, CN3, CN2] ! Side 2 = South (xi2 = -1) = [CN1, CN2, CN6, CN5] ! Side 3 = East (xi1 = 1) = [CN2, CN3, CN7, CN6] ! Side 4 = North (xi2 = 1) = [CN3, CN4, CN8, CN7] ! Side 5 = West (xi1 = -1) = [CN1, CN5, CN8, CN4] ! Side 6 = Top (xi3 = 1) = [CN5, CN6, CN7, CN8] ! ! In 3-D, corner nodes are order counter-clockwise (looking in the -xi3 direction) from ! bottom to top. ! ! CornerNode 1 = Bottom-South-West = (-1,-1,-1) ! CornerNode 2 = Bottom-South-East = ( 1,-1,-1) ! CornerNode 3 = Bottom-North-East = ( 1, 1,-1) ! CornerNode 4 = Bottom-North-West = (-1, 1,-1) ! CornerNode 5 = Top-South-West = (-1,-1, 1) ! CornerNode 6 = Top-South-East = ( 1,-1, 1) ! CornerNode 7 = Top-North-East = ( 1, 1, 1) ! CornerNode 8 = Top-North-West = (-1, 1, 1) ! ! ! Notes: ! * cornerNode attributes have not been implemented yet ! ! * For line segments, quads, and hexes, SELF uses Legendre-Gauss-Lobatto quadrature ! ! ! Connectivity information ! ! sideInfo(1:5,iSide,iEl) ! ! 1 - Side Type ! 2 - Global Side ID ! 3 - Neighbor Element ID ! 4 - 10*( neighbor local side ) + flip ! 5 - Boundary Condition ID ! ! ! ========================================================================= ! ! Side Ordering integer,parameter :: selfSide3D_Bottom = 1 integer,parameter :: selfSide3D_South = 2 integer,parameter :: selfSide3D_East = 3 integer,parameter :: selfSide3D_North = 4 integer,parameter :: selfSide3D_West = 5 integer,parameter :: selfSide3D_Top = 6 type,extends(SEMMesh) :: Mesh3D_t integer,pointer,dimension(:,:,:) :: sideInfo real(prec),pointer,dimension(:,:,:,:,:) :: nodeCoords integer,pointer,dimension(:,:) :: elemInfo integer,pointer,dimension(:,:,:,:) :: globalNodeIDs integer,pointer,dimension(:,:) :: CGNSCornerMap integer,pointer,dimension(:,:) :: sideMap integer,pointer,dimension(:,:) :: CGNSSideMap integer,pointer,dimension(:,:) :: BCType character(LEN=255),allocatable :: BCNames(:) ! 2:1 nonconforming (mortar) interface bookkeeping. ! ! A mortar interface joins one "big" element face to the four quarter-faces of its ! 2:1 refined neighbors ("small" faces). Faces that participate in a mortar carry ! sideInfo(3) = 0 (no conforming neighbor) and sideInfo(5) = 0 (no boundary ! condition), so all conforming-face machinery (SideExchange, RecalculateFlip, ! boundary-condition mapping) skips them; sideInfo(1) is set to the mortar index ! for reference. ! ! The big face is split into four sub-faces ("quadrants") indexed in the big ! face's trace coordinates: q = kx + 2*(ky-1) with kx, ky in {1,2} the ! half-interval index along the face's first and second coordinate ! (1 = [-1,0], 2 = [0,1]). ! ! mortarInfo(1:14,1:nMortars) is replicated on all ranks (element ids are global): ! 1 - big element id ; 2 - big local face id ! 3 - small element id on the sub-face covering big-face quadrant 1 ! 4 - 10*(small local face) + flip for that sub-face ! 5,6 - same for quadrant 2 ! 7,8 - same for quadrant 3 ! 9,10 - same for quadrant 4 ! 11:14 - global side ids of sub-faces 1..4 (used as MPI message tags) ! The flip maps big-face quadrant coordinates to the small face's own coordinates ! (the same receiver-to-donor convention as sideInfo(4); see MortarFaceMap). integer :: nMortars = 0 integer,pointer,dimension(:,:) :: mortarInfo => null() ! Anisotropic (2-subface) nonconforming interfaces. ! ! An anisotropic mortar joins one big face to two small faces that split it along ! one of its two trace coordinates only (the split coordinate) and span the whole ! interval [-1,1] of the other. This is the interface between a coarse column and two ! fine columns of an extruded mesh whose surface mesh has 2-D mortars. Its faces ! carry sideInfo(3) = 0 and sideInfo(5) = 0 like those of a 2x2 mortar, and ! sideInfo(1) = nMortars + a, so a nonzero sideInfo(1) still marks a mortar face. ! ! anisoMortarInfo(1:9,1:nAnisoMortars) is replicated on all ranks (global element ids): ! 1 - big element id ; 2 - big local face id ! 3 - small element id on the half h = 1 ([-1,0] of the split coordinate) ! 4 - 10*(small local face) + flip for that half ! 5,6 - same for the half h = 2 ([0,1]) ! 7,8 - global side ids of the two small faces (used as MPI message tags) ! 9 - the split coordinate of the big face, 1 or 2 (the face trace-coordinate ! convention of BoundaryInterp: faces 1,6 -> (xi1,xi2); 2,4 -> (xi1,xi3); ! 3,5 -> (xi2,xi3)) ! The flip maps big-face half coordinates to the small face's own coordinates, as for ! mortarInfo. The 1-D mortarR/mortarP of the half act on the split coordinate and the ! identity on the other. integer :: nAnisoMortars = 0 integer,pointer,dimension(:,:) :: anisoMortarInfo => null() ! Material tracking: every element has an integer material id ! indexing into materialNames. Single-material constructors and ! readers (HOPr, structured, periodic) leave nMaterials = 1 with ! the name "default". The HOHQMesh ISM-MM reader populates the ! table with the material strings from the .mesh file. integer :: nMaterials = 0 integer,allocatable :: elemMaterial(:) character(LEN=SELF_MESH_MATNAME_LENGTH),allocatable :: materialNames(:) ! Layer count of a mesh built by ExtrudedMesh, and 0 for every other mesh. ! ! An extruded mesh stores its elements column by column: local element ! e3 = (e2-1)*nLayers + k is layer k (k = 1 at the bottom) of the column above local ! element e2 of the surface Mesh2D it was built from, and the same rule maps global ! ids. A column therefore never straddles two ranks, and the surface element and ! layer of an element follow from its id alone: ! e2 = (e3-1)/nLayers + 1 , k = e3 - (e2-1)*nLayers integer :: nLayers = 0 contains procedure,public :: Init => Init_Mesh3D_t procedure,public :: Free => Free_Mesh3D_t procedure,public :: UpdateDevice => UpdateDevice_Mesh3D_t generic,public :: StructuredMesh => UniformStructuredMesh_Mesh3D_t procedure,private :: UniformStructuredMesh_Mesh3D_t procedure,public :: SimpleMortarMesh => SimpleMortarMesh_Mesh3D_t procedure,public :: SimpleAnisoMortarMesh => SimpleAnisoMortarMesh_Mesh3D_t generic,public :: PeriodicStructuredMesh => UniformPeriodicMesh_Mesh3D_t procedure,private :: UniformPeriodicMesh_Mesh3D_t procedure,public :: ExtrudedMesh => ExtrudedMesh_Mesh3D_t procedure,public :: Read_HOPr => Read_HOPr_Mesh3D_t procedure,public :: Read_HOHQMesh => Read_HOHQMesh_Mesh3D_t procedure,public :: ResetBoundaryConditionType => ResetBoundaryConditionType_Mesh3D_t procedure,public :: Write_Mesh => Write_Mesh3D_t procedure,public :: RecalculateFlip => RecalculateFlip_Mesh3D_t endtype Mesh3D_t integer,private :: CGNStoSELFflip(1:6,1:6,1:4) private :: GeometryNodes1D ! This table maps the primary side, secondary side, and CGNS flip values ! to indexing flips that are used in SELF. ! This table is used after reading in HOPr mesh information in "RecalculateFlip" ! SELF's flip indices correspond to the following scenarios ! ! 0 i2 = i1 j2 = j1 ! 1 i2 = N-i1 j2 = j1 ! 2 i2 = N-i1 j2 = N-j1 ! 3 i2 = i1 j2 = N-j1 ! 4 i2 = j1 j2 = i1 ! 5 i2 = N-j1 j2 = i1 ! 6 i2 = N-j1 j2 = N-i1 ! 7 i2 = j1 j2 = N-i1 ! data CGNStoSELFflip/ & 4,0,0,1,4,0, & 0,4,4,5,0,4, & 0,4,4,5,0,4, & 1,7,7,6,1,7, & 4,0,0,1,4,0, & 0,4,4,5,0,4, & 3,5,5,4,3,5, & 7,1,1,0,7,1, & 7,1,1,0,7,1, & 4,0,0,1,4,0, & 3,5,5,4,3,5, & 7,1,1,0,7,1, & 6,2,2,3,6,2, & 2,6,6,7,2,6, & 2,6,6,7,2,6, & 3,5,5,4,3,5, & 6,2,2,3,6,2, & 2,6,6,7,2,6, & 1,7,7,6,1,7, & 5,3,3,2,5,3, & 5,3,3,2,5,3, & 6,2,2,3,6,2, & 1,7,7,6,1,7, & 5,3,3,2,5,3/ ! Half-interval index along the big face's first (kx) and second (ky) trace ! coordinate for each mortar sub-face quadrant q = kx + 2*(ky-1). The 1-D mortar ! restriction/projection matrices (Lagrange % mortarR/mortarP) are applied as the ! tensor product R_kx x R_ky (P_kx x P_ky) on sub-face q. integer,parameter,public :: mortarQuadKx(1:4) = [1,2,1,2] integer,parameter,public :: mortarQuadKy(1:4) = [1,1,2,2] contains pure subroutine MortarFaceMap(i,j,N,flip,i2,j2) !! Maps receiver-face node indices (i,j) to donor-face node indices (i2,j2) for !! each of the eight SELF face flips (see the CGNStoSELFflip table above; indices !! here are 1-based, so the 0-based rule i2 = N-i1 reads i2 = N+2-i). Mortar !! staging uses this to reorient a small face's trace into the big face's !! coordinates and to scatter restricted traces back. implicit none integer,intent(in) :: i,j,N,flip integer,intent(out) :: i2,j2 select case(flip) case(0) i2 = i j2 = j case(1) i2 = N+2-i j2 = j case(2) i2 = N+2-i j2 = N+2-j case(3) i2 = i j2 = N+2-j case(4) i2 = j j2 = i case(5) i2 = N+2-j j2 = i case(6) i2 = N+2-j j2 = N+2-i case default ! flip == 7 i2 = j j2 = N+2-i endselect endsubroutine MortarFaceMap subroutine Init_Mesh3D_t(this,nGeo,nElem,nSides,nNodes,nBCs) implicit none class(Mesh3D_t),intent(inout) :: this integer,intent(in) :: nGeo integer,intent(in) :: nElem integer,intent(in) :: nSides integer,intent(in) :: nNodes integer,intent(in) :: nBCs this%nElem = nElem this%nGlobalElem = nElem this%nGeo = nGeo this%nSides = nSides this%nNodes = nNodes this%nCornerNodes = 0 this%nUniqueSides = 0 this%nUniqueNodes = 0 this%nBCs = nBCs allocate(this%elemInfo(1:6,1:nElem)) allocate(this%sideInfo(1:5,1:6,1:nElem)) allocate(this%nodeCoords(1:3,1:nGeo+1,1:nGeo+1,1:nGeo+1,1:nElem)) allocate(this%globalNodeIDs(1:nGeo+1,1:nGeo+1,1:nGeo+1,1:nElem)) allocate(this%CGNSCornerMap(1:3,1:8)) allocate(this%CGNSSideMap(1:4,1:6)) allocate(this%sideMap(1:4,1:6)) allocate(this%BCType(1:4,1:nBCs)) allocate(this%BCNames(1:nBCs)) ! Default material table: a single "default" material covers all ! elements. Readers that carry material information overwrite ! these allocations with the per-file table. this%nMaterials = 1 allocate(this%elemMaterial(1:nElem)) allocate(this%materialNames(1:1)) this%elemMaterial = 1 this%materialNames(1) = "default" ! Create lookup tables to assist with connectivity generation this%CGNSCornerMap(1:3,1) = (/1,1,1/) ! Bottom-South-West this%CGNSCornerMap(1:3,2) = (/nGeo+1,1,1/) ! Bottom-South-East this%CGNSCornerMap(1:3,3) = (/nGeo+1,nGeo+1,1/) ! Bottom-North-East this%CGNSCornerMap(1:3,4) = (/1,nGeo+1,1/) ! Bottom-North-West this%CGNSCornerMap(1:3,5) = (/1,1,nGeo+1/) ! Top-South-West this%CGNSCornerMap(1:3,6) = (/nGeo+1,1,nGeo+1/) ! Top-South-East this%CGNSCornerMap(1:3,7) = (/nGeo+1,nGeo+1,nGeo+1/) ! Top-North-East this%CGNSCornerMap(1:3,8) = (/1,nGeo+1,nGeo+1/) ! Top-North-West ! Maps from local corner node id to CGNS side this%CGNSSideMap(1:4,1) = (/1,4,3,2/) this%CGNSSideMap(1:4,2) = (/1,2,6,5/) this%CGNSSideMap(1:4,3) = (/2,3,7,6/) this%CGNSSideMap(1:4,4) = (/3,4,8,7/) this%CGNSSideMap(1:4,5) = (/1,5,8,4/) this%CGNSSideMap(1:4,6) = (/5,6,7,8/) ! Sidemap traverses each face so that the normal ! formed by the right hand rule is the coordinate ! positive pointing normal. For east,north,and top ! this is an outward facing normal. ! For bottom, south, and west, the normal is inward ! facing. this%sideMap(1:4,1) = (/1,2,3,4/) ! Bottom this%sideMap(1:4,2) = (/1,2,6,5/) ! South this%sideMap(1:4,3) = (/2,3,7,6/) ! East this%sideMap(1:4,4) = (/4,3,7,8/) ! North this%sideMap(1:4,5) = (/1,4,8,5/) ! West this%sideMap(1:4,6) = (/5,6,7,8/) ! Top endsubroutine Init_Mesh3D_t subroutine Free_Mesh3D_t(this) implicit none class(Mesh3D_t),intent(inout) :: this this%nElem = 0 this%nSides = 0 this%nNodes = 0 this%nCornerNodes = 0 this%nUniqueSides = 0 this%nUniqueNodes = 0 this%nBCs = 0 deallocate(this%elemInfo) deallocate(this%sideInfo) deallocate(this%nodeCoords) deallocate(this%globalNodeIDs) deallocate(this%CGNSCornerMap) deallocate(this%sideMap) deallocate(this%CGNSSideMap) deallocate(this%BCType) deallocate(this%BCNames) if(associated(this%mortarInfo)) then deallocate(this%mortarInfo) this%mortarInfo => null() endif this%nMortars = 0 if(associated(this%anisoMortarInfo)) then deallocate(this%anisoMortarInfo) this%anisoMortarInfo => null() endif this%nAnisoMortars = 0 if(allocated(this%elemMaterial)) deallocate(this%elemMaterial) if(allocated(this%materialNames)) deallocate(this%materialNames) this%nMaterials = 0 this%nLayers = 0 call this%decomp%Free() endsubroutine Free_Mesh3D_t subroutine UpdateDevice_Mesh3D_t(this) implicit none class(Mesh3D_t),intent(inout) :: this if(.false.) this%nElem = this%nElem ! CPU stub; suppress unused-dummy-argument warning endsubroutine UpdateDevice_Mesh3D_t subroutine ResetBoundaryConditionType_Mesh3D_t(this,bcid) !! This method can be used to reset all of the boundary elements !! boundary condition type to the desired value. !! !! Note that ALL physical boundaries will be set to have this boundary !! condition implicit none class(Mesh3D_t),intent(inout) :: this integer,intent(in) :: bcid ! Local integer :: iSide,iEl,e2 logical :: skipMortars skipMortars = this%nMortars+this%nAnisoMortars > 0 do iEl = 1,this%nElem do iSide = 1,6 e2 = this%sideInfo(3,iSide,iEl) ! Mortar faces carry sideInfo(3) = 0 but are interior faces, with sideInfo(1) ! holding the mortar index; they never take a physical boundary condition. The ! HOPr reader copies sideInfo(1) verbatim from the file, where it is the HOPr ! side type and may be nonzero on an ordinary face, so sideInfo(1) is only a ! mortar marker on a mesh that actually carries mortars. Same test as the ! unmapped-boundary scan in SELF_DGModel{2,3}D_t. if(e2 == 0 .and. .not.(skipMortars .and. this%sideInfo(1,iSide,iEl) /= 0)) then this%sideInfo(5,iSide,iEl) = bcid endif enddo enddo call this%UpdateDevice() endsubroutine ResetBoundaryConditionType_Mesh3D_t subroutine RecalculateFlip_Mesh3D_t(this) implicit none class(Mesh3D_t),intent(inout) :: this ! Local integer :: e1 integer :: s1 integer :: e2 integer :: s2 integer :: cgnsFlip,selfFlip do e1 = 1,this%nElem do s1 = 1,6 e2 = this%sideInfo(3,s1,e1) s2 = this%sideInfo(4,s1,e1)/10 cgnsFlip = this%sideInfo(4,s1,e1)-s2*10 if(e2 /= 0) then selfFlip = CGNStoSELFflip(s2,s1,cgnsFlip) this%sideInfo(4,s1,e1) = 10*s2+selfFlip endif enddo enddo endsubroutine RecalculateFlip_Mesh3D_t pure function elementid(i,j,k,ti,tj,tk,nxpertile,nypertile,nzpertile, & ntilex,ntiley,ntilez) result(eid) implicit none integer,intent(in) :: i,j,k integer,intent(in) :: ti,tj,tk integer,intent(in) :: nxpertile,nypertile,nzpertile integer,intent(in) :: ntilex,ntiley,ntilez integer :: eid eid = i+nxpertile*(j-1+nypertile*(k-1+nzpertile*( & ti-1+ntilex*(tj-1+ntiley*(tk-1))))) if(.false.) eid = eid+ntilez ! suppress unused-dummy-argument warning endfunction elementid pure function gid2eid(gx,gy,gz,nxpertile,nypertile,nzpertile, & ntilex,ntiley,ntilez) result(eid) !! Map a global element position (gx,gy,gz), with gx in [1,nX] etc., to the !! tile-ordered element id used throughout the structured mesh. This is the !! inverse of the (i,ti) decomposition: gx = i + nxpertile*(ti-1). implicit none integer,intent(in) :: gx,gy,gz integer,intent(in) :: nxpertile,nypertile,nzpertile integer,intent(in) :: ntilex,ntiley,ntilez integer :: eid integer :: i,j,k,ti,tj,tk ti = (gx-1)/nxpertile+1; i = gx-(ti-1)*nxpertile tj = (gy-1)/nypertile+1; j = gy-(tj-1)*nypertile tk = (gz-1)/nzpertile+1; k = gz-(tk-1)*nzpertile eid = elementid(i,j,k,ti,tj,tk,nxpertile,nypertile,nzpertile, & ntilex,ntiley,ntilez) endfunction gid2eid subroutine UniformStructuredMesh_Mesh3D_t(this,nxPerTile,nyPerTile,nzPerTile, & nTileX,nTileY,nTileZ,dx,dy,dz,bcids,comm) !! !! Create a structured mesh and store it in SELF's unstructured mesh format. !! The mesh is created in tiles of size (tnx,tny,tnz). Tiling is used to determine !! the element ordering. !! !! !! Input !! - this : Fresh/empty mesh2d_t object !! - nxPerTile : The number of elements in the x direction within a tile !! - nyPerTile : The number of elements in the y direction within a tile !! - nzPerTile : The number of elements in the z direction within a tile !! - nTileX : The number of tiles in the x direction !! - nTileY : The number of tiles in the y direction !! - nTileZ : The number of tiles in the z direction !! - dx : Element width in the x-direction !! - dy : Element width in the y-direction !! - dz : Element width in the z-direction !! - bcids(1:6) : Boundary condition flags for the south, east, north, and west sides of the domain !! - comm (optional) : Externally managed MPI communicator (Fortran integer handle). When !! provided, the caller owns MPI initialization/finalization. !! !! Output !! - this : mesh2d_t object with vertices, faces, and element information !! !! Total number of elements in the x-direction is nX = nxPerTile*nTileX !! Total number of elements in the y-direction is nY = nyPerTile*nTileY !! !! Length of the domain in the x-direction is Lx = dx*nX !! Length of the domain in the y-direction is Ly = dy*nY !! 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 integer,intent(in) :: bcids(1:6) integer,intent(in),optional :: comm ! Local integer :: nX,nY,nZ,nGeo,nBCs integer :: nGlobalElem integer :: nUniqueSides integer :: nUniqueNodes integer :: nLocalElems integer :: nLocalSides integer :: nLocalNodes 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,s1,s2 integer :: nfaces call this%decomp%init(comm) nX = nTileX*nxPerTile nY = nTileY*nyPerTile nZ = nTileZ*nzPerTile nGeo = 1 ! Force the geometry to be linear nBCs = 6 ! Force the number of boundary conditions to 4 nGlobalElem = nX*nY*nZ nUniqueSides = (nX+1)*nY*nZ+(nY+1)*nX*nZ+(nZ+1)*nX*nY 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)) ! Row 1 is the HOPr side type, which SELF repurposes as the mortar index; a ! structured mesh has no mortars, so it must read as 0 rather than as whatever ! the allocation happened to pick up. The mortar and periodic builders already ! zero the whole table for the same reason. sideInfo(1,1:6,1:nGlobalElem) = 0 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 ! sideInfo(1:5,iSide,iEl) ! 1 - Side Type (currently unused in SELF) ! 2 - Global Side ID (Used for message passing) ! 3 - Neighbor Element ID ! 4 - 10*( neighbor local side ) + flip ! 5 - Boundary Condition ID nfaces = 0 do tk = 1,nTileZ do tj = 1,nTileY do ti = 1,nTileX do k = 1,nzPerTile do j = 1,nyPerTile do i = 1,nxPerTile iel = elementid(i,j,k,ti,tj,tk, & nxpertile,nypertile,nzpertile, & ntilex,ntiley,ntilez) ! bottom, iside=1 s1 = 1 s2 = 6 if(k == 1) then ! bottom most part of the tile if(tk == 1) then ! bottom most tile nfaces = nfaces+1 sideinfo(2,s1,iel) = nfaces sideinfo(3,s1,iel) = 0 ! Neigbor element (null, boundary condition) sideinfo(4,s1,iel) = 0 ! Neighbor side id (null, boundary condition) sideinfo(5,s1,iel) = bcids(s1) ! Boundary condition id; set from the user input else ! interior tile !neighbor element is the top most element in the tile beneath e2 = elementid(i,j,nzpertile,ti,tj,tk-1, & nxpertile,nypertile,nzpertile, & ntilex,ntiley,ntilez) sideinfo(2,s1,iel) = sideInfo(2,s2,e2) ! Copy the face id from neighbor sideinfo(3,s1,iel) = e2 sideinfo(4,s1,iel) = 10*s2 ! Neighbor side id sideinfo(5,s1,iel) = 0 ! Boundary condition id; (null, interior face) endif else ! interior to the tile !neighbor element is in the same tile, but beneath e2 = elementid(i,j,k-1,ti,tj,tk, & nxpertile,nypertile,nzpertile, & ntilex,ntiley,ntilez) sideinfo(2,s1,iel) = sideInfo(2,s2,e2) ! Copy the face id from neighbor sideinfo(3,s1,iel) = e2 sideinfo(4,s1,iel) = 10*s2 ! Neighbor side id sideinfo(5,s1,iel) = 0 ! Boundary condition id; (null, interior face) endif ! south, iside=2 s1 = 2 s2 = 4 ! Neighbor side is north (4) if(j == 1) then ! southern most part of the tile if(tj == 1) then ! southern most tile nfaces = nfaces+1 sideinfo(2,s1,iel) = nfaces sideinfo(3,s1,iel) = 0 ! Neigbor element (null, boundary condition) sideinfo(4,s1,iel) = 0 ! Neighbor side id (null, boundary condition) sideinfo(5,s1,iel) = bcids(s1) ! Boundary condition id; eastern boundary set from the user input else ! interior tile !neighbor element is northernmost element in the tile to the south e2 = elementid(i,nypertile,k,ti,tj-1,tk, & nxpertile,nypertile,nzpertile, & ntilex,ntiley,ntilez) sideinfo(2,s1,iel) = sideInfo(2,s2,e2) ! Copy the face id from neighbor sideinfo(3,s1,iel) = e2 ! Neigbor element sideinfo(4,s1,iel) = 10*s2 ! Neighbor side id sideinfo(5,s1,iel) = 0 ! Boundary condition id; (null, interior face) endif else ! interior to the tile !neighbor element is in the same tile, to the south e2 = elementid(i,j-1,k,ti,tj,tk, & nxpertile,nypertile,nzpertile, & ntilex,ntiley,ntilez) sideinfo(2,s1,iel) = sideInfo(2,s2,e2) ! Copy the face id from neighbor sideinfo(3,s1,iel) = e2 ! Neigbor element sideinfo(4,s1,iel) = 10*s2 ! Neighbor side id sideinfo(5,s1,iel) = 0 ! Boundary condition id; (null, interior face) endif ! east, iside=3 s1 = 3 s2 = 5 ! neighbor side id is west (5) ! East faces are always new faces, due to the way we are traversing the grid nfaces = nfaces+1 sideinfo(2,s1,iel) = nfaces if(i == nxPerTile) then ! eastern most part of the tile if(ti == nTileX) then ! eastern most tile sideinfo(3,s1,iel) = 0 ! Neigbor element (null, boundary condition) sideinfo(4,s1,iel) = 0 ! Neighbor side id (null, boundary condition) sideinfo(5,s1,iel) = bcids(s1) ! Boundary condition id; else ! interior tile !neighbor element is westernmost element in tile to the east e2 = elementid(1,j,k,ti+1,tj,tk, & nxpertile,nypertile,nzpertile, & ntilex,ntiley,ntilez) sideinfo(3,s1,iel) = e2 ! Neigbor element sideinfo(4,s1,iel) = 10*s2 ! Neighbor side id sideinfo(5,s1,iel) = 0 ! Boundary condition id; (null, interior face) endif else ! interior to the tile !neighbor element is in the same tile, to the east e2 = elementid(i+1,j,k,ti,tj,tk, & nxpertile,nypertile,nzpertile, & ntilex,ntiley,ntilez) sideinfo(3,s1,iel) = e2 ! Neigbor element sideinfo(4,s1,iel) = 10*s2 ! Neighbor side id sideinfo(5,s1,iel) = 0 ! Boundary condition id; (null, interior face) endif ! north, iside=4 s1 = 4 s2 = 2 ! neighbor side is south (2) ! North faces are always new faces, due to the way we are traversing the grid nfaces = nfaces+1 sideinfo(2,s1,iel) = nfaces if(j == nyPerTile) then ! northern most part of the tile if(tj == nTileY) then ! northern most tile sideinfo(3,s1,iel) = 0 ! Neigbor element (null, boundary condition) sideinfo(4,s1,iel) = 0 ! Neighbor side id (null, boundary condition) sideinfo(5,s1,iel) = bcids(s1) ! Boundary condition id; set from the user input else ! interior tile, but northern most face of the tile !neighbor element is the southernmost element in the tile to the north e2 = elementid(i,1,k,ti,tj+1,tk, & nxpertile,nypertile,nzpertile, & ntilex,ntiley,ntilez) sideinfo(3,s1,iel) = e2 ! Neigbor element sideinfo(4,s1,iel) = 10*s2 ! Neighbor side id sideinfo(5,s1,iel) = 0 ! Boundary condition id; (null, interior face) endif else ! interior to the tile !neighbor element is the tile to the north e2 = elementid(i,j+1,k,ti,tj,tk, & nxpertile,nypertile,nzpertile, & ntilex,ntiley,ntilez) sideinfo(3,s1,iel) = e2 ! Neigbor element sideinfo(4,s1,iel) = 10*s2 ! Neighbor side id sideinfo(5,s1,iel) = 0 ! Boundary condition id; (null, interior face) endif ! west, iside=5 s1 = 5 s2 = 3 ! neighbor side id is east (3) if(i == 1) then ! western most part of the tile if(ti == 1) then ! western most tile nfaces = nfaces+1 sideinfo(2,s1,iel) = nfaces sideinfo(3,s1,iel) = 0 ! Neigbor element (null, boundary condition) sideinfo(4,s1,iel) = 0 ! Neighbor side id (null, boundary condition) sideinfo(5,s1,iel) = bcids(s1) ! Boundary condition id else ! interior tile, but western most face of the tile !neighbor element is the easternmost element in the tile to the west e2 = elementid(nxperTile,j,k,ti-1,tj,tk, & nxpertile,nypertile,nzpertile, & ntilex,ntiley,ntilez) sideinfo(2,s1,iel) = sideInfo(2,s2,e2) ! Copy the face id from neighbor's east face sideinfo(3,s1,iel) = e2 sideinfo(4,s1,iel) = 10*s2 ! Neighbor side id - neighbor to the west, east side (2) sideinfo(5,s1,iel) = 0 ! Boundary condition id; (null, interior face) endif else ! interior to the tile !neighbor element is the element to the west in the same tile e2 = elementid(i-1,j,k,ti,tj,tk, & nxpertile,nypertile,nzpertile, & ntilex,ntiley,ntilez) sideinfo(2,s1,iel) = sideInfo(2,s2,e2) ! Copy the face id from neighbor's east face sideinfo(3,s1,iel) = e2 sideinfo(4,s1,iel) = 10*s2 ! Neighbor side id - neighbor to the west, east side (2) sideinfo(5,s1,iel) = 0 ! Boundary condition id; (null, interior face) endif ! top, iside=6 s1 = 6 s2 = 1 ! neighbor side is bottom (1) ! Top faces are always new faces, due to the way we are traversing the grid nfaces = nfaces+1 sideinfo(2,s1,iel) = nfaces if(k == nzPerTile) then ! top most part of the tile if(tk == nTileZ) then ! top most tile sideinfo(3,s1,iel) = 0 ! Neigbor element (null, boundary condition) sideinfo(4,s1,iel) = 0 ! Neighbor side id (null, boundary condition) sideinfo(5,s1,iel) = bcids(s1) ! Boundary condition id; set from the user input else ! interior tile, but top most face of the tile !neighbor element is the bottom-most element in the tile above e2 = elementid(i,j,1,ti,tj,tk+1, & nxpertile,nypertile,nzpertile, & ntilex,ntiley,ntilez) sideinfo(3,s1,iel) = e2 ! Neigbor element sideinfo(4,s1,iel) = 10*s2 ! Neighbor side id sideinfo(5,s1,iel) = 0 ! Boundary condition id; (null, interior face) endif else ! interior to the tile !neighbor element is the tile above e2 = elementid(i,j,k+1,ti,tj,tk, & nxpertile,nypertile,nzpertile, & ntilex,ntiley,ntilez) sideinfo(3,s1,iel) = e2 ! Neigbor element, inside same tile, to the north sideinfo(4,s1,iel) = 10*s2 ! Neighbor side id - neighbor to the north, south side (1) sideinfo(5,s1,iel) = 0 ! Boundary condition id; (null, interior face) endif enddo enddo enddo 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 UniformStructuredMesh_Mesh3D_t subroutine SimpleMortarMesh_Mesh3D_t(this,dx,bcids,comm,flips) !! !! Create the smallest 3D 2:1 nonconforming (mortar) mesh: one "big" element of size !! 2*dx x 2*dx x 2*dx whose east face is shared with the faces of four "small" !! dx x dx x dx elements stacked 2x2 in (y,z). The mesh is conforming everywhere !! except at the single mortar interface. !! !! z = 2*dx ________________ ______ ______ !! | | e4 | e5 | (view of the x = 2*dx plane; !! | |______|______| e2..e5 sit in the big face's !! | e1 | e2 | e3 | quadrants 1..4) !! |________________|______|______| !! x = 0 2*dx 3*dx !! !! Input !! - this : Fresh/empty Mesh3D_t object !! - dx : Edge length of the small elements; the big element has edge length 2*dx !! - bcids(1:6) : Boundary condition flags for the bottom, south, east, north, !! west, and top sides of the domain (the side ordering of !! StructuredMesh) !! - flips(1:4) (optional) : requested face flip (0..7) for the small element in !! each big-face quadrant. Each small element is rigidly rotated so !! that its face on the mortar interface meets the big face with !! the requested flip; flips 0,2,5,7 place the small element's west !! face on the interface, flips 1,3,4,6 its east face (the flip !! fixes the face chirality). Defaults to 0 (no rotation). !! !! Output !! - this : Mesh3D_t object with five elements and one mortar interface !! !! All connectivity (conforming interior faces, boundary faces, and the mortar !! table) is derived by sampling the bilinear face maps and matching physical !! positions, so the emitted sideInfo/mortarInfo are correct by construction for !! every rotation. Element geometry is trilinear (nGeo=1). With domain !! decomposition on two or more ranks the mortar interface straddles the rank !! boundary. !! implicit none class(Mesh3D_t),intent(out) :: this real(prec),intent(in) :: dx integer,intent(in) :: bcids(1:6) integer,intent(in),optional :: comm integer,intent(in),optional :: flips(1:4) ! Local integer,parameter :: nGlobalElem = 5 ! Signed local axis directions (1=x,2=y,3=z) of the rotated small elements for ! each requested flip; column f holds (d1,d2,d3) such that the face meeting the ! big east face does so with flip f. The donor face is West (5) when d1 = +x and ! East (3) when d1 = -x. integer,parameter :: flipAxes(1:3,0:7) = reshape([ & 1,2,3, & ! flip 0 : West donor -1,-2,3, & ! flip 1 : East donor 1,-2,-3, & ! flip 2 : West donor -1,2,-3, & ! flip 3 : East donor -1,3,2, & ! flip 4 : East donor 1,-3,2, & ! flip 5 : West donor -1,-3,-2, & ! flip 6 : East donor 1,3,-2],[3,8]) ! flip 7 : West donor real(prec) :: nodeCoords(1:3,1:2,1:2,1:2,1:nGlobalElem) integer :: sideInfo(1:5,1:6,1:nGlobalElem) integer :: minfo(1:14) integer :: localFlips(1:4) integer,parameter :: invFlip(0:7) = [0,1,2,3,4,7,6,5] real(prec) :: center(1:3),axis(1:3,1:3) real(prec) :: tol integer :: q,f,i,j,k,d integer :: eA,eB,sA,sB,fAB integer :: e1,e2 integer :: gid,nUniqueSides integer :: nLocalElems integer :: nGeo,nBCs logical :: matched call this%decomp%init(comm) nGeo = 1 ! Trilinear element geometry nBCs = 6 tol = 1.0e-3_prec*dx localFlips = 0 if(present(flips)) then localFlips = flips do q = 1,4 if(localFlips(q) < 0 .or. localFlips(q) > 7) then print*,"SimpleMortarMesh_Mesh3D_t : requested flip out of range 0..7" stop 1 endif enddo endif ! Element 1 (big) : [0,2dx]^3, axis aligned do k = 1,2 do j = 1,2 do i = 1,2 nodeCoords(1:3,i,j,k,1) = [real(i-1,prec),real(j-1,prec),real(k-1,prec)]* & 2.0_prec*dx enddo enddo enddo ! Elements 2..5 (small) : [2dx,3dx] x quadrant q of the big east face, each ! rigidly rotated per the requested flip do q = 1,4 f = localFlips(q) center(1) = 2.5_prec*dx center(2) = (real(mortarQuadKx(q),prec)-0.5_prec)*dx center(3) = (real(mortarQuadKy(q),prec)-0.5_prec)*dx do d = 1,3 axis(1:3,d) = 0.0_prec axis(abs(flipAxes(d,f)),d) = real(sign(1,flipAxes(d,f)),prec) enddo do k = 1,2 do j = 1,2 do i = 1,2 nodeCoords(1:3,i,j,k,1+q) = center+ & 0.5_prec*dx*(real(2*i-3,prec)*axis(1:3,1)+ & real(2*j-3,prec)*axis(1:3,2)+ & real(2*k-3,prec)*axis(1:3,3)) enddo enddo enddo enddo sideInfo = 0 ! Conforming interior faces: match every pair of element faces by sampling the ! bilinear face maps under each of the eight flips. do eA = 1,nGlobalElem do eB = eA+1,nGlobalElem do sA = 1,6 if(sideInfo(3,sA,eA) /= 0) cycle do sB = 1,6 call MatchFaces3D(nodeCoords(:,:,:,:,eA),sA, & nodeCoords(:,:,:,:,eB),sB,tol,matched,fAB) if(matched) then sideInfo(3,sA,eA) = eB sideInfo(4,sA,eA) = 10*sB+fAB sideInfo(3,sB,eB) = eA sideInfo(4,sB,eB) = 10*sA+invFlip(fAB) exit endif enddo enddo enddo enddo ! Mortar interface: big element east face against each small element. The donor ! face and flip are recovered by matching each quadrant of the big face. minfo = 0 minfo(1) = 1 minfo(2) = selfSide3D_East do q = 1,4 matched = .false. do sB = 1,6 do f = 0,7 if(QuadrantMatches3D(nodeCoords(:,:,:,:,1),selfSide3D_East,q, & nodeCoords(:,:,:,:,1+q),sB,f,tol)) then minfo(2*q+1) = 1+q minfo(2*q+2) = 10*sB+f matched = .true. exit endif enddo if(matched) exit enddo if(.not. matched) then print*,"SimpleMortarMesh_Mesh3D_t : failed to match mortar sub-face",q stop 1 endif enddo ! Physical boundary conditions: any face that is neither conforming nor on the ! mortar interface lies on one of the six domain boundary planes; identify the ! plane from the face centroid. do eA = 1,nGlobalElem do sA = 1,6 if(sideInfo(3,sA,eA) /= 0) cycle if(eA == 1 .and. sA == selfSide3D_East) cycle matched = .false. do q = 1,4 if(eA == minfo(2*q+1) .and. sA == minfo(2*q+2)/10) matched = .true. enddo if(matched) cycle ! small mortar face sideInfo(5,sA,eA) = DomainBoundaryId3D(nodeCoords(:,:,:,:,eA),sA,dx,tol,bcids) enddo enddo ! Global side ids: conforming pairs share one id; the big mortar face shares the ! id of sub-face 1; sub-faces 2..4 get their own ids (MPI message tags). gid = 0 do eA = 1,nGlobalElem do sA = 1,6 if(sideInfo(2,sA,eA) /= 0) cycle if(eA == 1 .and. sA == selfSide3D_East) cycle ! assigned with the sub-faces matched = .false. do q = 1,4 if(eA == minfo(2*q+1) .and. sA == minfo(2*q+2)/10) then matched = .true. exit endif enddo if(matched) cycle ! assigned with the sub-faces below gid = gid+1 sideInfo(2,sA,eA) = gid eB = sideInfo(3,sA,eA) if(eB /= 0) then sB = sideInfo(4,sA,eA)/10 sideInfo(2,sB,eB) = gid endif enddo enddo do q = 1,4 gid = gid+1 minfo(10+q) = gid sideInfo(2,minfo(2*q+2)/10,minfo(2*q+1)) = gid sideInfo(1,minfo(2*q+2)/10,minfo(2*q+1)) = 1 ! mortar index if(q == 1) then sideInfo(2,selfSide3D_East,1) = gid sideInfo(1,selfSide3D_East,1) = 1 ! mortar index endif enddo nUniqueSides = gid ! Domain decomposition. The message count upper bound is oversized relative to ! nUniqueSides to accommodate the per-variable (and per-direction) mortar and ! conforming side messages on this small mesh. call this%decomp%GenerateDecomposition(nGlobalElem,64*nUniqueSides) e1 = this%decomp%offsetElem(this%decomp%rankId+1)+1 e2 = this%decomp%offsetElem(this%decomp%rankId+2) nLocalElems = e2-e1+1 call this%Init(nGeo,nLocalElems,nLocalElems*6,nLocalElems*8,nBCs) this%nUniqueSides = nUniqueSides this%quadrature = UNIFORM this%BCType = 0 this%elemInfo = 0 this%globalNodeIDs = 0 this%nodeCoords(1:3,1:2,1:2,1:2,1:nLocalElems) = nodeCoords(1:3,1:2,1:2,1:2,e1:e2) this%sideInfo(1:5,1:6,1:nLocalElems) = sideInfo(1:5,1:6,e1:e2) ! The mortar table is replicated on all ranks; element ids are global this%nMortars = 1 allocate(this%mortarInfo(1:14,1:1)) this%mortarInfo(1:14,1) = minfo(1:14) call this%UpdateDevice() endsubroutine SimpleMortarMesh_Mesh3D_t subroutine SimpleAnisoMortarMesh_Mesh3D_t(this,dx,bcids,comm,flips,splitCoord) !! !! Create the smallest 3D anisotropic (2-subface) nonconforming mesh: one "big" !! element whose east face is shared with the faces of two "small" dx x dx x dx !! elements that split it along one face coordinate only. The mesh is conforming !! everywhere except at the single mortar interface. !! !! With splitCoord = 1 the big element is [0,2dx] x [0,2dx] x [0,dx] and the small !! elements sit at y in [0,dx] and [dx,2dx]; the big east face's first trace !! coordinate (xi2, along y) is split. With splitCoord = 2 the big element is !! [0,2dx] x [0,dx] x [0,2dx] and the small elements are stacked in z, splitting the !! second trace coordinate (xi3). The small elements occupy x in [2dx,3dx]. !! !! y = 2dx ________________ ______ !! | | e3 | (splitCoord = 1, view along z; e2 and e3 !! | e1 |______| sit on the halves h = 1 and 2) !! | | e2 | !! |________________|______| !! x = 0 2dx 3dx !! !! Input !! - this : Fresh/empty Mesh3D_t object !! - dx : Edge length of the small elements (> 0) !! - bcids(1:6) : Boundary condition flags for the bottom, south, east, north, !! west, and top sides of the domain (the side ordering of !! StructuredMesh) !! - flips(1:2) (optional) : requested face flip (0..7) for the small element on !! each half; each small element is rigidly rotated about its !! centre as in SimpleMortarMesh. Defaults to 0. !! - splitCoord (optional) : 1 or 2, the big-face trace coordinate that is split. !! Defaults to 1. !! !! Output !! - this : Mesh3D_t object with three elements and one anisotropic mortar !! (nMortars = 0, nAnisoMortars = 1) !! !! As in SimpleMortarMesh, all connectivity is derived by matching sampled points of !! the bilinear face maps, so sideInfo and anisoMortarInfo are correct by !! construction for every rotation. Element geometry is trilinear (nGeo=1). On two !! ranks the mortar interface straddles the rank boundary. !! implicit none class(Mesh3D_t),intent(out) :: this real(prec),intent(in) :: dx integer,intent(in) :: bcids(1:6) integer,intent(in),optional :: comm integer,intent(in),optional :: flips(1:2) integer,intent(in),optional :: splitCoord ! Local integer,parameter :: nGlobalElem = 3 ! Signed local axis directions of the rotated small elements; see SimpleMortarMesh. integer,parameter :: flipAxes(1:3,0:7) = reshape([ & 1,2,3, & ! flip 0 : West donor -1,-2,3, & ! flip 1 : East donor 1,-2,-3, & ! flip 2 : West donor -1,2,-3, & ! flip 3 : East donor -1,3,2, & ! flip 4 : East donor 1,-3,2, & ! flip 5 : West donor -1,-3,-2, & ! flip 6 : East donor 1,3,-2],[3,8]) ! flip 7 : West donor integer,parameter :: invFlip(0:7) = [0,1,2,3,4,7,6,5] real(prec) :: nodeCoords(1:3,1:2,1:2,1:2,1:nGlobalElem) integer :: sideInfo(1:5,1:6,1:nGlobalElem) integer :: minfo(1:9) integer :: localFlips(1:2),split real(prec) :: center(1:3),axis(1:3,1:3),bigExtent(1:3),domainExtent(1:3) real(prec) :: tol integer :: h,f,i,j,k,d integer :: eA,eB,sA,sB,fAB integer :: e1,e2 integer :: gid,nUniqueSides integer :: nLocalElems integer :: nGeo,nBCs logical :: matched call this%decomp%init(comm) nGeo = 1 ! Trilinear element geometry nBCs = 6 tol = 1.0e-3_prec*dx localFlips = 0 if(present(flips)) then localFlips = flips do h = 1,2 if(localFlips(h) < 0 .or. localFlips(h) > 7) then print*,"SimpleAnisoMortarMesh_Mesh3D_t : requested flip out of range 0..7" stop 1 endif enddo endif split = 1 if(present(splitCoord)) split = splitCoord if(split /= 1 .and. split /= 2) then print*,"SimpleAnisoMortarMesh_Mesh3D_t : splitCoord must be 1 or 2" stop 1 endif ! Element 1 (big) : 2dx along x and along the split direction, dx along the other if(split == 1) then bigExtent = [2.0_prec*dx,2.0_prec*dx,dx] else bigExtent = [2.0_prec*dx,dx,2.0_prec*dx] endif domainExtent = bigExtent domainExtent(1) = 3.0_prec*dx do k = 1,2 do j = 1,2 do i = 1,2 nodeCoords(1:3,i,j,k,1) = [real(i-1,prec),real(j-1,prec),real(k-1,prec)]*bigExtent enddo enddo enddo ! Elements 2,3 (small) : [2dx,3dx] x half h of the big east face, each rigidly ! rotated per the requested flip do h = 1,2 f = localFlips(h) center(1) = 2.5_prec*dx center(2) = 0.5_prec*dx center(3) = 0.5_prec*dx center(1+split) = (real(h,prec)-0.5_prec)*dx do d = 1,3 axis(1:3,d) = 0.0_prec axis(abs(flipAxes(d,f)),d) = real(sign(1,flipAxes(d,f)),prec) enddo do k = 1,2 do j = 1,2 do i = 1,2 nodeCoords(1:3,i,j,k,1+h) = center+ & 0.5_prec*dx*(real(2*i-3,prec)*axis(1:3,1)+ & real(2*j-3,prec)*axis(1:3,2)+ & real(2*k-3,prec)*axis(1:3,3)) enddo enddo enddo enddo sideInfo = 0 ! Conforming interior faces do eA = 1,nGlobalElem do eB = eA+1,nGlobalElem do sA = 1,6 if(sideInfo(3,sA,eA) /= 0) cycle do sB = 1,6 call MatchFaces3D(nodeCoords(:,:,:,:,eA),sA, & nodeCoords(:,:,:,:,eB),sB,tol,matched,fAB) if(matched) then sideInfo(3,sA,eA) = eB sideInfo(4,sA,eA) = 10*sB+fAB sideInfo(3,sB,eB) = eA sideInfo(4,sB,eB) = 10*sA+invFlip(fAB) exit endif enddo enddo enddo enddo ! Mortar interface: the big element's east face against each small element minfo = 0 minfo(1) = 1 minfo(2) = selfSide3D_East minfo(9) = split do h = 1,2 matched = .false. do sB = 1,6 do f = 0,7 if(HalfMatches3D(nodeCoords(:,:,:,:,1),selfSide3D_East,h,split, & nodeCoords(:,:,:,:,1+h),sB,f,tol)) then minfo(2*h+1) = 1+h minfo(2*h+2) = 10*sB+f matched = .true. exit endif enddo if(matched) exit enddo if(.not. matched) then print*,"SimpleAnisoMortarMesh_Mesh3D_t : failed to match mortar half",h stop 1 endif enddo ! Physical boundary conditions on the faces that are neither conforming nor mortar do eA = 1,nGlobalElem do sA = 1,6 if(sideInfo(3,sA,eA) /= 0) cycle if(eA == 1 .and. sA == selfSide3D_East) cycle matched = .false. do h = 1,2 if(eA == minfo(2*h+1) .and. sA == minfo(2*h+2)/10) matched = .true. enddo if(matched) cycle ! small mortar face sideInfo(5,sA,eA) = DomainBoxBoundaryId3D(nodeCoords(:,:,:,:,eA),sA,domainExtent,tol,bcids) enddo enddo ! Global side ids: conforming pairs share one id; the big mortar face shares the ! id of half 1; the two halves get their own ids (MPI message tags). gid = 0 do eA = 1,nGlobalElem do sA = 1,6 if(sideInfo(2,sA,eA) /= 0) cycle if(eA == 1 .and. sA == selfSide3D_East) cycle matched = .false. do h = 1,2 if(eA == minfo(2*h+1) .and. sA == minfo(2*h+2)/10) matched = .true. enddo if(matched) cycle gid = gid+1 sideInfo(2,sA,eA) = gid eB = sideInfo(3,sA,eA) if(eB /= 0) then sB = sideInfo(4,sA,eA)/10 sideInfo(2,sB,eB) = gid endif enddo enddo do h = 1,2 gid = gid+1 minfo(6+h) = gid sideInfo(2,minfo(2*h+2)/10,minfo(2*h+1)) = gid sideInfo(1,minfo(2*h+2)/10,minfo(2*h+1)) = 1 ! nMortars + anisotropic mortar index if(h == 1) then sideInfo(2,selfSide3D_East,1) = gid sideInfo(1,selfSide3D_East,1) = 1 endif enddo nUniqueSides = gid call this%decomp%GenerateDecomposition(nGlobalElem,64*nUniqueSides) e1 = this%decomp%offsetElem(this%decomp%rankId+1)+1 e2 = this%decomp%offsetElem(this%decomp%rankId+2) nLocalElems = e2-e1+1 call this%Init(nGeo,nLocalElems,nLocalElems*6,nLocalElems*8,nBCs) this%nUniqueSides = nUniqueSides this%quadrature = UNIFORM this%BCType = 0 this%elemInfo = 0 this%globalNodeIDs = 0 this%nodeCoords(1:3,1:2,1:2,1:2,1:nLocalElems) = nodeCoords(1:3,1:2,1:2,1:2,e1:e2) this%sideInfo(1:5,1:6,1:nLocalElems) = sideInfo(1:5,1:6,e1:e2) ! The mortar table is replicated on all ranks; element ids are global this%nAnisoMortars = 1 allocate(this%anisoMortarInfo(1:9,1:1)) this%anisoMortarInfo(1:9,1) = minfo(1:9) call this%UpdateDevice() endsubroutine SimpleAnisoMortarMesh_Mesh3D_t pure function BilinearFacePoint3D(elemCoords,s,u,v) result(p) !! Evaluates the bilinear map of face s of a trilinear (nGeo=1) element at the !! face coordinates (u,v) in [-1,1]^2, using the face trace-coordinate convention !! of BoundaryInterp: faces 1,6 -> (xi1,xi2); faces 2,4 -> (xi1,xi3); !! faces 3,5 -> (xi2,xi3). implicit none real(prec),intent(in) :: elemCoords(1:3,1:2,1:2,1:2) integer,intent(in) :: s real(prec),intent(in) :: u,v real(prec) :: p(1:3) ! Local real(prec) :: c(1:3,1:2,1:2) real(prec) :: wu(1:2),wv(1:2) integer :: i,j select case(s) case(selfSide3D_Bottom) c = elemCoords(1:3,1:2,1:2,1) case(selfSide3D_South) c = elemCoords(1:3,1:2,1,1:2) case(selfSide3D_East) c = elemCoords(1:3,2,1:2,1:2) case(selfSide3D_North) c = elemCoords(1:3,1:2,2,1:2) case(selfSide3D_West) c = elemCoords(1:3,1,1:2,1:2) case default ! selfSide3D_Top c = elemCoords(1:3,1:2,1:2,2) endselect wu = [0.5_prec*(1.0_prec-u),0.5_prec*(1.0_prec+u)] wv = [0.5_prec*(1.0_prec-v),0.5_prec*(1.0_prec+v)] p = 0.0_prec do j = 1,2 do i = 1,2 p = p+wu(i)*wv(j)*c(1:3,i,j) enddo enddo endfunction BilinearFacePoint3D pure subroutine FlipFaceCoords3D(u,v,flip,u2,v2) !! Coordinate form of MortarFaceMap: maps receiver-face coordinates (u,v) to !! donor-face coordinates (u2,v2) for each of the eight SELF face flips. implicit none real(prec),intent(in) :: u,v integer,intent(in) :: flip real(prec),intent(out) :: u2,v2 select case(flip) case(0) u2 = u v2 = v case(1) u2 = -u v2 = v case(2) u2 = -u v2 = -v case(3) u2 = u v2 = -v case(4) u2 = v v2 = u case(5) u2 = -v v2 = u case(6) u2 = -v v2 = -u case default ! flip == 7 u2 = v v2 = -u endselect endsubroutine FlipFaceCoords3D subroutine MatchFaces3D(coordsA,sA,coordsB,sB,tol,matched,flip) !! Determines whether face sA of element A coincides with face sB of element B, !! and if so with which flip: face_B(F_flip(u,v)) == face_A(u,v) at the sampled !! points. Bilinear faces are pinned by their four corners plus one interior !! point. implicit none real(prec),intent(in) :: coordsA(1:3,1:2,1:2,1:2) integer,intent(in) :: sA real(prec),intent(in) :: coordsB(1:3,1:2,1:2,1:2) integer,intent(in) :: sB real(prec),intent(in) :: tol logical,intent(out) :: matched integer,intent(out) :: flip ! Local real(prec),parameter :: us(1:5) = [-1.0_prec,1.0_prec,-1.0_prec,1.0_prec,0.25_prec] real(prec),parameter :: vs(1:5) = [-1.0_prec,-1.0_prec,1.0_prec,1.0_prec,-0.5_prec] real(prec) :: pA(1:3),pB(1:3) real(prec) :: u2,v2 integer :: f,n logical :: ok matched = .false. flip = 0 do f = 0,7 ok = .true. do n = 1,5 pA = BilinearFacePoint3D(coordsA,sA,us(n),vs(n)) call FlipFaceCoords3D(us(n),vs(n),f,u2,v2) pB = BilinearFacePoint3D(coordsB,sB,u2,v2) if(maxval(abs(pA-pB)) > tol) then ok = .false. exit endif enddo if(ok) then matched = .true. flip = f return endif enddo endsubroutine MatchFaces3D function QuadrantMatches3D(coordsBig,sBig,q,coordsSmall,sSmall,flip,tol) result(matched) !! Determines whether face sSmall of the small element coincides with quadrant q !! of face sBig of the big element under the given flip: !! face_small(F_flip(u,v)) == face_big(quadrant_q(u,v)) at the sampled points. implicit none real(prec),intent(in) :: coordsBig(1:3,1:2,1:2,1:2) integer,intent(in) :: sBig,q real(prec),intent(in) :: coordsSmall(1:3,1:2,1:2,1:2) integer,intent(in) :: sSmall,flip real(prec),intent(in) :: tol logical :: matched ! Local real(prec),parameter :: us(1:5) = [-1.0_prec,1.0_prec,-1.0_prec,1.0_prec,0.25_prec] real(prec),parameter :: vs(1:5) = [-1.0_prec,-1.0_prec,1.0_prec,1.0_prec,-0.5_prec] real(prec) :: pBig(1:3),pSmall(1:3) real(prec) :: uq,vq,u2,v2 integer :: n matched = .true. do n = 1,5 ! Big-face coordinates of the quadrant point uq = 0.5_prec*(us(n)+real(2*mortarQuadKx(q)-3,prec)) vq = 0.5_prec*(vs(n)+real(2*mortarQuadKy(q)-3,prec)) pBig = BilinearFacePoint3D(coordsBig,sBig,uq,vq) call FlipFaceCoords3D(us(n),vs(n),flip,u2,v2) pSmall = BilinearFacePoint3D(coordsSmall,sSmall,u2,v2) if(maxval(abs(pBig-pSmall)) > tol) then matched = .false. return endif enddo endfunction QuadrantMatches3D function HalfMatches3D(coordsBig,sBig,h,split,coordsSmall,sSmall,flip,tol) result(matched) !! Determines whether face sSmall of the small element coincides with half h of !! face sBig of the big element, split along big-face trace coordinate split (1 or !! 2), under the given flip: face_small(F_flip(u,v)) == face_big(half_h(u,v)) at !! the sampled points. implicit none real(prec),intent(in) :: coordsBig(1:3,1:2,1:2,1:2) integer,intent(in) :: sBig,h,split real(prec),intent(in) :: coordsSmall(1:3,1:2,1:2,1:2) integer,intent(in) :: sSmall,flip real(prec),intent(in) :: tol logical :: matched ! Local real(prec),parameter :: us(1:5) = [-1.0_prec,1.0_prec,-1.0_prec,1.0_prec,0.25_prec] real(prec),parameter :: vs(1:5) = [-1.0_prec,-1.0_prec,1.0_prec,1.0_prec,-0.5_prec] real(prec) :: pBig(1:3),pSmall(1:3) real(prec) :: uh,vh,u2,v2 integer :: n matched = .true. do n = 1,5 ! Big-face coordinates of the half's point uh = us(n) vh = vs(n) if(split == 1) then uh = 0.5_prec*(us(n)+real(2*h-3,prec)) else vh = 0.5_prec*(vs(n)+real(2*h-3,prec)) endif pBig = BilinearFacePoint3D(coordsBig,sBig,uh,vh) call FlipFaceCoords3D(us(n),vs(n),flip,u2,v2) pSmall = BilinearFacePoint3D(coordsSmall,sSmall,u2,v2) if(maxval(abs(pBig-pSmall)) > tol) then matched = .false. return endif enddo endfunction HalfMatches3D function DomainBoxBoundaryId3D(elemCoords,s,extent,tol,bcids) result(bcid) !! Identifies which of the six planes of the box [0,extent(1)] x [0,extent(2)] x !! [0,extent(3)] contains face s of the element, and returns the corresponding !! boundary condition id (StructuredMesh side ordering). implicit none real(prec),intent(in) :: elemCoords(1:3,1:2,1:2,1:2) integer,intent(in) :: s real(prec),intent(in) :: extent(1:3),tol integer,intent(in) :: bcids(1:6) integer :: bcid ! Local real(prec) :: c(1:3) c = BilinearFacePoint3D(elemCoords,s,0.0_prec,0.0_prec) if(abs(c(3)) < tol) then bcid = bcids(1) ! bottom, z = 0 elseif(abs(c(2)) < tol) then bcid = bcids(2) ! south, y = 0 elseif(abs(c(1)-extent(1)) < tol) then bcid = bcids(3) ! east elseif(abs(c(2)-extent(2)) < tol) then bcid = bcids(4) ! north elseif(abs(c(1)) < tol) then bcid = bcids(5) ! west, x = 0 elseif(abs(c(3)-extent(3)) < tol) then bcid = bcids(6) ! top else print*,"DomainBoxBoundaryId3D : face centroid is not on a domain boundary plane" stop 1 endif endfunction DomainBoxBoundaryId3D function DomainBoundaryId3D(elemCoords,s,dx,tol,bcids) result(bcid) !! Identifies which of the six domain boundary planes of the SimpleMortarMesh !! contains face s of the element, and returns the corresponding boundary !! condition id. The domain is [0,3dx] x [0,2dx] x [0,2dx] with the region !! x > 2dx, above/below the small elements, outside the mesh; the x = 3dx plane !! is the domain east. implicit none real(prec),intent(in) :: elemCoords(1:3,1:2,1:2,1:2) integer,intent(in) :: s real(prec),intent(in) :: dx,tol integer,intent(in) :: bcids(1:6) integer :: bcid ! Local real(prec) :: c(1:3) c = BilinearFacePoint3D(elemCoords,s,0.0_prec,0.0_prec) if(abs(c(3)) < tol) then bcid = bcids(1) ! bottom, z = 0 elseif(abs(c(2)) < tol) then bcid = bcids(2) ! south, y = 0 elseif(abs(c(1)-3.0_prec*dx) < tol) then bcid = bcids(3) ! east, x = 3dx elseif(abs(c(2)-2.0_prec*dx) < tol) then bcid = bcids(4) ! north, y = 2dx elseif(abs(c(1)) < tol) then bcid = bcids(5) ! west, x = 0 elseif(abs(c(3)-2.0_prec*dx) < tol) then bcid = bcids(6) ! top, z = 2dx else print*,"DomainBoundaryId3D : face centroid is not on a domain boundary plane" stop 1 endif endfunction DomainBoundaryId3D subroutine UniformPeriodicMesh_Mesh3D_t(this,nxPerTile,nyPerTile,nzPerTile, & nTileX,nTileY,nTileZ,dx,dy,dz,comm) !! !! 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 integer,intent(in),optional :: comm ! 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(comm) 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 subroutine GeometryNodes1D(nGeo,quadrature,nodes) !! Returns the nGeo+1 one-dimensional geometry node positions, in [-1,1], of a mesh !! whose node coordinates are sampled on `quadrature` points. These are the points !! GenerateFromMesh interpolates the node coordinates from (see Lagrange_t%Init), so !! a coordinate sampled here is reproduced exactly by the geometry interpolant. implicit none integer,intent(in) :: nGeo integer,intent(in) :: quadrature real(prec),intent(out) :: nodes(1:nGeo+1) ! Local real(prec) :: weights(1:nGeo+1) if(quadrature == GAUSS .or. quadrature == GAUSS_LOBATTO) then call LegendreQuadrature(nGeo,nodes,weights,quadrature) elseif(quadrature == CHEBYSHEV_GAUSS .or. quadrature == CHEBYSHEV_GAUSS_LOBATTO) then call ChebyshevQuadrature(nGeo,nodes,weights,quadrature) elseif(quadrature == UNIFORM) then nodes = UniformPoints(-1.0_prec,1.0_prec,0,nGeo) else error stop "GeometryNodes1D : unrecognized mesh quadrature type" endif endsubroutine GeometryNodes1D subroutine ExtrudedMesh_Mesh3D_t(this,mesh2d,zLevels,bottomBC,topBC) !! Builds a layered hexahedral mesh by extruding a surface mesh downward through !! nz = size(zLevels,3)-1 layers bounded by the level surfaces z_0 < z_1 < ... < z_nz. !! !! Element k (k = 1..nz) of the column above surface element Q is the map !! !! x(xi,eta,zeta) = ( X(xi,eta), Y(xi,eta), Z(xi,eta,zeta) ) !! !! where (X,Y) is Q's isoparametric map of degree nGeo and, at every geometry node !! (xi_i,eta_j), !! !! Z(xi_i,eta_j,zeta) = (1-zeta)/2 z_{k-1}(xi_i,eta_j) + (1+zeta)/2 z_k(xi_i,eta_j). !! !! The horizontal coordinates do not depend on zeta, so lines of constant (xi,eta) !! are vertical, and Z is linear in zeta within a layer. The 3-D mesh keeps the 2-D !! nGeo and quadrature: its node coordinates are sampled on the same nGeo+1 points !! in all three directions, so they form the usual (1:3,1:nGeo+1,1:nGeo+1,1:nGeo+1, !! 1:nElem) array and the linear vertical map is represented exactly. The two !! endpoint weights are exact in floating point, so the top face of layer k and the !! bottom face of layer k+1 are sampled from the same z_k and agree bitwise. !! !! Input !! - this : Fresh/empty Mesh3D_t object. The result is an ordinary Mesh3D, so it !! works with SEMHex, the mapped data types and the DG models unchanged. !! - mesh2d : The surface mesh, already built and decomposed (any Mesh2D: structured, !! periodic structured, HOPr, HOHQMesh or an adapted quadtree mesh, with !! nGeo >= 1). It may carry 2:1 mortar edges (nMortars > 0); see "Mortars" below. !! - zLevels(1:nGeo+1,1:nGeo+1,0:nz,1:mesh2d%nElem) : Level heights z_k at the !! geometry nodes of each rank-local surface element, in the length unit of the !! mesh coordinates, positive upward. Required: nz >= 1, z_k strictly increasing !! in k at every node, and z_nz = 0 exactly (a flat top at the reference sea !! surface). The levels must also be continuous across lateral faces, i.e. two !! elements must give the same z_k at a shared node. That holds whenever the !! levels are computed pointwise from a continuous bathymetry sampled at the !! geometry nodes, as ExtrudedSigmaLevels and ExtrudedHybridLevels do; it is not !! checked here, because the neighbor of a lateral face may be on another rank. !! - bottomBC : Boundary condition id of the bottom face of every column (z = z_0) !! - topBC : Boundary condition id of the top face of every column (z = 0) !! !! Output !! - this : Mesh3D_t object with nLayers = nz and nz*mesh2d%nElem local elements !! !! Ordering and decomposition. Local element e3 = (e2-1)*nz + k is layer k of local !! surface element e2, and global ids follow the same rule. The 3-D decomposition is !! the 2-D one expanded by nz (offsetElem3D = nz*offsetElem2D), so a column is never !! split across ranks and the local surface element e2 of an extruded element is !! element e2 of the caller's mesh2d. The surface mesh is not copied; keep it alive if !! column operators need its geometry. !! !! Connectivity. The 2-D sides 1..4 (south, east, north, west) become the lateral !! faces 2..5 of the 3-D element, and faces 1 and 6 are the bottom and top. A lateral !! face inherits the neighbor, neighbor side and boundary id of its 2-D edge, in the !! same layer. Its trace coordinates are (edge coordinate, zeta) on both sides !! (BoundaryInterp: faces 2,4 -> (xi1,xi3), faces 3,5 -> (xi2,xi3)) and zeta runs !! upward in both elements, so the SELF face flip equals the 2-D edge flip: 0 maps to !! face flip 0 (i2 = i, j2 = j) and 1 maps to face flip 1 (i2 = N-i, j2 = j). Bottom !! and top faces between layers are conforming with flip 0. !! !! Global side ids are unique across ranks and serve as MPI message tags. The lateral !! face of layer k on the 2-D edge with id g has id sign(g)*((|g|-1)*nz + k), so the !! HOPr master/slave sign of the 2-D id is preserved. The horizontal faces of the !! column above global surface element E are numbered after all lateral ones, !! nz*nUniqueSides2D + (E-1)*(nz+1) + l for level l-1 (l = 1..nz+1), and !! !! nUniqueSides = nz*nUniqueSides2D + nGlobalElem2D*(nz+1). !! !! Mortars. Every 2-D mortar edge m becomes nz anisotropic (2-subface) mortars, one !! per layer: anisotropic mortar a = (m-1)*nz + k joins the big lateral face of layer k !! to the two small lateral faces of layer k. The face's first trace coordinate is the !! edge coordinate, so it is the split coordinate (anisoMortarInfo(9,a) = 1), the 2-D !! sub-edge h becomes half h, and the 2-D edge flip becomes the face flip as for a !! conforming lateral face. The sub-face global ids follow the lateral rule above !! applied to the 2-D sub-edge ids, and the faces carry sideInfo(1) = a. The result has !! nMortars = 0 and nAnisoMortars = nz*mesh2d%nMortars. !! !! The node ids (globalNodeIDs) and elemInfo are zeroed: the solver does not use !! them, and the flips are set directly rather than through RecalculateFlip. implicit none class(Mesh3D_t),intent(out) :: this class(Mesh2D_t),intent(in) :: mesh2d real(prec),intent(in) :: zLevels(1:,1:,0:,1:) integer,intent(in) :: bottomBC integer,intent(in) :: topBC ! Local integer :: nz,nGeo,nElem2D,nElem3D,nRanks integer :: nGlobalElem2D,nUniqueSides2D,nUniqueSides3D,nHorizontalBase integer :: e2,e3,k,s,i,j,nk integer :: E2global,E2neighbor,g,s2,flip,m,a,h integer,allocatable :: offsetElem3D(:) real(prec) :: zeta(1:mesh2d%nGeo+1) nGeo = mesh2d%nGeo nElem2D = mesh2d%nElem nz = size(zLevels,3)-1 ! ---- Validate the surface mesh and the levels ---- if(.not. mesh2d%decomp%initialized) then error stop "ExtrudedMesh : the surface mesh has no domain decomposition; build it first" endif if(nz < 1) then error stop "ExtrudedMesh : zLevels must carry at least two levels (nz >= 1)" endif if(size(zLevels,1) /= nGeo+1 .or. size(zLevels,2) /= nGeo+1) then error stop "ExtrudedMesh : zLevels must be dimensioned (1:nGeo+1,1:nGeo+1,0:nz,1:nElem2D)" endif if(size(zLevels,4) /= nElem2D) then error stop "ExtrudedMesh : the last dimension of zLevels must equal the local surface element count" endif do e2 = 1,nElem2D do j = 1,nGeo+1 do i = 1,nGeo+1 if(zLevels(i,j,nz,e2) /= 0.0_prec) then error stop "ExtrudedMesh : the top level zLevels(:,:,nz,:) must be exactly 0" endif do k = 1,nz if(.not.(zLevels(i,j,k,e2) > zLevels(i,j,k-1,e2))) then error stop "ExtrudedMesh : zLevels must be strictly increasing in the level index" endif enddo enddo enddo enddo nRanks = mesh2d%decomp%nRanks nGlobalElem2D = mesh2d%decomp%offsetElem(nRanks+1) if(nElem2D /= mesh2d%decomp%offsetElem(mesh2d%decomp%rankId+2)- & mesh2d%decomp%offsetElem(mesh2d%decomp%rankId+1)) then ! The columns inherit the surface partition, so the surface mesh must hold ! exactly the elements its decomposition assigns to this rank. error stop "ExtrudedMesh : the surface mesh does not hold the elements its decomposition assigns to this rank" endif nUniqueSides2D = mesh2d%nUniqueSides do e2 = 1,nElem2D do s = 1,4 g = mesh2d%sideInfo(2,s,e2) if(g == 0 .or. abs(g) > nUniqueSides2D) then error stop "ExtrudedMesh : every surface edge needs a global side id in 1..nUniqueSides" endif if(mesh2d%sideInfo(3,s,e2) > 0) then flip = mesh2d%sideInfo(4,s,e2)-10*(mesh2d%sideInfo(4,s,e2)/10) if(flip /= 0 .and. flip /= 1) then error stop "ExtrudedMesh : a surface edge flip must be 0 or 1" endif endif enddo enddo nElem3D = nz*nElem2D nHorizontalBase = nz*nUniqueSides2D nUniqueSides3D = nHorizontalBase+nGlobalElem2D*(nz+1) ! ---- Decomposition: the surface partition, expanded by nz ---- call this%decomp%Init(comm=mesh2d%decomp%mpiComm) allocate(offsetElem3D(1:nRanks+1)) offsetElem3D(1:nRanks+1) = nz*mesh2d%decomp%offsetElem(1:nRanks+1) call this%decomp%GenerateDecompositionFromOffsets(offsetElem3D,nUniqueSides3D) deallocate(offsetElem3D) call this%Init(nGeo,nElem3D,6*nElem3D,nElem3D*(nGeo+1)**3,mesh2d%nBCs) this%nGlobalElem = nz*nGlobalElem2D this%nUniqueSides = nUniqueSides3D ! GLOBAL side id space on every rank (the MPI tag stride) this%quadrature = mesh2d%quadrature this%nLayers = nz if(mesh2d%nBCs > 0) then this%BCType(1:4,1:mesh2d%nBCs) = mesh2d%BCType(1:4,1:mesh2d%nBCs) this%BCNames(1:mesh2d%nBCs) = mesh2d%BCNames(1:mesh2d%nBCs) endif if(allocated(mesh2d%materialNames) .and. allocated(mesh2d%elemMaterial)) then deallocate(this%materialNames) allocate(this%materialNames(1:mesh2d%nMaterials)) this%nMaterials = mesh2d%nMaterials this%materialNames(1:mesh2d%nMaterials) = mesh2d%materialNames(1:mesh2d%nMaterials) do e2 = 1,nElem2D do k = 1,nz this%elemMaterial((e2-1)*nz+k) = mesh2d%elemMaterial(e2) enddo enddo endif this%globalNodeIDs = 0 ! node ids are unused by the solver (flips are set directly) this%elemInfo = 0 ! ---- Geometry: vertical lines of constant (xi,eta), linear in zeta per layer ---- call GeometryNodes1D(nGeo,mesh2d%quadrature,zeta) if(mesh2d%quadrature == GAUSS_LOBATTO .or. mesh2d%quadrature == CHEBYSHEV_GAUSS_LOBATTO .or. & mesh2d%quadrature == UNIFORM) then ! The end nodes lie on the faces. Pin them to -1 and 1 so that the end weights ! below are exactly 0 and 1 and each face is sampled from one level alone. zeta(1) = -1.0_prec zeta(nGeo+1) = 1.0_prec endif do e2 = 1,nElem2D do k = 1,nz e3 = (e2-1)*nz+k do nk = 1,nGeo+1 do j = 1,nGeo+1 do i = 1,nGeo+1 this%nodeCoords(1,i,j,nk,e3) = mesh2d%nodeCoords(1,i,j,e2) this%nodeCoords(2,i,j,nk,e3) = mesh2d%nodeCoords(2,i,j,e2) this%nodeCoords(3,i,j,nk,e3) = 0.5_prec*(1.0_prec-zeta(nk))*zLevels(i,j,k-1,e2)+ & 0.5_prec*(1.0_prec+zeta(nk))*zLevels(i,j,k,e2) enddo enddo enddo enddo enddo ! ---- Connectivity ---- ! sideInfo(1:5,iSide,iEl) ! 1 - Side Type (0, or the anisotropic mortar index on a mortar face) ! 2 - Global Side ID ! 3 - Neighbor Element ID (global; 0 on a domain boundary) ! 4 - 10*( neighbor local side ) + flip ! 5 - Boundary Condition ID (0 on an interior face) this%sideInfo = 0 do e2 = 1,nElem2D E2global = mesh2d%decomp%offsetElem(mesh2d%decomp%rankId+1)+e2 do k = 1,nz e3 = (e2-1)*nz+k ! Lateral faces: 2-D side s is 3-D face s+1 (south, east, north, west) do s = 1,4 g = mesh2d%sideInfo(2,s,e2) this%sideInfo(2,s+1,e3) = sign(1,g)*((abs(g)-1)*nz+k) E2neighbor = mesh2d%sideInfo(3,s,e2) if(E2neighbor > 0) then s2 = mesh2d%sideInfo(4,s,e2)/10 flip = mesh2d%sideInfo(4,s,e2)-10*s2 this%sideInfo(3,s+1,e3) = (E2neighbor-1)*nz+k this%sideInfo(4,s+1,e3) = 10*(s2+1)+flip endif this%sideInfo(5,s+1,e3) = mesh2d%sideInfo(5,s,e2) ! On a mesh with mortars, sideInfo(1) of a mortar edge is its 2-D mortar index if(mesh2d%nMortars > 0 .and. mesh2d%sideInfo(1,s,e2) /= 0) then this%sideInfo(1,s+1,e3) = (mesh2d%sideInfo(1,s,e2)-1)*nz+k endif enddo ! Bottom face, on level k-1 of the column this%sideInfo(2,selfSide3D_Bottom,e3) = nHorizontalBase+(E2global-1)*(nz+1)+k if(k == 1) then this%sideInfo(5,selfSide3D_Bottom,e3) = bottomBC else this%sideInfo(3,selfSide3D_Bottom,e3) = (E2global-1)*nz+k-1 this%sideInfo(4,selfSide3D_Bottom,e3) = 10*selfSide3D_Top endif ! Top face, on level k of the column this%sideInfo(2,selfSide3D_Top,e3) = nHorizontalBase+(E2global-1)*(nz+1)+k+1 if(k == nz) then this%sideInfo(5,selfSide3D_Top,e3) = topBC else this%sideInfo(3,selfSide3D_Top,e3) = (E2global-1)*nz+k+1 this%sideInfo(4,selfSide3D_Top,e3) = 10*selfSide3D_Bottom endif enddo enddo ! ---- Anisotropic mortars: one per 2-D mortar edge and layer ---- ! mortarInfo of the surface mesh is replicated and holds global element ids, so the ! extruded table is too. if(mesh2d%nMortars > 0) then this%nAnisoMortars = nz*mesh2d%nMortars allocate(this%anisoMortarInfo(1:9,1:this%nAnisoMortars)) do m = 1,mesh2d%nMortars do k = 1,nz a = (m-1)*nz+k this%anisoMortarInfo(1,a) = (mesh2d%mortarInfo(1,m)-1)*nz+k this%anisoMortarInfo(2,a) = mesh2d%mortarInfo(2,m)+1 do h = 1,2 s2 = mesh2d%mortarInfo(2*h+2,m)/10 flip = mesh2d%mortarInfo(2*h+2,m)-10*s2 this%anisoMortarInfo(2*h+1,a) = (mesh2d%mortarInfo(2*h+1,m)-1)*nz+k this%anisoMortarInfo(2*h+2,a) = 10*(s2+1)+flip g = mesh2d%mortarInfo(6+h,m) this%anisoMortarInfo(6+h,a) = sign(1,g)*((abs(g)-1)*nz+k) enddo this%anisoMortarInfo(9,a) = 1 enddo enddo endif call this%UpdateDevice() endsubroutine ExtrudedMesh_Mesh3D_t subroutine ExtrudedHybridLevels(h,depths,sigma,zLevels) !! Fills the level heights for ExtrudedMesh with nzFlat = size(depths) flat !! interfaces near the surface over nzS = size(sigma)-1 terrain-following layers. !! With d_0 = 0 and d_F = depths(nzFlat) (d_F = 0 when nzFlat = 0), the nz = nzS+nzFlat !! levels at every geometry node are !! !! z_k = -( (1+sigma_k) d_F - sigma_k h ) , k = 0..nzS !! z_{nzS+m} = -d_{nzFlat-m} , m = 1..nzFlat !! !! so z_0 = -h, z_nzS = -d_F, z_nz = 0, and the layers between z = -d_F and the !! bottom follow the terrain with the stretching sigma. The form used for z_k makes !! z_0 = -h and z_nzS = -d_F exact in floating point. With nzFlat = 0 this is the !! sigma coordinate z_k = sigma_k h (see ExtrudedSigmaLevels); with a constant h it !! gives flat z-levels throughout. !! !! The levels are evaluated pointwise from h, so they are continuous across lateral !! faces exactly when h is: two elements that share a geometry node must be given the !! same h there, as they are when h is sampled from one bathymetry function at the !! mesh's node coordinates. !! !! Input !! - h(1:nGeo+1,1:nGeo+1,1:nElem2D) : Water depth at the surface geometry nodes, in !! the length unit of the mesh, with h > d_F (> 0) everywhere !! - depths(1:nzFlat) : Depths of the flat interfaces below the surface, same unit, !! 0 < d_1 < ... < d_nzFlat < min(h). May be empty. !! - sigma(0:nzS) : Terrain-following stretching, nzS >= 1, sigma(0) = -1, !! sigma(nzS) = 0, strictly increasing !! !! Output !! - zLevels(1:nGeo+1,1:nGeo+1,0:nz,1:nElem2D) : Level heights, nz = nzS+nzFlat implicit none real(prec),intent(in) :: h(1:,1:,1:) real(prec),intent(in) :: depths(1:) real(prec),intent(in) :: sigma(0:) real(prec),intent(out) :: zLevels(1:,1:,0:,1:) ! Local integer :: nzFlat,nzS,nz,i,j,k,m,e real(prec) :: dF nzFlat = size(depths) nzS = size(sigma)-1 nz = nzS+nzFlat if(nzS < 1) then error stop "ExtrudedHybridLevels : sigma must carry at least two values (one terrain-following layer)" endif if(sigma(0) /= -1.0_prec .or. sigma(nzS) /= 0.0_prec) then error stop "ExtrudedHybridLevels : sigma must run from sigma(0) = -1 to sigma(nzS) = 0" endif do k = 1,nzS if(.not.(sigma(k) > sigma(k-1))) then error stop "ExtrudedHybridLevels : sigma must be strictly increasing" endif enddo if(nzFlat > 0) then if(.not.(depths(1) > 0.0_prec)) then error stop "ExtrudedHybridLevels : the flat interface depths must be positive" endif do m = 2,nzFlat if(.not.(depths(m) > depths(m-1))) then error stop "ExtrudedHybridLevels : the flat interface depths must be strictly increasing" endif enddo dF = depths(nzFlat) else dF = 0.0_prec endif if(size(zLevels,1) /= size(h,1) .or. size(zLevels,2) /= size(h,2) .or. & size(zLevels,3) /= nz+1 .or. size(zLevels,4) /= size(h,3)) then error stop "ExtrudedHybridLevels : zLevels must be dimensioned (1:nGeo+1,1:nGeo+1,0:nz,1:nElem2D)" endif if(size(h) > 0) then if(.not.(minval(h) > dF)) then error stop "ExtrudedHybridLevels : the depth h must exceed the deepest flat interface (and 0)" endif endif do e = 1,size(h,3) do j = 1,size(h,2) do i = 1,size(h,1) do k = 0,nzS zLevels(i,j,k,e) = -((1.0_prec+sigma(k))*dF-sigma(k)*h(i,j,e)) enddo do m = 1,nzFlat-1 zLevels(i,j,nzS+m,e) = -depths(nzFlat-m) enddo zLevels(i,j,nz,e) = 0.0_prec enddo enddo enddo endsubroutine ExtrudedHybridLevels subroutine ExtrudedSigmaLevels(h,sigma,zLevels) !! Fills the level heights for ExtrudedMesh with terrain-following (sigma) levels, !! !! z_k = sigma_k h , k = 0..nz , !! !! at every geometry node, so z_0 = -h and z_nz = 0. This is ExtrudedHybridLevels with !! no flat interfaces. The levels are continuous across lateral faces when h is (see !! ExtrudedHybridLevels). !! !! Input !! - h(1:nGeo+1,1:nGeo+1,1:nElem2D) : Water depth at the surface geometry nodes, in !! the length unit of the mesh, h > 0 !! - sigma(0:nz) : nz >= 1, sigma(0) = -1, sigma(nz) = 0, strictly increasing !! !! Output !! - zLevels(1:nGeo+1,1:nGeo+1,0:nz,1:nElem2D) : Level heights implicit none real(prec),intent(in) :: h(1:,1:,1:) real(prec),intent(in) :: sigma(0:) real(prec),intent(out) :: zLevels(1:,1:,0:,1:) ! Local real(prec) :: noDepths(1:0) call ExtrudedHybridLevels(h,noDepths,sigma,zLevels) endsubroutine ExtrudedSigmaLevels subroutine Read_HOPr_Mesh3D_t(this,meshFile,comm) ! From https://www.hopr-project.org/externals/Meshformat.pdf, Algorithm 6 implicit none class(Mesh3D_t),intent(out) :: this character(*),intent(in) :: meshFile integer,intent(in),optional :: comm ! Local integer(HID_T) :: fileId integer(HID_T) :: offset(1:2),gOffset(1) integer :: nGlobalElem integer :: firstElem integer :: firstNode integer :: firstSide integer :: nLocalElems integer :: nLocalNodes integer :: nLocalSides integer :: nUniqueSides integer :: nGeo,nBCs integer :: eid,lsid,iSide integer :: i,j,k,nid integer,dimension(:,:),allocatable :: hopr_elemInfo integer,dimension(:,:),allocatable :: hopr_sideInfo real(prec),dimension(:,:),allocatable :: hopr_nodeCoords integer,dimension(:),allocatable :: hopr_globalNodeIDs integer,dimension(:,:),allocatable :: bcType call this%decomp%init(comm) if(this%decomp%mpiEnabled) then call Open_HDF5(meshFile,H5F_ACC_RDONLY_F,fileId,this%decomp%mpiComm) else call Open_HDF5(meshFile,H5F_ACC_RDONLY_F,fileId) endif call ReadAttribute_HDF5(fileId,'nElems',nGlobalElem) call ReadAttribute_HDF5(fileId,'Ngeo',nGeo) call ReadAttribute_HDF5(fileId,'nBCs',nBCs) call ReadAttribute_HDF5(fileId,'nUniqueSides',nUniqueSides) ! Read BCType allocate(bcType(1:4,1:nBCs)) if(this%decomp%mpiEnabled) then offset(:) = 0 call ReadArray_HDF5(fileId,'BCType',bcType,offset) else call ReadArray_HDF5(fileId,'BCType',bcType) endif ! Read local subarray of ElemInfo call this%decomp%GenerateDecomposition(nGlobalElem,nUniqueSides) firstElem = this%decomp%offsetElem(this%decomp%rankId+1)+1 nLocalElems = this%decomp%offsetElem(this%decomp%rankId+2)- & this%decomp%offsetElem(this%decomp%rankId+1) ! Allocate Space for hopr_elemInfo! allocate(hopr_elemInfo(1:6,1:nLocalElems)) if(this%decomp%mpiEnabled) then offset = (/0,firstElem-1/) call ReadArray_HDF5(fileId,'ElemInfo',hopr_elemInfo,offset) else call ReadArray_HDF5(fileId,'ElemInfo',hopr_elemInfo) endif ! Read local subarray of NodeCoords and GlobalNodeIDs firstNode = hopr_elemInfo(5,1)+1 nLocalNodes = hopr_elemInfo(6,nLocalElems)-hopr_elemInfo(5,1) ! Allocate Space for hopr_nodeCoords and hopr_globalNodeIDs ! allocate(hopr_nodeCoords(1:3,1:nLocalNodes),hopr_globalNodeIDs(1:nLocalNodes)) if(this%decomp%mpiEnabled) then offset = (/0,firstNode-1/) call ReadArray_HDF5(fileId,'NodeCoords',hopr_nodeCoords,offset) gOffset = (/firstNode-1/) call ReadArray_HDF5(fileId,'GlobalNodeIDs',hopr_globalNodeIDs,gOffset) else call ReadArray_HDF5(fileId,'NodeCoords',hopr_nodeCoords) call ReadArray_HDF5(fileId,'GlobalNodeIDs',hopr_globalNodeIDs) endif ! Read local subarray of SideInfo firstSide = hopr_elemInfo(3,1)+1 nLocalSides = hopr_elemInfo(4,nLocalElems)-hopr_elemInfo(3,1) ! Allocate space for hopr_sideInfo allocate(hopr_sideInfo(1:5,1:nLocalSides)) if(this%decomp%mpiEnabled) then offset = (/0,firstSide-1/) call ReadArray_HDF5(fileId,'SideInfo',hopr_sideInfo,offset) else call ReadArray_HDF5(fileId,'SideInfo',hopr_sideInfo) endif call Close_HDF5(fileID) ! ---- Done reading 3-D Mesh information ---- ! ! Load hopr data into mesh data structure call this%Init(nGeo,nLocalElems,nLocalSides,nLocalNodes,nBCs) ! Copy data from local arrays into this this%elemInfo = hopr_elemInfo this%nUniqueSides = nUniqueSides this%quadrature = UNIFORM ! Grab the node coordinates do eid = 1,this%nElem do k = 1,nGeo+1 do j = 1,nGeo+1 do i = 1,nGeo+1 nid = i+(nGeo+1)*(j-1+(nGeo+1)*(k-1+(nGeo+1)*(eid-1))) this%nodeCoords(1:3,i,j,k,eid) = hopr_nodeCoords(1:3,nid) this%globalNodeIDs(i,j,k,eid) = hopr_globalNodeIDs(nid) enddo enddo enddo enddo iSide = 0 do eid = 1,this%nElem do lsid = 1,6 iSide = iSide+1 this%sideInfo(1:5,lsid,eid) = hopr_sideInfo(1:5,iSide) enddo enddo call this%RecalculateFlip() deallocate(hopr_elemInfo,hopr_nodeCoords,hopr_globalNodeIDs,hopr_sideInfo) call this%UpdateDevice() endsubroutine Read_HOPr_Mesh3D_t subroutine Read_HOHQMesh_Mesh3D_t(this,meshFile,comm) !! Reader for HOHQMesh 3-D (hexahedral) text mesh files in the ISM !! and ISM-MM formats, as written by HOHQMesh's `WriteISMHexMeshFile` !! (`Source/3DSource/Mesh3DOutputMethods.f90`). Unlike the 2-D !! writer, the 3-D writer emits NO format header line: the first !! line is always the count line "nNodes nElems polyOrder". The two !! variants differ only in the per-element corner-node line: !! * ISM : 8 corner-node ids !! * ISM-MM : 8 corner-node ids followed by a material-name string !! The variant is auto-detected from the presence of the 9th token. !! !! Each element block contains, in order: the corner-node line, a !! line of 6 boundary-face flags, a (polyOrder+1)^2 block of face !! points (x,y,z per line, inner index first) for every flagged !! face, and a line of 6 boundary-condition names ("---" marks an !! interior face). Face points are sampled at Chebyshev-Gauss- !! Lobatto points, so the resulting mesh has !! `quadrature = CHEBYSHEV_GAUSS_LOBATTO`. !! !! HOHQMesh numbers hex corners with nodes 1-4 on the bottom face !! (counter-clockwise) and 5-8 above them, which matches SELF's !! CGNS corner convention exactly. HOHQMesh face numbering !! (1=south, 2=north, 3=bottom, 4=east, 5=top, 6=west; see !! `FaceFromVolume`) is remapped to SELF's side ordering !! (1=bottom, 2=south, 3=east, 4=north, 5=west, 6=top). Face-point !! grids are written with the two on-face volume axes in natural !! order, which coincides with SELF's boundary index convention, !! so no reorientation of the face data is required. !! !! Element interior nodes are reconstructed by transfinite (Coons) !! interpolation of the six face grids; unflagged faces are !! bilinear patches of their corner nodes. Element-face !! connectivity is not present in the format, so neighbors are !! reconstructed by matching the sorted corner-node ids of each !! face across elements, and the side "flip" is computed by !! matching the corner-node orderings of the paired faces. !! Boundary names populate this%BCNames and sideInfo(5,...) carries !! the 1-based index into that table (0 for interior faces). !! Material names (ISM-MM) populate this%materialNames and !! this%elemMaterial; plain ISM meshes keep the single "default" !! material. implicit none class(Mesh3D_t),intent(out) :: this character(*),intent(in) :: meshFile integer,intent(in),optional :: comm ! Local integer :: iUnit integer :: ios integer :: nNodesFile,nElemFile,polyOrder integer :: nGeo,ng1 integer :: i,j,k,e,f,s,iSide,l,p integer :: cornerIDs(1:8) integer :: probe8(1:8) integer :: hohqFlag(1:6) integer :: quadA(1:4),quadB(1:4) integer :: sortedA(1:4) integer :: ePair,sPair,flip integer :: bcIdx integer :: matIdx integer :: nBCsLocal integer :: nMatsLocal integer :: nFaceRecords integer :: hashKey,bucket,probe integer :: hashSize integer,allocatable :: hashHead(:) integer,allocatable :: hashNext(:) integer,allocatable :: pairKey(:,:) ! sorted corner ids, 4 x 6*nElem integer,allocatable :: pairElem(:),pairSide(:) integer,allocatable :: ismCorners(:,:) ! 8 x nElem integer,allocatable :: ismBCid(:,:) ! 6 x nElem (SELF side ordering) integer,allocatable :: ismFlag(:,:) ! 6 x nElem (SELF side ordering) integer,allocatable :: ismMat(:) ! nElem real(prec),allocatable :: nodeXYZ(:,:) ! 3 x nNodesFile real(prec),allocatable :: faceCurve(:,:,:,:,:) ! 3, ng1, ng1, 6, nElem real(prec) :: xyz(1:3) character(LEN=512) :: lineBuf character(LEN=SELF_MESH_MATNAME_LENGTH) :: matName character(LEN=255) :: bdyNames(1:6) character(LEN=255),allocatable :: BCNamesLocal(:) character(LEN=SELF_MESH_MATNAME_LENGTH),allocatable :: matNamesLocal(:) logical :: isISM_MM ! Map from HOHQMesh hex face id (1=south, 2=north, 3=bottom, ! 4=east, 5=top, 6=west) to SELF local side id integer,parameter :: hohqToSelfSide(1:6) = [2,4,1,3,6,5] ! Corner permutations realised by SELF's eight face flips: for a ! side pair (s1,s2) with flip p, corner l of side s1 (in face-index ! order, see sideMap) coincides with corner flipPerm(l,p) of side ! s2. Flips 0-3 are rotations/reflections keeping the face axes, ! flips 4-7 transpose them (see ApplyFlip/SideExchange kernels). integer,parameter :: flipPerm(1:4,0:7) = reshape( & [1,2,3,4, & 2,1,4,3, & 3,4,1,2, & 4,3,2,1, & 1,4,3,2, & 2,3,4,1, & 3,2,1,4, & 4,1,2,3],[4,8]) call this%decomp%init(comm) open(newunit=iUnit,file=trim(meshFile),status='old',action='read', & form='formatted',iostat=ios) if(ios /= 0) then print*,__FILE__//' : Failed to open '//trim(meshFile) stop 1 endif print*,__FILE__//' : Reading HOHQMesh mesh from '//trim(meshFile) ! ---- 1. Count line (3-D ISM files carry no format header) ---- read(iUnit,*) nNodesFile,nElemFile,polyOrder nGeo = polyOrder ng1 = nGeo+1 print*,__FILE__//' : nNodes = ',nNodesFile,' nElem = ',nElemFile, & ' polyOrder = ',polyOrder ! ---- 2. Read all node coordinates ---- allocate(nodeXYZ(1:3,1:nNodesFile)) do i = 1,nNodesFile read(iUnit,*) xyz(1:3) nodeXYZ(1:3,i) = xyz(1:3) enddo ! ---- 3. Per-element block ---- allocate(ismCorners(1:8,1:nElemFile)) allocate(ismFlag(1:6,1:nElemFile)) allocate(ismBCid(1:6,1:nElemFile)) allocate(ismMat(1:nElemFile)) allocate(faceCurve(1:3,1:ng1,1:ng1,1:6,1:nElemFile)) faceCurve = 0.0_prec ! Boundary-name table built incrementally nBCsLocal = 0 allocate(BCNamesLocal(1:16)) BCNamesLocal = "" ! Material-name table built incrementally nMatsLocal = 0 allocate(matNamesLocal(1:8)) matNamesLocal = "" isISM_MM = .false. nFaceRecords = 0 do e = 1,nElemFile ! Corner-node line; the trailing material-name token (ISM-MM) ! is detected by probing for a 9th list item after the 8 ! integer tokens. read(iUnit,'(A)') lineBuf read(lineBuf,*) cornerIDs(1:8) read(lineBuf,*,iostat=ios) probe8,matName if(ios /= 0) matName = "" if(e == 1) then isISM_MM = (matName /= "") print*,__FILE__//' : Format = ',merge("ISM-MM","ISM ",isISM_MM) endif if(matName /= "") then ! lookup/insert material name matIdx = 0 do k = 1,nMatsLocal if(trim(matNamesLocal(k)) == trim(matName)) then matIdx = k; exit endif enddo if(matIdx == 0) then nMatsLocal = nMatsLocal+1 if(nMatsLocal > size(matNamesLocal)) call grow_string_table(matNamesLocal) matNamesLocal(nMatsLocal) = matName matIdx = nMatsLocal endif ismMat(e) = matIdx else ismMat(e) = 1 endif ismCorners(1:8,e) = cornerIDs ! Face flags and face-point grids, in HOHQMesh face order read(iUnit,*) hohqFlag(1:6) do f = 1,6 s = hohqToSelfSide(f) ismFlag(s,e) = hohqFlag(f) if(hohqFlag(f) == 1) then nFaceRecords = nFaceRecords+1 do j = 1,ng1 do i = 1,ng1 read(iUnit,*) xyz(1:3) faceCurve(1:3,i,j,s,e) = xyz(1:3) enddo enddo endif enddo ! Boundary-condition names, in HOHQMesh face order read(iUnit,*) bdyNames(1:6) do f = 1,6 s = hohqToSelfSide(f) if(trim(adjustl(bdyNames(f))) == "---") then ismBCid(s,e) = 0 else ! lookup/insert bdy name bcIdx = 0 do i = 1,nBCsLocal if(trim(BCNamesLocal(i)) == trim(adjustl(bdyNames(f)))) then bcIdx = i; exit endif enddo if(bcIdx == 0) then nBCsLocal = nBCsLocal+1 if(nBCsLocal > size(BCNamesLocal)) call grow_bc_table(BCNamesLocal) BCNamesLocal(nBCsLocal) = trim(adjustl(bdyNames(f))) bcIdx = nBCsLocal endif ismBCid(s,e) = bcIdx endif enddo enddo close(iUnit) print*,__FILE__//' : n curved faces = ',nFaceRecords, & ' n materials = ',nMatsLocal,' n boundary names = ',nBCsLocal ! At least one BC slot must exist so that Init allocates BCNames/BCType nBCsLocal = max(nBCsLocal,1) ! Set up the domain decomposition arrays (elemToRank, offsetElem, ! request/stat slots) using the file's element count. This is ! required for SideExchange even in the serial single-rank case. call this%decomp%GenerateDecomposition(nElemFile,6*nElemFile) ! ---- 4. Allocate SELF Mesh3D_t and populate ---- call this%Init(nGeo,nElemFile,6*nElemFile,nElemFile*ng1**3,nBCsLocal) this%nUniqueSides = 0 ! filled below after face matching this%quadrature = CHEBYSHEV_GAUSS_LOBATTO ! HOHQMesh face samples are at CGL points this%BCType = 0 do i = 1,nBCsLocal if(BCNamesLocal(i) /= "") then this%BCNames(i) = BCNamesLocal(i) else this%BCNames(i) = "unused" endif enddo ! Replace the default single-material table with the parsed one deallocate(this%materialNames) this%nMaterials = max(nMatsLocal,1) allocate(this%materialNames(1:this%nMaterials)) if(nMatsLocal == 0) then this%materialNames(1) = "default" else this%materialNames(1:nMatsLocal) = matNamesLocal(1:nMatsLocal) endif this%elemMaterial = ismMat ! Place corner nodes and run transfinite interpolation per element do e = 1,nElemFile call build_nodeCoords_for_hex(this,e,nGeo,ismCorners,ismFlag,faceCurve,nodeXYZ) ! Synthesize globalNodeIDs: stamp the eight corners with their ! file IDs and leave interior IDs as 0 (interior nodes are ! private to the element under our tensor product layout). this%globalNodeIDs(:,:,:,e) = 0 do l = 1,8 i = this%CGNSCornerMap(1,l) j = this%CGNSCornerMap(2,l) k = this%CGNSCornerMap(3,l) this%globalNodeIDs(i,j,k,e) = ismCorners(l,e) enddo ! Pack elemInfo with simple placeholders; SELF's 3D path does not ! depend on the HOPR-style offset fields when the mesh comes from ! a non-HOPR reader. this%elemInfo(1,e) = 0 this%elemInfo(2,e) = ismMat(e) ! Zone = material id this%elemInfo(3,e) = 6*(e-1) this%elemInfo(4,e) = 6*e this%elemInfo(5,e) = ng1**3*(e-1) this%elemInfo(6,e) = ng1**3*e enddo ! ---- 5. Build sideInfo via face corner matching ---- ! Each face is keyed by its four sorted corner-node ids; two local ! faces with the same key are the two sides of an interior face. allocate(pairKey(1:4,1:6*nElemFile)) allocate(pairElem(1:6*nElemFile),pairSide(1:6*nElemFile)) do e = 1,nElemFile do iSide = 1,6 do l = 1,4 quadA(l) = ismCorners(this%sideMap(l,iSide),e) enddo call sort4(quadA,sortedA) i = iSide+6*(e-1) pairKey(1:4,i) = sortedA pairElem(i) = e pairSide(i) = iSide enddo enddo ! Hash chain by smallest corner node id hashSize = max(nNodesFile,6*nElemFile)+1 allocate(hashHead(0:hashSize-1),hashNext(1:6*nElemFile)) hashHead = 0 hashNext = 0 do i = 1,6*nElemFile hashKey = pairKey(1,i) bucket = modulo(hashKey,hashSize) hashNext(i) = hashHead(bucket) hashHead(bucket) = i enddo this%sideInfo = 0 this%nUniqueSides = 0 do e = 1,nElemFile do iSide = 1,6 i = iSide+6*(e-1) ePair = 0; sPair = 0 bucket = modulo(pairKey(1,i),hashSize) probe = hashHead(bucket) do while(probe /= 0) if(probe /= i) then if(all(pairKey(1:4,probe) == pairKey(1:4,i))) then ePair = pairElem(probe) sPair = pairSide(probe) exit endif endif probe = hashNext(probe) enddo flip = 0 if(ePair /= 0) then ! Determine the flip by matching the corner-node orderings of ! the paired faces (face-index corner order per sideMap). do l = 1,4 quadA(l) = ismCorners(this%sideMap(l,iSide),e) quadB(l) = ismCorners(this%sideMap(l,sPair),ePair) enddo flip = -1 do p = 0,7 if(quadB(flipPerm(1,p)) == quadA(1) .and. & quadB(flipPerm(2,p)) == quadA(2) .and. & quadB(flipPerm(3,p)) == quadA(3) .and. & quadB(flipPerm(4,p)) == quadA(4)) then flip = p exit endif enddo if(flip < 0) then print*,__FILE__//' : Inconsistent face corner ordering between elements ', & e,' and ',ePair stop 1 endif endif this%sideInfo(3,iSide,e) = ePair this%sideInfo(4,iSide,e) = 10*sPair+flip this%sideInfo(5,iSide,e) = ismBCid(iSide,e) ! Allocate a globalSideID (count each shared face once) if(ePair == 0 .or. e < ePair) then this%nUniqueSides = this%nUniqueSides+1 this%sideInfo(2,iSide,e) = this%nUniqueSides else this%sideInfo(2,iSide,e) = this%sideInfo(2,sPair,ePair) endif enddo enddo deallocate(hashHead,hashNext,pairKey,pairElem,pairSide) deallocate(nodeXYZ,ismCorners,ismFlag,ismBCid,ismMat,faceCurve) deallocate(BCNamesLocal,matNamesLocal) call this%UpdateDevice() contains subroutine sort4(a,b) !! Ascending insertion sort of four integers implicit none integer,intent(in) :: a(1:4) integer,intent(out) :: b(1:4) integer :: m,n,tmp b = a do m = 2,4 tmp = b(m) n = m-1 do while(n >= 1) if(b(n) <= tmp) exit b(n+1) = b(n) n = n-1 enddo b(n+1) = tmp enddo endsubroutine sort4 subroutine grow_bc_table(tbl) implicit none character(LEN=255),allocatable,intent(inout) :: tbl(:) character(LEN=255),allocatable :: tmp(:) integer :: oldSize oldSize = size(tbl) allocate(tmp(1:2*oldSize)) tmp(1:oldSize) = tbl(1:oldSize) tmp(oldSize+1:) = "" call move_alloc(tmp,tbl) endsubroutine grow_bc_table subroutine grow_string_table(tbl) implicit none character(LEN=SELF_MESH_MATNAME_LENGTH),allocatable,intent(inout) :: tbl(:) character(LEN=SELF_MESH_MATNAME_LENGTH),allocatable :: tmp(:) integer :: oldSize oldSize = size(tbl) allocate(tmp(1:2*oldSize)) tmp(1:oldSize) = tbl(1:oldSize) tmp(oldSize+1:) = "" call move_alloc(tmp,tbl) endsubroutine grow_string_table endsubroutine Read_HOHQMesh_Mesh3D_t subroutine build_nodeCoords_for_hex(mesh,e,nGeo,allCorners,allFlags,allFaces,nodeXYZ) !! Fill mesh%nodeCoords(:,:,:,:,e) for one hexahedral element by !! transfinite (Coons) interpolation of its six face grids. Faces !! flagged in the mesh file use the file's face-point grids; the !! remaining faces are bilinear patches of their four corner nodes !! evaluated at Chebyshev-Gauss-Lobatto parametric coordinates. !! Edge curves are extracted from the face grids (preferring a !! flagged face when an edge borders one flagged and one bilinear !! face) and the standard Boolean-sum formula !! x = Pxi + Peta + Pzeta - Pxi*Peta - Pxi*Pzeta - Peta*Pzeta !! + Pxi*Peta*Pzeta !! combines face, edge, and corner contributions. For an element !! with all-straight faces this reduces to trilinear interpolation !! of the eight corners. implicit none class(Mesh3D_t),intent(inout) :: mesh integer,intent(in) :: e integer,intent(in) :: nGeo integer,intent(in) :: allCorners(:,:) ! 8 x nElem corner-node ids integer,intent(in) :: allFlags(:,:) ! 6 x nElem face flags (SELF side order) real(prec),intent(in) :: allFaces(:,:,:,:,:) ! 3 x nGeo+1 x nGeo+1 x 6 x nElem real(prec),intent(in) :: nodeXYZ(:,:) ! Local integer :: ng1 integer :: i,j,k,s,l integer :: cornerIDs(1:8) integer :: flag(1:6) real(prec) :: P(1:3,1:8) real(prec) :: face(1:3,1:nGeo+1,1:nGeo+1,1:6) real(prec) :: exEdge(1:3,1:nGeo+1,1:4) ! xi-directed edges: bs, bn, ts, tn real(prec) :: eyEdge(1:3,1:nGeo+1,1:4) ! eta-directed edges: bw, be, tw, te real(prec) :: ezEdge(1:3,1:nGeo+1,1:4) ! zeta-directed edges: sw, se, nw, ne real(prec) :: xi(0:nGeo),wts(0:nGeo) real(prec) :: u(1:nGeo+1) real(prec) :: a,b,c real(prec) :: xfp(1:3),xep(1:3),xcp(1:3) ng1 = nGeo+1 cornerIDs = allCorners(1:8,e) flag = allFlags(1:6,e) do l = 1,8 P(1:3,l) = nodeXYZ(1:3,cornerIDs(l)) enddo if(nGeo == 1) then ! Pure trilinear element with the 8 corners do l = 1,8 i = mesh%CGNSCornerMap(1,l) j = mesh%CGNSCornerMap(2,l) k = mesh%CGNSCornerMap(3,l) mesh%nodeCoords(1:3,i,j,k,e) = P(1:3,l) enddo return endif ! Chebyshev-Gauss-Lobatto parametric coordinates on [-1,1] and ! their [0,1] image used in the blending weights. We use SELF's ! quadrature routine to keep node placement consistent with the ! rest of the code. call ChebyshevQuadrature(nGeo,xi,wts,CHEBYSHEV_GAUSS_LOBATTO) do i = 1,ng1 u(i) = 0.5_prec*(xi(i-1)+1.0_prec) enddo ! Face grids: curved faces from the file, straight faces as ! bilinear patches of their four corners (face-index corner order ! per sideMap; all faces use the natural on-face volume axes). do s = 1,6 if(flag(s) == 1) then face(1:3,1:ng1,1:ng1,s) = allFaces(1:3,1:ng1,1:ng1,s,e) else do j = 1,ng1 do i = 1,ng1 a = u(i) b = u(j) face(1:3,i,j,s) = (1.0_prec-a)*(1.0_prec-b)*P(1:3,mesh%sideMap(1,s))+ & a*(1.0_prec-b)*P(1:3,mesh%sideMap(2,s))+ & a*b*P(1:3,mesh%sideMap(3,s))+ & (1.0_prec-a)*b*P(1:3,mesh%sideMap(4,s)) enddo enddo endif enddo ! Edge curves extracted from the face grids. Each edge borders two ! faces; a flagged (curved) face is preferred so that straight-face ! bilinear patches never override curved edge data. When both ! bordering faces are straight, the two restrictions coincide (the ! straight line between the edge's corner nodes). SELF side ids: ! 1=bottom, 2=south, 3=east, 4=north, 5=west, 6=top. do i = 1,ng1 ! xi-directed edges if(flag(1) == 1 .or. flag(2) /= 1) then ! bottom-south: bottom or south exEdge(1:3,i,1) = face(1:3,i,1,1) else exEdge(1:3,i,1) = face(1:3,i,1,2) endif if(flag(1) == 1 .or. flag(4) /= 1) then ! bottom-north: bottom or north exEdge(1:3,i,2) = face(1:3,i,ng1,1) else exEdge(1:3,i,2) = face(1:3,i,1,4) endif if(flag(6) == 1 .or. flag(2) /= 1) then ! top-south: top or south exEdge(1:3,i,3) = face(1:3,i,1,6) else exEdge(1:3,i,3) = face(1:3,i,ng1,2) endif if(flag(6) == 1 .or. flag(4) /= 1) then ! top-north: top or north exEdge(1:3,i,4) = face(1:3,i,ng1,6) else exEdge(1:3,i,4) = face(1:3,i,ng1,4) endif ! eta-directed edges if(flag(1) == 1 .or. flag(5) /= 1) then ! bottom-west: bottom or west eyEdge(1:3,i,1) = face(1:3,1,i,1) else eyEdge(1:3,i,1) = face(1:3,i,1,5) endif if(flag(1) == 1 .or. flag(3) /= 1) then ! bottom-east: bottom or east eyEdge(1:3,i,2) = face(1:3,ng1,i,1) else eyEdge(1:3,i,2) = face(1:3,i,1,3) endif if(flag(6) == 1 .or. flag(5) /= 1) then ! top-west: top or west eyEdge(1:3,i,3) = face(1:3,1,i,6) else eyEdge(1:3,i,3) = face(1:3,i,ng1,5) endif if(flag(6) == 1 .or. flag(3) /= 1) then ! top-east: top or east eyEdge(1:3,i,4) = face(1:3,ng1,i,6) else eyEdge(1:3,i,4) = face(1:3,i,ng1,3) endif ! zeta-directed edges if(flag(2) == 1 .or. flag(5) /= 1) then ! south-west: south or west ezEdge(1:3,i,1) = face(1:3,1,i,2) else ezEdge(1:3,i,1) = face(1:3,1,i,5) endif if(flag(2) == 1 .or. flag(3) /= 1) then ! south-east: south or east ezEdge(1:3,i,2) = face(1:3,ng1,i,2) else ezEdge(1:3,i,2) = face(1:3,1,i,3) endif if(flag(4) == 1 .or. flag(5) /= 1) then ! north-west: north or west ezEdge(1:3,i,3) = face(1:3,1,i,4) else ezEdge(1:3,i,3) = face(1:3,ng1,i,5) endif if(flag(4) == 1 .or. flag(3) /= 1) then ! north-east: north or east ezEdge(1:3,i,4) = face(1:3,ng1,i,4) else ezEdge(1:3,i,4) = face(1:3,ng1,i,3) endif enddo ! Boolean-sum transfinite interpolation. u,v,w in [0,1] are the ! blending parameters along xi, eta, zeta. do k = 1,ng1 c = u(k) do j = 1,ng1 b = u(j) do i = 1,ng1 a = u(i) ! Face contributions (Pxi + Peta + Pzeta) xfp = (1.0_prec-a)*face(1:3,j,k,5)+a*face(1:3,j,k,3)+ & (1.0_prec-b)*face(1:3,i,k,2)+b*face(1:3,i,k,4)+ & (1.0_prec-c)*face(1:3,i,j,1)+c*face(1:3,i,j,6) ! Edge corrections (Pxi*Peta + Pxi*Pzeta + Peta*Pzeta) xep = (1.0_prec-a)*(1.0_prec-b)*ezEdge(1:3,k,1)+ & a*(1.0_prec-b)*ezEdge(1:3,k,2)+ & a*b*ezEdge(1:3,k,4)+ & (1.0_prec-a)*b*ezEdge(1:3,k,3)+ & (1.0_prec-a)*(1.0_prec-c)*eyEdge(1:3,j,1)+ & a*(1.0_prec-c)*eyEdge(1:3,j,2)+ & a*c*eyEdge(1:3,j,4)+ & (1.0_prec-a)*c*eyEdge(1:3,j,3)+ & (1.0_prec-b)*(1.0_prec-c)*exEdge(1:3,i,1)+ & b*(1.0_prec-c)*exEdge(1:3,i,2)+ & b*c*exEdge(1:3,i,4)+ & (1.0_prec-b)*c*exEdge(1:3,i,3) ! Corner contributions (Pxi*Peta*Pzeta) xcp = (1.0_prec-a)*(1.0_prec-b)*(1.0_prec-c)*P(1:3,1)+ & a*(1.0_prec-b)*(1.0_prec-c)*P(1:3,2)+ & a*b*(1.0_prec-c)*P(1:3,3)+ & (1.0_prec-a)*b*(1.0_prec-c)*P(1:3,4)+ & (1.0_prec-a)*(1.0_prec-b)*c*P(1:3,5)+ & a*(1.0_prec-b)*c*P(1:3,6)+ & a*b*c*P(1:3,7)+ & (1.0_prec-a)*b*c*P(1:3,8) mesh%nodeCoords(1:3,i,j,k,e) = xfp-xep+xcp enddo enddo enddo endsubroutine build_nodeCoords_for_hex subroutine Write_Mesh3D_t(this,meshFile) ! Writes mesh output in HOPR format (serial only) implicit none class(Mesh3D_t),intent(inout) :: this character(*),intent(in) :: meshFile ! Local integer(HID_T) :: fileId call Open_HDF5(meshFile,H5F_ACC_RDWR_F,fileId) call WriteAttribute_HDF5(fileId,'nElems',this%nElem) call WriteAttribute_HDF5(fileId,'Ngeo',this%nGeo) call WriteAttribute_HDF5(fileId,'nBCs',this%nBCs) call WriteArray_HDF5(fileId,'BCType',this%bcType) call WriteArray_HDF5(fileId,'ElemInfo',this%elemInfo) ! Read local subarray of NodeCoords and GlobalNodeIDs call WriteArray_HDF5(fileId,'NodeCoords',this%nodeCoords) call WriteArray_HDF5(fileId,'GlobalNodeIDs',this%globalNodeIDs) ! Read local subarray of SideInfo call WriteArray_HDF5(fileId,'SideInfo',this%sideInfo) call Close_HDF5(fileID) endsubroutine Write_Mesh3D_t endmodule SELF_Mesh_3D_t