# Gridgen Glyph Script for generating grids around airfoils.
# Airfoils are generated in DreeseCODE's DesignFOIL.
#
# Written by: Pointwise (glyph@pointwise.com)
#
# Assumptions:
# If only one curve in file, airfoil has closed TE.
# If two curves in file, second curve closes TE.
# Grid is in XZ plane.
# O-grid topology
# Outer boundary shape is undefined (is defined by the grid method).
# Units are in inches.
#

# Initialization
gg::memClear
gg::aswSet GENERIC -dim 2
gg::defReset
gg::tolReset
gg::dispViewReset
gg::updatePolicy DISPLAY_AND_INPUT

# Directory from which script is run
set cwd [file dirname [info script]]

# USER PARAMETERS
# Set the file (from DesignFOIL) containing the airfoil shape.
set fname "NACA23012_gridgen.dat"
# Number of grid points around airfoil
set af_dim 101
# Minimum grid spacing as percentage of chord
set min_ds_pct 0.2
# First step size as percentage of minimum grid spacing
set ds1_pct 20.0
# Outer boundary distance as percentage of chord
set ff_pct 500.0
# Maximum number of extrusion steps
set num_steps 101
# Minimum number of grids points on a blunt TE
set min_te_pts 5
# CFD Solver
# PLOT3D, CFD++, CGNS-STRUCT, COBALT, FLUENT, STAR-CCM+, WIND
set cfd_solver "PLOT3D"
# CFD solver file name
# Note: Some solvers automatically add extensions; others don't.
if { $cfd_solver == "CFD++" } {
   set solver_file "CFD++_airfoil"
} elseif { $cfd_solver == "CGNS-STRUCT" } {
   set solver_file "CGNS_airfoil.cgns"
} elseif { $cfd_solver == "COBALT" } {
   set solver_file "COBALT_airfoil.inp"
} elseif { $cfd_solver == "FLUENT" } {
   set solver_file "FLUENT_airfoil.cas"
} elseif { $cfd_solver == "STAR-CCM+" } {
   set solver_file "STAR-CCM+_airfoil.inp"
} elseif { $cfd_solver == "WIND" } {
   set solver_file "WIND_airfoil.cgd"
} else {
   set cfd_solver "PLOT3D"
   set solver_file "PLOT3D_airfoil.x"
}

# Import the database.
set DB [gg::dbImport [file join $cwd $fname]]
set N_db [llength $DB]
if { $N_db == 1 } {
   # The airfoil consists of a single curve.
   set topo 1
} elseif { $N_db == 3 } {
   # There are three entities: two curves and a group.
   # The airfoil has two curves; second closes TE.
   set topo 2
} else {
   puts "BUG: Wrong number of entities in file."
}
set af_db [lindex $DB 0]
if { $topo == 2 } { 
   set te_db [lindex $DB 1] 
}

# Airfoil is in XZ plane; reorient the view accordingly.
gg::dispViewRot X -90
gg::dbDisp ALL -isolines {2 2}

# Find the LE (min X) and TE (max X) points on the airfoil.
set IJLimits [gg::dbGetIJLimits $af_db]
set I_min [expr int([lindex $IJLimits 0])]
set I_max [expr int([lindex $IJLimits 2])]
set LE_i [expr $I_min - 1]
set LE_pt [list 9999.0 0 0]
set TE_i [expr $I_max + 1]
set TE_pt [list -9999.0 0 0]
for {set i $I_min} {$i <= $I_max} {incr i} {
   set UV [gg::dbIJToUV  [list $i 1 $af_db]]
   set pt [gg::dbUVToXYZ $UV]
   if { [lindex $pt 0] < [lindex $LE_pt 0] } {
      set LE_i $i
      set LE_pt $pt
   } 
   if { [lindex $pt 0] > [lindex $TE_pt 0] } {
      set TE_i $i
      set TE_pt $pt
   } 
}

# Compute chord as distance between LE and TE points.
set dx [expr [lindex $TE_pt 0] - [lindex $LE_pt 0]]
set dy [expr [lindex $TE_pt 1] - [lindex $LE_pt 1]]
set dz [expr [lindex $TE_pt 2] - [lindex $LE_pt 2]]
set D2 [expr $dx * $dx + $dy * $dy + $dz * $dz]
set chord [expr sqrt($D2)]

