vtf-logo

fsi/sfc-amroc/VibratingPanel/src/physbd3.f

!-----------------------------------------------------------------------
! Physical boundary conditions for 3d Euler equations.
! Reflecting walls at all sides.
! Physical boundary conditions
! Interface:
!   mx,my,mz  := shape of grid 
!
!   u(,) := grid function
!
!   lb(3) := lower bound for grid
!   ub(3) := upper bound for grid
!   lbbnd(3) := lower bound for boundary region
!   ubnd(3) := upper bound for boundary region
!   shapebnd(3) := shape of boundary region 
!   xc(3) := lower left corner of grid
!   dx(3) := grid spacing
!   dir := at which side of the grid is the boundary?
!   bnd(,2,3) := lower left and upper right corner of global grid and 
!      of mb-1 internal boundary regions 
!
! Copyright (C) 2002 Ralf Deiterding
! Brandenburgische Universitaet Cottbus
!
! Copyright (C) 2003-2007 California Institute of Technology
! Ralf Deiterding, ralf@amroc.net
!
!-----------------------------------------------------------------------

      subroutine physbd(u,mx,my,mz,lb,ub,lbbnd,ubbnd,shapebnd,
     &     xc,dx,dir,bnd,mb,time,meqn)

      implicit none
      include  "cuser.i"

      integer   mx, my, mz, meqn, mb, dir
      integer   lb(3), ub(3), lbbnd(3), ubbnd(3), shapebnd(3)
      double precision u(meqn,mx,my,mz), xc(3), dx(3), bnd(mb,2,3), time

!      Local variables
      integer   i, j, k, imin, imax, jmin, jmax, kmin, kmax, m 
      integer   stride, getindx, isym, jsym, ksym
      integer   isx(5), isy(5), isz(5)
c
      data isx /1,-1,1,1,1/
      data isy /1,1,-1,1,1/
      data isz /1,1,1,-1,1/
      
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!      See definition of member-function extents() in BBox.h 
!      for calculation of stride
         
      stride = (ub(1) - lb(1))/(mx-1)
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!     Find working domain
         
      imin = getindx(max(lbbnd(1), lb(1)), lb(1), stride)
      imax = getindx(min(ubbnd(1), ub(1)), lb(1), stride)
      
      jmin = getindx(max(lbbnd(2), lb(2)), lb(2), stride)
      jmax = getindx(min(ubbnd(2), ub(2)), lb(2), stride)
      
      kmin = getindx(max(lbbnd(3), lb(3)), lb(3), stride)
      kmax = getindx(min(ubbnd(3), ub(3)), lb(3), stride)
      
      if(imax .gt. mx .or. jmax .gt. my .or. kmax .gt. mz .or. 
     &     imin .lt. 1 .or. jmin .lt. 1 .or. kmin .lt. 1) then
         write(0,*)'INDEX ERROR in physbd'
      end if

      go to (100,200,300,400,500,600) dir+1

!        Left Side --- Inflow/Outflow
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 100  continue
      do 110 i = imax, imin, -1
         if (xc(1)+(i-0.5d0)*dx(1).lt.bnd(1,1,1)) then
            do 120 k = kmin, kmax
               do 120 j = jmin, jmax
                  do 120 m = 1, meqn
                     u(m,i,j,k) = u(m,i+1,j,k)
 120        continue
         endif
 110  continue
      return

!        Right Side --- Inflow/Outflow or Reflection
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 200  continue
      do 210 i = imin, imax
         if (xc(1)+(i-0.5d0)*dx(1).gt.bnd(1,2,1)) then
            if (ibnd.eq.1) then 
               do 220 k = kmin, kmax
                  do 220 j = jmin, jmax
                     do 220 m = 1, meqn
                        u(m,i,j,k) = u(m,i-1,j,k)
 220           continue
            else
               isym = 2*imin-1-i
               do 230 k = kmin, kmax
                  do 230 j = jmin, jmax
                     do 230 m = 1, meqn
                        u(m,i,j,k) = u(m,isym,j,k)*isx(m)
 230           continue
            endif
         endif
 210  continue
      return
         
!        Bottom Side --- Reflection
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 300  continue
      do 310 j = jmax, jmin, -1
         if (xc(2)+(j-0.5d0)*dx(2).lt.bnd(1,1,2)) then
            jsym = 2*jmax+1-j
            do 320 k = kmin, kmax
               do 320 i = imin, imax
                  do 320 m = 1, meqn
                     u(m,i,j,k) = u(m,i,jsym,k)*isy(m)
 320        continue
         endif
 310  continue
      return

!        Top Side --- Reflection
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 400  continue
      do 410 j = jmin, jmax
         if (xc(2)+(j-0.5d0)*dx(2).gt.bnd(1,2,2)) then
            jsym = 2*jmin-1-j
            do 420 k = kmin, kmax
               do 420 i = imin, imax
                  do 420 m = 1, meqn
                     u(m,i,j,k) = u(m,i,jsym,k)*isy(m)
 420        continue
         endif
 410  continue
      return

!        Front side --- Reflection
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 500  continue
      do 510 k = kmax, kmin, -1
         if (xc(3)+(k-0.5d0)*dx(3).lt.bnd(1,1,3)) then
            ksym = 2*kmax+1-k
            do 520 j = jmin, jmax
               do 520 i = imin, imax
                  do 520 m = 1, meqn
                     u(m,i,j,k) = u(m,i,j,ksym)*isz(m)
 520        continue
         endif
 510  continue
      return

!        Back Side --- Reflection
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 600  continue
      do 610 k = kmin, kmax
         if (xc(3)+(k-0.5d0)*dx(3).gt.bnd(1,2,3)) then
            ksym = 2*kmin-1-k
            do 620 j = jmin, jmax
               do 620 i = imin, imax
                  do 620 m = 1, meqn
                     u(m,i,j,k) = u(m,i,j,ksym)*isz(m)
 620        continue
         endif
 610  continue
      return

      end