#
# Copyright 2005 (c) Pointwise, Inc.
# All rights reserved.
# 
# This sample Gridgen script is not supported by Pointwise, Inc.
# It is provided freely for demonstration purposes only.  
# SEE THE WARRANTY DISCLAIMER AT THE BOTTOM OF THIS FILE.
#
# set GLF [file rootname [info script]]
set GLF "   quadWarp"

# convert radians to degrees
proc rad2deg { r } {
   set d [expr $r * 180 / 3.1415927]
   return $d
}

# get domain's dimensions
proc getDomainDimension { d } {
   gg::domReport $d i STRUCTURE
   set dim $i(dimensions)
   return $dim
}

# Compute the warp of the given quad cell.
# Warp is the angle between normal vectors of the two triangles that
# are obtained by diagonalizing the quad.
# First, compute warp relative to one diagonal.
# Second, compute warp relative to other diagonal.
# Return greater of the two.
# (i,j+1)   (i+1,j+1)
# C---------D
# |         |
# |         |
# |         |
# A---------B
# (i,j)  (i+1,j)
proc warp { d i j } {
   set T [gg::tolGP]
   set ij [list $i $j]
   set A [gg::domGetPt $d $ij]
   set ij [list [expr $i + 1] $j]
   set B [gg::domGetPt $d $ij]
   set ij [list $i [expr $j + 1] ]
   set C [gg::domGetPt $d $ij]
   set ij [list [expr $i + 1] [expr $j + 1] ]
   set D [gg::domGetPt $d $ij]
  
   # diagonal CB
   # triangle ABC
   set E1 [ggu::vec3Sub $B $A]
   set E2 [ggu::vec3Sub $C $A]
   set X1 [ggu::vec3Cross $E1 $E2]
   set L1 [ggu::vec3Length $X1]
   # triangle DCB
   set E1 [ggu::vec3Sub $D $C]
   set E2 [ggu::vec3Sub $D $B]
   set X2 [ggu::vec3Cross $E1 $E2]
   set L2 [ggu::vec3Length $X2]
   if { $L1 > $T && $L2 > $T } {
      # compute angle
      set DP [ggu::vec3Dot $X1 $X2]
      set cosTheta [expr $DP / $L1 / $L2]
      if { $cosTheta > 1 } {
         set cosTheta 1.0 
      } elseif { $cosTheta < -1 } {
         set cosTheta -1.0
      }
      set w_rad [expr acos($cosTheta) ]
      set w1_deg [rad2deg $w_rad]
   } else {
      set w1_deg 0
   }

   # diagonal AD
   # triangle BDA
   set E1 [ggu::vec3Sub $B $A]
   set E2 [ggu::vec3Sub $D $B]
   set X1 [ggu::vec3Cross $E1 $E2]
   set L1 [ggu::vec3Length $X1]
   # triangle CAD
   set E1 [ggu::vec3Sub $D $C]
   set E2 [ggu::vec3Sub $C $A]
   set X2 [ggu::vec3Cross $E1 $E2]
   set L2 [ggu::vec3Length $X2]
   if { $L1 > $T && $L2 > $T } {
      # compute angle
      set DP [ggu::vec3Dot $X1 $X2]
      set cosTheta [expr $DP / $L1 / $L2]
      if { $cosTheta > 1 } {
         set cosTheta 1.0 
      } elseif { $cosTheta < -1 } {
         set cosTheta -1.0
      }
      set w_rad [expr acos($cosTheta) ]
      set w2_deg [rad2deg $w_rad]
   } else {
      set w2_deg 0
   }

   # return greater of the two angles
   if { $w1_deg > $w2_deg } {
      set w_deg $w1_deg 
   } else {
      set w_deg $w2_deg
   }
   return $w_deg
}

# user picks domains
if { [catch { set D [gg::dispPick DOMAIN \
                   -title "quadWarp" \
                   -message "Pick domains for quad warp." \
                   -description DIMENSIONS \
                   -type STRUCTURED] } ] } {
   puts "$GLF: Aborting because there aren't any structured domains."
   exit
}
set n_D [llength $D]
if { $n_D == 0 } {
   puts "$GLF: Aborting because no domains were picked."
   exit
}

# loop through each domain, compute warp of each cell,
# track min/max warp per domain
# track min/max warp for all picked domains
puts " "
puts "script: quadWarp"
puts "------------------------------------------------"
puts "         Q U A D   C E L L S      W A R P (deg) "
puts "domain      #i      #j   Total      min      max"
puts "------  ------  ------  ------  -------  -------"
set total_quads 0
set g_max 0.0
set g_min 360.0
foreach dom $D {
   set dnum [lindex [gg::domIdToNum $dom] 0]
   set w_max 0.0
   set w_min 360.0
   set dim [getDomainDimension $dom]
   set Imax [lindex $dim 0]
   set Jmax [lindex $dim 1]
   for {set j 1} {$j < $Jmax} {incr j} {
      for {set i 1} {$i < $Imax} {incr i} {
         set w [warp $dom $i $j]
	 if { $w < $w_min } {
            set w_min $w
         }
	 if { $w > $w_max } {
            set w_max $w
         }
      }
   }
   # write data for current domain
   set n_quads [expr [expr $Imax - 1] * [expr $Jmax - 1]]
   puts [format "%6d  %6d  %6d  %6d  %7.4f  %7.4f" \
         $dnum [expr $Imax - 1] [expr $Jmax - 1] $n_quads $w_min $w_max]
   set total_quads [expr $total_quads + $n_quads]
   if { $w_min < $g_min } {
      set g_min $w_min
   }
   if { $w_max > $g_max } {
      set g_max $w_max
   }
}
puts "------------------------------------------------"
puts [format "%6d  %14s  %6d  %7.4f  %7.4f" \
      $n_D " " $total_quads $g_min $g_max]
puts "------------------------------------------------"

#
# DISCLAIMER:
# TO THE MAXIMUM EXTENT PERMITTED BY APPLICABLE LAW, POINTWISE DISCLAIMS
# ALL WARRANTIES, EITHER EXPRESS OR IMPLIED, INCLUDING, BUT NOT LIMITED
# TO, IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
# PURPOSE, WITH REGARD TO THIS SCRIPT.  TO THE MAXIMUM EXTENT PERMITTED 
# BY APPLICABLE LAW, IN NO EVENT SHALL POINTWISE BE LIABLE TO ANY PARTY 
# FOR ANY SPECIAL, INCIDENTAL, INDIRECT, OR CONSEQUENTIAL DAMAGES 
# WHATSOEVER (INCLUDING, WITHOUT LIMITATION, DAMAGES FOR LOSS OF 
# BUSINESS INFORMATION, OR ANY OTHER PECUNIARY LOSS) ARISING OUT OF THE 
# USE OF OR INABILITY TO USE THIS SCRIPT EVEN IF POINTWISE HAS BEEN 
# ADVISED OF THE POSSIBILITY OF SUCH DAMAGES AND REGARDLESS OF THE 
# FAULT OR NEGLIGENCE OF POINTWISE.
#