! Copyright (c) 2012-2014 Simon Zimny <s.zimny@grs-sim.de> ! Copyright (c) 2012-2013 Manuel Hasert <m.hasert@grs-sim.de> ! Copyright (c) 2012-2016, 2019, 2021 Kannan Masilamani <kannan.masilamani@uni-siegen.de> ! Copyright (c) 2013 Melven Zoellner <yameta@freenet.de> ! Copyright (c) 2013-2019 Harald Klimach <harald.klimach@uni-siegen.de> ! Copyright (c) 2013, 2016 Nikhil Anand <nikhil.anand@uni-siegen.de> ! Copyright (c) 2014 Kartik Jain <kartik.jain@uni-siegen.de> ! Copyright (c) 2014-2016 Jiaxing Qi <jiaxing.qi@uni-siegen.de> ! Copyright (c) 2015-2016 Peter Vitt <peter.vitt2@uni-siegen.de> ! Copyright (c) 2016 Tobias Schneider <tobias1.schneider@student.uni-siegen.de> ! Copyright (c) 2016 Verena Krupp <verena.krupp@uni-siegen.de> ! Copyright (c) 2017 Neda Ebrahimi Pour <neda.epour@uni-siegen.de> ! Copyright (c) 2018 Raphael Haupt <Raphael.Haupt@student.uni-siegen.de> ! ! 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. ! ! 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. ! ****************************************************************************** ! !> This module contains subroutines needed to create and update subTrees from !! corresponding global trees. !! !! The subroutine tem_create_subTree_of creates a subTree from a global tree !! using an array of shapes. This functionality is used when certain operations !! (e.g. tracking, source terms) shall be executed only on a subset of the fluid !! tree (or another tree). !! module tem_subTree_module ! include treelm modules use mpi use env_module, only: long_k, globalMaxLevels, pathLen, & & long_k_mpi, labelLen use treelmesh_module, only: treelmesh_type use tem_bc_prop_module, only: tem_bc_prop_type use tem_aux_module, only: tem_abort use tem_tools_module, only: tem_horizontalSpacer use tem_debug_module, only: tem_debug_type, main_debug use tem_dyn_array_module, only: dyn_intArray_type, init, append, destroy use tem_grow_array_module, only: grw_intArray_type, init, append, destroy use tem_property_module, only: gather_property, prp_hasQVal use tem_prophead_module, only: tem_prop_countelems use tem_shape_module, only: tem_shape_type, tem_global_shape, & & tem_geometrical_shape, tem_property_shape,& & tem_local_shape, tem_boundary_shape, & & tem_level_shape, & & tem_shape_subTreeFromGeomInters use tem_geometry_module, only: tem_setEffBoundingBox, tem_baryOfId use tem_topology_module, only: tem_levelOf, tem_firstIdAtLevel, & & tem_lastIdAtLevel use tem_param_module, only: qQQQ, qDir use tem_construction_module, only: tem_levelDesc_type use tem_subTree_type_module, only: tem_subTree_type, & & tem_dump_subTree use tem_stencil_module, only: tem_stencilHeader_type use tem_logging_module, only: logUnit, tem_log, tem_toStr use tem_pointData_module, only: tem_grwPoints_type, destroy implicit none private public :: tem_updatePropertyBits public :: tem_create_subTree_of public :: tem_write_debugMesh public :: tem_create_tree_from_sub public :: tem_subTree_from !> interface for copying the property bits from the global tree or the !! level descriptor interface tem_copyPropertyBits module procedure tem_copyPropertyBitsFromLevelDesc module procedure tem_copyPropertyBitsFromTree end interface contains ! ****************************************************************************** ! !> Update the property of subTree with the ones from inTree, if something !! changed update the logical meshChange. !! subroutine tem_updatePropertyBits( inTree, subTree ) ! --------------------------------------------------------------------------- !> tree to get information from type(treelmesh_type), intent(in) :: inTree !> tree to pass information to type(tem_subTree_type), intent(inout) :: subTree ! --------------------------------------------------------------------------- integer :: iElem integer :: elemPos ! --------------------------------------------------------------------------- do iElem = 1, subTree%nElems ! get the position in the global tree elemPos = subTree%map2global(iElem) ! test the individual property bits if ( inTree%ElemPropertyBits(elempos) & & /= subTree%ElemPropertyBits(iElem) ) then ! update the propertyBits in the subTree subTree%ElemPropertyBits(iElem) = inTree%ElemPropertyBits(elemPos) subTree%global%meshChange = .true. end if end do end subroutine tem_updatePropertyBits ! ****************************************************************************** ! ! ****************************************************************************** ! !> Copy the properties of the level decriptor to the ones in subTree. !! subroutine tem_copyPropertyBitsFromLevelDesc( levelDesc, subTree ) ! --------------------------------------------------------------------------- !> level descriptor including all elements (fluid, ghost, halo) type(tem_levelDesc_type ),intent(in) :: levelDesc(:) !> tree to pass information to type(tem_subTree_type), intent(inout) :: subTree ! --------------------------------------------------------------------------- integer :: iLevel integer :: nElems integer :: elemCounter ! --------------------------------------------------------------------------- allocate( subTree%ElemPropertyBits( subTree%nElems )) elemCounter = 0 do iLevel = 1, size(levelDesc) nElems = size(levelDesc(iLevel)%property(:)) subTree%ElemPropertyBits(elemCounter+1 : elemCounter+nElems) = & & levelDesc(iLevel)%property(:) elemCounter = elemCounter + nElems end do end subroutine tem_copyPropertyBitsFromLevelDesc ! ****************************************************************************** ! ! ****************************************************************************** ! !> Copy the properties of inTree to the ones in subTree. !! subroutine tem_copyPropertyBitsFromTree( inTree, subTree ) ! --------------------------------------------------------------------------- !> tree to get information from type(treelmesh_type), intent(in) :: inTree !> tree to pass information to type(tem_subTree_type), intent(inout) :: subTree ! --------------------------------------------------------------------------- integer :: iElem integer :: elemPos integer :: iProp ! --------------------------------------------------------------------------- ! allocate the bitfield allocate(subTree%elemPropertyBits(subTree%nElems)) ! initialize the bitfield with 0 subTree%elemPropertyBits = 0 do iElem = 1, subTree%nElems ! get the position in the global tree elemPos = subTree%map2global(iElem) ! copy the propertyBits in the subTree subTree%ElemPropertyBits(iElem) = inTree%ElemPropertyBits(elemPos) end do ! update the property information allocate(subTree%property(subTree%global%nProperties)) do iProp=1,subTree%global%nProperties call gather_Property( Property = subTree%Property(iProp), & & Header = subTree%global%Property(iProp), & & BitField = subTree%ElemPropertyBits, & & comm = subTree%global%comm ) end do end subroutine tem_copyPropertyBitsFromTree ! ****************************************************************************** ! ! ****************************************************************************** ! !> This subroutine creates a subtree based on a provided map or list of !! treeIDs (in case a local shape is used) to the corresponding tree. !! Only processes in comm will be involved. subroutine tem_subTree_from( me, map2global, treeID, comm, dirname, grwPnts ) ! --------------------------------------------------------------------------- !> subTree to be created from list of elements (map2global) type(tem_subTree_type), intent(inout) :: me !> position of the treeID in the global treeID list integer, optional, intent(in) :: map2global(:) !> list of treeIDs only use this in case a local shape is set integer(kind=long_k), optional, intent(in) :: treeID(:) !> mpi communicator to use, defaults to the one in me%global%comm if not !! specified integer, intent(in), optional :: comm !> directory to store the mesh in. is taken to be me%global%dirname if not !! specified character(len=*), intent(in), optional :: dirname !> array of point vaues that neeeds to be stored in the subtree type(tem_grwPoints_type), intent(in), optional :: grwPnts ! --------------------------------------------------------------------------- integer :: commloc integer(kind=long_k) :: offset, nElems, nPoints integer :: ierror integer :: nElemsList ! mpi variables integer :: commsize, rank, iPnt ! --------------------------------------------------------------------------- if (present(dirname)) then me%global%dirname = dirname end if if (present(comm)) then commloc = comm else commloc = me%global%comm end if ! allocate and store the points for get point tracking if ! applicable if (present(grwPnts)) then me%nPoints = grwPnts%coordX%nVals allocate(me%points(me%nPoints,3)) do iPnt = 1, me%nPoints me%points(:, 1) = grwPnts%coordX%val(1:me%nPoints) me%points(:, 2) = grwPnts%coordY%val(1:me%nPoints) me%points(:, 3) = grwPnts%coordZ%val(1:me%nPoints) end do end if ! determine number of elements in the provided map2global nElemsList = 0 if(present(map2global) .and. .not. present(treeID))then nElemsList = size(map2global) allocate(me%map2global( nElemsList )) ! copy the map2global me%map2global = map2global else if(present(treeID) .and. .not. present(map2global))then nElemsList = size(treeID) allocate(me%treeID( nElemsList )) ! copy the list of treeID me%treeID = treeID ! in case a treeID is set we have to set useLocalMesh to true me%useLocalMesh = .true. else write(logUnit(1),*)'error: either none or both of treeID or map2global ' write(logUnit(1),*)' are passed to tem_subTree_from. '// & & ' only one of them is allowed' call tem_abort() end if ! get size of communicator and my rank in it call mpi_comm_size(commloc, commsize, ierror) call mpi_comm_rank(commloc, rank, ierror) me%global%nparts = commsize me%global%mypart = rank ! assign treeID list to the new mesh. sort first offset = 0 nElems = int(nElemsList, kind=long_k) ! determine offset for each process (get nElems of ranks before) call MPI_Exscan( nElems, offset, 1, long_k_mpi, mpi_sum, commloc, ierror ) me%elemOffset = offset me%nElems = nElemsList ! last process offset + its nElems is total number of elems nElems = offset + nElems ! distribute among new communicator call MPI_Bcast( nElems, 1, long_k_mpi, commsize-1, commloc, ierror ) ! store in global number of elements me%global%nElems = nElems if (present(grwPnts)) then nPoints = me%nPoints call MPI_Reduce(nPoints, me%glob_nPoints, 1, long_k_mpi, mpi_sum, & & 0, commloc, ierror) end if end subroutine tem_subTree_from ! ****************************************************************************** ! ! ****************************************************************************** ! !> Create a subTree from a given inTree and an array of shapes. !! !! All elements to the inTree%treeID are collected into the subTree%map2global !! array and a sub-communicator with participating processes is created. !! Additionally how many and which elements exist on my local process and are !! requested from the shapes is identified. !! subroutine tem_create_subTree_of( inTree, subTree, inShape, levelDesc, & & bc_prop, storePnts, stencil, prefix ) ! --------------------------------------------------------------------------- !> Global mesh from which the elements are identified and then stored to type(treelmesh_type), intent(in) :: inTree !> new mesh type(tem_subTree_type), intent(out) :: subTree !> shape objects on which to work type(tem_shape_type ),intent(in) :: inShape(:) !> optional level descriptor needed for local shape type(tem_levelDesc_type ), optional, intent(in) :: levelDesc(:) !> bc property which is used to identify elements belong to certain BCs type( tem_bc_prop_type ), optional, intent(in) :: bc_prop !> To store space points in subTree logical, optional, intent(in) :: storePnts !> stencil used to find bcID on certain links type( tem_stencilHeader_type ), optional, intent(in) :: stencil !> prefix for the subTree character(len=*), optional, intent(in) :: prefix ! --------------------------------------------------------------------------- ! local growing array for the map to the global tree type(dyn_intArray_type) :: local_map2global type(tem_grwPoints_type) :: local_grwPnts integer, allocatable :: map2globalTemp(:) integer(kind=long_k), allocatable :: tmp_treeID(:) integer :: iShape, iElem ! mpi variables for parallel writing of the log-wise trees integer :: color, iError ! if the local rank has part of the tracking object logical :: participateInMesh, globalParticipateInMesh ! if any element be found for this iShape logical :: foundAny integer :: local_countElems( globalMaxLevels ) integer :: local_countPnts logical :: local_storePnts ! -------------------------------------------------------------------------- if (present(storePnts)) then local_storePnts = storePnts else local_storePnts = .false. end if subTree%useGlobalMesh = .false. subTree%useLocalMesh = .false. subTree%created_new_comm = .false. ! Prepare output structure ! init here the sub-tree mesh subTree%global = inTree%global ! initialize the number of elements on this partition with 0 subTree%nElems = 0 subTree%global%nElems = 0 if ( present( prefix )) subTree%global%dirname = trim(prefix) ! Reset counter for actual elements for this shape object ! @todo: indeed only sum(local_countElems) is used later. ! thus no need to calculate countElems in following routines. local_countElems = 0 local_countPnts = 0 call init( local_map2global, length = 128) do iShape = 1, size( inShape ) select case( inShape(iShape)%shapeID ) ! Use elements belonging to a canonical n-dimensional shape case( tem_geometrical_shape ) call tem_log(5, 'iShape '//trim(tem_toStr(iShape)) & & //' is geometrical shape') call tem_shape_subTreeFromGeomInters( me = inShape(iShape), & & inTree = inTree, & & countElems = local_countElems, & & countPoints = local_countPnts, & & grwPnts = local_grwPnts, & & storePnts = local_storePnts, & & map2global = local_map2global ) ! Use all elements for this shape case( tem_global_shape ) call tem_log(5, 'iShape '//trim(tem_toStr(iShape))//' is a global shape') ! Global shape means, all elements are part of the shape. ! so nothing to do, just continue with the next shape local_countElems = 1 subTree%useGlobalMesh = .true. exit ! Only use elements with a certain property case( tem_property_shape ) call tem_log(5, 'iShape '//trim(tem_toStr(iShape))//' is a property shape') call tem_shape_initPropElements( inShape(iShape)%propBits, inTree, & & local_countElems, local_map2global ) ! Only use elements belong to certain boundaries case( tem_boundary_shape ) call tem_log(5, 'iShape '//trim(tem_toStr(iShape))//' is a boundary shape') if (present(bc_prop) .and. present(stencil)) then call tem_shape_findElemByBCLabels( bcLabels = inShape(iShape) & & %bcLabels, & & bc_prop = bc_prop, & & foundAny = foundAny, & & map2global = local_map2Global, & & bcIDs = subTree%bc_ID, & & stencil = stencil ) if ( foundAny ) local_countElems = 1 else call tem_abort('In tem_create_subTree_of: Stencil or bc_prop not' & & //'passed') end if ! elements on the local partition case( tem_local_shape ) call tem_log(5, 'iShape '//trim(tem_toStr(iShape))//' is a local shape') if(present(levelDesc)) then subTree%useLocalMesh = .true. call tem_shape_initLocal( levelDesc, tmp_treeID ) else write(logUnit(1),*)'ERROR:levelDesc is not passed to '// & & 'tem_creat_subTree_of' write(logUnit(1),*)' levelDesc is required to create '// & & 'subTree for tem_local_shape' call tem_abort() end if local_countElems = 1 ! update the communicator subTree%global%comm = mpi_comm_self exit case( tem_level_shape ) call tem_log(5, 'iShape '//trim(tem_toStr(iShape))//' is a level shape') call tem_shape_initByLevels( inTree, inShape(iShape)%minLevel, & & inShape(iShape)%maxLevel, & & local_countElems, local_map2global ) end select ! shapeID end do ! iShape if( .not. subTree%useGlobalMesh ) then subTree%global%predefined = '' end if ! Exchange information about the mesh ! mpi_exscan: track( iLog )%tree%offset ! Now create a new communicator with ranks that own parts of the sub-mesh ! This is done with the color participateInMesh = ( sum( local_countElems ) > 0 .or. local_countPnts > 0 ) ! abort if no elements are found for the given shape call MPI_ALLREDUCE( participateInMesh, globalParticipateInMesh, 1, & & MPI_LOGICAL, MPI_LOR, inTree%global%comm, iError ) if (.not. globalParticipateInMesh) then write(logUnit(1),*) 'Error: No elements found for '// & & 'defined shape '//trim(prefix)//' in the mesh' call tem_abort() end if if( participateInMesh ) then color = 1 else ! Not participating, set color to undefined color = MPI_UNDEFINED end if ! Input rank as a key to keep sorting as in global communicator ! reorder process ranks into tracking communicator if( .not. subTree%useLocalMesh ) then call mpi_comm_split( inTree%global%comm, color, inTree%global%myPart, & & subTree%global%comm, iError ) if( iError .ne. mpi_success ) then write(logUnit(1),*)' There was an error splitting the communicator' write(logUnit(1),*)' for the tracking subset. Aborting...' call tem_abort() end if subTree%created_new_comm = .true. end if if( participateInMesh ) then call mpi_comm_size( subTree%global%comm, subTree%global%nParts, iError ) call mpi_comm_rank( subTree%global%comm, subTree%global%myPart, iError ) if( .not. subTree%useGlobalmesh) then if( subTree%useLocalMesh ) then call tem_subTree_from( me = subTree, & & treeID = tmp_treeID ) ! update the property list of the newly defined tree according to the ! the local shape (simply copy the bits) call tem_copyPropertyBits( levelDesc = levelDesc, subTree = subTree ) else ! NOT useLocalMesh ! Copy the created dynamic local_tIDlist to a long integer array allocate( map2globalTemp( local_map2global%nVals )) ! Assign list of links to the temporary list, first sorting do iElem = 1, local_map2global%nVals map2globalTemp( iElem ) = local_map2global%val( & & local_map2global%sorted( iElem )) end do call tem_subTree_from( me = subTree, & & map2global = map2globalTemp, & & grwPnts = local_grwPnts ) call destroy( me = local_grwPnts ) ! copy property bits from the global tree call tem_copyPropertyBits( inTree = inTree, subTree = subTree ) deallocate( map2globalTemp ) end if ! useLocalMesh with level desc ! Calculate the effective extents of the tracking sub-tree bounding box call tem_setEffBoundingBox( subTree, inTree ) else ! useGlobalMesh ! assign the global properties to the subTree subTree%nElems = inTree%nElems subTree%global%nElems = inTree%global%nElems end if ! not global mesh? if (subTree%global%myPart == 0) then if (subTree%global%nElems == 1 .and. allocated(subTree%map2global) ) then write(*,"(A,3F12.6)") 'Barycenter of tracking element: ', & & tem_baryOfId(inTree, inTree%treeID(subTree%map2global(1))) end if end if !subTree root else ! NOT participate in mesh subTree%nElems = 0 subTree%global%nElems = 0 allocate( subTree%map2global(0) ) subTree%nPoints = 0 subTree%glob_nPoints = 0_long_k allocate( subTree%points(0,0) ) end if ! Participate in mesh? write(logUnit(5),"(A,I0)") 'Found corresponding nElems: ', & & subTree%global%nElems call destroy( me = local_map2global ) call tem_horizontalSpacer( fUnit= logUnit(5) ) end subroutine tem_create_subTree_of ! ****************************************************************************** ! ! ****************************************************************************** ! !> Create newtree out of intree by restricting to the elements of subtree. !! !! The new mesh will have no properties subroutine tem_create_tree_from_sub(intree, subtree, newtree, keep_props) !> The tree on which the subtree is defined. type(treelmesh_type), intent(in) :: intree !> Subtree describing the part of the mesh to create a new mesh from. type(tem_subtree_type), intent(in) :: subtree !> Resulting new tree with the elements selected by subtree from newtree. type(treelmesh_type), intent(out) :: newtree !> Flag to indicate whether to keep properties from intree also in newtree. !! !! If this is true, the properties will be copied from the intree to the !! newtree. An actual copy is done, as we can not rely on the pointer !! targets in intree to exist further on. !! Default is `.false.`, which means all properties will be dropped and !! newtree will have no properties at all. logical, optional, intent(in) :: keep_props logical :: withprop integer(kind=long_k) :: nNewElems integer :: iProp integer :: iError withprop = .false. if (present(keep_props)) withprop = keep_props newtree%global%maxlevel = intree%global%maxlevel newtree%global%minlevel = intree%global%minlevel newtree%global%origin = intree%global%origin newtree%global%BoundingCubeLength = intree%global%BoundingCubeLength if (subtree%useGlobalMesh) then ! Copy complete tree, but ignore properties. newtree = intree nullify(newtree%global%property) nullify(newtree%property) else newtree%nelems = subtree%nElems nullify(newtree%global%property) nullify(newtree%property) newtree%global%comm = subtree%global%comm newtree%global%nparts = subtree%global%nparts newtree%global%myPart = subtree%global%myPart allocate(newtree%treeID(newtree%nelems)) allocate(newtree%ElemPropertyBits(newtree%nelems)) allocate(newtree%Part_First(newtree%global%nparts)) allocate(newtree%Part_Last(newtree%global%nparts)) newtree%treeID = intree%treeID(subtree%map2global) if (withprop) then newtree%ElemPropertyBits = intree%ElemPropertyBits(subtree%map2global) else newtree%ElemPropertyBits = 0_long_k end if nNewElems = int(newtree%nElems, kind=long_k) ! Overall number of elements in the new mesh and offsets. call MPI_Exscan(nNewelems, newtree%ElemOffset, 1, long_k_mpi, & & MPI_SUM, newtree%global%comm, iError ) newtree%global%nElems = newtree%ElemOffset+nNewElems call MPI_Bcast( newtree%global%nElems, 1, long_k_mpi, & & newtree%global%nParts-1, & & newtree%global%comm, iError ) call MPI_Allgather( newtree%treeID(1), 1, long_k_mpi, & & newtree%Part_First, 1, long_k_mpi, & & newtree%global%comm, iError ) call MPI_Allgather( newtree%treeID(newtree%nElems), 1, long_k_mpi, & & newtree%Part_Last, 1, long_k_mpi, & & newtree%global%comm, iError ) end if if (withprop) then newtree%global%nProperties = intree%global%nProperties allocate(newtree%global%property(newtree%global%nProperties)) allocate(newtree%property(newtree%global%nProperties)) newtree%global%property = intree%global%property do iProp=1,newtree%global%nProperties ! In the new mesh there may be a different number of elements with ! this property, recount them and update the header information ! accordingly. call tem_prop_countelems( me = newtree%global & & %Property(iProp), & & elempropertybits = newtree%Elempropertybits, & & comm = newtree%global%comm ) ! Now create the process local information on the property. call gather_Property( Property = newtree%Property(iProp), & & Header = newtree%global%Property(iProp), & & BitField = newtree%ElemPropertyBits, & & comm = newtree%global%comm ) end do else newtree%global%nProperties = 0 allocate(newtree%global%property(0)) allocate(newtree%property(0)) end if end subroutine tem_create_tree_from_sub ! ****************************************************************************** ! ! ****************************************************************************** ! !> This routine identify all the elements in inTree that has a certain property !! bit, save their positions in inTree into array: map2global, !! save the number of these elements into level wise array: countElems !! (e.g. for shape kind='property'). !! subroutine tem_shape_initPropElements( propBits, inTree, countElems, map2global ) ! --------------------------------------------------------------------------- !> shape objects on which to work integer( kind=long_k ), intent(in) :: propBits !> Global mesh from which the elements are identified and then stored to type( treelmesh_type ), intent(in) :: inTree !> How many elements there will be for each level in the track integer, intent( out ) :: countElems( globalMaxLevels ) !> growing array. Elements positions in inTree%treeID type(dyn_intArray_type), intent(inout) :: map2global ! --------------------------------------------------------------------------- integer(kind=long_k) :: elemProp, match integer :: tLevel, dPos, iElem logical :: wasAdded ! --------------------------------------------------------------------------- write(logUnit(1),*) 'Initializing shapes of elements that have a '// & & 'certain property' ! Loop over all elements in inTree do iElem = 1, inTree%nElems ! get its property bits elemProp = inTree%elemPropertyBits( iElem ) ! check whether its property match = iand( propBits, elemProp ) if( match > 0_long_k ) then ! Append to treeID list (note that already existing ones are ! omitted) call append( me = map2global, & & pos = dPos, & & val = iElem, & & wasAdded = wasAdded ) ! Count up if it was added if( wasAdded ) then tLevel = tem_levelOf( inTree%treeID(iElem) ) countElems( tLevel ) = countElems( tLevel ) + 1 end if ! wasAdded end if ! match > 0 end do ! iElem end subroutine tem_shape_initPropElements ! ****************************************************************************** ! ! ****************************************************************************** ! !> This subroutine collects the treeIDs of all elements in the level !! descriptor on the local partition (e.g. used in debug mesh, no general !! option for now). !! subroutine tem_shape_initLocal( levelDesc, treeID ) ! --------------------------------------------------------------------------- !> level descriptor including all elements (fluid, ghost, halo) type(tem_levelDesc_type ),intent(in) :: levelDesc(:) !> temporary array of treeIDs integer(kind=long_k), allocatable, intent(out) :: treeID(:) ! --------------------------------------------------------------------------- integer :: offset integer :: nTotalElems integer :: iLevel integer :: iElem ! --------------------------------------------------------------------------- write(logUnit(1),*) 'Initializing shape with all elements on this partition' ! get the total number of elements in the level descriptor nTotalElems = sum(levelDesc(:)%nElems) ! allocate the array of treeIDs allocate(treeID(nTotalElems)) offset = 0 ! copy the treeIDs to the array do iLevel = 1,size(levelDesc) do iElem = 1, size( levelDesc( iLevel )%total ) treeID( offset + iElem ) = levelDesc(iLevel)%total(iElem) end do offset = offset + size(levelDesc( iLevel )%total) end do end subroutine tem_shape_initLocal ! ****************************************************************************** ! ! ****************************************************************************** ! !> Write the complete mesh including fluid, ghosts and halo elements to disk !! !! Define in the [[tem_debug_module:tem_load_debug]] table as !!```lua !! debug = { debugMode = true, !! debugMesh = true } !!``` !! subroutine tem_write_debugMesh( globtree, levelDesc, debug, myPart ) ! --------------------------------------------------------------------------- !> mesh to locate the point in type(treelmesh_type), intent(in) :: globtree !> current level descriptor type(tem_levelDesc_type ),intent(in) :: levelDesc(:) !> debug info type(tem_debug_type ), optional, intent(inout) :: debug !> Partition to use on the calling process (= MPI Rank in comm) integer, intent(in) :: myPart ! --------------------------------------------------------------------------- type( tem_shape_type ) :: inShape(1) character(len=pathLen) :: dirname type( tem_debug_type ) :: local_debug ! local subTree that is dumped to disc type( tem_subTree_type ) :: subTree ! dummy variable which is required by tem_create_subTree_of but indeed not used type( tem_bc_prop_type ) :: bc_prop ! dummy variable which is required by tem_create_subTree_of but indeed not used type( tem_stencilHeader_type ) :: stencil ! --------------------------------------------------------------------------- if( present( debug ))then local_debug = debug else local_debug = main_debug end if if( trim( local_debug%debugMesh ) .ne. '')then write(logUnit(1),*)' Writing debug mesh to disk ' inShape(1)%shapeID = tem_local_shape write(dirname ,'(a,i6.6,a)')trim( local_debug%debugMesh ),myPart,'_' call tem_create_subTree_of( inTree = globtree, & & subTree = subTree, & & inShape = inShape, & & levelDesc = levelDesc, & & bc_prop = bc_prop, & & stencil = stencil, & & prefix = trim( dirname )) call tem_dump_subTree( subTree, globtree, root_only = .false. ) write(logUnit(1),*) 'Done writing debugMesh' endif ! output debugmesh end subroutine tem_write_debugMesh ! ****************************************************************************** ! ! ****************************************************************************** ! !> This routine identify elements that belong to certain bounaries. !! Labels of required boundaries are given by bcLabels. !! bc_prop contains boudnary_ID of all local elements. !! Firstly, bcLabels are converted into bcIDs. !! Then all elements in bc_prop are looped over to check if it matches one of !! the required bcID. If match, its position is save in map2global. !! Number of elements found on each level is saved in countElems. subroutine tem_shape_findElemByBCLabels( bcLabels, bc_prop, & & foundAny, map2global, bcIDs, stencil) ! --------------------------------------------------------------------------- !> bcLabels character(len=labelLen), intent(in) :: bcLabels(:) !> bc property type( tem_bc_prop_type ), intent(in) :: bc_prop !> if any element be identified logical, intent(out) :: foundAny !> dynamic array. Elements positions in bc_prop%property%elemID type(dyn_intArray_type), intent(inout) :: map2global !> id of boundary condition to be tracked integer, allocatable, intent(out) :: bcIDs(:) !> stencil required to get useful links type( tem_stencilHeader_type ), optional, intent(in) :: stencil ! --------------------------------------------------------------------------- integer :: dPos, iElem, iBC, nBCs, iBCtype, posInTree, QQN logical :: wasAdded, found integer, allocatable :: map(:) ! --------------------------------------------------------------------------- if( present( stencil ) ) then allocate( map( size(stencil%map) ) ) map = stencil%map qqn = stencil%QQN else allocate( map(qQQQ) ) map = qDir qqn = qQQQ end if ! write(*,*) 'inside tem_shape_findElemByBCLabels routine' foundAny = .false. ! convert bcLabels to bcIDs nBCs = size(bcLabels) allocate( bcIDs( nBCs ) ) ! write(*,*) 'size of bcLabels: ', nBCs ! loop over all required bcLabels do iBC = 1, nBCs found = .false. ! loop over all BCs in the mesh do iBCtype = 1, bc_prop%nBCtypes if ( trim(bcLabels(iBC)) == bc_prop%bc_label( iBCtype ) ) then bcIDs( iBC ) = iBCtype found = .true. ! write(*,*) 'required bc label: '//trim(bcLabels(iBC)) ! write(*,*) 'its bcID: ', bcIDs(iBC) exit end if end do if ( .not. found ) then write(logUnit(1),*) 'Required BC label: '//trim(bcLabels(iBC)) write(logUnit(1),*) 'can not be found in given mesh!' stop end if end do ! iBC ! write(*,*) 'finished find bcIDs' ! write(*,*) 'nElems of bc_prop%property:', bc_prop%property%nElems ! Loop over all element with boundary property do iElem = 1, bc_prop%property%nElems do iBC = 1, nBCs if ( any(int(bc_prop%boundary_ID( map(:qqn), iElem )) == bcIDs(iBC) ) ) then posInTree = bc_prop%property%elemID(iElem) ! Append to treeID list (note that already existing ones are ! omitted) call append( me = map2global, & & pos = dPos, & & val = posInTree, & & wasAdded = wasAdded ) ! Count up if it was added if( wasAdded ) then ! tLevel = tem_levelOf( inTree%treeID(posInTree) ) ! countElems( tLevel ) = countElems( tLevel ) + 1 foundAny = .true. ! write(*,*) 'Found a element on level: ', tLevel exit ! continue with next element end if ! wasAdded end if ! boundary_ID == bcIDs end do end do ! iElem end subroutine tem_shape_findElemByBCLabels ! ****************************************************************************** ! ! ****************************************************************************** ! subroutine tem_shape_initByLevels( inTree, minLevel, maxLevel, countElems, & & map2global ) ! --------------------------------------------------------------------------- !> Global mesh from which the elements are identified and then stored to type( treelmesh_type ), intent(in) :: inTree !> level range of target elements integer, intent(in) :: minLevel, maxLevel !> How many elements there will be for each level in the track integer, intent(out) :: countElems( globalMaxLevels ) !> growing array. Elements positions in inTree%treeID type(dyn_intArray_type), intent(inout) :: map2global ! --------------------------------------------------------------------------- integer(kind=long_k) :: myID, minID, maxID integer :: tLevel, dPos, iElem, loc_min, loc_max logical :: wasAdded ! --------------------------------------------------------------------------- loc_min = minLevel loc_max = maxLevel if ( minLevel > maxLevel ) then ! take inverse loc_min = maxLevel loc_max = minLevel end if if ( minLevel < 1 ) loc_min = 1 if ( maxLevel > globalMaxLevels ) loc_max = globalMaxLevels call tem_log(3, 'Initializing shapes by elements between level '& & //trim(tem_toStr(loc_min))//' and '//trim(tem_toStr(loc_max)) ) ! the treeID range is the first ID on min level and the last ID on max level minID = tem_firstIdAtLevel( loc_min ) maxID = tem_lastIdAtLevel( loc_max ) ! Loop over all elements in inTree do iElem = 1, inTree%nElems myID = inTree%treeID(iElem) if( (myID >= minID) .and. (myID <= maxID) ) then ! Append to treeID list (note that already existing ones are ! omitted) call append( me = map2global, & & pos = dPos, & & val = iElem, & & wasAdded = wasAdded ) ! Count up if it was added if( wasAdded ) then tLevel = tem_levelOf( inTree%treeID(iElem) ) countElems( tLevel ) = countElems( tLevel ) + 1 end if ! wasAdded end if end do ! iElem end subroutine tem_shape_initByLevels ! ****************************************************************************** ! end module tem_subTree_module ! ****************************************************************************** !