1 subroutine eqn_plane(x, iBC) 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 25 use pointer_data 26 include "common.h" 27 include "mpif.h" 28 include "auxmpi.h" 29c 30 dimension x(numnp,nsd), nrin(numnp), 31 & xA(nsd), xB(nsd), xC(nsd), xD(nsd), fourth(nsd), 32 & erreur(nshl), iBC(numnp), xtmp(nsd) 33 34 integer temp, tempb, etmp 35 36 integer, allocatable :: ien(:,:) 37 integer, allocatable :: ienb(:,:) 38 39 40c 41c .... find all nodes on inlet plane 42c 43 j = 0 44 do i=1,numnp 45 if(btest(iBC(i),11)) then 46 j = j + 1 47 nen1(j) = i 48 endif 49 enddo 50 npin = j 51 52c 53c .... find one element and its vertices on inlet plane 54c 55 do iblk=1,nelblb 56 iel=lcblkb(1,iblk) 57 nenl=lcblkb(5,iblk) 58 nenlb=lcblkb(6,iblk) 59 nshl=lcblkb(9,iblk) 60 nshlb=lcblkb(10,iblk) 61 npro=lcblkb(1,iblk+1)-iel 62 63 64 allocate (ienb(npro,nshl)) 65 66 ienb(:,:)=mienb(iblk)%p(:,:) 67 do i=1,npro 68 inum=i 69 tempb=0 70 do k=1,nshlb 71 do j=1,npin 72 if (ienb(i,k).eq.nen1(j)) then 73 tempb=tempb+1 74 endif 75 enddo 76 enddo 77 if (tempb.eq.3) then 78 xA(:) = x(ienb(inum,1),:) 79 xB(:) = x(ienb(inum,2),:) 80 xC(:) = x(ienb(inum,3),:) 81 xD(:) = x(ienb(inum,4),:) 82 deallocate (ienb) 83 goto 92 84 endif 85 enddo 86 deallocate (ienb) 87 enddo 88 92 continue 89 90c 91c .... find normal to inlet plane 92c 93 xnrml = (xB(2) - xA(2)) * (xC(3) - xA(3)) 94 & - (xB(3) - xA(3)) * (xC(2) - xA(2)) 95 96 ynrml = (xB(3) - xA(3)) * (xC(1) - xA(1)) 97 & - (xB(1) - xA(1)) * (xC(3) - xA(3)) 98 99 znrml = (xB(1) - xA(1)) * (xC(2) - xA(2)) 100 & - (xB(2) - xA(2)) * (xC(1) - xA(1)) 101 102 tmp = xnrml*xnrml + ynrml*ynrml + znrml*znrml 103 tmp = sqrt(tmp) 104 105 xnrml = xnrml/tmp 106 ynrml = ynrml/tmp 107 znrml = znrml/tmp 108 109 fourth(:) = xD(:) - xA(:) 110 scal = xnrml*fourth(1) + ynrml*fourth(2) + znrml*fourth(3) 111 112 if (scal .lt. 0.0) then 113 xnrml = -1.0*xnrml 114 ynrml = -1.0*ynrml 115 znrml = -1.0*znrml 116 endif 117 118c 119c .... find the equation of internal plane 120c equation of plane given by: 121c (xn)x + (yn)y + (zn)z = a 122c 123 aI = xnrml*xA(1)+ynrml*xA(2)+znrml*xA(3) 124c aR = xnrml*xA(1)+ynrml*xA(2)+znrml*(xA(3)+rcydist) 125 aR = aI + plandist 126 if (thetag.eq.0.0) then 127 sang = 1.0 128 else 129 angle = atan2(znrml,sqrt(xnrml*xnrml+ynrml*ynrml)) 130 sang = sin(angle) 131 endif 132 133c 134c .... blocking elements cutting the recycle plane 135c 136 itmp = lcblk(1,nelblk+1) 137 allocate (ien2D(itmp,4)) 138 nelint=0 139 do iblk=1,nelblk 140 iel=lcblk(1,iblk) 141 nenl=lcblk(5,iblk) 142 nshl=lcblk(10,iblk) 143 npro=lcblk(1,iblk+1)-iel 144 145 allocate (ien(npro,nshl)) 146 147 ien(:,:)=mien(iblk)%p(:,:) 148 149 do i=1,npro 150 inum=iel+i-1 151 152 temp=0 153 do k=1,nshl 154 erreur(k) = aR - xnrml*x(ien(i,k),1) 155 & - ynrml*x(ien(i,k),2) - znrml*x(ien(i,k),3) 156 enddo 157 158 do j=1,nshl 159 do l=j,nshl 160 if (erreur(j)*erreur(l).le. 0.0) then 161 nelint=nelint+1 162 ien2D(nelint,:) = ien(i,:) 163 goto 95 164 endif 165 enddo 166 enddo 167 95 continue 168 enddo 169 deallocate (ien) 170 171 enddo 172 173 174c 175c .... For each node of inlet plane find the corresponding element 176c and local coordinates on recycle plane 177c 178 allocate (xintl(nelint,nshl,nsd)) 179 do i = 1, nelint 180 xintl(i,:,1) = x(ien2D(i,:),1) 181 xintl(i,:,2) = x(ien2D(i,:),2) 182 xintl(i,:,3) = x(ien2D(i,:),3) 183 enddo 184 185 if (thetag .eq. 0.0) then 186 call renum_cart(x) 187 else 188 call renum_cyl(x) 189 endif 190 191 return 192 end 193