Fill mesh%nodeCoords(:,:,:,:,e) for one hexahedral element by transfinite (Coons) interpolation of its six face grids. Faces flagged in the mesh file use the file's face-point grids; the remaining faces are bilinear patches of their four corner nodes evaluated at Chebyshev-Gauss-Lobatto parametric coordinates. Edge curves are extracted from the face grids (preferring a flagged face when an edge borders one flagged and one bilinear face) and the standard Boolean-sum formula x = Pxi + Peta + Pzeta - PxiPeta - PxiPzeta - PetaPzeta + PxiPeta*Pzeta combines face, edge, and corner contributions. For an element with all-straight faces this reduces to trilinear interpolation of the eight corners.
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| class(Mesh3D_t), | intent(inout) | :: | mesh | |||
| integer, | intent(in) | :: | e | |||
| integer, | intent(in) | :: | nGeo | |||
| integer, | intent(in) | :: | allCorners(:,:) | |||
| integer, | intent(in) | :: | allFlags(:,:) | |||
| real(kind=prec), | intent(in) | :: | allFaces(:,:,:,:,:) | |||
| real(kind=prec), | intent(in) | :: | nodeXYZ(:,:) |
subroutine build_nodeCoords_for_hex(mesh,e,nGeo,allCorners,allFlags,allFaces,nodeXYZ)
!! Fill mesh%nodeCoords(:,:,:,:,e) for one hexahedral element by
!! transfinite (Coons) interpolation of its six face grids. Faces
!! flagged in the mesh file use the file's face-point grids; the
!! remaining faces are bilinear patches of their four corner nodes
!! evaluated at Chebyshev-Gauss-Lobatto parametric coordinates.
!! Edge curves are extracted from the face grids (preferring a
!! flagged face when an edge borders one flagged and one bilinear
!! face) and the standard Boolean-sum formula
!! x = Pxi + Peta + Pzeta - Pxi*Peta - Pxi*Pzeta - Peta*Pzeta
!! + Pxi*Peta*Pzeta
!! combines face, edge, and corner contributions. For an element
!! with all-straight faces this reduces to trilinear interpolation
!! of the eight corners.
implicit none
class(Mesh3D_t),intent(inout) :: mesh
integer,intent(in) :: e
integer,intent(in) :: nGeo
integer,intent(in) :: allCorners(:,:) ! 8 x nElem corner-node ids
integer,intent(in) :: allFlags(:,:) ! 6 x nElem face flags (SELF side order)
real(prec),intent(in) :: allFaces(:,:,:,:,:) ! 3 x nGeo+1 x nGeo+1 x 6 x nElem
real(prec),intent(in) :: nodeXYZ(:,:)
! Local
integer :: ng1
integer :: i,j,k,s,l
integer :: cornerIDs(1:8)
integer :: flag(1:6)
real(prec) :: P(1:3,1:8)
real(prec) :: face(1:3,1:nGeo+1,1:nGeo+1,1:6)
real(prec) :: exEdge(1:3,1:nGeo+1,1:4) ! xi-directed edges: bs, bn, ts, tn
real(prec) :: eyEdge(1:3,1:nGeo+1,1:4) ! eta-directed edges: bw, be, tw, te
real(prec) :: ezEdge(1:3,1:nGeo+1,1:4) ! zeta-directed edges: sw, se, nw, ne
real(prec) :: xi(0:nGeo),wts(0:nGeo)
real(prec) :: u(1:nGeo+1)
real(prec) :: a,b,c
real(prec) :: xfp(1:3),xep(1:3),xcp(1:3)
ng1 = nGeo+1
cornerIDs = allCorners(1:8,e)
flag = allFlags(1:6,e)
do l = 1,8
P(1:3,l) = nodeXYZ(1:3,cornerIDs(l))
enddo
if(nGeo == 1) then
! Pure trilinear element with the 8 corners
do l = 1,8
i = mesh%CGNSCornerMap(1,l)
j = mesh%CGNSCornerMap(2,l)
k = mesh%CGNSCornerMap(3,l)
mesh%nodeCoords(1:3,i,j,k,e) = P(1:3,l)
enddo
return
endif
! Chebyshev-Gauss-Lobatto parametric coordinates on [-1,1] and
! their [0,1] image used in the blending weights. We use SELF's
! quadrature routine to keep node placement consistent with the
! rest of the code.
call ChebyshevQuadrature(nGeo,xi,wts,CHEBYSHEV_GAUSS_LOBATTO)
do i = 1,ng1
u(i) = 0.5_prec*(xi(i-1)+1.0_prec)
enddo
! Face grids: curved faces from the file, straight faces as
! bilinear patches of their four corners (face-index corner order
! per sideMap; all faces use the natural on-face volume axes).
do s = 1,6
if(flag(s) == 1) then
face(1:3,1:ng1,1:ng1,s) = allFaces(1:3,1:ng1,1:ng1,s,e)
else
do j = 1,ng1
do i = 1,ng1
a = u(i)
b = u(j)
face(1:3,i,j,s) = (1.0_prec-a)*(1.0_prec-b)*P(1:3,mesh%sideMap(1,s))+ &
a*(1.0_prec-b)*P(1:3,mesh%sideMap(2,s))+ &
a*b*P(1:3,mesh%sideMap(3,s))+ &
(1.0_prec-a)*b*P(1:3,mesh%sideMap(4,s))
enddo
enddo
endif
enddo
! Edge curves extracted from the face grids. Each edge borders two
! faces; a flagged (curved) face is preferred so that straight-face
! bilinear patches never override curved edge data. When both
! bordering faces are straight, the two restrictions coincide (the
! straight line between the edge's corner nodes). SELF side ids:
! 1=bottom, 2=south, 3=east, 4=north, 5=west, 6=top.
do i = 1,ng1
! xi-directed edges
if(flag(1) == 1 .or. flag(2) /= 1) then ! bottom-south: bottom or south
exEdge(1:3,i,1) = face(1:3,i,1,1)
else
exEdge(1:3,i,1) = face(1:3,i,1,2)
endif
if(flag(1) == 1 .or. flag(4) /= 1) then ! bottom-north: bottom or north
exEdge(1:3,i,2) = face(1:3,i,ng1,1)
else
exEdge(1:3,i,2) = face(1:3,i,1,4)
endif
if(flag(6) == 1 .or. flag(2) /= 1) then ! top-south: top or south
exEdge(1:3,i,3) = face(1:3,i,1,6)
else
exEdge(1:3,i,3) = face(1:3,i,ng1,2)
endif
if(flag(6) == 1 .or. flag(4) /= 1) then ! top-north: top or north
exEdge(1:3,i,4) = face(1:3,i,ng1,6)
else
exEdge(1:3,i,4) = face(1:3,i,ng1,4)
endif
! eta-directed edges
if(flag(1) == 1 .or. flag(5) /= 1) then ! bottom-west: bottom or west
eyEdge(1:3,i,1) = face(1:3,1,i,1)
else
eyEdge(1:3,i,1) = face(1:3,i,1,5)
endif
if(flag(1) == 1 .or. flag(3) /= 1) then ! bottom-east: bottom or east
eyEdge(1:3,i,2) = face(1:3,ng1,i,1)
else
eyEdge(1:3,i,2) = face(1:3,i,1,3)
endif
if(flag(6) == 1 .or. flag(5) /= 1) then ! top-west: top or west
eyEdge(1:3,i,3) = face(1:3,1,i,6)
else
eyEdge(1:3,i,3) = face(1:3,i,ng1,5)
endif
if(flag(6) == 1 .or. flag(3) /= 1) then ! top-east: top or east
eyEdge(1:3,i,4) = face(1:3,ng1,i,6)
else
eyEdge(1:3,i,4) = face(1:3,i,ng1,3)
endif
! zeta-directed edges
if(flag(2) == 1 .or. flag(5) /= 1) then ! south-west: south or west
ezEdge(1:3,i,1) = face(1:3,1,i,2)
else
ezEdge(1:3,i,1) = face(1:3,1,i,5)
endif
if(flag(2) == 1 .or. flag(3) /= 1) then ! south-east: south or east
ezEdge(1:3,i,2) = face(1:3,ng1,i,2)
else
ezEdge(1:3,i,2) = face(1:3,1,i,3)
endif
if(flag(4) == 1 .or. flag(5) /= 1) then ! north-west: north or west
ezEdge(1:3,i,3) = face(1:3,1,i,4)
else
ezEdge(1:3,i,3) = face(1:3,ng1,i,5)
endif
if(flag(4) == 1 .or. flag(3) /= 1) then ! north-east: north or east
ezEdge(1:3,i,4) = face(1:3,ng1,i,4)
else
ezEdge(1:3,i,4) = face(1:3,ng1,i,3)
endif
enddo
! Boolean-sum transfinite interpolation. u,v,w in [0,1] are the
! blending parameters along xi, eta, zeta.
do k = 1,ng1
c = u(k)
do j = 1,ng1
b = u(j)
do i = 1,ng1
a = u(i)
! Face contributions (Pxi + Peta + Pzeta)
xfp = (1.0_prec-a)*face(1:3,j,k,5)+a*face(1:3,j,k,3)+ &
(1.0_prec-b)*face(1:3,i,k,2)+b*face(1:3,i,k,4)+ &
(1.0_prec-c)*face(1:3,i,j,1)+c*face(1:3,i,j,6)
! Edge corrections (Pxi*Peta + Pxi*Pzeta + Peta*Pzeta)
xep = (1.0_prec-a)*(1.0_prec-b)*ezEdge(1:3,k,1)+ &
a*(1.0_prec-b)*ezEdge(1:3,k,2)+ &
a*b*ezEdge(1:3,k,4)+ &
(1.0_prec-a)*b*ezEdge(1:3,k,3)+ &
(1.0_prec-a)*(1.0_prec-c)*eyEdge(1:3,j,1)+ &
a*(1.0_prec-c)*eyEdge(1:3,j,2)+ &
a*c*eyEdge(1:3,j,4)+ &
(1.0_prec-a)*c*eyEdge(1:3,j,3)+ &
(1.0_prec-b)*(1.0_prec-c)*exEdge(1:3,i,1)+ &
b*(1.0_prec-c)*exEdge(1:3,i,2)+ &
b*c*exEdge(1:3,i,4)+ &
(1.0_prec-b)*c*exEdge(1:3,i,3)
! Corner contributions (Pxi*Peta*Pzeta)
xcp = (1.0_prec-a)*(1.0_prec-b)*(1.0_prec-c)*P(1:3,1)+ &
a*(1.0_prec-b)*(1.0_prec-c)*P(1:3,2)+ &
a*b*(1.0_prec-c)*P(1:3,3)+ &
(1.0_prec-a)*b*(1.0_prec-c)*P(1:3,4)+ &
(1.0_prec-a)*(1.0_prec-b)*c*P(1:3,5)+ &
a*(1.0_prec-b)*c*P(1:3,6)+ &
a*b*c*P(1:3,7)+ &
(1.0_prec-a)*b*c*P(1:3,8)
mesh%nodeCoords(1:3,i,j,k,e) = xfp-xep+xcp
enddo
enddo
enddo
endsubroutine build_nodeCoords_for_hex