! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// !
!
! Maintainers : support@fluidnumerics.com
! Official Repository : https://github.com/FluidNumerics/self/
!
! Copyright © 2024 Fluid Numerics LLC
!
! Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
!
! 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
!
! 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in
!    the documentation and/or other materials provided with the distribution.
!
! 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from
!    this software without specific prior written permission.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS “AS IS” AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
! HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
! LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
! THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
! THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
!
! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// !

module SELF_MeshRefinement_3D
!! Mesh-level h-refinement operations for 3-D hexahedral meshes, the direct analogue of
!! SELF_MeshRefinement_2D.
!!
!! This module assembles a refined Mesh3D_t from the element-local primitives in
!! SELF_RefinementPrimitives_3D. The currently supported operation is uniform refinement -
!! splitting every element 2:1 in each direction into eight children - which keeps the
!! mesh conforming (no hanging faces) and therefore needs neither the mortar machinery nor
!! 2:1 balancing. Adaptive (flagged) refinement with hanging faces and mortar generation
!! lives in SELF_AdaptiveMesh_3D. Uniform refinement is restricted to a single rank here,
!! and guards against multi-rank input.

  use SELF_Constants
  use SELF_Lagrange
  use SELF_Mesh_3D
  use SELF_RefinementPrimitives_3D

  implicit none

contains

  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

endmodule SELF_MeshRefinement_3D
