!******************************************************************************** !* Educationally-Designed Unstructured 3D (EDU3D) Code !* !* --- EDU3D tetgrid_cube_ptb !* !* This program generates a irregular 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 7 (2019). !* ! 07-24-19: Added .su2 and .vtk files. ! !* Ver.5 (November 2015) !* Ver.6 !* 10-06-17: Added a subroutine that computes the edge-based data, which are !* used in the edge-based discretization, and verifies them. !* Ver.6.1 !* 10-07-17: Made a bit more efficeint. !* !* !* This F90 program is written and made available for an educational purpose. !* !* This file may be updated in future. !* !* Katate Masatsuka, July 2019. http://www.cfdbooks.com !******************************************************************************** !******************************************************************************** !* Input : nx, ny, nz ! Number of nodes in each coordinate direction !* adjustable_parameter ! Parameter for nodal perturbation !* Larger value gives more perturbation. !* !* Output: tetgrid.ugrid ! Grid file (in .ugrid format) !* tetgrid.mapbc ! Boundary condition file !* tetgrid_boundary_tecplot.dat !Boundary grid file for viewing by Tecplot !* tetgrid.su2 ! SU2 grid file !* tetgrid.vtk ! VTK file !* !* Optional output: !* tetgrid_tecplot.dat ! Volume 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. Make your own. !* !* Note: Connectivity is fixed. Only nodal coordinates are perturbed. !* !* Note: This program checks negative volumes created by perturbation, and !* tries to fix them by reducing the perturbation. !* !******************************************************************************** program tetgrid_cube_perturbed 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, ii ! Variables used for perturbing the nodal coordiantes. real(dp) :: rn, adjustable_parameter integer, allocatable, dimension(:) :: n_node2tet integer, allocatable, dimension(:,:) :: node2tet real(dp) :: volume, x1,x2,x3,x4, y1,y2,y3,y4, z1,z2,z3,z4 real(dp) :: x0, y0, z0, vol_max, vol_min, vol_ave, vol_tot logical :: negative_volume_exists integer :: n_negative_vol_detected, inode, it ! Change this to .true. if you wish to write out a Tecplot file for the entire grid. ! By default, this program only writes a Tecplot file for the boudnary. logical :: write_tecplot_volume_file = .false. character(80) :: filename_su2, filename_vtk integer :: nb character(80), dimension(: ), pointer :: bnames !Boundary names integer :: nprs, nhex, nquad integer , allocatable, dimension(:,:) :: prs, hex, quad !******************************************************************************* ! 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 ! ! !******************************************************************************* !------------------------------------------------ ! Input parameters: ! (1)Number of points in each direction. nx = 15 ny = 15 nz = 15 ! (2)A parameter for nodal perturbation (larger value for larger perturbation). adjustable_parameter = 1.2_dp write(*,*) write(*,*) "--- Input parameters ----------------------" write(*,*) write(*,*) " nx = ", nx write(*,*) " ny = ", ny write(*,*) " nz = ", nz write(*,*) write(*,*) " For nodal perturbation:" write(*,*) " adjustable_parameter = ", adjustable_parameter write(*,*) write(*,*) "-------------------------------------------" write(*,*) !------------------------------------------------ ! 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(*,*) write(*,*) " A regular grid generated: " write(*,*) " - Number of nodes = ", nnodes write(*,*) " - Number of tetrahedra = ", ntet write(*,*) " - Number of tria faces = ", ntria write(*,*) !******************************************************************************* ! Construct data for tetrahedra that share a node. !******************************************************************************* write(*,*) write(*,*) " Constructing data for tetrahedra that share a node....." write(*,*) allocate( n_node2tet(nnodes) ) ! Count the number of tetrahedra sharing the same node: n_node2tet = 0 do i = 1, ntet do k = 1, 4 n_node2tet( tet(i,k) ) = n_node2tet( tet(i,k) ) + 1 !# of tets around the node tet(i,k). end do end do ! Allocate the array to store the neighboring tetrahedra at each node. ! Use the maximum n_node2tet. allocate( node2tet(nnodes, maxval(n_node2tet) ) ) ! Construct the list of tetrahedra around each node. n_node2tet = 0 do i = 1, ntet do k = 1, 4 n_node2tet( tet(i,k) ) = n_node2tet( tet(i,k) ) + 1 !# of tets around the node tet(i,k). node2tet( tet(i,k), n_node2tet( tet(i,k) ) ) = i !Add tet i to the node tet(i,k). end do end do !******************************************************************************* ! Perturb the nodal coordinates. !******************************************************************************* write(*,*) write(*,*) " Start perturbing the nodal coordinates....." write(*,*) ! Randomly perturb the nodal coordinates. ! Check the volume. If negative, reduce the perturbation and try again. n_negative_vol_detected = 0 do i = 1, nx if ( mod(i,10) == 0 ) write(*,*) " Perturbing the nodal coordinates....." do j = 1, ny do k = 1, nz inode = ijk2n(i,j,k) x0 = xyz(inode, 1) y0 = xyz(inode, 2) z0 = xyz(inode, 3) negative_volume_exists = .true. do while (negative_volume_exists) call random_number(rn) if (i/=1 .and. i/=nx) & xyz(inode, 1) = x0 + (rn-0.5_dp)*dx * adjustable_parameter call random_number(rn) if (j/=1 .and. j/=ny) & xyz(inode, 2) = y0 + (rn-0.5_dp)*dy * adjustable_parameter call random_number(rn) if (k/=1 .and. k/=nz) & xyz(inode, 3) = z0 + (rn-0.5_dp)*dz * adjustable_parameter !Compute the volume of all surrounding tetrahedra. negative_volume_exists = .false. do ii = 1, n_node2tet(inode) it = node2tet(inode,ii) x1 = xyz(tet(it,1),1); x2 = xyz(tet(it,2),1); x3 = xyz(tet(it,3),1); x4 = xyz(tet(it,4),1); y1 = xyz(tet(it,1),2); y2 = xyz(tet(it,2),2); y3 = xyz(tet(it,3),2); y4 = xyz(tet(it,4),2); z1 = xyz(tet(it,1),3); z2 = xyz(tet(it,2),3); z3 = xyz(tet(it,3),3); z4 = xyz(tet(it,4),3); volume = tet_volume(x1,x2,x3,x4, y1,y2,y3,y4, z1,z2,z3,z4) if (volume < 0.0_dp) negative_volume_exists = .true. end do !If negative volume is detected, recude the perturbation parameter and try again. if (negative_volume_exists) then adjustable_parameter = adjustable_parameter*0.9_dp n_negative_vol_detected = n_negative_vol_detected + 1 !write(*,*) " Negative volume created. Adjusting a parameter..." endif end do end do end do end do write(*,*) write(*,*) " # of negative volume fix = ", n_negative_vol_detected write(*,*) write(*,*) write(*,*) " Finished perturbing the nodal coordinates....." write(*,*) !******************************************************************************* ! Final volume check !******************************************************************************* vol_max = -1.0_dp vol_min = 1.0e+15_dp vol_ave = 0.0_dp do i = 1, ntet x1 = xyz(tet(i,1),1); x2 = xyz(tet(i,2),1); x3 = xyz(tet(i,3),1); x4 = xyz(tet(i,4),1); y1 = xyz(tet(i,1),2); y2 = xyz(tet(i,2),2); y3 = xyz(tet(i,3),2); y4 = xyz(tet(i,4),2); z1 = xyz(tet(i,1),3); z2 = xyz(tet(i,2),3); z3 = xyz(tet(i,3),3); z4 = xyz(tet(i,4),3); volume = tet_volume(x1,x2,x3,x4, y1,y2,y3,y4, z1,z2,z3,z4) vol_max = max(vol_max,volume) vol_min = min(vol_min,volume) vol_ave = vol_ave + volume end do vol_tot = vol_ave vol_ave = vol_tot / real(ntet,dp) if (vol_min > 0.0_dp) then write(*,*) write(*,*) " Negative volume successfully avoided." write(*,*) endif write(*,*) write(*,*) " Minimum tetrahedral volume = ", vol_min ! <- Should be positive! write(*,*) " Maximum tetrahedral volume = ", vol_max write(*,*) " Average tetrahedral volume = ", vol_ave write(*,*) " Total tetrahedral volume = ", vol_tot, " (which should be 1.0)" write(*,*) if (vol_min < 0.0_dp) then write(*,*) " Tetrahedra with negative volume created... Stop." stop endif !******************************************************************************* ! Write a UGRID file !******************************************************************************* write(*,*) write(*,*) " Writing a UGRID file (tetgrid.ugrid)...." write(*,*) 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. !******************************************************************************* write(*,*) write(*,*) " Writing a BC file (tetgrid.mapbc)...." write(*,*) open(unit=3, file="tetgrid.mapbc", status="unknown", iostat=os) write(3,*) 6, " !Number of boundary parts (boundary conditions)" 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. !******************************************************************************* if (write_tecplot_volume_file) then write(*,*) write(*,*) " Writing a Tecplot volume file (tetgrid_tecplot.dat)...." write(*,*) 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) else write(*,*) " write_tecplot_volume_file = ", write_tecplot_volume_file write(*,*) " OK, skip the Tecplot volume-grid file." endif !******************************************************************************* ! Write a Tecplot boundary grid file for viewing. Just for viewing. !****************************************************************************** write(*,*) write(*,*) " Writing a Tecplot boundary file (tetgrid_boundary_tecplot.dat)...." write(*,*) 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" call compute_edge_based_data(xyz(:,1),xyz(:,2),xyz(:,3),nnodes,ntet,tet,ntria,tria) !------------------------------------------------------------------------------- ! ! Preparation for calling vtk and su2. ! !------------------------------------------------------------------------------- !This is a pure tet grid. So, no prism and hex. nprs = 0 allocate( prs(1,1) ) nhex = 0 allocate( hex(1,1) ) nquad = 0 allocate( quad(1,1) ) !Boundary information nb = 6 allocate( bnames(nb) ) bnames(1) = "xmin" bnames(2) = "xmax" bnames(3) = "ymin" bnames(4) = "ymax" bnames(5) = "zmin" bnames(6) = "zmax" !------------------------------------------------------------------------------- ! ! Generate a .vtk file. ! !------------------------------------------------------------------------------- filename_vtk = "tetgrid.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, & nprs,prs, nhex,hex ) !------------------------------------------------------------------------------- ! ! Generate a .su2 file for SU2 solver. ! !------------------------------------------------------------------------------- filename_su2 = "tetgrid.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, nprs,prs, & nhex,hex, ntria,tria, nquad,quad, nb, bnames ) !------------------------------------------------------------------------------- ! !------------------------------------------------------------------------------- stop contains !******************************************************************************* ! Compute the volume of a tetrahedron defined by 4 vertices: ! ! (x1,y1,z1), (x2,y2,z2), (x3,y3,z3), (x4,y4,z4), ! ! which are ordered as follows: ! ! 1 ! o ! /| . ! / | . ! / | . ! / | . ! 2 o----|-------o 3 ! \ | . ! \ | . ! \ | . ! \|. ! o ! 4 ! ! Note: Volume = volume integral of 1 = 1/3 * volume integral of div(x,y,z) dV ! = surface integral of (x,y,z)*dS ! = sum of [ (xc,yc,zc)*area_vector ] over triangular faces. ! ! where the last step is exact because (x,y,z) vary linearly over the triangle. ! There are other ways to compute the volume, of course. ! !******************************************************************************* function tet_volume(x1,x2,x3,x4, y1,y2,y3,y4, z1,z2,z3,z4) implicit none integer , parameter :: dp = selected_real_kind(15) !Double precision !Input real(dp), intent(in) :: x1,x2,x3,x4, y1,y2,y3,y4, z1,z2,z3,z4 !Output real(dp) :: tet_volume real(dp) :: xc, yc, zc real(dp), dimension(3) :: area integer :: ix=1, iy=2, iz=3 tet_volume = 0.0_dp ! Triangle 1-3-2 !Centroid of the triangular face xc = (x1+x3+x2)/3.0_dp yc = (y1+y3+y2)/3.0_dp zc = (z1+z3+z2)/3.0_dp !Outward normal surface vector area = triangle_area_vector(x1,x3,x2, y1,y3,y2, z1,z3,z2) tet_volume = tet_volume + ( xc*area(ix) + yc*area(iy) + zc*area(iz) ) ! Triangle 1-4-3 !Centroid of the triangular face xc = (x1+x4+x3)/3.0_dp yc = (y1+y4+y3)/3.0_dp zc = (z1+z4+z3)/3.0_dp !Outward normal surface vector area = triangle_area_vector(x1,x4,x3, y1,y4,y3, z1,z4,z3) tet_volume = tet_volume + ( xc*area(ix) + yc*area(iy) + zc*area(iz) ) ! Triangle 1-2-4 !Centroid of the triangular face xc = (x1+x2+x4)/3.0_dp yc = (y1+y2+y4)/3.0_dp zc = (z1+z2+z4)/3.0_dp !Outward normal surface vector area = triangle_area_vector(x1,x2,x4, y1,y2,y4, z1,z2,z4) tet_volume = tet_volume + ( xc*area(ix) + yc*area(iy) + zc*area(iz) ) ! Triangle 2-3-4 !Centroid of the triangular face xc = (x2+x3+x4)/3.0_dp yc = (y2+y3+y4)/3.0_dp zc = (z2+z3+z4)/3.0_dp !Outward normal surface vector area = triangle_area_vector(x2,x3,x4, y2,y3,y4, z2,z3,z4) tet_volume = tet_volume + ( xc*area(ix) + yc*area(iy) + zc*area(iz) ) tet_volume = tet_volume / 3.0_dp end function tet_volume !******************************************************************************* ! Compute the area of a triangle in 3D defined by 3 vertices: ! ! (x1,y1,z1), (x2,y2,z2), (x3,y3,z3), ! ! which is assumed to be ordered clockwise. ! ! 1 2 ! o------------o ! \ . ! \ . ---------> ! \ . ! \ . ! o ! 3 ! ! Note: Area is a vector based on the right-hand rule: ! when wrapping the right hand around the triangle with the fingers in the ! direction of the vertices [1,2,3], the thumb points in the positive ! direction of the area. ! ! Note: Area vector is computed as the cross product of edge vectors [31] and [32]. ! !******************************************************************************* function triangle_area_vector(x1,x2,x3, y1,y2,y3, z1,z2,z3) result(area_vector) implicit none integer , parameter :: dp = selected_real_kind(15) !Double precision !Input real(dp), intent(in) :: x1,x2,x3, y1,y2,y3, z1,z2,z3 !Output real(dp), dimension(3) :: area_vector integer :: ix=1, iy=2, iz=3 area_vector(ix) = 0.5_dp*( (y1-y3)*(z2-z3)-(z1-z3)*(y2-y3) ) area_vector(iy) = 0.5_dp*( (z1-z3)*(x2-x3)-(x1-x3)*(z2-z3) ) area_vector(iz) = 0.5_dp*( (x1-x3)*(y2-y3)-(y1-y3)*(x2-x3) ) end function triangle_area_vector !******************************************************************************* ! Compute the data required for the edeg-based scheme, and check the data. ! ! References: ! ! [1] H. Nishikawa and Y. Liu, " Accuracy-Preserving Source Term Quadrature for ! Third-Order Edge-Based Discretization", Journal of Computational Physics, ! Volume 344, September 2017, Pages 595-622. ! https://doi.org/10.1016/j.jcp.2017.04.075 ! http://hiroakinishikawa.com/My_papers/nishikawa_liu_jcp2017_preprint.pdf ! ! [2] H. Nishikawa, "Beyond Interface Gradient: A General Principle for ! Constructing Diffusion Schemes", AIAA Paper 2010-5093, 40th AIAA Fluid ! Dynamics Conference and Exhibit, Chicago, 2010. ! http://hiroakinishikawa.com/My_papers/nishikawa_AIAA-2010-5093.pdf ! (NOTE: See Appendix B for dual volumes and directed area vectors, boundary.) ! !******************************************************************************* subroutine compute_edge_based_data(x,y,z,nnodes,ntet,tet,nbtria,btria) implicit none integer , parameter :: dp = selected_real_kind(15) !Double precision real(dp), parameter :: zero = 0.0_dp real(dp), parameter :: half = 0.5_dp real(dp), parameter :: one = 1.0_dp real(dp), parameter :: third = 1.0_dp/3.0_dp real(dp), parameter :: fourth = 1.0_dp/4.0_dp real(dp), parameter :: sixth = 1.0_dp/6.0_dp ! Input integer , intent(in) :: nnodes, ntet, nbtria real(dp), dimension(nnodes) , intent(in) :: x real(dp), dimension(nnodes) , intent(in) :: y real(dp), dimension(nnodes) , intent(in) :: z integer , dimension(ntet ,4) , intent(in) :: tet integer , dimension(nbtria,3), intent(in) :: btria ! Dual volume at nodes real(dp), allocatable, dimension(:) :: dual_vol1 real(dp), allocatable, dimension(:) :: dual_vol2 ! Directed area vector at edges. real(dp), allocatable, dimension(:,:) :: directed_area_vector ! Node neighbors of nodes. integer , allocatable, dimension(:) :: nnghbrs integer , allocatable, dimension(:,:) :: nghbr integer , allocatable, dimension(:,:) :: edge_number ! Edge array. integer :: nedges integer , allocatable, dimension(:,:) :: edge ! Cell centroid coordinates. real(dp) :: xm , ym , zm real(dp) :: xc , yc , zc real(dp) :: xcl, ycl, zcl real(dp) :: xcr, ycr, zcr integer :: i_edge integer :: k1, k2 integer :: n1, n2 integer :: kv1, kv2, kv3 integer :: v1, v2, v3 ! Local pointers integer , dimension(4,3) :: loc_nghbr integer , dimension(6,2) :: loc_edge integer , dimension(6,2,3) :: loc_face integer :: left, rght ! Left and right face of an edge in tet. integer :: max_nnghbrs integer :: i, j, k, m integer :: ix, iy, iz logical :: found real(dp), dimension(3) :: edge_vector real(dp), dimension(3) :: directed_area_contribution real(dp) :: dxjk_dot_njk ! Two dual face normals of edge in a tetrahedron. real(dp), dimension(3) :: normal1, normal2 ! Boundary normal real(dp), dimension(3) :: normal_b ! For data check. real(dp), allocatable, dimension(:) :: diff_dual real(dp) :: total_vol1 real(dp) :: total_vol2 real(dp), allocatable, dimension(:,:) :: sum_directed_area_vector real(dp) :: x1,x2,x3,x4, y1,y2,y3,y4, z1,z2,z3,z4 write(*,*) write(*,*) " ---------------------------------------------------------------------" write(*,*) write(*,*) " Edge based data construction and check." write(*,*) left = 1 rght = 2 ix = 1 iy = 2 iz = 3 !------------------------------------------------------------------------------------ !------------------------------------------------------------------------------------ ! Preparation (1): ! ! Define some local pointers within a tetrahedron ordered as below. ! ! 4 ! o ! /| . ! / | . ! / | . ! / | . ! 3 o----|-------o 2 ! \ | . ! \ | . ! \ | . ! \|. ! o ! 1 ! !Neighbor nodes of a given node within a tetrahedron. loc_nghbr(1,1) = 2 !nghbr1 of vertex 1. loc_nghbr(1,2) = 3 !nghbr2 of vertex 1. loc_nghbr(1,3) = 4 !nghbr3 of vertex 1. loc_nghbr(2,1) = 1 !nghbr1 of vertex 2. loc_nghbr(2,2) = 4 !nghbr2 of vertex 2. loc_nghbr(2,3) = 3 !nghbr3 of vertex 2. loc_nghbr(3,1) = 1 !nghbr1 of vertex 3. loc_nghbr(3,2) = 2 !nghbr2 of vertex 3. loc_nghbr(3,3) = 4 !nghbr3 of vertex 3. loc_nghbr(4,1) = 1 !nghbr1 of vertex 4. loc_nghbr(4,2) = 3 !nghbr2 of vertex 4. loc_nghbr(4,3) = 2 !nghbr3 of vertex 4. !Local edges in a tetrahedron: pointing from node1 to node2. loc_edge(1,1) = 1 !node1 of local edge 1 loc_edge(1,2) = 2 !node2 of local edge 1 loc_edge(2,1) = 2 !node1 of local edge 2 loc_edge(2,2) = 3 !node2 of local edge 2 loc_edge(3,1) = 3 !node1 of local edge 3 loc_edge(3,2) = 1 !node2 of local edge 3 loc_edge(4,1) = 1 !node1 of local edge 4 loc_edge(4,2) = 4 !node2 of local edge 4 loc_edge(5,1) = 2 !node1 of local edge 5 loc_edge(5,2) = 4 !node2 of local edge 5 loc_edge(6,1) = 3 !node1 of local edge 6 loc_edge(6,2) = 4 !node2 of local edge 6 !Local adjacent faces of a local edge in a tetrahedron. ! - Left and right face seen from outside of the tetrahedron w.r.t. the edge orientation. ! - Nodes ordered in each face to point outward. loc_face(1,left,1) = 1 !left face of edge 1 (edge nodes: 1-2) loc_face(1,left,2) = 2 !left face of edge 1 (edge nodes: 1-2) loc_face(1,left,3) = 4 !left face of edge 1 (edge nodes: 1-2) loc_face(1,rght,1) = 1 !rght face of edge 1 (edge nodes: 1-2) loc_face(1,rght,2) = 3 !rght face of edge 1 (edge nodes: 1-2) loc_face(1,rght,3) = 2 !rght face of edge 1 (edge nodes: 1-2) loc_face(2,left,1) = 2 !left face of edge 2 (edge nodes: 2-3) loc_face(2,left,2) = 3 !left face of edge 2 (edge nodes: 2-3) loc_face(2,left,3) = 4 !left face of edge 2 (edge nodes: 2-3) loc_face(2,rght,1) = 2 !rght face of edge 2 (edge nodes: 2-3) loc_face(2,rght,2) = 1 !rght face of edge 2 (edge nodes: 2-3) loc_face(2,rght,3) = 3 !rght face of edge 2 (edge nodes: 2-3) loc_face(3,left,1) = 1 !left face of edge 3 (edge nodes: 3-1) loc_face(3,left,2) = 4 !left face of edge 3 (edge nodes: 3-1) loc_face(3,left,3) = 3 !left face of edge 3 (edge nodes: 3-1) loc_face(3,rght,1) = 1 !rght face of edge 3 (edge nodes: 3-1) loc_face(3,rght,2) = 3 !rght face of edge 3 (edge nodes: 3-1) loc_face(3,rght,3) = 2 !rght face of edge 3 (edge nodes: 3-1) loc_face(4,left,1) = 1 !left face of edge 4 (edge nodes: 1-4) loc_face(4,left,2) = 4 !left face of edge 4 (edge nodes: 1-4) loc_face(4,left,3) = 3 !left face of edge 4 (edge nodes: 1-4) loc_face(4,rght,1) = 1 !rght face of edge 4 (edge nodes: 1-4) loc_face(4,rght,2) = 2 !rght face of edge 4 (edge nodes: 1-4) loc_face(4,rght,3) = 4 !rght face of edge 4 (edge nodes: 1-4) loc_face(5,left,1) = 1 !left face of edge 5 (edge nodes: 2-4) loc_face(5,left,2) = 2 !left face of edge 5 (edge nodes: 2-4) loc_face(5,left,3) = 4 !left face of edge 5 (edge nodes: 2-4) loc_face(5,rght,1) = 2 !rght face of edge 5 (edge nodes: 2-4) loc_face(5,rght,2) = 3 !rght face of edge 5 (edge nodes: 2-4) loc_face(5,rght,3) = 4 !rght face of edge 5 (edge nodes: 2-4) loc_face(6,left,1) = 2 !left face of edge 6 (edge nodes: 3-4) loc_face(6,left,2) = 3 !left face of edge 6 (edge nodes: 3-4) loc_face(6,left,3) = 4 !left face of edge 6 (edge nodes: 3-4) loc_face(6,rght,1) = 1 !rght face of edge 6 (edge nodes: 3-4) loc_face(6,rght,2) = 4 !rght face of edge 6 (edge nodes: 3-4) loc_face(6,rght,3) = 3 !rght face of edge 6 (edge nodes: 3-4) !------------------------------------------------------------------------------------ !------------------------------------------------------------------------------------ ! Preparation (2): ! ! Construct a list of neighbor nodes of each node. allocate( nnghbrs(nnodes) ) nnghbrs = 0 allocate( nghbr(nnodes, 50 ) ) !Assume nnghbrs < 51. nghbr = 0 !------------------------------------------------- ! Find edge-connected neghbr nodes within a tetrahedron. do i = 1, ntet !------------------------------------------------- ! Loop over the nodes of tetrahedron i. do k = 1, 4 n1 = tet(i,k) !------------------------------------------------- ! Loop over the rest of nodes. do j = 1, 3 n2 = tet(i,loc_nghbr(k,j)) if ( nnghbrs(n1) == 0 ) then nnghbrs(n1) = 1 nghbr(n1,1) = n2 else !Check if the neighbor is already added. found = .false. do m = 1, nnghbrs(n1) if (n2 == nghbr(n1,m)) then found = .true. exit endif end do !If not added yet, add it to the list. if (.not.found) then nnghbrs(n1) = nnghbrs(n1) + 1 nghbr(n1,nnghbrs(n1)) = n2 endif endif end do !------------------------------------------------- ! End of Loop over the rest of nodes. end do !------------------------------------------------- ! End of Loop over the nodes of tetrahedron i. end do ! End of Find edge-connected neghbr nodes from tets. !------------------------------------------------- !------------------------------------------------------------------------------------ !------------------------------------------------------------------------------------ ! Preparation (3): ! ! Define the edges. Edge points from node1 to node2, ! where node2 > node1 (i.e., smaller node number to larger node number). !This will be used to identify the edge number from two node numbers. max_nnghbrs = maxval(nnghbrs) allocate( edge_number(nnodes, max_nnghbrs ) ) edge_number = 0 !--------------------------------------- ! First, count the number of total edges. nedges = 0 do i = 1, nnodes do k = 1, nnghbrs(i) if ( i < nghbr(i,k) ) then nedges = nedges + 1 endif end do end do !--------------------------------------- ! Allocate and fill the edge array. allocate(edge(nedges,2)) nedges = 0 do i = 1, nnodes do k = 1, nnghbrs(i) if ( i < nghbr(i,k) ) then nedges = nedges + 1 edge(nedges,1) = i edge(nedges,2) = nghbr(i,k) edge_number(i,k) = nedges endif end do end do !------------------------------------------------------------------------------------ !------------------------------------------------------------------------------------ ! Now, we can compute the directed area vector at each edge and the dual volume at ! each node. !Dual volume is at nodes, and computed here in two different ways for checking. allocate( dual_vol1(nnodes)) allocate( dual_vol2(nnodes)) !Directed area vector is at edge. allocate(directed_area_vector(nedges,3)) !Initialize them. dual_vol1 = zero dual_vol2 = zero directed_area_vector = zero !Compute them via a loop over tetrahedra. tet_loop : do i = 1, ntet !Compute the tetra centroid coordinates. xc = fourth*( x(tet(i,1)) + x(tet(i,2)) + x(tet(i,3)) + x(tet(i,4)) ) yc = fourth*( y(tet(i,1)) + y(tet(i,2)) + y(tet(i,3)) + y(tet(i,4)) ) zc = fourth*( z(tet(i,1)) + z(tet(i,2)) + z(tet(i,3)) + z(tet(i,4)) ) local_edge_loop : do j = 1, 6 ! Each tetra has 6 local edges. !Local edge orientation: k1 -> k2. k1 = loc_edge(j,1) ! Local node number of node1 of edge j. k2 = loc_edge(j,2) ! Local node number of node2 of edge j. !Global edge orientation not necessarily n1 -> n2. It is min(n1,n2)-> max(n1,n2). n1 = tet(i,k1) ! Global node number of node1 of edge j. n2 = tet(i,k2) ! Global node number of node2 of edge j. !Edge midpoint coordinates. xm = half*( x(n1) + x(n2) ) ym = half*( y(n1) + y(n2) ) zm = half*( z(n1) + z(n2) ) !Centroid coordinates of the left face of the edge j. kv1 = loc_face(j,left,1) kv2 = loc_face(j,left,2) kv3 = loc_face(j,left,3) v1 = tet(i,kv1) v2 = tet(i,kv2) v3 = tet(i,kv3) xcl = third*( x(v1) + x(v2) + x(v3) ) ycl = third*( y(v1) + y(v2) + y(v3) ) zcl = third*( z(v1) + z(v2) + z(v3) ) !Centroid coordinates of the right face of the edge j. kv1 = loc_face(j,rght,1) kv2 = loc_face(j,rght,2) kv3 = loc_face(j,rght,3) v1 = tet(i,kv1) v2 = tet(i,kv2) v3 = tet(i,kv3) xcr = third*( x(v1) + x(v2) + x(v3) ) ycr = third*( y(v1) + y(v2) + y(v3) ) zcr = third*( z(v1) + z(v2) + z(v3) ) !------------------------------------------------------------------------------- ! Dual face normal vectors. Here, we compute them on the assumption that the edge ! orientation is n1 -> n2. The sum is the contributon to directed area vector. !------------------------------------------------------------------------------- ! Dual face 1: triangle with (edge midpoint)-(tetra centroid)-(left face centrpoid). normal1 = triangle_area_vector(xm,xc,xcl, ym,yc,ycl, zm,zc,zcl) ! Dual face 2: triangle with (edge midpoint)-(right face centrpoid)-(tetra centroid). normal2 = triangle_area_vector(xm,xcr,xc, ym,ycr,yc, zm,zcr,zc) ! Sum of the two. directed_area_contribution = normal1 + normal2 !------------------------------------------------------------------------------- ! To compute the dual volume by ! dual volume = sum_over_edges [directed area vector]*[edge vector]/6 in 3D. ! Accumulate the contribution to the dual volume. !------------------------------------------------------------------------------- ! Compute the inner product: directed_area_vector * edge_vector. ! This is independent of the edge orientation. So, keep assuming n1 -> n2, ! for consistency with the directed area vector contribution. edge_vector(ix) = x(n2) - x(n1) edge_vector(iy) = y(n2) - y(n1) edge_vector(iz) = z(n2) - z(n1) dxjk_dot_njk = dot_product( directed_area_contribution, edge_vector ) dual_vol1( n1) = dual_vol1( n1) + sixth * dxjk_dot_njk dual_vol1( n2) = dual_vol1( n2) + sixth * dxjk_dot_njk !------------------------------------------------------------------------------- ! To compute the directed-area vector of the edge. ! Add the contribution to the directed-area vector. !------------------------------------------------------------------------------- i_edge = 0 ! Case(1): Edge orientations match: local k1 -> k2, global n1 -> n2. if (n1 < n2) then !Find the global edge number for n1 -> n2. do k = 1, nnghbrs(n1) if (nghbr(n1,k) == n2) i_edge = edge_number(n1,k) end do if (i_edge == 0) stop ! Error !Add the contribution to the directed area vector to 'i_edge'. directed_area_vector(i_edge,:) = directed_area_vector(i_edge,:) + directed_area_contribution ! Case(2): Edge orientations are opposite: local k1 -> k2, global n2 -> n1. else !Find the global edge number for n2 -> n1. do k = 1, nnghbrs(n2) if (nghbr(n2,k) == n1) i_edge = edge_number(n2,k) end do if (i_edge == 0) stop ! Error !Add the contribution to the directed area vector to 'i_edge' with minus sign. directed_area_vector(i_edge,:) = directed_area_vector(i_edge,:) - directed_area_contribution endif end do local_edge_loop !---------------------------------------------------------------------------------- ! Alternative: Compute the dual volume as a sum of 1/4 of the tetrahedral volume. x1 = x(tet(i,1)) x2 = x(tet(i,2)) x3 = x(tet(i,3)) x4 = x(tet(i,4)) y1 = y(tet(i,1)) y2 = y(tet(i,2)) y3 = y(tet(i,3)) y4 = y(tet(i,4)) z1 = z(tet(i,1)) z2 = z(tet(i,2)) z3 = z(tet(i,3)) z4 = z(tet(i,4)) dual_vol2(tet(i,:)) = dual_vol2(tet(i,:)) + fourth*tet_volume(x1,x2,x3,x4, y1,y2,y3,y4, z1,z2,z3,z4) !---------------------------------------------------------------------------------- end do tet_loop ! Edge-based data construction is complete. Now, we verify them below. !------------------------------------------------------------------------------------ !------------------------------------------------------------------------------------ ! Verify the data. ! 1. Check the dual volume sum, which must be equal to the total volume. ! Also, check if the dual volumes computed in diffrent ways match. allocate(diff_dual(nnodes)) total_vol1 = zero total_vol2 = zero do i = 1, nnodes total_vol1 = total_vol1 + dual_vol1(i) total_vol2 = total_vol2 + dual_vol2(i) diff_dual(i) = dual_vol1(i) - dual_vol2(i) end do write(*,*) write(*,*) " sum of the dual_vol1 = ", total_vol1 write(*,*) " sum of the dual_vol2 = ", total_vol2 write(*,*) " Total domain volume = ", one write(*,*) write(*,*) " Max of sum_at_node (dual_vol1-dual_vol2) ", maxval(diff_dual) write(*,*) " Min of sum_at_node (dual_vol1-dual_vol2) ", minval(diff_dual) write(*,*) ! 2. Check the directed area: Sum must be zero at all nodes. !Accumulate the directed area at nodes: all outward from the node. allocate(sum_directed_area_vector(nnodes,3)) sum_directed_area_vector = zero edge_loop : do j = 1, nedges n1 = edge(j,1) n2 = edge(j,2) !directed_area_vector points from n1 to n2. So, add to n1. sum_directed_area_vector(n1,:) = sum_directed_area_vector(n1,:) + directed_area_vector(j,:) !directed_area_vector points from n1 to n2. So, add to n2 with minus sign. sum_directed_area_vector(n2,:) = sum_directed_area_vector(n2,:) - directed_area_vector(j,:) end do edge_loop !Close them with boundary outward normals. boundary_element : do i = 1, nbtria !Three nodes of the triangle i. v1 = btria(i,1) v2 = btria(i,2) v3 = btria(i,3) !Boundary triangle is oriented inward. So, reverse the order to compute the outward vector. normal_b = triangle_area_vector( x(v3),x(v2),x(v1), y(v3),y(v2),y(v1), z(v3),z(v2),z(v1) ) !Dual area for each node is exactly 1/3 of the triangle area. sum_directed_area_vector(v1,:) = sum_directed_area_vector(v1,:) + third*normal_b sum_directed_area_vector(v2,:) = sum_directed_area_vector(v2,:) + third*normal_b sum_directed_area_vector(v3,:) = sum_directed_area_vector(v3,:) + third*normal_b end do boundary_element ! Print the min and max over all nodes in the grid, which must be zero. write(*,*) write(*,*) " Max of sum_at_node(directed area vector(x)) = ", maxval(sum_directed_area_vector(:,ix)) write(*,*) " Max of sum_at_node(directed area vector(y)) = ", maxval(sum_directed_area_vector(:,iy)) write(*,*) " Max of sum_at_node(directed area vector(z)) = ", maxval(sum_directed_area_vector(:,iz)) write(*,*) write(*,*) " Min of sum_at_node(directed area vector(x)) = ", minval(sum_directed_area_vector(:,ix)) write(*,*) " Min of sum_at_node(directed area vector(y)) = ", minval(sum_directed_area_vector(:,iy)) write(*,*) " Min of sum_at_node(directed area vector(z)) = ", minval(sum_directed_area_vector(:,iz)) write(*,*) write(*,*) write(*,*) " ---------------------------------------------------------------------" end subroutine compute_edge_based_data !******************************************************************************* !******************************************************************************* ! 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 !******************************************************************************** end program tetgrid_cube_perturbed