overset.F90

Go to the documentation of this file.

module oversetdata
 
    use constants
    use adtdata, only: adttype
    use block, only: fringetype
    use kdtree2_module, only: kdtree2
#ifndef USE_TAPENADE
#include <petsc/finclude/petsc.h>
    use petsc
    implicit none
 
#endif
    ! Helper dataType for communicated overset grid points. This data
    ! structure mirrros the blockType structure in block.F90, but only
    ! contains minimum amount of information required for computing
    ! overset connectivities.
 
    ! Store the coordinates from a block that will need to be searched.
 
    ! A simple generic sparse matrix storage container for storing the
    ! (sparse) overlap structure of an overset mesh
    type csrmatrix
        integer(kind=intType) :: nrow, ncol, nnz, nnzlocal
        integer(kind=intType), dimension(:), pointer :: colind, rowptr
        real(kind=realtype), dimension(:), pointer :: data
        integer(Kind=intType), dimension(:), pointer :: assignedproc
 
        ! Flag if this matrix is allocated
        logical :: allocated = .false.
 
    end type csrmatrix
 
    ! This derived type contains sufficient information to perfom ADT
    ! donor searches.
    type oversetblock
 
        ! Sizes for the block
        integer(kind=intType) :: il, jl, kl
 
        ! Cluster this block belongs to
        integer(kind=intType) :: cluster
 
        ! This is the cell volume of the donor
        real(kind=realtype), dimension(:, :), pointer :: qualdonor
 
        ! Connectivity for the ADT
        integer(kind=intType), dimension(:, :), pointer :: hexaconn
 
        ! Coordinates for the ADT
        real(kind=realtype), dimension(:, :), pointer :: xadt
 
        ! Copy of global cell
        integer(kind=intType), dimension(:, :, :), pointer :: globalcell
 
        ! Whether or not a cell is "near" a wall
        integer(kind=intType), dimension(:, :, :), pointer :: nearwall
 
        ! Whether or not a cell is not possible to be a donor. Ie. a forceRecv
        integer(kind=intType), dimension(:, :, :), pointer :: invaliddonor
 
        ! Minimum volume for this block
        real(kind=realtype) :: minvol
 
        ! The ADT for this block
        type(adttype) :: adt
 
        ! The processor for this block
        integer(kind=intType) :: proc
 
        ! And the local block index
        integer(kind=intType) :: block
 
        ! Buffer space for sending/receiving the block
        real(kind=realtype), dimension(:), allocatable :: rbuffer
        integer(kind=intType), dimension(:), allocatable :: ibuffer
 
        ! Flag if this block got allocated
        logical :: allocated = .false.
 
        ! Flag if the real/int Buffers are ready after receiving info
        logical :: realbufferready = .false.
        logical :: intbufferready = .false.
 
    end type oversetblock
 
    type oversetfringe
 
        ! The processor where this set of fringes came from
        integer(kind=intType) :: proc
 
        ! The block number of these fringes on processor 'proc'
        integer(kind=intType) :: block
 
        ! Sizes
        integer(kind=intType) :: il, jl, kl, nx, ny, nz
 
        ! Cluster this set of fringes belongs to
        integer(kind=intType) :: cluster
 
        ! Buffer space for sending/receiving the fringes
        real(kind=realtype), dimension(:), allocatable :: rbuffer
        integer(kind=intType), dimension(:), allocatable :: ibuffer
 
        ! These are the coordinates of the points we are searching for
        real(kind=realtype), dimension(:, :), allocatable :: x
 
        ! These are the coordinate of its wall surface pt if applicable
        real(kind=realtype), dimension(:, :), allocatable :: xseed
 
        ! These are the indices of the wall surfaces
        integer(kind=intType), dimension(:), allocatable :: wallind
 
        ! Flag specifying if this cell is next to a wall or not. 1 for
        ! next to wall, 0 otherwise.
        integer(kind=intType), dimension(:), allocatable :: iswall
 
        ! This is where we will store all the potential donors that have
        ! been found for this set of fringes
        integer(kind=intType) :: ndonor = 0
        integer(kind=intType), dimension(:, :), pointer :: fringeintbuffer => null()
        real(kind=realtype), dimension(:, :), pointer :: fringerealbuffer => null()
 
        ! Flag if this set of fringes got allocated
        logical :: allocated = .false.
 
        ! Flag if the real/int Buffers are ready after receiving info
        logical :: realbufferready = .false.
        logical :: intbufferready = .false.
 
        ! The number of actual fringes that need to communicated.
        integer(kind=intType) :: fringereturnsize = 0
 
    end type oversetfringe
 
    type oversetwall
 
        ! Sizes
        integer(kind=intType) :: il, jl, kl
        integer(kind=intType) :: nnodes = 0
        integer(kind=intType) :: ncells = 0
        integer(kind=intType) :: maxcells = 0
        integer(kind=intType) :: cluster = 0
 
        ! Buffer space for sending/receiving the fringes
        real(kind=realtype), dimension(:), allocatable :: rbuffer
        integer(kind=intType), dimension(:), allocatable :: ibuffer
 
        ! Surface nodes used to build the tree:
        real(kind=realtype), dimension(:, :), pointer :: x => null()
 
        ! Only primal mesh cell centers for dual mesh. Only allocated as
        ! needed.
        real(kind=realtype), dimension(:, :), pointer :: xprimalcen => null()
 
        ! Surface nonal used for determining if point is "underneath" the
        ! surface.
        real(kind=realtype), dimension(:, :), pointer :: norm => null()
 
        ! Local estimate of surface error
        real(kind=realtype), dimension(:), pointer :: delta => null()
 
        ! Connectivity for the surface
        integer(kind=intType), dimension(:, :), pointer :: conn => null()
 
        ! ind: Global node index for nodes
        integer(kind=intType), dimension(:), pointer :: ind => null()
 
        ! indPrimal: Global node index. Only temporarly used to store
        ! index for strictly primal cells on dual mesh.
        integer(kind=intType), dimension(:), pointer :: indprimal => null()
 
        ! indCell: Global cell index for wall cells
        integer(kind=intType), dimension(:), pointer :: indcell
 
        ! Blanking values for Nodes
        integer(kind=intType), dimension(:), allocatable :: iblank
        integer(kind=intType), dimension(:), allocatable :: cellptr
 
        ! Node to element array
        integer(kind=intType), dimension(:, :), allocatable :: nte
 
        ! The ADT for this block's wall(s)
        type(adttype) :: adt
 
        ! This KDTree for this block's wall
        type(kdtree2), pointer :: tree
 
        ! Flag if the real/int Buffers are ready after receiving info
        logical :: realbufferready = .false.
        logical :: intbufferready = .false.
 
        ! Flag if this wall got allocated
        logical :: allocated = .false.
 
    end type oversetwall
 
    type oversetstring
        ! This is a generic type defining a string list. It may be used
        ! as both a "parent" or a "child".
 
        ! Sizes
        integer(kind=intType) :: nnodes, nelems
 
        ! My String's index
        integer(kind=intType) :: myid
 
        ! Whether or no this string is periodic
        logical :: isperiodic = .false.
 
        ! Whether or no this string is a pocket
        logical :: ispocket = .false.
 
        ! --------------------------------------------------------------------
        ! Node Data: The actual physical node locations, unit surface normal,
        ! perpNormal and mesh size. x is from index 1:3, normal from 4:6,
        ! perpNormal form 7:9 and h is index 10. This pointer gets allocated.
        real(kind=realtype), dimension(:, :), pointer :: nodedata => null()
 
        ! Pointer for physical node location. Points to nodeData
        real(kind=realtype), dimension(:, :), pointer :: x => null()
 
        ! Pointer for nodal unit normal. Points to nodeData
        real(kind=realtype), dimension(:, :), pointer :: norm => null()
 
        ! Pointer for nodal unit perpendicual in-plane normal. Points to nodeData
        real(kind=realtype), dimension(:, :), pointer :: perpnorm => null()
 
        ! Pointer for nodal element size. Points to nodeData
        real(kind=realtype), dimension(:), pointer :: h => null()
 
        ! --------------------------------------------------------------------
        ! Integer Node Data: This stores the global node index (into
        ! xVec), the cluster ID of the node, and the family ID of the node
        integer(kind=intType), dimension(:, :), pointer :: intnodedata => null()
 
        ! The orignal nodal index. Size nNodes. Pointer into intNodeData
        integer(kind=intType), dimension(:), pointer :: ind => null()
 
        ! The cluster the node came from. Pointer into intNodeData
        integer(kind=intType), dimension(:), pointer :: cluster => null()
 
        ! The family the node came from. Pointer into intNodeData
        integer(kind=intType), dimension(:), pointer :: family => null()
        ! --------------------------------------------------------------------
 
        ! The connectivity of the nodes forming 1D bar elements. Size (2, nElems)
        integer(kind=intType), dimension(:, :), pointer :: conn => null()
 
        ! The index of my node numbers on my parent
        integer(kind=intType), dimension(:), pointer :: pnodes => null()
 
        ! The index of my elements numbers on my parent
        integer(kind=intType), dimension(:), pointer :: pelems => null()
 
        ! The string ID and index of my nodes on a split substing
        integer(kind=intType), dimension(:, :), pointer :: cnodes => null()
 
        ! The cloest string ID of each node *AND* the node index on the
        ! other string. Size (2, nNodes)
        integer(kind=intType), dimension(:, :), pointer :: otherid => null()
 
        ! The inverse of the connectivity node to elem array. Size (5,
        ! nNodes). First index is the number of elements, other 4 entries
        ! are the up to 4 possible element neighbours.
        integer(kind=intType), dimension(:, :), pointer :: nte => null()
 
        ! Two buffer used for storing element indices while creating
        ! chains. Size (2, nElem)
        integer(kind=intType), dimension(:, :), pointer :: substr => null()
 
        ! The sizes of the two substrings
        integer(kind=intType), dimension(2) :: nsubstr
 
        ! A array to keep track of the number of elements
        ! "consumed" during chain searches or during zipping.
        integer(kind=intType), dimension(:), pointer :: elemused => null()
 
        ! Array to keep track of nodes used to contruct string pairs for
        ! crossZipping.
        integer(kind=intType), dimension(:), pointer :: xzipnodeused => null()
 
        ! The KD tree for this string for performing fast seaches.
        !type(tree_master_record), pointer :: tree
        type(kdtree2), pointer :: tree => null()
 
        ! Pointer to the parent string
        type(oversetstring), pointer :: p => null()
 
        ! Pointer for next string for a linked list
        type(oversetstring), pointer :: next => null()
 
        ! List of all all directed edges.
        type(oversetedge), pointer, dimension(:) :: edges => null()
 
        ! nEdges: The number of new edges added due to triangles.
        integer(kind=intTYpe) :: nedges
 
        ! List of all computed triangles
        integer(kind=intType), dimension(:, :), pointer :: tris => null()
 
        ! Number of trianges
        integer(kind=intType) :: ntris
 
        ! surfCellID(1:nTris)
        ! Global cellID of the primal cell containing the triangle centroid
        integer(kind=intType), dimension(:), pointer :: surfcellid
 
    end type oversetstring
 
    type oversetedge
        ! Simple data structure representing a directed edge from n1->n2
        integer(kind=intType) :: n1, n2
    end type oversetedge
 
    interface operator(<=)
        module procedure lessequaledgetype
    end interface operator(<=)
 
    interface operator(<)
        module procedure lessedgetype
    end interface operator(<)
 
    type pocketedge
        ! Simple data structure representing a directed edge from n1->n2
        ! Similar to oversetEdge, but introduced to do different type
        ! of sort on edges
        integer(kind=intType) :: n1, n2
    end type pocketedge
 
    interface operator(<=)
        module procedure lessequalpocketedgen2
    end interface operator(<=)
 
    interface operator(<)
        module procedure lesspocketedgen2
    end interface operator(<)
 
    interface operator(==)
        module procedure equalpocketedgen2
    end interface operator(==)
 
    type zippermesh
 
        ! Data required for zipper mesh surface integer for a particular BCGroup
        integer(kind=intType), dimension(:, :), allocatable :: conn
        integer(kind=intType), dimension(:), allocatable :: fam, indices
        logical :: allocated = .false.
