!******************************************************************************* ! Educationally-Designed Unstructured 2D (EDU2D) Code ! ! --- EDU2D-50YEN-Tria-Grids ! ! This is a specialized grid generation code. ! ! ! This code generates various 2D triangular grids: ! ! ---------------------------------------------------------- ! Type Name ! ---------------------------------------------------------- ! 1. Sector : sector ! 2. Disk : disk ! 3. 50-yen : 50yen ! 4. Inverted 50-yen : i50yen ! 5. Scaled-Inverted 50-yen : si50yen ! 6. Perturbed-Scaled-Inverted 50-yen : psi50yen ! 7. Swapped-Perturbed-Scaled-Inverted 50-yen : spsi50yen ! ---------------------------------------------------------- ! ! For each case, it generates the following data files: ! ! - Grid file for EDU2D codes (.grid) ! - Tecplot file for viewing (.dat) ! - Sample boundary condition file. (.bcmap) ! - A file containing structured indices. (.k) ! - VTK file. (.vtk) ! - SU2 grid file. (.su2) ! ! ! This is Version 3 (October 2021). ! ! - 10-01-2032: Fixed a bug about nodes2 array. It didn't work with gfortran; it does now. ! ! - 07-23-2019: Added .su2 and .vkt files. ! ! - A bug fixed in the swapping part. Skip if the two elements shearing an edge ! involve empty neighbors. <- Version 0 (July 2015). ! ! ! Note: The final grid, e.g., si50y-grid, may be useful for CFD simulations over a cylinder ! because it has much less nodes on the outer boundary than on the inner cylinder, ! i.e., saving lots of computational time by having less nodes in the outer field ! while keeping isotropic node distribution (unlike grids generated regularly in ! the polar coordinate system -> the same number of nodes on inner and outer boundaries). ! Example: 384 nodes around the inner body (circle) and 62 nodes in the radial direction ! towards the outer boundary (circle) gives 384x62=23808 nodes in total in the case ! of a regular grid in the polar coordinates; while this code results in a grid with ! 12474 nodes (with 384 nodes in the inner body, 12 nodes on the outer boundary, and ! 62 nodes in the radial direction). ! ! ! Note: Structured indices (k1,k2,k3) may be used to regularly coarsen the grid for multigrid. ! Coarsening is not performed in this progarm. ! ! Note: You can Google image-search "50 yen". ! ! ! The program may be updated in future. ! ! ------------------------------------------------------------------------------ ! ! Input: ! ! smoothing_type !(!) Specify a type of grid smoothing (0, 1, or 2) ! nr_gs !(2) Number of division of the generating sector. ! nr_gs_in !(3) Specify the location of a hole to create a 50-yen domain. ! sf !(4) Stretching factor ! spacing_min !(5) Minimum grid spacing desired. ! ! Output: ! ! (1) sector - sector.grid : Grid file for use in EDU2D codes ! sector.dat : Tecplot file for viewing ! sector.bcmap : Boundary condition file (sample) ! sector.k : Structured indices ! ! (2) disk - disk.grid : Grid file for use in EDU2D codes ! disk.dat : Tecplot file for viewing ! disk.bcmap : Boundary condition file (sample) ! disk.k : Structured indices ! ! (3) 50yen - 50yen.grid : Grid file for use in EDU2D codes ! 50yen.dat : Tecplot file for viewing ! 50yen.bcmap : Boundary condition file (sample) ! 50yen.k : Structured indices ! ! (4) i50yen - i50yen.grid : Grid file for use in EDU2D codes ! i50yen.dat : Tecplot file for viewing ! i50yen.bcmap : Boundary condition file (sample) ! i50yen.k : Structured indices ! ! (5) si50yen - si50yen.grid : Grid file for use in EDU2D codes ! si50yen.dat : Tecplot file for viewing ! si50yen.bcmap : Boundary condition file (sample) ! si50yen.k : Structured indices ! ! (6) psi50yen - psi50yen.grid : Grid file for use in EDU2D codes ! psi50yen.dat : Tecplot file for viewing ! psi50yen.bcmap : Boundary condition file (sample) ! psi50yen.k : Structured indices ! ! (7) spsi50yen - spsi50yen.grid : Grid file for use in EDU2D codes ! spsi50yen.dat : Tecplot file for viewing ! spsi50yen.bcmap : Boundary condition file (sample) ! spsi50yen.k : Structured indices ! ! ! Note: Grid (7) can be highly and really really irregular. ! I can imagine many CFD codes would fail on this grid. ! ! Note: I have used Grid (5), si50yen, for various test cases, e.g., ! a laminar flow over a cylinder, and I like it. You may also find it ! useful for a similar purpose or some other purposes. ! ! Note: The resulting grids are good grids (plese don't blame them; even Grid (7) ). ! If your CFD code doesn't work with it, it means that the numerical ! scheme and the solver used are not robust and so they need to be improved. ! ! ! ! ! written by Dr. Katate Masatsuka (info[at]cfdbooks.com), ! ! the author of useful CFD books, "I do like CFD" (http://www.cfdbooks.com). ! ! This F90 code is written and made available for an educational purpose. ! This file may be updated in future. ! ! Katate Masatsuka, http://www.cfdbooks.com !******************************************************************************* program edu2d_50yen_tria_grids implicit none !--------------------------- ! Parameters integer , parameter :: dp = selected_real_kind(P=15) real(dp), parameter :: zero = 0.0_dp real(dp), parameter :: one = 1.0_dp real(dp), parameter :: two = 2.0_dp real(dp), parameter :: three = 3.0_dp real(dp), parameter :: half = 0.5_dp real(dp), parameter :: pi = 3.14159265358979323846_dp !--------------------------- ! Custom data types type tria_data integer, dimension(3) :: v !Vertices (nodes) of the triangle integer, dimension(3) :: e !Element neighbors integer :: type !Type of triangle (upward, downward, cylinder) end type tria_data type node_data_xy integer :: gnode !Global node number real(dp) :: x, y !Coordinates in xy-plane. integer, dimension(6) :: elm !Neighbor elements (triangles) integer :: nelms !Actual number of neighbor elements, which !is less than 6 for the grid considered here. end type node_data_xy type node_data real(dp) :: x, y !Coordinates in xy space end type node_data type bnode_data integer, dimension(:), pointer :: node integer :: nnodes end type bnode_data !--------------------------- ! Input variables integer :: smoothing_type !(!) Specify a type of grid smoothing (0, 1, or 2) integer :: nr_gs !(2) Number of division of the generating sector. integer :: nr_gs_in !(3) Specify the location of a hole to create a 50-yen domain. real(dp) :: sf !(4) Stretching factor !--------------------------- ! Output File names character(80) :: filename_sector = "sector.dat" character(80) :: filename_sector_k = "sector.k" character(80) :: filename_sector_grid = "sector.grid" character(80) :: filename_sector_bcmap = "sector.bcmap" character(80) :: filename_sector_su2 = "sector.su2" character(80) :: filename_sector_vtk = "sector.vtk" character(80) :: filename_disk = "disk.dat" character(80) :: filename_disk_k = "disk.k" character(80) :: filename_disk_grid = "disk.grid" character(80) :: filename_disk_bcmap = "disk.bcmap" character(80) :: filename_disk_su2 = "disk.su2" character(80) :: filename_disk_vtk = "disk.vtk" character(80) :: filename_50yen = "50yen.dat" character(80) :: filename_50yen_k = "50yen.k" character(80) :: filename_50yen_grid = "50yen.grid" character(80) :: filename_50yen_bcmap = "50yen.bcmap" character(80) :: filename_50yen_su2 = "50yen.su2" character(80) :: filename_50yen_vtk = "50yen.vtk" character(80) :: filename_i50yen = "i50yen.dat" character(80) :: filename_i50yen_k = "i50yen.k" character(80) :: filename_i50yen_grid = "i50yen.grid" character(80) :: filename_i50yen_bcmap = "i50yen.bcmap" character(80) :: filename_i50yen_su2 = "i50yen.su2" character(80) :: filename_i50yen_vtk = "i50yen.vtk" character(80) :: filename_si50yen = "si50yen.dat" character(80) :: filename_si50yen_k = "si50yen.k" character(80) :: filename_si50yen_grid = "si50yen.grid" character(80) :: filename_si50yen_bcmap = "si50yen.bcmap" character(80) :: filename_si50yen_su2 = "si50yen.su2" character(80) :: filename_si50yen_vtk = "si50yen.vtk" character(80) :: filename_psi50yen = "psi50yen.dat" character(80) :: filename_psi50yen_k = "psi50yen.k" character(80) :: filename_psi50yen_grid = "psi50yen.grid" character(80) :: filename_psi50yen_bcmap = "psi50yen.bcmap" character(80) :: filename_psi50yen_su2 = "psi50yen.su2" character(80) :: filename_psi50yen_vtk = "psi50yen.vtk" character(80) :: filename_spsi50yen = "spsi50yen.dat" character(80) :: filename_spsi50yen_k = "spsi50yen.k" character(80) :: filename_spsi50yen_grid = "spsi50yen.grid" character(80) :: filename_spsi50yen_bcmap = "spsi50yen.bcmap" character(80) :: filename_spsi50yen_su2 = "spsi50yen.su2" character(80) :: filename_spsi50yen_vtk = "spsi50yen.vtk" !--------------------------- ! Local variables integer :: os !IO constant integer :: i,j,k,inode integer :: ntrias_gs !Number of triangles in the generating sector integer :: nnodes_gs !Number of nodes in the generating sector real(dp) :: r_gs !Raidus of the generating sector real(dp) :: dr_gs !Spacing along the side of the generating sector real(dp) :: r2, dtheta, theta type(node_data_xy), dimension(:), pointer :: nodes1, nodes2, node_gs integer , dimension(:), pointer :: k1_gs, k2_gs, k3_gs, kr_gs type(tria_data) , dimension(:), allocatable :: tria_gs type(bnode_data) , dimension(:), allocatable :: bnode_gs integer , dimension(:), allocatable :: node_map integer :: nnodes_disk, ntrias_disk type(node_data_xy), dimension(:), pointer :: node_disk integer , dimension(:), pointer :: k1_disk integer , dimension(:), pointer :: k2_disk integer , dimension(:), pointer :: k3_disk integer , dimension(:), pointer :: kr_disk type(tria_data) , dimension(:), allocatable :: tria_disk type(bnode_data) , dimension(:), allocatable :: bnode_disk integer :: nnodes_50yen, ntrias_50yen type(node_data_xy), dimension(:), pointer :: node_50yen integer , dimension(:), pointer :: k1_50yen integer , dimension(:), pointer :: k2_50yen integer , dimension(:), pointer :: k3_50yen integer , dimension(:), pointer :: kr_50yen type(tria_data) , dimension(:), allocatable :: tria_50yen type(bnode_data) , dimension(:), allocatable :: bnode_50yen integer , dimension(:), allocatable :: f2c integer , dimension(:), allocatable :: f2c_tria integer :: nnodes_i50yen, ntrias_i50yen type(node_data_xy), dimension(:), pointer :: node_i50yen integer , dimension(:), pointer :: k1_i50yen integer , dimension(:), pointer :: k2_i50yen integer , dimension(:), pointer :: k3_i50yen integer , dimension(:), pointer :: kr_i50yen type(tria_data) , dimension(:), allocatable :: tria_i50yen type(bnode_data) , dimension(:), allocatable :: bnode_i50yen integer :: i_smoothing, max_smoothing, v(3) real(dp) :: smoothing_factor, dxy_norm real(dp) , dimension(:,:), allocatable :: dxy real(dp) :: r, s, rmin, rmax, w, cp, spacing_min, spacing real(dp) :: rin, rout, nx, ny, area integer :: mglevels integer :: v1, v2, v3, v4, negative_volumes, n_attempts, n_attempts_max logical :: all_dirichlet, cp_min_reached integer, dimension(:,:), allocatable :: nghbr integer, dimension(:) , allocatable :: nnghbrs real(dp) :: dx, dy, rn integer :: ne, kn, e1, e2, e3, e4, ke1, ke2, ke3, ke4, n_swapped integer, dimension(3,2) :: km real(dp) :: area124, area234 character(80), dimension(3) :: bnames all_dirichlet = .true. n_attempts_max = 15 !******************************************************************************* ! 0. Selecet the type of grid smoothing !******************************************************************************* 10 write(*,*) write(*,*) " Input - Smoothing on the disk triangulation (e.g., Try 0)" write(*,*) " 0 : No smoothing" write(*,*) " 1 : Smoothing" write(*,*) " 2 : Constrained smoothing (some grid lines are fixed)" write(*,*) write(*,*) " smoothing_type = ?" read(*,*) smoothing_type if ( smoothing_type/=0 .and. & smoothing_type/=1 .and. & smoothing_type/=2 ) then write(*,*) " Invalid input : ", smoothing_type go to 10 endif write(*,*) write(*,*) " smoothing_type = ", smoothing_type write(*,*) !******************************************************************************* ! 1. Systematic triangulation of the generating sector (isotropic). ! ! This is a sector with the central angle 60 degrees. ! It is called here the generating sector. ! It is located in the xy-plane. ! We first triangulate this, and use it to build a triangulation of the ! whole circle (the disk) later. ! ! ____________ ! \/\/\/\/\/\/ This is an example corresponding to nr_gs = 6. ! \/\/\/\/\/ ! \/\/\/\/ y ^ ! \/\/\/ | ! \/\/ | ! \/ ----> x ! ! NOTE: The number of triangules = 1 + 3 + 5 + ... + 2*nr_gs-1 = nr_gs*nr_gs ! The number of ndoes = 1 + 2 + 3 + ... + nr_gs+1 ! ! NOTE: The length of a side is set to be 0.5 inside the program: r_gs = half. ! !******************************************************************************* write(*,*) "*****************************************************************" write(*,*) " 1. Generate a triangulated sector" write(*,*) "*****************************************************************" write(*,*) "Input - Division of the generating sector: Positive integer. (e.g., Try 64) " write(*,*) " nr_gs = ?" read(*,*) nr_gs k = 0 do k = k + 1 if (mod(nr_gs,2**(k-1))==0) then mglevels = k else write(*,*) write(*,*) " Maximum multigrid level is ", mglevels, " just FYI..." write(*,*) exit endif end do r_gs = half ! Radius of the isotropic triangle dr_gs = r_gs/real(nr_gs,dp) ! Uniform refinement nnodes_gs = (nr_gs+1)*((nr_gs+1)+1)/2 ntrias_gs = nr_gs**2 allocate(node_gs(nnodes_gs)) allocate(tria_gs(ntrias_gs)) allocate(k1_gs(nnodes_gs)) allocate(k2_gs(nnodes_gs)) allocate(k3_gs(nnodes_gs)) allocate(kr_gs(nnodes_gs)) allocate(bnode_gs(3)) ! 3 boundaries allocate(bnode_gs(1)%node(nr_gs+1)) allocate(bnode_gs(2)%node(nr_gs+1)) allocate(bnode_gs(3)%node(nr_gs+1)) nnodes_gs = 0 ntrias_gs = 0 triangulate : do i = 1, nr_gs ! r1 = dr_gs * real(i-1,dp) ! Current arc r2 = dr_gs * real(i,dp) ! Next arc ! Nodes on the current arc (r = r1) call my_alloc_ndxy_ptr(nodes1,i) if (i == 1) then nodes1(1)%x = zero nodes1(1)%y = zero nodes1(1)%gnode = 1 nnodes_gs = 1 node_gs(1)%x = zero node_gs(1)%y = zero k1_gs(1) = 0 k2_gs(1) = 0 k3_gs(1) = - ( k1_gs(1) + k2_gs(1) ) else do k = 1, i nodes1(k)%x = nodes2(k)%x nodes1(k)%y = nodes2(k)%y nodes1(k)%gnode = nodes2(k)%gnode end do endif allocate(nodes2(1)) !This is necessary to initialize nodes2. ! Nodes on the next arc (r = r2): New nodes call my_alloc_ndxy_ptr(nodes2,i+1) dtheta = (pi/three) / real(i,dp) do k = 1, i+1 theta = dtheta * real(k-1,dp) nodes2(k)%x = r2 * cos(theta) nodes2(k)%y = r2 * sin(theta) nnodes_gs = nnodes_gs + 1 nodes2(k)%gnode = nnodes_gs node_gs(nnodes_gs)%x = nodes2(k)%x node_gs(nnodes_gs)%y = nodes2(k)%y k1_gs(nnodes_gs) = i + (1 - k) k2_gs(nnodes_gs) = k - 1 k3_gs(nnodes_gs) = - ( k1_gs(nnodes_gs) + k2_gs(nnodes_gs) ) end do ! Triangulate the region between nodes1 and nodes2. ! Right to left ! Note: There are Type 1 and 2, but in this code, it doesn't matter. ! (It'd be good to distinguish them if we wish to build a 3D grid ! from the triangulation.) ! Type 1 triangle ! ! nodes2(k+1) nodes2(k) ! 1 2 ! o-----------o ! \ / ! \ / ! \ / ! \ / ! \ / ! o ! 3 ! nodes1(k) downward_tria : do k = 1, i ntrias_gs = ntrias_gs + 1 tria_gs(ntrias_gs)%v(1) = nodes2(k+1)%gnode tria_gs(ntrias_gs)%v(2) = nodes2(k )%gnode tria_gs(ntrias_gs)%v(3) = nodes1(k )%gnode tria_gs(ntrias_gs)%type = 1 end do downward_tria ! Type 2 triangle ! ! nodes2(k+1) ! 3 ! o ! / \ ! / \ ! / \ ! / \ ! / \ ! o-----------o ! 2 1 ! nodes1(k+1) nodes1(k) if (i > 1) then upward_tria : do k = 1, i-1 ntrias_gs = ntrias_gs + 1 tria_gs(ntrias_gs)%v(1) = nodes1(k )%gnode tria_gs(ntrias_gs)%v(2) = nodes1(k+1)%gnode tria_gs(ntrias_gs)%v(3) = nodes2(k+1)%gnode tria_gs(ntrias_gs)%type = 2 end do upward_tria endif end do triangulate ! In fact, the ordering needs to be reversed.. do i = 1, ntrias_gs v1 = tria_gs(i)%v(1) v2 = tria_gs(i)%v(2) v3 = tria_gs(i)%v(3) tria_gs(i)%v(1) = v3 tria_gs(i)%v(2) = v2 tria_gs(i)%v(3) = v1 end do ! So, now I have a triangulation defined by tria_gs and node_gs. ! Number of triangles = ntrias_gs: tria_gs(:)%v(1:3) ! Number of nodes = nnodes_gs: node_gs(:)%x, node_gs(:)%y ! ! Construct a boundary data. bnode_gs(1)%nnodes = 0 do i = nnodes_gs, 1, -1 !Left boundary if ( k1_gs(i) == 0 ) then bnode_gs(1)%nnodes = bnode_gs(1)%nnodes + 1 bnode_gs(1)%node(bnode_gs(1)%nnodes) = i endif end do bnode_gs(2)%nnodes = 0 bnode_gs(3)%nnodes = 0 do i = 1, nnodes_gs !Right boundary if ( k2_gs(i) == 0 ) then bnode_gs(2)%nnodes = bnode_gs(2)%nnodes + 1 bnode_gs(2)%node(bnode_gs(2)%nnodes) = i endif !Top boundary (circumferential) if ( k3_gs(i) == -nr_gs ) then bnode_gs(3)%nnodes = bnode_gs(3)%nnodes + 1 bnode_gs(3)%node(bnode_gs(3)%nnodes) = i endif end do call write_k_file( filename_sector_k, k1_gs, k2_gs, k3_gs, nnodes_gs) call k_radial(kr_gs, k1_gs,k2_gs,k3_gs,nnodes_gs) call write_tecplot_file(filename_sector, node_gs, k1_gs, k2_gs, k3_gs, kr_gs, & nnodes_gs, tria_gs, ntrias_gs ) call write_bcmap_file( filename_sector_bcmap, 3, all_dirichlet) call write_grid_file( filename_sector_grid,node_gs,nnodes_gs,tria_gs,ntrias_gs, 3,bnode_gs) call write_vtk_file(filename_sector_vtk, nnodes_gs,node_gs,ntrias_gs,tria_gs ) bnames(1) = "LEFT" bnames(2) = "RIGHT" bnames(3) = "TOP" call write_su2_file(filename_sector_su2, nnodes_gs,node_gs,ntrias_gs,tria_gs, 3,bnode_gs,bnames ) !******************************************************************************* ! 2. Generate a triangulated disk. ! ! Now, rotate and copy the sector triangulation onto 5 places to form ! a triangulation of a whole disk (6 patches in total). ! Generate new data: tria_disk and node_disk ! Number of triangles = ntrias_disk ! Number of nodes = nnodes_disk !******************************************************************************* write(*,*) "*****************************************************************" write(*,*) " 2. Generate a triangulated disk." write(*,*) "*****************************************************************" write(*,*) nnodes_disk = nnodes_gs*6 !<- More than enough (by the # of overlapping ndoes) ntrias_disk = ntrias_gs*6 allocate(node_disk(nnodes_disk)) allocate(tria_disk(ntrias_disk)) allocate(k1_disk(nnodes_disk)) allocate(k2_disk(nnodes_disk)) allocate(k3_disk(nnodes_disk)) allocate(kr_disk(nnodes_disk)) allocate(bnode_disk(1)) ! 1 boundary allocate(bnode_disk(1)%node(6*nr_gs+1)) ! nnodes_disk = 0 ntrias_disk = 0 ! Copy the data from the generating-sector data to the disk data. ! Copy the node data ! This is the first one: theta = 0 to 60 degrees do i = 1, nnodes_gs node_disk(i)%x = node_gs(i)%x node_disk(i)%y = node_gs(i)%y k1_disk(i) = k1_gs(i) k2_disk(i) = k2_gs(i) k3_disk(i) = k3_gs(i) end do nnodes_disk = nnodes_gs ! Copy the triangle data ! This is the first one: theta = 0 to 60 degrees do i = 1, ntrias_gs tria_disk(i)%v = tria_gs(i)%v tria_disk(i)%type = tria_gs(i)%type end do ntrias_disk = ntrias_gs ! Now generate other parts of the disk: i=1,2,3,4,5 ! 1. theta = 60 to 120 degrees ! 2. theta = 120 to 180 degrees ! 3. theta = 180 to 240 degrees ! 4. theta = 240 to 300 degrees ! 5. theta = 300 to 360 degrees ! Each part has (n+1)*(n+2)/2 nodes: 1+2+3+...+(nr_gs+1) = (n_gs+1)*(n_gs+2)/2 allocate(node_map((nr_gs + 1)*(nr_gs + 2)/2)) new_sectors : do i = 1, 5 ! (1)Generate new nodes and construct a node map (one to one) ! inode = local node numbering for node_map(:) node_map(1) = 1 !Node at the origin do k = 2, nr_gs + 1 !Origin to outer boundary do j = 1, k !Right to left inode = (k-1)*k/2 + j !Local node number = Right-most node + increment(1,2,..,k) ! Right side nodes are existing nodes: Left side nodes of the previous sector if (j==1) then if (i == 1) then node_map(inode) = (k-1)*k/2 + k !Left-most node of the original sector elseif (i == 2) then node_map(inode) = nnodes_gs + (k-1)*k/2 !Left-most node of the second sector else node_map(inode) = nnodes_gs + (nnodes_gs-(nr_gs+1))*(i-2) + (k-1)*k/2 endif ! Left side of the last one is the right side of the original sector elseif (i==5 .and. j==k) then node_map(inode) = (k-1)*k/2 + 1 !Right-most node of the original sector ! New nodes: Rotate the original nodes by theta = i*pi/3 (i times 60 degrees). else theta = (pi/three) * real(i,dp) nnodes_disk = nnodes_disk + 1 node_map(inode) = nnodes_disk node_disk(nnodes_disk)%x = cos(theta)*node_gs(inode)%x - sin(theta)*node_gs(inode)%y node_disk(nnodes_disk)%y = sin(theta)*node_gs(inode)%x + cos(theta)*node_gs(inode)%y if (i==1) then k1_disk(nnodes_disk) = -k2_gs(inode) k2_disk(nnodes_disk) = -k3_gs(inode) k3_disk(nnodes_disk) = -k1_gs(inode) elseif(i==2) then k1_disk(nnodes_disk) = k3_gs(inode) k2_disk(nnodes_disk) = k1_gs(inode) k3_disk(nnodes_disk) = k2_gs(inode) elseif(i==3) then k1_disk(nnodes_disk) = -k1_gs(inode) k2_disk(nnodes_disk) = -k2_gs(inode) k3_disk(nnodes_disk) = -k3_gs(inode) elseif(i==4) then k1_disk(nnodes_disk) = k2_gs(inode) k2_disk(nnodes_disk) = k3_gs(inode) k3_disk(nnodes_disk) = k1_gs(inode) elseif(i==5) then k1_disk(nnodes_disk) = -k3_gs(inode) k2_disk(nnodes_disk) = -k1_gs(inode) k3_disk(nnodes_disk) = -k2_gs(inode) endif endif end do end do ! (2)Generate triangles on the new sector. do k = 1, ntrias_gs ntrias_disk = ntrias_disk + 1 tria_disk(ntrias_disk)%v = node_map(tria_gs(k)%v) tria_disk(ntrias_disk)%type = tria_gs(k)%type end do end do new_sectors !-------------------------------------------------------------------------------- ! (If requested) Apply smoothing on the disk triangulation if (smoothing_type /= 0) then write(*,*) "-----------------------------------------------------" write(*,*) " Applying smoothing..." write(*,*) "-----------------------------------------------------" allocate(dxy(nnodes_disk,2)) smoothing_factor = 0.05_dp max_smoothing = 1000 smoothing : do i_smoothing = 1, max_smoothing dxy = zero !Initialize the changes ! Accumulate the changes by looping over triangles do i = 1, ntrias_disk v = tria_disk(i)%v ! x-coordinate dxy(v(1),1) = dxy(v(1),1) + ( node_disk(v(2))%x - node_disk(v(1))%x ) dxy(v(1),1) = dxy(v(1),1) + ( node_disk(v(3))%x - node_disk(v(1))%x ) dxy(v(2),1) = dxy(v(2),1) + ( node_disk(v(1))%x - node_disk(v(2))%x ) dxy(v(2),1) = dxy(v(2),1) + ( node_disk(v(3))%x - node_disk(v(2))%x ) dxy(v(3),1) = dxy(v(3),1) + ( node_disk(v(1))%x - node_disk(v(3))%x ) dxy(v(3),1) = dxy(v(3),1) + ( node_disk(v(2))%x - node_disk(v(3))%x ) ! y-coordinate dxy(v(1),2) = dxy(v(1),2) + ( node_disk(v(2))%y - node_disk(v(1))%y ) dxy(v(1),2) = dxy(v(1),2) + ( node_disk(v(3))%y - node_disk(v(1))%y ) dxy(v(2),2) = dxy(v(2),2) + ( node_disk(v(1))%y - node_disk(v(2))%y ) dxy(v(2),2) = dxy(v(2),2) + ( node_disk(v(3))%y - node_disk(v(2))%y ) dxy(v(3),2) = dxy(v(3),2) + ( node_disk(v(1))%y - node_disk(v(3))%y ) dxy(v(3),2) = dxy(v(3),2) + ( node_disk(v(2))%y - node_disk(v(3))%y ) end do ! Make changes to each node except the boundary nodes. dxy_norm = -1000000.0_dp do i=1, nnodes_disk if ( abs( sqrt(node_disk(i)%x**2 + node_disk(i)%y**2) - r_gs ) < 1.0e-14 ) then else ! Constrained smoothing: skip nodes in the sector boundaries. if (smoothing_type == 2) then if ( abs(atan2(node_disk(i)%y, node_disk(i)%x) - pi/three) < 1.0e-13 ) cycle if ( abs(atan2(node_disk(i)%y, node_disk(i)%x) + pi/three) < 1.0e-13 ) cycle if ( abs(atan2(node_disk(i)%y, node_disk(i)%x) - two*pi/three) < 1.0e-13 ) cycle if ( abs(atan2(node_disk(i)%y, node_disk(i)%x) + two*pi/three) < 1.0e-13 ) cycle if ( abs(node_disk(i)%y) < 1.0e-13 ) cycle endif node_disk(i)%x = node_disk(i)%x + smoothing_factor * dxy(i,1) node_disk(i)%y = node_disk(i)%y + smoothing_factor * dxy(i,2) ! L_inf norm of changes scaled by the typical mesh spacing, dr_gs. dxy_norm = max( dxy_norm, abs(dxy(i,1)**2 + dxy(i,2)**2)/dr_gs ) endif end do ! Exit if converged if ( dxy_norm < 1.0e-04) then write(*,*) " Smoothing converged at ", i_smoothing exit smoothing elseif (i_smoothing == max_smoothing) then write(*,*) " Smoothing didn't converge... ", " dxy_norm = ", dxy_norm endif end do smoothing deallocate(dxy) write(*,*) "-----------------------------------------------------" endif !-------------------------------------------------------------------------------- ! Now, at this point, we have a triangulation of a disk defined by ! tria_disk and node_disk. ! Number of triangles = ntrias_disk ! Number of nodes = nnodes_disk call k_radial(kr_disk, k1_disk,k2_disk,k3_disk,nnodes_disk) ! Construct a boundary data. bnode_disk(1)%nnodes = 0 do i = 1, nnodes_disk if ( kr_disk(i) == nr_gs ) then bnode_disk(1)%nnodes = bnode_disk(1)%nnodes + 1 bnode_disk(1)%node(bnode_disk(1)%nnodes) = i endif end do !Add the first one at the end for a closed boundary. bnode_disk(1)%nnodes = bnode_disk(1)%nnodes + 1 bnode_disk(1)%node(bnode_disk(1)%nnodes) = bnode_disk(1)%node(1) call write_k_file( filename_disk_k, k1_disk, k2_disk, k3_disk, nnodes_disk) call write_tecplot_file(filename_disk, node_disk, k1_disk, k2_disk, k3_disk, kr_disk, & nnodes_disk, tria_disk, ntrias_disk ) call write_bcmap_file( filename_disk_bcmap, 1, all_dirichlet) call write_grid_file( filename_disk_grid, node_disk, nnodes_disk, tria_disk, & ntrias_disk, 1,bnode_disk ) call write_vtk_file(filename_disk_vtk, nnodes_disk,node_disk,ntrias_disk,tria_disk ) bnames(1) = "OUTER" call write_su2_file(filename_disk_su2, nnodes_disk,node_disk,ntrias_disk,tria_disk, 1,bnode_disk,bnames ) if (nr_gs == 1) then write(*,*) " cannot generate 50-yen grids... Stop" write(*,*) " [ nr_gs must be greater than 1 to go further. ]" stop endif deallocate(k1_gs,k2_gs,k3_gs,kr_gs,node_gs,tria_gs,bnode_gs) !******************************************************************************* ! 3. Generate a triangulated 50-yen grid. !******************************************************************************* write(*,*) "*****************************************************************" write(*,*) " 3. Generate a triangulated 50-yen grid." write(*,*) "*****************************************************************" write(*,*) write(*,*) " Input - Inner boundary location index (e.g., Try 2) : " write(*,*) " [ Integer greater than 1 and smaller than ", nr_gs, "]" write(*,*) " [ It leads to 6*(input value) nodes along the inner boundary. ]" write(*,*) " nr_gs_in = ? " read(*,*) nr_gs_in if (nr_gs_in >= nr_gs .or. nr_gs_in < 0) then write(*,*) " Invalid input: ", nr_gs_in, " Stop..." stop endif nnodes_50yen = nnodes_disk !<- More than enough (by the # of overlapping ndoes) ntrias_50yen = ntrias_disk allocate(node_50yen(nnodes_50yen)) allocate(tria_50yen(ntrias_50yen)) allocate(k1_50yen(nnodes_50yen)) allocate(k2_50yen(nnodes_50yen)) allocate(k3_50yen(nnodes_50yen)) allocate(kr_50yen(nnodes_50yen)) allocate( f2c(nnodes_disk)) allocate( f2c_tria(ntrias_disk)) f2c = 0 f2c_tria = 0 allocate(bnode_50yen(2)) ! 2 boundaries: allocate(bnode_50yen(1)%node(6*nr_gs+1)) ! <- Inner circle allocate(bnode_50yen(2)%node(6*nr_gs+1)) ! <- outer circle nnodes_50yen = 0 do k = nr_gs, nr_gs_in, -1 do i = 1, nnodes_disk if ( kr_disk(i) == k ) then nnodes_50yen = nnodes_50yen + 1 f2c(i) = nnodes_50yen endif end do end do do i = 1, nnodes_disk if (f2c(i) > 0) then node_50yen(f2c(i))%x = node_disk(i)%x node_50yen(f2c(i))%y = node_disk(i)%y k1_50yen(f2c(i)) = k1_disk(i) k2_50yen(f2c(i)) = k2_disk(i) k3_50yen(f2c(i)) = k3_disk(i) kr_50yen(f2c(i)) = kr_disk(i) end if end do ntrias_50yen = 0 do i = 1, ntrias_disk if ( f2c( tria_disk(i)%v(1) ) > 0 .and. & f2c( tria_disk(i)%v(2) ) > 0 .and. & f2c( tria_disk(i)%v(3) ) > 0 ) then ntrias_50yen = ntrias_50yen + 1 f2c_tria(i) = ntrias_50yen endif end do do i = 1, ntrias_disk if (f2c_tria(i) > 0) then tria_50yen(f2c_tria(i))%v = f2c(tria_disk(i)%v) endif end do ! Construct boundary data. bnode_50yen(1)%nnodes = 0 do i = 1, 6*nr_gs bnode_50yen(1)%nnodes = bnode_50yen(1)%nnodes + 1 bnode_50yen(1)%node(bnode_50yen(1)%nnodes) = i end do bnode_50yen(1)%nnodes = bnode_50yen(1)%nnodes + 1 bnode_50yen(1)%node(bnode_50yen(1)%nnodes) = bnode_50yen(1)%node(1) bnode_50yen(2)%nnodes = 0 do i = nnodes_50yen, nnodes_50yen-(6*nr_gs), -1 if (kr_50yen(i) == nr_gs_in) then bnode_50yen(2)%nnodes = bnode_50yen(2)%nnodes + 1 bnode_50yen(2)%node(bnode_50yen(2)%nnodes) = i endif end do bnode_50yen(2)%nnodes = bnode_50yen(2)%nnodes + 1 bnode_50yen(2)%node(bnode_50yen(2)%nnodes) = bnode_50yen(2)%node(1) call write_k_file( filename_50yen_k, k1_50yen, k2_50yen, k3_50yen, nnodes_50yen) call write_tecplot_file(filename_50yen, node_50yen, k1_50yen, k2_50yen, k3_50yen, & kr_50yen, nnodes_50yen, tria_50yen, ntrias_50yen ) call write_bcmap_file(filename_50yen_bcmap, 2, .false.) call write_grid_file(filename_50yen_grid, node_50yen, nnodes_50yen, tria_50yen, & ntrias_50yen, 2,bnode_50yen ) call write_vtk_file(filename_50yen_vtk, nnodes_50yen,node_50yen,ntrias_50yen,tria_50yen) bnames(1) = "INNER" bnames(2) = "OUTER" call write_su2_file(filename_50yen_su2, nnodes_50yen,node_50yen,ntrias_50yen, & tria_50yen, 2,bnode_50yen,bnames ) deallocate(k1_disk,k2_disk,k3_disk,kr_disk,node_disk,tria_disk,bnode_disk) !******************************************************************************* ! 4. Generate an inverted 50-yen grid. ! NOTE: This grid may have negative-volume triangles near the outer boundary. ! Adjust (increase) the stretching factor 'sf' to avoid it. !******************************************************************************* write(*,*) "*****************************************************************" write(*,*) " 4. Generate an inverted 50-yen grid." write(*,*) "*****************************************************************" write(*,*) write(*,*) " Input - Stretching factor (>1.0; e.g., Try 5.0):" write(*,*) " sf = ?" read(*,*) sf nnodes_i50yen = nnodes_50yen !<- More than enough (by the # of overlapping nodes) ntrias_i50yen = ntrias_50yen allocate(node_i50yen(nnodes_i50yen)) allocate(tria_i50yen(ntrias_i50yen)) allocate(k1_i50yen(nnodes_i50yen)) allocate(k2_i50yen(nnodes_i50yen)) allocate(k3_i50yen(nnodes_i50yen)) allocate(kr_i50yen(nnodes_i50yen)) ! 3 boundary segments: Inner circle, left and right portions of the outer circle. ! Note: Here, we split the outer circle into two parts: Left and right in order ! to be able to treat the left as inflow and the right as outflow in a CFD solver. allocate(bnode_i50yen(3)) allocate(bnode_i50yen(1)%node(6*nr_gs+1)) ! Inner circle allocate(bnode_i50yen(2)%node(6*nr_gs+1)) ! Inflow (Left) allocate(bnode_i50yen(3)%node(6*nr_gs+1)) ! Outflow (Right) !Find out the minimum and maximum distances of the nodes from the origin. rmin = 1.0e+12_dp rmax = -1.0_dp do i = 1, nnodes_50yen r = sqrt( node_50yen(i)%x**2 + node_50yen(i)%y**2 ) rmin = min(r, rmin) rmax = max(r, rmax) end do cp = 0.07_dp !Fixed value here. n_attempts = 0 200 continue do i = 1, nnodes_50yen r = sqrt( node_50yen(i)%x**2 + node_50yen(i)%y**2 ) nx = node_50yen(i)%x/r ny = node_50yen(i)%y/r s = (r-rmax)/(rmin-rmax) !Stretching introduced here. w = ( cp*s + (one-exp(sf*s))/(one-exp(sf)) )/(cp+one) r = ( w*rmax + (one-w)*rmin ) node_i50yen(i)%x = r*nx node_i50yen(i)%y = r*ny k1_i50yen(i) = k1_50yen(i) k2_i50yen(i) = k2_50yen(i) k3_i50yen(i) = k3_50yen(i) kr_i50yen(i) = kr_50yen(i) end do negative_volumes = 0 ntrias_i50yen = ntrias_50yen do i = 1, ntrias_50yen tria_i50yen(i)%v(3) = tria_50yen(i)%v(1) tria_i50yen(i)%v(2) = tria_50yen(i)%v(2) tria_i50yen(i)%v(1) = tria_50yen(i)%v(3) area = 0.5_dp*( node_i50yen(tria_i50yen(i)%v(1))%x * & (node_i50yen(tria_i50yen(i)%v(2))%y-node_i50yen(tria_i50yen(i)%v(3))%y) & + node_i50yen(tria_i50yen(i)%v(2))%x * & (node_i50yen(tria_i50yen(i)%v(3))%y-node_i50yen(tria_i50yen(i)%v(1))%y) & + node_i50yen(tria_i50yen(i)%v(3))%x * & (node_i50yen(tria_i50yen(i)%v(1))%y-node_i50yen(tria_i50yen(i)%v(2))%y) ) if (area < 0.0_dp) negative_volumes = negative_volumes + 1 end do if (negative_volumes > 0) then if (n_attempts == 0) write(*,*) sf = sf*1.25_dp !Increased by 25% n_attempts = n_attempts+ 1 if (n_attempts > n_attempts_max+1) then write(*,*) " Tried to fix negative volume issue times, but failed..." write(*,*) " Try again with a larger stretching factor sf. Stop." stop endif write(*,'(a45,f10.3,a11,i5)') " --- Trying to fix negative volume with sf = ", & sf, " : Attempt ", n_attempts go to 200 else write(*,*) endif ! Construct boundary data. bnode_i50yen(1)%nnodes = 0 do i = bnode_50yen(1)%nnodes, 1, -1 bnode_i50yen(1)%nnodes = bnode_i50yen(1)%nnodes + 1 bnode_i50yen(1)%node(bnode_i50yen(1)%nnodes) = bnode_50yen(1)%node(i) end do bnode_i50yen(2)%nnodes = 0 bnode_i50yen(3)%nnodes = 1 do i = bnode_50yen(2)%nnodes-1, 1, -1 !Left outer boundary if ( node_50yen(bnode_50yen(2)%node(i))%x < zero) then bnode_i50yen(2)%nnodes = bnode_i50yen(2)%nnodes + 1 bnode_i50yen(2)%node(bnode_i50yen(2)%nnodes) = bnode_50yen(2)%node(i) !Right outer boundary else !First, collect the nodes below x-axis. if ( node_50yen(bnode_50yen(2)%node(i))%y < -1.0e-16_dp) then bnode_i50yen(3)%nnodes = bnode_i50yen(3)%nnodes + 1 bnode_i50yen(3)%node(bnode_i50yen(3)%nnodes) = bnode_50yen(2)%node(i) endif endif end do do i = bnode_50yen(2)%nnodes-1, 1, -1 if ( node_50yen(bnode_50yen(2)%node(i))%x < zero) then else !Right outer boundary !First, now collect the nodes above and on the x-axis. if ( node_50yen(bnode_50yen(2)%node(i))%y > -1.0e-16_dp) then bnode_i50yen(3)%nnodes = bnode_i50yen(3)%nnodes + 1 bnode_i50yen(3)%node(bnode_i50yen(3)%nnodes) = bnode_50yen(2)%node(i) endif endif end do !Add the first node to the right boundary bnode_i50yen(3)%node(1) = bnode_i50yen(2)%node(bnode_i50yen(2)%nnodes) !Add the last node to the left boundary bnode_i50yen(3)%nnodes = bnode_i50yen(3)%nnodes + 1 bnode_i50yen(3)%node(bnode_i50yen(3)%nnodes) = bnode_i50yen(2)%node(1) call write_k_file( filename_i50yen_k, k1_i50yen, k2_i50yen, k3_i50yen, nnodes_i50yen) call write_tecplot_file(filename_i50yen, node_i50yen, k1_i50yen, k2_i50yen, k3_i50yen, & kr_i50yen, nnodes_i50yen, tria_i50yen, ntrias_i50yen ) call write_bcmap_file(filename_i50yen_bcmap, 3, .false.) call write_grid_file(filename_i50yen_grid, node_i50yen, nnodes_i50yen, tria_i50yen, & ntrias_i50yen, 3,bnode_i50yen ) call write_vtk_file(filename_i50yen_vtk, nnodes_i50yen,node_i50yen,ntrias_i50yen,tria_i50yen) bnames(1) = "INNER" bnames(2) = "OUTER_LEFT" bnames(3) = "OUTER_RIGHT" call write_su2_file(filename_i50yen_su2, nnodes_i50yen,node_i50yen,ntrias_i50yen, & tria_i50yen, 3,bnode_i50yen,bnames ) deallocate(k1_50yen,k2_50yen,k3_50yen,kr_50yen,node_50yen,tria_50yen,bnode_50yen) !******************************************************************************* ! 5. Scale the inverted 50-yen grid. ! ! Here, we simply re-define the domain. ! Distance from the origin to the inner boundary = rin ! Distance from the origin to the outer boundary = rout ! ! Note: rin and rout are specified inside the progarm: ! ! rin = 0.5_dp ! rout = 50.0_dp ! ! Change these values to generate a 50-yen grid with desired distances. !******************************************************************************* write(*,*) "*****************************************************************" write(*,*) " 5. Scale the inverted 50-yen grid." write(*,*) "*****************************************************************" write(*,*) write(*,*) " Input - Minimum spacing desired, spacing_min (e.g., Try 0.001) :" write(*,*) " spacing_min = ?" read(*,*) spacing_min write(*,*) write(*,*) " spacing_min = ", spacing_min write(*,*) rin = 0.5_dp rout = 50.0_dp write(*,*) write(*,*) " Inner boundary is now located at r = rin = ", rin write(*,*) " Outer boundary is now located at r = rout = ", rout write(*,*) " [ You can change these values inside the code. ] " write(*,*) !------------------------------------------------------------------------------------ ! Adjust parameters (cp and possibly sf) to meet the desrired minimum spacing. r2 = 0.0_dp do i = 1, nnodes_i50yen if (kr_i50yen(i) == nr_gs-1) then r2 = sqrt( node_i50yen(i)%x**2 + node_i50yen(i)%y**2 ) endif end do cp = 0.05_dp !Just an initial guess. cp_min_reached = .false. k = 0 !------------------------------------------------------------------ adjust_cp : do !------------------------------------------------------------------ k = k + 1 if (cp < 1.0e-08_dp) then cp_min_reached = .true. sf = sf * 1.02_dp endif if (k > 1500) then write(*,*) " Sorry, I can't adjust the parameter cp... Stop" stop endif s = (r2-rmin)/(rmax-rmin) w = ( cp*s + (one-exp(sf*s))/(one-exp(sf)) )/(cp+one) r = ( w*rout + (one-w)*rin ) spacing = abs(r-rin) if( abs( spacing - spacing_min ) < 0.1_dp*spacing_min ) then !write(*,*) " Successfully adjusted. Exit." exit adjust_cp endif if (.not.cp_min_reached ) then if( spacing > spacing_min ) then cp = cp * 0.95_dp !write(*,*) " Decrease and try again: cp = ", cp, k else cp = cp * 1.05_dp !write(*,*) " Increase and try again: cp = ", cp, k endif endif !------------------------------------------------------------------ end do adjust_cp !------------------------------------------------------------------ write(*,*) " ---> Adjusted cp = ", cp if (cp_min_reached) write(*,*) " ---> Adjusted sf = ", sf write(*,*) !------------------------------------------------------------------ ! Now scale the grid with the determined parameters. do i = 1, nnodes_i50yen r = sqrt( node_i50yen(i)%x**2 + node_i50yen(i)%y**2 ) nx = node_i50yen(i)%x/r ny = node_i50yen(i)%y/r !Note: kr decreases from nr_gs to nr_gs_in towards the outer boundary. if (kr_i50yen(i) == nr_gs) then !Inner circle node_i50yen(i)%x = rin * nx node_i50yen(i)%y = rin * ny elseif (kr_i50yen(i) == nr_gs_in) then !Outer circle node_i50yen(i)%x = rout * nx node_i50yen(i)%y = rout * ny else s = (r-rmin)/(rmax-rmin) !Adjusted distance w = ( cp*s + (one-exp(sf*s))/(one-exp(sf)) )/(cp+one) r = ( w*rout + (one-w)*rin ) node_i50yen(i)%x = r * nx node_i50yen(i)%y = r * ny endif end do !If you wish, you can enable the following to perform grid smoothing to the resulting grid. ! [To enable this, simply rewrite '0 == 1' as '1 == 1'.] if (0 == 1) then call direct_smoothing(node_i50yen,nnodes_i50yen,tria_i50yen,ntrias_i50yen, 2,bnode_i50yen) endif call write_k_file( filename_si50yen_k, k1_i50yen, k2_i50yen, k3_i50yen, nnodes_i50yen) call write_tecplot_file(filename_si50yen , node_i50yen, k1_i50yen, k2_i50yen, k3_i50yen, & kr_i50yen, nnodes_i50yen, tria_i50yen, ntrias_i50yen ) call write_bcmap_file( filename_si50yen_bcmap, 3, .false.) call write_grid_file( filename_si50yen_grid, node_i50yen, nnodes_i50yen, tria_i50yen, & ntrias_i50yen, 3,bnode_i50yen ) call write_vtk_file(filename_si50yen_vtk, nnodes_i50yen,node_i50yen, & ntrias_i50yen,tria_i50yen ) bnames(1) = "INNER" bnames(2) = "OUTER_LEFT" bnames(3) = "OUTER_RIGHT" call write_su2_file(filename_si50yen_su2, nnodes_i50yen,node_i50yen,ntrias_i50yen, & tria_i50yen, 3,bnode_i50yen,bnames ) !******************************************************************************* ! 6. Perturb the scaled inverted 50-yen grid to generate an irregular grid. ! !******************************************************************************* write(*,*) "*****************************************************************" write(*,*) " 6. Perturb the scaled-inverted 50-yen grid." write(*,*) "*****************************************************************" write(*,*) ! Construct neighbor data allocate(nnghbrs(nnodes_i50yen)) !Note: We always have 6 neighbors in the interior node by construction ! for the grid considered here. But set the max to be 12 because ! we don't mind double-counting the neighbors here. allocate( nghbr(nnodes_i50yen,12)) nnghbrs = 0 nghbr = 0 ! Collect node-neighbor data. ! Here, we double count the interior neighbors, but we don't care. ! It is OK to have the same node twice for our purpose. do i = 1, ntrias_i50yen v = tria_i50yen(i)%v nnghbrs(v(1)) = nnghbrs(v(1)) + 1 nghbr(v(1),nnghbrs(v(1))) = v(2) nnghbrs(v(1)) = nnghbrs(v(1)) + 1 nghbr(v(1),nnghbrs(v(1))) = v(3) nnghbrs(v(2)) = nnghbrs(v(2)) + 1 nghbr(v(2),nnghbrs(v(2))) = v(1) nnghbrs(v(2)) = nnghbrs(v(2)) + 1 nghbr(v(2),nnghbrs(v(2))) = v(3) nnghbrs(v(3)) = nnghbrs(v(3)) + 1 nghbr(v(3),nnghbrs(v(3))) = v(2) nnghbrs(v(3)) = nnghbrs(v(3)) + 1 nghbr(v(3),nnghbrs(v(3))) = v(1) end do ! Now perturb nodes one by one by shrinking a randomly chosen edge. do i = 1, nnodes_i50yen !Skip boundary nodes (Do not perturb boundary nodes). if ( kr_i50yen(i) < nr_gs .and. kr_i50yen(i) > nr_gs_in ) then !Randomly choose a neighbor: k call random_number(rn) k = int( min( max(10.0_dp*rn,1.0_dp), 6.0_dp) ) ! so that k is in [1,6]. !Direction to the chosen neighbor. dx = node_i50yen(nghbr(i,k))%x - node_i50yen(i)%x dy = node_i50yen(nghbr(i,k))%y - node_i50yen(i)%y !Move the node i in that direction at 25% distance. node_i50yen(i)%x = node_i50yen(i)%x + 0.25_dp*dx node_i50yen(i)%y = node_i50yen(i)%y + 0.25_dp*dy endif end do call write_k_file( filename_psi50yen_k, k1_i50yen, k2_i50yen, k3_i50yen, nnodes_i50yen) call write_tecplot_file(filename_psi50yen , node_i50yen, k1_i50yen, k2_i50yen, k3_i50yen, & kr_i50yen, nnodes_i50yen, tria_i50yen, ntrias_i50yen ) call write_bcmap_file( filename_psi50yen_bcmap, 3, .false.) call write_grid_file( filename_psi50yen_grid, node_i50yen, nnodes_i50yen, tria_i50yen, & ntrias_i50yen, 3,bnode_i50yen ) call write_vtk_file(filename_psi50yen_vtk, nnodes_i50yen,node_i50yen, & ntrias_i50yen,tria_i50yen ) bnames(1) = "INNER" bnames(2) = "OUTER_LEFT" bnames(3) = "OUTER_RIGHT" call write_su2_file(filename_psi50yen_su2, nnodes_i50yen,node_i50yen,ntrias_i50yen, & tria_i50yen, 3,bnode_i50yen,bnames ) !******************************************************************************* ! 7. Swap diagonals in the perturbed-scaled-inverted 50-yen grid. ! Do it randomly! !******************************************************************************* write(*,*) "*****************************************************************" write(*,*) " 7. Swap diagonals in the perturbed-scaled-inverted 50-yen grid." write(*,*) "*****************************************************************" write(*,*) ! Construct node-to-element pointers: ! ! Example: A stencil centered at node i. ! Node i is shared by the elements [2,3,6,10,23,33,55,101]. ! ! o-------o-------------o ! / . 6 | 55 . | ! / 33 . | . 3 | ! o----------o-------------o ! \ 23 /i\ . 101 | ! \ / \ . | ! \ / 10 \ 2 . | ! o----------o---------o ! ! The information is stored as follows: ! ! node_i50yen(i)%nelms = 8 ! node_i50yen(i)%elm(1:8) = [2,3,6,10,23,33,55,101]. ! ! Initialization do i = 1, nnodes_i50yen node_i50yen(i)%nelms = 0 end do ! Add each element to the list at each node in that element. do i = 1, ntrias_i50yen v = tria_i50yen(i)%v do k = 1, 3 node_i50yen(v(k))%nelms = node_i50yen(v(k))%nelms + 1 node_i50yen(v(k))%elm(node_i50yen(v(k))%nelms) = i end do end do ! Construct element-neighbor data: do i = 1, ntrias_i50yen v = tria_i50yen(i)%v tria_i50yen(i)%e = 0 !=0 if there is no neighbor element. !Its element neighbors must be contained in the list of surrounding elements !of the nodes in that element. And so, we find neighbors by looping over the !surrounding elements of each node: v(k), k=1,2,3. !Loop over element-vertices: 1, 2, 3. do k = 1, 3 !Loop over the element-nghbors of the element vertex v(k): do j = 1, node_i50yen(v(k))%nelms ne = node_i50yen(v(k))%elm(j) v1 = tria_i50yen(ne)%v(1) v2 = tria_i50yen(ne)%v(2) v3 = tria_i50yen(ne)%v(3) !Let's find the element-neighbor: !Note: Neighbors are ordered as follows: ! ! /\ ! / \ ! / \ ! / \ ! / e1 \ ! v3/__________\v2 ! /\ /\ ! / \ i / \ ! / \ / \ ! / \ / \ ! / e2 \ / e3 \ ! /__________\/__________\ ! v1 ! !Neighbor e3: across v(1)-v(2) if ( v2 == v(1) .and. v1 == v(2) ) then tria_i50yen(i)%e(3) = ne elseif ( v3 == v(1) .and. v2 == v(2) ) then tria_i50yen(i)%e(3) = ne elseif ( v1 == v(1) .and. v3 == v(2) ) then tria_i50yen(i)%e(3) = ne else endif !Neighbor e1: across v(2)-v(3) if ( v2 == v(2) .and. v1 == v(3) ) then tria_i50yen(i)%e(1) = ne elseif ( v3 == v(2) .and. v2 == v(3) ) then tria_i50yen(i)%e(1) = ne elseif ( v1 == v(2) .and. v3 == v(3) ) then tria_i50yen(i)%e(1) = ne else endif !Neighbor e2: across v(3)-v(1) if ( v2 == v(3) .and. v1 == v(1) ) then tria_i50yen(i)%e(2) = ne elseif ( v3 == v(3) .and. v2 == v(1) ) then tria_i50yen(i)%e(2) = ne elseif ( v1 == v(3) .and. v3 == v(1) ) then tria_i50yen(i)%e(2) = ne else endif end do end do end do !Permutation matrix: ! k km(k,1) km(k,2) ! km(k,2) ------------------------- ! /\ 1 2 3 ! / \ 2 3 1 ! / \ 3 1 2 ! / \ ! / e1 \ ! k/__________\km(k,1) km(1,1) = 2 km(1,2) = 3 km(2,1) = 3 km(2,2) = 1 km(3,1) = 1 km(3,2) = 2 !We are now ready to swap the diagonals! n_swapped = 0 do i = 1, ntrias_i50yen v = tria_i50yen(i)%v e_nghbr : do k = 1, 3 ne = tria_i50yen(i)%e(k) if (ne == 0) cycle !Skip if no element neighbor. if (ne < i) cycle !Skip lower-numbered neighbors (just feel like; you don't have to). !Swap elements, i and ne, if (random_number-1/2) is positive!! call random_number(rn) if ( (rn-0.5_dp) > 0.0_dp ) then ! We try to swap elements i and ne. ! /\ ! / \ ! / \ ! / \ ! / e3 \ ! v4/__________\v3__________ ! /\ /\ / ! / \ ne / \ e2 / ! / \ / \ / ! / \ / \ / ! / e4 \ / i \ / ! /__________\/__________\/ ! v1 /V2 ! \ e1 / ! \ / ! \ / ! \ / ! \/ ! v3 = v( km(k,1) ) v2 = v( k ) v1 = v( km(k,2) ) e1 = tria_i50yen(i)%e( km(k,1) ) if (e1 == 0) cycle ! Let's skip if e1 doesn't exist. ke1 = 0 do j = 1, 3 if ( tria_i50yen(e1)%e(j) == i) ke1 = j end do e2 = tria_i50yen(i)%e( km(k,2) ) if (e2 == 0) cycle ! Let's skip if e2 doesn't exist. ke2 = 0 do j = 1, 3 if ( tria_i50yen(e2)%e(j) == i) ke2 = j end do kn = 0 do j = 1, 3 if ( tria_i50yen(ne)%e(j) == i) kn = j end do v4 = tria_i50yen(ne)%v( kn ) e3 = tria_i50yen(ne)%e( km(kn,1) ) if (e3 == 0) cycle ! Let's skip if e3 doesn't exist. ke3 = 0 do j = 1, 3 if ( tria_i50yen(e3)%e(j) == ne) ke3 = j end do e4 = tria_i50yen(ne)%e( km(kn,2) ) if (e4 == 0) cycle ! Let's skip if e4 doesn't exist. ke4 = 0 do j = 1, 3 if ( tria_i50yen(e4)%e(j) == ne) ke4 = j end do !But swap only if it leads to positive volumes. area124 = node_i50yen(v1)%x*( node_i50yen(v2)%y - node_i50yen(v4)%y ) & + node_i50yen(v2)%x*( node_i50yen(v4)%y - node_i50yen(v1)%y ) & + node_i50yen(v4)%x*( node_i50yen(v1)%y - node_i50yen(v2)%y ) area234 = node_i50yen(v2)%x*( node_i50yen(v3)%y - node_i50yen(v4)%y ) & + node_i50yen(v3)%x*( node_i50yen(v4)%y - node_i50yen(v2)%y ) & + node_i50yen(v4)%x*( node_i50yen(v2)%y - node_i50yen(v3)%y ) if ( area124 > 0.0_dp .and. area234 > 0.0_dp ) then !OK, now SWAP! ! We now swap and obtain the following grid: ! /\ ! / \ ! / \ ! / \ ! / e3 \ ! v4/__________\v3__________ ! /\ . \ / ! / \ . ne \ e2 / ! / \ . \ / ! / \ . \ / ! / e4 \ i . \ / ! /__________\___________\/ ! v1 /V2 ! \ e1 / ! \ / ! \ / ! \ / ! \/ ! !New element i tria_i50yen( i)%v(1) = v1 tria_i50yen( i)%v(2) = v2 tria_i50yen( i)%v(3) = v4 tria_i50yen( i)%e(1) = ne tria_i50yen( i)%e(2) = e4 tria_i50yen( i)%e(3) = e1 !New element ne tria_i50yen(ne)%v(1) = v2 tria_i50yen(ne)%v(2) = v3 tria_i50yen(ne)%v(3) = v4 tria_i50yen(ne)%e(1) = e3 tria_i50yen(ne)%e(2) = i tria_i50yen(ne)%e(3) = e2 !Update the element neighbor information in the neighbor elements: e1,e2,e3,e4, tria_i50yen(e1)%e(ke1) = i tria_i50yen(e2)%e(ke2) = ne tria_i50yen(e3)%e(ke3) = ne tria_i50yen(e4)%e(ke4) = i n_swapped = n_swapped + 1 !Let's skip other neighbors and move onto the next element. exit e_nghbr endif endif end do e_nghbr end do write(*,*) write(*,*) " The number of swapped edges = ", n_swapped write(*,*) call write_k_file( filename_spsi50yen_k, k1_i50yen, k2_i50yen, k3_i50yen, nnodes_i50yen) call write_tecplot_file(filename_spsi50yen , node_i50yen, k1_i50yen, k2_i50yen, k3_i50yen, & kr_i50yen, nnodes_i50yen, tria_i50yen, ntrias_i50yen ) call write_bcmap_file( filename_spsi50yen_bcmap, 3, .false.) call write_grid_file( filename_spsi50yen_grid, node_i50yen, nnodes_i50yen, tria_i50yen, & ntrias_i50yen, 3,bnode_i50yen ) call write_vtk_file(filename_spsi50yen_vtk, nnodes_i50yen,node_i50yen, & ntrias_i50yen,tria_i50yen ) bnames(1) = "INNER" bnames(2) = "OUTER_LEFT" bnames(3) = "OUTER_RIGHT" call write_su2_file(filename_spsi50yen_su2, nnodes_i50yen,node_i50yen,ntrias_i50yen, & tria_i50yen, 3,bnode_i50yen,bnames ) deallocate(k1_i50yen,k2_i50yen,k3_i50yen,kr_i50yen,node_i50yen,tria_i50yen,bnode_i50yen) !******************************************************************************* ! That's all for now... !******************************************************************************* write(*,*) write(*,*) "Grid generation successfuly completed. Bye." stop contains !******************************************************************************* ! This subroutine applies smoothing to the grid provided. !******************************************************************************* subroutine direct_smoothing(node,nnodes,tria,ntrias,nb,bnode) implicit none type(node_data_xy), dimension(:), intent(inout) :: node integer , intent(in) :: nnodes, ntrias type(tria_data) , dimension(:), intent(in) :: tria integer , intent(in) :: nb type(bnode_data) , dimension(:), intent(in) :: bnode real(dp), dimension(nnodes,2) :: dxy integer :: i_smoothing, max_smoothing, i, j real(dp) :: smoothing_factor, r integer , dimension(3) :: v integer , dimension(nnodes) :: bmark smoothing_factor = 0.05_dp max_smoothing = 1000 ! Construct boundary marker: bmark = 0 for the interior nodes, = 1 for boundary nodes. ! We'll use this to skip and thus avoid smoothing the boundary nodes. bmark = 0 do i = 1, nb do j = 1, bnode(i)%nnodes bmark(bnode(i)%node(j)) = 1 end do end do !---------------------------------------------------------------------------- ! Smoothing loop begins here. smoothing : do i_smoothing = 1, max_smoothing dxy = zero !Initialize the changes ! Accumulate the changes by looping over triangles do i = 1, ntrias v = tria(i)%v ! x-coordinate dxy(v(1),1) = dxy(v(1),1) + ( node(v(2))%x - node(v(1))%x ) dxy(v(1),1) = dxy(v(1),1) + ( node(v(3))%x - node(v(1))%x ) dxy(v(2),1) = dxy(v(2),1) + ( node(v(1))%x - node(v(2))%x ) dxy(v(2),1) = dxy(v(2),1) + ( node(v(3))%x - node(v(2))%x ) dxy(v(3),1) = dxy(v(3),1) + ( node(v(1))%x - node(v(3))%x ) dxy(v(3),1) = dxy(v(3),1) + ( node(v(2))%x - node(v(3))%x ) ! y-coordinate dxy(v(1),2) = dxy(v(1),2) + ( node(v(2))%y - node(v(1))%y ) dxy(v(1),2) = dxy(v(1),2) + ( node(v(3))%y - node(v(1))%y ) dxy(v(2),2) = dxy(v(2),2) + ( node(v(1))%y - node(v(2))%y ) dxy(v(2),2) = dxy(v(2),2) + ( node(v(3))%y - node(v(2))%y ) dxy(v(3),2) = dxy(v(3),2) + ( node(v(1))%y - node(v(3))%y ) dxy(v(3),2) = dxy(v(3),2) + ( node(v(2))%y - node(v(3))%y ) end do ! Make changes to each node except the boundary nodes. dxy_norm = -1000000.0_dp do i=1, nnodes if ( bmark(i) == 0 ) then node(i)%x = node(i)%x + smoothing_factor * dxy(i,1) node(i)%y = node(i)%y + smoothing_factor * dxy(i,2) ! L_inf norm of relative changes r = sqrt( node(i)%x**2 + node(i)%y**2 ) dxy_norm = max( dxy_norm, abs(dxy(i,1)**2 + dxy(i,2)**2)/r ) endif end do ! Exit if converged if ( dxy_norm < 1.0e-04) then write(*,*) " Smoothing converged at ", i_smoothing exit smoothing elseif (i_smoothing == max_smoothing) then write(*,*) " Smoothing didn't converge... ", " dxy_norm = ", dxy_norm endif end do smoothing ! Smoothing loop ends here. !---------------------------------------------------------------------------- end subroutine direct_smoothing !******************************************************************************* ! This subroutine writes a Tecplot file. !******************************************************************************* subroutine write_tecplot_file(filename,node,k1,k2,k3,kr,nnodes,tria,ntrias) implicit none character(80) , intent(in) :: filename type(node_data_xy), dimension(:), intent(in) :: node integer , dimension(:), intent(in) :: k1, k2, k3, kr integer , intent(in) :: nnodes, ntrias type(tria_data) , dimension(:), intent(in) :: tria open(unit=1, file=filename, status="unknown", iostat=os) write(1,*) 'TITLE = "grid"' write(1,*) 'VARIABLES = "x","y","k1","k2","k3","kr"' write(1,*) 'ZONE N=', nnodes,',E=', ntrias,' , ET=triangle, F=FEPOINT' ! Nodes do i = 1, nnodes write(1,'(2ES27.15,4i12)') node(i)%x, node(i)%y, k1(i),k2(i),k3(i),kr(i) end do ! Triangles do i = 1, ntrias write(1,'(3I10)') tria(i)%v(1), tria(i)%v(2), tria(i)%v(3) end do close(1) write(*,*) " --- Tecplot file written -> ", trim(filename) end subroutine write_tecplot_file !******************************************************************************* ! This subroutine writes a k file. !****************************************************************************** subroutine write_k_file(filename,k1,k2,k3,nnodes) implicit none integer , intent(in) :: nnodes character(80) , intent(in) :: filename integer, dimension(nnodes), intent(in) :: k1, k2, k3 open(unit=10, file=filename, status="unknown", iostat=os) write(10,*) nnodes do i = 1, nnodes write(10,'(4i15)') i, k1(i),k2(i),k3(i) end do close(10) write(*,*) " --- .k file written -> ", trim(filename) end subroutine write_k_file !******************************************************************************* ! This subroutine writes a bcmap file. ! Note: This is just an example. !****************************************************************************** subroutine write_bcmap_file(filename,nb, all_dirichlet) implicit none character(80), intent(in) :: filename integer , intent(in) :: nb logical , intent(in) :: all_dirichlet if (all_dirichlet) then open(unit=1, file=filename, status="unknown") write(1,*) ' Boundary Segment Boundary Condition' do i = 1, nb write(1,*) i, " dirichlet" end do close(1) return endif if (nb == 2) then open(unit=1, file=filename, status="unknown") write(1,*) ' Boundary Segment Boundary Condition' write(1,*) ' 1 "viscous_wall" ' write(1,*) ' 2 "freestream" ' close(1) elseif (nb == 3) then open(unit=1, file=filename, status="unknown") write(1,*) ' Boundary Segment Boundary Condition' write(1,*) ' 1 "viscous_wall" ' write(1,*) ' 2 "freestream" ' write(1,*) ' 3 "outflow_back_pressure"' close(1) else ! No such case is considered in this code. endif write(*,*) " --- .bcmap file written -> ", trim(filename) end subroutine write_bcmap_file !******************************************************************************* ! k-index in the radial direction. !****************************************************************************** subroutine k_radial(kr, k1,k2,k3,nnodes) implicit none integer , intent( in) :: nnodes integer, dimension(nnodes), intent( in) :: k1, k2, k3 integer, dimension(nnodes), intent(out) :: kr integer :: i integer :: sk12p, sk23p, sk31p integer :: sk12m, sk23m, sk31m do i = 1, nnodes sk12p = ( 1 + int_sign( k1(i)*k2(i) ) )/2 sk23p = ( 1 + int_sign( k2(i)*k3(i) ) )/2 sk31p = ( 1 + int_sign( k3(i)*k1(i) ) )/2 sk12m = ( 1 - int_sign( k1(i)*k2(i) ) )/2 sk23m = ( 1 - int_sign( k2(i)*k3(i) ) )/2 sk31m = ( 1 - int_sign( k3(i)*k1(i) ) )/2 if (sk12p > 0) then sk23p = 0 sk23m = 1 sk31p = 0 sk31m = 1 elseif (sk23p > 0) then sk12p = 0 sk12m = 1 sk31p = 0 sk31m = 1 elseif (sk31p > 0) then sk12p = 0 sk12m = 1 sk23p = 0 sk23m = 1 endif kr(i) = sk31m*sk23m*sk12p*abs(k3(i)) + & sk31m*sk12m*sk23p*abs(k1(i)) + & sk23m*sk12m*sk31p*abs(k2(i)) end do end subroutine k_radial !******************************************************************************** ! This subroutine writes a grid file to be read by a solver. ! NOTE: Unlike the tecplot file, this files contains boundary info. !******************************************************************************** subroutine write_grid_file(filename,node,nnodes,tria,ntrias, nb,bnode) implicit none character(80) , intent(in) :: filename type(node_data_xy), dimension(:), intent(in) :: node integer , intent(in) :: nnodes, ntrias type(tria_data) , dimension(:), intent(in) :: tria integer , intent(in) :: nb type(bnode_data) , dimension(:), intent(in) :: bnode integer :: os, negative_volumes real(dp) :: area negative_volumes = 0 !-------------------------------------------------------------------------------- open(unit=1, file=filename, status="unknown", iostat=os) !-------------------------------------------------------------------------------- ! Grid size: # of nodes, # of triangles, # of quadrilaterals write(1,*) nnodes, ntrias, 0 !-------------------------------------------------------------------------------- ! Node data do i = 1, nnodes write(1,*) node(i)%x, node(i)%y end do !-------------------------------------------------------------------------------- ! Triangle connectivity do i = 1, ntrias write(1,*) tria(i)%v(1), tria(i)%v(2), tria(i)%v(3) !Check the triangle volume: area = 0.5_dp*( node(tria(i)%v(1))%x*(node(tria(i)%v(2))%y-node(tria(i)%v(3))%y) & + node(tria(i)%v(2))%x*(node(tria(i)%v(3))%y-node(tria(i)%v(1))%y) & + node(tria(i)%v(3))%x*(node(tria(i)%v(1))%y-node(tria(i)%v(2))%y) ) if (area < 0.0_dp) negative_volumes = negative_volumes + 1 end do !-------------------------------------------------------------------------------- ! Boundary data: ! ! The number of boundary segments write(1,*) nb ! The number of nodes in each segment: do i = 1, nb write(1,*) bnode(i)%nnodes end do write(1,*) ! List of boundary nodes do i = 1, nb do j = 1, bnode(i)%nnodes write(1,*) bnode(i)%node(j) end do write(1,*) end do !-------------------------------------------------------------------------------- close(1) write(*,*) " --- .grid file written -> ", trim(filename) write(*,*) write(*,*) ">>>>> Grid generated ... " write(*,*) " nodes = ", nnodes write(*,*) " triangles = ", ntrias write(*,*) write(*,*) " Boundary segments = ", nb do i = 1, nb write(*,*) " Boundary # = ", i, ": nodes = ", bnode(i)%nnodes end do write(*,*) if (negative_volumes > 0) then write(*,*) " !!!!!!!!!!!!!!!!!!!!!!!!!!" write(*,*) " Negative volume found..... " write(*,*) " # of triangles with negative volume = ", negative_volumes write(*,*) " --> Try stronger stretching: Increase the factor sf and try again." write(*,*) " !!!!!!!!!!!!!!!!!!!!!!!!!!" else write(*,*) " All triangles have positive volume. Good." endif write(*,*) end subroutine write_grid_file !******************************************************************************** !******************************************************************************** !* This subroutine is useful to expand or shrink type(node_data_xy) arrays. !* !* Array, a, will be allocated if the requested dimension is 1 (i.e., n=1) !* expanded to the requested dimension, n, if n > dim(a). !* shrinked to the requested dimension, n, if n < dim(a). !* !******************************************************************************** subroutine my_alloc_ndxy_ptr(a,n) implicit none integer, intent(in ) :: n type(node_data_xy), dimension(:), pointer :: a integer :: i type(node_data_xy), dimension(:), pointer :: temp if (n <= 0) then write(*,*) "my_alloc_ndxy_ptr received non-positive dimension. Stop." stop endif ! If initial, allocate and return if (.not.(associated(a))) then allocate(a(n)) return endif ! If reallocation, create a pointer with a target of new dimension. allocate(temp(n)) temp(n)%gnode = 0 temp(n)%x = zero temp(n)%y = zero ! (1) Expand the array dimension if ( n > size(a) ) then do i = 1, size(a) temp(i)%gnode = a(i)%gnode temp(i)%x = a(i)%x temp(i)%y = a(i)%y end do ! (2) Shrink the array dimension: the extra data, a(n+1:size(x)), discarded. else do i = 1, n temp(i)%gnode = a(i)%gnode temp(i)%x = a(i)%x temp(i)%y = a(i)%y end do endif ! Destroy the target of a ! deallocate(a) ! Re-assign the pointer a => temp return end subroutine my_alloc_ndxy_ptr !******************************************************************************** !******************************************************************************** ! This is just to find the sign of an integer. function int_sign(n) implicit none integer, intent(in) :: n integer :: int_sign if (n >= 0) then int_sign = 1 else int_sign = -1 endif end function int_sign !******************************************************************************** !******************************************************************************* ! This subroutine writes a su2 grid file. ! ! Note: Nodes -> i = 0,1,2,...; Elements -> i = 0,1,2,... ! ! ! Identifier: ! Line 3 ! Triangle 5 ! Quadrilateral 9 ! Tetrahedral 10 ! Hexahedral 12 ! Prism 13 ! Pyramid 14 ! ! !******************************************************************************* subroutine write_su2_file(filename, nnodes,node, ntria,tria, nb,bnode,bnames ) character(80), intent(in) :: filename type(node_data_xy), dimension(:), intent(in) :: node integer , intent(in) :: nnodes, ntria type(tria_data) , dimension(:), intent(in) :: tria integer , intent(in) :: nb type(bnode_data) , dimension(:), intent(in) :: bnode character(80) , dimension(:), intent(in) :: bnames !Local variables integer :: i, ib, j, os write(*,*) write(*,*) " Writing .su2 file: ", trim(filename) open(unit=7, file=filename, status="unknown", iostat=os) write(7,*) "%" write(7,*) "% Problem dimension" write(7,*) "%" write(7,5) 2 5 format('NDIME= ',i12) write(7,*) "%" write(7,*) "% Inner element connectivity" write(7,10) ntria 10 format('NELEM= ',i12) !------------------------------------------------------------------------- ! Elements do i = 1, ntria write(7,'(4i20)') 5, tria(i)%v(1)-1, tria(i)%v(2)-1, tria(i)%v(3)-1 end do !-------------------------------------------------------------------------- ! Nodes write(7,*) "%" write(7,*) "% Node coordinates" write(7,*) "%" write(7,20) nnodes 20 format('NPOIN= ', i12) ! Nodes do i = 1, nnodes write(7,'(2es26.15)') node(i)%x, node(i)%y end do !-------------------------------------------------------------------------- ! Boundary 30 format('NMARK= ',i12) 40 format('MARKER_TAG= ',a) 50 format('MARKER_ELEMS= ', i12) write(7,*) "%" write(7,*) "% Boundary elements" write(7,*) "%" write(7,30) nb !# of boundary parts. do ib = 1, nb !------------------------- !Just to print on screen write(*,*) write(*,40) trim(bnames(ib)) !ib-th boundary-part name, e.g., "farfield". write(*,50) bnode(ib)%nnodes-1 !# of boundary elements (edges) !------------------------- write(7,40) trim(bnames(ib)) !ib-th boundary-part name, e.g., "farfield". write(7,50) bnode(ib)%nnodes-1 !# of boundary elements (edges) do j = 1, bnode(ib)%nnodes-1 write(7,'(3i20)') 3, bnode(ib)%node(j)-1, bnode(ib)%node(j+1)-1 end do end do close(7) end subroutine write_su2_file !******************************************************************************** !******************************************************************************* ! This subroutine writes a .vtk file for the grid whose name is defined by ! filename_vtk. ! ! 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,node, ntria,tria ) implicit none character(80), intent(in) :: filename type(node_data_xy), dimension(:), intent(in) :: node integer , intent(in) :: nnodes, ntria type(tria_data) , dimension(:), intent(in) :: tria !Local variables integer :: i, j, os !------------------------------------------------------------------------------ !------------------------------------------------------------------------------ !------------------------------------------------------------------------------ write(*,*) write(*,*) ' Writing .vtk file = ', trim(filename) write(*,*) !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)') node(j)%x, node(j)%y, zero end do !--------------------------------------------------------------------------- ! Cell information. !CELLS: # of total cells (tria+quad), total size of the cell list. write(8,'(a,i12,i12)') 'CELLS ', ntria, (3+1)*ntria ! Note: The latter is the number of integer values written below as data. ! 4 for triangles (# of vertices + 3 vertices), and ! 5 for quads (# of vertices + 4 vertices). !--------------------------------- ! 2.1 List of triangles (counterclockwise vertex ordering) if (ntria > 0) then ! (# of vertices = 3), 3 vertices in counterclockwise do i = 1, ntria write(8,'(a,4i12)') '3', tria(i)%v(1)-1, tria(i)%v(2)-1, tria(i)%v(3)-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 ', ntria !Triangle is classified as the cell type 5 in the .vtk format. if (ntria > 0) then do i = 1, ntria write(8,'(i3)') 5 end do endif !-------------------------------------------------------------------------------- ! NOTE: Commented out because there are no solution data. This part should be ! uncommented if this is used in a solver or if solution data are available ! and one would like to plot them. ! ! Field data (e.g., density, pressure, velocity)are added here for visualization. ! write(8,*) 'POINT_DATA ',nnodes ! ! FIELD dataName # of arrays (variables to plot) <--This is a commnet. ! write(8,*) 'FIELD FlowField ', 4 ! write(8,*) 'Density ', 1 , nnodes, ' double' ! do j = 1, nnodes ! write(8,'(es25.15)') w(j,1) ! end do ! write(8,*) 'X-velocity ', 1 , nnodes, ' double' ! do j = 1, nnodes ! write(8,'(es25.15)') w(j,2) ! end do ! write(8,*) 'Y-veloiity ', 1 , nnodes, ' double' ! do j = 1, nnodes ! write(8,'(es25.15)') w(j,3) ! end do ! write(8,*) 'Pressure ', 1 , nnodes, ' double' ! do j = 1, nnodes ! write(8,'(es25.15)') w(j,4) ! end do !--------------------------------------------------------------------------- !Close the output file. <--This is a commnet. close(8) end subroutine write_vtk_file !******************************************************************************** end program edu2d_50yen_tria_grids