!******************************************************************************** !* This program generates a prismatic 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). !* !* This is Version 3 (2011). !* !* This F90 program is written and made available for an educational purpose. !* It may be useful for learning how to create unstructured grid data. !* !* 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: przgrid.ugrid ! Grid file (in .ugrid format) !* przgrid.mapbc ! Boundary condition file !* przgrid_tecplot.dat ! Volume grid file for viewing by Tecplot !* przgrid_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: przgrid.ugrid and przgrid.mapbc are the files to be read by a flow solver !* for running a simulation. !* !* Note: BC numbers in przgrid.mapbc must be chosen from available numbers in the !* flow solver. !* !******************************************************************************** program przgrid_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(:,:) :: prz, tria, quad ! Prz/tria/quad connectivities ! Number of nodes in each coordinate direction integer :: nx, ny, nz ! Number of prisms, triangles and quadrilaterals(boundary faces) integer :: nprz, ntria, nquad ! 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 prz element and quad boundary connectivity data: ! Prz: 6 nodes that define the hexahedron. ! Tria: 3 nodes that define the tria face. ! Quad: 4 nodes that define the quad face. ! Allocate the arrays with known dimensions. nprz = 2 * (nx-1)*(ny-1)*(nz-1) ntria = 2 * 2*(nx-1)*(ny-1) !Triangles on top and bottom nquad = 2*(ny-1)*(nz-1) + 2*(nz-1)*(nx-1)!Quadrilaterals on sides allocate( prz(nprz,8), tria(ntria,4), quad(nquad,5) ) ! Reset nhex and nquad. nprz = 0 ntria = 0 nquad = 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) ! ! ! Break it into 2 prisms: ! ! Prism #1 Prism #2 ! ! (i,j+1,k+1) 8 (i,j+1,k+1) 8----------------------7 (i+1,j+1,k+1) ! /.\ |\ /| ! / . \ | \ / | ! / . \ | \ / | ! / . \ | \ / | ! / . \ | \ / | ! / . \ | \ / | ! / . \ | \ / | ! / . \ | \ / | ! / . \ | \ / | ! / . \ | \ / | ! / . \ | \/ | !(i,j,k+1) 5----------------------6 (i+1,j,k+1) | 6 (i+1,j,k+1) ! ! . ! | ! ! ! ! . ! | ! ! ! ! . ! | ! ! ! | 4 | 4..........|...........3 ! | (i,j+1,k). . | (i,j+1,k) \ | /(i+1,j+1,k) ! | . . | \ | / ! | . . | \ | / ! | . . | \ | / ! | . . | \ | / ! | . . | \ | / ! | . . | \ | / ! | . . | \ | / ! | . . | \ | / ! | . .| \| / ! |. | |/ ! 1----------------------2 2 ! (i,j,k) (i+1,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 2 prisms as indicated by the lines of 'p' ! int he above figure: ! ! Prism #1 nprz = nprz + 1 prz(nprz,1) = hex_loc(1) prz(nprz,2) = hex_loc(2) prz(nprz,3) = hex_loc(4) prz(nprz,4) = hex_loc(5) prz(nprz,5) = hex_loc(6) prz(nprz,6) = hex_loc(8) ! Prism #2 nprz = nprz + 1 prz(nprz,1) = hex_loc(4) prz(nprz,2) = hex_loc(2) prz(nprz,3) = hex_loc(3) prz(nprz,4) = hex_loc(8) prz(nprz,5) = hex_loc(6) prz(nprz,6) = hex_loc(7) !--- 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) : Quad face if (i == 1) then nquad = nquad + 1 quad(nquad,1) = hex_loc(1) quad(nquad,2) = hex_loc(4) quad(nquad,3) = hex_loc(8) quad(nquad,4) = hex_loc(5) quad(nquad,5) = 1 !<------ Boundary group number ! 2. Right boundary face (x=1) : Quad face elseif (i == nx-1) then nquad = nquad + 1 quad(nquad,1) = hex_loc(3) quad(nquad,2) = hex_loc(2) quad(nquad,3) = hex_loc(6) quad(nquad,4) = hex_loc(7) quad(nquad,5) = 2 !<------ Boundary group number endif ! 3. Front boundary face (y=0) : Quad face if (j == 1) then nquad = nquad + 1 quad(nquad,1) = hex_loc(2) quad(nquad,2) = hex_loc(1) quad(nquad,3) = hex_loc(5) quad(nquad,4) = hex_loc(6) quad(nquad,5) = 3 !<------ Boundary group number ! 4. Back boundary face (y=1) : Quad face elseif (j == ny-1) then nquad = nquad + 1 quad(nquad,1) = hex_loc(4) quad(nquad,2) = hex_loc(3) quad(nquad,3) = hex_loc(7) quad(nquad,4) = hex_loc(8) quad(nquad,5) = 4 !<------ Boundary group number endif ! 5. Bottom boundary face (z=0) : Tria face 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) : Tria face 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 prisms = ", nprz write(*,*) " Number of tria faces = ", ntria write(*,*) " Number of quad faces = ", nquad !******************************************************************************* ! Write a UGRID file !******************************************************************************* open(unit=2, file="przgrid.ugrid", status="unknown", iostat=os) ! #nodes, #tri_faces, #quad_faces, #tetra, #pyr, #prz, #hex write(2,'(7i10)') nnodes, ntria, nquad, 0, 0, nprz, 0 write(2, *) ! Nodal coordinates do i = 1, nnodes write(2,'(3es27.15)') (xyz(i,j), j=1,3) end do ! Tria boundary faces do i = 1, ntria write(2,'(3i10)') tria(i,1), tria(i,2), tria(i,3) end do ! Quad boundary faces do i = 1, nquad write(2,'(4i10)') quad(i,1), quad(i,2), quad(i,3), quad(i,4) end do ! Face tag: Boundary group number do i = 1, ntria write(2,'(i10)') tria(i,4) end do do i = 1, nquad write(2,'(i10)') quad(i,5) end do ! Prisms do i = 1, nprz write(2,'(6i10)') prz(i,1), prz(i,2), prz(i,3), prz(i,4), prz(i,5), prz(i,6) 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="przgrid.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="przgrid_tecplot.dat", status="unknown", iostat=os) write(1,*) 'title = "przmatic grid"' write(1,*) 'variables = "x","y","z"' write(1,*) 'zone n=', nnodes,',e=', nprz,' , et=brick, f=fepoint' do i = 1, nnodes write(1,'(3es27.15)') (xyz(i,j),j=1,3) end do do i = 1, nprz write(1,'(8i10)') prz(i,1), prz(i,2), prz(i,3), prz(i,3), & prz(i,4), prz(i,5), prz(i,6), prz(i,6) end do close(1) !******************************************************************************* ! Write a Tecplot boundary grid file for viewing. Just for viewing. !****************************************************************************** open(unit=4, file="przgrid_boundary_tecplot.dat", status="unknown", iostat=os) write(4,*) 'title = "przmatic grid boundary"' write(4,*) 'variables = "x","y","z"' write(4,*) 'zone n=', nnodes,',e=', ntria+nquad,', et=quadrilateral, f=fepoint' do i = 1, nnodes write(4,'(3es27.15)') (xyz(i,j),j=1,3) end do do i = 1, ntria write(4,'(4i10)') tria(i,1), tria(i,2), tria(i,3), tria(i,1) end do do i = 1, nquad write(4,'(4i10)') quad(i,1), quad(i,2), quad(i,3), quad(i,4) end do close(4) !******************************************************************************* write(*,*) write(*,*) " Grid files have been successfully generated:" write(*,*) " --- przgrid.ugrid (boundary faces are ordered pointing inward.)" write(*,*) " --- przgrid.mapbc" write(*,*) write(*,*) " Tecplot files have been successfully generated:" write(*,*) " --- przgrid_tecplot.dat" write(*,*) " --- przgrid_boundary_tecplot.dat" stop end program przgrid_cube