Create the smallest 3D anisotropic (2-subface) nonconforming mesh: one "big" element whose east face is shared with the faces of two "small" dx x dx x dx elements that split it along one face coordinate only. The mesh is conforming everywhere except at the single mortar interface.
With splitCoord = 1 the big element is [0,2dx] x [0,2dx] x [0,dx] and the small elements sit at y in [0,dx] and [dx,2dx]; the big east face's first trace coordinate (xi2, along y) is split. With splitCoord = 2 the big element is [0,2dx] x [0,dx] x [0,2dx] and the small elements are stacked in z, splitting the second trace coordinate (xi3). The small elements occupy x in [2dx,3dx].
y = 2dx ________________ ______
| | e3 | (splitCoord = 1, view along z; e2 and e3
| e1 |______| sit on the halves h = 1 and 2)
| | e2 |
|________________|______|
x = 0 2dx 3dx
Input - this : Fresh/empty Mesh3D_t object - dx : Edge length of the small elements (> 0) - 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:2) (optional) : requested face flip (0..7) for the small element on each half; each small element is rigidly rotated about its centre as in SimpleMortarMesh. Defaults to 0. - splitCoord (optional) : 1 or 2, the big-face trace coordinate that is split. Defaults to 1.
Output - this : Mesh3D_t object with three elements and one anisotropic mortar (nMortars = 0, nAnisoMortars = 1)
As in SimpleMortarMesh, all connectivity is derived by matching sampled points of the bilinear face maps, so sideInfo and anisoMortarInfo are correct by construction for every rotation. Element geometry is trilinear (nGeo=1). On two 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:2) | ||
| integer, | intent(in), | optional | :: | splitCoord |
subroutine SimpleAnisoMortarMesh_Mesh3D_t(this,dx,bcids,comm,flips,splitCoord)
!!
!! Create the smallest 3D anisotropic (2-subface) nonconforming mesh: one "big"
!! element whose east face is shared with the faces of two "small" dx x dx x dx
!! elements that split it along one face coordinate only. The mesh is conforming
!! everywhere except at the single mortar interface.
!!
!! With splitCoord = 1 the big element is [0,2dx] x [0,2dx] x [0,dx] and the small
!! elements sit at y in [0,dx] and [dx,2dx]; the big east face's first trace
!! coordinate (xi2, along y) is split. With splitCoord = 2 the big element is
!! [0,2dx] x [0,dx] x [0,2dx] and the small elements are stacked in z, splitting the
!! second trace coordinate (xi3). The small elements occupy x in [2dx,3dx].
!!
!! y = 2dx ________________ ______
!! | | e3 | (splitCoord = 1, view along z; e2 and e3
!! | e1 |______| sit on the halves h = 1 and 2)
!! | | e2 |
!! |________________|______|
!! x = 0 2dx 3dx
!!
!! Input
!! - this : Fresh/empty Mesh3D_t object
!! - dx : Edge length of the small elements (> 0)
!! - 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:2) (optional) : requested face flip (0..7) for the small element on
!! each half; each small element is rigidly rotated about its
!! centre as in SimpleMortarMesh. Defaults to 0.
!! - splitCoord (optional) : 1 or 2, the big-face trace coordinate that is split.
!! Defaults to 1.
!!
!! Output
!! - this : Mesh3D_t object with three elements and one anisotropic mortar
!! (nMortars = 0, nAnisoMortars = 1)
!!
!! As in SimpleMortarMesh, all connectivity is derived by matching sampled points of
!! the bilinear face maps, so sideInfo and anisoMortarInfo are correct by
!! construction for every rotation. Element geometry is trilinear (nGeo=1). On two
!! 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:2)
integer,intent(in),optional :: splitCoord
! Local
integer,parameter :: nGlobalElem = 3
! Signed local axis directions of the rotated small elements; see SimpleMortarMesh.
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
integer,parameter :: invFlip(0:7) = [0,1,2,3,4,7,6,5]
real(prec) :: nodeCoords(1:3,1:2,1:2,1:2,1:nGlobalElem)
integer :: sideInfo(1:5,1:6,1:nGlobalElem)
integer :: minfo(1:9)
integer :: localFlips(1:2),split
real(prec) :: center(1:3),axis(1:3,1:3),bigExtent(1:3),domainExtent(1:3)
real(prec) :: tol
integer :: h,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 h = 1,2
if(localFlips(h) < 0 .or. localFlips(h) > 7) then
print*,"SimpleAnisoMortarMesh_Mesh3D_t : requested flip out of range 0..7"
stop 1
endif
enddo
endif
split = 1
if(present(splitCoord)) split = splitCoord
if(split /= 1 .and. split /= 2) then
print*,"SimpleAnisoMortarMesh_Mesh3D_t : splitCoord must be 1 or 2"
stop 1
endif
! Element 1 (big) : 2dx along x and along the split direction, dx along the other
if(split == 1) then
bigExtent = [2.0_prec*dx,2.0_prec*dx,dx]
else
bigExtent = [2.0_prec*dx,dx,2.0_prec*dx]
endif
domainExtent = bigExtent
domainExtent(1) = 3.0_prec*dx
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)]*bigExtent
enddo
enddo
enddo
! Elements 2,3 (small) : [2dx,3dx] x half h of the big east face, each rigidly
! rotated per the requested flip
do h = 1,2
f = localFlips(h)
center(1) = 2.5_prec*dx
center(2) = 0.5_prec*dx
center(3) = 0.5_prec*dx
center(1+split) = (real(h,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+h) = 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
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: the big element's east face against each small element
minfo = 0
minfo(1) = 1
minfo(2) = selfSide3D_East
minfo(9) = split
do h = 1,2
matched = .false.
do sB = 1,6
do f = 0,7
if(HalfMatches3D(nodeCoords(:,:,:,:,1),selfSide3D_East,h,split, &
nodeCoords(:,:,:,:,1+h),sB,f,tol)) then
minfo(2*h+1) = 1+h
minfo(2*h+2) = 10*sB+f
matched = .true.
exit
endif
enddo
if(matched) exit
enddo
if(.not. matched) then
print*,"SimpleAnisoMortarMesh_Mesh3D_t : failed to match mortar half",h
stop 1
endif
enddo
! Physical boundary conditions on the faces that are neither conforming nor mortar
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 h = 1,2
if(eA == minfo(2*h+1) .and. sA == minfo(2*h+2)/10) matched = .true.
enddo
if(matched) cycle ! small mortar face
sideInfo(5,sA,eA) = DomainBoxBoundaryId3D(nodeCoords(:,:,:,:,eA),sA,domainExtent,tol,bcids)
enddo
enddo
! Global side ids: conforming pairs share one id; the big mortar face shares the
! id of half 1; the two halves 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
matched = .false.
do h = 1,2
if(eA == minfo(2*h+1) .and. sA == minfo(2*h+2)/10) matched = .true.
enddo
if(matched) cycle
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 h = 1,2
gid = gid+1
minfo(6+h) = gid
sideInfo(2,minfo(2*h+2)/10,minfo(2*h+1)) = gid
sideInfo(1,minfo(2*h+2)/10,minfo(2*h+1)) = 1 ! nMortars + anisotropic mortar index
if(h == 1) then
sideInfo(2,selfSide3D_East,1) = gid
sideInfo(1,selfSide3D_East,1) = 1
endif
enddo
nUniqueSides = gid
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%nAnisoMortars = 1
allocate(this%anisoMortarInfo(1:9,1:1))
this%anisoMortarInfo(1:9,1) = minfo(1:9)
call this%UpdateDevice()
endsubroutine SimpleAnisoMortarMesh_Mesh3D_t