calculator.py 8.1 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220
  1. r"""This module is used by vtkPythonCalculator. It encapsulates the logic
  2. implemented by the vtkPythonCalculator to operate on datasets to compute
  3. derived quantities.
  4. """
  5. try:
  6. import numpy as np
  7. except ImportError:
  8. raise RuntimeError("'numpy' module is not found. numpy is needed for this functionality to work. ")
  9. import paraview
  10. import vtkmodules.numpy_interface.dataset_adapter as dsa
  11. from vtkmodules.numpy_interface.algorithms import *
  12. # -- this will import vtkMultiProcessController and vtkMPI4PyCommunicator
  13. from paraview.vtk import vtkDataObject, vtkDoubleArray, vtkSelectionNode, vtkSelection, vtkStreamingDemandDrivenPipeline
  14. from paraview.modules import vtkPVVTKExtensionsFiltersPython
  15. import sys
  16. if sys.version_info >= (3,):
  17. xrange = range
  18. def get_arrays(attribs, controller=None):
  19. """Returns a 'dict' referring to arrays in dsa.DataSetAttributes or
  20. dsa.CompositeDataSetAttributes instance.
  21. When running in parallel, this method will ensure that arraynames are
  22. reduced across all ranks and for any arrays missing on the local process, a
  23. NoneArray will be added to the returned dictionary. This ensures that
  24. expressions evaluate without issues due to missing arrays on certain ranks.
  25. """
  26. if not isinstance(attribs, dsa.DataSetAttributes) and \
  27. not isinstance(attribs, dsa.CompositeDataSetAttributes):
  28. raise ValueError(
  29. "Argument must be DataSetAttributes or CompositeDataSetAttributes.")
  30. arrays = dict()
  31. for key in attribs.keys():
  32. varname = paraview.make_name_valid(key)
  33. arrays[varname] = attribs[key]
  34. # If running in parallel, ensure that the arrays are synced up so that
  35. # missing arrays get NoneArray assigned to them avoiding any unnecessary
  36. # errors when evaluating expressions.
  37. if controller is None and vtkMultiProcessController is not None:
  38. controller = vtkMultiProcessController.GetGlobalController()
  39. if controller and controller.IsA("vtkMPIController") and controller.GetNumberOfProcesses() > 1:
  40. from mpi4py import MPI
  41. comm = vtkMPI4PyCommunicator.ConvertToPython(controller.GetCommunicator())
  42. rank = comm.Get_rank()
  43. # reduce the array names across processes to ensure arrays missing on
  44. # certain ranks are handled correctly.
  45. arraynames = list(arrays) # get keys from the arrays as a list.
  46. # gather to root and then broadcast
  47. # I couldn't get Allgather/Allreduce to work properly with strings.
  48. gathered_names = comm.gather(arraynames, root=0)
  49. # gathered_names is a list of lists.
  50. if rank == 0:
  51. result = set()
  52. for alist in gathered_names:
  53. for val in alist: result.add(val)
  54. gathered_names = [x for x in result]
  55. arraynames = comm.bcast(gathered_names, root=0)
  56. for name in arraynames:
  57. if name not in arrays:
  58. arrays[name] = dsa.NoneArray
  59. return arrays
  60. def pointIsNear(locations, distance, inputs):
  61. array = vtkDoubleArray()
  62. array.SetNumberOfComponents(3)
  63. array.SetNumberOfTuples(len(locations))
  64. for i in range(len(locations)):
  65. array.SetTuple(i, locations[i])
  66. node = vtkSelectionNode()
  67. node.SetFieldType(vtkSelectionNode.POINT)
  68. node.SetContentType(vtkSelectionNode.LOCATIONS)
  69. node.GetProperties().Set(vtkSelectionNode.EPSILON(), distance)
  70. node.SetSelectionList(array)
  71. from paraview.vtk.vtkFiltersExtraction import vtkLocationSelector
  72. selector = vtkLocationSelector()
  73. selector.SetInsidednessArrayName("vtkInsidedness")
  74. selector.Initialize(node)
  75. inputDO = inputs[0].VTKObject
  76. outputDO = inputDO.NewInstance()
  77. outputDO.CopyStructure(inputDO)
  78. output = dsa.WrapDataObject(outputDO)
  79. if outputDO.IsA('vtkCompositeDataSet'):
  80. it = inputDO.NewIterator()
  81. it.InitTraversal()
  82. while not it.IsDoneWithTraversal():
  83. outputDO.SetDataSet(it, inputDO.GetDataSet(it).NewInstance())
  84. it.GoToNextItem()
  85. selector.Execute(inputDO, outputDO)
  86. return output.PointData.GetArray('vtkInsidedness')
  87. def cellContainsPoint(inputs, locations):
  88. array = vtkDoubleArray()
  89. array.SetNumberOfComponents(3)
  90. array.SetNumberOfTuples(len(locations))
  91. for i in range(len(locations)):
  92. array.SetTuple(i, locations[i])
  93. node = vtkSelectionNode()
  94. node.SetFieldType(vtkSelectionNode.CELL)
  95. node.SetContentType(vtkSelectionNode.LOCATIONS)
  96. node.SetSelectionList(array)
  97. from paraview.vtk.vtkFiltersExtraction import vtkLocationSelector
  98. selector = vtkLocationSelector()
  99. selector.SetInsidednessArrayName("vtkInsidedness")
  100. selector.Initialize(node)
  101. inputDO = inputs[0].VTKObject
  102. outputDO = inputDO.NewInstance()
  103. outputDO.CopyStructure(inputDO)
  104. output = dsa.WrapDataObject(outputDO)
  105. if outputDO.IsA('vtkCompositeDataSet'):
  106. it = inputDO.NewIterator()
  107. it.InitTraversal()
  108. while not it.IsDoneWithTraversal():
  109. outputDO.SetDataSet(it, inputDO.GetDataSet(it).NewInstance())
  110. it.GoToNextItem()
  111. selector.Execute(inputDO, outputDO)
  112. return output.CellData.GetArray('vtkInsidedness')
  113. def compute(inputs, expression, ns=None):
  114. # build the locals environment used to eval the expression.
  115. mylocals = dict()
  116. if ns:
  117. mylocals.update(ns)
  118. mylocals["inputs"] = inputs
  119. try:
  120. mylocals["points"] = inputs[0].Points
  121. except AttributeError:
  122. pass
  123. finalRet = None
  124. for subEx in expression.split(' and '):
  125. retVal = eval(subEx, globals(), mylocals)
  126. if finalRet is None:
  127. finalRet = retVal
  128. else:
  129. finalRet = dsa.VTKArray([a & b for a, b in zip(finalRet, retVal)])
  130. return finalRet
  131. def get_data_time(self, do, ininfo):
  132. dinfo = do.GetInformation()
  133. if dinfo and dinfo.Has(do.DATA_TIME_STEP()):
  134. t = dinfo.Get(do.DATA_TIME_STEP())
  135. else:
  136. t = None
  137. key = vtkStreamingDemandDrivenPipeline.TIME_STEPS()
  138. t_index = None
  139. if ininfo.Has(key):
  140. tsteps = [ininfo.Get(key, x) for x in range(ininfo.Length(key))]
  141. try:
  142. t_index = tsteps.index(t)
  143. except ValueError:
  144. pass
  145. return (t, t_index)
  146. def execute(self, expression):
  147. """
  148. **Internal Method**
  149. Called by vtkPythonCalculator in its RequestData(...) method. This is not
  150. intended for use externally except from within
  151. vtkPythonCalculator::RequestData(...).
  152. """
  153. # Add inputs.
  154. inputs = []
  155. for index in range(self.GetNumberOfInputConnections(0)):
  156. # wrap all input data objects using vtkmodules.numpy_interface.dataset_adapter
  157. wdo_input = dsa.WrapDataObject(self.GetInputDataObject(0, index))
  158. t, t_index = get_data_time(self, wdo_input.VTKObject, self.GetInputInformation(0, index))
  159. wdo_input.time_value = wdo_input.t_value = t
  160. wdo_input.time_index = wdo_input.t_index = t_index
  161. inputs.append(wdo_input)
  162. # Setup output.
  163. output = dsa.WrapDataObject(self.GetOutputDataObject(0))
  164. if self.GetCopyArrays():
  165. for attr in range(vtkDataObject.NUMBER_OF_ATTRIBUTE_TYPES):
  166. if inputs[0].HasAttributes(attr):
  167. inputAttribute = inputs[0].GetAttributes(attr)
  168. output.GetAttributes(attr).PassData(inputAttribute)
  169. # get a dictionary for arrays in the dataset attributes. We pass that
  170. # as the variables in the eval namespace for compute.
  171. variables = get_arrays(inputs[0].GetAttributes(self.GetArrayAssociation()))
  172. variables.update({"time_value": inputs[0].time_value,
  173. "t_value": inputs[0].t_value,
  174. "time_index": inputs[0].time_index,
  175. "t_index": inputs[0].t_index})
  176. retVal = compute(inputs, expression, ns=variables)
  177. if retVal is not None:
  178. if hasattr(retVal, "Association") and retVal.Association is not None:
  179. output.GetAttributes(retVal.Association).append(retVal, self.GetArrayName())
  180. else:
  181. # if somehow the association was removed we
  182. # fall back to the input array association
  183. output.GetAttributes(self.GetArrayAssociation()).append(retVal, self.GetArrayName())