xref: /phasta/phSolver/common/eqn_plane.f (revision 595995161822a203c8467e0e4a253d7bd7d6df32)

      subroutine eqn_plane(x, iBC)
c
c----------------------------------------------------------------------
c This subroutine finds all nodes that are on the inlet plane and all
c nodes that are on the recycle plane; it also blocks elements that
c are on recycle plane; all nodes are also stored with cylindrical
c coordinates; find all father nodes for recycle plane (i.e. for
c theta = -theta given)
c
c input:
c  x      (numnp,nsd)           : node coordinates
c
c output:
c  xcyl   (numnp,nsd)           : node cylindrical coordinates
c  ien2D  (npro, nshl)		: connectivity array for recycle plane
c				  assuming tethraheadral elements, i.e.
c				  triangular elements on face
c
c
c  Elaine Bohr
c  June 2002
c----------------------------------------------------------------------
c
       use spebc
       use pointer_data
       include "common.h"
       include "mpif.h"
       include "auxmpi.h"
c
        dimension x(numnp,nsd),       nrin(numnp),
     &  	  xA(nsd),  xB(nsd),  xC(nsd), xD(nsd), fourth(nsd),
     &		  erreur(nshl), iBC(numnp), xtmp(nsd)

        integer  temp, tempb, etmp

        integer, allocatable :: ien(:,:)
        integer, allocatable :: ienb(:,:)


c
c .... find all nodes on inlet plane
c
        j = 0
        do i=1,numnp
	  if(btest(iBC(i),11)) then
	    j = j + 1
	    nen1(j) = i
	  endif
	enddo
	npin = j

c
c .... find one element and its vertices on inlet plane
c
        do iblk=1,nelblb
          iel=lcblkb(1,iblk)
	  nenl=lcblkb(5,iblk)
          nenlb=lcblkb(6,iblk)
	  nshl=lcblkb(9,iblk)
          nshlb=lcblkb(10,iblk)
	  npro=lcblkb(1,iblk+1)-iel


	  allocate (ienb(npro,nshl))

	  ienb(:,:)=mienb(iblk)%p(:,:)
	  do i=1,npro
	    inum=i
	    tempb=0
	    do k=1,nshlb
	      do j=1,npin
	        if (ienb(i,k).eq.nen1(j)) then
	          tempb=tempb+1
	        endif
	      enddo
	    enddo
	    if (tempb.eq.3) then
	      xA(:) = x(ienb(inum,1),:)
	      xB(:) = x(ienb(inum,2),:)
	      xC(:) = x(ienb(inum,3),:)
	      xD(:) = x(ienb(inum,4),:)
	      deallocate (ienb)
	      goto 92
	    endif
	  enddo
	  deallocate (ienb)
        enddo
 92     continue

c
c .... find normal to inlet plane
c
	xnrml = (xB(2) - xA(2)) * (xC(3) - xA(3))
     &        - (xB(3) - xA(3)) * (xC(2) - xA(2))

     	ynrml = (xB(3) - xA(3)) * (xC(1) - xA(1))
     &        - (xB(1) - xA(1)) * (xC(3) - xA(3))

     	znrml = (xB(1) - xA(1)) * (xC(2) - xA(2))
     &        - (xB(2) - xA(2)) * (xC(1) - xA(1))

     	tmp = xnrml*xnrml + ynrml*ynrml + znrml*znrml
	tmp = sqrt(tmp)

	xnrml = xnrml/tmp
     	ynrml = ynrml/tmp
     	znrml = znrml/tmp

	fourth(:) = xD(:) - xA(:)
	scal = xnrml*fourth(1) + ynrml*fourth(2) + znrml*fourth(3)

	if (scal .lt. 0.0) then
	  xnrml = -1.0*xnrml
	  ynrml = -1.0*ynrml
	  znrml = -1.0*znrml
	endif

c
c .... find the equation of internal plane
c      equation of plane given by:
c      (xn)x + (yn)y + (zn)z = a
c
        aI = xnrml*xA(1)+ynrml*xA(2)+znrml*xA(3)
c	aR = xnrml*xA(1)+ynrml*xA(2)+znrml*(xA(3)+rcydist)
	aR = aI + plandist
	if (thetag.eq.0.0) then
	  sang = 1.0
	else
	  angle = atan2(znrml,sqrt(xnrml*xnrml+ynrml*ynrml))
	  sang = sin(angle)
	endif

c
c .... blocking elements cutting the recycle plane
c
	itmp = lcblk(1,nelblk+1)
	allocate (ien2D(itmp,4))
       	nelint=0
       	do iblk=1,nelblk
          iel=lcblk(1,iblk)
	  nenl=lcblk(5,iblk)
	  nshl=lcblk(10,iblk)
	  npro=lcblk(1,iblk+1)-iel

	  allocate (ien(npro,nshl))

	  ien(:,:)=mien(iblk)%p(:,:)

	  do i=1,npro
	    inum=iel+i-1

	    temp=0
	    do k=1,nshl
	      erreur(k) = aR - xnrml*x(ien(i,k),1)
     &		- ynrml*x(ien(i,k),2) - znrml*x(ien(i,k),3)
     	    enddo

	    do j=1,nshl
	      do l=j,nshl
	        if (erreur(j)*erreur(l).le. 0.0) then
 	          nelint=nelint+1
		  ien2D(nelint,:) = ien(i,:)
		  goto 95
	        endif
	      enddo
	    enddo
 95	    continue
	  enddo
	  deallocate (ien)

        enddo


c
c .... For each node of inlet plane find the corresponding element
c      and local coordinates on recycle plane
c
	allocate (xintl(nelint,nshl,nsd))
	do i = 1, nelint
              xintl(i,:,1) = x(ien2D(i,:),1)
	      xintl(i,:,2) = x(ien2D(i,:),2)
	      xintl(i,:,3) = x(ien2D(i,:),3)
        enddo

	if (thetag .eq. 0.0) then
	  call renum_cart(x)
	else
	  call renum_cyl(x)
	endif

        return
        end