Posted to tcl by ahype18 at Thu Jan 08 22:08:32 GMT 2015view raw

  1. # Routine that extracts a VOI based on the sections and integrates a volume based on
  2. # a the extracted volume, a cutoff and the elevations planes
  3. #
  4. set showvol 1
  5. set savestrgrd 0
  6. toplevel .contour
  7. frame .contour.buttons -borderwidth 3 -relief ridge
  8. pack .contour.buttons -side top -padx 8 -pady 5 -ipadx 3 -ipady 1
  9.  
  10. frame .contour.buttons.f1
  11. label .contour.buttons.l1 -text "Thresh Upper:" -width $lwidth -anchor w
  12. entry .contour.buttons.r1 -textvariable volcontvalU -width $rwidth
  13. pack .contour.buttons.f1
  14. pack .contour.buttons.l1 .contour.buttons.r1 -in .contour.buttons.f1 -side left -padx $xpad -pady $ypad
  15.  
  16. frame .contour.buttons.f1a
  17. label .contour.buttons.l1a -text "Thresh Lower:" -width $lwidth -anchor w
  18. entry .contour.buttons.r1a -textvariable volcontvalL -width $rwidth
  19. pack .contour.buttons.f1a
  20. pack .contour.buttons.l1a .contour.buttons.r1a -in .contour.buttons.f1a -side left -padx $xpad -pady $ypad
  21.  
  22. frame .contour.buttons.f1b
  23. label .contour.buttons.l1b -text "BVWirr:" -width $lwidth -anchor w
  24. entry .contour.buttons.r1b -textvariable bvwirr -width $rwidth
  25. pack .contour.buttons.f1b
  26. pack .contour.buttons.l1b .contour.buttons.r1b -in .contour.buttons.f1b -side left -padx $xpad -pady $ypad
  27.  
  28. frame .contour.buttons.f1c
  29. label .contour.buttons.l1c -text "Bo:" -width $lwidth -anchor w
  30. entry .contour.buttons.r1c -textvariable shrinkage -width $rwidth
  31. pack .contour.buttons.f1c
  32. pack .contour.buttons.l1c .contour.buttons.r1c -in .contour.buttons.f1c -side left -padx $xpad -pady $ypad
  33.  
  34. frame .contour.buttons.f2
  35. button .contour.buttons.f2.b1 -text "Contour" -width 10 -command contour
  36. pack .contour.buttons.f2
  37. pack .contour.buttons.f2.b1 -side left -padx 4 -pady 4
  38.  
  39. frame .contour.save -borderwidth 3 -relief ridge
  40. pack .contour.save -side top -padx 8 -pady 5 -ipadx 3 -ipady 1
  41.  
  42. frame .contour.save.f1
  43. label .contour.save.l1 -text "Save Poly:" -width $lwidth -anchor w
  44. entry .contour.save.r1 -textvariable volfile -width $rwidth
  45. pack .contour.save.f1
  46. pack .contour.save.l1 .contour.save.r1 -in .contour.save.f1 -side left -padx $xpad -pady $ypad
  47.  
  48. frame .contour.save.f2
  49. button .contour.save.f2.b1 -text "Save" -width 10 -command volsave
  50. pack .contour.save.f2
  51. pack .contour.save.f2.b1 -side left -padx 4 -pady 4
  52.  
  53. frame .contour.save.f1b
  54. label .contour.save.l1b -text "Save Map:" -width $lwidth -anchor w
  55. entry .contour.save.r1b -textvariable mapfile -width $rwidth
  56. pack .contour.save.f1b
  57. pack .contour.save.l1b .contour.save.r1b -in .contour.save.f1b -side left -padx $xpad -pady $ypad
  58.  
  59. frame .contour.save.f2b
  60. button .contour.save.f2b.b1 -text "Save" -width 10 -command mapsave
  61. pack .contour.save.f2b
  62. pack .contour.save.f2b.b1 -side left -padx 4 -pady 4
  63.  
  64. frame .contour.save.f1c
  65. label .contour.save.l1c -text "Save StrGrd:" -width $lwidth -anchor w
  66. entry .contour.save.r1c -textvariable strgrdfile -width $rwidth
  67. pack .contour.save.f1c
  68. pack .contour.save.l1c .contour.save.r1c -in .contour.save.f1c -side left -padx $xpad -pady $ypad
  69.  
  70. frame .contour.save.f2c
  71. button .contour.save.f2c.b1 -text "Save" -width 10 -command strgrdsave
  72. pack .contour.save.f2c
  73. pack .contour.save.f2c.b1 -side left -padx 4 -pady 4
  74.  
  75. radiobutton .contour.buttons.rb1 -variable threshtype -text "Higher" -value "1"
  76. pack .contour.buttons.rb1 -side left
  77. radiobutton .contour.buttons.rb2 -variable threshtype -text "Between" -value "0"
  78. pack .contour.buttons.rb2 -side left
  79. radiobutton .contour.buttons.rb3 -variable threshtype -text "Lower" -value "-1"
  80. pack .contour.buttons.rb3 -side left
  81.  
  82. checkbutton .contour.buttons.cb1 -variable showvol -text "Show Vol" -command showvol
  83. pack .contour.buttons.cb1 -side top -anchor w
  84.  
  85. text .contour.t -height 20 -width 60 -background white
  86. pack .contour.t
  87.  
  88. set volcontvalU 0.0
  89. set volcontvalL 0.0
  90. set bvwirr 0.06
  91. set shrinkage 1.1
  92. set volrunflg 0
  93. set volfile "VolPoly.vtk"
  94. set mapfile "Map.dat"
  95. set strgrdfile "StrGrd-Extract.vtk"
  96.  
  97. proc contour {} {
  98. global volcontvalU volcontvalL status dx dy gx1 gx2 gy1 gy2 z1 z2 elev1 elev2 volrunflg bvwirr shrinkage sumvol cellvol
  99. global grdindex sumnetthick grdsize volxdim volxmin volymin threshtype passcut sumpor ptscal sw sumhc zgriddim sumthick
  100. set status "Contouring" ; update
  101. .contour.t delete 1.0 end
  102. if {$volrunflg>0} {
  103. contGeom Delete
  104. p2c Delete
  105. thresh Delete
  106. connect Delete
  107. transConn Delete
  108. lut2 Delete
  109. geomFilt Delete
  110. connTriang Delete
  111. connSmooth Delete
  112. connNormals Delete
  113. moldMapper Delete
  114. ren1 RemoveActor moldActor
  115. moldActor Delete
  116. }
  117. vtkExtractGrid contGeom
  118. contGeom SetInput [reader GetOutput]
  119. contGeom SetVOI $gx1 $gx2 $gy1 $gy2 $z1 $z2
  120. contGeom Update
  121. set voldim [[contGeom GetOutput] GetDimensions]
  122. set volxdim [lindex $voldim 0]
  123. set volydim [lindex $voldim 1]
  124. set volzdim [lindex $voldim 2]
  125. .contour.t insert end "GX1= $gx1\n"
  126. .contour.t insert end "GX2= $gx2\n"
  127. .contour.t insert end "GY1= $gy1\n"
  128. .contour.t insert end "GY2= $gy2\n"
  129. .contour.t insert end "Zmin= $z1\n"
  130. .contour.t insert end "Zmax= $z2\n"
  131. .contour.t insert end "Elev1= $elev1\n"
  132. .contour.t insert end "Elev2= $elev2\n"
  133. .contour.t insert end "Xdim= $volxdim\n"
  134. .contour.t insert end "Ydim= $volydim\n"
  135. .contour.t insert end "Zdim= $volzdim\n"
  136. set ncells [expr $volxdim*$volydim*$volzdim]
  137. .contour.t insert end "nCells= $ncells\n"
  138. set passcut 0
  139. set sumpor 0.0
  140. set sumhc 0.0
  141. set avgpor 0.0
  142. set sumthick 0.0
  143. set sumvol 0.0
  144. set grdsize [expr $volxdim*$volydim]
  145. set repor 0.0
  146. for {set i4 0} {$i4 < $grdsize} {incr i4} {
  147. set sumnetthick($i4) 0.0
  148. }
  149. set volpt0 [[contGeom GetOutput] GetPoint 0]
  150. set volxmin [lindex $volpt0 0]
  151. set volymin [lindex $volpt0 1]
  152. for {set i5 0} {$i5 < $ncells} {incr i5} {
  153. set calck [expr $i5/$grdsize]
  154. set grdindex [expr $i5-$calck*$grdsize]
  155. set volpt [[contGeom GetOutput] GetPoint $i5]
  156. set volz [lindex $volpt 2]
  157. if {$calck < [expr $volzdim-1]} {
  158. set i5b [expr $i5 + $grdsize]
  159. } else {
  160. set i5b [expr $i5 - $grdsize]
  161. }
  162. set volptb [[contGeom GetOutput] GetPoint $i5b]
  163. set volzb [lindex $volptb 2]
  164. set zgriddim [expr abs($volz - $volzb)]
  165. set cellvol [expr $dx*$dy*$zgriddim]
  166. set ptscal [[[[contGeom GetOutput] GetPointData] GetScalars] GetValue $i5]
  167.  
  168. if {$volz>=$elev1 && $volz<=$elev2} {
  169. if {$threshtype==1 && $ptscal>= $volcontvalU} {
  170. sumcells
  171. } elseif {$threshtype==0 && $ptscal>=$volcontvalL && $ptscal<=$volcontvalU} {
  172. sumcells
  173. } elseif {$threshtype==-1 && $ptscal<=$volcontvalL} {
  174. sumcells
  175. }
  176. } else {
  177.  
  178. eval [[[[contGeom GetOutput] GetPointData] GetScalars] SetValue $i5 $repor]
  179. }
  180. }
  181.  
  182. #set avgpor [expr $sumpor/$sumvol]
  183. #set porevolbbls [expr $sumpor*0.178108]
  184. #set hcvolbbls [expr $sumhc*0.178108]
  185. #.contour.t insert end "PassCont= $passcut\n"
  186. #.contour.t insert end "AvgPor= $avgpor\n"
  187. #.contour.t insert end "PoreVol= $sumpor\n"
  188. #.contour.t insert end "PoreVol(bbl)= $porevolbbls\n"
  189. #.contour.t insert end "HCVol= $sumhc\n"
  190. #.contour.t insert end "HCVol(bbl)= $hcvolbbls\n"
  191. update
  192.  
  193. vtkPointDataToCellData p2c
  194. p2c SetInput [contGeom GetOutput]
  195. vtkThreshold thresh
  196. thresh SetInput [contGeom GetOutput]
  197. if {$threshtype==1} {thresh ThresholdByUpper $volcontvalU}
  198. if {$threshtype==0} {thresh ThresholdBetween $volcontvalL $volcontvalU}
  199. if {$threshtype==-1} {thresh ThresholdByLower $volcontvalL}
  200. vtkConnectivityFilter connect
  201. connect SetInput [thresh GetOutput]
  202. connect SetExtractionModeToAllRegions
  203. connect ColorRegionsOn
  204. connect Update
  205. set scalrange2 [[connect GetOutput] GetScalarRange]
  206. set scalmin2 [lindex $scalrange2 0]
  207. set scalmax2 [lindex $scalrange2 1]
  208. .debug.t insert end "Scalmin2= $scalmin2\n"
  209. .debug.t insert end "Scalmax2= $scalmax2\n" ; update
  210. vtkLookupTable lut2
  211. lut2 SetHueRange 0.667 0.0
  212. lut2 SetTableRange $scalmin2 $scalmax2
  213. lut2 Build
  214.  
  215. vtkGeometryFilter geomFilt
  216. geomFilt SetInput [connect GetOutput]
  217. vtkTriangleFilter connTriang
  218. connTriang SetInput [geomFilt GetOutput]
  219. vtkSmoothPolyDataFilter connSmooth
  220. connSmooth SetInput [connTriang GetOutput]
  221. connSmooth BoundarySmoothingOn
  222.  
  223. vtkTransformFilter transConn
  224. transConn SetInput [connSmooth GetOutput]
  225. transConn SetTransform trans1
  226. transConn Update
  227. vtkPolyDataNormals connNormals
  228. connNormals SetInput [transConn GetOutput]
  229. vtkPolyDataMapper moldMapper
  230. moldMapper SetInput [connNormals GetOutput]
  231. moldMapper ScalarVisibilityOn
  232. moldMapper SetLookupTable lut2
  233. moldMapper SetScalarRange $scalmin2 $scalmax2
  234. vtkActor moldActor
  235. moldActor SetMapper moldMapper
  236. set volrunflg 1
  237. showvol
  238. set status "Ready" ; update
  239. }
  240.  
  241. proc sumcells {} {
  242. global passcut sumpor ptscal sw bvwirr sumhc shrinkage sumnetthick grdindex zgriddim sumthick sumvol cellvol
  243. incr passcut
  244. set sumthick [expr $sumthick+$zgriddim]
  245. set sumvol [expr $sumvol+$cellvol]
  246. set sumpor [expr $sumpor+$ptscal*$cellvol]
  247. if {$ptscal <= 0.0} {
  248. set sw 1.0
  249. } else {
  250. set sw [expr $bvwirr/$ptscal]
  251. }
  252. set sumhc [expr $sumhc+$cellvol*$ptscal*(1-$sw)/$shrinkage]
  253. set sumnetthick($grdindex) [expr $sumnetthick($grdindex)+$zgriddim] ; # feet meeting cutoff
  254. # set sumnetthick($grdindex) [expr $sumnetthick($grdindex)+$zgriddim*(1.0-$ptscal)] ; # net sand thIckness
  255. # set sumnetthick($grdindex) [expr $sumnetthick($grdindex)+$zgriddim*$ptscal] ; # porosity thickness
  256.  
  257. }
  258.  
  259. proc showvol {} {
  260. global showvol
  261. set status "Displaying Volume" ; update
  262. set existsVol [[ren1 GetActors] IsItemPresent moldActor]
  263. if {$showvol==0 && $existsVol>0} {ren1 RemoveActor moldActor}
  264. if {$showvol==1 && $existsVol<=0} {ren1 AddActor moldActor}
  265. renWin Render
  266.  
  267. if {$showvol==0 && $existsVol>0} {
  268. ren1 RemoveActor moldActor
  269. } elseif {$showvol==1 && $existsVol<=0} {
  270. ren1 AddActor moldActor
  271. }
  272. }
  273.  
  274. set volsaveflg 0
  275. proc volsave {} {
  276. global status volsaveflg volfile
  277. set status "Saving Vol Poly" ; update
  278. if {$volsaveflg>0} {writeVolPoly Delete}
  279. vtkPolyDataWriter writeVolPoly
  280. writeVolPoly SetInput [connSmooth GetOutput]
  281. writeVolPoly SetFileName data/$volfile
  282. writeVolPoly SetFileTypeToBinary
  283. writeVolPoly Update
  284. set volsaveflg 1
  285. set status "Ready" ; update
  286. }
  287.  
  288. proc mapsave {} {
  289. global status sumnetthick mapfile grdsize volxdim dx dy volxmin volymin
  290. set status "Saving Map Grid" ; update
  291. set mapFileId [open $mapfile w]
  292. for {set i 0} {$i < $grdsize} {incr i} {
  293. # set calcgrdj [expr double($i/$volxdim)]
  294. # set grdj [expr double(int($calcgrdj)*$dy)+$volymin]
  295. # set calcgrdi [expr $calcgrdj*double($volxdim)]
  296. # set grdi [expr double(($i-$calcgrdi)*$dx)+$volxmin]
  297. set grdpt [[contGeom GetOutput] GetPoint $i]
  298. set grdx [lindex $grdpt 0]
  299. set grdy [lindex $grdpt 1]
  300. puts $mapFileId "$grdx $grdy $sumnetthick($i)"
  301. }
  302. close $mapFileId
  303. set status "Ready" ; update
  304. }
  305.  
  306. set savestrgrdflg 0
  307. proc strgrdsave {} {
  308. global status strgrdfile savestrgrdflg gx1 gx2 gy1 gy2 z1 z2
  309. set status "Saving Structured Grid" ; update
  310. if {$savestrgrdflg==1} {
  311. writeStrgrd Delete
  312. }
  313. vtkExtractGrid contGeom
  314. contGeom SetInput [reader GetOutput]
  315. contGeom SetVOI $gx1 $gx2 $gy1 $gy2 $z1 $z2
  316. contGeom Update
  317. vtkStructuredGridWriter writeStrgrd
  318. writeStrgrd SetInput [contGeom GetOutput]
  319. writeStrgrd SetFileName $strgrdfile
  320. # writeStrgrd SetFileTypeToBinary
  321. writeStrgrd Update
  322. set savestrgrdflg 1
  323. contGeom Delete
  324. set status "Ready" ; update
  325. }
  326.