!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!                                                                              |
!              ||===\\                                                         | 
!              ||    \\                                                        |
!              ||     ||   //==\\   ||  ||   //==||  ||/==\\                   |
!              ||     ||  ||    ||  ||  ||  ||   ||  ||    ||                  |
!              ||    //   ||    ||  ||  ||  ||   ||  ||                        |
!              ||===//     \\==//    \\==\\  \\==\\  ||                        |
!                                                                              |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!                                                                              |
!              DEFINE_BC_3Dpunch    Apr. 2007                                  |
!                                                                              |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|

subroutine define_bc_3Dpunch (nnode,kfix,kfixt,x,y,z,u,v,w,temp,vo) 

!------------------------------------------------------------------------------|
!(((((((((((((((( Purpose of the routine  ))))))))))))))))))))))))))))))))))))))
!------------------------------------------------------------------------------|
! this routine assigns the boundary condition in the case of a square, 
! rectangular, or circular punch

!------------------------------------------------------------------------------|
!((((((((((((((((  declaration of the subroutine arguments  ))))))))))))))))))))
!------------------------------------------------------------------------------|

use definitions

implicit none

integer nnode
integer kfix(nnode*3)
integer kfixt(nnode)
double precision x(nnode),y(nnode),z(nnode)
double precision u(nnode),v(nnode),w(nnode)
double precision temp(nnode)
type (void) vo

!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine internal variables )))))))))))))
!------------------------------------------------------------------------------|

integer i
double precision eps,square_size,rect_x_size,rect_y_size,circle_radius
double precision imposed_velocity
logical circle_punch
logical square_punch
logical rect_punch
logical rough

!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|

rough=.true.

square_punch = .true.
square_size = 0.15d0

rect_punch=.false.
rect_x_size =  0.06d0
rect_y_size =  0.08d0

circle_punch=.false.
circle_radius = 0.05d0

imposed_velocity=1.05d0

eps=1.d-10

do i=1,nnode
   if (x(i).lt.eps) then
      kfix((i-1)*3+1)=1 ; u(i)=0.d0
      kfix((i-1)*3+2)=1 ; v(i)=0.d0
      kfix((i-1)*3+3)=1 ; w(i)=0.d0
   endif
   if (x(i).gt.1.d0-eps) then
      kfix((i-1)*3+1)=1 ; u(i)=0.d0
      kfix((i-1)*3+2)=1 ; v(i)=0.d0
      kfix((i-1)*3+3)=1 ; w(i)=0.d0
   endif
   if (y(i).lt.eps) then
      kfix((i-1)*3+1)=1 ; u(i)=0.d0
      kfix((i-1)*3+2)=1 ; v(i)=0.d0
      kfix((i-1)*3+3)=1 ; w(i)=0.d0
   endif
   if (y(i).gt.1.d0-eps) then
      kfix((i-1)*3+1)=1 ; u(i)=0.d0
      kfix((i-1)*3+2)=1 ; v(i)=0.d0
      kfix((i-1)*3+3)=1 ; w(i)=0.d0
   endif
   if (z(i).lt.eps) then
      if (square_punch) then
         if ((x(i)-(.5d0-.5d0*square_size))*(x(i)-(.5d0+.5d0*square_size)).le.0.d0 .and. &
             (y(i)-(.5d0-.5d0*square_size))*(y(i)-(.5d0+.5d0*square_size)).le.0.d0 ) then 
            if (rough) then
               kfix((i-1)*3+1)=1 ; u(i)=0.d0
               kfix((i-1)*3+2)=1 ; v(i)=0.d0
            end if
            kfix((i-1)*3+3)=1 ; w(i)=imposed_velocity
         endif
      endif
      if (rect_punch) then
         if ((x(i)-(.5d0-.5d0*rect_x_size))*(x(i)-(.5d0+.5d0*rect_x_size)).le.0.d0 .and. &
             (y(i)-(.5d0-.5d0*rect_y_size))*(y(i)-(.5d0+.5d0*rect_y_size)).le.0.d0 ) then 
            if (rough) then
               kfix((i-1)*3+1)=1 ; u(i)=0.d0
               kfix((i-1)*3+2)=1 ; v(i)=0.d0
            end if
            kfix((i-1)*3+3)=1 ; w(i)=imposed_velocity
         endif
      endif
      if (circle_punch) then
         if (sqrt((x(i)-.5d0)**2+(y(i)-.5d0)**2).le.circle_radius) then
            if (rough) then
               kfix((i-1)*3+1)=1 ; u(i)=0.d0
               kfix((i-1)*3+2)=1 ; v(i)=0.d0
            end if
            kfix((i-1)*3+3)=1 ; w(i)=imposed_velocity
         endif
      endif
      kfixt(i)=1
      temp(i)=1.d0
   endif
   if (z(i).gt.1.d0-eps) then
      kfix((i-1)*3+1)=1 ; u(i)=0.d0
      kfix((i-1)*3+2)=1 ; v(i)=0.d0
      kfix((i-1)*3+3)=1 ; w(i)=0.d0
      kfixt(i)=1
      temp(i)=0.d0
   endif
   if (.not.vo%influid(i)) then
      kfixt(i)=1
      temp(i)=0.d0
   endif
enddo

return
end

!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------