Hi nek,

I wrote before in terms of findpts subroutine and tried some code. Unfortunately, I failed to get it right. May I attach my code here and please help me to sort it.

I want to map the given set of arbitrary points xyz(ldim,n) below to the local coordinates in the mesh (velocity grid). Currently, I stop at how to use findpts results to return ix, iy, iz and iel (element number). Can anyone help me finish it?

subroutine interp_p(ielll,ixx,iyy,izz,xyz,n) !find elememt and processor for points xyz

include 'SIZE'

include 'TOTAL'

real xyz(ldim,n), r(n), s(n), s(n)

logical ifjac,ifpts

integer rcodee(n), ielll(n), nidd(n)

integer nfail, inn

real dist(lpart)

real rst(lpart*ldim)

parameter(nmax=lpart,nfldmax=ldim)

common /rv_intp/ pts(ldim*nmax)

common /iv_intp/ ihandle

integer icalld,e

save icalld

data icalld /0/

nfail = 0

nxyz = nx1*ny1*nz1

ntot = nxyz*nelt

if (n.gt.nmax) call exitti ('ABORT: interp_p() n > nmax!$',n)

if (nelgt.ne.nelgv) call exitti

$ ('ABORT: interp_p() nelgt.ne.nelgv not yet supported!$',nelgv)

do i=1,n ! ? not moving -> save?

pts(i) = xyz(1,i)

pts(i + n) = xyz(2,i)

if (if3d) pts(i + n*2) = xyz(3,i)

enddo

if (icalld.eq.0) then ! interpolation setup !? intpts_done(ih_intp_v)?

icalld = 1

tolin = 1.e-8

call intpts_setup(tolin,ihandle)

endif

nflds = ndim ! number of fields to interpolate

! interpolate

ifjac = .true. ! output transpose (of Jacobian)

ifpts = .true. ! find points

if(nio.eq.0) write(6,*) 'call findpts'

call fintpts(ihandle,rcodee,1,ielll,1,rst,ndim,dist,1,pts(1),1,pts(n+1)

$ ,1,pts(2*n+1),1,n)

do inn = 1, n ! check return code

if(rcodee(inn) .eq. 1) then

nfail = nfail + 1

if(nfail .le. 5) write(6,'(a,1p4e15.7)')

$ ' WARNING: point on boundary or outside the mesh xy[z]d^2: ',

$ (pts(inn+k*n),k=0,ndim-1),dist(inn)

endif

elseif(rcodee(inn) .eq. 2) then

nfail = nfail + 1

if(nfail .le. 5) write(6,'(a,1p3e15.7)')

$ ' WARNING: point not within mesh xy[z]: !',

$ (pts(inn+k*n),k=0,ndim-1)

endif

enddo

Kind regards,

Jian