[Nek5000-users] findpts

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

Check out the hemi example. This will show you how to use the interpolation tool.

> On 9 Sep 2018, at 13:27, "nek5000-users at lists.mcs.anl.gov" <nek5000-users at lists.mcs.anl.gov> wrote:
> 
> 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
> 
> _______________________________________________
> 
> Nek5000-users mailing list
> 
> Nek5000-users at lists.mcs.anl.gov
> 
> https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users
>