stl_preproc_des_mod.f

Go to the documentation of this file.

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!                                                                      !
!  Module name: stl_functions_des                                      !
!  Author: Rahul Garg                                 Date: 24-Oct-13  !
!                                                                      !
!  Purpose: This module containd routines for geometric interaction    !
!  required for STL files.                                             !
!                                                                      !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
      MODULE stl_preproc_des

      IMPLICIT NONE

! Use this module only to define functions and subroutines.
      CONTAINS

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!                                                                      !
!  Subroutine: DES_STL_PREPROCESSING                                   !
!  Author: Rahul Garg                                 Date: 24-Oct-13  !
!                                                                      !
!  Purpose:                                                            !
!                                                                      !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
      SUBROUTINE des_stl_preprocessing

! Flag to for STL defined geometry
      use cutcell, only: use_stl
! Number of facets from STL files, (plus DES generated)
      use stl, only: n_facets, n_facets_des
! Start/End position of different STLs
      use stl, only: stl_start, stl_end
! All STLS
      use stl, only: all_stl
! STLs read from geometry files
      use stl, only: base_stl
! STLs for user specified walls (NSW, PSW, FSW)
      use stl, only: bcwalls_stl
! STLs for impermeable surfaces
      use stl, only: imprmbl_stl
! STLs for default walls
      use stl, only: default_stl

      use stl_dbg_des
      use error_manager

      IMPLICIT NONE

! Pre-procssing for the des in order to assign facets to grid cells.
      WRITE(err_msg,"('Pre-Processing geometry for DES.')")
      CALL flush_err_msg(header=.false., footer=.false.)

! Process the STL files
      n_facets_des = merge(n_facets, 0, use_stl)
! Store the Start/End of the base STLs from geometry files
      stl_start(base_stl)=1;   stl_end(base_stl)=n_facets_des

! Process stair-step geometries
      CALL convert_bc_walls_to_stl
! Process stair-step geometries
      CALL convert_impermeable_is_to_stl
! Process default walls
      CALL convert_default_walls_to_stl

! Bin the STL to the DES grid.
      CALL bin_facets_to_dg

! Some functions for debugging.
!      CALL STL_DBG_WRITE_FACETS(BASE_STL)
!      CALL STL_DBG_WRITE_FACETS(BCWALLS_STL)
!      CALL STL_DBG_WRITE_FACETS(IMPRMBL_STL)
!      CALL STL_DBG_WRITE_FACETS(DEFAULT_STL)
!      CALL STL_DBG_WRITE_FACETS(ALL_STL)
!      CALL STL_DBG_WRITE_STL_FROM_DG(STL_TYPE=BASE_STL)

! Pre-procssing for the des in order to assign facets to grid cells.
      WRITE(err_msg,"('DES geometry pre-processing complete.')")
      CALL flush_err_msg(header=.false., footer=.false.)

      RETURN
      END SUBROUTINE des_stl_preprocessing


!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!                                                                      !
!  Subroutine: BIN_FACETS_TO_DG                                        !
!  Author: Rahul Garg                                 Date: 24-Oct-13  !
!                                                                      !
!  Purpose:                                                            !
!                                                                      !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
      Subroutine bin_facets_to_dg

      use desgrid, only: dg_ijksize2
      use desgrid, only: dg_iend2, dg_istart2
      use desgrid, only: dg_jend2, dg_jstart2
      use desgrid, only: dg_kend2, dg_kstart2
      use desgrid, only: dg_dxinv, dg_dyinv, dg_dzinv

      use stl, only: facets_at_dg

      use geometry, only: xlength, ylength, zlength, do_k
      use stl, only: n_facets_des
      use stl, only: vertex

      use stl, only: tol_stl
      use param1, only: zero, one

      use desgrid, only: dg_funijk
      use desgrid, only: iofpos, jofpos, kofpos
      use desgrid, only: dg_is_on_mype_plus1layers

      IMPLICIT NONE

! DES Grid cell index.
      INTEGER :: IJK, IJK2
! Loop counters:
      INTEGER :: I1, I2, II  ! X-axis
      INTEGER :: J1, J2, JJ  ! Y-axis
      INTEGER :: K1, K2, KK  ! Z-axis
      INTEGER :: NN          ! STLs
! Generic accumulator
      INTEGER :: COUNT_FAC

