| 1 | #define INTP_HMAX 20
|
|---|
| 2 |
|
|---|
| 3 | c-----------------------------------------------------------------------
|
|---|
| 4 | subroutine interp_setup(ih,tolin,nmsh,nelm)
|
|---|
| 5 | c
|
|---|
| 6 | c input:
|
|---|
| 7 | c tolin ... tolerance newton solve (use 0 for default)
|
|---|
| 8 | c nmsh ... polynomial order for mesh (use 0 to fallback to lx1-1)
|
|---|
| 9 | c nelm ... number of local mesh elements (typically nelt)
|
|---|
| 10 | c
|
|---|
| 11 | c output:
|
|---|
| 12 | c ih ... handle
|
|---|
| 13 | c
|
|---|
| 14 |
|
|---|
| 15 | include 'SIZE'
|
|---|
| 16 | include 'INPUT'
|
|---|
| 17 | include 'GEOM'
|
|---|
| 18 |
|
|---|
| 19 | common /nekmpi/ nidd,npp,nekcomm,nekgroup,nekreal
|
|---|
| 20 |
|
|---|
| 21 | common /intp_h/ ih_intp(2,INTP_HMAX)
|
|---|
| 22 | common /intp/ tol
|
|---|
| 23 |
|
|---|
| 24 | data ihcounter /0/
|
|---|
| 25 | save ihcounter
|
|---|
| 26 |
|
|---|
| 27 | real xmi, ymi, zmi
|
|---|
| 28 | common /CBXMI/ xmi(lx1*ly1*lz1*lelt),
|
|---|
| 29 | $ ymi(lx1*ly1*lz1*lelt),
|
|---|
| 30 | $ zmi(lx1*ly1*lz1*lelt)
|
|---|
| 31 |
|
|---|
| 32 | real w(2*lx1**3)
|
|---|
| 33 |
|
|---|
| 34 | tol = max(5e-13,tolin)
|
|---|
| 35 | npt_max = 128
|
|---|
| 36 | bb_t = 0.01
|
|---|
| 37 |
|
|---|
| 38 | if (nio.eq.0)
|
|---|
| 39 | $ write(6,*) 'call intp_setup ','tol=', tol
|
|---|
| 40 |
|
|---|
| 41 | if (nmsh.ge.1 .and. nmsh.lt.lx1-1) then
|
|---|
| 42 | if (nio.eq.0) write(6,*) 'Nmsh for findpts:',nmsh
|
|---|
| 43 | nxi = nmsh+1
|
|---|
| 44 | nyi = nxi
|
|---|
| 45 | nzi = 1
|
|---|
| 46 | if (if3d) nzi = nxi
|
|---|
| 47 | do ie = 1,nelm
|
|---|
| 48 | call map_m_to_n(xmi((ie-1)*nxi*nyi*nzi + 1),nxi,xm1(1,1,1,ie)
|
|---|
| 49 | $ ,lx1,if3d,w,size(w))
|
|---|
| 50 | call map_m_to_n(ymi((ie-1)*nxi*nyi*nzi + 1),nyi,ym1(1,1,1,ie)
|
|---|
| 51 | $ ,ly1,if3d,w,size(w))
|
|---|
| 52 | if (if3d)
|
|---|
| 53 | $ call map_m_to_n(zmi((ie-1)*nxi*nyi*nzi + 1),nzi,zm1(1,1,1,ie)
|
|---|
| 54 | $ ,lz1,if3d,w,size(w))
|
|---|
| 55 | enddo
|
|---|
| 56 |
|
|---|
| 57 | n = nelm*nxi*nyi*nzi
|
|---|
| 58 | call fgslib_findpts_setup(ih_intp1,nekcomm,npp,ldim,
|
|---|
| 59 | & xm1,ym1,zm1,nx1,ny1,nz1,
|
|---|
| 60 | & nelm,nx1,ny1,nz1,bb_t,nelm+2,nelm+2,
|
|---|
| 61 | & npt_max,tol)
|
|---|
| 62 |
|
|---|
| 63 | call fgslib_findpts_setup(ih_intp2,nekcomm,npp,ldim,
|
|---|
| 64 | $ xmi,ymi,zmi,nxi,nyi,nzi,
|
|---|
| 65 | $ nelm,2*nxi,2*nyi,2*nzi,bb_t,n,n,
|
|---|
| 66 | $ npt_max,tol)
|
|---|
| 67 | else
|
|---|
| 68 | n = nelm*nx1*ny1*nz1
|
|---|
| 69 | call fgslib_findpts_setup(ih_intp1,nekcomm,npp,ldim,
|
|---|
| 70 | & xm1,ym1,zm1,nx1,ny1,nz1,
|
|---|
| 71 | & nelm,2*nx1,2*ny1,2*nz1,bb_t,n,n,
|
|---|
| 72 | & npt_max,tol)
|
|---|
| 73 | ih_intp2 = ih_intp1
|
|---|
| 74 | endif
|
|---|
| 75 |
|
|---|
| 76 | ihcounter = ihcounter + 1
|
|---|
| 77 | ih = ihcounter
|
|---|
| 78 | if (ih .gt. INTP_HMAX)
|
|---|
| 79 | $ call exitti('Maximum number of handles exceeded!$',INTP_HMAX)
|
|---|
| 80 | ih_intp(1,ih) = ih_intp1
|
|---|
| 81 | ih_intp(2,ih) = ih_intp2
|
|---|
| 82 |
|
|---|
| 83 | return
|
|---|
| 84 | end
|
|---|
| 85 | c-----------------------------------------------------------------------
|
|---|
| 86 | subroutine interp_nfld(out,fld,nfld,xp,yp,zp,n,iwk,rwk,nmax,
|
|---|
| 87 | $ iflp,ih)
|
|---|
| 88 | c
|
|---|
| 89 | c output:
|
|---|
| 90 | c out ... interpolation value(s) dim (n,nfld)
|
|---|
| 91 | c
|
|---|
| 92 | c input:
|
|---|
| 93 | c fld ... source field(s)
|
|---|
| 94 | c nfld ... number of fields
|
|---|
| 95 | c xp,yp,zp ... interpolation points
|
|---|
| 96 | c n ... number of points dim(xp,yp,zp)
|
|---|
| 97 | c iwk ... integer working array dim(nmax,3)
|
|---|
| 98 | c rwk ... real working array dim(nmax,ldim+1)
|
|---|
| 99 | c nmax ... leading dimension of iwk and rwk
|
|---|
| 100 | c iflp ... locate interpolation points (proc,el,r,s,t)
|
|---|
| 101 | c ih ... handle
|
|---|
| 102 | c
|
|---|
| 103 | include 'SIZE'
|
|---|
| 104 |
|
|---|
| 105 | common /intp/ tol
|
|---|
| 106 | common /intp_h/ ih_intp(2,INTP_HMAX)
|
|---|
| 107 |
|
|---|
| 108 | real fld(*),out(*)
|
|---|
| 109 | real xp(*),yp(*),zp(*)
|
|---|
| 110 |
|
|---|
| 111 | logical iflp
|
|---|
| 112 |
|
|---|
| 113 | real rwk(nmax,*)
|
|---|
| 114 | integer iwk(nmax,*)
|
|---|
| 115 |
|
|---|
| 116 | integer nn(2)
|
|---|
| 117 | logical ifot
|
|---|
| 118 |
|
|---|
| 119 | if (ih.lt.1 .or. ih.gt.INTP_HMAX)
|
|---|
| 120 | $ call exitti('invalid intp handle!$',ih)
|
|---|
| 121 |
|
|---|
| 122 | ih_intp1 = ih_intp(1,ih)
|
|---|
| 123 | ih_intp2 = ih_intp(2,ih)
|
|---|
| 124 |
|
|---|
| 125 | ifot = .false. ! transpose output field
|
|---|
| 126 |
|
|---|
| 127 | if (nio.eq.0 .and. loglevel.gt.2)
|
|---|
| 128 | $ write(6,*) 'call intp_nfld', ih, ih_intp1, ih_intp2
|
|---|
| 129 |
|
|---|
| 130 | if (n.gt.nmax) then
|
|---|
| 131 | write(6,*)
|
|---|
| 132 | & 'ABORT: n>nmax in intp_nfld', n, nmax
|
|---|
| 133 | call exitt
|
|---|
| 134 | endif
|
|---|
| 135 |
|
|---|
| 136 | ! locate points (iel,iproc,r,s,t)
|
|---|
| 137 | nfail = 0
|
|---|
| 138 | if(iflp) then
|
|---|
| 139 | if(nio.eq.0 .and. loglevel.gt.2) write(6,*) 'call findpts'
|
|---|
| 140 | call fgslib_findpts(ih_intp2,
|
|---|
| 141 | & iwk(1,1),1,
|
|---|
| 142 | & iwk(1,3),1,
|
|---|
| 143 | & iwk(1,2),1,
|
|---|
| 144 | & rwk(1,2),ldim,
|
|---|
| 145 | & rwk(1,1),1,
|
|---|
| 146 | & xp,1,
|
|---|
| 147 | & yp,1,
|
|---|
| 148 | & zp,1,n)
|
|---|
| 149 | do in=1,n
|
|---|
| 150 | ! check return code
|
|---|
| 151 | if(iwk(in,1).eq.1) then
|
|---|
| 152 | if(rwk(in,1).gt.10*tol) then
|
|---|
| 153 | nfail = nfail + 1
|
|---|
| 154 | if (nfail.le.5) write(6,'(a,1p4e15.7)')
|
|---|
| 155 | & ' WARNING: point on boundary or outside the mesh xy[z]d^2: ',
|
|---|
| 156 | & xp(in),yp(in),zp(in),rwk(in,1)
|
|---|
| 157 | endif
|
|---|
| 158 | elseif(iwk(in,1).eq.2) then
|
|---|
| 159 | nfail = nfail + 1
|
|---|
| 160 | if (nfail.le.5) write(6,'(a,1p3e15.7)')
|
|---|
| 161 | & ' WARNING: point not within mesh xy[z]: !',
|
|---|
| 162 | & xp(in),yp(in),zp(in)
|
|---|
| 163 | endif
|
|---|
| 164 | enddo
|
|---|
| 165 | endif
|
|---|
| 166 |
|
|---|
| 167 | ! evaluate inut field at given points
|
|---|
| 168 | ltot = lelt*lx1*ly1*lz1
|
|---|
| 169 | do ifld = 1,nfld
|
|---|
| 170 | iin = (ifld-1)*ltot + 1
|
|---|
| 171 | iout = (ifld-1)*n + 1
|
|---|
| 172 | is_out = 1
|
|---|
| 173 | if(ifot) then ! transpose output
|
|---|
| 174 | iout = ifld
|
|---|
| 175 | is_out = nfld
|
|---|
| 176 | endif
|
|---|
| 177 | call fgslib_findpts_eval(ih_intp1,out(iout),is_out,
|
|---|
| 178 | & iwk(1,1),1,
|
|---|
| 179 | & iwk(1,3),1,
|
|---|
| 180 | & iwk(1,2),1,
|
|---|
| 181 | & rwk(1,2),ldim,n,
|
|---|
| 182 | & fld(iin))
|
|---|
| 183 | enddo
|
|---|
| 184 |
|
|---|
| 185 | nn(1) = iglsum(n,1)
|
|---|
| 186 | nn(2) = iglsum(nfail,1)
|
|---|
| 187 | if(nio.eq.0) then
|
|---|
| 188 | if(nn(2).gt.0 .or. loglevel.gt.2) write(6,1) nn(1),nn(2)
|
|---|
| 189 | 1 format(' total number of points = ',i12,/,' failed = '
|
|---|
| 190 | & ,i12,/,' done :: intp_nfld')
|
|---|
| 191 | endif
|
|---|
| 192 |
|
|---|
| 193 | return
|
|---|
| 194 | end
|
|---|
| 195 | c-----------------------------------------------------------------------
|
|---|
| 196 | subroutine interp_free(ih)
|
|---|
| 197 |
|
|---|
| 198 | common /intp_h/ ih_intp(2,INTP_HMAX)
|
|---|
| 199 |
|
|---|
| 200 | ih_intp1 = ih_intp(1,ih)
|
|---|
| 201 | ih_intp2 = ih_intp(2,ih)
|
|---|
| 202 |
|
|---|
| 203 | call fgslib_findpts_free(ih_intp1)
|
|---|
| 204 | call fgslib_findpts_free(ih_intp2)
|
|---|
| 205 |
|
|---|
| 206 | return
|
|---|
| 207 | end
|
|---|
| 208 | c-----------------------------------------------------------------------
|
|---|
| 209 | subroutine intp_setup(tolin)
|
|---|
| 210 |
|
|---|
| 211 | call exitti('intp_setup is deprecated, see release notes!!$',1)
|
|---|
| 212 |
|
|---|
| 213 | return
|
|---|
| 214 | end
|
|---|
| 215 | c-----------------------------------------------------------------------
|
|---|
| 216 | subroutine intp_do(fldout,fldin,nfld,xp,yp,zp,n,iwk,rwk,nmax,iflp)
|
|---|
| 217 |
|
|---|
| 218 | call exitti('intp_do is deprecated, see release notes!!$',1)
|
|---|
| 219 |
|
|---|
| 220 | return
|
|---|
| 221 | end
|
|---|
| 222 | c-----------------------------------------------------------------------
|
|---|
| 223 | subroutine intp_free()
|
|---|
| 224 |
|
|---|
| 225 | call exitti('intp_free is deprecated, see release notes!!$',1)
|
|---|
| 226 |
|
|---|
| 227 | return
|
|---|
| 228 | end
|
|---|
