<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>