! Maximum and minimum extents of the indexed STL
      DOUBLE PRECISION:: X1,Y1,Z1
      DOUBLE PRECISION:: X2,Y2,Z2

! Allocate the data storage array.
      IF(.not.allocated(facets_at_dg)) &
         allocate(facets_at_dg(dg_ijksize2))

      facets_at_dg(:)%COUNT = 0

      DO nn = 1,n_facets_des

         x1 = minval(vertex(1:3,1,nn))
         x2 = maxval(vertex(1:3,1,nn))
         y1 = minval(vertex(1:3,2,nn))
         y2 = maxval(vertex(1:3,2,nn))
         z1 = minval(vertex(1:3,3,nn))
         z2 = maxval(vertex(1:3,3,nn))

         i1 = dg_iend2
         i2 = dg_istart2
         IF(x2>=-tol_stl .AND. x1<=xlength+tol_stl) THEN
            i1 = max(iofpos(x1)-1, dg_istart2)
            i2 = min(iofpos(x2)+1, dg_iend2)
         ENDIF

         j1 = dg_jend2
         j2 = dg_jstart2
         IF(y2>=-tol_stl .AND. y1<=ylength+tol_stl) THEN
            j1 = max(jofpos(y1)-1, dg_jstart2)
            j2 = min(jofpos(y2)+1, dg_jend2)
         ENDIF

         k1 = dg_kend2
         k2 = dg_kstart2
         IF(do_k) THEN
            IF(z2>=-tol_stl .AND. z1<=zlength+tol_stl) THEN
               k1 = max(kofpos(z1)-1, dg_kstart2)
               k2 = min(kofpos(z2)+1, dg_kend2)
            ENDIF
         ENDIF

         DO kk=k1,k2
         DO jj=j1,j2
         DO ii=i1,i2
            IF(dg_is_on_mype_plus1layers(ii,jj,kk)) THEN
               ijk = dg_funijk(ii,jj,kk)
               CALL add_facet_for_des(ii,jj,kk,ijk,nn)
            ENDIF
         ENDDO
         ENDDO
         ENDDO

      ENDDO

      RETURN
      END SUBROUTINE bin_facets_to_dg




!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!                                                                      !
!  Subroutine: ADD_FACET_FOR_DES                                       !
!  Author: Rahul Garg                                  Date: 24-Oct-13 !
!                                                                      !
!  Purpose: Add facets to DES grid cells.                              !
!                                                                      !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
      SUBROUTINE add_facet_for_des(I,J,K,IJK,N)

      use geometry, only: do_k

      use desgrid, only: dg_dxinv, dg_xstart, dg_istart1
      use desgrid, only: dg_dyinv, dg_ystart, dg_jstart1
      use desgrid, only: dg_dzinv, dg_zstart, dg_kstart1

      use discretelement, only: max_radius

      use stl, only: vertex

      use stl_functions_des, only: tri_box_overlap

      use param1, only: zero, half, one
      use error_manager

      IMPLICIT NONE

! DES grid index and facet index
      INTEGER, INTENT(IN) :: I,J,K,IJK, N

! Center of DES grid cell and half size. Note that a buffer is added to
! the half size to make the cell appear a little larger. This ensures
! that paricles near the edge 'see' STLs that are nearby but do not
! directly intersect the DES grid cell contain the particle center.
      DOUBLE PRECISION :: CENTER(3), HALFSIZE(3)
! Flag: STL intersects the DES grid cell
      LOGICAL :: OVERLAP
! DES grid cell dimensions
      DOUBLE PRECISION :: lDX, lDY, lDZ
! Buffer to ensure all particle-STL collisions are captured.
      DOUBLE PRECISION :: BUFFER
! Legacy variable - should be removed
      INTEGER ::  CURRENT_COUNT

      buffer = 1.1d0*max_radius

      ldx = one/dg_dxinv
      ldy = one/dg_dyinv
      ldz = one/dg_dzinv

      center(1) = dg_xstart + (dble(i-dg_istart1)+half)*ldx
      halfsize(1) = half*ldx + buffer

      center(2) = dg_ystart + (dble(j-dg_jstart1)+half)*ldy
      halfsize(2) = half*ldy + buffer

      IF(do_k)THEN
         center(3) = dg_zstart + (dble(k-dg_kstart1)+half)*ldz
         halfsize(3) = half*ldz + buffer
      ELSE
         center(3) = half*ldz
         halfsize(3) = half*ldz
      ENDIF

      CALL tri_box_overlap(center, halfsize, vertex(:,:,n), overlap)

      IF(overlap) CALL add_facet(ijk, n)

      RETURN
      END SUBROUTINE add_facet_for_des



