API tools faq
paste
Advertisement
Guest User

Path Lines

a guest
Apr 28th, 2025
21
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 32.11 KB | None | 0 0
raw download clone embed print report
  1. from collections import defaultdict
  2. from paraview.util.vtkAlgorithm import *
  3. import numpy as np
  4. import vtk
  5.  
  6. from vtkmodules.numpy_interface import algorithms as algs
  7. from vtkmodules.vtkCommonDataModel import vtkPolyData
  8. from vtkmodules.util.vtkAlgorithm import VTKPythonAlgorithmBase
  9. from vtkmodules.numpy_interface import dataset_adapter as dsa
  10. from vtkmodules.vtkCommonTransforms import vtkTransform
  11. from vtkmodules.vtkFiltersGeneral import vtkTransformFilter
  12. import random
  13. from vtkmodules.vtkFiltersCore import vtkAppendPolyData
  14.  
  15. from vtkmodules import vtkCommonCore
  16. import sys
  17.  
  18.  
  19. @smproxy.filter(label="New Pathlines")
  20. @smproperty.input(name="Input")
  21. @smproperty.xml("""
  22. <OutputPort index="0" name='Transform_Z_0.01'/>
  23. <OutputPort index="1" name='Landmass'/>
  24. <OutputPort index="2" name='Water'/>
  25. <OutputPort index="3" name='Velocity_Field'/>
  26. <OutputPort index="4" name='Seed'/>
  27. <OutputPort index="5" name='ParticleTracer' />
  28. <OutputPort index="6" name='Pathlines' />
  29. <Hints>
  30. <ShowInMenu category="pyParaOcean Filters"/>
  31. </Hints>
  32. """)
  33.  
  34. class Pathlines(VTKPythonAlgorithmBase):
  35. # the rest of the code here is unchanged.
  36. def __init__(self):
  37. print("---__init__ Called---")
  38. from vtkmodules.vtkCommonExecutionModel import vtkTrivialProducer
  39.  
  40. self.z1 = 70
  41. self.z2 = 72
  42.  
  43. self.salinity_variable_name = "so"
  44. self.temperature_variable_name = "thetao"
  45.  
  46. #variables set through UI elements
  47. self.inputVectorStr= "" # default value of the input vector field for vector weighted seeding
  48. self.inputScalarStr= "" # default value of scalar
  49. self.seedingStrategy = 1 # default value of seeding strategy
  50. self.numSeeds = 300 # default number of seeds
  51. self.calculated = False
  52.  
  53. self.start_time_set = False
  54.  
  55. self.particle_tracer = vtk.vtkParticleTracer()
  56. self.particle_tracer.SetStaticSeeds(False)
  57. self.particle_tracer.SetIgnorePipelineTime(False)
  58.  
  59.  
  60. self.seeds = vtkPolyData()
  61. self.seed_producer = vtkTrivialProducer()
  62. self.seeded = False
  63.  
  64. super().__init__(nInputPorts = 1, nOutputPorts = 7)
  65.  
  66. # textbox to enter flow field vector
  67. @smproperty.stringvector(name ="Enter Flow Vector Components (comma separated)",
  68. default_values='uo,vo')
  69. def SetVectorField(self, name):
  70. self.inputVectorStr = name.split(",")
  71. self.Modified()
  72.  
  73. # dropdown to select seeding strategy
  74. @smproperty.xml("""
  75. <IntVectorProperty
  76. name="Seeding Strategy"
  77. command="SetSeedingStrategy"
  78. number_of_elements="1"
  79. default_values="1">
  80. <EnumerationDomain name="enum">
  81. <Entry value="1" text="Uniform"/>
  82. <Entry value="2" text="Chosen Scalar's Magnitude"/>
  83. <Entry value="3" text="Curl"/>
  84. <Entry value="4" text="Vorticity"/>
  85. <Entry value="5" text="Okubo-Weiss"/>
  86. <Entry value="6" text="Flow Magnitude"/>
  87.  
  88. </EnumerationDomain>
  89. <Documentation>
  90. This property indicates which seeding strategy will be used.
  91. </Documentation>
  92. </IntVectorProperty>
  93. """)
  94. def SetSeedingStrategy(self, n):
  95. self.seedingStrategy = n
  96. self.Modified()
  97.  
  98. def get_slice(self, pdi, depth):
  99. plane = vtk.vtkPlane()
  100. plane.SetOrigin(0, 0, depth)
  101. plane.SetNormal(0, 0, 1)
  102. cutter = vtk.vtkCutter()
  103. cutter.SetCutFunction(plane)
  104. cutter.SetInputData(pdi)
  105. cutter.Update()
  106.  
  107. return cutter
  108.  
  109. @smproperty.stringvector(name = "Salinity Variable Name", default_values = "so")
  110. def SetSalinityName(self, x):
  111. self.salinity_variable_name = x
  112. self.Modified()
  113.  
  114. @smproperty.stringvector(name = "Temperature Variable Name", default_values = "thetao")
  115. def SetTemperatureName(self, x):
  116. self.temperature_variable_name = x
  117. self.Modified()
  118.  
  119. def Tranform_z(self, data):
  120. t = vtkTransform()
  121. t.Scale(1, 1, 0.01)
  122.  
  123. tf = vtkTransformFilter()
  124. tf.SetTransform(t)
  125. tf.SetInputData(data)
  126. tf.Update()
  127.  
  128. return tf
  129.  
  130. # textbox to enter number of seeding points
  131. @smproperty.intvector(name ="Number of seeds", default_values=5000)
  132. def SetNumSeeds(self, n):
  133. self.numSeeds = n
  134. self.Modified()
  135.  
  136. #TODO: deprecate this helper function
  137. def point_generate(self,input0,id,n):
  138. ip = vtkPolyData()
  139. ps1 = vtk.vtkPointSource()
  140. pt = input0.GetPoint(int(id))
  141. ps1.SetCenter(pt)
  142. ps1.SetNumberOfPoints(n)
  143. ps1.SetRadius(0.4)
  144. ps1.Update()
  145. ip.ShallowCopy(ps1.GetOutput())
  146. return ip
  147.  
  148. def FillOutputPortInformation(self, port, info):
  149. print("---FillOutputPortInformation Called---")
  150. if port in [0, 1, 2, 3]:
  151. info.Set(vtk.vtkDataObject.DATA_TYPE_NAME(), "vtkUnstructuredGrid")
  152. elif port in [4, 5, 6]:
  153. info.Set(vtk.vtkDataObject.DATA_TYPE_NAME(), "vtkPolyData")
  154.  
  155. return 1
  156.  
  157. def RequestInformation(self, request, inInfo, outInfo):
  158. return 1
  159.  
  160. def RequestUpdateExtent(self, request, inInfo, outInfo):
  161. return 1
  162.  
  163. def RequestDataObject(self, request, inInfo, outInfo):
  164. print("---RequestDataObject Called---")
  165. from vtkmodules.vtkCommonDataModel import vtkUnstructuredGrid, vtkRectilinearGrid
  166.  
  167. for i in range(self.GetNumberOfOutputPorts()):
  168. if i in [0, 1, 2, 3]:
  169. outInfo.GetInformationObject(i).Set(vtk.vtkDataObject.DATA_OBJECT(), vtkUnstructuredGrid())
  170. elif i in [4, 5, 6]:
  171. outInfo.GetInformationObject(i).Set(vtk.vtkDataObject.DATA_OBJECT(), vtkPolyData())
  172.  
  173. return 1
  174.  
  175. def RequestData(self, request, inInfo, outInfo):
  176. print("---RequestData Called---")
  177. from vtkmodules.vtkCommonDataModel import vtkRectilinearGrid
  178. from vtkmodules.vtkCommonDataModel import vtkUnstructuredGrid, vtkDataSet
  179. from vtkmodules.vtkFiltersCore import vtkThreshold
  180. from vtkmodules.vtkFiltersGeneral import vtkTransformFilter, vtkRectilinearGridToPointSet
  181. from vtkmodules.vtkFiltersFlowPaths import vtkParticleTracer
  182. from vtkmodules.vtkCommonExecutionModel import vtkStreamingDemandDrivenPipeline as SDSP
  183. from vtkmodules.vtkCommonExecutionModel import vtkTrivialProducer
  184.  
  185. # Get input/output objects
  186. input_data = vtkRectilinearGrid.GetData(inInfo[0])
  187. output_transformed_z = vtkUnstructuredGrid.GetData(outInfo, 0)
  188. output_data_land = vtkUnstructuredGrid.GetData(outInfo, 1)
  189. output_data_water = vtkUnstructuredGrid.GetData(outInfo, 2)
  190. output_velocity_field = dsa.WrapDataObject(vtkDataSet.GetData(outInfo, 3))
  191. output_seeds = vtkPolyData.GetData(outInfo, 4)
  192. particle_tracer_object = vtkPolyData.GetData(outInfo, 5)
  193. pathlines = vtkPolyData.GetData(outInfo, 6)
  194.  
  195. if output_data_land is None or output_data_water is None or output_transformed_z is None:
  196. raise RuntimeError("Output object is None. ParaView did not allocate it.")
  197.  
  198. # Convert to vtkPointSet using vtkRectilinearGridToPointSet
  199. converter = vtkRectilinearGridToPointSet()
  200. converter.SetInputData(input_data)
  201. converter.Update()
  202. pointset = converter.GetOutput()
  203.  
  204. # Apply Transform
  205. transform = vtkTransform()
  206. transform.Scale(1.0, 1.0, 0.01) # example translation
  207.  
  208. transform_filter = vtkTransformFilter()
  209. transform_filter.SetTransform(transform)
  210. transform_filter.SetInputData(pointset)
  211. transform_filter.Update()
  212. transformed = transform_filter.GetOutput()
  213.  
  214. # Apply Threshold
  215. threshold_land = vtkThreshold()
  216. threshold_land.SetInputData(transformed)
  217. threshold_land.SetInputArrayToProcess(0, 0, 0, 0, self.salinity_variable_name)
  218. threshold_land.SetLowerThreshold(-100)
  219. threshold_land.SetUpperThreshold(0.0)
  220. threshold_land.SetThresholdFunction(threshold_land.THRESHOLD_BETWEEN)
  221. threshold_land.Update()
  222.  
  223. threshold_water = vtkThreshold()
  224. threshold_water.SetInputData(transformed)
  225. threshold_water.SetInputArrayToProcess(0, 0, 0, 0, self.salinity_variable_name)
  226. threshold_water.SetLowerThreshold(0.0)
  227. threshold_water.SetUpperThreshold(100.0)
  228. threshold_water.SetThresholdFunction(threshold_water.THRESHOLD_BETWEEN)
  229. threshold_water.Update()
  230.  
  231. # Output
  232. output_data_land.ShallowCopy(threshold_land.GetOutput())
  233. output_data_water.ShallowCopy(threshold_water.GetOutput())
  234.  
  235. input_data_np = dsa.WrapDataObject(vtkDataSet.GetData(inInfo[0]))
  236. output_data_water_np = dsa.WrapDataObject(output_data_water)
  237.  
  238. # Put vector field in velocity
  239. vec = []
  240. for i in self.inputVectorStr:
  241. vec.append(output_data_water_np.PointData[i])
  242. try:
  243. velocity = algs.make_vector(vec[0],vec[1],vec[2])
  244. except:
  245. velocity = algs.make_vector(vec[0],vec[1],np.zeros(len(vec[0])))
  246.  
  247. #print("velocity type: ", type(velocity))
  248. #print("velocity: ", velocity)
  249.  
  250. # Compute jacobian
  251. J = algs.gradient(velocity) # Jacobian
  252. ux, uy, uz = J[:,0,0], J[:,0,1], J[:,0,2]
  253. vx, vy, vz = J[:,1,0], J[:,1,1], J[:,1,2]
  254. wx, wy, wz = J[:,2,0], J[:,2,1], J[:,2,2]
  255.  
  256. # Compute seeding distribution
  257. criterion = None
  258. if self.seedingStrategy == 1:
  259. criterion = None
  260. elif self.seedingStrategy == 2:
  261. scalarMag = algs.abs(output_data_water_np.PointData[self.inputScalarStr])
  262. criterion = scalarMag
  263. elif self.seedingStrategy == 3:
  264. curlSqrMag = (wy - vz)**2 + (uz - wx)**2 + (vx - uy)**2
  265. criterion = curlSqrMag
  266. elif self.seedingStrategy == 4:
  267. vorSqrMag = (vx - uy)**2
  268. criterion = vorSqrMag
  269. elif self.seedingStrategy == 5:
  270. okuboWeiss = np.abs((ux - vy)**2 + (vx + uy)**2 - (vx - uy)**2)
  271. criterion = okuboWeiss
  272. elif self.seedingStrategy == 6:
  273. flowMag = algs.mag(velocity)
  274. criterion = flowMag
  275. else:
  276. criterion = None
  277.  
  278. # Mask unviable points
  279. validPointIds = np.where((vx - uy)**2 > 0)[0]
  280. if criterion is not None:
  281. criterion = criterion[validPointIds]
  282.  
  283. # Sample from set of valid points in domain
  284. seedPointIds = random.choices(validPointIds, weights=criterion, k=self.numSeeds)
  285. newSeedPoints = vtkAppendPolyData()
  286.  
  287. for idval in seedPointIds:
  288. newSeedPoints.AddInputData(self.point_generate(output_data_water_np, idval, 1))
  289. newSeedPoints.Update()
  290.  
  291. # Populate output ports for velocity and seeds
  292. output_velocity_field.ShallowCopy(output_data_water_np.VTKObject)
  293. output_velocity_field.PointData.append(velocity, "velocity")
  294.  
  295. output_seeds.ShallowCopy(newSeedPoints.GetOutput())
  296. if not self.seeded:
  297. self.seeds.ShallowCopy(newSeedPoints.GetOutput())
  298.  
  299. # -------------------------------------------------------------------- ParticleTracer --------------------------------------------------------------------
  300. print("----------------------------------")
  301. # Get input info object from main input (index 0)
  302. input_info = self.GetExecutive().GetInputInformation(0, 0)
  303. output_info = output_velocity_field.GetInformation()
  304. output_info_seeds = self.seeds.GetInformation()
  305.  
  306. # Forward TIME_STEPS
  307. if input_info.Has(SDSP.TIME_STEPS()):
  308. time_steps = input_info.Get(SDSP.TIME_STEPS())
  309. output_info.Set(SDSP.TIME_STEPS(), time_steps, len(time_steps))
  310. output_info_seeds.Set(SDSP.TIME_STEPS(), time_steps, len(time_steps))
  311.  
  312. # Forward TIME_RANGE
  313. if input_info.Has(SDSP.TIME_RANGE()):
  314. time_range = input_info.Get(SDSP.TIME_RANGE())
  315. output_info.Set(SDSP.TIME_RANGE(), time_range, 2)
  316. output_info_seeds.Set(SDSP.TIME_RANGE(), time_range, 2)
  317.  
  318. # Forward UPDATE_TIME_STEP (optional, sometimes needed)
  319. if input_info.Has(SDSP.UPDATE_TIME_STEP()):
  320. current_time = input_info.Get(SDSP.UPDATE_TIME_STEP())
  321. output_info.Set(SDSP.UPDATE_TIME_STEP(), current_time)
  322. output_info.Set(vtk.vtkDataObject.DATA_TIME_STEP(), current_time)
  323.  
  324. output_info_seeds.Set(SDSP.UPDATE_TIME_STEP(), current_time)
  325. output_info_seeds.Set(vtk.vtkDataObject.DATA_TIME_STEP(), current_time)
  326.  
  327. print("------")
  328. if not self.seeded:
  329. self.seeded = True
  330. self.seed_producer.SetOutput(self.seeds)
  331. print("seed_producer.SetOutput")
  332. print("------")
  333.  
  334. print("------")
  335. if output_velocity_field.GetInformation().Has(SDSP.TIME_STEPS()):
  336. time_steps = output_velocity_field.GetInformation().Get(SDSP.TIME_STEPS())
  337. print("TIME_STEPS in output_velocity_field:", time_steps)
  338. else:
  339. print("No TIME_STEPS in output_velocity_field")
  340.  
  341. if output_velocity_field.GetInformation().Has(SDSP.TIME_RANGE()):
  342. time_steps = output_velocity_field.GetInformation().Get(SDSP.TIME_RANGE())
  343. print("TIME_RANGE in output_velocity_field:", time_steps)
  344. else:
  345. print("No TIME_RANGE in output_velocity_field")
  346.  
  347. if output_velocity_field.GetInformation().Has(SDSP.UPDATE_TIME_STEP()):
  348. time_steps = output_velocity_field.GetInformation().Get(SDSP.UPDATE_TIME_STEP())
  349. print("UPDATE_TIME_STEP in output_velocity_field:", time_steps)
  350. else:
  351. print("No UPDATE_TIME_STEP in output_velocity_field")
  352.  
  353. if output_velocity_field.GetInformation().Has(vtk.vtkDataObject.DATA_TIME_STEP()):
  354. time_steps = output_velocity_field.GetInformation().Get(vtk.vtkDataObject.DATA_TIME_STEP())
  355. print("DATA_TIME_STEP in output_velocity_field:", time_steps)
  356. else:
  357. print("No DATA_TIME_STEP in output_velocity_field")
  358.  
  359. velocity_field_producer = vtkTrivialProducer()
  360. velocity_field_producer.SetOutput(output_velocity_field.VTKObject)
  361. velocity_field_producer.Update()
  362.  
  363. info = output_velocity_field.VTKObject.GetInformation()
  364.  
  365. time = info.Get(SDSP.TIME_STEPS())
  366. velocity_field_producer.GetOutputInformation(0).Set(SDSP.TIME_STEPS(), time, len(time_steps))
  367. self.seed_producer.GetOutputInformation(0).Set(SDSP.TIME_STEPS(), time, len(time_steps))
  368.  
  369. time = info.Get(SDSP.TIME_RANGE())
  370. velocity_field_producer.GetOutputInformation(0).Set(SDSP.TIME_RANGE(), time, 2)
  371. self.seed_producer.GetOutputInformation(0).Set(SDSP.TIME_RANGE(), time, len(time_steps))
  372.  
  373. time = info.Get(SDSP.UPDATE_TIME_STEP())
  374. velocity_field_producer.GetOutputInformation(0).Set(SDSP.UPDATE_TIME_STEP(), time)
  375. self.seed_producer.GetOutputInformation(0).Set(SDSP.UPDATE_TIME_STEP(), time)
  376. velocity_field_producer.GetOutputInformation(0).Set(vtk.vtkDataObject.DATA_TIME_STEP(), current_time)
  377. self.seed_producer.GetOutputInformation(0).Set(vtk.vtkDataObject.DATA_TIME_STEP(), time)
  378. # ------
  379.  
  380. print("------")
  381. velocity_field_producer_output = velocity_field_producer.GetOutputDataObject(0)
  382.  
  383. point_data = velocity_field_producer_output.GetPointData()
  384. velocity_array = point_data.GetArray("velocity")
  385. print(f"velocity_field_producer_output: {velocity_array}:")
  386.  
  387.  
  388.  
  389. print("------")
  390. if self.seeds.GetInformation().Has(SDSP.TIME_STEPS()):
  391. time_steps = self.seeds.GetInformation().Get(SDSP.TIME_STEPS())
  392. print("TIME_STEPS in self.seeds:", time_steps)
  393. else:
  394. print("No TIME_STEPS in self.seeds")
  395.  
  396. if self.seeds.GetInformation().Has(SDSP.TIME_RANGE()):
  397. time_steps = self.seeds.GetInformation().Get(SDSP.TIME_RANGE())
  398. print("TIME_RANGE in self.seeds:", time_steps)
  399. else:
  400. print("No TIME_RANGE in output_velocity_field")
  401.  
  402. if self.seeds.GetInformation().Has(SDSP.UPDATE_TIME_STEP()):
  403. time_steps = self.seeds.GetInformation().Get(SDSP.UPDATE_TIME_STEP())
  404. print("UPDATE_TIME_STEP in self.seeds:", time_steps)
  405. else:
  406. print("No UPDATE_TIME_STEP in self.seeds")
  407.  
  408. if self.seeds.GetInformation().Has(vtk.vtkDataObject.DATA_TIME_STEP()):
  409. time_steps = self.seeds.GetInformation().Get(vtk.vtkDataObject.DATA_TIME_STEP())
  410. print("DATA_TIME_STEP in self.seeds:", time_steps)
  411. else:
  412. print("No DATA_TIME_STEP in self.seeds")
  413.  
  414. print("------")
  415. time_steps = 0
  416. if velocity_field_producer.GetOutputInformation(0).Has(SDSP.TIME_STEPS()):
  417. time_steps = velocity_field_producer.GetOutputInformation(0).Get(SDSP.TIME_STEPS())
  418. print("TIME_STEPS in velocity_field_producer: ", time_steps)
  419. else:
  420. print("No TIME_STEPS in velocity_field_producer")
  421.  
  422. if velocity_field_producer.GetOutputInformation(0).Has(SDSP.TIME_RANGE()):
  423. time_steps = velocity_field_producer.GetOutputInformation(0).Get(SDSP.TIME_RANGE())
  424. print("TIME_RANGE in velocity_field_producer: ", time_steps)
  425. else:
  426. print("No TIME_RANGE in velocity_field_producer")
  427.  
  428. if velocity_field_producer.GetOutputInformation(0).Has(SDSP.UPDATE_TIME_STEP()):
  429. time_steps = velocity_field_producer.GetOutputInformation(0).Get(SDSP.UPDATE_TIME_STEP())
  430. print("UPDATE_TIME_STEP in velocity_field_producer: ", time_steps)
  431. else:
  432. print("No UPDATE_TIME_STEP in velocity_field_producer")
  433.  
  434. if velocity_field_producer.GetOutputInformation(0).Has(vtk.vtkDataObject.DATA_TIME_STEP()):
  435. time_steps = velocity_field_producer.GetOutputInformation(0).Get(vtk.vtkDataObject.DATA_TIME_STEP())
  436. print("DATA_TIME_STEP in velocity_field_producer: ", time_steps)
  437. else:
  438. print("No DATA_TIME_STEP in velocity_field_producer")
  439.  
  440.  
  441. print("------")
  442. time_steps = 0
  443. if self.seed_producer.GetOutputInformation(0).Has(SDSP.TIME_STEPS()):
  444. time_steps = self.seed_producer.GetOutputInformation(0).Get(SDSP.TIME_STEPS())
  445. print("TIME_STEPS in self.seed_producer: ", time_steps)
  446. else:
  447. print("No TIME_STEPS in self.seed_producer")
  448.  
  449. if self.seed_producer.GetOutputInformation(0).Has(SDSP.TIME_RANGE()):
  450. time_steps = self.seed_producer.GetOutputInformation(0).Get(SDSP.TIME_RANGE())
  451. print("TIME_RANGE in self.seed_producer: ", time_steps)
  452. else:
  453. print("No TIME_RANGE in self.seed_producer")
  454.  
  455. if self.seed_producer.GetOutputInformation(0).Has(SDSP.UPDATE_TIME_STEP()):
  456. time_steps = self.seed_producer.GetOutputInformation(0).Get(SDSP.UPDATE_TIME_STEP())
  457. print("UPDATE_TIME_STEP in self.seed_producer: ", time_steps)
  458. else:
  459. print("No UPDATE_TIME_STEP in self.seed_producer")
  460.  
  461. if self.seed_producer.GetOutputInformation(0).Has(vtk.vtkDataObject.DATA_TIME_STEP()):
  462. time_steps = self.seed_producer.GetOutputInformation(0).Get(vtk.vtkDataObject.DATA_TIME_STEP())
  463. print("DATA_TIME_STEP in self.seed_producer: ", time_steps)
  464. else:
  465. print("No DATA_TIME_STEP in self.seed_producer")
  466.  
  467. print("------")
  468.  
  469. num_points = self.seeds.GetNumberOfPoints()
  470. print(f"Number of seed points: {num_points}")
  471.  
  472. pd = output_velocity_field.VTKObject.GetPointData()
  473. velocity_array = pd.GetArray("velocity")
  474. print("------")
  475. if velocity_array:
  476. pd.SetVectors(velocity_array)
  477. print("Set 'velocity' as active vectors.")
  478. else:
  479. print("No array named 'velocity' found!")
  480.  
  481.  
  482.  
  483.  
  484. self.particle_tracer.SetInputConnection(0, velocity_field_producer.GetOutputPort())
  485.  
  486. self.particle_tracer.SetInputConnection(1, self.seed_producer.GetOutputPort())
  487.  
  488.  
  489.  
  490.  
  491. velPortInfo = outInfo.GetInformationObject(3)
  492.  
  493.  
  494. tracerExec = self.particle_tracer.GetExecutive()
  495. inVelInfo = tracerExec.GetInputInformation(0, 0)
  496. input_info_tracer_1 = tracerExec.GetInputInformation(1, 0)
  497.  
  498.  
  499. if output_velocity_field.GetInformation().Has(SDSP.TIME_STEPS()):
  500. ts = output_velocity_field.GetInformation().Get(SDSP.TIME_STEPS())
  501.  
  502. inVelInfo.Set(SDSP.TIME_STEPS(), ts, len(ts))
  503. input_info_tracer_1.Set(SDSP.TIME_STEPS(), ts, len(ts))
  504.  
  505. if output_velocity_field.GetInformation().Has(SDSP.TIME_RANGE()):
  506. tr = output_velocity_field.GetInformation().Get(SDSP.TIME_RANGE())
  507.  
  508. inVelInfo.Set(SDSP.TIME_RANGE(), tr, 2)
  509. input_info_tracer_1.Set(SDSP.TIME_RANGE(), tr, 2)
  510.  
  511. if output_velocity_field.GetInformation().Has(SDSP.UPDATE_TIME_STEP()):
  512. ut = output_velocity_field.GetInformation().Get(SDSP.UPDATE_TIME_STEP())
  513.  
  514. inVelInfo.Set(SDSP.UPDATE_TIME_STEP(), ut)
  515. inVelInfo.Set(vtk.vtkDataObject.DATA_TIME_STEP(), ut)
  516.  
  517. input_info_tracer_1.Set(SDSP.UPDATE_TIME_STEP(), ut)
  518. input_info_tracer_1.Set(vtk.vtkDataObject.DATA_TIME_STEP(), ut)
  519.  
  520. self.particle_tracer.SetInputArrayToProcess(0, 0, 0, vtk.vtkDataObject.FIELD_ASSOCIATION_POINTS, "velocity")
  521.  
  522. print("------")
  523. input_info_tracer_0 = self.particle_tracer.GetInputInformation(0, 0)
  524. input_info_tracer_1 = self.particle_tracer.GetInputInformation(1, 0)
  525.  
  526. if input_info_tracer_0.Has(SDSP.TIME_STEPS()) and input_info_tracer_1.Has(SDSP.TIME_STEPS()):
  527. time_steps = input_info_tracer_0.Get(SDSP.TIME_STEPS())
  528. print("TIME_STEPS in input_info_tracer: ", time_steps)
  529. else:
  530. print("No TIME_STEP in input_info_tracer")
  531.  
  532. if input_info_tracer_0.Has(SDSP.TIME_RANGE()) and input_info_tracer_1.Has(SDSP.TIME_RANGE()):
  533. time_steps = input_info_tracer_0.Get(SDSP.TIME_RANGE())
  534. print("TIME_RANGE in input_info_tracer: ", time_steps)
  535. else:
  536. print("No TIME_RANGE in input_info_tracer")
  537.  
  538. if input_info_tracer_0.Has(SDSP.UPDATE_TIME_STEP()) and input_info_tracer_1.Has(SDSP.UPDATE_TIME_STEP()):
  539. time_steps = input_info_tracer_0.Get(SDSP.UPDATE_TIME_STEP())
  540. print("UPDATE_TIME_STEP in input_info_tracer: ", time_steps)
  541. else:
  542. print("No UPDATE_TIME_STEP in input_info_tracer")
  543.  
  544. if input_info_tracer_0.Has(vtk.vtkDataObject.DATA_TIME_STEP()) and input_info_tracer_1.Has(vtk.vtkDataObject.DATA_TIME_STEP()):
  545. time_steps = input_info_tracer_0.Get(vtk.vtkDataObject.DATA_TIME_STEP())
  546. print("DATA_TIME_STEP in input_info_tracer: ", time_steps)
  547. else:
  548. print("No DATA_TIME_STEP in input_info_tracer")
  549.  
  550. print("---")
  551. print("REQUEST_UPDATE_TIME: ", SDSP.REQUEST_UPDATE_TIME())
  552. # suppose 'alg' is any vtkAlgorithm (e.g. your reader or trivial producer)
  553. # info = tracer.GetExecutive().GetOutputInformation(0)
  554. # or if you have the producer: info = producer.GetOutputInformation(0)
  555.  
  556. # request time T
  557. # info.Set(SDSP.UPDATE_TIME_STEP(), time_steps)
  558.  
  559. # output_info = particle_tracer.GetInformation()
  560.  
  561. # time_steps = input_info_tracer.Get(SDSP.TIME_STEPS())
  562. # output_info.Set(SDSP.TIME_STEPS(), time_steps, len(time_steps))
  563.  
  564. # time_range = input_info_tracer.Get(SDSP.TIME_RANGE())
  565. # output_info.Set(SDSP.TIME_RANGE(), time_range, 2)
  566.  
  567. # update_time_step = input_info_tracer.Get(SDSP.UPDATE_TIME_STEP())
  568. # output_info.Set(SDSP.UPDATE_TIME_STEP(), update_time_step)
  569.  
  570. # data_time_step = input_info_tracer.Get(vtk.vtkDataObject.DATA_TIME_STEP())
  571. # output_info.Set(vtk.vtkDataObject.DATA_TIME_STEP(), data_time_step)
  572.  
  573. # self.particle_tracer.SetInitialIntegrationStep(0.1) # Adjust as necessary
  574. # self.particle_tracer.SetMaximumNumberOfSteps(1000)
  575. # print("------")
  576. # input_data = self.particle_tracer.GetInputDataObject(0, 0)
  577. # if input_data is not None:
  578. # vectors = input_data.GetPointData().GetVectors()
  579. # if vectors:
  580. # print("Velocity field found:", vectors.GetName())
  581. # print("Number of vectors:", vectors.GetNumberOfTuples())
  582. # else:
  583. # print("No velocity vectors found in point data.")
  584. # else:
  585. # print("No input dataset connected.")
  586.  
  587. # This works if you've connected your seeds to port 1 (seeds port)
  588. seeds = self.particle_tracer.GetInputDataObject(1, 0)
  589. if seeds is not None:
  590. print("Seeds present.")
  591. print("Number of seed points:", seeds.GetNumberOfPoints())
  592. else:
  593. print("No seed points connected.")
  594. print("------")
  595. print("BEFORE Update(): Num Particles:", self.particle_tracer.GetOutput().GetNumberOfPoints())
  596. print("BEFORE Update(): Current Time(Port-0):", self.particle_tracer.GetInputInformation(0, 0).Get(SDSP.UPDATE_TIME_STEP()))
  597. print("BEFORE Update(): Current Time(Port-1):", self.particle_tracer.GetInputInformation(1, 0).Get(SDSP.UPDATE_TIME_STEP()))
  598. print("BEFORE Update(): Current Time Object(Port-0):", self.particle_tracer.GetInputInformation(0, 0))
  599. print("BEFORE Update(): Current Time Object(Port-1):", self.particle_tracer.GetInputInformation(1, 0))
  600.  
  601. rk4 = vtk.vtkRungeKutta4()
  602. #self.particle_tracer.SetIntegrator(rk4)
  603. #self.particle_tracer.SetIntegratorType(vtkParticleTracer.RUNGE_KUTTA2)
  604. self.particle_tracer.SetInterpolatorType(vtkParticleTracer.INTERPOLATOR_WITH_DATASET_POINT_LOCATOR)
  605.  
  606. if not self.start_time_set:
  607. self.start_time_set = True
  608. self.particle_tracer.SetStartTime(self.particle_tracer.GetInputInformation(0, 0).Get(SDSP.UPDATE_TIME_STEP()))
  609.  
  610. print("------")
  611. print("GetIntegrator: ", self.particle_tracer.GetIntegrator())
  612. print("GetIntegratorType: ", self.particle_tracer.GetIntegratorType())
  613. print("GetStartTime: ", self.particle_tracer.GetStartTime())
  614. print("GetIgnorePipelineTime: ", self.particle_tracer.GetIgnorePipelineTime())
  615.  
  616. self.particle_tracer.InitializeObjectBase()
  617. self.particle_tracer.Modified()
  618. # Pass time to the tracer
  619. #tracer.SetCurrentTime(input_info_tracer_0.Has(SDSP.UPDATE_TIME_STEP()))
  620. self.particle_tracer.Update()
  621.  
  622. # current_time_step = self.particle_tracer.GetInputInformation(0, 0).Get(SDSP.UPDATE_TIME_STEP())
  623. # current_time_step = current_time_step + 24
  624. # next_time_step = current_time_step # Or however you want to calculate the next time step
  625. # self.particle_tracer.GetInputInformation(0, 0).Set(SDSP.UPDATE_TIME_STEP(), next_time_step)
  626.  
  627. print("------")
  628. print("AFTER Update(): Num Particles:", self.particle_tracer.GetOutput().GetNumberOfPoints())
  629. print("AFTER Update(): Current Time(Port-0):", self.particle_tracer.GetInputInformation(0, 0).Get(SDSP.UPDATE_TIME_STEP()))
  630. print("AFTER Update(): Current Time(Port-1):", self.particle_tracer.GetInputInformation(1, 0).Get(SDSP.UPDATE_TIME_STEP()))
  631. print("AFTER Update(): Current Time Object(Port-0):", self.particle_tracer.GetInputInformation(0, 0))
  632. print("AFTER Update(): Current Time Object(Port-1):", self.particle_tracer.GetInputInformation(1, 0))
  633. print("------")
  634.  
  635. # Check number of particles
  636. print(f"Particles after update: {self.particle_tracer.GetOutput().GetNumberOfPoints()}")
  637.  
  638. # Ensure ParticleAge is being correctly updated
  639. particle_age_array = self.particle_tracer.GetOutput().GetPointData().GetArray('ParticleAge')
  640. if particle_age_array:
  641. print("ParticleAge values after update:", particle_age_array)
  642.  
  643. # for i in range(self.particle_tracer.GetOutput().GetNumberOfPoints()):
  644. # age = particle_age_array.GetValue(i)
  645. # print("age: ", age)
  646.  
  647. print("------")
  648. output = self.particle_tracer.GetOutput()
  649. if output is not None and output.GetNumberOfPoints() > 0:
  650. print("Particles generated:", output.GetNumberOfPoints())
  651. else:
  652. print("No particles generated.")
  653.  
  654. particle_tracer_object.ShallowCopy(self.particle_tracer.GetOutput())
  655. print("------")
  656. # After Copying Tracer's output
  657. if output_info.Has(SDSP.TIME_STEPS()):
  658. time_steps = output_info.Get(SDSP.TIME_STEPS())
  659. print("TIME_STEP in particle_tracer: ", time_steps)
  660. else:
  661. print("No TIME_STEP in particle_tracer")
  662.  
  663. if output_info.Has(SDSP.TIME_RANGE()):
  664. time_steps = output_info.Get(SDSP.TIME_RANGE())
  665. print("TIME_RANGE in particle_tracer: ", time_steps)
  666. else:
  667. print("No TIME_RANGE in particle_tracer")
  668.  
  669. if output_info.Has(SDSP.UPDATE_TIME_STEP()):
  670. time_steps = output_info.Get(SDSP.UPDATE_TIME_STEP())
  671. print("UPDATE_TIME_STEP in particle_tracer: ", time_steps)
  672. else:
  673. print("No UPDATE_TIME_STEP in particle_tracer")
  674.  
  675. if output_info.Has(vtk.vtkDataObject.DATA_TIME_STEP()):
  676. time_steps = output_info.Get(vtk.vtkDataObject.DATA_TIME_STEP())
  677. print("DATA_TIME_STEP in particle_tracer: ", time_steps)
  678. else:
  679. print("No DATA_TIME_STEP in particle_tracer")
  680.  
  681.  
  682.  
  683. # ------------------------------------------------------------------------------------------------------------------------------------------------------
  684.  
  685. # ---------------------------------------------------- TemporalParticlesToPathlines --------------------------------------------------------------------
  686.  
  687. pathline_filter = vtk.vtkTemporalPathLineFilter()
  688.  
  689. #pathline_filter.SetInputConnection(self.particle_tracer.GetOutputPort())
  690. pathline_filter.SetInputData(particle_tracer_object)
  691.  
  692. # pathline_filter.SetMaskPoints(1)
  693. # pathline_filter.SetMaxTrackLength(200)
  694. # pathline_filter.SetIdChannelArray("ParticleId")
  695. # 3. Set parameters (optional, depending on what you're tracking)
  696. #pathline_filter.SetMaskArray("Mask") # optional: if you want to mask/filter some particles
  697. pathline_filter.SetIdChannelArray("ParticleId") # required to track individual particles
  698. #pathline_filter.SetTimeChannelArray("IntegrationTime") # required to trace particle over time
  699. #pathline_filter.SetKeepFinalPoints(1) # optional: keeps the final points at the end of pathlines
  700. pathline_filter.SetMaxTrackLength(2500) # -1 means no limit
  701. pathline_filter.SetMaxStepDistance(0.0, 0.0, 0.0) # 0 means automatic
  702.  
  703. pathline_filter.Update()
  704.  
  705. print("Num Lines:", pathline_filter.GetOutput().GetNumberOfLines())
  706.  
  707. # 4. Update the filter to produce the output
  708. # pathline_filter.Update()
  709.  
  710. pathlines.ShallowCopy(pathline_filter.GetOutput())
  711.  
  712. # Ensure connectivity is available
  713. # pathlines.GetCellData().AddArray(pathline_filter.GetOutput().GetCellData().GetArray("ParticleId"))
  714.  
  715.  
  716.  
  717. # ------------------------------------------------------------------------------------------------------------------------------------------------------
  718.  
  719. # 'Hide' Water
  720. empty = vtkUnstructuredGrid()
  721. output_data_water.ShallowCopy(empty)
  722. #output_velocity_field.ShallowCopy(empty)
  723.  
  724. # 'Hide' Seeds
  725. #empty = vtkPolyData()
  726. #output_seeds.ShallowCopy(empty)
  727.  
  728. return 1
  729.  
  730.  
  731.  
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement
Public Pastes
Advertisement
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!