<div dir="ltr"><div dir="ltr"><div dir="ltr"><div dir="ltr">Dear Paul,<div><br></div><div>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). </div><div><br></div><div>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.</div><div><br></div><div>Would you please help me to look at the subroutines below? Thanks so much.</div><div><br></div><div>I just took the case with one point as example. Results are</div><div><br></div><div>call interp_pp(piel,pproc,pix,piy,piz,ppos0,lpart)</div><div><br></div><div> ppx:  -9.995177957692734E-002<br></div><div><div> pix:   0.523597350230634</div><div> ===</div><div> ppy:  -0.232351957858704</div><div> piy:   0.605225220747395</div><div> ===</div><div> ppz:    6.08294103920131</div><div> piz:   4.940656458412465E-324</div><div> ===</div><div> piel:   3.695611030892524E-321</div><div> pproc:   1.432790372939615E-322</div></div><div><br></div><div>*ppx ppy ppy are local point's coordinates to be interpolated  <br></div><div><br></div><div><div>C-----------------------------------------------------------------------</div><div>      subroutine interp_p(fieldin,nfld,pts,n,filedout,ifot,ifpts,ih,</div><div>     $                    elid,proc,rst)</div><div><br></div><div>c     based on intpts()</div><div><br></div><div><br></div><div>c     in:</div><div>c        fieldin ... input field(s) to interpolate (lelt*lxyz,nfld)</div><div>c        nfld    ... number of fields in fieldin</div><div>c        pts     ... packed list of interpolation points</div><div>c                    pts:=[x(1)...x(n),y(1)...y(n),z(1)...z(n)]</div><div>c        n       ... local number of interpolation points</div><div>c        ifto    ... transpose output (n,nfld)^T = (nfld,n) </div><div>c        itpts   ... find interpolation points  </div><div>c        ih      ... interpolation handle  </div><div>c     out:  </div><div>c        fieldout ... packed list of interpolated values (n,nfld)</div><div>c        elid     ... element on remote processor in which the point was found</div><div>c        proc     ... remote processor on which the point was found</div><div>c        rst      ... parametric coordinates for point</div><div><br></div><div>      include 'SIZE'</div><div>      include 'TOTAL'</div><div><br></div><div>      real fieldin(1), fieldout(1)</div><div>      real pts(1)</div><div><span style="white-space:pre">   </span>  </div><div>      real dist(lpart) ! squared distance</div><div>      real    rst(lpart*ldim)</div><div>      integer rcode(lpart),elid(lpart),proc(lpart)</div><div><br></div><div>      integer nn(2)</div><div>      logical ifot,ifpts</div><div>      ifjac  = .true.           ! output transpose (of Jacobian)</div><div>      ifpts  = .true.            ! find points</div><div>      if(nio.eq.0) write(6,*) 'call intpts in interp_p(1)'</div><div>      if(n.gt.lpart) then</div><div>        write(6,*)</div><div>     $   'ABORT: intpts() n>lpart, increase lpart in SIZE', n, lpart</div><div>        call exitt</div><div>        call exitt</div><div>      endif</div><div><br></div><div>c     ! locate points (iel,iproc,r,s,t)</div><div>      nfail = 0</div><div>      if(ifpts) then</div><div>         if(nio.eq.0) write(6,*) 'call findpts in interp_p(2)'</div><div>         call findpts(ih,rcode,1,</div><div>     $                  proc,1,</div><div>     $                  elid,1,</div><div>     $                  rst,ndim,</div><div>     $                  dist,1,</div><div>     $                  pts(    1),1,</div><div>     $                  pts(  n+1),1,</div><div>     $                  pts(2*n+1),1,n)</div><div>         do in=1,n ! check return code </div><div>         if(rcode(in).eq.1) then</div><div>             if(dist(in).gt.1e-12) 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(in+k*n),k=0,ndim-1),dist(in)</div><div>             endif</div><div>             elseif(rcode(in).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(in+k*n),k=0,ndim-1)</div><div>                 endif</div><div>         enddo</div><div>      endif</div><div><br></div><div>c     ! evaluate inut field at given points</div><div>      ltot = lelt*lx1*ly1*lz1</div><div>      do ifld = 1,nfld</div><div>         iin    = (ifld-1)*ltot + 1</div><div>         iout   = (ifld-1)*n + 1</div><div>         is_out = 1</div><div>         if(ifot) then ! transpose output</div><div>             iout   = ifld</div><div>             is_out = nfld </div><div>         endif</div><div>         call findpts_eval(ih,fieldout(iout),is_out,</div><div>     $                         rcode,1,</div><div>     $                     proc,1,</div><div>     $                     elid,1,</div><div>     $                     rst,ndim,n,</div><div>     $                     fieldin(iin))</div><div>      enddo</div><div>      nn(1) = iglsum(n,1)</div><div>      nn(2) = iglsum(nfail,1)</div><div>      if(nio.eq.0) then</div><div>         write(6,1) nn(1),nn(2)</div><div>1     format('   total number of points = ',i12,/,'   failed = '</div><div>     $     ,i12,/,' done :: intpts')</div><div>      endif</div><div>      return</div><div>      end</div><div>c-----------------------------------------------------------------------</div><div>      subroutine interp_pp(elidp,procp,rr,ss,tt,xyz,n) !evaluate elid,rst for xyz</div><div><br></div><div>c     based on interp_v()</div><div><br></div><div>      include 'SIZE'</div><div>      include 'TOTAL'</div><div><br></div><div>      real    rstp(lpart*ldim),rr(lpart),ss(lpart),tt(lpart)</div><div>      integer elidp(lpart),procp(lpart)</div><div>      real xyz(ldim,n)</div><div>      logical ifjac,ifpts</div><div><br></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>      common /outtmp/ wrk(lx1*ly1*lz1*lelt,nfldmax)</div><div><br></div><div>      integer icalld,e</div><div>      save    icalld</div><div>      data    icalld /0/</div><div><br></div><div>      nxyz  = nx1*ny1*nz1</div><div>      ntot  = nxyz*nelt</div><div><br></div><div>      if (n.gt.nmax) call exitti ('ABORT: interp_v() n > nmax!$',n)</div><div>      </div><div>      if (nelgt.ne.nelgv) call exitti</div><div>     $   ('ABORT: interp_v() 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>      ! pack working array</div><div>      call opcopy(wrk(1,1),wrk(1,2),wrk(1,3),vx,vy,vz)</div><div><br></div><div>      ! interpolate</div><div>      ifjac  = .true.           ! output transpose (of Jacobian)</div><div>      ifpts  = .true.            ! find points</div><div>      call interp_p(wrk,nflds,pts,n,uvw,ifjac,ifpts,ihandle</div><div>     $             ,elidp,procp,rstp)</div><div><span style="white-space:pre">      </span> </div><div>      do i=1,n</div><div>         rr(i) = rstp(i+n)</div><div>         ss(i) = rstp(i+2*n)</div><div>         tt(i) = rstp(i+3*n)</div><div>      enddo</div><div><br></div><div>      return</div><div>      end</div></div><div><br></div><div><br></div><div>Kind regards,</div><div><br></div><div>Jian</div></div></div></div></div>