tcltk vtk optimize

# Routine that extracts a VOI based on the sections and integrates a volume based on
# a the extracted volume, a cutoff and the elevations planes
#
set showvol 1
set savestrgrd 0
toplevel .contour
frame .contour.buttons -borderwidth 3 -relief ridge
pack .contour.buttons -side top -padx 8 -pady 5 -ipadx 3 -ipady 1

frame .contour.buttons.f1
label .contour.buttons.l1 -text "Thresh Upper:" -width $lwidth -anchor w
entry .contour.buttons.r1 -textvariable volcontvalU -width $rwidth
pack .contour.buttons.f1
pack .contour.buttons.l1 .contour.buttons.r1 -in .contour.buttons.f1 -side left -padx $xpad -pady $ypad

frame .contour.buttons.f1a
label .contour.buttons.l1a -text "Thresh Lower:" -width $lwidth -anchor w
entry .contour.buttons.r1a -textvariable volcontvalL -width $rwidth
pack .contour.buttons.f1a
pack .contour.buttons.l1a .contour.buttons.r1a -in .contour.buttons.f1a -side left -padx $xpad -pady $ypad

frame .contour.buttons.f1b
label .contour.buttons.l1b -text "BVWirr:" -width $lwidth -anchor w
entry .contour.buttons.r1b -textvariable bvwirr -width $rwidth
pack .contour.buttons.f1b
pack .contour.buttons.l1b .contour.buttons.r1b -in .contour.buttons.f1b -side left -padx $xpad -pady $ypad

frame .contour.buttons.f1c
label .contour.buttons.l1c -text "Bo:" -width $lwidth -anchor w
entry .contour.buttons.r1c -textvariable shrinkage -width $rwidth
pack .contour.buttons.f1c
pack .contour.buttons.l1c .contour.buttons.r1c -in .contour.buttons.f1c -side left -padx $xpad -pady $ypad

frame .contour.buttons.f2
button .contour.buttons.f2.b1 -text "Contour" -width 10 -command contour
pack .contour.buttons.f2
pack .contour.buttons.f2.b1 -side left -padx 4 -pady 4

frame .contour.save -borderwidth 3 -relief ridge
pack .contour.save -side top -padx 8 -pady 5 -ipadx 3 -ipady 1

frame .contour.save.f1
label .contour.save.l1 -text "Save Poly:" -width $lwidth -anchor w
entry .contour.save.r1 -textvariable volfile -width $rwidth
pack .contour.save.f1
pack .contour.save.l1 .contour.save.r1 -in .contour.save.f1 -side left -padx $xpad -pady $ypad

frame .contour.save.f2
button .contour.save.f2.b1 -text "Save" -width 10 -command volsave
pack .contour.save.f2
pack .contour.save.f2.b1 -side left -padx 4 -pady 4

frame .contour.save.f1b
label .contour.save.l1b -text "Save Map:" -width $lwidth -anchor w
entry .contour.save.r1b -textvariable mapfile -width $rwidth
pack .contour.save.f1b
pack .contour.save.l1b .contour.save.r1b -in .contour.save.f1b -side left -padx $xpad -pady $ypad

frame .contour.save.f2b
button .contour.save.f2b.b1 -text "Save" -width 10 -command mapsave
pack .contour.save.f2b
pack .contour.save.f2b.b1 -side left -padx 4 -pady 4

frame .contour.save.f1c
label .contour.save.l1c -text "Save StrGrd:" -width $lwidth -anchor w
entry .contour.save.r1c -textvariable strgrdfile -width $rwidth
pack .contour.save.f1c
pack .contour.save.l1c .contour.save.r1c -in .contour.save.f1c -side left -padx $xpad -pady $ypad

frame .contour.save.f2c
button .contour.save.f2c.b1 -text "Save" -width 10 -command strgrdsave
pack .contour.save.f2c
pack .contour.save.f2c.b1 -side left -padx 4 -pady 4

radiobutton .contour.buttons.rb1 -variable threshtype -text "Higher" -value "1"
pack .contour.buttons.rb1 -side left
radiobutton .contour.buttons.rb2 -variable threshtype -text "Between" -value "0"
pack .contour.buttons.rb2 -side left
radiobutton .contour.buttons.rb3 -variable threshtype -text "Lower" -value "-1"
pack .contour.buttons.rb3 -side left

checkbutton .contour.buttons.cb1 -variable showvol -text "Show Vol" -command showvol
pack .contour.buttons.cb1 -side top -anchor w

text .contour.t -height 20 -width 60 -background white
pack .contour.t

set volcontvalU 0.0
set volcontvalL 0.0
set bvwirr 0.06
set shrinkage 1.1
set volrunflg 0
set volfile "VolPoly.vtk"
set mapfile "Map.dat"
set strgrdfile "StrGrd-Extract.vtk"