!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!                                                                      !
!  Subroutine: ADD_FACET                                               !
!                                                                      !
!                                                                      !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
      SUBROUTINE add_facet(IJK, FACET_ID)

      use stl, only: vertex
      use stl, only: facets_at_dg
      use param1, only: zero

      implicit none

      INTEGER, INTENT(IN) :: IJK, facet_id

      INTEGER, ALLOCATABLE :: int_tmp(:)
      DOUBLE PRECISION, ALLOCATABLE :: real_tmp(:)

      INTEGER :: lSIZE, II
      DOUBLE PRECISION :: smallest_extent, min_temp, max_temp


      IF(facets_at_dg(ijk)%COUNT > 0) THEN

         DO ii=1, facets_at_dg(ijk)%COUNT
            IF(facet_id == facets_at_dg(ijk)%ID(ii)) RETURN
         ENDDO

         facets_at_dg(ijk)%COUNT = facets_at_dg(ijk)%COUNT+1

         lsize = size(facets_at_dg(ijk)%ID)
         IF(facets_at_dg(ijk)%COUNT +1> lsize) THEN
            allocate(int_tmp(2*lsize)); int_tmp=0
            int_tmp(1:lsize) = facets_at_dg(ijk)%ID(1:lsize)
            call move_alloc(int_tmp,facets_at_dg(ijk)%ID)

            allocate(int_tmp(2*lsize)); int_tmp=0
            int_tmp(1:lsize) = facets_at_dg(ijk)%DIR(1:lsize)
            call move_alloc(int_tmp, facets_at_dg(ijk)%DIR)

            allocate(real_tmp(2*lsize)); real_tmp=zero
            real_tmp(1:lsize) = facets_at_dg(ijk)%MIN(1:lsize)
            call move_alloc(real_tmp, facets_at_dg(ijk)%MIN)

            allocate(real_tmp(2*lsize)); real_tmp=zero
            real_tmp(1:lsize) = facets_at_dg(ijk)%MAX(1:lsize)
            call move_alloc(real_tmp, facets_at_dg(ijk)%MAX)
         ENDIF

      ELSE
         facets_at_dg(ijk)%COUNT = 1
         IF(.not.allocated(facets_at_dg(ijk)%ID)) &
            allocate(facets_at_dg(ijk)%ID(4))
         IF(.not.allocated(facets_at_dg(ijk)%DIR)) &
            allocate(facets_at_dg(ijk)%DIR(4))
         IF(.not.allocated(facets_at_dg(ijk)%MIN)) &
            allocate(facets_at_dg(ijk)%MIN(4))
         IF(.not.allocated(facets_at_dg(ijk)%MAX)) &
            allocate(facets_at_dg(ijk)%MAX(4))
      ENDIF

      facets_at_dg(ijk)%ID(facets_at_dg(ijk)%COUNT) = facet_id

      smallest_extent = huge(0.0)

      DO ii=1,3
         min_temp = minval(vertex(:,ii,facet_id))
         max_temp = maxval(vertex(:,ii,facet_id))
         IF(abs(max_temp - min_temp) < smallest_extent ) THEN
            facets_at_dg(ijk)%DIR(facets_at_dg(ijk)%COUNT) = ii
            facets_at_dg(ijk)%MIN(facets_at_dg(ijk)%COUNT) = min_temp
            facets_at_dg(ijk)%MAX(facets_at_dg(ijk)%COUNT) = max_temp
            smallest_extent = abs(max_temp - min_temp)
         ENDIF
      ENDDO

      RETURN
      END SUBROUTINE add_facet



