Create the smallest 3D 2:1 nonconforming (mortar) mesh: one "big" element of size 2dx x 2dx x 2*dx whose east face is shared with the faces of four "small" dx x dx x dx elements stacked 2x2 in (y,z). The mesh is conforming everywhere except at the single mortar interface.
z = 2*dx ________________ ______ ______
| | e4 | e5 | (view of the x = 2*dx plane;
| |______|______| e2..e5 sit in the big face's
| e1 | e2 | e3 | quadrants 1..4)
|________________|______|______|
x = 0 2*dx 3*dx
Input - this : Fresh/empty Mesh3D_t object - dx : Edge length of the small elements; the big element has edge length 2*dx - bcids(1:6) : Boundary condition flags for the bottom, south, east, north, west, and top sides of the domain (the side ordering of StructuredMesh) - flips(1:4) (optional) : requested face flip (0..7) for the small element in each big-face quadrant. Each small element is rigidly rotated so that its face on the mortar interface meets the big face with the requested flip; flips 0,2,5,7 place the small element's west face on the interface, flips 1,3,4,6 its east face (the flip fixes the face chirality). Defaults to 0 (no rotation).
Output - this : Mesh3D_t object with five elements and one mortar interface
All connectivity (conforming interior faces, boundary faces, and the mortar table) is derived by sampling the bilinear face maps and matching physical positions, so the emitted sideInfo/mortarInfo are correct by construction for every rotation. Element geometry is trilinear (nGeo=1). With domain decomposition on two or more ranks the mortar interface straddles the rank boundary.
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| class(Mesh3D_t), | intent(out) | :: | this | |||
| real(kind=prec), | intent(in) | :: | dx | |||
| integer, | intent(in) | :: | bcids(1:6) | |||
| integer, | intent(in), | optional | :: | comm | ||
| integer, | intent(in), | optional | :: | flips(1:4) |
subroutine SimpleMortarMesh_Mesh3D_t(this,dx,bcids,comm,flips)
!!
!! Create the smallest 3D 2:1 nonconforming (mortar) mesh: one "big" element of size
!! 2*dx x 2*dx x 2*dx whose east face is shared with the faces of four "small"
!! dx x dx x dx elements stacked 2x2 in (y,z). The mesh is conforming everywhere
!! except at the single mortar interface.
!!
!! z = 2*dx ________________ ______ ______
!! | | e4 | e5 | (view of the x = 2*dx plane;
!! | |______|______| e2..e5 sit in the big face's
!! | e1 | e2 | e3 | quadrants 1..4)
!! |________________|______|______|
!! x = 0 2*dx 3*dx
!!
!! Input
!! - this : Fresh/empty Mesh3D_t object
!! - dx : Edge length of the small elements; the big element has edge length 2*dx
!! - bcids(1:6) : Boundary condition flags for the bottom, south, east, north,
!! west, and top sides of the domain (the side ordering of
!! StructuredMesh)
!! - flips(1:4) (optional) : requested face flip (0..7) for the small element in
!! each big-face quadrant. Each small element is rigidly rotated so
!! that its face on the mortar interface meets the big face with
!! the requested flip; flips 0,2,5,7 place the small element's west
!! face on the interface, flips 1,3,4,6 its east face (the flip
!! fixes the face chirality). Defaults to 0 (no rotation).
!!
!! Output
!! - this : Mesh3D_t object with five elements and one mortar interface
!!
!! All connectivity (conforming interior faces, boundary faces, and the mortar
!! table) is derived by sampling the bilinear face maps and matching physical
!! positions, so the emitted sideInfo/mortarInfo are correct by construction for
!! every rotation. Element geometry is trilinear (nGeo=1). With domain
!! decomposition on two or more ranks the mortar interface straddles the rank
!! boundary.
!!
implicit none
class(Mesh3D_t),intent(out) :: this
real(prec),intent(in) :: dx
integer,intent(in) :: bcids(1:6)
integer,intent(in),optional :: comm
integer,intent(in),optional :: flips(1:4)
! Local
integer,parameter :: nGlobalElem = 5
! Signed local axis directions (1=x,2=y,3=z) of the rotated small elements for
! each requested flip; column f holds (d1,d2,d3) such that the face meeting the
! big east face does so with flip f. The donor face is West (5) when d1 = +x and
! East (3) when d1 = -x.
integer,parameter :: flipAxes(1:3,0:7) = reshape([ &
1,2,3, & ! flip 0 : West donor
-1,-2,3, & ! flip 1 : East donor
1,-2,-3, & ! flip 2 : West donor
-1,2,-3, & ! flip 3 : East donor
-1,3,2, & ! flip 4 : East donor
1,-3,2, & ! flip 5 : West donor
-1,-3,-2, & ! flip 6 : East donor
1,3,-2],[3,8]) ! flip 7 : West donor
real(prec) :: nodeCoords(1:3,1:2,1:2,1:2,1:nGlobalElem)
integer :: sideInfo(1:5,1:6,1:nGlobalElem)
integer :: minfo(1:14)
integer :: localFlips(1:4)
integer,parameter :: invFlip(0:7) = [0,1,2,3,4,7,6,5]
real(prec) :: center(1:3),axis(1:3,1:3)
real(prec) :: tol
integer :: q,f,i,j,k,d
integer :: eA,eB,sA,sB,fAB
integer :: e1,e2
integer :: gid,nUniqueSides
integer :: nLocalElems
integer :: nGeo,nBCs
logical :: matched
call this%decomp%init(comm)
nGeo = 1 ! Trilinear element geometry
nBCs = 6
tol = 1.0e-3_prec*dx
localFlips = 0
if(present(flips)) then
localFlips = flips
do q = 1,4
if(localFlips(q) < 0 .or. localFlips(q) > 7) then
print*,"SimpleMortarMesh_Mesh3D_t : requested flip out of range 0..7"
stop 1
endif
enddo
endif
! Element 1 (big) : [0,2dx]^3, axis aligned
do k = 1,2
do j = 1,2
do i = 1,2
nodeCoords(1:3,i,j,k,1) = [real(i-1,prec),real(j-1,prec),real(k-1,prec)]* &
2.0_prec*dx
enddo
enddo
enddo
! Elements 2..5 (small) : [2dx,3dx] x quadrant q of the big east face, each
! rigidly rotated per the requested flip
do q = 1,4
f = localFlips(q)
center(1) = 2.5_prec*dx
center(2) = (real(mortarQuadKx(q),prec)-0.5_prec)*dx
center(3) = (real(mortarQuadKy(q),prec)-0.5_prec)*dx
do d = 1,3
axis(1:3,d) = 0.0_prec
axis(abs(flipAxes(d,f)),d) = real(sign(1,flipAxes(d,f)),prec)
enddo
do k = 1,2
do j = 1,2
do i = 1,2
nodeCoords(1:3,i,j,k,1+q) = center+ &
0.5_prec*dx*(real(2*i-3,prec)*axis(1:3,1)+ &
real(2*j-3,prec)*axis(1:3,2)+ &
real(2*k-3,prec)*axis(1:3,3))
enddo
enddo
enddo
enddo
sideInfo = 0
! Conforming interior faces: match every pair of element faces by sampling the
! bilinear face maps under each of the eight flips.
do eA = 1,nGlobalElem
do eB = eA+1,nGlobalElem
do sA = 1,6
if(sideInfo(3,sA,eA) /= 0) cycle
do sB = 1,6
call MatchFaces3D(nodeCoords(:,:,:,:,eA),sA, &
nodeCoords(:,:,:,:,eB),sB,tol,matched,fAB)
if(matched) then
sideInfo(3,sA,eA) = eB
sideInfo(4,sA,eA) = 10*sB+fAB
sideInfo(3,sB,eB) = eA
sideInfo(4,sB,eB) = 10*sA+invFlip(fAB)
exit
endif
enddo
enddo
enddo
enddo
! Mortar interface: big element east face against each small element. The donor
! face and flip are recovered by matching each quadrant of the big face.
minfo = 0
minfo(1) = 1
minfo(2) = selfSide3D_East
do q = 1,4
matched = .false.
do sB = 1,6
do f = 0,7
if(QuadrantMatches3D(nodeCoords(:,:,:,:,1),selfSide3D_East,q, &
nodeCoords(:,:,:,:,1+q),sB,f,tol)) then
minfo(2*q+1) = 1+q
minfo(2*q+2) = 10*sB+f
matched = .true.
exit
endif
enddo
if(matched) exit
enddo
if(.not. matched) then
print*,"SimpleMortarMesh_Mesh3D_t : failed to match mortar sub-face",q
stop 1
endif
enddo
! Physical boundary conditions: any face that is neither conforming nor on the
! mortar interface lies on one of the six domain boundary planes; identify the
! plane from the face centroid.
do eA = 1,nGlobalElem
do sA = 1,6
if(sideInfo(3,sA,eA) /= 0) cycle
if(eA == 1 .and. sA == selfSide3D_East) cycle
matched = .false.
do q = 1,4
if(eA == minfo(2*q+1) .and. sA == minfo(2*q+2)/10) matched = .true.
enddo
if(matched) cycle ! small mortar face
sideInfo(5,sA,eA) = DomainBoundaryId3D(nodeCoords(:,:,:,:,eA),sA,dx,tol,bcids)
enddo
enddo
! Global side ids: conforming pairs share one id; the big mortar face shares the
! id of sub-face 1; sub-faces 2..4 get their own ids (MPI message tags).
gid = 0
do eA = 1,nGlobalElem
do sA = 1,6
if(sideInfo(2,sA,eA) /= 0) cycle
if(eA == 1 .and. sA == selfSide3D_East) cycle ! assigned with the sub-faces
matched = .false.
do q = 1,4
if(eA == minfo(2*q+1) .and. sA == minfo(2*q+2)/10) then
matched = .true.
exit
endif
enddo
if(matched) cycle ! assigned with the sub-faces below
gid = gid+1
sideInfo(2,sA,eA) = gid
eB = sideInfo(3,sA,eA)
if(eB /= 0) then
sB = sideInfo(4,sA,eA)/10
sideInfo(2,sB,eB) = gid
endif
enddo
enddo
do q = 1,4
gid = gid+1
minfo(10+q) = gid
sideInfo(2,minfo(2*q+2)/10,minfo(2*q+1)) = gid
sideInfo(1,minfo(2*q+2)/10,minfo(2*q+1)) = 1 ! mortar index
if(q == 1) then
sideInfo(2,selfSide3D_East,1) = gid
sideInfo(1,selfSide3D_East,1) = 1 ! mortar index
endif
enddo
nUniqueSides = gid
! Domain decomposition. The message count upper bound is oversized relative to
! nUniqueSides to accommodate the per-variable (and per-direction) mortar and
! conforming side messages on this small mesh.
call this%decomp%GenerateDecomposition(nGlobalElem,64*nUniqueSides)
e1 = this%decomp%offsetElem(this%decomp%rankId+1)+1
e2 = this%decomp%offsetElem(this%decomp%rankId+2)
nLocalElems = e2-e1+1
call this%Init(nGeo,nLocalElems,nLocalElems*6,nLocalElems*8,nBCs)
this%nUniqueSides = nUniqueSides
this%quadrature = UNIFORM
this%BCType = 0
this%elemInfo = 0
this%globalNodeIDs = 0
this%nodeCoords(1:3,1:2,1:2,1:2,1:nLocalElems) = nodeCoords(1:3,1:2,1:2,1:2,e1:e2)
this%sideInfo(1:5,1:6,1:nLocalElems) = sideInfo(1:5,1:6,e1:e2)
! The mortar table is replicated on all ranks; element ids are global
this%nMortars = 1
allocate(this%mortarInfo(1:14,1:1))
this%mortarInfo(1:14,1) = minfo(1:14)
call this%UpdateDevice()
endsubroutine SimpleMortarMesh_Mesh3D_t