Build outMesh (a conforming-or-mortar Mesh2D_t) from a 2:1-balanced forest. baseMesh is the mesh the forest was initialised from (supplies BC metadata and the communicator; on nRanks > 1 the forest must be rank-replicated so every rank emits identical global tables). The forest must already be balanced (MaxLevelJump <= 1); EmitMesh does not mutate it.
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| type(QuadTreeMesh2D), | intent(in) | :: | forest | |||
| type(Mesh2D), | intent(in) | :: | baseMesh | |||
| type(Mesh2D), | intent(out) | :: | outMesh |
subroutine EmitMesh(forest,baseMesh,outMesh)
!! Build outMesh (a conforming-or-mortar Mesh2D_t) from a 2:1-balanced forest. baseMesh is
!! the mesh the forest was initialised from (supplies BC metadata and the communicator; on
!! nRanks > 1 the forest must be rank-replicated so every rank emits identical global
!! tables). The forest must already be balanced (MaxLevelJump <= 1); EmitMesh does not
!! mutate it.
implicit none
type(QuadTreeMesh2D),intent(in) :: forest
type(Mesh2D),intent(in) :: baseMesh
type(Mesh2D),intent(out) :: outMesh
! Local
integer :: nEl,nGeo,nBCs,li,s,node,nbr,ns,nf,e,ne,k
integer :: m,nMortar,gid,gidA,gidB,t1,t2,c1,c2,es1,es2
integer :: eFirst,eLast,nLocal
type(Lagrange) :: geomInterp
integer,allocatable :: leafIdx(:)
integer,allocatable :: si(:,:,:)
integer,allocatable :: minfo(:,:)
real(prec),allocatable :: coords(:,:,:)
if(forest%MaxLevelJump() > 1) then
print*,__FILE__,':',__LINE__, &
' : Error : EmitMesh requires a 2:1-balanced forest; call Balance2to1 first.'
stop 1
endif
nEl = forest%nLeaves
nGeo = forest%nGeo
nBCs = baseMesh%nBCs
call geomInterp%Init(nGeo,forest%quadrature,nGeo,forest%quadrature)
! node id -> emitted element id (leaf-list order); 0 for non-leaf nodes.
allocate(leafIdx(1:forest%nNodes))
leafIdx = 0
do li = 1,nEl
leafIdx(forest%leaf(li)) = li
enddo
! ---- Classify every leaf face; build sideInfo and the mortar table ----
allocate(si(1:5,1:4,1:nEl))
si = 0
allocate(minfo(1:8,1:4*nEl)) ! upper bound: at most one mortar per leaf face
nMortar = 0
gid = 0
do li = 1,nEl
node = forest%leaf(li)
e = li
do s = 1,4
if(si(2,s,e) /= 0) cycle ! already filled (conforming partner, or small side of a mortar)
call forest%FaceNeighbor(node,s,nbr,ns,nf)
if(nbr == 0) then
! Physical domain boundary.
gid = gid+1
si(2,s,e) = gid
si(5,s,e) = forest%rootBC(s,forest%rootElem(node))
elseif(forest%child(1,nbr) == 0) then
! Neighbour is a leaf.
ne = leafIdx(nbr)
if(forest%level(nbr) == forest%level(node)) then
! Conforming same-level interior side; assign a shared global id to both.
gid = gid+1
si(2,s,e) = gid
si(3,s,e) = ne
si(4,s,e) = 10*ns+nf
si(2,ns,ne) = gid
si(3,ns,ne) = e
si(4,ns,ne) = 10*s+nf
else
! Neighbour is one level coarser: this is a SMALL side; its big side fills it later.
cycle
endif
else
! Neighbour node is internal (finer) -> this leaf is the BIG side of a 2:1 mortar.
nMortar = nMortar+1
m = nMortar
! The two small elements are the finer neighbour's children on its side ns.
! Big edge coordinate [-1,0] (sub-edge 1) maps to neighbour sub-position t1, and
! [0,1] (sub-edge 2) to t2, reversed when the shared face has flip 1.
if(nf == 0) then
t1 = 1; t2 = 2
else
t1 = 2; t2 = 1
endif
c1 = forest%child(childOfSide(t1,ns),nbr)
c2 = forest%child(childOfSide(t2,ns),nbr)
es1 = leafIdx(c1)
es2 = leafIdx(c2)
gid = gid+1; gidA = gid ! sub-edge 1 (shared by the big side and small 1)
gid = gid+1; gidB = gid ! sub-edge 2
! Big side.
si(1,s,e) = m
si(2,s,e) = gidA
! Small sides (both on neighbour local side ns).
si(1,ns,es1) = m
si(2,ns,es1) = gidA
si(1,ns,es2) = m
si(2,ns,es2) = gidB
minfo(1,m) = e
minfo(2,m) = s
minfo(3,m) = es1
minfo(4,m) = 10*ns+nf
minfo(5,m) = es2
minfo(6,m) = 10*ns+nf
minfo(7,m) = gidA
minfo(8,m) = gidB
endif
enddo
enddo
! ---- Allocate and populate the output mesh (fresh contiguous decomposition) ----
! Initialize on the base mesh's communicator so MPI is reused (not re-initialized) and the
! process-wide live-decomposition count stays correct across mesh lifetimes. The
! decomposition is regenerated over the (global) leaf list, and this rank stores only its
! contiguous slice eFirst:eLast, exactly as the built-in mesh constructors do.
call outMesh%decomp%Init(comm=baseMesh%decomp%mpiComm)
call outMesh%decomp%GenerateDecomposition(nEl,64*max(gid,1))
eFirst = outMesh%decomp%offsetElem(outMesh%decomp%rankId+1)+1
eLast = outMesh%decomp%offsetElem(outMesh%decomp%rankId+2)
nLocal = eLast-eFirst+1
call outMesh%Init(nGeo,nLocal,4*nLocal,4*nLocal,nBCs)
outMesh%nGlobalElem = nEl
outMesh%nUniqueSides = gid ! GLOBAL side count on every rank (the MPI tag stride)
outMesh%quadrature = forest%quadrature
! Leaf geometry (rank-local leaves only).
allocate(coords(1:2,1:nGeo+1,1:nGeo+1))
do li = eFirst,eLast
call forest%LeafCoords(li,geomInterp,coords)
outMesh%nodeCoords(1:2,1:nGeo+1,1:nGeo+1,li-eFirst+1) = coords(1:2,1:nGeo+1,1:nGeo+1)
enddo
! Local slice of the global side table; sideInfo(3) keeps GLOBAL neighbour element ids,
! which is what SideExchange consumes (locality decided through decomp%elemToRank).
outMesh%sideInfo(1:5,1:4,1:nLocal) = si(1:5,1:4,eFirst:eLast)
outMesh%globalNodeIDs = 0 ! node ids are unused by the solver (flips are set directly)
! Boundary-condition metadata (replicated on every rank).
if(nBCs > 0) then
outMesh%BCType(1:4,1:nBCs) = baseMesh%BCType(1:4,1:nBCs)
do k = 1,nBCs
outMesh%BCNames(k) = baseMesh%BCNames(k)
enddo
endif
! Material table: each leaf inherits its root element's material (rootMaterial is global
! on the forest, so this works for any decomposition of the emitted mesh).
outMesh%nMaterials = baseMesh%nMaterials
if(allocated(outMesh%materialNames)) deallocate(outMesh%materialNames)
allocate(outMesh%materialNames(1:baseMesh%nMaterials))
outMesh%materialNames(1:baseMesh%nMaterials) = baseMesh%materialNames(1:baseMesh%nMaterials)
do li = eFirst,eLast
outMesh%elemMaterial(li-eFirst+1) = forest%rootMaterial(forest%rootElem(forest%leaf(li)))
enddo
! Mortar table.
outMesh%nMortars = nMortar
if(associated(outMesh%mortarInfo)) deallocate(outMesh%mortarInfo)
if(nMortar > 0) then
allocate(outMesh%mortarInfo(1:8,1:nMortar))
outMesh%mortarInfo(1:8,1:nMortar) = minfo(1:8,1:nMortar)
else
outMesh%mortarInfo => null()
endif
deallocate(leafIdx,si,minfo,coords)
call geomInterp%Free()
call outMesh%UpdateDevice()
endsubroutine EmitMesh