!******************************************************************************** ! This program generates a grid for a 3D rectangular domain with a bump. ! ! 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 4 (2019). ! ! -- 12/09/19 ! Added an option to generate a triangular bump. ! Write boundary parts separetely in the tecplot boundary file. ! ! ! NOTE: The file name is not updated to preserve existing links to the code. ! ! ! This F90 program is written and made available for an educational purpose. ! ! This file may be updated in future. ! ! ! Katate Masatsuka, December 2019. http://www.cfdbooks.com !******************************************************************************** !------------------------------------------------------------------------------- !------------------------------------------------------------------------------- !------------------------------------------------------------------------------- !------------------------------------------------------------------------------- !------------------------------------------------------------------------------- !------------------------------------------------------------------------------- ! ! Input parameter module ! !------------------------------------------------------------------------------- ! ! Sample input file: 'input.nml' to generate a grid. ! ------------------------------------------------------ ! &input_parameters ! bump_type = 1 !1=circular arc, 2=smooth sine, 3=triangular ! bump_on_top = F !T=Bump place on top, F=on the bottom ! prism_on_top = F !T=Prism/Tet, F=Tet/Prism (<-Top/Bottom) ! separate_bump = T !T = Sym|Bump|Sym, F=Bottom ! xmin = -2.0 !Left end of the domain ! xmax = 3.0 !Right end of the domain ! ymin = 0.0 !y-plane facing me in xz-plane. ! ymax = 1.0 !y-plane facing away from me in xz-plane. ! zmin = 0.0 !Bottom ! zmax = 1.0 !Top ! xle = 0.0 !Leading edge location of the bump ! xte = 1.0 !Trailing edge location of the bump ! tria_height = 1.0 !Height of triangular bump ! nx1 = 32 !# of cells along x-drection in the left region. ! nx2 = 24 !# of cells along a bump. ! nx3 = 32 !# of cells along x-drection in the right region. ! ny = 3 !# of cells along y-drection ! nz1 = 16 !# of cells in the top part in z-direction ! nz2 = 18 !# of cells in the bottom part in z-direction ! circular_bump_radius = 2.5 !Circular bump radius. The larger the thinner. ! sine_bump_zmax = 0.05 !Sine_bump_zmax ! stretching_factor_z = 2.0 !Stretching in z: +ve cluster to bottom; -ve to top. ! generate_ugrid_file = T !T = Write a .ugrid file.(required by the coarsening program) ! ugrid_file_unformatted = T !T = unformatted .ugrid/.ufmt, F = formatted .ugrid/.p3d ! generate_tec_file_b = T !T = Write a boundary grid (Tecplot) ! generate_tec_file_v = F !T = Write a volume grid (Tecplot) ! generate_su2grid_file = T !T = Write .su2 grid file. ! generate_vtk_file = T !T = Write .vtk grid file. ! / ! ------------------------------------------------------ ! ! Note: No need to specify all namelist variables. ! In the above, those not shown are given their default values ! as defined below. ! !------------------------------------------------------------------------------- !------------------------------------------------------------------------------- !------------------------------------------------------------------------------- module input_parameter_module implicit none integer , parameter :: dp = selected_real_kind(P=15) public !---------------------------- ! Default input values !---------------------------- !---------------------------- ! Type of bump: ! 1 = circular arc ! 2 = smooth sine (Generarization of sine^4 at https://turbmodels.larc.nasa.gov/bump.html) ! 3 = triangular bump ! ! Note: If xle=0.3, xte=1.2, and sine_bump_zmax = 0.05, then the bump is the same ! as the one described in https://turbmodels.larc.nasa.gov/bump.html. integer :: bump_type = 1 !---------------------------- ! Place the bump on the top part (=.true.) ! or on the bottom (=.false.) ! logical :: bump_on_top = .false. !---------------------------- ! Prism layers on the top part (=.true.) ! or on the bottom (=.false.) for the mixed grid case: ! nz1 /= 0 and nz2 /= 0. ! logical :: prism_on_top = .false. !---------------------------- ! To treat the bump as a separate boundary. ! T = Top/Bottom is split into symmetry, bump, symmetry ! F = Top/Bottom is treated as a single boundary ! ! Note: Top if bump_on_top = T ! Bottom if bump_on_top = F ! logical :: separate_bump = .false. !---------------------------- ! Domain real(dp) :: xmin = -2.0_dp !Left end of the domain real(dp) :: xmax = 3.0_dp !Right end of the domain real(dp) :: ymin = 0.0_dp !y-plane facing me in xz-plane. real(dp) :: ymax = 1.0_dp !y-plane facing away from me in xz-plane. real(dp) :: zmin = 0.0_dp !Bottom real(dp) :: zmax = 1.0_dp !Top !---------------------------- ! Leading edge location of the bump, ! and trailing edge location of the bump. real(dp) :: xle = 0.0_dp real(dp) :: xte = 1.0_dp !---------------------------- ! Height of triangular bump ! This is used only when bump_type = 3. real(dp) :: tria_height = 0.1_dp !---------------------------- ! ! integer :: nx1 = 32 !# of cells along x-drection in the left region. integer :: nx2 = 24 !# of cells along a bump. integer :: nx3 = 32 !# of cells along x-drection in the right region. integer :: ny = 3 !# of cells along y-drection integer :: nz1 = 16 !# of cells in the bottom part in z-direction integer :: nz2 = 18 !# of cells in the top part in z-direction ! Example: To generate a pure tetrahedral grid (prism_on_top=.true. ) ! or a pure prismatic grid (prism_on_top=.false.). ! nz1 = 34 ! nz2 = 0 ! Example: To generate a pure prismatic grid (prism_on_top=.true. ) ! or a pure tetrahedral grid (prism_on_top=.false.). ! nz1 = 0 ! nz2 = 34 !---------------------------- ! Circular bump radius. The larger the thinner. real(dp) :: circular_bump_radius = 2.5_dp !---------------------------- ! Maximum thickness of a sine^4 bump, ! which is set to be 0.05 in https://turbmodels.larc.nasa.gov/bump.html real(dp) :: sine_bump_zmax = 0.05_dp !---------------------------- ! Exponential stretching in z-direction ! stretching_factor_z > 1.0 -> Cluster towards the bottom. ! stretching_factor_z < -1.0 -> Cluster towards the top. real(dp) :: stretching_factor_z = 2.0_dp !---------------------------- ! generate_ugrid_file = T to write .ugrid file (required by the coarsening program). ! F not to write. ! >>> [=T Required by the coarsening program] logical :: generate_ugrid_file = .true. !---------------------------- ! ugrid_file_unformatted = T: unformatted, F: formatted logical :: ugrid_file_unformatted = .true. !---------------------------- ! generate_tec_file_b = T: Write a Tecplot file for boundary. logical :: generate_tec_file_b = .false. !---------------------------- ! generate_tec_file_v = T: Write a Tecplot file for the entire volume grid. logical :: generate_tec_file_v = .false. !---------------------------- ! generate_su2grid_file = T to write .su2 file ! F not to write. logical :: generate_su2grid_file = .false. !---------------------------- ! generate_vtk_file = T to write .vtk file ! F not to write. logical :: generate_vtk_file = .false. !---------------------------- ! End of Default input values !---------------------------- !------------------------------------------------------------------------ ! List of all input variables: !------------------------------------------------------------------------ namelist / input_parameters / & bump_type , & bump_on_top, & prism_on_top, & separate_bump, & xmin, & xmax, & ymin, & ymax, & zmin, & zmax, & xle, & xte, & tria_height, & nx1, & nx2, & nx3, & ny , & nz1, & nz2, & circular_bump_radius, & sine_bump_zmax, & stretching_factor_z, & generate_ugrid_file, & ugrid_file_unformatted, & generate_tec_file_b, & generate_tec_file_v, & generate_su2grid_file, & generate_vtk_file !------------------------------------------------------------------------ ! End of List of all input variables: !------------------------------------------------------------------------ contains !***************************************************************************** !* Read input_parameters in the input file: file name = namelist_file !***************************************************************************** subroutine read_nml_input_parameters(namelist_file) implicit none character(9), intent(in) :: namelist_file integer :: os write(*,*) "**************************************************************" write(*,*) " List of namelist variables and their values" write(*,*) open(unit=10,file=namelist_file,form='formatted',status='old',iostat=os) read(unit=10,nml=input_parameters) write(*,nml=input_parameters) ! Print the namelist variables. close(10) if ( bump_type==3 .and. mod(nx2,2)/=0 ) then write(*,*) " For a triangular bump, nx2 must be an even number... Stop." stop endif end subroutine read_nml_input_parameters end module input_parameter_module !------------------------------------------------------------------------------- !------------------------------------------------------------------------------- !------------------------------------------------------------------------------- ! ! End of input parameter module ! !------------------------------------------------------------------------------- !------------------------------------------------------------------------------- !------------------------------------------------------------------------------- !------------------------------------------------------------------------------- !------------------------------------------------------------------------------- !------------------------------------------------------------------------------- !******************************************************************************** ! Input : See "Input parameters" section in the code below. ! ! Output: bump3d.ugrid ! Grid file (in .ugrid format) ! bump3d.mapbc ! Boundary condition file ! bump3d.su2 ! Grid file (in .su2 format) ! bump3d.vtk ! Grid file (in .vtk format) ! bump3d_tecplot.dat ! Volume grid file for viewing by Tecplot ! bump3d_boundary_tecplot.dat !Boundary grid file for viewing by Tecplot ! !******************************************************************************** program edu3d_bump3d use input_parameter_module implicit none real(dp), parameter :: one = 1.0_dp real(dp), parameter :: pi = 3.141592653589793238_dp ! 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, tet, tria, quad ! Element connectivities ! Number of nodes in each coordinate direction integer :: nx, nz ! Number of prisms, tetrahedra, triangles, quadrilaterals integer :: nprz, ntet, ntria, nquad ! Number of nodes integer :: nnodes ! Local numbering of a hexahedron integer :: hex_loc(8) ! Local variables integer :: i, j, k, os real(dp) :: zmaxf, zminf real(dp) :: xc, dzc real(dp) :: sf real(dp) :: smin, smax, s, xm character(80) :: filename_tec_b character(80) :: filename_su2, filename_vtk integer :: nb character(80), dimension(: ), pointer :: bnames !Boundary names integer :: nhex integer , allocatable, dimension(:,:) :: hex !******************************************************************************* ! Unit cube to be gridded. ! ! z ! ^ Circuler or sine^4 ! | <-----------> ! | ------------------_____________----------------------- zmax ! | | | ! | | upper part of the domain | ! | |. . . . . . . . . . . . . . . . . . . . . . . . . . | ! | | lower part of the domain | ! | | | ! | ------------------------------------------------------ zmin ! | . . ! -------------------------------------------------------------->x ! xmin xle xte xmax ! <----- nx1 -----> <-- nx2 --> <-------- nx3 --------> ( <- # of cells ) ! ! e.g., nx1=5, nx2=4 ! ! o---o--o--o---o--o--o--o--o--o--o--o-........ ! 1 2 3 4 5 6 7 8 9 10 ... ! nx1+1 nx1+nx2+1 ! !******************************************************************************* !-------------------------------------------------------------------- ! Read the input file named "input.nml": write(*,*) "Reading the input file: input.nml..... " write(*,*) call read_nml_input_parameters('input.nml') write(*,*) !-------------------------------------------------------------------- nx = nx1 + nx2 + nx3 + 1 nz = nz1 + nz2 write(*,*) write(*,*) " nx = ", nx write(*,*) " ny = ", ny write(*,*) " nz = ", nz write(*,*) !-------------------------------------------------------------------- ! Allocate structured arrays allocate( x(nx+1,ny+1,nz+1), y(nx+1,ny+1,nz+1), z(nx+1,ny+1,nz+1) ) !-------------------------------------------------------------------- ! Let's begin by generating a 2D grid on the (x,z)-plane at y=0 j = 1 !<- at y=0 do i = 1, nx+1 do k = 1, nz+1 !---------------------------------------------------------------------- ! Left side of the bump if ( i < nx1+1 ) then x(i,j,k) = xmin + real(i-1,dp) * (xle-xmin)/real(nx1,dp) !Uniform spacing !---------------------------------------------------------------------- ! Over the bump elseif ( i >= nx1+1 .and. i < nx1+nx2+1 ) then x(i,j,k) = xle + real(i-(nx1+1),dp)*(xte-xle)/real(nx2,dp) !Uniform spacing !---------------------------------------------------------------------- ! Right side of the bump else x(i,j,k) = xte + real(i-(nx1+nx2+1),dp)*(xmax-xte)/real(nx3,dp) !Uniform spacing endif !---------------------------------------------------------------------- !-------------------------------------------------------------------------------- !-------------------------------------------------------------------------------- ! (1)Circular bump !-------------------------------------------------------------------------------- !-------------------------------------------------------------------------------- type_of_bump : if (bump_type == 1) then !---------------------------------------------------------------------- ! Bump is inserted here. The bump is a circular arc. !---------------------------------------------------------------------- if (bump_on_top) then !------------------------------------------- ! Adjust the zmax over the bump: Bump on top. ! x-coodrinate of the center of a circle. xc = (xte+xle)/2.0_dp ! displacement of the z-coodinate of the center of a circle. dzc = sqrt( circular_bump_radius**2 - ( (xte-xle)/2.0_dp )**2 ) if (i > nx1 .and. i < nx1+nx2+1) then zmaxf = (zmax + dzc) - sqrt( circular_bump_radius**2 - (x(i,j,k)-xc)**2 ) !Circular arc. else zmaxf = zmax !Original zmax (no bump) endif !Compute z-coordinates wit the modified zmax. z(i,j,k) = zmin + (zmaxf-zmin)*real(k-1,dp)/real(nz,dp) !(-) bump at top !------------------------------------------- else !------------------------------------------- ! Adjust the zmin over the bump: Bump on bottom. ! x-coodrinate of the center of a circle. xc = (xte+xle)/2.0_dp ! displacement of the z-coodinate of the center of a circle. dzc = sqrt( circular_bump_radius**2 - ( (xte-xle)/2.0_dp )**2 ) if (i > nx1 .and. i < nx1+nx2+1) then zminf = zmin - ( dzc - sqrt( circular_bump_radius**2 - (x(i,j,k)-xc)**2 ) ) !Circular arc. else zminf = zmin !Original zmax (no bump) endif !Compute z-coordinates wit the modified zmax. z(i,j,k) = zminf + (zmax-zminf)*real(k-1,dp)/real(nz,dp) !(-) bump at top !------------------------------------------- endif !-------------------------------------------------------------------------------- !-------------------------------------------------------------------------------- ! (2)Smooth sine^4 bump: ! If xle=0.3, xte=1.2, and sine_bump_zmax = 0.05, then the bump is the same ! as the one described in https://turbmodels.larc.nasa.gov/bump.html. !-------------------------------------------------------------------------------- !-------------------------------------------------------------------------------- elseif (bump_type == 2) then if (bump_on_top) then !------------------------------------------- ! Adjust the zmax over the bump: Bump on top. if (i > nx1 .and. i < nx1+nx2+1) then zmaxf = zmax - sine_bump_zmax*( sin( pi*(x(i,j,k)-xle)/(xte-xle) ) )**4 else zmaxf = zmax !Original zmax (no bump) endif !Compute z-coordinates wit the modified zmax. z(i,j,k) = zmin + (zmaxf-zmin)*real(k-1,dp)/real(nz,dp) !(-) bump at top !------------------------------------------- else !------------------------------------------- ! Adjust the zmin over the bump: Bump on bottom. if (i > nx1 .and. i < nx1+nx2+1) then zminf = zmin + sine_bump_zmax*( sin( pi*(x(i,j,k)-xle)/(xte-xle) ) )**4 else zminf = zmin !Original zmax (no bump) endif !Compute z-coordinates wit the modified zmax. z(i,j,k) = zminf + (zmax-zminf)*real(k-1,dp)/real(nz,dp) !(-) bump at top !------------------------------------------- endif !-------------------------------------------------------------------------------- !-------------------------------------------------------------------------------- ! (3)Triangular bump: ! Triangle is defined wirh xle, xte, and tria_height. !-------------------------------------------------------------------------------- !-------------------------------------------------------------------------------- elseif (bump_type == 3) then xm = 0.5_dp*(xle+xte) !Make sure the x-coordinate is orginally equal to xm at the midpoint of the bump. if (i == nx1+1+nx2/2) x(i,j,k) = xm if (bump_on_top) then !------------------------------------------- ! Adjust the zmax over the bump: Bump on top. if (i > nx1 .and. i < nx1+1+nx2/2) then zmaxf = zmax - tria_height*(x(i,j,k)-xle)/(xm-xle) elseif (i >= nx1+1+nx2/2 .and. i < nx1+nx2+1) then zmaxf = zmax - tria_height*(x(i,j,k)-xte)/(xm-xte) else zmaxf = zmax !Original zmax (no bump) endif !Compute z-coordinates with the modified zmax. z(i,j,k) = zmin + (zmaxf-zmin)*real(k-1,dp)/real(nz,dp) !(-) bump at top !------------------------------------------- else !------------------------------------------- ! Adjust the zmin over the bump: Bump on bottom. if (i > nx1 .and. i < nx1+1+nx2/2) then zminf = zmin + tria_height*(x(i,j,k)-xle)/(xm-xle) elseif (i >= nx1+1+nx2/2 .and. i < nx1+nx2+1) then zminf = zmin + tria_height*(x(i,j,k)-xte)/(xm-xte) else zminf = zmin !Original zmax (no bump) endif !Compute z-coordinates with the modified zmax. z(i,j,k) = zminf + (zmax-zminf)*real(k-1,dp)/real(nz,dp) !(-) bump at top !------------------------------------------- endif else write(*,*) " Invalid input: bump_type = ", bump_type write(*,*) " bump_type must be 1 or 2.... Stop." stop endif type_of_bump end do end do !-------------------------------------------------------------------- ! Stretching in z-direction. sf = stretching_factor_z j = 1 !<- at y=0 do i = 1, nx+1 smin = z(i,j, 1) smax = z(i,j,nz+1) do k = 1, nz+1 s = (z(i,j,k)-smin)/(smax-smin) ! Transform z into s = [0,1]. s = (one-exp(sf*s))/(one-exp(sf)) ! Stretching in [0,1] z(i,j,k) = s*(smax-smin) + smin ! Transform back to z. end do end do !-------------------------------------------------------------------- ! Generate points in the 3D domain. do j = 1, ny+1 do i = 1, nx+1 do k = 1, nz+1 x(i,j,k) = x(i,1,k) y(i,j,k) = ymin + (ymax-ymin)*real(j-1,dp)/real(ny,dp) !Uniform spacing in y. z(i,j,k) = z(i,1,k) end do end do end do ! For smooth transition from prismatic layer to the tetrahedral region, ! you may want to apply some stretching in z-direction in the prismatic ! layer. I'll leave it to you. !-------------------------------------------------------------------- ! Construct a node array ! (1)Construct a map first. allocate(ijk2n(nx+1,ny+1,nz+1)) ! (i,j,k) to node map nnodes = 0 do i = 1, nx+1 do j = 1, ny+1 do k = 1, nz+1 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+1 do j = 1, ny+1 do k = 1, nz+1 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: ! Prz: 6 nodes that define the prism. ! Quad: 4 nodes that define the quad face. ! 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 tetrahedra ! - Quad faces are divided into 2 triangles. if (prism_on_top) then !Prisms on top / Tets on the bottom. nprz = 2 * nx*ny* nz2 ntet = 6 * nx*ny*(nz-nz2) ntria = 4 * ( nx*ny + ny*(nz-nz2) + nx*(nz-nz2) ) nquad = 2 * ( nz2*ny + nz2*nx ) else !Tets on top / Prisms on the bottom. nprz = 2 * nx*ny* nz1 ntet = 6 * nx*ny*(nz-nz1) ntria = 4 * ( nx*ny + ny*(nz-nz1) + nx*(nz-nz1) ) nquad = 2 * ( nz1*ny + nz1*nx ) endif write(*,*) " nz1 = ", nz1 write(*,*) write(*,*) " prism = ", nprz write(*,*) " tetra = ", ntet write(*,*) " tria = ", ntria write(*,*) " quad = ", nquad write(*,*) allocate( prz(max(nprz,1),8) ) allocate( tet(max(ntet,1),4) ) allocate( tria(ntria,5), quad(nquad,5) ) !-------------------------------------------------------------------- ! Reset ntet and ntria, nhex and nquad. nprz = 0 ntet = 0 ntria = 0 nquad = 0 !-------------------------------------------------------------------- ! Loop over (i,j,k) and construct connectivity data. k_loop : do k = 1, nz i_loop : do i = 1, nx j_loop : do j = 1, ny ! 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 numbers 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 defined by the ! numbers {1,2,3,4,5,6,7,8}. !--- Hexahedron defined by the nodes 1265-4378 !-------------------------------------------------------------------------- ! Lower Part: Divide the hexahedra into 2 prisms as follows. !-------------------------------------------------------------------------- if ( (.not.prism_on_top .and. k <= nz1) .or. (prism_on_top .and. k > nz1) ) then ! 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 hexa is adjacent to the boundary: ! Note: Nodes are ordered so that the boundary face points INWARD. ! ! 1. Left boundary face (x=xmin) : 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 endif ! 2. Right boundary face (x=xmax) : Quad face if (i == nx) 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=ymin) : 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 endif ! 4. Back boundary face (y=ymax) : Quad face if (j == ny) 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=zmin) : 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 endif ! 6. Top boundary face (z=zmax) : Tria face if (k == nz) 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 else !-------------------------------------------------------------------------- ! Upper Part: 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 hexa is adjacent to the boundary: ! Note: Nodes are ordered so that the boundary face points INWARD. ! ! 1. Left boundary face (x=xmin) 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 endif ! 2. Right boundary face (x=xmax) if (i == nx) 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=ymin) 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 endif ! 4. Back boundary face (y=ymax) if (j == ny) 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=zmin) 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 endif ! 6. Top boundary face (z=zmax) if (k == nz) 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 endif end do j_loop end do i_loop end do k_loop !******************************************************************************* ! Separate the boundary if requested. ! ! Up to this point, the domain boundary is split into 6 parts: ! ! Left = 1 ! Right = 2 ! Front = 3 ! Back = 4 ! Bottom = 5 ! Top = 6 ! ! Note: We only have triangles on top and bottom boundaries. ! ! We split the bottom or top into left/bump/right or in the case of triangular ! bump, left/left_slop/right_slope/right. ! !******************************************************************************* !---------------------------------------------------------- ! Triangular bump triangular_bump : if (bump_type == 3) then !------------------------------------------------------------------------- !(1)Bump on the top !------------------------------------------------------------------------- if (bump_on_top) then !Split triangles on the top ( tria(i,4)=6 ) do i = 1, ntria if (tria(i,4) == 6) then xc = ( xyz(tria(i,1),1) + xyz(tria(i,2),1) + xyz(tria(i,3),1) )/3.0_dp !Left of the bump if (xc < xle) then tria(i,4) = 6 !Over the left of triangle elseif (xc > xle .and. xc < xm) then tria(i,4) = 7 !Over the right of triangle elseif (xc > xm .and. xc < xte) then tria(i,4) = 8 !Right of the bump else tria(i,4) = 9 endif endif end do !------------------------------------------------------------------------- !(2)Bump on the bottom !------------------------------------------------------------------------- else do i = 1, ntria if (tria(i,4) == 6) tria(i,4) = 16 !Temporarily save off (triangles on top) end do !Split triangles on the bottom ( tria(i,4)=5 ) do i = 1, ntria if (tria(i,4) == 5) then xc = ( xyz(tria(i,1),1) + xyz(tria(i,2),1) + xyz(tria(i,3),1) )/3.0_dp !Left of the bump if (xc < xle) then tria(i,4) = 6 !Over the left of triangle elseif (xc > xle .and. xc < xm) then tria(i,4) = 7 !Over the right of triangle elseif (xc > xm .and. xc < xte) then tria(i,4) = 8 !Right of the bump else tria(i,4) = 9 endif endif end do do i = 1, ntria if (tria(i,4) == 16) tria(i,4) = 5 !Re-assign (triangles on top) end do endif !------------------------------------------------------------------------- ! !------------------------------------------------------------------------- !---------------------------------------------------------- ! Others else if (separate_bump) then !------------------------------------------------------------------------- !(1)Bump on the top !------------------------------------------------------------------------- if (bump_on_top) then !Split triangles on the top ( tria(i,4)=6 ) do i = 1, ntria if (tria(i,4) == 6) then xc = ( xyz(tria(i,1),1) + xyz(tria(i,2),1) + xyz(tria(i,3),1) )/3.0_dp !Left of the bump if (xc < xle) then tria(i,4) = 6 !Over the bump elseif (xc > xle .and. xc < xte) then tria(i,4) = 7 !Right of the bump else tria(i,4) = 8 endif endif end do !------------------------------------------------------------------------- !(2)Bump on the bottom !------------------------------------------------------------------------- else do i = 1, ntria if (tria(i,4) == 6) tria(i,4) = 16 !Temporarily save off (triangles on top) end do !Split triangles on the bottom ( tria(i,4)=5 ) do i = 1, ntria if (tria(i,4) == 5) then xc = ( xyz(tria(i,1),1) + xyz(tria(i,2),1) + xyz(tria(i,3),1) )/3.0_dp !Left of the bump if (xc < xle) then tria(i,4) = 6 !Over the bump elseif (xc > xle .and. xc < xte) then tria(i,4) = 7 !Right of the bump else tria(i,4) = 8 endif endif end do do i = 1, ntria if (tria(i,4) == 16) tria(i,4) = 5 !Re-assign (triangles on top) end do endif !------------------------------------------------------------------------- ! !------------------------------------------------------------------------- endif endif triangular_bump !---------------------------------------------------------- !******************************************************************************* write(*,*) write(*,*) " Nodes = ", nnodes write(*,*) " Prisms = ", nprz write(*,*) " Tetrahedra = ", ntet write(*,*) " Quadrilaterals = ", nquad write(*,*) " Triangles = ", ntria !******************************************************************************* ! Write a UGRID file !******************************************************************************* write_ugrid : if (generate_ugrid_file) then !------------------------------------------------------------------------------- ! Unformatted !------------------------------------------------------------------------------- if (ugrid_file_unformatted) then if ( big_endian_io(9999) ) then open(unit=2, file="bump3d.b8.ugrid", form='unformatted',access="stream",& status='unknown', iostat=os ) else open(unit=2, file="bump3d.lb8.ugrid", form='unformatted',access="stream",& status='unknown', iostat=os ) end if ! #nodes, #tri_faces, #quad_faces, #tetra, #pyr, #prz, #hex write(2) nnodes, ntria, nquad, ntet, 0, nprz, 0 ! Nodal coordinates do i = 1, nnodes write(2) (xyz(i,j), j=1,3) end do ! Triangular boundary faces do i = 1, ntria write(2) tria(i,1), tria(i,2), tria(i,3) end do ! Quad boundary faces do i = 1, nquad write(2) 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) tria(i,4) end do do i = 1, nquad write(2) quad(i,5) end do ! Tetrahedra do i = 1, ntet write(2) tet(i,1), tet(i,2), tet(i,3), tet(i,4) end do ! Prisms do i = 1, nprz write(2) prz(i,1), prz(i,2), prz(i,3), prz(i,4), prz(i,5), prz(i,6) end do close(2) !------------------------------------------------------------------------------- ! Formatted !------------------------------------------------------------------------------- else open(unit=2, file="bump3d.ugrid", form='formatted',access="stream",& status='unknown', iostat=os ) write(2,'(7i20)') nnodes, ntria, nquad, ntet, 0, nprz, 0 ! Nodes do i = 1, nnodes write(2,'(3es26.15)') (xyz(i,j), j=1,3) end do ! Triangular faces = ntri if (ntria > 0) then do i = 1, ntria write(2,'(3i20)') tria(i,1), tria(i,2), tria(i,3) end do endif ! Quad faces = nquad if (nquad > 0) then do i = 1, nquad write(2,'(4i20)') quad(i,1), quad(i,2), quad(i,3), quad(i,4) end do endif ! Face tag if (ntria > 0) then do i = 1, ntria write(2,'(i10)') tria(i,4) end do endif if (nquad > 0) then do i = 1, nquad write(2,'(i10)') quad(i,5) end do endif ! tet if (ntet > 0) then do i = 1, ntet write(2,'(4i20)') tet(i,1), tet(i,2), tet(i,3), tet(i,4) end do endif ! Prism if (nprz > 0) then do i = 1, nprz write(2,'(6i20)') prz(i,1), prz(i,2), prz(i,3), & prz(i,4), prz(i,5), prz(i,6) end do endif endif endif write_ugrid !******************************************************************************* ! Write a boundary condition map file: boundary marks ! Note: Set appropriate boundary condition numbers in this file. !******************************************************************************* open(unit=3, file="bump3d.mapbc", status="unknown", iostat=os) if (bump_type == 3) then !Boundary information nb = 9 allocate( bnames(nb) ) if (bump_on_top) then write(3,'(a)') "9 !Number of boundary parts (boundary conditions)" write(3,'(a)') "1 5050 xmin !Viscous wall in FUN3D" write(3,'(a)') "2 5051 xmax !Outflow" write(3,'(a)') "3 6662 ymin !Y-Symmetry" write(3,'(a)') "4 6662 ymax !Y-Symmetry" write(3,'(a)') "5 4000 zmin !Wall" write(3,'(a)') "6 6663 left_sym !Z-Symmetry" write(3,'(a)') "7 4000 bump_left !Left part of triangular bump" write(3,'(a)') "8 4000 bump_right !Right part of triangular bump" write(3,'(a)') "9 6663 right_sym !Z-Symmetry" bnames(1) = "xmin" bnames(2) = "xmax" bnames(3) = "ymin" bnames(4) = "ymax" bnames(5) = "zmin" bnames(6) = "left_sym" bnames(7) = "bump_left" bnames(8) = "bump_right" bnames(9) = "right_sym" else write(3,'(a)') "9 !Number of boundary parts (boundary conditions)" write(3,'(a)') "1 5050 xmin !Viscous wall in FUN3D" write(3,'(a)') "2 5051 xmax !Outflow" write(3,'(a)') "3 6662 ymin !Y-Symmetry" write(3,'(a)') "4 6662 ymax !Y-Symmetry" write(3,'(a)') "5 4000 zmax !Wall" write(3,'(a)') "6 6663 left_sym !Z-Symmetry" write(3,'(a)') "7 4000 bump_left !Left part of triangular bump" write(3,'(a)') "8 4000 bump_right !Right part of triangular bump" write(3,'(a)') "9 6663 right_sym !Z-Symmetry" bnames(1) = "xmin" bnames(2) = "xmax" bnames(3) = "ymin" bnames(4) = "ymax" bnames(5) = "zmax" bnames(6) = "left_sym" bnames(7) = "bump_left" bnames(8) = "bump_right" bnames(9) = "right_sym" endif else if (separate_bump) then !Boundary information nb = 8 allocate( bnames(nb) ) if (bump_on_top) then write(3,'(a)') "8 !Number of boundary parts (boundary conditions)" write(3,'(a)') "1 5050 xmin !Viscous wall in FUN3D" write(3,'(a)') "2 5051 xmax !Outflow" write(3,'(a)') "3 6662 ymin !Y-Symmetry" write(3,'(a)') "4 6662 ymax !Y-Symmetry" write(3,'(a)') "5 4000 zmin !Wall" write(3,'(a)') "6 6663 left_sym !Z-Symmetry" write(3,'(a)') "7 4000 bump !Wall" write(3,'(a)') "8 6663 right_sym !Z-Symmetry" bnames(1) = "xmin" bnames(2) = "xmax" bnames(3) = "ymin" bnames(4) = "ymax" bnames(5) = "zmin" bnames(6) = "left_sym" bnames(7) = "bump" bnames(8) = "right_sym" else write(3,'(a)') "8 !Number of boundary parts (boundary conditions)" write(3,'(a)') "1 5050 xmin !Viscous wall in FUN3D" write(3,'(a)') "2 5051 xmax !Outflow" write(3,'(a)') "3 6662 ymin !Y-Symmetry" write(3,'(a)') "4 6662 ymax !Y-Symmetry" write(3,'(a)') "5 4000 zmax !Wall" write(3,'(a)') "6 6663 left_sym !Z-Symmetry" write(3,'(a)') "7 4000 bump !Wall" write(3,'(a)') "8 6663 right_sym !Z-Symmetry" bnames(1) = "xmin" bnames(2) = "xmax" bnames(3) = "ymin" bnames(4) = "ymax" bnames(5) = "zmax" bnames(6) = "left_sym" bnames(7) = "bump" bnames(8) = "right_sym" endif else !Boundary information nb = 6 allocate( bnames(nb) ) write(3,'(a)') "6 !Number of boundary parts (boundary conditions)" write(3,'(a)') "1 5050 xmin !Viscous wall in FUN3D" write(3,'(a)') "2 5051 xmax !Outflow" write(3,'(a)') "3 6662 ymin !Y-Symmetry" write(3,'(a)') "4 6662 ymax !Y-Symmetry" write(3,'(a)') "5 4000 zmin !Wall" write(3,'(a)') "6 4000 zmax !Wall" bnames(1) = "xmin" bnames(2) = "xmax" bnames(3) = "ymin" bnames(4) = "ymax" bnames(5) = "zmin" bnames(6) = "zmax" endif endif close(3) !******************************************************************************* ! Write a Tecplot volume grid file for viewing. Just for viewing. !******************************************************************************* write_volume_tec_file : if (generate_tec_file_v) then open(unit=1, file="bump3d_tecplot.dat", status="unknown", iostat=os) write(1,*) 'title = "tet grid"' write(1,*) 'variables = "x","y","z"' ! Tetra Zone write(1,*) 'zone n=', nnodes,',e=', ntet,' , et=tetrahedron, f=fepoint' do i = 1, nnodes write(1,'(3ES20.10)') (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 ! Prism zone write(1,*) 'zone n=', nnodes,',e=', nprz,' , et=brick, f=fepoint' do i = 1, nnodes write(1,'(3es20.10)') (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) endif write_volume_tec_file !******************************************************************************* ! Write a Tecplot boundary grid file for viewing. Just for viewing. ! Boundary parts are written separately. !******************************************************************************* if (generate_tec_file_b) then filename_tec_b = "bump3d_boundary_tecplot.dat" call write_tecplot_boundary_file(filename_tec_b, nnodes, & xyz(:,1),xyz(:,2),xyz(:,3), & ntria,tria, nquad,quad, nb, bnames ) endif !------------------------------------------------------------------------------- ! ! Preparation for calling vtk and su2. ! !------------------------------------------------------------------------------- !No hex. nhex = 0 allocate( hex(1,1) ) !------------------------------------------------------------------------------- ! ! Generate a .vtk file. ! !------------------------------------------------------------------------------- write_vtk_grid : if (generate_vtk_file) then filename_vtk = "bump3d.vtk" write(*,*) write(*,*) " Writing .vtk file... ", trim(filename_vtk) write(*,*) call write_vtk_file(filename_vtk, nnodes,xyz(:,1),xyz(:,2),xyz(:,3), ntet,tet, & nprz,prz, nhex,hex ) endif write_vtk_grid !------------------------------------------------------------------------------- ! ! Generate a .su2 file for SU2 solver. ! !------------------------------------------------------------------------------- write_su2_grid : if (generate_su2grid_file) then filename_su2 = "bump3d.su2" write(*,*) write(*,*) " Writing .su2 file... ", trim(filename_su2) write(*,*) call write_su2grid_file(filename_su2, nnodes,xyz(:,1),xyz(:,2),xyz(:,3), ntet,tet, nprz,prz, & nhex,hex, ntria,tria, nquad,quad, nb, bnames ) endif write_su2_grid !------------------------------------------------------------------------------- ! !------------------------------------------------------------------------------- write(*,*) write(*,*) " Grid files have been successfully generated:" if (generate_ugrid_file) write(*,*) " --- bump3d.ugrid (boundary faces are ordered pointing inward.)" if (generate_su2grid_file) write(*,*) " --- bump3d.su2" if (generate_vtk_file) write(*,*) " --- bump3d.vtk" write(*,*) " --- bump3d.mapbc" write(*,*) write(*,*) " Tecplot files have been successfully generated:" if (generate_tec_file_v) write(*,*) " --- bump3d_tecplot.dat" if (generate_tec_file_b) write(*,*) " --- bump3d_boundary_tecplot.dat" write(*,*) stop contains !******************************************************************************** ! Find out big_endian_io. !******************************************************************************** function big_endian_io( opt_unit ) integer, intent(in) :: opt_unit logical :: big_endian_io ! one-byte integer integer, parameter :: i1 = selected_int_kind(2) ! two-byte integer integer, parameter :: i2 = selected_int_kind(4) integer(i1) :: byte_one, byte_two ! 00000000 00000001 big-endian binary integer(i2) :: two_byte_int = 1_i2 open(opt_unit,status='scratch',form='unformatted') write( opt_unit) two_byte_int rewind(opt_unit) read( opt_unit) byte_one, byte_two close(opt_unit) big_endian_io = ( byte_one == 0 .and. byte_two == 1 ) end function big_endian_io !******************************************************************************** !******************************************************************************* ! This subroutine writes a su2 grid file. ! ! Note: Nodes -> i = 0,1,2,...; Elements -> i = 0,1,2,... ! !******************************************************************************* subroutine write_su2grid_file(filename, nnodes,xp,yp,zp, ntet,tet, nprs,prs, & nhex,hex, ntria,tria, nquad,quad, & nb, bnames ) character(80), intent(in) :: filename integer , intent(in) :: nnodes integer , intent(in) :: ntet, nprs, nhex integer , intent(in) :: ntria, nquad real(dp) , dimension(: ), intent(in) :: xp, yp, zp integer , dimension(:,:), intent(in) :: tria integer , dimension(:,:), intent(in) :: quad integer , dimension(:,:), intent(in) :: tet integer , dimension(:,:), intent(in) :: prs integer , dimension(:,:), intent(in) :: hex integer , intent(in) :: nb character(80), dimension(:) , intent(in) :: bnames integer :: k, itag, i open(unit=7, file=filename, status="unknown", iostat=os) write(7,*) "%" write(7,*) "% Problem dimension" write(7,*) "%" write(7,5) 3 5 format('NDIME= ',i12) write(7,*) "%" write(7,*) "% Inner element connectivity" k = ntet + nprs + nhex write(7,10) k 10 format('NELEM= ',i12) k = 0 !------------------------------------------------------------------------- ! Elements ! tet if (ntet > 0) then do i = 1, ntet write(7,'(6i20)') 10, tet(i,1)-1, tet(i,2)-1, tet(i,3)-1, tet(i,4)-1, k k = k + 1 end do endif ! Prism: Orietation is reversed (See VTK format). if (nprs > 0) then do i = 1, nprs write(7,'(8i20)') 13, prs(i,3)-1, prs(i,2)-1, prs(i,1)-1, & prs(i,6)-1, prs(i,5)-1, prs(i,4)-1, k k = k + 1 end do endif ! Hex if (nhex > 0) then do i = 1, nhex write(7,'(10i20)') 12, hex(i,1)-1, hex(i,2)-1, hex(i,3)-1, hex(i,4)-1, & hex(i,5)-1, hex(i,6)-1, hex(i,7)-1, hex(i,8)-1, k k = k + 1 end do endif write(*,*) " --- elm check", ntet + nprs + nhex, k !-------------------------------------------------------------------------- ! Nodes write(7,*) "%" write(7,*) "% Node coordinates" write(7,*) "%" write(7,20) nnodes 20 format('NPOIN= ', i12) k = 0 ! Nodes do i = 1, nnodes write(7,'(3es26.15,i20)') xp(i), yp(i), zp(i) k = k + 1 end do write(*,*) " --- node check", nnodes, k write(*,*) !-------------------------------------------------------------------------- ! Boundary write(7,*) "%" write(7,*) "% Boundary elements" write(7,*) "%" write(7,30) nb 30 format('NMARK= ',i12) 40 format('MARKER_TAG= ',a) 50 format('MARKER_ELEMS= ', i12) boundary_loop : do itag = 1, nb write(7,40) trim( bnames(itag) ) write(*, *) ' MARKER_TAG = ', trim( bnames(itag) ) k = 0 ! Triangular faces = ntri if (ntria > 0) then do i = 1, ntria if ( tria(i,4) == itag ) k = k + 1 end do endif ! Quad faces = nquad if (nquad > 0) then do i = 1, nquad if ( quad(i,5) == itag ) k = k + 1 end do endif write(7,50) k ! Triangular faces = ntri if (ntria > 0) then do i = 1, ntria if ( tria(i,4) == itag ) write(7,'(4i20)') 5, tria(i,1)-1, tria(i,2)-1, tria(i,3)-1 end do endif ! Quad faces = nquad if (nquad > 0) then do i = 1, nquad if ( quad(i,5) == itag ) write(7,'(5i20)') 9, quad(i,1)-1, quad(i,2)-1, quad(i,3)-1, quad(i,4)-1 end do endif end do boundary_loop !-------------------------------------------------------------------------- !-------------------------------------------------------------------------- close(7) end subroutine write_su2grid_file !******************************************************************************** !******************************************************************************* ! This subroutine writes a .vtk file for the grid whose name is defined by ! filename_vtk. ! ! Identifier: ! Line 3 ! Triangle 5 ! Quadrilateral 9 ! Tetrahedral 10 ! Hexahedral 12 ! Prism 13 ! Pyramid 14 ! ! Note: This version is only for tet, prs, and hex. Need to add pyramids and others ! if needed. ! ! Use Paraview to read .vtk and visualize it. https://www.paraview.org ! ! Search in Google for 'vkt format' to learn .vtk file format. !******************************************************************************* subroutine write_vtk_file(filename, nnodes,xp,yp,zp, ntet,tet, nprs,prs, nhex,hex ) implicit none character(80), intent(in) :: filename integer , intent(in) :: nnodes real(dp) , dimension(: ), intent(in) :: xp, yp, zp integer , intent(in) :: ntet, nprs, nhex integer , dimension(:,:), intent(in) :: tet integer , dimension(:,:), intent(in) :: prs integer , dimension(:,:), intent(in) :: hex !Local variables integer :: i, j, os !------------------------------------------------------------------------------ !------------------------------------------------------------------------------ !------------------------------------------------------------------------------ !Open the output file. open(unit=8, file=filename, status="unknown", iostat=os) !--------------------------------------------------------------------------- ! Header information write(8,'(a)') '# vtk DataFile Version 3.0' write(8,'(a)') filename write(8,'(a)') 'ASCII' write(8,'(a)') 'DATASET UNSTRUCTURED_GRID' !--------------------------------------------------------------------------- ! Nodal information ! ! Note: These nodes i=1,nnodes are interpreted as i=0,nnodes-1 in .vtk file. ! So, later below, the connectivity list for tria and quad will be ! shifted by -1. write(8,*) 'POINTS ', nnodes, ' double' do j = 1, nnodes write(8,'(3es25.15)') xp(j), yp(j), zp(j) end do !--------------------------------------------------------------------------- ! Cell information. !CELLS: # of total cells (ntet+nprs+nhex), total size of the cell list. write(8,'(a,i12,i12)') 'CELLS ',ntet+nprs+nhex, (4+1)*ntet + (6+1)*nprs + (8+1)*nhex ! Note: The latter is the number of integer values written below as data. ! 5 for tets (# of vertices + 4 vertices), ! 7 for prisms (# of vertices + 6 vertices), ! 9 for hexa (# of vertices + 8 vertices). !--------------------------------- ! 2.1 List of tets if (ntet > 0) then do i = 1, ntet write(8,'(a,4i12)') '4', tet(i,1)-1, tet(i,2)-1, tet(i,3)-1 ! -1 since VTK reads the nodes as 0,1,2,3,..., not 1,2,3,.. end do endif !--------------------------------- ! 2.2 List of prisms if (nprs > 0) then do i = 1, nprs write(8,'(a,6i12)') '6', prs(i,3)-1, prs(i,2)-1, prs(i,1)-1, & prs(i,6)-1, prs(i,5)-1, prs(i,4)-1 ! -1 since VTK reads the nodes as 0,1,2,3,..., not 1,2,3,.. end do endif !--------------------------------- ! 2.3 List of hexa if (nhex > 0) then do i = 1, nhex write(8,'(a,8i12)') '8', hex(i,1)-1, hex(i,2)-1, hex(i,3)-1, hex(i,4)-1, & hex(i,5)-1, hex(i,6)-1, hex(i,7)-1, hex(i,8)-1 ! -1 since VTK reads the nodes as 0,1,2,3,..., not 1,2,3,.. end do endif !--------------------------------------------------------------------------- ! Cell type information. !# of all cells write(8,'(a,i11)') 'CELL_TYPES ', ntet+nprs+nhex !Tetrahedron is classified as the cell type 10 in the .vtk format. if (ntet > 0) then do i = 1, ntet write(8,'(i3)') 10 end do endif !Prism is classified as the cell type 13 in the .vtk format. if (nprs > 0) then do i = 1, nprs write(8,'(i3)') 13 end do endif !Hexahedron is classified as the cell type 12 in the .vtk format. if (nhex > 0) then do i = 1, nhex write(8,'(i3)') 12 end do endif !--------------------------------------------------------------------------- close(8) end subroutine write_vtk_file !******************************************************************************** !******************************************************************************* ! ! This subroutine writes a Tecplot file for boundaries. ! !******************************************************************************* subroutine write_tecplot_boundary_file(filename_tecplot_b, nnodes,x,y,z, & ntria,tria, nquad,quad, nb, bnames ) Implicit none character(80), intent(in) :: filename_tecplot_b integer , intent(in) :: nnodes integer , intent(in) :: ntria, nquad real(dp) , dimension(: ), intent(in) :: x, y, z integer , dimension(:,:), intent(in) :: tria integer , dimension(:,:), intent(in) :: quad integer , intent(in) :: nb character(80), dimension(:) , intent(in) :: bnames !Input character(80) :: boundary_name integer :: nelms integer :: nnodes_loc, i, j, k, ib, os integer, dimension(:,:), allocatable :: g2l integer, dimension(: ), allocatable :: l2g write(*,*) write(*,*) ' Tecplot boundary file = ', trim(filename_tecplot_b) write(*,*) !Allocate global-to-local and local-to-global arrays. allocate(g2l(nnodes,2)) allocate(l2g(nnodes )) open(unit=7, file=filename_tecplot_b, status="unknown", iostat=os) write(7,*) 'TITLE = "GRID"' write(7,*) 'VARIABLES = "x","y","z"' !Write boundaries separately. !--------------------------------------------------------------------------- ! Triangular boundaries. tria_boundary_loop : do ib = 1, nb nelms = 0 nnodes_loc = 0 g2l = -1 l2g = 0 !To save memory, we write the boundary data based on the data !local to each boundary. Local numbers are assigned to each !node on the boundary, which is accessed by g2l (global to local). do i = 1, ntria if ( tria(i,4) == ib ) then nelms = nelms + 1 do k = 1, 3 if (g2l(tria(i,k),2) == -1) then nnodes_loc = nnodes_loc + 1 g2l(tria(i,k),1) = nnodes_loc g2l(tria(i,k),2) = 100 !Flag to indicate the node is recorded. l2g(nnodes_loc) = tria(i,k) endif end do endif end do !So, there are nnodes_loc nodes and nelms elements on this boundary (ib). if (nelms > 0) then boundary_name = trim(bnames(ib)) // " (triangles)" write(7,*) 'ZONE T="', trim(boundary_name), '" N=', nnodes_loc, & ', E=', nelms, ' , ET=quadrilateral, F=FEPOINT' !Write only the local nodes. do j = 1, nnodes_loc i = l2g(j) write(7,'(3es25.15)') x(i), y(i), z(i) end do !Write only the local elements. do i = 1, ntria if ( tria(i,4) == ib ) write(7,'(4i10)') g2l(tria(i,1),1), g2l(tria(i,2),1), & g2l(tria(i,3),1), g2l(tria(i,3),1) end do endif end do tria_boundary_loop !--------------------------------------------------------------------------- !--------------------------------------------------------------------------- ! quad boundaries. if (nquad > 0) then quad_boundary_loop : do ib = 1, nb nelms = 0 nnodes_loc = 0 g2l = -1 l2g = 0 !To save memory, we write the boundary data based on the data !local to each boundary. Local numbers are assigned to each !node on the boundary, which is accessed by g2l (global to local). do i = 1, nquad if ( quad(i,5) == ib ) then nelms = nelms + 1 do k = 1, 4 if (g2l(quad(i,k),2) == -1) then nnodes_loc = nnodes_loc + 1 g2l(quad(i,k),1) = nnodes_loc g2l(quad(i,k),2) = 100 !Flag to indicate the node is recorded. l2g(nnodes_loc) = quad(i,k) endif end do endif end do !So, there are nnodes_loc nodes and nelms elements on this boundary (ib). if (nelms > 0) then boundary_name = trim(bnames(ib)) // " (quads)" write(7,*) 'ZONE T="', trim(boundary_name), '" N=', nnodes_loc, & ', E=', nelms, ' , ET=quadrilateral, F=FEPOINT' !Write only the local nodes. do j = 1, nnodes_loc i = l2g(j) write(7,'(3es25.15)') x(i), y(i), z(i) end do !Write only the local elements. do i = 1, nquad if ( quad(i,5) == ib ) write(7,'(4i10)') g2l(quad(i,1),1), g2l(quad(i,2),1), & g2l(quad(i,3),1), g2l(quad(i,4),1) end do endif end do quad_boundary_loop endif !--------------------------------------------------------------------------- close(7) end subroutine write_tecplot_boundary_file !******************************************************************************** end program edu3d_bump3d