#ifndef USE_TAPENADE
        vecscatter :: scatter
        vec :: localval
#endif
    end type zippermesh
 
    ! This is the flattened list of the fringes next to the wall that we
    !  have actually found donors for.
    ! tmpFringePtr is only used if we need to realloc.
    type(fringetype), dimension(:), pointer :: localwallfringes, wallfringes, tmpfringeptr
    integer(kind=intType) :: nlocalwallfringe, nwallfringe
 
    ! These are the master overlap matrices
    type(csrmatrix), dimension(:, :), allocatable, target :: overlapmatrix
 
    ! Some additional helper stuff
    integer(kind=intType), dimension(:), allocatable :: ndomproc, cumdomproc
    integer(kind=intType) :: ndomtotal
    integer(kind=intType) :: nclusters
    integer(kind=intType), dimension(:), allocatable :: clusters
    real(kind=realtype), dimension(:), allocatable :: clusterareas
    real(kind=realtype), dimension(:), allocatable :: clustermarchdist
    type(oversetwall), dimension(:), allocatable, target :: clusterwalls
 
    ! Flag specifying if overset is present in mesh
    logical :: oversetpresent
 
    ! Zipper meshes
    type(zippermesh), dimension(nFamExchange), target :: zippermeshes
 
    ! Static arrays for doing timings
    real(kind=realtype), dimension(iTotal) :: tstart
    real(kind=realtype), dimension(iTotal) :: oversettimes
 
