oversetInitialization.F90

Go to the documentation of this file.

module oversetinitialization
 
contains
 
    subroutine initializestatus(level, sps)
 
        ! This subroutine initializes the status variable for use in the
        ! overset process
 
        use constants
        use blockpointers, only: ndom, il, jl, kl, ib, jb, kb, status
        use oversetutilities, only: setiscompute, setiswalldonor, setisdonor, setisreceiver
        use utils, only: setpointers
        implicit none
 
        ! Input params
        integer(kind=intType), intent(in) :: level, sps
 
        ! Working parameters
        integer(kind=intType) :: nn, i, j, k
 
        do nn = 1, ndom
            call setpointers(nn, level, sps)
 
            ! Now loop over the owned cells and set the isCompute flag to
            ! true.
 
            do k = 2, kl
                do j = 2, jl
                    do i = 2, il
                        status(i, j, k) = 0
                        call setiscompute(status(i, j, k), .true.)
 
                        ! Additional initialization for full overset update mode
                        call setiswalldonor(status(i, j, k), .false.)
                        call setisdonor(status(i, j, k), .false.)
                        call setisreceiver(status(i, j, k), .false.)
 
                    end do
                end do
            end do
        end do
 
    end subroutine initializestatus
 
    subroutine reinitializestatus(level, sps)
 
        ! This subroutine reinitializes the status variable. However, if
        ! cell is a wallDonor, that information is kept.
 
        use constants
        use blockpointers, only: ndom, il, jl, kl, ib, jb, kb, status
        use oversetutilities, only: setiscompute, iswalldonor, setiswalldonor
        use utils, only: setpointers
        implicit none
 
        ! Input params
        integer(kind=intType), intent(in) :: level, sps
 
        ! Working parameters
        integer(kind=intType) :: nn, i, j, k
        logical :: wDonor
        do nn = 1, ndom
            call setpointers(nn, level, sps)
 
            ! Now loop over the owned cells and set the isCompute flag to
            ! true.
 
            do k = 2, kl
                do j = 2, jl
                    do i = 2, il
                        wdonor = iswalldonor(status(i, j, k))
                        status(i, j, k) = 0
                        call setiscompute(status(i, j, k), .true.)
                        call setiswalldonor(status(i, j, k), wdonor)
                    end do
                end do
            end do
        end do
    end subroutine reinitializestatus
 
    subroutine initializeoblock(oBlock, nn, level, sps)
 
        ! This routine allocates the data for the supplied oBlock using the
        !  data currently in blockPointers
        use constants
        use oversetdata, only: oversetblock, clusters, cumdomproc
        use inputoverset, only: backgroundvolscale, nearwalldist, useoversetwallscaling
        use blockpointers, only: x, globalcell, il, jl, kl, ib, jb, kb, &
                                 ie, je, ke, vol, iblank, xseed, forcedrecv, nbkglobal, si, sj, sk
        use cgnsgrid, only: cgnsdoms
        use adtbuild, only: buildserialhex
        use communication, only: myid
        use stencils, only: visc_drdw_stencil, n_visc_drdw
        use utils, only: mynorm2
        use oversetutilities, only: wallsonblock
        implicit none
 
        ! Input Params
        type(oversetblock), intent(inout) :: oBlock
        integer(kind=intType) :: nn, level, sps, kk
 
        ! Working paramters
        integer(kind=intType) :: i, j, k, mm, nADT, nHexa, planeOffset
        integer(kind=intType) :: iStart, iEnd, jStart, jEnd, kStart, kEnd
        real(kind=realtype) :: factor, frac, dist, xp(3), aspect(3), fact
        integer(kind=intType) :: i_stencil, ii, jj, iii
        logical :: wallsPresent
        logical, allocatable, dimension(:, :, :) :: nearWallTmp
 
        ! Set all the sizes for this block.
        oblock%il = il
        oblock%jl = jl
        oblock%kl = kl
 
        oblock%proc = myid
        oblock%block = nn
        oblock%cluster = clusters(cumdomproc(myid) + nn)
        call wallsonblock(wallspresent)
 
        ! Do the reset of the allocs
        allocate ( &
            oblock%qualDonor(1, ie * je * ke), &
            oblock%globalCell(0:ib, 0:jb, 0:kb), &
            oblock%invalidDonor(1:ie, 1:je, 1:ke))
 
        ! Invalid Donor array is simply if the cell is a forced receiver
        ! or not.
        do k = 1, ke
            do j = 1, je
                do i = 1, ie
                    oblock%invalidDonor(i, j, k) = forcedrecv(i, j, k)
                end do
            end do
        end do
 
        ! Compute the qualDonor depending on if we have a wall block or not.
        mm = 0
        do k = 1, ke
            do j = 1, je
                do i = 1, ie
                    mm = mm + 1
                    if (wallspresent) then
 
                        ! Modify based on aspect ratio of the cell in the
                        ! k-direcion. High aspect ratioin BL
                        aspect = one
                        if (useoversetwallscaling) then
                            if (cgnsdoms(nbkglobal)%viscousDir(1)) &
                                aspect(1) = (half * (mynorm2(si(i - 1, j, k, :)) + &
                                                     mynorm2(si(i, j, k, :)))) / vol(i, j, k)
                            if (cgnsdoms(nbkglobal)%viscousDir(2)) &
                                aspect(2) = (half * (mynorm2(sj(i, j - 1, k, :)) + &
                                                     mynorm2(sj(i, j, k, :)))) / vol(i, j, k)
                            if (cgnsdoms(nbkglobal)%viscousDir(3)) &
                                aspect(3) = (half * (mynorm2(sk(i, j, k - 1, :)) + &
                                                     mynorm2(sk(i, j, k, :)))) / vol(i, j, k)
                        end if
                        fact = min(aspect(1) * aspect(2) * aspect(3), 100.0_realtype)
 
                        oblock%qualDonor(1, mm) = (vol(i, j, k)**third) / fact
                    else
                        oblock%qualDonor(1, mm) = (backgroundvolscale * vol(i, j, k))**third
                    end if
 
                    ! Account for explicit scaling of quality
                    oblock%qualDonor(1, mm) = oblock%qualDonor(1, mm) * cgnsdoms(nbkglobal)%priority
                end do
            end do
        end do
 
        !Copy over global cell
        oblock%globalCell = globalcell
 
        ! Now setup the data for the ADT
        nhexa = il * jl * kl
        nadt = ie * je * ke
 
        allocate (oblock%xADT(3, nadt), oblock%hexaConn(8, nhexa))
        ! Fill up the xADT using cell centers (dual mesh)
        mm = 0
 
        ! Allocate the nearWall
        allocate (oblock%nearWall(1:il, 1:jl, 1:kl))
        oblock%nearWall = 0
 
        allocate (nearwalltmp(1:ie, 1:je, 1:ke))
        nearwalltmp = .false.
 
        do k = 1, ke
            do j = 1, je
                do i = 1, ie
                    mm = mm + 1
                    xp = eighth * ( &
                         x(i - 1, j - 1, k - 1, :) + &
                         x(i, j - 1, k - 1, :) + &
                         x(i - 1, j, k - 1, :) + &
                         x(i, j, k - 1, :) + &
                         x(i - 1, j - 1, k, :) + &
                         x(i, j - 1, k, :) + &
                         x(i - 1, j, k, :) + &
                         x(i, j, k, :))
                    oblock%xADT(:, mm) = xp
 
                    ! Determine if this point is near wall. Note that the
                    ! boundary halos sill have xSeed as "large" so these won't
                    ! be flagged as nearWall. We will account for this below.
                    dist = mynorm2(xp - xseed(i, j, k, :))
                    if (dist < nearwalldist) then
                        nearwalltmp(i, j, k) = .true.
                    end if
                end do
            end do
        end do
 
        ! Now finally set the nearwall for the dual mesh cells. It is
        ! considered a near wall if all "nodes" of the dual mesh cell are
        ! also near wall. Have to be carful not to count boundary halos
        ! since they do not have nearWallTmp Values.
 
        do k = 1, kl
            do j = 1, jl
                do i = 1, il
                    if ( &
                        (nearwalltmp(i, j, k) .or. globalcell(i, j, k) < 0) .and. &
                        (nearwalltmp(i + 1, j, k) .or. globalcell(i + 1, j, k) < 0) .and. &
                        (nearwalltmp(i, j + 1, k) .or. globalcell(i, j + 1, k) < 0) .and. &
                        (nearwalltmp(i + 1, j + 1, k) .or. globalcell(i + 1, j + 1, k) < 0) .and. &
                        (nearwalltmp(i, j, k + 1) .or. globalcell(i, j, k + 1) < 0) .and. &
                        (nearwalltmp(i + 1, j, k + 1) .or. globalcell(i + 1, j, k + 1) < 0) .and. &
                        (nearwalltmp(i, j + 1, k + 1) .or. globalcell(i, j + 1, k + 1) < 0) .and. &
                        (nearwalltmp(i + 1, j + 1, k + 1) .or. globalcell(i + 1, j + 1, k + 1) < 0)) then
                        oblock%nearWall(i, j, k) = 1
                    end if
                end do
            end do
        end do
 
        deallocate (nearwalltmp)
        mm = 0
        ! These are the 'elements' of the dual mesh.
        planeoffset = ie * je
        do k = 2, ke
            do j = 2, je
                do i = 2, ie
                    mm = mm + 1
                    oblock%hexaConn(1, mm) = (k - 2) * planeoffset + (j - 2) * ie + (i - 2) + 1
                    oblock%hexaConn(2, mm) = oblock%hexaConn(1, mm) + 1
                    oblock%hexaConn(3, mm) = oblock%hexaConn(2, mm) + ie
                    oblock%hexaConn(4, mm) = oblock%hexaConn(3, mm) - 1
 
                    oblock%hexaConn(5, mm) = oblock%hexaConn(1, mm) + planeoffset
                    oblock%hexaConn(6, mm) = oblock%hexaConn(2, mm) + planeoffset
                    oblock%hexaConn(7, mm) = oblock%hexaConn(3, mm) + planeoffset
                    oblock%hexaConn(8, mm) = oblock%hexaConn(4, mm) + planeoffset
                end do
            end do
        end do
 
        ! Call the custom build routine -- Serial only, only Hexa volumes,
        ! we supply our own ADT Type
 
        call buildserialhex(nhexa, nadt, oblock%xADT, oblock%hexaConn, oblock%ADT)
 
        ! Flag this block as being allocated
        oblock%allocated = .true.
 
    end subroutine initializeoblock
 
    subroutine initializeofringes(oFringe, nn, famList)
 
        ! This subroutine initializes the fringe information for the given
        ! block, level and spectral instance. It is assumed that
        ! blockPointers are already set.
        use constants
        use communication, only: myid
        use blockpointers
        use oversetdata, only: oversetfringe, clusters, cumdomproc
        use stencils, only: visc_drdw_stencil, n_visc_drdw
        use inputoverset, only: backgroundvolscale
        use sorting, only: faminlist
        use oversetutilities, only: wallsonblock, windindex, unwindindex
        implicit none
 
        ! Input Params
        type(oversetfringe), intent(inout) :: oFringe
        integer(kind=intType), intent(in) :: nn
        integer(kind=intType), intent(in), dimension(:) :: famList
        ! Working Params
        integer(kind=intTYpe) :: i, j, k, mm, iDim, ii, jj, kk, iii, jjj, myI, myJ, myK
        integer(kind=intTYpe) :: iStart, iEnd, jStart, jEnd, kStart, kEnd
        logical :: wallsPresent
        integer(kind=intType) :: i_stencil
        real(kind=realtype) :: dist, frac, xp(3)
        ! Check if we have walls:
        call wallsonblock(wallspresent)
 
        ! Set the sizes for the oFringe and allocate the required space.
        ofringe%il = il
        ofringe%jl = jl
        ofringe%kl = kl
        ofringe%nx = nx
        ofringe%ny = ny
        ofringe%nz = nz
        ofringe%block = nn
        ofringe%cluster = clusters(cumdomproc(myid) + nn)
        ofringe%proc = myid
 
        mm = nx * ny * nz
        allocate (ofringe%x(3, mm))
        allocate (ofringe%xSeed(3, mm))
        allocate (ofringe%wallInd(mm))
        allocate (ofringe%isWall(mm))
        ofringe%isWall = 0
        ofringe%xSeed = large
        ofringe%wallInd = 0
 
        ! Assume each cell will get just one donor. It's just a guess, it
        ! will be expanded if necessary so the exact value doesn't matter.
        allocate (ofringe%fringeIntBuffer(5, mm), ofringe%fringeRealBuffer(4, mm))
 
        ofringe%nDonor = 0
        ! Now loop over the actual compute cells, setting the cell center
        ! value 'x', the volume and flag these cells as compute
        ii = 0
        do k = 2, kl
            do j = 2, jl
                do i = 2, il
                    ii = ii + 1
                    do idim = 1, 3
                        ofringe%x(idim, ii) = eighth * ( &
                                              x(i - 1, j - 1, k - 1, idim) + &
                                              x(i, j - 1, k - 1, idim) + &
                                              x(i - 1, j, k - 1, idim) + &
                                              x(i, j, k - 1, idim) + &
                                              x(i - 1, j - 1, k, idim) + &
                                              x(i, j - 1, k, idim) + &
                                              x(i - 1, j, k, idim) + &
                                              x(i, j, k, idim))
                    end do
                    ofringe%xSeed(:, ii) = xseed(i, j, k, :)
                    ofringe%wallInd(ii) = wallind(i, j, k)
                    ofringe%fringeIntBuffer(4, ii) = nn
                    ofringe%fringeIntBuffer(5, ii) = windindex(i, j, k, il, jl, kl)
 
                end do
            end do
        end do
 
        ! We also need to flag a single layer of cells next a wall
        ! boundary condition as being "isWall". This information is
        ! necessary to be able to determine the "wall donors" which are
        ! the flood seeds.
 
        do mm = 1, nbocos
            select case (bcfaceid(mm))
            case (imin)
                istart = 2; iend = 2;
                jstart = bcdata(mm)%inBeg + 1; jend = bcdata(mm)%inEnd
                kstart = bcdata(mm)%jnBeg + 1; kend = bcdata(mm)%jnEnd
            case (imax)
                istart = il; iend = il;
                jstart = bcdata(mm)%inBeg + 1; jend = bcdata(mm)%inEnd
                kstart = bcdata(mm)%jnBeg + 1; kend = bcdata(mm)%jnEnd
            case (jmin)
                istart = bcdata(mm)%inBeg + 1; iend = bcdata(mm)%inEnd
                jstart = 2; jend = 2;
                kstart = bcdata(mm)%jnBeg + 1; kend = bcdata(mm)%jnEnd
            case (jmax)
                istart = bcdata(mm)%inBeg + 1; iend = bcdata(mm)%inEnd
                jstart = jl; jend = jl;
                kstart = bcdata(mm)%jnBeg + 1; kend = bcdata(mm)%jnEnd
            case (kmin)
                istart = bcdata(mm)%inBeg + 1; iend = bcdata(mm)%inEnd
                jstart = bcdata(mm)%jnBeg + 1; jend = bcdata(mm)%jnEnd
                kstart = 2; kend = 2;
            case (kmax)
                istart = bcdata(mm)%inBeg + 1; iend = bcdata(mm)%inEnd
                jstart = bcdata(mm)%jnBeg + 1; jend = bcdata(mm)%jnEnd
                kstart = kl; kend = kl;
            end select
 
            faminclude: if (faminlist(bcdata(mm)%famID, famlist)) then
                do k = kstart, kend
                    do j = jstart, jend
                        do i = istart, iend
                            ! Recompute the index
                            ii = (k - 2) * nx * ny + (j - 2) * nx + (i - 2) + 1
                            ofringe%isWall(ii) = bcfaceid(mm)
                        end do
                    end do
                end do
            end if faminclude
        end do ! BocoLoop
 
        ! Flag this set of fringes as being allocated
        ofringe%allocated = .true.
    end subroutine initializeofringes
 
    subroutine initializeosurf(famList, oSurf, dualMesh, cluster)
 
        ! This routine builds the ADT tree for any wall surfaces for the
        ! block currently being pointed to by block Pointers.
        use constants
        use oversetdata, only: oversetwall
        use blockpointers, only: nbocos, bcdata, bcfaceid, il, jl, kl, &
                                 bcfaceid, x, bctype, righthanded, nbklocal
        use adtbuild, only: buildserialquad
        use kdtree2_module, onlY: kdtree2_create
        use oversetpackingroutines, only: getwallsize
        use sorting, only: faminlist
        implicit none
 
        ! Input Params
        integer(kind=intType), intent(in), dimension(:) :: famList
        type(oversetwall), intent(inout) :: oSurf
        logical, intent(in) :: dualMesh
        integer(kind=intType), intent(in) :: cluster
 
        ! Working paramters
        integer(kind=intType) :: i, j, k, n, ii, jj, jjj, mm, ni, nj, nodeCount
        integer(kind=intType) :: iBeg, iEnd, jBeg, jEnd, nNodes, maxCells, nCells, iNode
        logical :: regularOrdering
 
        ! Set all the sizes for this block.
        osurf%il = il
        osurf%jl = jl
        osurf%kl = kl
 
        call getwallsize(famlist, nnodes, maxcells, dualmesh)
 
        osurf%nNodes = nnodes
        osurf%maxCells = maxcells
        osurf%cluster = cluster
        ! Allocate space for the x array and connectivity array. cellPtr is
        ! larger than necessary.
        allocate (osurf%x(3, nnodes), osurf%conn(4, maxcells), &
                  osurf%cellPtr(maxcells), osurf%iBlank(maxcells), &
                  osurf%delta(nnodes), osurf%nte(4, nnodes))
 
        ii = 0 ! Cumulative node counter
        jj = 0 ! Cumulative cell counter (with iblanks)
        jjj = 0 !Cumulative cell counter (without iblanks)
        nodecount = 0
 
        do mm = 1, nbocos
            faminclude: if (faminlist(bcdata(mm)%famID, famlist)) then
 
                select case (bcfaceid(mm))
                case (imin, jmax, kmin)
                    regularordering = .true.
                case default
                    regularordering = .false.
                end select
 
                ! Now this can be reversed *again* if we have a block that
                ! is left handed.
                if (.not. righthanded) then
                    regularordering = .not. (regularordering)
                end if
 
                ! THIS IS SUPER IMPORTANT: It is absolutely critical that the
                ! wall be built *FROM THE DUAL MESH!!* IT WILL NOT WORK IF YOU
                ! USE THE PRIMAL MESH! The -1 for the node ranges below gives
                ! the extra '1' node for the mesh formed from the dual cells.
 
                dualcheck: if (dualmesh) then
                    jbeg = bcdata(mm)%jnBeg - 1; jend = bcdata(mm)%jnEnd
                    ibeg = bcdata(mm)%inBeg - 1; iend = bcdata(mm)%inEnd
                    ! Now fill up the point array
                    do j = jbeg, jend
                        do i = ibeg, iend
                            ii = ii + 1
                            select case (bcfaceid(mm))
                            case (imin)
                                osurf%x(:, ii) = fourth * (x(1, i, j, :) + x(1, i + 1, j, :) + &
                                                           x(1, i, j + 1, :) + x(1, i + 1, j + 1, :))
                            case (imax)
                                osurf%x(:, ii) = fourth * (x(il, i, j, :) + x(il, i + 1, j, :) + &
                                                           x(il, i, j + 1, :) + x(il, i + 1, j + 1, :))
                            case (jmin)
                                osurf%x(:, ii) = fourth * (x(i, 1, j, :) + x(i + 1, 1, j, :) + &
                                                           x(i, 1, j + 1, :) + x(i + 1, 1, j + 1, :))
                            case (jmax)
                                osurf%x(:, ii) = fourth * (x(i, jl, j, :) + x(i + 1, jl, j, :) + &
                                                           x(i, jl, j + 1, :) + x(i + 1, jl, j + 1, :))
                            case (kmin)
                                osurf%x(:, ii) = fourth * (x(i, j, 1, :) + x(i + 1, j, 1, :) + &
                                                           x(i, j + 1, 1, :) + x(i + 1, j + 1, 1, :))
                            case (kmax)
                                osurf%x(:, ii) = fourth * (x(i, j, kl, :) + x(i + 1, j, kl, :) + &
                                                           x(i, j + 1, kl, :) + x(i + 1, j + 1, kl, :))
                            end select
                        end do
                    end do
 
                    ! Fill up the conn array. Note that don't take the
                    ! surface`normal direction (in or out) or the cell
                    ! handed-ness into account...it is not necessary since we
                    ! are just getting distance to the wall, which is
                    ! independent of the orientation.
 
                    ni = iend - ibeg + 1
                    nj = jend - jbeg + 1
                    do j = 0, nj - 2
                        do i = 0, ni - 2
                            jj = jj + 1
                            if (regularordering) then
                                osurf%conn(1, jj) = nodecount + (j) * ni + i + 1 ! n1
                                osurf%conn(2, jj) = nodecount + (j) * ni + i + 2 ! n2
                                osurf%conn(3, jj) = nodecount + (j + 1) * ni + i + 2 ! n3
                                osurf%conn(4, jj) = nodecount + (j + 1) * ni + i + 1 ! n4
                            else
                                osurf%conn(1, jj) = nodecount + (j) * ni + i + 1 ! n1
                                osurf%conn(2, jj) = nodecount + (j + 1) * ni + i + 1 ! n4
                                osurf%conn(3, jj) = nodecount + (j + 1) * ni + i + 2 ! n3
                                osurf%conn(4, jj) = nodecount + (j) * ni + i + 2 ! n2
                            end if
                        end do
                    end do
                    nodecount = nodecount + ni * nj
 
                    ! We don't care about iBlank, cellPtr or delta for the dual
                    ! mesh
                    osurf%iBlank = 1
                    osurf%cellPtr = 0
                    osurf%delta = zero
                else ! Using the primal mesh
                    jbeg = bcdata(mm)%jnBeg; jend = bcdata(mm)%jnEnd
                    ibeg = bcdata(mm)%inBeg; iend = bcdata(mm)%inEnd
 
                    ! Now fill up the point array. Owned node loop.
                    do j = jbeg, jend
                        do i = ibeg, iend
                            ii = ii + 1
                            select case (bcfaceid(mm))
                            case (imin)
                                osurf%x(:, ii) = x(1, i, j, :)
                            case (imax)
                                osurf%x(:, ii) = x(il, i, j, :)
                            case (jmin)
                                osurf%x(:, ii) = x(i, 1, j, :)
                            case (jmax)
                                osurf%x(:, ii) = x(i, jl, j, :)
                            case (kmin)
                                osurf%x(:, ii) = x(i, j, 1, :)
                            case (kmax)
                                osurf%x(:, ii) = x(i, j, kl, :)
                            end select
                            osurf%delta(ii) = bcdata(mm)%deltaNode(i, j)
                        end do
                    end do
 
                    ! Fill up the conn array being careful to *only* adding
                    ! cells that are not already blanked.
                    ni = iend - ibeg + 1
                    nj = jend - jbeg + 1
                    do j = 0, nj - 2
                        do i = 0, ni - 2
                            jjj = jjj + 1
                            osurf%iBlank(jjj) = bcdata(mm)%iblank(ibeg + i + 1, jbeg + j + 1)
                            if (osurf%iBlank(jjj) == 1) then
                                jj = jj + 1
                                if (regularordering) then
                                    osurf%conn(1, jj) = nodecount + (j) * ni + i + 1 ! n1
                                    osurf%conn(2, jj) = nodecount + (j) * ni + i + 2 ! n2
                                    osurf%conn(3, jj) = nodecount + (j + 1) * ni + i + 2 ! n3
                                    osurf%conn(4, jj) = nodecount + (j + 1) * ni + i + 1 ! n4
                                else
                                    osurf%conn(1, jj) = nodecount + (j) * ni + i + 1 ! n1
                                    osurf%conn(2, jj) = nodecount + (j + 1) * ni + i + 1 ! n4
                                    osurf%conn(3, jj) = nodecount + (j + 1) * ni + i + 2 ! n3
                                    osurf%conn(4, jj) = nodecount + (j) * ni + i + 2 ! n2
                                end if
                                osurf%cellPtr(jj) = jjj
                            end if
                        end do
                    end do
                    nodecount = nodecount + ni * nj
                end if dualcheck
            end if faminclude
        end do
 
        ! Set the actual number of cells
        osurf%nCells = jj
 
        ! Build the tree itself.
        call buildserialquad(osurf%nCells, nnodes, osurf%x, osurf%conn, osurf%ADT)
 
        ! Build the KDTree
        if (osurf%nNodes > 0) then
            osurf%tree => kdtree2_create(osurf%x)
        end if
 
        ! Build the inverse of the connectivity, the nodeToElem array.
        osurf%nte = 0
        do i = 1, osurf%nCells
            do j = 1, 4
                n = osurf%conn(j, i)
                inner: do k = 1, 4
                    if (osurf%nte(k, n) == 0) then
                        osurf%nte(k, n) = i
                        exit inner
                    end if
                end do inner
            end do
        end do
 
        ! Flag this wall as being allocated
        osurf%allocated = .true.
 
    end subroutine initializeosurf
 
end module oversetinitialization