!-----------------------------------------------------------------------
! Physical boundary conditions for 3d Euler equations.
! Reflecting walls at all sides.
! Interface:
! mx,my := shape of grid function
!
! u(,) := grid function
!
! lb(2) := lower bound for grid
! ub(2) := upper bound for grid
! lbbnd(2) := lower bound for boundary region
! ubnd(2) := upper bound for boundary region
! shapebnd(2) := shape of boundary region
! xc(2) := lower left corner of grid
! dx(2) := grid spacing
! dir := at which side of the grid is the boundary?
! bnd(,2,2) := lower left and upper right corner of global grid and
! of mb-1 internal boundary regions
!
! A general implementation of boundary conditions is required, because
! - a boundary region can have a width of 1 or 2 ghost cells.
! - corners of a grid might be part of two boundaries and it is NOT
! guaranteed that ALL cells of a particular corner region are
! contained in both boundaries.
!-----------------------------------------------------------------------
subroutine physbd(u,mx,my,mz,lb,ub,lbbnd,ubbnd,shapebnd,
& xc,dx,dir,bnd,mb,time,meqn)
implicit none
include "cuser.i"
integer meqn, mx, my, mz, 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(10), isy(10), isz(10)
c
data isx /1,-1,1,1,1,1,1,1,1,1/
data isy /1,1,-1,1,1,1,1,1,1,1/
data isz /1,1,1,-1,1,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 --- periodic
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
100 continue
return
! Right Side --- periodic
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
200 continue
return
! Bottom Side --- periodic
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
300 continue
return
! Top Side --- periodic
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
400 continue
return
! Front side --- periodic
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
500 continue
return
! Back Side --- periodic
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
600 continue
return
end