contains
    ! ==============================
    ! Operator overloading functions
    ! ==============================
    logical function lessequaledgetype(e1, e2)
        !
        !   lessEqualEdgeType returns .True. if e1<=e2 and .False. otherwise.
        !   Compared on the directed edge node indices n1 and n2.
        !   First compare wrt averaged node indices data. If same averaged
        !   node data, then compare wrt increasing or decreasing node indices.
 
        implicit none
 
        ! Input
        type(oversetedge), intent(in) :: e1, e2
 
        ! Local variables
        integer(kind=intType) :: nsum1, nsum2, ndiff1, ndiff2
 
        ! Compare the averaged (or just sum of) node indices values.
        ! Positive sign for increasing order of node indices.
        nsum1 = e1%n1 + e1%n2
        nsum2 = e2%n1 + e2%n2
        ndiff1 = e1%n2 - e1%n1
        ndiff2 = e2%n2 - e2%n1
 
        ! Compare based on averaged node indices values
        if (abs(nsum1) < abs(nsum2)) then
            lessequaledgetype = .true.
            return
        else if (abs(nsum1) > abs(nsum2)) then
            lessequaledgetype = .false.
            return
        end if
 
        if (abs(nsum1) /= abs(nsum2)) &
            stop ' *** Error in lessEqualEdgeType ***'
 
        ! Compare based on edge nodes difference
        if (abs(ndiff1) < abs(ndiff2)) then
            lessequaledgetype = .true.
            return
        else if (abs(ndiff1) > abs(ndiff2)) then
            lessequaledgetype = .false.
            return
        end if
 
        ! here abs(ndiff1) == abs(ndiff2) and
        ! abs(nsum1)== abs(nsum2), so same edge
        if (ndiff1 < ndiff2) then
            lessequaledgetype = .true.
            return
        else if (ndiff1 > ndiff2) then
            lessequaledgetype = .false.
            return
        end if
 
        ! here ndiff1 == ndiff2, hence .True.
 
        lessequaledgetype = .true.
 
    end function lessequaledgetype
 
    ! ---------------------------------
    logical function lessedgetype(e1, e2)
        !
        !   lessEdgeType returns .True. if e1<e2 and .False. otherwise.
        !   Compared on the directed edge node indices n1 and n2.
        !   First compare wrt averaged node indices data. If same averaged
        !   node data, then compare wrt increasing or decreasing node indices.
 
        implicit none
 
        ! Input
        type(oversetedge), intent(in) :: e1, e2
 
        ! Local variables
        integer(kind=intType) :: nsum1, nsum2, ndiff1, ndiff2
 
        ! Compare the averaged (or just sum of) node indices values.
        ! Positive sign for increasing order of node indices.
        nsum1 = e1%n1 + e1%n2
        nsum2 = e2%n1 + e2%n2
        ndiff1 = e1%n2 - e1%n1
        ndiff2 = e2%n2 - e2%n1
 
        ! Compare based on averaged node indices values
        if (abs(nsum1) < abs(nsum2)) then
            lessedgetype = .true.
            return
        else if (abs(nsum1) > abs(nsum2)) then
            lessedgetype = .false.
            return
        end if
 
        if (abs(nsum1) /= abs(nsum2)) &
            stop ' *** Error in lessEdgeType ***'
 
        ! Compare based on edge nodes difference
        if (abs(ndiff1) < abs(ndiff2)) then
            lessedgetype = .true.
            return
        else if (abs(ndiff1) > abs(ndiff2)) then
            lessedgetype = .false.
            return
        end if
 
        ! here abs(ndiff1) == abs(ndiff2) and
        ! abs(nsum1)== abs(nsum2), so same edge
        if (ndiff1 < ndiff2) then
            lessedgetype = .true.
            return
        else if (ndiff1 > ndiff2) then
            lessedgetype = .false.
            return
        end if
 
        ! here ndiff1 == ndiff2, hence .False.
 
        lessedgetype = .false.
 
    end function lessedgetype
 
    logical function lessequalpocketedgen2(e1, e2)
        !
        !   lessEqualPocketEdgeN2 returns .True. if e1%n2<=e2%n2 and .False.
        !   otherwise. Compared on the directed edge node indices n1 and n2.
 
        implicit none
 
        ! Input
        type(pocketedge), intent(in) :: e1, e2
 
        if (e1%n2 < e2%n2) then
            lessequalpocketedgen2 = .true.
            return
        else if (e1%n2 > e2%n2) then
            lessequalpocketedgen2 = .false.
            return
        end if
 
        ! Here e1%n2==e2%n2, so edges are equal, hence .True.
        lessequalpocketedgen2 = .true.
 
    end function lessequalpocketedgen2
 
    logical function lesspocketedgen2(e1, e2)
        !
        !   lessPocketEdgeN2 returns .True. if e1%N2<e2%N2 and .False.
        !   otherwise. Compared on the directed edge node indices n1 and n2.
 
        implicit none
 
        ! Input
        type(pocketedge), intent(in) :: e1, e2
 
        if (e1%N2 < e2%N2) then
            lesspocketedgen2 = .true.
            return
        else if (e1%N2 > e2%N2) then
            lesspocketedgen2 = .false.
            return
        end if
 
        ! Here e1%N2==e2%n2, so edges are equal, hence .False.
        lesspocketedgen2 = .false.
 
    end function lesspocketedgen2
 
    logical function equalpocketedgen2(e1, e2)
        !
        !   EqualPocketEdgeN2 returns .True. if e1%N2==e2%N2 and .False.
        !   otherwise. Compared on the directed edge node indices n1 and n2.
 
        implicit none
 
        ! Input
        type(pocketedge), intent(in) :: e1, e2
 
        if (e1%N2 == e2%N2) then
            equalpocketedgen2 = .true.
            return
        else
            equalpocketedgen2 = .false.
            return
        end if
 
    end function equalpocketedgen2
end module oversetdata