UniformRefineMesh Subroutine

  subroutine UniformRefineMesh(meshIn,meshOut)
    !! Produce meshOut, a uniformly refined copy of meshIn: every element is split into four
    !! children. Child geometry is generated by exact isoparametric subdivision of the parent
    !! geometry (so the refined mesh represents the identical curved domain), and connectivity /
    !! flips are inherited deterministically from the base mesh. meshOut is a fully-formed,
    !! conforming Mesh2D_t with 4x the elements, ready for geometry generation and time stepping.
    !!
    !! Serial only (pending AMR Stage 5); a multi-rank meshIn is rejected.
    implicit none
    type(Mesh2D),intent(in) :: meshIn
    type(Mesh2D),intent(out) :: meshOut
    ! Local
    integer :: nGeo,nElem,nElemR,nBCs
    integer :: p,c,k,eBase,nodeOffset,nUniqueSidesRef
    type(Lagrange) :: geomInterp
    integer,allocatable :: baseCorner(:,:)
    integer,allocatable :: refSideInfo(:,:,:)
    integer,allocatable :: refCorner(:,:)
    real(prec),allocatable :: childCoords(:,:,:,:)

    if(meshIn%decomp%nRanks > 1) then
      print*,__FILE__,':',__LINE__, &
        ' : Error : UniformRefineMesh is serial-only pending AMR Stage 5 (MPI repartitioning).'
      stop 1
    endif

    nGeo = meshIn%nGeo
    nElem = meshIn%nElem
    nElemR = 4*nElem
    nBCs = meshIn%nBCs

    ! Geometry interpolant over the mesh's geometry nodes (degree nGeo, mesh quadrature nodes).
    call geomInterp%Init(nGeo,meshIn%quadrature,nGeo,meshIn%quadrature)

    ! Extract base element corner node ids (SW,SE,NE,NW) from globalNodeIDs.
    allocate(baseCorner(1:4,1:nElem))
    do p = 1,nElem
      do k = 1,4
        baseCorner(k,p) = meshIn%globalNodeIDs(meshIn%CGNSCornerMap(1,k), &
                                               meshIn%CGNSCornerMap(2,k),p)
      enddo
    enddo
    nodeOffset = maxval(baseCorner)

    ! Build refined connectivity and corner ids (pure integer, deterministic).
    allocate(refSideInfo(1:5,1:4,1:nElemR))
    allocate(refCorner(1:4,1:nElemR))
    call RefineConnectivity(nElem,meshIn%sideInfo,baseCorner,nodeOffset,meshIn%nUniqueSides, &
                            refSideInfo,refCorner,nUniqueSidesRef)

    ! Set up the (serial) decomposition for the output mesh. Initialize it on the input mesh's
    ! communicator: MPI is already up, so AcquireMPI reuses that communicator without a second
    ! mpi_init and registers the new decomposition with the process-wide lifecycle counter (so
    ! freeing meshOut later will not pull MPI out from under meshIn). Multi-rank refinement is AMR
    ! Stage 5.
    call meshOut%decomp%Init(comm=meshIn%decomp%mpiComm)
    call meshOut%decomp%GenerateDecomposition(nElemR,nUniqueSidesRef)
    call meshOut%Init(nGeo,nElemR,4*nElemR,4*nElemR,nBCs)
    meshOut%nUniqueSides = nUniqueSidesRef
    meshOut%quadrature = meshIn%quadrature

    ! Child geometry by isoparametric subdivision.
    allocate(childCoords(1:2,1:nGeo+1,1:nGeo+1,1:4))
    do p = 1,nElem
      eBase = 4*(p-1)
      call SubdivideNodeCoords(geomInterp,nGeo,meshIn%nodeCoords(1:2,1:nGeo+1,1:nGeo+1,p), &
                               childCoords)
      do c = 1,4
        meshOut%nodeCoords(1:2,1:nGeo+1,1:nGeo+1,eBase+c) = childCoords(1:2,1:nGeo+1,1:nGeo+1,c)
      enddo
    enddo

    ! Connectivity.
    meshOut%sideInfo(1:5,1:4,1:nElemR) = refSideInfo(1:5,1:4,1:nElemR)

    ! Corner global node ids (interior geometry-node ids are private/zero, as elsewhere in SELF).
    meshOut%globalNodeIDs = 0
    do p = 1,nElemR
      do k = 1,4
        meshOut%globalNodeIDs(meshOut%CGNSCornerMap(1,k),meshOut%CGNSCornerMap(2,k),p) = &
          refCorner(k,p)
      enddo
    enddo

    ! Boundary-condition metadata carries over unchanged (side BC ids are already in sideInfo).
    if(nBCs > 0) then
      meshOut%BCType(1:4,1:nBCs) = meshIn%BCType(1:4,1:nBCs)
      do k = 1,nBCs
        meshOut%BCNames(k) = meshIn%BCNames(k)
      enddo
    endif

    ! Material table: children inherit their parent's material id.
    meshOut%nMaterials = meshIn%nMaterials
    if(allocated(meshOut%materialNames)) deallocate(meshOut%materialNames)
    allocate(meshOut%materialNames(1:meshIn%nMaterials))
    meshOut%materialNames(1:meshIn%nMaterials) = meshIn%materialNames(1:meshIn%nMaterials)
    do p = 1,nElem
      do c = 1,4
        meshOut%elemMaterial(4*(p-1)+c) = meshIn%elemMaterial(p)
      enddo
    enddo

    deallocate(baseCorner,refSideInfo,refCorner,childCoords)
    call geomInterp%Free()

    call meshOut%UpdateDevice()

  endsubroutine UniformRefineMesh