proc contour {} {
  global volcontvalU volcontvalL status dx dy gx1 gx2 gy1 gy2 z1 z2 elev1 elev2 volrunflg bvwirr shrinkage sumvol cellvol
  global grdindex sumnetthick grdsize volxdim volxmin volymin threshtype passcut sumpor ptscal sw sumhc zgriddim sumthick
  set status "Contouring" ; update
  .contour.t delete 1.0 end
  if {$volrunflg>0} {
    contGeom Delete
    p2c Delete
    thresh Delete
    connect Delete
    transConn Delete
    lut2 Delete
    geomFilt Delete
    connTriang Delete
    connSmooth Delete
    connNormals Delete
    moldMapper Delete
    ren1 RemoveActor moldActor
    moldActor Delete
  }  
  vtkExtractGrid contGeom
    contGeom SetInput [reader GetOutput]
    contGeom SetVOI $gx1 $gx2 $gy1 $gy2 $z1 $z2
    contGeom Update
  set voldim [[contGeom GetOutput] GetDimensions]
  set volxdim [lindex $voldim 0]
  set volydim [lindex $voldim 1]
  set volzdim [lindex $voldim 2]
  .contour.t insert end "GX1= $gx1\n"
  .contour.t insert end "GX2= $gx2\n"
  .contour.t insert end "GY1= $gy1\n"
  .contour.t insert end "GY2= $gy2\n"
  .contour.t insert end "Zmin= $z1\n"
  .contour.t insert end "Zmax= $z2\n"
  .contour.t insert end "Elev1= $elev1\n"
  .contour.t insert end "Elev2= $elev2\n"
  .contour.t insert end "Xdim= $volxdim\n"
  .contour.t insert end "Ydim= $volydim\n"
  .contour.t insert end "Zdim= $volzdim\n"
  set ncells [expr $volxdim*$volydim*$volzdim]
  .contour.t insert end "nCells= $ncells\n"
  set passcut 0
  set sumpor 0.0
  set sumhc 0.0
  set avgpor 0.0
  set sumthick 0.0
  set sumvol 0.0
  set grdsize [expr $volxdim*$volydim]
  set repor 0.0
  for {set i4 0} {$i4 < $grdsize} {incr i4} {
    set sumnetthick($i4) 0.0
  }
  set volpt0 [[contGeom GetOutput] GetPoint 0]
  set volxmin [lindex $volpt0 0]
  set volymin [lindex $volpt0 1]
  for {set i5 0} {$i5 < $ncells} {incr i5} {
    set calck [expr $i5/$grdsize]
    set grdindex [expr $i5-$calck*$grdsize]
    set volpt [[contGeom GetOutput] GetPoint $i5]
    set volz [lindex $volpt 2]
    if {$calck < [expr $volzdim-1]} {
      set i5b [expr $i5 + $grdsize]
    } else {
      set i5b [expr $i5 - $grdsize]
    }
    set volptb [[contGeom GetOutput] GetPoint $i5b]
    set volzb [lindex $volptb 2]
    set zgriddim [expr abs($volz - $volzb)]
    set cellvol [expr $dx*$dy*$zgriddim]      
    set ptscal [[[[contGeom GetOutput] GetPointData] GetScalars] GetValue $i5]
    
if {$volz>=$elev1 && $volz<=$elev2} {
  if {$threshtype==1 && $ptscal>= $volcontvalU} {
    sumcells
  } elseif {$threshtype==0 && $ptscal>=$volcontvalL && $ptscal<=$volcontvalU} {
    sumcells
  } elseif {$threshtype==-1 && $ptscal<=$volcontvalL} {
    sumcells
  }
} else {

  eval [[[[contGeom GetOutput] GetPointData] GetScalars] SetValue $i5 $repor]
}
  }
  
  #set avgpor [expr $sumpor/$sumvol]
  #set porevolbbls [expr $sumpor*0.178108]
  #set hcvolbbls [expr $sumhc*0.178108]
  #.contour.t insert end "PassCont= $passcut\n"
  #.contour.t insert end "AvgPor= $avgpor\n"
  #.contour.t insert end "PoreVol= $sumpor\n"
  #.contour.t insert end "PoreVol(bbl)= $porevolbbls\n"
  #.contour.t insert end "HCVol= $sumhc\n"
  #.contour.t insert end "HCVol(bbl)= $hcvolbbls\n"
  update    
  
  vtkPointDataToCellData p2c
    p2c SetInput [contGeom GetOutput]
  vtkThreshold thresh
    thresh SetInput [contGeom GetOutput]
    if {$threshtype==1} {thresh ThresholdByUpper $volcontvalU}
    if {$threshtype==0} {thresh ThresholdBetween $volcontvalL $volcontvalU}
    if {$threshtype==-1} {thresh ThresholdByLower $volcontvalL}
  vtkConnectivityFilter connect
    connect SetInput [thresh GetOutput]
    connect SetExtractionModeToAllRegions
    connect ColorRegionsOn
    connect Update
  set scalrange2 [[connect GetOutput] GetScalarRange]
  set scalmin2 [lindex $scalrange2 0]
  set scalmax2 [lindex $scalrange2 1]
  .debug.t insert end "Scalmin2= $scalmin2\n"
  .debug.t insert end "Scalmax2= $scalmax2\n" ; update    
  vtkLookupTable lut2
    lut2 SetHueRange 0.667 0.0
    lut2 SetTableRange $scalmin2 $scalmax2
    lut2 Build

  vtkGeometryFilter geomFilt
    geomFilt SetInput [connect GetOutput]
  vtkTriangleFilter connTriang
    connTriang SetInput [geomFilt GetOutput]
  vtkSmoothPolyDataFilter connSmooth
    connSmooth SetInput [connTriang GetOutput]
    connSmooth BoundarySmoothingOn

  vtkTransformFilter transConn
    transConn SetInput [connSmooth GetOutput]
    transConn SetTransform trans1
    transConn Update  
  vtkPolyDataNormals connNormals
    connNormals SetInput [transConn GetOutput]  
  vtkPolyDataMapper moldMapper
    moldMapper SetInput [connNormals GetOutput]
    moldMapper ScalarVisibilityOn
    moldMapper SetLookupTable lut2
    moldMapper SetScalarRange $scalmin2 $scalmax2
  vtkActor moldActor
    moldActor SetMapper moldMapper
  set volrunflg 1
  showvol 
  set status "Ready" ; update
}