!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!                                                                      !
!  Subroutine: CONVERT_BC_WALLS_TO_STL                                 !
!  Author: J.Musser                                   Date: 03-Nov-15  !
!                                                                      !
!  Purpose: Convert user specified walls to STLs for particle-wall     !
!  collision detection.                                                !
!                                                                      !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
      Subroutine convert_bc_walls_to_stl

      use geometry, only: zlength, do_k

      use bc, only: bc_defined, bc_type_enum, free_slip_wall, no_slip_wall, par_slip_wall
      use bc, only: bc_i_w, bc_i_e
      use bc, only: bc_j_s, bc_j_n
      use bc, only: bc_k_b, bc_k_t

      use stl, only: n_facets_des
      use stl, only: stl_start, stl_end, bcwalls_stl

      use discretelement, only: xe, yn, zt

      use param, only: dimension_bc
      USE param1, only: zero

      IMPLICIT NONE

! Loop counter.
      INTEGER :: BCV
! Extents of the BC region with respect to the fluid grid.
      DOUBLE PRECISION :: lXw, lXe, lYs, lYn, lZb, lZt

      stl_start(bcwalls_stl)=n_facets_des+1

      DO bcv=1, dimension_bc
         IF(.NOT.bc_defined(bcv)) cycle

         IF(bc_type_enum(bcv) == free_slip_wall .OR.   &
            bc_type_enum(bcv) == no_slip_wall   .OR.   &
            bc_type_enum(bcv) == par_slip_wall) THEN

            lxw = xe(bc_i_w(bcv)-1); lxe = xe(bc_i_e(bcv))
            lys = yn(bc_j_s(bcv)-1); lyn = yn(bc_j_n(bcv))
            IF(do_k) THEN
               lzb = zt(bc_k_b(bcv)-1); lzt = zt(bc_k_t(bcv))
            ELSE
               lzb = zero; lzt = zlength
            ENDIF
            CALL generate_stl_box(lxw, lxe, lys, lyn, lzb, lzt)
         ENDIF
      ENDDO
      stl_end(bcwalls_stl)=n_facets_des

      RETURN
      END SUBROUTINE convert_bc_walls_to_stl


!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!                                                                      !
!  Subroutine: CONVERT_IMPERMEABLE_IS_TO_STL                           !
!  Author: J.Musser                                   Date: 03-Nov-15  !
!                                                                      !
!  Purpose: Convert user specified impermeable surfaces to STLs for    !
!  particle-wall collision detection.                                  !
!                                                                      !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
      SUBROUTINE convert_impermeable_is_to_stl

      use geometry, only: do_k, zlength

      use is, only: is_defined, is_type
      use is, only: is_i_w, is_i_e
      use is, only: is_j_s, is_j_n
      use is, only: is_k_b, is_k_t

      use stl, only: n_facets_des
      use stl, only: stl_start, stl_end, imprmbl_stl
      use discretelement, only: xe, yn, zt

      use param, only: dimension_is
      USE param1, only: zero

      use error_manager

      IMPLICIT NONE

! Loop counter.
      INTEGER :: ISV
! Extents of the BC region with respect to the fluid grid.
      DOUBLE PRECISION :: lXw, lXe, lYs, lYn, lZb, lZt

      stl_start(imprmbl_stl)=n_facets_des+1

      DO isv=1, dimension_is
         IF(.NOT.is_defined(isv)) cycle

         IF(trim(is_type(isv)) == 'IMPERMEABLE') THEN

            lxw = xe(is_i_w(isv)-1); lxe = xe(is_i_e(isv))
            lys = yn(is_j_s(isv)-1); lyn = yn(is_j_n(isv))
            IF(do_k) THEN
               lzb = zt(is_k_b(isv)-1); lzt = zt(is_k_t(isv))
            ELSE
               lzb = zero; lzt = zlength
            ENDIF

            CALL generate_stl_box(lxw, lxe, lys, lyn, lzb, lzt)
         ELSE
            CALL init_err_msg('CONVERT_IMPERMEABLE_IS_TO_STL')
            WRITE(err_msg,1000) isv, is_type(isv)
            CALL flush_err_msg(abort=.true.)
         ENDIF
      ENDDO

      stl_end(imprmbl_stl)=n_facets_des

 1000 FORMAT("Error 1000: DES simulations do not support the ",/       &
         'specified IS TYPE:',/3x,'IS: ',i3,/3x,'IS_TYPE=',a)

      RETURN
      END SUBROUTINE convert_impermeable_is_to_stl


