#
# FixPoleDomain - a script to correct a bug in Gridgen V15.02 regarding
#   improper initialization of pole domains (SPR 3831).  This script
#   then attempts to force the proper initialization of a single pole
#   domain by creating a temporary database surface defining the pole
#   and projecting the domain to the database.  The database surface is
#   deleted when finished.  The script handles domains collapsed
#   to a point and domains collapsed to a line.  The latter case is
#   much more difficult, and complex line topologies may not be handled
#   correctly.  In such a case, simplify the non-pole edge of the domain
#   by joining connectors as much as possible.
#
#   By default, this script will prompt the user to pick a domain
#   interactively.  To run this in batch mode, define the variable
#   "FixPoleDomainBatch", source this file, and then call FixPoleDomain
#   with the pole domain.  For example:
#
#      set FixPoleDomainBatch 1
#      source [file join $thisScriptsDirectory FixPoleDomain.glf]
#      FixPoleDomain $poleDom
#
#   where "thisScriptsDirectory" is the directory containing this script,
#   and "poleDom" is the pole domain that did not get initialized
#   properly.
#
# Copyright 2004 (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.
#
package require PWI_Glyph 1.6.2

proc findNodes { con pt } {
  global FixPoleDomainChain
  set tol [gg::tolNode]
  for {set n 0} {$n < $FixPoleDomainChain(numNodes)} {incr n} {
    set v $FixPoleDomainChain(nodept-$n)
    if {$tol > [ggu::vec3Length [ggu::vec3Sub $v $pt]]} {
      lappend FixPoleDomainChain(nodecon-$n) $con
      set FixPoleDomainChain(nodecon-$n) \
        [lsort -unique $FixPoleDomainChain(nodecon-$n)]
      return $n
    }
  }
  # New node
  set n $FixPoleDomainChain(numNodes)
  set FixPoleDomainChain(nodept-$n) $pt
  set FixPoleDomainChain(nodecon-$n) $con
  incr FixPoleDomainChain(numNodes)
  return $n
}

