[Nek5000-users] findpts

[nek5000-users at lists.mcs.anl.gov]

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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/nek5000-users/attachments/20180909/a74ae995/attachment.html>