!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!                                                                      !
!  Subroutine: CONVERT_DEFAULT_WALLS_TO_STL                            !
!  Author: J.Musser                                   Date: 03-Nov-15  !
!                                                                      !
!  Purpose: Convert user specified walls to STLs for particle-wall     !
!  collision detection.                                                !
!                                                                      !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
      Subroutine convert_default_walls_to_stl

      USE geometry, only: do_k
      USE geometry, only: xlength, ylength, zlength
      use stl, only: vertex, norm_face
      use stl, only: n_facets_des
      use stl, only: stl_start, stl_end, default_stl

      use discretelement, only: des_periodic_walls_x
      use discretelement, only: des_periodic_walls_y
      use discretelement, only: des_periodic_walls_z

      USE param1, only: zero, one

      IMPLICIT NONE

      stl_start(default_stl)=n_facets_des+1

! West Face
      IF(.NOT.des_periodic_walls_x)THEN
         n_facets_des = n_facets_des+1
         vertex(1,:,n_facets_des) = (/zero, zero, zero/)
         vertex(2,:,n_facets_des) = (/zero, 2*ylength, zero/)
         vertex(3,:,n_facets_des) = (/zero, zero, 2*zlength/)
         norm_face(:,n_facets_des) = (/one, zero, zero/)

! East Face
         n_facets_des = n_facets_des+1
         vertex(1,:,n_facets_des) = (/xlength, zero, zero/)
         vertex(2,:,n_facets_des) = (/xlength, 2*ylength, zero/)
         vertex(3,:,n_facets_des) = (/xlength, zero, 2*zlength/)
         norm_face(:,n_facets_des) = (/-one, zero, zero/)
      ENDIF

! South Face
      IF(.NOT.des_periodic_walls_y)THEN
         n_facets_des = n_facets_des+1
         vertex(1,:,n_facets_des) = (/zero, zero, zero/)
         vertex(2,:,n_facets_des) = (/2*xlength, zero, zero/)
         vertex(3,:,n_facets_des) = (/zero, zero, 2*zlength/)
         norm_face(:,n_facets_des) = (/zero, one, zero/)

! North Face
         n_facets_des = n_facets_des+1
         vertex(1,:,n_facets_des) = (/zero, ylength, zero/)
         vertex(2,:,n_facets_des) = (/2*xlength, ylength, zero/)
         vertex(3,:,n_facets_des) = (/zero, ylength, 2*zlength/)
         norm_face(:,n_facets_des) = (/zero, -one, zero/)
      ENDIF

! Bottom Face
      IF(.NOT.des_periodic_walls_z .AND. do_k) THEN
         n_facets_des = n_facets_des+1
         vertex(1,:,n_facets_des) = (/zero, zero, zero/)
         vertex(2,:,n_facets_des) = (/2*xlength, zero, zero/)
         vertex(3,:,n_facets_des) = (/zero, 2*ylength, zero/)
         norm_face(:,n_facets_des) = (/zero, zero, one/)

! Top Face
         n_facets_des = n_facets_des+1
         vertex(1,:,n_facets_des) = (/zero, zero, zlength/)
         vertex(2,:,n_facets_des) = (/2*xlength, zero, zlength/)
         vertex(3,:,n_facets_des) = (/zero, 2*ylength, zlength/)
         norm_face(:,n_facets_des) = (/zero, zero, -one/)
      ENDIF

      stl_end(default_stl)=n_facets_des

      RETURN
      END SUBROUTINE convert_default_walls_to_stl

      END MODULE stl_preproc_des




!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!                                                                      !
!  Subroutine: GENERATE_STL_BOX                                        !
!  Author: J.Musser                                   Date: 03-Nov-15  !
!                                                                      !
!  Purpose: Given the six corners of a box, create the 12 STLs needed  !
!  to define the geometry.                                             !
!                                                                      !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
      SUBROUTINE generate_stl_box(pXw, pXe, pYs, pYn, pZb, pZt)

      use stl, only: vertex, norm_face
      use stl, only: n_facets_des

      use geometry, only: do_k

      use param1, only: zero, one

      IMPLICIT NONE

      DOUBLE PRECISION, INTENT(IN) :: pXw, pXe, pYs, pYn, pZb, pZt

