UniformRefineMesh Subroutine

  subroutine UniformRefineMesh(meshIn,meshOut)
    !! Produce meshOut, a uniformly refined copy of meshIn: every element is split into
    !! eight 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 Mesh3D_t with 8x the elements, ready for
    !! geometry generation and time stepping. Refined globalNodeIDs are left zero: the
    !! solver reads flips from sideInfo directly and no 3-D consumer requires them.
    !!
    !! Serial only (matching the 2-D restriction); a multi-rank meshIn is rejected.
    implicit none
    type(Mesh3D),intent(in) :: meshIn
    type(Mesh3D),intent(out) :: meshOut
    ! Local
    integer :: nGeo,nElem,nElemR,nBCs
    integer :: p,c,k,eBase,nUniqueSidesRef
    type(Lagrange) :: geomInterp
    integer,allocatable :: refSideInfo(:,:,:)
    real(prec),allocatable :: childCoords(:,:,:,:,:)

    if(meshIn%decomp%nRanks > 1) then
      print*,__FILE__,':',__LINE__, &
        ' : Error : UniformRefineMesh is serial-only; adaptive multi-rank refinement'// &
        ' goes through the octree forest and EmitMesh.'
      stop 1
    endif

    nGeo = meshIn%nGeo
    nElem = meshIn%nElem
    nElemR = 8*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)

    ! Build refined connectivity (pure integer, deterministic).
    allocate(refSideInfo(1:5,1:6,1:nElemR))
    call RefineConnectivity(nElem,meshIn%sideInfo,refSideInfo,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).
    call meshOut%decomp%Init(comm=meshIn%decomp%mpiComm)
    call meshOut%decomp%GenerateDecomposition(nElemR,nUniqueSidesRef)
    call meshOut%Init(nGeo,nElemR,6*nElemR,8*nElemR,nBCs)
    meshOut%nUniqueSides = nUniqueSidesRef
    meshOut%quadrature = meshIn%quadrature

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

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

    ! Global node ids are not propagated (see the routine docstring).
    meshOut%globalNodeIDs = 0
    meshOut%elemInfo = 0

    ! Boundary-condition metadata carries over unchanged (face 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,8
        meshOut%elemMaterial(8*(p-1)+c) = meshIn%elemMaterial(p)
      enddo
    enddo

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

    call meshOut%UpdateDevice()

  endsubroutine UniformRefineMesh