proc sumcells {} {
  global passcut sumpor ptscal sw bvwirr sumhc shrinkage sumnetthick grdindex zgriddim sumthick sumvol cellvol
  incr passcut
  set sumthick [expr $sumthick+$zgriddim]
  set sumvol [expr $sumvol+$cellvol]
  set sumpor [expr $sumpor+$ptscal*$cellvol]
  if {$ptscal <= 0.0} {
    set sw 1.0
  } else {
    set sw [expr $bvwirr/$ptscal]
  }
  set sumhc [expr $sumhc+$cellvol*$ptscal*(1-$sw)/$shrinkage]
  set sumnetthick($grdindex) [expr $sumnetthick($grdindex)+$zgriddim] ; # feet meeting cutoff
#  set sumnetthick($grdindex) [expr $sumnetthick($grdindex)+$zgriddim*(1.0-$ptscal)] ; # net sand thIckness
#  set sumnetthick($grdindex) [expr $sumnetthick($grdindex)+$zgriddim*$ptscal] ; # porosity thickness
  
}

proc showvol {} {
  global showvol
  set status "Displaying Volume" ; update
    set existsVol [[ren1 GetActors] IsItemPresent moldActor]
    if {$showvol==0 && $existsVol>0} {ren1 RemoveActor moldActor}
    if {$showvol==1 && $existsVol<=0} {ren1 AddActor moldActor}
    renWin Render
    
  if {$showvol==0 && $existsVol>0} {
    ren1 RemoveActor moldActor
  } elseif {$showvol==1 && $existsVol<=0} {
    ren1 AddActor moldActor
  }
}

set volsaveflg 0
proc volsave {} {
  global status volsaveflg volfile
  set status "Saving Vol Poly" ; update
  if {$volsaveflg>0} {writeVolPoly Delete}
  vtkPolyDataWriter writeVolPoly
    writeVolPoly SetInput [connSmooth GetOutput]
    writeVolPoly SetFileName data/$volfile
    writeVolPoly SetFileTypeToBinary
    writeVolPoly Update
  set volsaveflg 1    
  set status "Ready" ; update
}

proc mapsave {} {
  global status sumnetthick mapfile grdsize volxdim dx dy volxmin volymin
  set status "Saving Map Grid" ; update
  set mapFileId [open $mapfile w]
  for {set i 0} {$i < $grdsize} {incr i} {
#    set calcgrdj [expr double($i/$volxdim)]
#    set grdj [expr double(int($calcgrdj)*$dy)+$volymin]
#    set calcgrdi [expr $calcgrdj*double($volxdim)]
#    set grdi [expr double(($i-$calcgrdi)*$dx)+$volxmin]    
    set grdpt [[contGeom GetOutput] GetPoint $i]
    set grdx [lindex $grdpt 0]    
    set grdy [lindex $grdpt 1]
    puts $mapFileId "$grdx  $grdy  $sumnetthick($i)"
  }
  close $mapFileId    
  set status "Ready" ; update
}

set savestrgrdflg 0
proc strgrdsave {} {
  global status strgrdfile savestrgrdflg gx1 gx2 gy1 gy2 z1 z2 
  set status "Saving Structured Grid" ; update
  if {$savestrgrdflg==1} {
    writeStrgrd Delete
  }
  vtkExtractGrid contGeom
    contGeom SetInput [reader GetOutput]
    contGeom SetVOI $gx1 $gx2 $gy1 $gy2 $z1 $z2
    contGeom Update  
  vtkStructuredGridWriter writeStrgrd
   writeStrgrd SetInput [contGeom GetOutput]
   writeStrgrd SetFileName $strgrdfile
#   writeStrgrd SetFileTypeToBinary
   writeStrgrd Update
  set savestrgrdflg 1
  contGeom Delete
  set status "Ready" ; update
}