Posted to tcl by ahype18 at Fri Jan 09 02:32:39 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.  
  155. set timesteptimer [time {expr {$i5/$grdsize}}]
  156. .contour.t insert end "time= $timesteptimer\n" ; update
  157.  
  158. set grdindex [expr {$i5-$calck*$grdsize}]
  159.  
  160. set timesteptimer [time {expr {$i5-$calck*$grdsize}}]
  161. .contour.t insert end "time= $timesteptimer\n" ; update
  162.  
  163. set volpt [[contGeom GetOutput] GetPoint $i5]
  164.  
  165. set volz [lindex $volpt 2]
  166. if {$calck < [expr {$volzdim-1}]} {
  167.  
  168. set timesteptimer [time {expr {$i5 + $grdsize}}]
  169. .contour.t insert end "time= $timesteptimer\n" ; update
  170.  
  171. set i5b [expr {$i5 + $grdsize}]
  172. } else {
  173. set i5b [expr {$i5 - $grdsize}]
  174. }
  175. set volptb [[contGeom GetOutput] GetPoint $i5b]
  176. set volzb [lindex $volptb 2]
  177. set zgriddim [expr {abs($volz - $volzb)}]
  178.  
  179. set timesteptimer [time {expr {abs($volz - $volzb)}}]
  180. .contour.t insert end "time= $timesteptimer\n" ; update
  181.  
  182. set cellvol [expr {$dx*$dy*$zgriddim}]
  183.  
  184. set timesteptimer [time {expr {$dx*$dy*$zgriddim}}]
  185. .contour.t insert end "time= $timesteptimer\n" ; update
  186.  
  187. set ptscal [[[[contGeom GetOutput] GetPointData] GetScalars] GetValue $i5]
  188.  
  189. if {$volz>=$elev1 && $volz<=$elev2} {
  190. if {$threshtype==1 && $ptscal>= $volcontvalU} {
  191. sumcells
  192. } elseif {$threshtype==0 && $ptscal>=$volcontvalL && $ptscal<=$volcontvalU} {
  193. sumcells
  194. } elseif {$threshtype==-1 && $ptscal<=$volcontvalL} {
  195. sumcells
  196. }
  197. } else {
  198.  
  199. eval [[[[contGeom GetOutput] GetPointData] GetScalars] SetValue $i5 $repor]
  200. }
  201. }
  202.  
  203. #set avgpor [expr $sumpor/$sumvol]
  204. #set porevolbbls [expr $sumpor*0.178108]
  205. #set hcvolbbls [expr $sumhc*0.178108]
  206. #.contour.t insert end "PassCont= $passcut\n"
  207. #.contour.t insert end "AvgPor= $avgpor\n"
  208. #.contour.t insert end "PoreVol= $sumpor\n"
  209. #.contour.t insert end "PoreVol(bbl)= $porevolbbls\n"
  210. #.contour.t insert end "HCVol= $sumhc\n"
  211. #.contour.t insert end "HCVol(bbl)= $hcvolbbls\n"
  212. update
  213.  
  214. vtkPointDataToCellData p2c
  215. p2c SetInput [contGeom GetOutput]
  216. vtkThreshold thresh
  217. thresh SetInput [contGeom GetOutput]
  218. if {$threshtype==1} {thresh ThresholdByUpper $volcontvalU}
  219. if {$threshtype==0} {thresh ThresholdBetween $volcontvalL $volcontvalU}
  220. if {$threshtype==-1} {thresh ThresholdByLower $volcontvalL}
  221. vtkConnectivityFilter connect
  222. connect SetInput [thresh GetOutput]
  223. connect SetExtractionModeToAllRegions
  224. connect ColorRegionsOn
  225. connect Update
  226. set scalrange2 [[connect GetOutput] GetScalarRange]
  227. set scalmin2 [lindex $scalrange2 0]
  228. set scalmax2 [lindex $scalrange2 1]
  229. .debug.t insert end "Scalmin2= $scalmin2\n"
  230. .debug.t insert end "Scalmax2= $scalmax2\n" ; update
  231. vtkLookupTable lut2
  232. lut2 SetHueRange 0.667 0.0
  233. lut2 SetTableRange $scalmin2 $scalmax2
  234. lut2 Build
  235.  
  236. vtkGeometryFilter geomFilt
  237. geomFilt SetInput [connect GetOutput]
  238. vtkTriangleFilter connTriang
  239. connTriang SetInput [geomFilt GetOutput]
  240. vtkSmoothPolyDataFilter connSmooth
  241. connSmooth SetInput [connTriang GetOutput]
  242. connSmooth BoundarySmoothingOn
  243.  
  244. vtkTransformFilter transConn
  245. transConn SetInput [connSmooth GetOutput]
  246. transConn SetTransform trans1
  247. transConn Update
  248. vtkPolyDataNormals connNormals
  249. connNormals SetInput [transConn GetOutput]
  250. vtkPolyDataMapper moldMapper
  251. moldMapper SetInput [connNormals GetOutput]
  252. moldMapper ScalarVisibilityOn
  253. moldMapper SetLookupTable lut2
  254. moldMapper SetScalarRange $scalmin2 $scalmax2
  255. vtkActor moldActor
  256. moldActor SetMapper moldMapper
  257. set volrunflg 1
  258. showvol
  259. set status "Ready" ; update
  260. }
  261.  
  262. proc sumcells {} {
  263. global passcut sumpor ptscal sw bvwirr sumhc shrinkage sumnetthick grdindex zgriddim sumthick sumvol cellvol
  264. incr passcut
  265. set sumthick [expr {$sumthick+$zgriddim}]
  266. set sumvol [expr {$sumvol+$cellvol}]
  267. set sumpor [expr {$sumpor+$ptscal*$cellvol}]
  268. if {$ptscal <= 0.0} {
  269. set sw 1.0
  270. } else {
  271. set sw [expr {$bvwirr/$ptscal}]
  272. }
  273. set sumhc [expr {$sumhc+$cellvol*$ptscal*(1-$sw)/$shrinkage}]
  274. set sumnetthick($grdindex) [expr {$sumnetthick($grdindex)+$zgriddim}] ; # feet meeting cutoff
  275. # set sumnetthick($grdindex) [expr $sumnetthick($grdindex)+$zgriddim*(1.0-$ptscal)] ; # net sand thIckness
  276. # set sumnetthick($grdindex) [expr $sumnetthick($grdindex)+$zgriddim*$ptscal] ; # porosity thickness
  277.  
  278. }
  279.  
  280. proc showvol {} {
  281. global showvol
  282. set status "Displaying Volume" ; update
  283. set existsVol [[ren1 GetActors] IsItemPresent moldActor]
  284. if {$showvol==0 && $existsVol>0} {ren1 RemoveActor moldActor}
  285. if {$showvol==1 && $existsVol<=0} {ren1 AddActor moldActor}
  286. renWin Render
  287.  
  288. if {$showvol==0 && $existsVol>0} {
  289. ren1 RemoveActor moldActor
  290. } elseif {$showvol==1 && $existsVol<=0} {
  291. ren1 AddActor moldActor
  292. }
  293. }
  294.  
  295. set volsaveflg 0
  296. proc volsave {} {
  297. global status volsaveflg volfile
  298. set status "Saving Vol Poly" ; update
  299. if {$volsaveflg>0} {writeVolPoly Delete}
  300. vtkPolyDataWriter writeVolPoly
  301. writeVolPoly SetInput [connSmooth GetOutput]
  302. writeVolPoly SetFileName data/$volfile
  303. writeVolPoly SetFileTypeToBinary
  304. writeVolPoly Update
  305. set volsaveflg 1
  306. set status "Ready" ; update
  307. }
  308.  
  309. proc mapsave {} {
  310. global status sumnetthick mapfile grdsize volxdim dx dy volxmin volymin
  311. set status "Saving Map Grid" ; update
  312. set mapFileId [open $mapfile w]
  313. for {set i 0} {$i < $grdsize} {incr i} {
  314. # set calcgrdj [expr double($i/$volxdim)]
  315. # set grdj [expr double(int($calcgrdj)*$dy)+$volymin]
  316. # set calcgrdi [expr $calcgrdj*double($volxdim)]
  317. # set grdi [expr double(($i-$calcgrdi)*$dx)+$volxmin]
  318. set grdpt [[contGeom GetOutput] GetPoint $i]
  319. set grdx [lindex $grdpt 0]
  320. set grdy [lindex $grdpt 1]
  321. puts $mapFileId "$grdx $grdy $sumnetthick($i)"
  322. }
  323. close $mapFileId
  324. set status "Ready" ; update
  325. }
  326.  
  327. set savestrgrdflg 0
  328. proc strgrdsave {} {
  329. global status strgrdfile savestrgrdflg gx1 gx2 gy1 gy2 z1 z2
  330. set status "Saving Structured Grid" ; update
  331. if {$savestrgrdflg==1} {
  332. writeStrgrd Delete
  333. }
  334. vtkExtractGrid contGeom
  335. contGeom SetInput [reader GetOutput]
  336. contGeom SetVOI $gx1 $gx2 $gy1 $gy2 $z1 $z2
  337. contGeom Update
  338. vtkStructuredGridWriter writeStrgrd
  339. writeStrgrd SetInput [contGeom GetOutput]
  340. writeStrgrd SetFileName $strgrdfile
  341. # writeStrgrd SetFileTypeToBinary
  342. writeStrgrd Update
  343. set savestrgrdflg 1
  344. contGeom Delete
  345. set status "Ready" ; update
  346. }
  347.