! West Face
      n_facets_des = n_facets_des+1
      vertex(1,:,n_facets_des) = (/pxw, pys, pzb/)
      vertex(2,:,n_facets_des) = (/pxw, pyn, pzb/)
      vertex(3,:,n_facets_des) = (/pxw, pyn, pzt/)
      norm_face(:,n_facets_des) = (/-one, zero, zero/)

      IF(do_k)THEN
         n_facets_des = n_facets_des+1
         vertex(1,:,n_facets_des) = (/pxw, pys, pzt/)
         vertex(2,:,n_facets_des) = (/pxw, pys, pzb/)
         vertex(3,:,n_facets_des) = (/pxw, pyn, pzt/)
         norm_face(:,n_facets_des) = (/-one, zero, zero/)
      ENDIF

! East Face
      n_facets_des = n_facets_des+1
      vertex(3,:,n_facets_des) = (/pxe, pys, pzb/)
      vertex(2,:,n_facets_des) = (/pxe, pyn, pzb/)
      vertex(1,:,n_facets_des) = (/pxe, pyn, pzt/)
      norm_face(:,n_facets_des) = (/one, zero, zero/)

      IF(do_k) THEN
         n_facets_des = n_facets_des+1
         vertex(3,:,n_facets_des) = (/pxe, pys, pzt/)
         vertex(2,:,n_facets_des) = (/pxe, pys, pzb/)
         vertex(1,:,n_facets_des) = (/pxe, pyn, pzt/)
         norm_face(:,n_facets_des) = (/one, zero, zero/)
      ENDIF

! South Face
      n_facets_des = n_facets_des+1
      vertex(1,:,n_facets_des) = (/pxw, pys, pzb/)
      vertex(2,:,n_facets_des) = (/pxe, pys, pzb/)
      vertex(3,:,n_facets_des) = (/pxw, pys, pzt/)
      norm_face(:,n_facets_des) = (/zero, -one, zero/)

      IF(do_k) THEN
         n_facets_des = n_facets_des+1
         vertex(1,:,n_facets_des) = (/pxe, pys, pzt/)
         vertex(2,:,n_facets_des) = (/pxe, pys, pzb/)
         vertex(3,:,n_facets_des) = (/pxw, pys, pzb/)
         norm_face(:,n_facets_des) = (/zero, -one, zero/)
      ENDIF

! North Face
      n_facets_des = n_facets_des+1
      vertex(3,:,n_facets_des) = (/pxw, pyn, pzb/)
      vertex(2,:,n_facets_des) = (/pxe, pyn, pzb/)
      vertex(1,:,n_facets_des) = (/pxw, pyn, pzt/)
      norm_face(:,n_facets_des) = (/zero, one, zero/)

      IF(do_k) THEN
         n_facets_des = n_facets_des+1
         vertex(3,:,n_facets_des) = (/pxe, pyn, pzt/)
         vertex(2,:,n_facets_des) = (/pxe, pyn, pzb/)
         vertex(1,:,n_facets_des) = (/pxw, pyn, pzb/)
         norm_face(:,n_facets_des) = (/zero, one, zero/)
      ENDIF

! Bottom Face
      IF(do_k)THEN
         n_facets_des = n_facets_des+1
         vertex(1,:,n_facets_des) = (/pxw, pys, pzb/)
         vertex(2,:,n_facets_des) = (/pxe, pys, pzb/)
         vertex(3,:,n_facets_des) = (/pxe, pyn, pzb/)
         norm_face(:,n_facets_des) = (/zero, zero, -one/)

         n_facets_des = n_facets_des+1
         vertex(1,:,n_facets_des) = (/pxe, pyn, pzb/)
         vertex(2,:,n_facets_des) = (/pxw, pyn, pzb/)
         vertex(3,:,n_facets_des) = (/pxw, pys, pzb/)
         norm_face(:,n_facets_des) = (/zero, zero, -one/)

! Top Face
         n_facets_des = n_facets_des+1
         vertex(3,:,n_facets_des) = (/pxw, pys, pzb/)
         vertex(2,:,n_facets_des) = (/pxe, pys, pzb/)
         vertex(1,:,n_facets_des) = (/pxe, pyn, pzb/)
         norm_face(:,n_facets_des) = (/zero, zero, one/)

         n_facets_des = n_facets_des+1
         vertex(3,:,n_facets_des) = (/pxe, pyn, pzb/)
         vertex(2,:,n_facets_des) = (/pxw, pyn, pzb/)
         vertex(1,:,n_facets_des) = (/pxw, pys, pzb/)
         norm_face(:,n_facets_des) = (/zero, zero, one/)
      ENDIF

      RETURN
      END SUBROUTINE generate_stl_box