proc FixPoleDomain { dom } {
  global FixPoleDomainChain
  # Get point from connector
  gg::domReport $dom data [list STRUCTURE REFERENCE]
  set tol [gg::tolGP]

  # Make sure this is a pole domain
  set nonPoleCons [list]
  foreach con $data(refCons) {
    set len [gg::conGetLength $con]
    if {$tol < $len} {
      lappend nonPoleCons $con
    }
  }

  # Check to see if there are any pole connectors
  if {[llength $nonPoleCons] == [llength $data(refCons)]} {
    puts "There does not appear to be any pole connectors in this domain."
    gg::abort
  # Handle case of point singularity domain
  } elseif {[llength $nonPoleCons] == 0} {
    # Get the singular point
    set con [lindex $data(refCons) 0]
    set pt [gg::conGetPt $con 1]
    # Create a temporary db surface to project to.  Create a square in the
    # XY plane through the point with a side length of 200 times the node
    # tolerance.
    set len [expr {100.0 * [gg::tolNode]} ]
    gg::dbCurveBegin -type 3D_Line
      set a [ggu::vec3Add $pt [ggu::vec3Scale [list -1.0 -1.0 0.0] $len]]
      gg::dbCurveAddPt $a
      set b [ggu::vec3Add $pt [ggu::vec3Scale [list -1.0 1.0 0.0] $len]]
      gg::dbCurveAddPt $b
    set db1 [gg::dbCurveEnd]
    gg::dbCurveBegin -type 3D_Line
      set a [ggu::vec3Add $pt [ggu::vec3Scale [list 1.0 -1.0 0.0] $len]]
      gg::dbCurveAddPt $a
      set b [ggu::vec3Add $pt [ggu::vec3Scale [list 1.0 1.0 0.0] $len]]
      gg::dbCurveAddPt $b
    set db2 [gg::dbCurveEnd]
    gg::dbSurfBegin -type RULED
      gg::dbSurfAddEnt $db1
      gg::dbSurfAddEnt $db2
    set db [gg::dbSurfEnd]
    # Delete the curves used to create the surface
    gg::dbDelete [list $db1 $db2]
    # Projecting the domain along with the connectors causes the domain to
    # be initialized with the correct values.
    gg::domProject $dom -type CLOSEST_PT -db $db
    # Delete the surface
    gg::dbDelete $db -force

    # Verify points were initialized correctly
    foreach {imax jmax} $data(dimensions) {break}
    for {set j 1} {$j <= $jmax} {incr j} {
      for {set i 1} {$i <= $imax} {incr i} {
	set v [gg::domGetPt $dom [list $i $j]]
	if {$tol < [ggu::vec3Length [ggu::vec3Sub $v $pt]]} {
	  error "The selected domain did not initialize properly"
	}
      }
    }
    set msg "The selected domain successfully collapsed to $pt."
  # Handle case of line singularity
  } else {
    # Map connectors to nodes
    set FixPoleDomainChain(numNodes) 0
    for {set c 0} {$c < [llength $nonPoleCons]} {incr c} {
      set con [lindex $nonPoleCons $c]
      set FixPoleDomainChain(con-$con-nodes) [list \
	[findNodes $con [gg::conGetPt $con -arc 0.0]] \
	[findNodes $con [gg::conGetPt $con -arc 1.0]] \
      ]
    }
    # Find the first node that is referenced only once.  If not such node
    # exists, either there's a tolerance problem or the non-pole edge loops
    # back on itself.
    foreach n [array names FixPoleDomainChain {nodecon-*}] {
      if {[llength $FixPoleDomainChain($n)] == 1} {
	set startNode [string range $n 8 end]
	break
      }
    }
    if {![info exists startNode]} {
      error "Could not determine edge linkage (1)"
    }
    # Starting with the first node, build up the list of points making up
    # the edge.
    set ptList [list $FixPoleDomainChain(nodept-$startNode)]
    set availableCons $nonPoleCons
    set lastNode $startNode
    set numAvail [llength $availableCons]
    while {$numAvail > 0} {
      # Loop through all of the available connectors looking for the connector
      # (should be only 1) that references the last node.
      for {set c 0} {$c < $numAvail} {incr c} {
	set con [lindex $availableCons $c]
	# Get this connector's nodes
	foreach {n0 n1} $FixPoleDomainChain(con-${con}-nodes) {break}
	# See if a node matches
	if {$n0 == $lastNode || $n1 == $lastNode} {
	  set dim [gg::conDim $con]
	  # If the first node matches, add the points in order (skipping the
	  # matching end point since it's already in the list).
	  if {$n0 == $lastNode} {
	    for {set i 2} {$i <= $dim} {incr i} {
	      lappend ptList [gg::conGetPt $con $i]
	    }
	    set lastNode $n1
	  # The last node matches, so add the points in reverse order (skipping
	  # the matching end point since it's already in the list).
	  } else {
	    for {set i [expr {$dim - 1}]} {$i > 0} {incr i -1} {
	      lappend ptList [gg::conGetPt $con $i]
	    }
	    set lastNode $n0
	  }
	  # Remove this connector from the list of available connectors.
	  set availableCons [lreplace $availableCons $c $c]
	  # Break out of the connector loop
	  break
	}
      }
      # If the list of available connectors did not shrink, that means we
      # did not find a node match.  This should not happen, but if it does,
      # it either means there is a tolerance problem or an odd topology
      # was encountered.
      if {$numAvail == [llength $availableCons]} {
	error "Could not determine edge linkage (2)"
      }
      # Update the number of available connectors
      set numAvail [llength $availableCons]
    }
    # Now find the limits of the points to determine the plane in which the
    # database will be created.
    set minPt [list 1000000.0 1000000.0 1000000.0]
    set maxPt [list -1000000.0 -1000000.0 -1000000.0]
    # Have to have at least 2 points.
    if {2 > [llength $ptList]} {
      error "Could not determine edge linkage (3)"
    }
    # Create the first curve.  This curve should follow the non-pole edge
    # of the domain.  Update the limits as we go.
    gg::dbCurveBegin -type 3D_LINE
      foreach pt $ptList {
	for {set i 0} {$i < 3} {incr i} {
	  set x [lindex $pt $i]
	  if {$x < [lindex $minPt $i]} {
	    set minPt [lreplace $minPt $i $i $x]
	  }
	  if {$x > [lindex $maxPt $i]} {
	    set maxPt [lreplace $maxPt $i $i $x]
	  }
	}
	gg::dbCurveAddPt $pt
      }
    set db1 [gg::dbCurveEnd]
    # We have the limits, so now determine the smallest extent along an
    # axial direction.  This will be the direction in which the edge curve
    # is swept.
    set delta [ggu::vec3Sub $maxPt $minPt]
    set dir [list 1.0 0.0 0.0]
    set i 0
    if {[lindex $delta 1] < [lindex $delta $i]} {
      set i 1
    }
    if {[lindex $delta 2] < [lindex $delta $i]} {
      set i 2
    }
    set dir [lreplace [list 0.0 0.0 0.0] $i $i [expr {100.0 * [gg::tolNode]}]]
    # Create the sweep curve
    gg::dbCurveBegin
      set pt [lindex $ptList 0]
      gg::dbCurveAddPt $pt
      gg::dbCurveAddPt [ggu::vec3Add $pt $dir]
    set db2 [gg::dbCurveEnd]
    # Create the sweep surface
    gg::dbSurfBegin -type SWEEP
      gg::dbSurfAddEnt $db1
      gg::dbSurfAddEnt $db2
    set db [gg::dbSurfEnd]
    # Delete the temporary curves
    gg::dbDelete [list $db1 $db2] -force
    # Project to the sweep surface
    gg::domProject $dom -type CLOSEST_PT -db $db
    # At this point, the edges of the domain should be right.  The interior
    # points are probably collapsed to a point since they were at the origin
    # to begin with.  By running the TFI solver, the interior points get updated
    # based on the now correct edge points.
    gg::domTFISolverRun $dom
    # Delete the sweep surface
    gg::dbDelete $db -force
    # Verify results by checking that the max and min area is zero
    set maxArea [expr {[gg::tolNode] * [gg::tolNode]}]
    gg::domExamine $dom data AREA -minimum 0.0 -maximum $maxArea
    if {$data(aboveRange) > 0 || $data(belowRange) > 0} {
      error "The selected domain did not initialize properly"
    }
    set msg "The selected domain successfully collapsed to a line."
  }
  return $msg
}

# To run this in batch mode, define the variable "FixPoleDomainBatch", 
# source this file, and then call FixPoleDomain with the pole domain.
# By default, this script will prompt the user to pick a domain
# interactively.
if {![info exists FixPoleDomainBatch]} {
  # Pick the pole domain
  set dom [gg::dispPick DOMAIN -multiple FALSE -type STRUCTURED \
    -title "PICK POLE DOMAIN" \
    -message "Pick the pole domain to be initialized"]
  puts "\n[FixPoleDomain $dom]"
  puts "\nIf the domain was previously constrained to a database, you will need"
  puts "to re-project the domain back to the desired database entities.\n"
}

#
# 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.
#