Refine all elements defined in subtree by one level in the original mesh, and create a new mesh.
Orig_mesh needs to be a properly defined treelmesh, and subtree an accompanying subtree in that mesh. orig_bcs needs to be the boundary conditions accompanying the original mesh. New_mesh will be the resulting mesh after refining the elements at elempos. New_bcs will be the resulting boundary conditions for that new mesh. Per default, the elements from the orig_mesh that are not covered by the subtree are still part of new_mesh, just not refined. To change this behavior and include only the refined elements in the new mesh, set restrict_to_sub to true.
The boundary properties are the only ones, that will be inherited to the new mesh. All other properties will get lost!!!
Type | Intent | Optional | Attributes | Name | ||
---|---|---|---|---|---|---|
type(treelmesh_type), | intent(in) | :: | orig_mesh |
The original mesh to be refined. |
||
type(tem_BC_prop_type), | intent(in) | :: | orig_bcs |
Boundary conditions for the original mesh. |
||
type(tem_subTree_type), | intent(in) | :: | subtree |
(Process local) positions of elements to refine. |
||
integer, | intent(in) | :: | ndims |
Number of dimensions for the refinement The dimensionality can restrict the elements to create. |
||
type(treelmesh_type), | intent(out) | :: | new_mesh |
Newly created refined mesh. |
||
type(tem_BC_prop_type), | intent(out) | :: | new_bcs |
Boundary conditions for the new mesh. |
||
logical, | intent(in), | optional | :: | restrict_to_sub |
Flag to indicate, wether only refined elements should be put into the new mesh. |
subroutine tem_refine_global_subtree( orig_mesh, orig_bcs, subtree, ndims, &
& new_mesh, new_bcs, restrict_to_sub )
!> The original mesh to be refined.
type(treelmesh_type), intent(in) :: orig_mesh
!> Boundary conditions for the original mesh.
type(tem_BC_prop_type), intent(in) :: orig_bcs
!> (Process local) positions of elements to refine.
type(tem_subtree_type), intent(in) :: subtree
!> Number of dimensions for the refinement
!!
!! The dimensionality can restrict the elements to create.
integer, intent(in) :: ndims
!> Newly created refined mesh.
type(treelmesh_type), intent(out) :: new_mesh
!> Boundary conditions for the new mesh.
type(tem_BC_prop_type), intent(out) :: new_bcs
!> Flag to indicate, wether only refined elements should be put into the
!! new mesh.
logical, optional, intent(in) :: restrict_to_sub
! -------------------------------------------------------------------- !
logical :: restrict
logical :: refine_all
integer :: iParent
integer :: nParents
integer :: parentpos
integer :: childpos
integer :: newpos
integer :: subpos
integer :: iError
integer :: bcprop
integer :: iBC_elem
integer :: iElem
integer :: iChild
integer :: iProp
integer :: iSide
integer(kind=long_k) :: nNewElems
integer(kind=long_k) :: nBCElems
integer(kind=long_k) :: childID_off
integer(kind=long_k) :: child_bcid(orig_bcs%nSides,2**ndims)
integer, allocatable :: bc_child(:)
type(grw_long2dArray_type) :: newbcid
logical :: has_boundary
! -------------------------------------------------------------------- !
refine_all = subtree%useGlobalMesh
if (present(restrict_to_sub)) then
restrict = restrict_to_sub
else
restrict = .false.
end if
nParents = subtree%nElems
! (local) Number of new elements:
if (restrict) then
nNewElems = nParents*(2**ndims)
else
nNewElems = nParents*(2**ndims-1) + orig_mesh%nElems
end if
new_mesh%nElems = int(nNewElems)
new_mesh%global = orig_mesh%global
! Need to create new properties for the new mesh
nullify(new_mesh%global%Property)
nullify(new_mesh%Property)
! Find the boundary property
bcprop = 0
do iProp=1,orig_mesh%global%nProperties
if (orig_mesh%global%Property(iProp)%bitpos == prp_hasBnd) then
bcprop = iProp
EXIT
end if
end do
! Only carry on boundary conditions
if (bcprop > 0) then
new_mesh%global%nProperties = 1
new_bcs%nSides = orig_bcs%nSides
new_bcs%nBCtypes = orig_bcs%nBCtypes
new_bcs%BC_label = orig_bcs%BC_label
call init( me = newbcid, &
& width = orig_bcs%nSides, &
& length = (ndims*2-1)*orig_bcs%property%nElems )
else
new_mesh%global%nProperties = 0
end if
! Adapt the level range in the new mesh.
new_mesh%global%maxlevel = new_mesh%global%maxlevel + 1
if (refine_all .or. restrict) then
new_mesh%global%minlevel = new_mesh%global%minlevel + 1
end if
! Overall number of elements in the new mesh and offsets.
call MPI_Exscan(nNewelems, new_mesh%ElemOffset, 1, long_k_mpi, &
& MPI_SUM, new_mesh%global%comm, iError )
! Initialize for rank 0
if (new_mesh%global%myPart == 0) new_mesh%ElemOffset = 0_long_k
new_mesh%global%nElems = new_mesh%ElemOffset+nNewElems
call MPI_Bcast( new_mesh%global%nElems, 1, long_k_mpi, &
& new_mesh%global%nParts-1, &
& new_mesh%global%comm, iError )
allocate(new_mesh%treeID(nNewElems))
allocate(new_mesh%ElemPropertyBits(nNewElems))
allocate(new_mesh%Part_First(orig_mesh%global%nParts))
allocate(new_mesh%Part_Last(orig_mesh%global%nParts))
newpos = 0
subpos = 1
iBC_elem = 0
! Need to iterate through all elements of the original mesh to properly
! account for the boundary elements.
origsall: do iParent=1,orig_mesh%nElems
has_boundary = btest( orig_mesh%ElemPropertyBits(iParent), &
& prp_hasbnd )
if (has_boundary) iBC_elem = iBC_elem + 1
if (refine_all) then
parentPos = iParent
else
if (nParents >= subpos) then
parentPos = subtree%map2global(subpos)
else
! Element not part of the subtree!
parentPos = -1
end if
end if
dorefine: if (parentPos == iParent) then
! Element in subtree to be refined.
childID_off = orig_mesh%treeID(iParent)*8_long_k
childpos = newpos ! store current position for boundaries
do iChild=1,2**ndims
newpos = newpos+1
new_mesh%treeID(newpos) = childID_off + iChild
new_mesh%ElemPropertyBits(newpos) = 0_long_k
end do
if (has_boundary) then
child_bcid = 0_long_k
do iSide=1,2*ndims ! Only iterate over direct neighbors
if (orig_bcs%boundary_ID(iSide, iBC_elem) > 0) then
call tem_eligibleChildren( eligible_child = bc_child, &
& direction = iSide )
do iChild=1,2**(ndims-1)
new_mesh%ElemPropertyBits(childpos+bc_child(iChild)) &
& = ibset(0_long_k, prp_hasbnd)
child_bcid(iSide, bc_child(iChild)) = orig_bcs &
& %boundary_ID( iSide, &
& iBC_elem )
end do
deallocate(bc_child)
end if
end do
do iChild=1,2**ndims
if ( btest(new_mesh%ElemPropertyBits(childpos+iChild), &
& prp_hasBnd) ) then
call append( me = newbcid, &
& val = child_bcid(:,iChild) )
end if
end do
end if
subpos = min(subpos + 1, subtree%nelems)
else dorefine
! Element NOT in subtree to be refined.
if (.not. restrict) then
! Add also non-refined elements to the new mesh.
newpos = newpos+1
new_mesh%treeID(newpos) = orig_mesh%treeID(iParent)
new_mesh%ElemPropertyBits(newpos) = 0_long_k
if (has_boundary) then
new_mesh%ElemPropertyBits(newpos) = ibset(0_long_k, prp_hasbnd)
call append( me = newbcid, &
& val = orig_bcs%boundary_ID(:,iBC_elem) )
end if
end if
end if dorefine
end do origsall
call MPI_Allgather( new_mesh%treeID(1), 1, long_k_mpi, &
& new_mesh%Part_First, 1, long_k_mpi, &
& new_mesh%global%comm, iError )
call MPI_Allgather( new_mesh%treeID(new_mesh%nElems), 1, long_k_mpi, &
& new_mesh%Part_Last, 1, long_k_mpi, &
& new_mesh%global%comm, iError )
! Set boundary property, if needed.
if (bcprop > 0) then
allocate(new_mesh%Property(1))
allocate(new_mesh%global%Property(1))
new_mesh%global%property(1)%label &
& = orig_mesh%global%Property(bcprop)%label
new_mesh%global%property(1)%bitpos &
& = orig_mesh%global%Property(bcprop)%bitpos
new_mesh%Property(1)%nElems = newbcid%nVals
nBCElems = int(newbcid%nVals, kind=long_k)
allocate(new_mesh%Property(1)%ElemID(newbcid%nVals))
iBC_elem = 0
do iElem=1,int(nNewElems)
if ( btest(new_mesh%ElemPropertyBits(iElem), prp_hasBnd) ) then
iBC_elem = iBC_elem + 1
new_mesh%Property(1)%ElemID(iBC_elem) = iElem
end if
end do
allocate(new_bcs%boundary_ID(new_bcs%nSides, newbcid%nVals))
new_bcs%property => new_mesh%Property(1)
new_bcs%header => new_mesh%global%Property(1)
new_bcs%boundary_ID = newbcid%val(:,:newbcid%nVals)
call destroy(newbcid)
call MPI_Exscan(nBCElems, new_bcs%property%Offset, 1, long_k_mpi, &
& MPI_SUM, new_mesh%global%comm, iError )
nBCElems = new_bcs%property%Offset+nBCElems
call MPI_Bcast(nBCElems, 1, long_k_mpi, new_mesh%global%nParts-1, &
& new_mesh%global%comm, iError )
new_bcs%header%nElems = nBCElems
else
if (associated(new_mesh%property)) &
& deallocate(new_mesh%property)
if (associated(new_mesh%global%property)) &
& deallocate(new_mesh%global%property)
allocate(new_mesh%Property(0))
allocate(new_mesh%global%Property(0))
end if
end subroutine tem_refine_global_subtree