1 subroutine renum_cyl(x) 2c 3c---------------------------------------------------------------------- 4c This subroutine finds all nodes that are on the inlet plane and all 5c nodes that are on the recycle plane; it also blocks elements that 6c are on recycle plane; all nodes are also stored with cylindrical 7c coordinates; find all father nodes for recycle plane (i.e. for 8c theta = -theta given) 9c 10c input: 11c x (numnp,nsd) : node coordinates 12c 13c output: 14c xcyl (numnp,nsd) : node cylindrical coordinates 15c ien2D (npro, nshl) : connectivity array for recycle plane 16c assuming tethraheadral elements, i.e. 17c triangular elements on face 18c 19c 20c Elaine Bohr 21c June 2002 22c---------------------------------------------------------------------- 23c 24 use spebc 25c use pointer_data 26 include "common.h" 27 include "mpif.h" 28 include "auxmpi.h" 29c 30 dimension x(numnp,nsd), nrin(numnp), nula(numnp), 31 & erreur(nshl), xtmp(nsd) 32 33 integer temp, etmp, s 34 35 36 real*8, allocatable :: xrtmp(:) 37 38 39c thetag = thetag/180.0*pi 40 41 42c 43c .... changing to cylindrical coordinate system for nodal point 44c 45 46 xcyl(:,1) = sqrt(x(:,1)*x(:,1) + x(:,2)*x(:,2)) 47 48 j = 0 49 do i=1,numnp 50 if ((x(i,1).eq.0).and.(x(i,2).eq.0)) then 51 j = j+1 52 nula(j) = i 53 else 54 xcyl(i,2) = atan2(x(i,2),x(i,1)) 55 endif 56 enddo 57 xcyl(:,3) = x(:,3) 58 59 60c 61c .... finding the minimum and maximum angles 62c 63 64 thmin = xcyl(nen1(1),2) 65 thmax = xcyl(nen1(1),2) 66 do i=2,npin 67 if((x(nen1(i),1).eq.0).and.(x(nen1(i),2).eq.0)) then 68 goto 10 69 else 70 if(thmin.gt.xcyl(nen1(i),2)) thmin = xcyl(nen1(i),2) 71 if(thmax.lt.xcyl(nen1(i),2)) thmax = xcyl(nen1(i),2) 72 endif 73 10 continue 74 enddo 75 76 do i=1,j 77 xcyl(nula(i),2) = thmin 78 enddo 79 80 81c 82c .... Finding nodes from inlet plane with same theta given 83c 84 85 j = 0 86 do i=1,npin 87 error1 = abs(xcyl(nen1(i),2) - thmin) 88 if(error1.le.0.00001) then 89 j = j + 1 90 nrin(j)=nen1(i) 91 endif 92 enddo 93 nfin = j 94 allocate(xrtmp(nfin)) 95 do i=1,nfin 96 xrtmp(i) = xcyl(nrin(i),1) 97 enddo 98 99c 100c .... Ordering nrin by decreasing radius 101c 102 j = 0 103 allocate (nrint(nfin)) 104 105 radinl = xcyl(nrin(1),1) 106 do i=1,nfin 107 if (radinl.lt.xcyl(nrin(i),1)) radinl = xcyl(nrin(i),1) 108 rmaxtemp=xrtmp(1) 109 itmp = 1 110 do k=2,nfin 111 if (rmaxtemp.le.xrtmp(k)) then 112 rmaxtemp = xrtmp(k) 113 itmp = k 114 endif 115 enddo 116 if (radcyl.ge.xrtmp(itmp)) then 117 j = j + 1 118 nrint(j)=nrin(itmp) 119 endif 120 xrtmp(itmp) = -1 121 enddo 122 nfint = j 123 124 deallocate(xrtmp) 125 126c 127c .... off wall coordinate for the inlet plane 128c 129 do i=1,npin 130 xynin(i) = (radinl - xcyl(nen1(i),1))/sang 131 enddo 132 133c 134c .... off wall coordinate for the virtual points on recycle plane 135c 136 137 do i=1,nfint 138 xyn(i) = (radcyl - xcyl(nrint(i),1))/sang 139 enddo 140 141c 142c .... Finding corresponding elements and local coordinates 143c on recycle plane for every arc with radius from nrint 144c 145 s = (thmax-thmin)*radcyl/ds 146 allocate (xsinfin(nfint,s+1,nsd)) 147 allocate (elcnfin(nfint,s+1)) 148 allocate (imax(nfint)) 149 150 do jj = 1, nfint 151 152 xts1 = x(nrint(jj),1) !*cos(xcyl(nrint(jj),2)+tolerence) 153 xts2 = x(nrint(jj),2) !*sin(xcyl(nrint(jj),2)+tolerence) 154 xts3 = x(nrint(jj),3) + plandist 155 call elem_search(xintl, xts1, xts2, xts3, 156 & xtmp(:), etmp, 1) 157 xsinfin(jj,1,:) = xtmp(:) 158 elcnfin(jj,1) = etmp 159 imax(jj) = (thmax-thmin)*xcyl(nrint(jj),1)/ds 160 do i=1,imax(jj) 161 if ( xcyl(nrint(jj),1) .eq. radcyl) then 162 xts1 = (xcyl(nrint(jj),1)-tolerence) 163 & *cos(1.0*i/imax(jj)*(thmax-thmin)+thmin) 164 xts2 = (xcyl(nrint(jj),1)-tolerence) 165 & *sin(1.0*i/imax(jj)*(thmax-thmin)+thmin) 166 else 167c reel=i*ds/xcyl(nrint(jj),1) 168 xts1 = xcyl(nrint(jj),1) 169 & *cos(1.0*i/imax(jj)*(thmax-thmin)+thmin) 170 xts2 = xcyl(nrint(jj),1) 171 & *sin(1.0*i/imax(jj)*(thmax-thmin)+thmin) 172 endif 173 xts3 = x(nrint(jj),3) + plandist 174 call elem_search(xintl, xts1, xts2, xts3, 175 & xtmp(:), etmp, 1) 176 xsinfin(jj,1+i,:) = xtmp(:) 177 elcnfin(jj,1+i) = etmp 178 enddo 179 enddo 180 181 return 182 end 183 184 185 subroutine renum_cart(x) 186c 187c---------------------------------------------------------------------- 188c This subroutine finds all nodes that are on the inlet plane and all 189c nodes that are on the recycle plane; it also blocks elements that 190c are on recycle plane; all nodes are also stored with cylindrical 191c coordinates; find all father nodes for recycle plane (i.e. for 192c theta = -theta given) 193c 194c input: 195c x (numnp,nsd) : node coordinates 196c 197c output: 198c xcyl (numnp,nsd) : node cylindrical coordinates 199c ien2D (npro, nshl) : connectivity array for recycle plane 200c assuming tethraheadral elements, i.e. 201c triangular elements on face 202c 203c 204c Elaine Bohr 205c July 2002 206c---------------------------------------------------------------------- 207c 208 use spebc 209c use pointer_data 210 include "common.h" 211 include "mpif.h" 212 include "auxmpi.h" 213c 214 dimension x(numnp,nsd), nyin(numnp), 215 & erreur(nshl), xtmp(nsd), yrtmp(numnp) 216 217 integer temp, etmp, s 218 219c 220c .... Finding nodes from inlet plane with minimal z value 221c 222 j = 0 223 zmin = x(nen1(1),3) 224 zmax = x(nen1(1),3) 225 do i=2,npin 226 if(x(nen1(i),3).lt.zmin) zmin = x(nen1(i),3) 227 if(x(nen1(i),3).gt.zmax) zmax = x(nen1(i),3) 228 enddo 229 230 do i=1,npin 231 if (x(nen1(i),3).eq.zmin) then 232 j = j + 1 233 nyin(j) = nen1(i) 234 endif 235 enddo 236 nfin = j 237 do i=1,nfin 238 yrtmp(i) = x(nyin(i),2) 239 enddo 240 241c 242c .... Ordering nyin by increasing y 243c 244 j = 0 245 allocate (nrint(nfin)) 246 247 do i=1,nfin 248 rmintemp=yrtmp(1) 249 itmp = 1 250 do k=2,nfin 251 if (rmintemp.ge.yrtmp(k)) then 252 rmintemp = yrtmp(k) 253 itmp = k 254 endif 255 enddo 256 j = j + 1 257 nrint(j)=nyin(itmp) 258 yrtmp(itmp) = 10000000 259 enddo 260 nfint = j 261 262c 263c .... y coordinate for the inlet plane 264c 265 do i=1,npin 266 xynin(i) = x(nen1(i),2) 267 enddo 268 269c 270c .... y coordinate for the virtual points on recycle plane 271c 272 273 do i=1,nfint 274 xyn(i) = x(nrint(i),2) 275 enddo 276 277c 278c .... Finding corresponding elements and local coordinates 279c on recycle plane for every arc with radius from nrint 280c 281 s = (zmax - zmin) / ds 282 allocate (xsinfin(nfint,s+1,nsd)) 283 allocate (elcnfin(nfint,s+1)) 284 allocate (imax(nfint)) 285 286 do jj = 1, nfint 287 288 xts1 = x(nrint(jj),1) + plandist 289 xts2 = x(nrint(jj),2) 290 xts3 = x(nrint(jj),3) 291 call elem_search(xintl, xts1, xts2, xts3, 292 & xtmp(:), etmp, 1) 293 xsinfin(jj,1,:) = xtmp(:) 294 elcnfin(jj,1) = etmp 295 imax(jj) = s 296 do i=1,imax(jj) 297 xts1 = x(nrint(jj),1) + plandist 298 xts2 = x(nrint(jj),2) 299 xts3 = x(nrint(jj),3) + i*ds 300 call elem_search(xintl, xts1, xts2, xts3, 301 & xtmp(:), etmp, 1) 302 xsinfin(jj,1+i,:) = xtmp(:) 303 elcnfin(jj,1+i) = etmp 304 enddo 305 enddo 306 307 return 308 end 309