!********************************************************************************
!* This program generates a tetrahedral grid in a unit cubic domain, and generates
!* grid and boundary condition files in an unstructured format (.ugrid).
!*
!*        written by Dr. Katate Masatsuka (info[at]cfdbooks.com),
!*
!* the author of useful CFD books, "I do like CFD" (http://www.cfdbooks.com).
!*
!* Tetrahedral grid is generated by dividing a hexehedron into 6 tetrahedara.
!* It is possible to divide it into 5 tetrahedra, but not implemented here.
!* You can modify this code to generate such a grid.
!*
!* This is Version 3 (2011).
!* 
!* This F90 program is written and made available for an educational purpose.
!*
!* This file may be updated in future.
!*
!* Katate Masatsuka, February 2011. http://www.cfdbooks.com
!********************************************************************************
!********************************************************************************
!* Input : nx, ny, nz    ! Number of nodes in each coordinate direction
!*
!* Output: tetgrid.ugrid ! Grid file (in .ugrid format)
!*         tetgrid.mapbc ! Boundary condition file
!*         tetgrid_tecplot.dat          !  Volume grid file for viewing by Tecplot
!*         tetgrid_boundary_tecplot.dat !Boundary grid file for viewing by Tecplot
!*
!* Note: Information on UGRID format can be found at
!*       http://www.simcenter.msstate.edu/docs/solidmesh/ugridformat.html
!*
!* Note: tetgrid.ugrid and tetgrid.mapbc are the files to be read by a flow solver
!*       for running a simulation.
!*
!* Note: BC numbers in tetgrid.mapbc must be chosen from available numbers in the
!*       flow solver.
!*
!********************************************************************************
 program tetgrid_cube

  implicit none
  integer , parameter :: dp = selected_real_kind(15) !Double precision

! Structured data
  real(dp), allocatable, dimension(:,:,:) :: x, y, z

! Mapping from (i,j,k) to the node number
  integer , allocatable, dimension(:,:,:) :: ijk2n

! Unstructured data: these data will be used to write output files.
  real(dp), allocatable, dimension(:,:)   :: xyz       ! Coordinates
  integer , allocatable, dimension(:,:)   :: tet, tria ! Tet/tria connectivities

! Number of nodes in each coordinate direction
  integer :: nx, ny, nz
! Number of tetrahedra and triangles(boundary faces)
  integer :: ntet, ntria
! Number of nodes
  integer :: nnodes
! Local numbering of a hexahedron
  integer :: hex_loc(8)
! Grid spacing in each coordinate direction
  real(dp) :: dx, dy, dz
! Local variables
  integer  :: i, j, k, os

!*******************************************************************************
! Unit cube to be gridded.
!        ______________  
!       /             /|      Boundary group numbers are assigned as below:
!      /             / |       1:  Boundary face on x = 0.0
! z=1 /____________ /  |       2:  Boundary face on x = 1.0
!     |  |         |   |       3:  Boundary face on y = 0.0
!     |  |         |   |       4:  Boundary face on y = 1.0
!     |  |         |   |       5:  Boundary face on z = 0.0
!     |  |y=1______|___|       6:  Boundary face on z = 1.0
!     |  /         |  /
!     | /          | /
!     |/___________|/
! x=y=z=0        x=1
!   
!
!*******************************************************************************
! Number of points in each direction.
  nx = 10
  ny = 10
  nz = 10

! Allocate structured arrays
  allocate( x(nx,ny,nz), y(nx,ny,nz), z(nx,ny,nz) )

! Grid spacing (uniform)

  dx = 1.0_dp / real(nx-1)
  dy = 1.0_dp / real(ny-1)
  dz = 1.0_dp / real(nz-1)

! Generate points.

  do i = 1, nx
   do j = 1, ny
    do k = 1, nz
     x(i,j,k) = dx*real(i-1)
     y(i,j,k) = dy*real(j-1)
     z(i,j,k) = dz*real(k-1)
    end do
   end do
  end do

! Construct a node array

! (1)Construct a map first.

    allocate(ijk2n(nx,ny,nz)) ! (i,j,k) to node map

     nnodes = 0

    do i = 1, nx
     do j = 1, ny
      do k = 1, nz
              nnodes = nnodes + 1
        ijk2n(i,j,k) = nnodes
      end do
     end do
    end do

! (2)Construct a node array using ijk2n(:,:,:)
!     xyz(i,1) = x at node i.
!     xyz(i,2) = y at node i.
!     xyz(i,3) = z at node i.

    allocate(xyz(nnodes,3))

    do i = 1, nx
     do j = 1, ny
      do k = 1, nz
       xyz(ijk2n(i,j,k), 1) = x(i,j,k)
       xyz(ijk2n(i,j,k), 2) = y(i,j,k)
       xyz(ijk2n(i,j,k), 3) = z(i,j,k)
      end do
     end do
    end do

! Nodes are now ordered from 1 to nnodes(=total # of nodes).
! The nodal coordinates are stored in a one-dimensional fashion:
! x_i = xyz(i,1), y_i = xyz(i,2), z_i = xyz(i,3), i=1,2,3,...,nnodes.

! So, we don't need x, y, z, in the (i,j,k) data format.
! Goodbye, x, y, z!

   deallocate(x,y,z)

! Now, construct tet element and tria boundary connectivity data:
!  Tet: 4 nodes that define the tetrahedron.
! Tria: 3 nodes that define the tria face.

! Allocate the arrays with known dimensions:
!  - Hexahedron is divided into 6 tetrahedara
!  - Quad faces are divided into 2 triangles.

  ntet  = 6 * (nx-1)*(ny-1)*(nz-1)
  ntria = 2 * ( 2*(nx-1)*(ny-1) + 2*(ny-1)*(nz-1) + 2*(nz-1)*(nx-1) )

  allocate( tet(ntet,4), tria(ntria,5) )

! Reset ntet and ntria.
   ntet = 0
  ntria = 0

! Loop over (i,j,k) and construct connectivity data.

   i_loop : do i = 1, nx-1
    j_loop : do j = 1, ny-1
     k_loop : do k = 1, nz-1

! At each step, we look at a hexahedron as shown below.
! We work with local numbering because it is convenient.
!
!  Local node numbering, {1,2,3,4,5,6,7,8}, is defined as follows.
!
!
!                         (i,j+1,k+1) 8----------------------7 (i+1,j+1,k+1)
!                                    /.                     /|
!                                   / .                    / |
!                                  /  .                   /  |
!                                 /   .                  /   |
!                                /    .                 /    |
!                               /     .                /     |
!                              /      .               /      |
!                             /       .              /       |
!                            /        .             /        |
!                           /         .            /         |
!                          /          .           /          |
!               (i,j,k+1) 5----------------------6 (i+1,j,k+1)
!                         !           .          !           !
!                         !           .          !           !
!                         !           .          !           !
!                         |           4..........|...........3
!                         | (i,j+1,k).           |          /(i+1,j+1,k)
!                         |         .            |         /
!                         |        .             |        /
!                         |       .              |       /
!                         |      .               |      /
!                         |     .                |     /
!                         |    .                 |    /
!                         |   .                  |   /
!                         |  .                   |  /
!                         | .                    | /
!                         |.                     |/
!                         1----------------------2
!                   (i,j,k)                     (i+1,j,k)
!

! Store the node nunmbers into hex_loc(1:8), so that
! we can look at the hex as defined by local numbers {1,2,3,4,5,6,7,8}.
!
!    Lower 4 vertices
     hex_loc(1) = ijk2n(i  , j  , k  )
     hex_loc(2) = ijk2n(i+1, j  , k  )
     hex_loc(3) = ijk2n(i+1, j+1, k  )
     hex_loc(4) = ijk2n(i  , j+1, k  )

!    Upper 4 vertices
     hex_loc(5) = ijk2n(i  , j  , k+1)
     hex_loc(6) = ijk2n(i+1, j  , k+1)
     hex_loc(7) = ijk2n(i+1, j+1, k+1)
     hex_loc(8) = ijk2n(i  , j+1, k+1)

! From here, we can think of the hex defned by the 
! numbers {1,2,3,4,5,6,7,8}.

!--- Hexahedron defined by the nodes 1265-4378
!    Divide the hexahedra into 6 tetrahedra as follows.

!    Tet #1
     ntet = ntet + 1
     tet(ntet,1) = hex_loc(1)
     tet(ntet,2) = hex_loc(4)
     tet(ntet,3) = hex_loc(5)
     tet(ntet,4) = hex_loc(6)

!    Tet #2
     ntet = ntet + 1
     tet(ntet,1) = hex_loc(4)
     tet(ntet,2) = hex_loc(8)
     tet(ntet,3) = hex_loc(5)
     tet(ntet,4) = hex_loc(6)


!    Tet #3
     ntet = ntet + 1
     tet(ntet,1) = hex_loc(3)
     tet(ntet,2) = hex_loc(2)
     tet(ntet,3) = hex_loc(6)
     tet(ntet,4) = hex_loc(4)
 
!    Tet #4
    ntet = ntet + 1
     tet(ntet,1) = hex_loc(3)
     tet(ntet,2) = hex_loc(6)
     tet(ntet,3) = hex_loc(7)
     tet(ntet,4) = hex_loc(4)


!    Tet #5
     ntet = ntet + 1
     tet(ntet,1) = hex_loc(1)
     tet(ntet,2) = hex_loc(2)
     tet(ntet,3) = hex_loc(4)
     tet(ntet,4) = hex_loc(6)

!    Tet #6
     ntet = ntet + 1
     tet(ntet,1) = hex_loc(4)
     tet(ntet,2) = hex_loc(7)
     tet(ntet,3) = hex_loc(8)
     tet(ntet,4) = hex_loc(6)

!--- Boundary faces if the hex is adjacent to the boundary:
!    Note: Nodes are ordered so that the boundary face points INWARD.
!
!    1. Left boundary face (x=0)
     if     (i == 1) then

              ntria = ntria + 1
      tria(ntria,1) = hex_loc(4)
      tria(ntria,2) = hex_loc(8)
      tria(ntria,3) = hex_loc(5)
      tria(ntria,4) = 1  !<------ Boundary group number

              ntria = ntria + 1
      tria(ntria,1) = hex_loc(4)
      tria(ntria,2) = hex_loc(5)
      tria(ntria,3) = hex_loc(1)
      tria(ntria,4) = 1  !<------ Boundary group number

!    2. Right boundary face (x=1)
     elseif (i == nx-1) then

              ntria = ntria + 1
      tria(ntria,1) = hex_loc(3)
      tria(ntria,2) = hex_loc(2)
      tria(ntria,3) = hex_loc(6)
      tria(ntria,4) = 2  !<------ Boundary group number

              ntria = ntria + 1
      tria(ntria,1) = hex_loc(3)
      tria(ntria,2) = hex_loc(6)
      tria(ntria,3) = hex_loc(7)
      tria(ntria,4) = 2  !<------ Boundary group number

     endif

!    3. Front boundary face (y=0)
     if     (j == 1) then

              ntria = ntria + 1
      tria(ntria,1) = hex_loc(1)
      tria(ntria,2) = hex_loc(5)
      tria(ntria,3) = hex_loc(6)
      tria(ntria,4) = 3  !<------ Boundary group number

              ntria = ntria + 1
      tria(ntria,1) = hex_loc(1)
      tria(ntria,2) = hex_loc(6)
      tria(ntria,3) = hex_loc(2)
      tria(ntria,4) = 3  !<------ Boundary group number

!    4. Back boundary face (y=1)
     elseif (j == ny-1) then

              ntria = ntria + 1
      tria(ntria,1) = hex_loc(4)
      tria(ntria,2) = hex_loc(3)
      tria(ntria,3) = hex_loc(7)
      tria(ntria,4) = 4  !<------ Boundary group number

              ntria = ntria + 1
      tria(ntria,1) = hex_loc(4)
      tria(ntria,2) = hex_loc(7)
      tria(ntria,3) = hex_loc(8)
      tria(ntria,4) = 4  !<------ Boundary group number

     endif

!    5. Bottom boundary face (z=0)
     if     (k == 1) then

              ntria = ntria + 1
      tria(ntria,1) = hex_loc(4)
      tria(ntria,2) = hex_loc(1)
      tria(ntria,3) = hex_loc(2)
      tria(ntria,4) = 5  !<------ Boundary group number

              ntria = ntria + 1
      tria(ntria,1) = hex_loc(4)
      tria(ntria,2) = hex_loc(2)
      tria(ntria,3) = hex_loc(3)
      tria(ntria,4) = 5  !<------ Boundary group number

!    6. Top boundary face (z=1)
     elseif (k == nz-1) then

              ntria = ntria + 1
      tria(ntria,1) = hex_loc(8)
      tria(ntria,2) = hex_loc(7)
      tria(ntria,3) = hex_loc(6)
      tria(ntria,4) = 6  !<------ Boundary group number

              ntria = ntria + 1
      tria(ntria,1) = hex_loc(8)
      tria(ntria,2) = hex_loc(6)
      tria(ntria,3) = hex_loc(5)
      tria(ntria,4) = 6  !<------ Boundary group number

     endif

    end do k_loop
   end do j_loop
  end do i_loop

 write(*,*) " Number of nodes      = ", nnodes
 write(*,*) " Number of tetrahedra = ", ntet
 write(*,*) " Number of tria faces = ", ntria
 
!*******************************************************************************
! Write a UGRID file
!*******************************************************************************
  open(unit=2, file="tetgrid.ugrid", status="unknown", iostat=os)

!                    #nodes, #tri_faces, #quad_faces, #tetra, #pyr, #prz, #hex
   write(2,'(7i10)') nnodes,      ntria,           0,   ntet,    0,    0,    0
   write(2, *)

! Nodal coordinates
  do i = 1, nnodes
   write(2,'(3es27.15)') (xyz(i,j), j=1,3)
  end do
 
! Triangular boundary faces
   do i = 1, ntria
    write(2,'(3i10)') tria(i,1), tria(i,2), tria(i,3)
   end do

! Face tag: Boundary group number
   do i = 1, ntria
    write(2,'(i10)') tria(i,4)
   end do

! Tetrahedra
   do i = 1, ntet
   write(2,'(4i10)') tet(i,1), tet(i,2), tet(i,3), tet(i,4)
   end do

  close(2)

!*******************************************************************************
! Write a boundary condition map file: boundary marks
! Note: Set appropriate boundary condition numbers in this file.
!*******************************************************************************
 open(unit=3, file="tetgrid.mapbc", status="unknown", iostat=os)

  write(3,*) "Boundary Group", "   BC Index"
  write(3,*) 1, 5000 !<--- These numbers must be chosen from available boundary
  write(3,*) 2, 5000 !     condition numbers in a flow solver that read this
  write(3,*) 3, 5000 !     file. Here, I just set 5000; I don't know what it is.
  write(3,*) 4, 5000 !
  write(3,*) 5, 5000 !
  write(3,*) 6, 5000 !

 close(3)

!*******************************************************************************
! Write a Tecplot volume grid file for viweing. Just for viewing.
!*******************************************************************************
 open(unit=1, file="tetgrid_tecplot.dat", status="unknown", iostat=os)

  write(1,*) 'title = "tet grid"'
  write(1,*) 'variables = "x","y","z"'
  write(1,*) 'zone  n=', nnodes,',e=', ntet,' , et=tetrahedron, f=fepoint'

  do i = 1, nnodes
    write(1,'(3es27.15)') (xyz(i,j),j=1,3)
  end do

  do i = 1, ntet
   write(1,'(4i10)') tet(i,1), tet(i,2), tet(i,3), tet(i,4)
  end do

 close(1)

!*******************************************************************************
! Write a Tecplot boundary grid file for viewing. Just for viewing.
!******************************************************************************

 open(unit=4, file="tetgrid_boundary_tecplot.dat", status="unknown", iostat=os)

  write(4,*) 'title = "tet grid boundary"'
  write(4,*) 'variables = "x","y","z"'
  write(4,*) 'zone  n=', nnodes,',e=', ntria,' , et=triangle, f=fepoint'

  do i = 1, nnodes
    write(4,'(3es27.15)') (xyz(i,j),j=1,3)
  end do

  do i = 1, ntria
   write(4,'(3i10)') tria(i,1), tria(i,2), tria(i,3)
  end do

 close(4)

!*******************************************************************************

  write(*,*) 
  write(*,*) " Grid files have been successfully generated:"
  write(*,*) " --- tetgrid.ugrid (boundary faces are ordered pointing inward.)"
  write(*,*) " --- tetgrid.mapbc"
  write(*,*) 
  write(*,*) " Tecplot files have been successfully generated:"
  write(*,*) " --- tetgrid_tecplot.dat"
  write(*,*) " --- tetgrid_boundary_tecplot.dat"

 stop

 end program tetgrid_cube