<div dir="ltr"><div dir="ltr"><div dir="ltr">Hi nek,<div><br></div><div>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.</div><div><br></div><div>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?</div><div><br></div><div><div>      subroutine interp_p(ielll,ixx,iyy,izz,xyz,n) !find elememt and processor for points xyz</div><div><br></div><div>      include 'SIZE'</div><div>      include 'TOTAL'</div><div><br></div><div>      real xyz(ldim,n), r(n), s(n), s(n)</div><div>      logical ifjac,ifpts</div><div><br></div><div>      integer rcodee(n), ielll(n), nidd(n)<span style="white-space:pre">      </span>  </div><div>      integer nfail, inn</div><div>      real dist(lpart)<span style="white-space:pre">       </span>  </div><div>      real rst(lpart*ldim)<span style="white-space:pre">     </span>  </div><div><br></div><div>      parameter(nmax=lpart,nfldmax=ldim) </div><div>      common /rv_intp/ pts(ldim*nmax)</div><div>      common /iv_intp/ ihandle</div><div><br></div><div>      integer icalld,e</div><div>      save    icalld</div><div>      data    icalld /0/</div><div><br></div><div>      nfail = 0</div><div>      nxyz  = nx1*ny1*nz1</div><div>      ntot  = nxyz*nelt</div><div><br></div><div>      if (n.gt.nmax) call exitti ('ABORT: interp_p() n > nmax!$',n)</div><div>      </div><div>      if (nelgt.ne.nelgv) call exitti</div><div>     $   ('ABORT: interp_p() nelgt.ne.nelgv not yet supported!$',nelgv)</div><div><br></div><div>      do i=1,n<span style="white-space:pre">                           </span>! ? not moving -> save?</div><div>         pts(i)     = xyz(1,i)</div><div>         pts(i + n) = xyz(2,i)</div><div>         if (if3d)  pts(i + n*2) = xyz(3,i)</div><div>      enddo</div><div><br></div><div>      if (icalld.eq.0) then<span style="white-space:pre">         </span>! interpolation setup<span style="white-space:pre">        </span>!? intpts_done(ih_intp_v)?</div><div>        icalld = 1</div><div>        tolin  = 1.e-8</div><div>        call intpts_setup(tolin,ihandle)</div><div>      endif</div><div><br></div><div>      nflds  = ndim ! number of fields to interpolate</div><div><br></div><div>      ! interpolate</div><div>      ifjac  = .true.           ! output transpose (of Jacobian)</div><div>      ifpts  = .true.            ! find points</div><div>      if(nio.eq.0) write(6,*) 'call findpts'</div><div>      call fintpts(ihandle,rcodee,1,ielll,1,rst,ndim,dist,1,pts(1),1,pts(n+1)</div><div>     $             ,1,pts(2*n+1),1,n)</div><div>      do inn = 1, n ! check return code</div><div>         if(rcodee(inn) .eq. 1) then</div><div>            nfail = nfail + 1</div><div>            if(nfail .le. 5) write(6,'(a,1p4e15.7)') </div><div>     $     ' WARNING: point on boundary or outside the mesh xy[z]d^2: ',</div><div>     $     (pts(inn+k*n),k=0,ndim-1),dist(inn)</div><div>            endif</div><div>         elseif(rcodee(inn) .eq. 2) then</div><div>            nfail = nfail + 1</div><div>            if(nfail .le. 5) write(6,'(a,1p3e15.7)') </div><div>     $     ' WARNING: point not within mesh xy[z]: !',</div><div>     $     (pts(inn+k*n),k=0,ndim-1)</div><div>         endif</div><div>      enddo</div></div><div><br></div><div><br></div><div>Kind regards,</div><div><br></div><div>Jian</div><div><br></div></div></div></div>