Locate each stored physical point inside the 2D SEMQuad geometry. On exit, elements(p) holds the (rank-local) element id and coordinates(p,:) holds the (s,t) reference coordinates, both for points that resolve. Points that resolve to no element are marked with elements(p) = 0.
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| class(Points_t), | intent(inout) | :: | this | |||
| type(SEMQuad), | intent(in) | :: | geometry |
subroutine LocatePoints_2D_Points_t(this,geometry)
!! Locate each stored physical point inside the 2D SEMQuad geometry. On
!! exit, elements(p) holds the (rank-local) element id and coordinates(p,:)
!! holds the (s,t) reference coordinates, both for points that resolve.
!! Points that resolve to no element are marked with elements(p) = 0.
implicit none
class(Points_t),intent(inout) :: this
type(SEMQuad),intent(in) :: geometry
! Local
integer :: p,iEl,iCand,nCand
integer :: nElem,N
integer :: ix,iy,nCellsX,nCellsY,nCellsTotal,cellId
real(prec) :: gMin(2),gMax(2),cellSize(2),slack
real(prec) :: bbDiag,padX,padY
real(prec) :: xTarget(2),xi(2)
real(prec),allocatable :: bbMin(:,:),bbMax(:,:)
integer,allocatable :: cellStart(:),cellElems(:)
integer :: totalEntries,ixLo,ixHi,iyLo,iyHi,k
logical :: converged
if(this%nDim /= 2) then
print*,"SELF_Points_t::LocatePoints (2D): nDim must be 2"
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 axis-aligned bounding boxes --------------------
allocate(bbMin(1:2,1:nElem),bbMax(1:2,1:nElem))
call BuildElementBBoxes_2D(geometry,N,nElem,bbMin,bbMax)
! --- Step 2: global bbox + spatial-hash grid sizing ---------------------
gMin(1) = minval(bbMin(1,:))
gMin(2) = minval(bbMin(2,:))
gMax(1) = maxval(bbMax(1,:))
gMax(2) = maxval(bbMax(2,:))
bbDiag = sqrt((gMax(1)-gMin(1))**2+(gMax(2)-gMin(2))**2)
slack = max(bbDiag*1.0e-8_prec,tiny(1.0_prec)*1.0e6_prec)
gMin = gMin-slack
gMax = gMax+slack
! Aim for ~1 element per cell along each axis. Inflate by 25% to keep
! buckets short. Floor at 1.
nCellsX = max(1,ceiling(sqrt(real(nElem,prec)*1.25_prec)))
nCellsY = nCellsX
nCellsTotal = nCellsX*nCellsY
cellSize(1) = (gMax(1)-gMin(1))/real(nCellsX,prec)
cellSize(2) = (gMax(2)-gMin(2))/real(nCellsY,prec)
if(cellSize(1) <= 0.0_prec) cellSize(1) = 1.0_prec
if(cellSize(2) <= 0.0_prec) cellSize(2) = 1.0_prec
! --- Step 3: build CSR-style cell->elements hash ------------------------
allocate(cellStart(0:nCellsTotal))
cellStart = 0
! Pass 1: count entries per cell.
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)
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)
do iy = iyLo,iyHi
do ix = ixLo,ixHi
cellId = ix+iy*nCellsX
cellStart(cellId+1) = cellStart(cellId+1)+1
enddo
enddo
enddo
! Prefix sum -> cellStart now indexes into cellElems.
do k = 1,nCellsTotal
cellStart(k) = cellStart(k)+cellStart(k-1)
enddo
totalEntries = cellStart(nCellsTotal)
allocate(cellElems(1:totalEntries))
! Pass 2: place entries.
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)
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)
do iy = iyLo,iyHi
do ix = ixLo,ixHi
cellId = ix+iy*nCellsX
cellStart(cellId) = cellStart(cellId)+1
cellElems(cellStart(cellId)) = iEl
enddo
enddo
enddo
! Restore cellStart to start-offsets (currently holds end-offsets).
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)
! Quick reject: outside global bbox.
if(xTarget(1) < gMin(1) .or. xTarget(1) > gMax(1) .or. &
xTarget(2) < gMin(2) .or. xTarget(2) > gMax(2)) cycle
ix = ClampCell((xTarget(1)-gMin(1))/cellSize(1),nCellsX)
iy = ClampCell((xTarget(2)-gMin(2))/cellSize(2),nCellsY)
cellId = ix+iy*nCellsX
nCand = cellStart(cellId+1)-cellStart(cellId)
do iCand = 1,nCand
iEl = cellElems(cellStart(cellId)+iCand)
! Cheap AABB reject (with slack).
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)
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
call NewtonInverse_2D(geometry,iEl,N,xTarget,xi,converged)
if(.not. converged) cycle
! Accept if reference coords lie in the (slightly inflated) bi-unit cube.
if(abs(xi(1)) <= 1.0_prec+1.0e-6_prec .and. &
abs(xi(2)) <= 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)))
exit
endif
enddo
enddo
! Fill per-point Lagrange basis cache: the basis values are a function of
! the reference coordinates only, so callers that sample many times reuse
! these without recomputing CalculateLagrangePolynomials per call.
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))
this%lS_cache = 0.0_prec
this%lT_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))
enddo
this%nCached = N
deallocate(bbMin,bbMax,cellStart,cellElems)
endsubroutine LocatePoints_2D_Points_t