SELF_MeshRefinement_2D.f90 Source File

! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// !
!
! 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_2D
!! Mesh-level h-refinement operations for 2-D quadrilateral meshes (AMR Stage 2).
!!
!! This module assembles a refined Mesh2D_t from the element-local primitives in
!! SELF_RefinementPrimitives_2D. The currently supported operation is uniform refinement -
!! splitting every element 2:1 in each direction into four children - which keeps the mesh
!! conforming (no hanging nodes) and therefore needs neither the mortar machinery nor 2:1
!! balancing. Adaptive (flagged) refinement with hanging nodes and mortar regeneration is AMR
!! Stage 4; dynamic MPI re-partitioning of the refined mesh is Stage 5. Uniform refinement is
!! therefore restricted to a single rank here, and guards against multi-rank input.

  use SELF_Constants
  use SELF_Lagrange
  use SELF_Mesh_2D
  use SELF_RefinementPrimitives_2D

  implicit none

contains

  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

endmodule SELF_MeshRefinement_2D