! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ! ! ! 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_QuadTreeMesh_2D !! Forest-of-quadtrees data structure for adaptive h-refinement of 2-D quadrilateral meshes !! (AMR Stage 2b). Each base-mesh element is the root of a quadtree; refinement replaces a leaf !! with four children (one per reference sub-quadrant, SELF child ordering 1=SW,2=SE,3=NE,4=NW), !! and coarsening merges a family of four leaf siblings back into their parent. The forest tracks !! refinement levels and the active leaf set; the physical geometry of any leaf is produced by !! repeated exact isoparametric subdivision of its root (SELF_RefinementPrimitives_2D). !! !! This module owns the adaptive *mesh mutation* only. Turning an adaptively refined (and hence !! generally nonconforming) forest into a solver-ready Mesh2D_t - face-neighbour queries, 2:1 !! balancing, hanging-node/mortar generation - is AMR Stage 4, and dynamic MPI re-partitioning is !! Stage 5. The intended driver loop each adaptation step is: !! !! indicator%Estimate(solution,ivar) ! AMR Stage 1 : per-leaf refine/keep/coarsen flags !! forest%AdaptFromFlags(indicator%flag) ! AMR Stage 2b : mutate the forest (this module) !! ... Stage 4 : balance + emit Mesh2D_t + mortars ; Stage 3 : transfer the solution ... !! !! Coarsening does not currently reclaim the storage of removed child nodes (the node arrays grow !! monotonically); the active leaf set is always recovered by traversal from the roots, so !! orphaned nodes are simply never revisited. Storage compaction is a future optimisation. use SELF_Constants use SELF_Lagrange use SELF_Mesh_2D use SELF_RefinementPrimitives_2D implicit none ! Per-leaf adaptation flags (chosen to match SELF_RefinementIndicator_2D's SELF_AMR_* codes). integer,parameter :: QUADTREE_COARSEN = -1 integer,parameter :: QUADTREE_KEEP = 0 integer,parameter :: QUADTREE_REFINE = 1 type :: QuadTreeMesh2D ! ---- Base-mesh (root) reference data ---- integer :: nGeo = 0 integer :: quadrature = 0 integer :: nRoots = 0 real(prec),allocatable :: rootCoords(:,:,:,:) ! (2,nGeo+1,nGeo+1,nRoots) ! Base-mesh face connectivity at the roots (conforming base assumed). rootNbr(s,r) is the ! base neighbour element of root r across local side s (0 = physical boundary), which is also ! that neighbour's root node id; rootNbrSide / rootFlip decode base sideInfo(4,s,r). integer,allocatable :: rootNbr(:,:) ! (4,nRoots) integer,allocatable :: rootNbrSide(:,:) ! (4,nRoots) integer,allocatable :: rootFlip(:,:) ! (4,nRoots) integer,allocatable :: rootBC(:,:) ! (4,nRoots) base boundary-condition id per side integer,allocatable :: rootMaterial(:) ! (nRoots) base material id per root element ! ---- Forest node storage (roots occupy node ids 1:nRoots) ---- integer :: nNodes = 0 integer :: capacity = 0 integer,allocatable :: level(:) ! refinement level (0 at roots) integer,allocatable :: parent(:) ! parent node id (0 at roots) integer,allocatable :: quadrant(:) ! child index within parent (1..4; 0 at roots) integer,allocatable :: rootElem(:) ! base element this node descends from integer,allocatable :: child(:,:) ! (4,capacity) child node ids (0 => leaf) ! ---- Active leaf set (rebuilt after every mutation) ---- integer :: nLeaves = 0 integer,allocatable :: leaf(:) ! (nLeaves) leaf node ids, root-major DFS order contains procedure,public :: Init => Init_QuadTreeMesh2D procedure,public :: InitGlobal => InitGlobal_QuadTreeMesh2D procedure,public :: Free => Free_QuadTreeMesh2D procedure,public :: RefineNode => RefineNode_QuadTreeMesh2D procedure,public :: AdaptFromFlags => AdaptFromFlags_QuadTreeMesh2D procedure,public :: RebuildLeaves => RebuildLeaves_QuadTreeMesh2D procedure,public :: LeafCoords => LeafCoords_QuadTreeMesh2D procedure,public :: MaxLevel => MaxLevel_QuadTreeMesh2D procedure,public :: FaceNeighbor => FaceNeighbor_QuadTreeMesh2D procedure,public :: Balance2to1 => Balance2to1_QuadTreeMesh2D procedure,public :: MaxLevelJump => MaxLevelJump_QuadTreeMesh2D procedure,private :: EnsureCapacity => EnsureCapacity_QuadTreeMesh2D endtype QuadTreeMesh2D contains subroutine Init_QuadTreeMesh2D(this,mesh) !! Initialize the forest with one root per base-mesh element (all leaves at level 0). The !! base geometry (node coordinates) is copied so leaf geometry can be regenerated after any !! amount of refinement without holding a reference to the mesh. Requires a single-rank !! mesh (a decomposed mesh only stores its local elements); a rank-replicated forest over a !! decomposed base is built by gathering the global tables and calling InitGlobal. implicit none class(QuadTreeMesh2D),intent(out) :: this type(Mesh2D),intent(in) :: mesh ! Local integer :: r,s integer,allocatable :: nbr(:,:),nbrSide(:,:),flip(:,:),bc(:,:),mat(:) if(mesh%decomp%nRanks > 1) then print*,__FILE__,':',__LINE__, & ' : Error : QuadTreeMesh2D%Init requires a single-rank mesh; gather the global base'// & ' tables and call InitGlobal for a decomposed base mesh.' stop 1 endif allocate(nbr(1:4,1:mesh%nElem),nbrSide(1:4,1:mesh%nElem)) allocate(flip(1:4,1:mesh%nElem),bc(1:4,1:mesh%nElem)) allocate(mat(1:mesh%nElem)) do r = 1,mesh%nElem do s = 1,4 nbr(s,r) = mesh%sideInfo(3,s,r) nbrSide(s,r) = mesh%sideInfo(4,s,r)/10 flip(s,r) = mod(mesh%sideInfo(4,s,r),10) bc(s,r) = mesh%sideInfo(5,s,r) enddo mat(r) = mesh%elemMaterial(r) enddo call this%InitGlobal(mesh%nElem,mesh%nGeo,mesh%quadrature, & mesh%nodeCoords,nbr,nbrSide,flip,bc,mat) deallocate(nbr,nbrSide,flip,bc,mat) endsubroutine Init_QuadTreeMesh2D subroutine InitGlobal_QuadTreeMesh2D(this,nRoots,nGeo,quadrature,rootCoords, & rootNbr,rootNbrSide,rootFlip,rootBC,rootMaterial) !! Initialize the forest directly from GLOBAL base-mesh tables (one root per global base !! element, all leaves at level 0). This is the initialization path for a rank-replicated !! forest over a decomposed base mesh (AMR Stage 5): every rank passes the same gathered !! tables and holds an identical forest. rootNbr carries global element ids (0 = physical !! boundary), rootNbrSide/rootFlip decode the base sideInfo(4) pairing, rootBC the base !! boundary-condition id per side, and rootMaterial the base material id per element. implicit none class(QuadTreeMesh2D),intent(out) :: this integer,intent(in) :: nRoots integer,intent(in) :: nGeo integer,intent(in) :: quadrature real(prec),intent(in) :: rootCoords(1:2,1:nGeo+1,1:nGeo+1,1:nRoots) integer,intent(in) :: rootNbr(1:4,1:nRoots) integer,intent(in) :: rootNbrSide(1:4,1:nRoots) integer,intent(in) :: rootFlip(1:4,1:nRoots) integer,intent(in) :: rootBC(1:4,1:nRoots) integer,intent(in) :: rootMaterial(1:nRoots) ! Local integer :: r this%nGeo = nGeo this%quadrature = quadrature this%nRoots = nRoots allocate(this%rootCoords(1:2,1:nGeo+1,1:nGeo+1,1:nRoots)) this%rootCoords(1:2,1:nGeo+1,1:nGeo+1,1:nRoots) = & rootCoords(1:2,1:nGeo+1,1:nGeo+1,1:nRoots) allocate(this%rootNbr(1:4,1:nRoots)) allocate(this%rootNbrSide(1:4,1:nRoots)) allocate(this%rootFlip(1:4,1:nRoots)) allocate(this%rootBC(1:4,1:nRoots)) allocate(this%rootMaterial(1:nRoots)) this%rootNbr(1:4,1:nRoots) = rootNbr(1:4,1:nRoots) this%rootNbrSide(1:4,1:nRoots) = rootNbrSide(1:4,1:nRoots) this%rootFlip(1:4,1:nRoots) = rootFlip(1:4,1:nRoots) this%rootBC(1:4,1:nRoots) = rootBC(1:4,1:nRoots) this%rootMaterial(1:nRoots) = rootMaterial(1:nRoots) ! Roots are the first nRoots nodes. this%capacity = max(4*this%nRoots,16) this%nNodes = this%nRoots allocate(this%level(1:this%capacity)) allocate(this%parent(1:this%capacity)) allocate(this%quadrant(1:this%capacity)) allocate(this%rootElem(1:this%capacity)) allocate(this%child(1:4,1:this%capacity)) this%level = 0 this%parent = 0 this%quadrant = 0 this%rootElem = 0 this%child = 0 do r = 1,this%nRoots this%rootElem(r) = r enddo call this%RebuildLeaves() endsubroutine InitGlobal_QuadTreeMesh2D subroutine Free_QuadTreeMesh2D(this) implicit none class(QuadTreeMesh2D),intent(inout) :: this if(allocated(this%rootCoords)) deallocate(this%rootCoords) if(allocated(this%rootNbr)) deallocate(this%rootNbr) if(allocated(this%rootNbrSide)) deallocate(this%rootNbrSide) if(allocated(this%rootFlip)) deallocate(this%rootFlip) if(allocated(this%rootBC)) deallocate(this%rootBC) if(allocated(this%rootMaterial)) deallocate(this%rootMaterial) if(allocated(this%level)) deallocate(this%level) if(allocated(this%parent)) deallocate(this%parent) if(allocated(this%quadrant)) deallocate(this%quadrant) if(allocated(this%rootElem)) deallocate(this%rootElem) if(allocated(this%child)) deallocate(this%child) if(allocated(this%leaf)) deallocate(this%leaf) this%nGeo = 0 this%nRoots = 0 this%nNodes = 0 this%capacity = 0 this%nLeaves = 0 endsubroutine Free_QuadTreeMesh2D subroutine EnsureCapacity_QuadTreeMesh2D(this,need) !! Grow the node arrays (amortized doubling) so at least `need` nodes fit. implicit none class(QuadTreeMesh2D),intent(inout) :: this integer,intent(in) :: need ! Local integer :: newCap integer,allocatable :: itmp(:),ctmp(:,:) if(need <= this%capacity) return newCap = this%capacity do while(newCap < need) newCap = 2*newCap enddo allocate(itmp(1:newCap)) itmp = 0; itmp(1:this%nNodes) = this%level(1:this%nNodes); call move_alloc(itmp,this%level) allocate(itmp(1:newCap)) itmp = 0; itmp(1:this%nNodes) = this%parent(1:this%nNodes); call move_alloc(itmp,this%parent) allocate(itmp(1:newCap)) itmp = 0; itmp(1:this%nNodes) = this%quadrant(1:this%nNodes); call move_alloc(itmp,this%quadrant) allocate(itmp(1:newCap)) itmp = 0; itmp(1:this%nNodes) = this%rootElem(1:this%nNodes); call move_alloc(itmp,this%rootElem) allocate(ctmp(1:4,1:newCap)) ctmp = 0; ctmp(1:4,1:this%nNodes) = this%child(1:4,1:this%nNodes); call move_alloc(ctmp,this%child) this%capacity = newCap endsubroutine EnsureCapacity_QuadTreeMesh2D subroutine RefineNode_QuadTreeMesh2D(this,node) !! Subdivide one leaf node into four children. Does nothing (with a warning) if the node is !! already refined. The caller is responsible for rebuilding the leaf list afterwards (or !! calling AdaptFromFlags, which does so). implicit none class(QuadTreeMesh2D),intent(inout) :: this integer,intent(in) :: node ! Local integer :: c,newid if(this%child(1,node) /= 0) then print*,__FILE__,':',__LINE__,' : Warning : RefineNode called on an already-refined node.' return endif call this%EnsureCapacity(this%nNodes+4) do c = 1,4 newid = this%nNodes+c this%level(newid) = this%level(node)+1 this%parent(newid) = node this%quadrant(newid) = c this%rootElem(newid) = this%rootElem(node) this%child(1:4,newid) = 0 enddo this%child(1,node) = this%nNodes+1 this%child(2,node) = this%nNodes+2 this%child(3,node) = this%nNodes+3 this%child(4,node) = this%nNodes+4 this%nNodes = this%nNodes+4 endsubroutine RefineNode_QuadTreeMesh2D subroutine AdaptFromFlags_QuadTreeMesh2D(this,flag) !! Mutate the forest from a per-leaf flag array (indexed over the current leaves in this%leaf !! order, e.g. the flag array produced by SELF_RefinementIndicator_2D): !! !! flag(i) = QUADTREE_REFINE (+1) -> subdivide leaf i !! flag(i) = QUADTREE_COARSEN (-1) -> merge leaf i's family if ALL four siblings are leaves !! and ALL are flagged COARSEN (standard coarsening rule) !! flag(i) = QUADTREE_KEEP (0) -> unchanged !! !! Refinement and coarsening are both resolved against the pre-adaptation leaf snapshot, so !! the two operations never interfere. The leaf list is rebuilt on return. implicit none class(QuadTreeMesh2D),intent(inout) :: this integer,intent(in) :: flag(1:this%nLeaves) ! Local integer :: i,c,p,node,nSnap integer,allocatable :: flagOfNode(:) integer,allocatable :: refineList(:) logical :: family nSnap = this%nLeaves ! Map the flag onto node ids so coarsening can test whole families by node. allocate(flagOfNode(1:this%nNodes)) flagOfNode = QUADTREE_KEEP allocate(refineList(1:nSnap)) refineList = 0 do i = 1,nSnap flagOfNode(this%leaf(i)) = flag(i) if(flag(i) == QUADTREE_REFINE) refineList(i) = this%leaf(i) enddo ! ---- Coarsening : merge families all of whose four leaf children are flagged COARSEN ---- do i = 1,nSnap if(flag(i) /= QUADTREE_COARSEN) cycle node = this%leaf(i) p = this%parent(node) if(p == 0) cycle ! a root cannot be coarsened ! Only act once per family (when processing its first child), and only if every child is a ! leaf flagged COARSEN. if(this%child(1,p) /= node) cycle family = .true. do c = 1,4 if(this%child(1,this%child(c,p)) /= 0) family = .false. ! child not a leaf if(flagOfNode(this%child(c,p)) /= QUADTREE_COARSEN) family = .false. enddo if(family) this%child(1:4,p) = 0 ! detach children -> p becomes a leaf again enddo ! ---- Refinement : subdivide flagged leaves that are still leaves ---- do i = 1,nSnap if(refineList(i) == 0) cycle node = refineList(i) if(this%child(1,node) == 0) call this%RefineNode(node) enddo deallocate(flagOfNode,refineList) call this%RebuildLeaves() endsubroutine AdaptFromFlags_QuadTreeMesh2D subroutine RebuildLeaves_QuadTreeMesh2D(this) !! Recompute the active leaf set by depth-first traversal from the roots. Traversal (rather !! than a scan of all nodes) is what makes orphaned nodes left behind by coarsening invisible. implicit none class(QuadTreeMesh2D),intent(inout) :: this ! Local integer :: r,n ! First pass: count leaves. n = 0 do r = 1,this%nRoots call count_leaves(r,n) enddo this%nLeaves = n if(allocated(this%leaf)) deallocate(this%leaf) allocate(this%leaf(1:max(n,1))) this%leaf = 0 ! Second pass: collect leaves in root-major DFS order. n = 0 do r = 1,this%nRoots call collect_leaves(r,n) enddo contains recursive subroutine count_leaves(node,cnt) integer,intent(in) :: node integer,intent(inout) :: cnt integer :: k if(this%child(1,node) == 0) then cnt = cnt+1 else do k = 1,4 call count_leaves(this%child(k,node),cnt) enddo endif endsubroutine count_leaves recursive subroutine collect_leaves(node,idx) integer,intent(in) :: node integer,intent(inout) :: idx integer :: k if(this%child(1,node) == 0) then idx = idx+1 this%leaf(idx) = node else do k = 1,4 call collect_leaves(this%child(k,node),idx) enddo endif endsubroutine collect_leaves endsubroutine RebuildLeaves_QuadTreeMesh2D function MaxLevel_QuadTreeMesh2D(this) result(mx) !! Highest refinement level among the active leaves. implicit none class(QuadTreeMesh2D),intent(in) :: this integer :: mx integer :: i mx = 0 do i = 1,this%nLeaves mx = max(mx,this%level(this%leaf(i))) enddo endfunction MaxLevel_QuadTreeMesh2D subroutine LeafCoords_QuadTreeMesh2D(this,leafIndex,geomInterp,coords) !! Physical geometry-node coordinates of leaf `leafIndex`, produced by repeated exact !! isoparametric subdivision of its root element along the quadtree path. geomInterp must be a !! degree-nGeo Lagrange interpolant on the mesh's geometry (quadrature) nodes - the same !! interpolant SELF_MeshRefinement_2D builds. Level 0 leaves return the root geometry directly. implicit none class(QuadTreeMesh2D),intent(in) :: this integer,intent(in) :: leafIndex type(Lagrange),intent(in) :: geomInterp real(prec),intent(out) :: coords(1:2,1:this%nGeo+1,1:this%nGeo+1) ! Local integer :: node,lvl,step integer,allocatable :: path(:) real(prec),allocatable :: cur(:,:,:),kids(:,:,:,:) node = this%leaf(leafIndex) lvl = this%level(node) allocate(cur(1:2,1:this%nGeo+1,1:this%nGeo+1)) if(lvl == 0) then cur(1:2,:,:) = this%rootCoords(1:2,:,:,this%rootElem(node)) coords = cur deallocate(cur) return endif ! Path of quadrant indices from root (step 1) down to the leaf (step lvl). allocate(path(1:lvl)) do step = lvl,1,-1 path(step) = this%quadrant(node) node = this%parent(node) enddo ! `node` is now the root. cur(1:2,:,:) = this%rootCoords(1:2,:,:,this%rootElem(this%leaf(leafIndex))) allocate(kids(1:2,1:this%nGeo+1,1:this%nGeo+1,1:4)) do step = 1,lvl call SubdivideNodeCoords(geomInterp,this%nGeo,cur,kids) cur(1:2,:,:) = kids(1:2,:,:,path(step)) enddo coords = cur deallocate(cur,kids,path) endsubroutine LeafCoords_QuadTreeMesh2D recursive subroutine FaceNeighbor_QuadTreeMesh2D(this,node,s,nbr,ns,nf) !! Find the equal-or-larger face neighbour of `node` across its local side s using the !! classic quadtree ascend/descend search. Returns: !! nbr - neighbour node id (0 = physical domain boundary). It is either a LEAF at any level !! <= level(node), or an INTERNAL node at exactly level(node) (meaning the shared face !! is subdivided on the neighbour side, i.e. finer neighbours exist). !! ns - the neighbour's local side that faces `node`. !! nf - the flip between the two shared edges (0 same direction, 1 reversed), inherited !! from the base-mesh face where the search crosses a root boundary. !! With this, a 2:1 hanging face is exactly "nbr is a leaf with level(nbr) = level(node)-1", !! and finer neighbours are exactly "nbr is internal". implicit none class(QuadTreeMesh2D),intent(in) :: this integer,intent(in) :: node,s integer,intent(out) :: nbr,ns,nf ! Local integer :: p,c,pnbr,ps,pf,t,tq if(this%level(node) == 0) then ! Cross a base-mesh face: neighbour root, its side, and the base flip. nbr = this%rootNbr(s,this%rootElem(node)) ns = this%rootNbrSide(s,this%rootElem(node)) nf = this%rootFlip(s,this%rootElem(node)) return endif p = this%parent(node) c = this%quadrant(node) if(qt_internal(s,c)) then ! Neighbour is the sibling on the other side of an interior face of the parent. nbr = this%child(qt_reflect(s,c),p) ns = qt_opposite(s) nf = 0 return endif ! Otherwise ascend: find the parent's neighbour across the same side, then descend. call FaceNeighbor_QuadTreeMesh2D(this,p,s,pnbr,ps,pf) if(pnbr == 0) then nbr = 0; ns = 0; nf = 0 return endif if(this%child(1,pnbr) == 0) then ! Parent's neighbour is a leaf (equal or larger than the parent) -> our larger neighbour. nbr = pnbr; ns = ps; nf = pf return endif ! Parent's neighbour is internal (at level(node)-1): descend one level to the child that ! borders `node`, matching sub-positions across the face through the flip. t = qt_subpos(c,s) if(pf == 0) then tq = t else tq = 3-t endif nbr = this%child(childOfSide(tq,ps),pnbr) ns = ps nf = pf endsubroutine FaceNeighbor_QuadTreeMesh2D subroutine Balance2to1_QuadTreeMesh2D(this) !! Enforce the 2:1 balance condition: no leaf face may separate elements differing by more !! than one refinement level. Iterates to a fixed point - in each sweep, any leaf whose !! equal-or-larger neighbour is a leaf two or more levels coarser triggers refinement of that !! coarser neighbour; refinement can ripple, so sweeps repeat until nothing changes. The leaf !! set is rebuilt on return. implicit none class(QuadTreeMesh2D),intent(inout) :: this ! Local integer :: li,s,node,nbr,ns,nf,nSnap integer,allocatable :: snap(:) logical :: changed do changed = .false. nSnap = this%nLeaves allocate(snap(1:nSnap)) snap(1:nSnap) = this%leaf(1:nSnap) do li = 1,nSnap node = snap(li) if(this%child(1,node) /= 0) cycle ! refined earlier in this sweep do s = 1,4 call this%FaceNeighbor(node,s,nbr,ns,nf) if(nbr /= 0) then if(this%child(1,nbr) == 0 .and. this%level(nbr) <= this%level(node)-2) then call this%RefineNode(nbr) changed = .true. endif endif enddo enddo deallocate(snap) call this%RebuildLeaves() if(.not. changed) exit enddo endsubroutine Balance2to1_QuadTreeMesh2D function MaxLevelJump_QuadTreeMesh2D(this) result(mx) !! Largest refinement-level difference across any leaf face (0 on a conforming or uniformly !! refined forest, 1 on a 2:1-balanced adaptive forest). Because FaceNeighbor returns the !! equal-or-larger neighbour, every level difference is observed from the finer leaf as a !! coarser leaf neighbour; internal (finer) neighbours contribute nothing from this side. implicit none class(QuadTreeMesh2D),intent(in) :: this integer :: mx ! Local integer :: li,s,node,nbr,ns,nf mx = 0 do li = 1,this%nLeaves node = this%leaf(li) do s = 1,4 call this%FaceNeighbor(node,s,nbr,ns,nf) if(nbr /= 0) then if(this%child(1,nbr) == 0) mx = max(mx,this%level(node)-this%level(nbr)) endif enddo enddo endfunction MaxLevelJump_QuadTreeMesh2D ! -------------------------------------------------------------------------------------------- ! ! Quadtree face-adjacency helpers (SELF child ordering 1=SW,2=SE,3=NE,4=NW; sides 1=S,2=E,3=N, ! 4=W). See SELF_RefinementPrimitives_2D for the quadrant/side conventions and childOfSide. ! -------------------------------------------------------------------------------------------- ! pure function qt_opposite(s) result(o) !! The local side directly across an element from side s. implicit none integer,intent(in) :: s integer :: o integer,parameter :: opp(1:4) = [3,4,1,2] o = opp(s) endfunction qt_opposite pure function qt_internal(s,c) result(isInternal) !! .true. if child c's side s is interior to its parent (its neighbour across s is a sibling). implicit none integer,intent(in) :: s,c logical :: isInternal select case(s) case(1); isInternal = (c == 3 .or. c == 4) ! South interior for the top children case(2); isInternal = (c == 1 .or. c == 4) ! East interior for the left children case(3); isInternal = (c == 1 .or. c == 2) ! North interior for the bottom children case(4); isInternal = (c == 2 .or. c == 3) ! West interior for the right children case default; isInternal = .false. endselect endfunction qt_internal pure function qt_reflect(s,c) result(rc) !! Sibling child index obtained by reflecting c across side s (swap the y-half for the !! horizontal faces S/N, the x-half for the vertical faces E/W). implicit none integer,intent(in) :: s,c integer :: rc integer,parameter :: vref(1:4) = [4,3,2,1] ! swap ay (s = 1,3) integer,parameter :: href(1:4) = [2,1,4,3] ! swap ax (s = 2,4) if(s == 1 .or. s == 3) then rc = vref(c) else rc = href(c) endif endfunction qt_reflect pure function qt_subpos(c,s) result(t) !! Sub-position (1 or 2, in the positive direction of side s) of child c along side s; !! 0 if c does not touch side s. implicit none integer,intent(in) :: c,s integer :: t if(childOfSide(1,s) == c) then t = 1 elseif(childOfSide(2,s) == c) then t = 2 else t = 0 endif endfunction qt_subpos endmodule SELF_QuadTreeMesh_2D