Read_HOHQMesh_Mesh3D_t Subroutine

  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
      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)
      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)
      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