[Nek5000-users] findpts

Fri Sep 14 16:03:31 CDT 2018

Dear Paul,

When I looked at the interp_v in hemi case, I didn't see it returns the
element number, processor number, and ix, iy and iz (ranging from 1 to nx1,
ny1, and nz1, respectively).

So for that purpose, I worte two subroutines which are similar to interp_v
and intpts to return the data above. But the results show those data are
not reasonable.

Would you please help me to look at the subroutines below? Thanks so much.

I just took the case with one point as example. Results are

call interp_pp(piel,pproc,pix,piy,piz,ppos0,lpart)

 ppx:  -9.995177957692734E-002
 pix:   0.523597350230634
 ===
 ppy:  -0.232351957858704
 piy:   0.605225220747395
 ===
 ppz:    6.08294103920131
 piz:   4.940656458412465E-324
 ===
 piel:   3.695611030892524E-321
 pproc:   1.432790372939615E-322

*ppx ppy ppy are local point's coordinates to be interpolated

C-----------------------------------------------------------------------
      subroutine interp_p(fieldin,nfld,pts,n,filedout,ifot,ifpts,ih,
     $                    elid,proc,rst)

c     based on intpts()

c     in:
c        fieldin ... input field(s) to interpolate (lelt*lxyz,nfld)
c        nfld    ... number of fields in fieldin
c        pts     ... packed list of interpolation points
c                    pts:=[x(1)...x(n),y(1)...y(n),z(1)...z(n)]
c        n       ... local number of interpolation points
c        ifto    ... transpose output (n,nfld)^T = (nfld,n)
c        itpts   ... find interpolation points
c        ih      ... interpolation handle
c     out:
c        fieldout ... packed list of interpolated values (n,nfld)
c        elid     ... element on remote processor in which the point was
found
c        proc     ... remote processor on which the point was found
c        rst      ... parametric coordinates for point

      include 'SIZE'
      include 'TOTAL'

      real fieldin(1), fieldout(1)
      real pts(1)

      real dist(lpart) ! squared distance
      real    rst(lpart*ldim)
      integer rcode(lpart),elid(lpart),proc(lpart)

      integer nn(2)
      logical ifot,ifpts
      ifjac  = .true.           ! output transpose (of Jacobian)
      ifpts  = .true.            ! find points
      if(nio.eq.0) write(6,*) 'call intpts in interp_p(1)'
      if(n.gt.lpart) then
        write(6,*)
     $   'ABORT: intpts() n>lpart, increase lpart in SIZE', n, lpart
        call exitt
        call exitt
      endif

c     ! locate points (iel,iproc,r,s,t)
      nfail = 0
      if(ifpts) then
         if(nio.eq.0) write(6,*) 'call findpts in interp_p(2)'
         call findpts(ih,rcode,1,
     $                  proc,1,
     $                  elid,1,
     $                  rst,ndim,
     $                  dist,1,
     $                  pts(    1),1,
     $                  pts(  n+1),1,
     $                  pts(2*n+1),1,n)
         do in=1,n ! check return code
         if(rcode(in).eq.1) then
             if(dist(in).gt.1e-12) 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(in+k*n),k=0,ndim-1),dist(in)
             endif
             elseif(rcode(in).eq.2) then
                 nfail = nfail + 1
                 if (nfail.le.5) write(6,'(a,1p3e15.7)')
     $  ' WARNING: point not within mesh xy[z]: !',
     $  (pts(in+k*n),k=0,ndim-1)
                 endif
         enddo
      endif

c     ! evaluate inut field at given points
      ltot = lelt*lx1*ly1*lz1
      do ifld = 1,nfld
         iin    = (ifld-1)*ltot + 1
         iout   = (ifld-1)*n + 1
         is_out = 1
         if(ifot) then ! transpose output
             iout   = ifld
             is_out = nfld
         endif
         call findpts_eval(ih,fieldout(iout),is_out,
     $                         rcode,1,
     $                     proc,1,
     $                     elid,1,
     $                     rst,ndim,n,
     $                     fieldin(iin))
      enddo
      nn(1) = iglsum(n,1)
      nn(2) = iglsum(nfail,1)
      if(nio.eq.0) then
         write(6,1) nn(1),nn(2)
1     format('   total number of points = ',i12,/,'   failed = '
     $     ,i12,/,' done :: intpts')
      endif
      return
      end
c-----------------------------------------------------------------------
      subroutine interp_pp(elidp,procp,rr,ss,tt,xyz,n) !evaluate elid,rst
for xyz

c     based on interp_v()

      include 'SIZE'
      include 'TOTAL'

      real    rstp(lpart*ldim),rr(lpart),ss(lpart),tt(lpart)
      integer elidp(lpart),procp(lpart)
      real xyz(ldim,n)
      logical ifjac,ifpts

      parameter(nmax=lpart,nfldmax=ldim)
      common /rv_intp/ pts(ldim*nmax)
      common /iv_intp/ ihandle
      common /outtmp/ wrk(lx1*ly1*lz1*lelt,nfldmax)

      integer icalld,e
      save    icalld
      data    icalld /0/

      nxyz  = nx1*ny1*nz1
      ntot  = nxyz*nelt

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

      if (nelgt.ne.nelgv) call exitti
     $   ('ABORT: interp_v() 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

      ! pack working array
      call opcopy(wrk(1,1),wrk(1,2),wrk(1,3),vx,vy,vz)

      ! interpolate
      ifjac  = .true.           ! output transpose (of Jacobian)
      ifpts  = .true.            ! find points
      call interp_p(wrk,nflds,pts,n,uvw,ifjac,ifpts,ihandle
     $             ,elidp,procp,rstp)

      do i=1,n
         rr(i) = rstp(i+n)
         ss(i) = rstp(i+2*n)
         tt(i) = rstp(i+3*n)
      enddo

      return
      end

Kind regards,

Jian
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/nek5000-users/attachments/20180914/2fa78c5e/attachment.html>