# Create a connector on the airfoil.
# Place a CP at the LE for eventual use as a BP.
gg::conBegin
   gg::segBegin -type DB_LINE
      set UV [gg::dbIJToUV [list $I_min 1 $af_db]]
      gg::segAddControlPt -db $UV
      set UV [gg::dbIJToUV [list $LE_i  1 $af_db]]
      gg::segAddControlPt -db $UV
      set UV [gg::dbIJToUV [list $I_max 1 $af_db]]
      gg::segAddControlPt -db $UV
   gg::segEnd
set af_con [gg::conEnd]
gg::conDim $af_con $af_dim

# Compute min. grid spacing relative to chord.
set chord [expr [gg::conGetLength $af_con] / 2]
set min_ds [expr $chord * $min_ds_pct / 100 ]

if { $topo == 2 } {
   # Create and dimension a connector on the blunt TE.
   set te_con [gg::conOnDBEnt $te_db]
   set te_length [gg::conGetLength $te_con]
   set te_dim [expr int($te_length / $min_ds) + 1]
   gg::conDim $te_con $te_dim
   if { $te_dim < $min_te_pts } {
      # if the number of grid points on the TE is less than the minimum
      # redimension to the minimum and recompute the minimum spacing
      set te_dim $min_te_pts
      gg::conDim $te_con $te_dim
      set min_ds [expr $te_length / 2]
   }
}

# Distribute grid points on the airfoil.
# Set a breakpoint at the LE.
gg::conSetBreakPt $af_con $LE_pt
for {set i 1} {$i<=2} {incr i} {
   gg::conBeginSpacing $af_con -sub $i $min_ds
   gg::conEndSpacing   $af_con -sub $i $min_ds
}
gg::dispConGPS TRUE

# Extrude a domain from the connectors.
# Set a limit on the farfield height.
# First step size is a fraction of TE spacing.
set ff_height [expr $chord * $ff_pct / 100]
set ds1 [expr $min_ds * $ds1_pct / 100]
if { $topo == 1 } {
   set af_edge [list $af_con]
} else {
   set af_edge [list $af_con $te_con]
}
gg::domExtrusionBegin $af_edge -edge
   gg::domExtrusionMode HYPERBOLIC
   gg::domExtrusionAtt -s_init $ds1 \
                       -march_plane {0 1 0} \
                       -kb_smoothing NOMINAL \
                       -stop_height $ff_height
   gg::domExtrusionStep -result ExtResult $num_steps
set af_dom [gg::domExtrusionEnd]
gg::dispConGPS FALSE

# Create a 2D block for application of BCs.
gg::blkBegin -type STRUCTURED
   gg::faceBegin
      gg::faceAddDom $af_dom
   gg::faceEnd
set af_blk [gg::blkEnd]

# Set CFD solver and boundary conditions.
if { $cfd_solver != "PLOT3D" } {
   gg::aswSet $cfd_solver -dim 2
   set options ""
   # The farfield connector should be the last in the list.
   gg::domReport $af_dom Report REFERENCE
   set last [expr [llength $Report(refCons)] -1]
   set ff_con [lindex $Report(refCons) $last]
   # Set the BC using the types specific to the chosen solver.
   if { $cfd_solver == "CFD++" } {
      gg::aswSetBC $af_edge "Wall"
      gg::aswSetBC $ff_con "Unspecified"
   } elseif { $cfd_solver == "CGNS-STRUCT" } {
      gg::aswSetBC $af_edge "Wall"
      gg::aswSetBC $ff_con "Unspecified"
   } elseif { $cfd_solver == "COBALT" } {
      gg::aswSetBC $af_edge "Solid Wall"
      gg::aswSetBC $ff_con "Farfield"
   } elseif { $cfd_solver == "FLUENT" } {
      gg::aswSetBC $af_edge "Wall"
      gg::aswSetBC $ff_con "Pressure Far Field"
   } elseif { $cfd_solver == "STAR-CCM+" } {
      gg::aswSetBC $af_edge "Wall"
      gg::aswSetBC $ff_con "Unspecified"
   } elseif { $cfd_solver == "WIND" } {
      gg::aswSetBC $af_edge "Viscous wall"
      gg::aswSetBC $ff_con "Freestream"
      set options "-units inches"
   } else {
      puts "Can't determine CFD solver."
   }
}

if { $cfd_solver == "PLOT3D" } {
   # Export a PLOT3D file (grid only, no BCs).
   gg::domExport $af_dom [file join $cwd $solver_file] -style PLOT3D
} else {
   # Export CFD solver file(s).  (Some solvers write more than 1 file.)
   eval "gg::aswExport [file join $cwd $solver_file] $options"
}
puts "Wrote file $solver_file"