Locate each stored physical point inside the 3D SEMHex geometry. On exit, elements(p) and coordinates(p,1:3) hold the located element id and reference (s,t,u) for points that resolve; elements(p) = 0 otherwise.
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| class(Points_t), | intent(inout) | :: | this | |||
| type(SEMHex), | intent(in) | :: | geometry |
subroutine LocatePoints_3D_Points_t(this,geometry)
!! Locate each stored physical point inside the 3D SEMHex geometry. On
!! exit, elements(p) and coordinates(p,1:3) hold the located element id
!! and reference (s,t,u) for points that resolve; elements(p) = 0
!! otherwise.
implicit none
class(Points_t),intent(inout) :: this
type(SEMHex),intent(in) :: geometry
! Local
integer :: p,iEl,iCand,nCand
integer :: nElem,N
integer :: ix,iy,iz,nCellsX,nCellsY,nCellsZ,nCellsTotal,cellId
real(prec) :: gMin(3),gMax(3),cellSize(3),slack,bbDiag
real(prec) :: padX,padY,padZ
real(prec) :: xTarget(3),xi(3)
real(prec),allocatable :: bbMin(:,:),bbMax(:,:)
integer,allocatable :: cellStart(:),cellElems(:)
integer :: totalEntries,ixLo,ixHi,iyLo,iyHi,izLo,izHi,k
real(prec) :: nC
logical :: converged
if(this%nDim /= 3) then
print*,"SELF_Points_t::LocatePoints (3D): nDim must be 3"
stop 1
endif
nElem = geometry%nElem
N = geometry%x%interp%N
this%elements = 0
this%coordinates = 0.0_prec
if(nElem == 0 .or. this%nPoints == 0) return
! --- Step 1: per-element AABBs ------------------------------------------
allocate(bbMin(1:3,1:nElem),bbMax(1:3,1:nElem))
call BuildElementBBoxes_3D(geometry,N,nElem,bbMin,bbMax)
! --- Step 2: global bbox + grid sizing ----------------------------------
gMin(1) = minval(bbMin(1,:))
gMin(2) = minval(bbMin(2,:))
gMin(3) = minval(bbMin(3,:))
gMax(1) = maxval(bbMax(1,:))
gMax(2) = maxval(bbMax(2,:))
gMax(3) = maxval(bbMax(3,:))
bbDiag = sqrt((gMax(1)-gMin(1))**2+(gMax(2)-gMin(2))**2+(gMax(3)-gMin(3))**2)
slack = max(bbDiag*1.0e-8_prec,tiny(1.0_prec)*1.0e6_prec)
gMin = gMin-slack
gMax = gMax+slack
nC = real(nElem,prec)*1.25_prec
nCellsX = max(1,ceiling(nC**(1.0_prec/3.0_prec)))
nCellsY = nCellsX
nCellsZ = nCellsX
nCellsTotal = nCellsX*nCellsY*nCellsZ
cellSize(1) = (gMax(1)-gMin(1))/real(nCellsX,prec)
cellSize(2) = (gMax(2)-gMin(2))/real(nCellsY,prec)
cellSize(3) = (gMax(3)-gMin(3))/real(nCellsZ,prec)
if(cellSize(1) <= 0.0_prec) cellSize(1) = 1.0_prec
if(cellSize(2) <= 0.0_prec) cellSize(2) = 1.0_prec
if(cellSize(3) <= 0.0_prec) cellSize(3) = 1.0_prec
! --- Step 3: CSR-style cell->elements hash ------------------------------
allocate(cellStart(0:nCellsTotal))
cellStart = 0
do iEl = 1,nElem
padX = max((bbMax(1,iEl)-bbMin(1,iEl))*1.0e-8_prec,slack)
padY = max((bbMax(2,iEl)-bbMin(2,iEl))*1.0e-8_prec,slack)
padZ = max((bbMax(3,iEl)-bbMin(3,iEl))*1.0e-8_prec,slack)
ixLo = ClampCell((bbMin(1,iEl)-padX-gMin(1))/cellSize(1),nCellsX)
ixHi = ClampCell((bbMax(1,iEl)+padX-gMin(1))/cellSize(1),nCellsX)
iyLo = ClampCell((bbMin(2,iEl)-padY-gMin(2))/cellSize(2),nCellsY)
iyHi = ClampCell((bbMax(2,iEl)+padY-gMin(2))/cellSize(2),nCellsY)
izLo = ClampCell((bbMin(3,iEl)-padZ-gMin(3))/cellSize(3),nCellsZ)
izHi = ClampCell((bbMax(3,iEl)+padZ-gMin(3))/cellSize(3),nCellsZ)
do iz = izLo,izHi
do iy = iyLo,iyHi
do ix = ixLo,ixHi
cellId = ix+nCellsX*(iy+nCellsY*iz)
cellStart(cellId+1) = cellStart(cellId+1)+1
enddo
enddo
enddo
enddo
do k = 1,nCellsTotal
cellStart(k) = cellStart(k)+cellStart(k-1)
enddo
totalEntries = cellStart(nCellsTotal)
allocate(cellElems(1:totalEntries))
do iEl = 1,nElem
padX = max((bbMax(1,iEl)-bbMin(1,iEl))*1.0e-8_prec,slack)
padY = max((bbMax(2,iEl)-bbMin(2,iEl))*1.0e-8_prec,slack)
padZ = max((bbMax(3,iEl)-bbMin(3,iEl))*1.0e-8_prec,slack)
ixLo = ClampCell((bbMin(1,iEl)-padX-gMin(1))/cellSize(1),nCellsX)
ixHi = ClampCell((bbMax(1,iEl)+padX-gMin(1))/cellSize(1),nCellsX)
iyLo = ClampCell((bbMin(2,iEl)-padY-gMin(2))/cellSize(2),nCellsY)
iyHi = ClampCell((bbMax(2,iEl)+padY-gMin(2))/cellSize(2),nCellsY)
izLo = ClampCell((bbMin(3,iEl)-padZ-gMin(3))/cellSize(3),nCellsZ)
izHi = ClampCell((bbMax(3,iEl)+padZ-gMin(3))/cellSize(3),nCellsZ)
do iz = izLo,izHi
do iy = iyLo,iyHi
do ix = ixLo,ixHi
cellId = ix+nCellsX*(iy+nCellsY*iz)
cellStart(cellId) = cellStart(cellId)+1
cellElems(cellStart(cellId)) = iEl
enddo
enddo
enddo
enddo
do k = nCellsTotal,1,-1
cellStart(k) = cellStart(k-1)
enddo
cellStart(0) = 0
! --- Step 4: per-point location -----------------------------------------
do p = 1,this%nPoints
xTarget(1) = this%x(p,1)
xTarget(2) = this%x(p,2)
xTarget(3) = this%x(p,3)
if(xTarget(1) < gMin(1) .or. xTarget(1) > gMax(1) .or. &
xTarget(2) < gMin(2) .or. xTarget(2) > gMax(2) .or. &
xTarget(3) < gMin(3) .or. xTarget(3) > gMax(3)) cycle
ix = ClampCell((xTarget(1)-gMin(1))/cellSize(1),nCellsX)
iy = ClampCell((xTarget(2)-gMin(2))/cellSize(2),nCellsY)
iz = ClampCell((xTarget(3)-gMin(3))/cellSize(3),nCellsZ)
cellId = ix+nCellsX*(iy+nCellsY*iz)
nCand = cellStart(cellId+1)-cellStart(cellId)
do iCand = 1,nCand
iEl = cellElems(cellStart(cellId)+iCand)
padX = max((bbMax(1,iEl)-bbMin(1,iEl))*1.0e-8_prec,slack)
padY = max((bbMax(2,iEl)-bbMin(2,iEl))*1.0e-8_prec,slack)
padZ = max((bbMax(3,iEl)-bbMin(3,iEl))*1.0e-8_prec,slack)
if(xTarget(1) < bbMin(1,iEl)-padX .or. xTarget(1) > bbMax(1,iEl)+padX) cycle
if(xTarget(2) < bbMin(2,iEl)-padY .or. xTarget(2) > bbMax(2,iEl)+padY) cycle
if(xTarget(3) < bbMin(3,iEl)-padZ .or. xTarget(3) > bbMax(3,iEl)+padZ) cycle
call NewtonInverse_3D(geometry,iEl,N,xTarget,xi,converged)
if(.not. converged) cycle
if(abs(xi(1)) <= 1.0_prec+1.0e-6_prec .and. &
abs(xi(2)) <= 1.0_prec+1.0e-6_prec .and. &
abs(xi(3)) <= 1.0_prec+1.0e-6_prec) then
this%elements(p) = iEl
this%coordinates(p,1) = min(1.0_prec,max(-1.0_prec,xi(1)))
this%coordinates(p,2) = min(1.0_prec,max(-1.0_prec,xi(2)))
this%coordinates(p,3) = min(1.0_prec,max(-1.0_prec,xi(3)))
exit
endif
enddo
enddo
! Fill per-point basis cache (see 2D variant for rationale).
if(associated(this%lS_cache)) deallocate(this%lS_cache)
if(associated(this%lT_cache)) deallocate(this%lT_cache)
if(associated(this%lU_cache)) deallocate(this%lU_cache)
this%lS_cache => null()
this%lT_cache => null()
this%lU_cache => null()
allocate(this%lS_cache(0:N,1:this%nPoints))
allocate(this%lT_cache(0:N,1:this%nPoints))
allocate(this%lU_cache(0:N,1:this%nPoints))
this%lS_cache = 0.0_prec
this%lT_cache = 0.0_prec
this%lU_cache = 0.0_prec
do p = 1,this%nPoints
if(this%elements(p) <= 0) cycle
this%lS_cache(:,p) = geometry%x%interp%CalculateLagrangePolynomials(this%coordinates(p,1))
this%lT_cache(:,p) = geometry%x%interp%CalculateLagrangePolynomials(this%coordinates(p,2))
this%lU_cache(:,p) = geometry%x%interp%CalculateLagrangePolynomials(this%coordinates(p,3))
enddo
this%nCached = N
deallocate(bbMin,bbMax,cellStart,cellElems)
endsubroutine LocatePoints_3D_Points_t