BuildGeometry_AMRController2D Subroutine

  subroutine BuildGeometry_AMRController2D(this,newMesh,plan,newGeom,nReused)
    !! Fill newGeom for the emitted mesh, reusing the previous epoch's geometry for every element
    !! that did not change and generating only the rest (AMR Stage 6c). nReused reports how many
    !! elements were copied rather than computed.
    !!
    !! Why the reuse is exact. The transfer plan marks a new leaf SELF_TRANSFER_COPY only when the
    !! walk depth is zero, i.e. when the new element and old element sourceElem(li) are the SAME
    !! forest node (see BuildTransferPlan). A leaf's mesh node coordinates come from LeafCoords,
    !! which is a pure function of the root element's coordinates, the leaf's level and its
    !! quadrant path, evaluated in a fixed order; root coordinates are never mutated and node ids,
    !! levels and quadrants are stable across forest mutations. So an unchanged leaf's coordinates
    !! are bit-identical between epochs, and because per-element geometry generation touches only
    !! that element's own coordinates, its whole geometry block is too.
    !!
    !! Multi-rank: sourceElem indexes the GLOBAL old element list while each rank holds only its
    !! own slice, so reuse additionally requires the source to be locally owned. That is a range
    !! test against the OLD decomposition, which is still valid here because Regrid has not run
    !! yet. On one rank every COPY element qualifies and the test is always true; on several ranks
    !! whatever migrated is simply regenerated. No communication is added either way.
    implicit none
    class(AMRController2D),intent(inout) :: this
    type(Mesh2D),intent(in) :: newMesh
    type(TransferPlan2D),intent(in) :: plan
    type(SEMQuad),intent(inout) :: newGeom
    integer,intent(out) :: nReused
    ! Local
    integer :: li,gi,src,nLocal,nGen,oldFirst,oldLast,eFirst,nGeo,k
    integer,allocatable :: srcIdx(:),dstIdx(:),genIdx(:)
    real(prec),allocatable :: genCoords(:,:,:,:)

    nLocal = newMesh%nElem
    nGeo = newMesh%nGeo
    eFirst = newMesh%decomp%offsetElem(newMesh%decomp%rankId+1)+1

    ! Rank-local range of the OLD element list, i.e. what activeGeom actually holds.
    oldFirst = this%activeMesh%decomp%offsetElem(this%activeMesh%decomp%rankId+1)+1
    oldLast = this%activeMesh%decomp%offsetElem(this%activeMesh%decomp%rankId+2)

    allocate(srcIdx(1:nLocal),dstIdx(1:nLocal),genIdx(1:nLocal))

    call ResolveGeomDebug()

    ! Diagnostic bypass: reproduce the pre-6c behaviour exactly (full regeneration on the target
    ! buffer) while keeping the persistent alternating buffers, to tell a defect in the
    ! incremental assembly apart from one in the buffer reuse itself.
    if(geomFull) then
      call newGeom%GenerateFromMesh(newMesh)
      nReused = 0
      deallocate(srcIdx,dstIdx,genIdx)
      return
    endif

    nReused = 0
    nGen = 0
    do li = 1,nLocal
      gi = eFirst+li-1 ! this element's index in the plan's global new-leaf arrays
      src = plan%sourceElem(gi)
      if(.not. geomNoReuse .and. plan%sourceKind(gi) == SELF_TRANSFER_COPY .and. &
         src >= oldFirst .and. src <= oldLast) then
        nReused = nReused+1
        srcIdx(nReused) = src-oldFirst+1 ! rank-local index into activeGeom
        dstIdx(nReused) = li
      else
        nGen = nGen+1
        genIdx(nGen) = li
      endif
    enddo

    ! Generate the changed elements, compacted, so the generation loops run over nGen elements
    ! instead of all of them. Their geometry is then scattered into place.
    if(nGen > 0) then
      allocate(genCoords(1:2,1:nGeo+1,1:nGeo+1,1:nGen))
      do k = 1,nGen
        genCoords(1:2,:,:,k) = newMesh%nodeCoords(1:2,:,:,genIdx(k))
      enddo

      if(.not. associated(this%genGeom)) allocate(this%genGeom)
      if(this%genGeom%nElem == 0) then
        call this%genGeom%Init(this%interp,nGen)
      else
        call this%genGeom%Resize(this%interp,nGen)
      endif
      call this%genGeom%GenerateFromNodeCoords(genCoords,nGeo,newMesh%quadrature,nGen)
      call this%genGeom%x%UpdateDevice()
      call this%genGeom%x%BoundaryInterp()
      call this%genGeom%x%UpdateHost()
      call this%genGeom%CalculateMetricTerms()

      do k = 1,nGen
        srcIdx(nReused+k) = k
        dstIdx(nReused+k) = genIdx(k)
      enddo
      call newGeom%CopyElements(this%genGeom,srcIdx(nReused+1:),dstIdx(nReused+1:),nGen)
      deallocate(genCoords)
    endif

    ! Carry the unchanged elements across from the previous epoch's geometry.
    if(nReused > 0) then
      call newGeom%CopyElements(this%activeGeom,srcIdx,dstIdx,nReused)
    endif

    call newGeom%UploadGeometry()

    if(geomVerify) call VerifyGeometry(this,newMesh,newGeom)

    deallocate(srcIdx,dstIdx,genIdx)

  endsubroutine BuildGeometry_AMRController2D