subroutine Read_HOPr_Mesh2D_t(this,meshFile,enableDomainDecomposition)
! From https://www.hopr-project.org/externals/Meshformat.pdf, Algorithm 6
! Adapted for 2D Mesh : Note that HOPR does not have 2D mesh output.
implicit none
class(Mesh2D_t),intent(out) :: this
character(*),intent(in) :: meshFile
logical,intent(in),optional :: enableDomainDecomposition
! Local
integer(HID_T) :: fileId
integer(HID_T) :: offset(1:2),gOffset(1)
integer :: nGlobalElem
integer :: firstElem
integer :: firstNode
integer :: firstSide
integer :: nLocalElems
integer :: nLocalNodes3D
integer :: nLocalSides3D
integer :: nUniqueSides3D
integer :: nLocalNodes2D
integer :: nLocalSides2D
integer :: nUniqueSides2D
integer :: nGeo,nBCs
integer :: eid,lsid,iSide
integer :: i,j,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
if(present(enableDomainDecomposition)) then
call this%decomp%init(enableDomainDecomposition)
else
call this%decomp%init(.false.)
endif
print*,__FILE__//' : Reading HOPr mesh from'//trim(meshfile)
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
print*,__FILE__//' : Loading mesh attributes'
call ReadAttribute_HDF5(fileId,'nElems',nGlobalElem)
call ReadAttribute_HDF5(fileId,'Ngeo',nGeo)
call ReadAttribute_HDF5(fileId,'nBCs',nBCs)
call ReadAttribute_HDF5(fileId,'nUniqueSides',nUniqueSides3D)
print*,__FILE__//' : N Global Elements = ',nGlobalElem
print*,__FILE__//' : Mesh geometry degree = ',nGeo
print*,__FILE__//' : N Boundary conditions = ',nBCs
print*,__FILE__//' : N Unique Sides (3D) = ',nUniqueSides3D
! 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
print*,__FILE__//' : Generating Domain Decomposition'
call this%decomp%GenerateDecomposition(nGlobalElem,nUniqueSides3D)
firstElem = this%decomp%offsetElem(this%decomp%rankId+1)+1
nLocalElems = this%decomp%offsetElem(this%decomp%rankId+2)- &
this%decomp%offsetElem(this%decomp%rankId+1)
print*,__FILE__//' : Rank ',this%decomp%rankId+1,' : element offset = ',firstElem
print*,__FILE__//' : Rank ',this%decomp%rankId+1,' : n_elements = ',nLocalElems
! 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
nLocalNodes3D = hopr_elemInfo(6,nLocalElems)-hopr_elemInfo(5,1)
! Allocate Space for hopr_nodeCoords and hopr_globalNodeIDs !
allocate(hopr_nodeCoords(1:3,nLocalNodes3D),hopr_globalNodeIDs(1:nLocalNodes3D))
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
nLocalSides3D = hopr_elemInfo(4,nLocalElems)-hopr_elemInfo(3,1)
! Allocate space for hopr_sideInfo
allocate(hopr_sideInfo(1:5,1:nLocalSides3D))
if(this%decomp%mpiEnabled) then
offset = (/0,firstSide-1/)
print*,__FILE__//' : Rank ',this%decomp%rankId+1,' Reading side information'
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 ---- !
! Now we need to convert from 3-D to 2-D !
nLocalSides2D = nLocalSides3D-2*nLocalElems
nUniqueSides2D = nUniqueSides3D-2*nGlobalElem ! Remove the "top" and "bottom" faces
nLocalNodes2D = nLocalNodes2D-nLocalElems*nGeo*(nGeo+1)**2 ! Remove the third dimension
print*,__FILE__//' : Rank ',this%decomp%rankId+1,' Allocating memory for mesh'
print*,__FILE__//' : Rank ',this%decomp%rankId+1,' n local sides : ',nLocalSides2D
call this%Init(nGeo,nLocalElems,nLocalSides2D,nLocalNodes2D,nBCs)
this%nUniqueSides = nUniqueSides2D ! Store the number of sides in the global mesh
! Copy data from local arrays into this
! elemInfo(1:6,iEl)
! 1 - Element Type
! 2 - Zone
! 3 - offset index for side array (not needed when all quads are assumed)
! 4 - last index for side array (not needed when all quads are assumed)
! 5 - offset index for node array (not needed when all quads are assumed)
! 6 - last index for node array (not needed when all quads are assumed)
this%elemInfo = hopr_elemInfo
this%quadrature = UNIFORM ! HOPr uses uniformly spaced points
! Grab the node coordinates (x and y only) from the "bottom" layer of the extruded mesh
do eid = 1,this%nElem
do j = 1,nGeo+1
do i = 1,nGeo+1
nid = i+(nGeo+1)*(j-1+(nGeo+1)*((nGeo+1)*(eid-1)))
this%nodeCoords(1:2,i,j,eid) = hopr_nodeCoords(1:2,nid)
this%globalNodeIDs(i,j,eid) = hopr_globalNodeIDs(nid)
enddo
enddo
enddo
! Grab the south, west, north, and south sides of the elements
! sideInfo(1:5,iSide,iEl)
!
! 1 - Side Type (currently unused in SELF)
! 2 - Global Side ID (Used for message passing. Don't need to change)
! 3 - Neighbor Element ID (Can stay the same)
! 4 - 10*( neighbor local side ) + flip (Need to recalculate flip)
! 5 - Boundary Condition ID (Can stay the same)
do eid = 1,this%nElem
do lsid = 1,4
! Calculate the 3-D side ID from the 2-D local side id and element ID
iSide = lsid+1+6*(eid-1)
this%sideInfo(1:5,lsid,eid) = hopr_sideInfo(1:5,iSide)
! Adjust the secondary side index for 2-D
this%sideInfo(4,lsid,eid) = this%sideInfo(4,lsid,eid)-10
enddo
enddo
call this%RecalculateFlip()
deallocate(hopr_elemInfo,hopr_nodeCoords,hopr_globalNodeIDs,hopr_sideInfo)
call this%UpdateDevice()
endsubroutine Read_HOPr_Mesh2D_t