Merge lp:~fluidity-core/fluidity/speedup-vtudiff into lp:fluidity

[Stephan Kramer]

This speeds up the vtudiff script run on vtus with multiple fields, by only setting up the VTK probe once instead of once per field. It does not change anything in the way fields are probed.

There is also a little bit of cleaning up:
* removing the option to call vtktools as a main: this really doesn't make sense and should be in a separate script
* removing the "field manipulation" methods. These were basically vectorial operations on fields but implemented in a rather inefficient way. I see no reason to not simply use numpy vector operations for that, so instead of:

vtu.SubFieldFromField("fieldName", array, "newFieldName")

just do:

field = vtu.GetField("fieldName")
vtu.AddField("newFieldName", field-array)

[lp:~fluidity-core/fluidity/speedup-vtudiff]

4308. By Stephan Kramer on 2014-01-21

For vtudiff: Also compute diff of cell-based fields (only works if both meshes are the same)

[Jon Hill (jon-hill)]

review: Approve

[Stephan Kramer (s-kramer)]

Thanks Jon, will merge once I have a green buildbot.

[lp:~fluidity-core/fluidity/speedup-vtudiff]

4309. By Stephan Kramer on 2014-02-17

Maintain the same sign convention in vtudiff as before:

With "vtudiff A.vtu B.vtu C.vtu" we get fieldC = fieldA-fieldB
Also don't try to subtract vtkGhostLevels field.

4310. By Stephan Kramer on 2014-02-19

Merge in trunk.

=== modified file 'manual/visualisation_and_diagnostics.tex'
--- manual/visualisation_and_diagnostics.tex	2013-04-12 10:31:41 +0000
+++ manual/visualisation_and_diagnostics.tex	2014-02-19 10:25:38 +0000
@@ -350,8 +350,6 @@
   method & arguments & use \\ \hline 
 \lstinline[language=Python]+AddField+ & \lstinline[language=Python]+name, array+ & Adds the values in \lstinline[language=Python]+array+ (the entries of which may have an arbitrary number of components) as a field called \lstinline[language=Python]+name+ \\ \hline
 %
-\lstinline[language=Python]+AddFieldtoField+ & \lstinline[language=Python]+fieldname,+ \lstinline[language=Python]+array,+ \lstinline[language=Python]+newFieldName=None+ & Adds the values in \lstinline[language=Python]+array+ to the field \lstinline[language=Python]+fieldname+. If \lstinline[language=Python]+newFieldName+ is specified then a new field with that name is created with the new values, otherwise the original field is replaced. \\ \hline
-%
 \lstinline[language=Python]+AddScalarField+ & \lstinline[language=Python]+name, array+ & Adds a scalar field called \lstinline[language=Python]+name+ using the values in \lstinline[language=Python]+array+ \\ \hline
 %
 \lstinline[language=Python]+AddVectorField+ & \lstinline[language=Python]+name, array+ & Adds a vector field called \lstinline[language=Python]+name+ using the values in \lstinline[language=Python]+array+ \\ \hline
@@ -365,10 +363,6 @@
 \lstinline[language=Python]+CellDataToPointData+ & & transforms all cell--wise fields to point--wise fields. All existing fields will remain, \\ \hline
 \lstinline[language=Python]+Crop+ & \lstinline[language=Python]+min_x, max_x,+ \lstinline[language=Python]+min_y, max_y+, \lstinline[language=Python]+min_z, max_z+ & Crops the edges defined by the bounding box given by the arguments \\ \hline
 %
-\lstinline[language=Python]+CrossFieldWithField+ &  \lstinline[language=Python]+fieldName,+ \lstinline[language=Python]+array,+ \lstinline[language=Python]+newFieldName=None,+ \lstinline[language=Python]+postMultiply=True+ & Calculates the cross product, ${\bf a} \times {\bf b}$, where ${\bf a}$ is the field \lstinline[language=Python]+fieldName+ and ${\bf b}$ is \lstinline[language=Python]+array+. If \lstinline[language=Python]+postMultiply$\neq$True+ then ${\bf b} \times {\bf a}$ will be calculated. If \lstinline[language=Python]+newFieldName+ is specified  then a new field with that name is created that takes the values of the cross product, otherwise the original field is replaced. \\ \hline  
-%
-\lstinline[language=Python]+DotFieldWithField+ & \lstinline[language=Python]+fieldName,+ \lstinline[language=Python]+array,+ \lstinline[language=Python]+newFieldName+ & Calculates the dot product of the field called \lstinline[language=Python]+fieldName+ and \lstinline[language=Python]+array+. If \lstinline[language=Python]+newFieldName+ is specified then a new field with that name is created that takes the values of the dot product, otherwise the original field is replaced. \\ \hline
-%
 \lstinline[language=Python]+GetCellPoints+ & \lstinline[language=Python]+id+ & Returns an array with the node numbers of the cell (mesh element) number \lstinline[language=Python]+id+. \\ \hline
 %
 \lstinline[language=Python]+GetCellVolume+ & \lstinline[language=Python]+id+ & Returns the volume of the cell (mesh element) with number \lstinline[language=Python]+id+ \\ \hline
@@ -404,18 +398,12 @@
 %
 \lstinline[language=Python]+IntegrateField+ & \lstinline[language=Python]+field+ & Returns the integral of the field called \lstinline[language=Python]+field+ assuming a linear representation on a tetrahedral mesh. \\ \hline
 %
-\lstinline[language=Python]+ManipulateField+ & \lstinline[language=Python]+fieldName,+ \lstinline[language=Python]+manipFunc,+ \lstinline[language=Python]+newFieldName=None+ & Generic field manipulation method. Applies the supplied manipulation function, \lstinline[language=Python]+manipFunc+, to the field called \lstinline[language=Python]+fieldName+. \lstinline[language=Python]+manipFunc+ must have form \lstinline[language=Python]+def manipFunc(field, index):+ $\ldots$ \lstinline[language=Python]+return fieldValAtIndex+. If \lstinline[language=Python]+newFieldName+ is specified then a new field with that name is created that takes the calculated values, otherwise the original field is replaced. \\ \hline
-%
-\lstinline[language=Python]+MatMulFieldWithField+ & \lstinline[language=Python]+fieldName,+ \lstinline[language=Python]+array,+ \lstinline[language=Python]+newFieldName=None,+ \lstinline[language=Python]+postMultiply=True+. & Multiplies two matrices $\bar{\bar{A}}\bar{\bar{B}}$ where $\bar{\bar{A}}$ is the field \lstinline[language=Python]+fieldName+, and $\bar{\bar{B}}$ is \lstinline[language=Python]+array+. If \lstinline[language=Python]+postMultiply+$\neq$\lstinline[language=Python]+True+ then $\bar{\bar{A}}\bar{\bar{B}}$ will be calculated. If \lstinline[language=Python]+newFieldName+ is specified  then a new field with that name is created that takes the values of the product, otherwise the original field is replaced. \\ \hline
-%
- \lstinline[language=Python]+ProbeData+ & \lstinline[language=Python]+coordinates,+ \lstinline[language=Python]+name+ & Returns an array of values of the field called \lstinline[language=Python]+name+ at the positions given in \lstinline[language=Python]+coordinates+. The values are calculated by interpolation of the field to the positions given. \lstinline[language=Python]+coordinates+ can be created using \lstinline[language=Python]+vtktools.arr()+ e.g. \lstinline[language=Python]+coordinates=vtktools.arr([[1,1,1],[1,1,2]])+. \\ \hline
+\lstinline[language=Python]+ProbeData+ & \lstinline[language=Python]+coordinates,+ \lstinline[language=Python]+name+ & Returns an array of values of the field called \lstinline[language=Python]+name+ at the positions given in \lstinline[language=Python]+coordinates+. The values are calculated by interpolation of the field to the positions given. \lstinline[language=Python]+coordinates+ can be created using \lstinline[language=Python]+vtktools.arr()+ e.g. \lstinline[language=Python]+coordinates=vtktools.arr([[1,1,1],[1,1,2]])+. \\ \hline
 %
 \lstinline[language=Python]+RemoveField+ & \lstinline[language=Python]+name+ & Removes the field called \lstinline[language=Python]+name+. \\ \hline 
 %
 \lstinline[language=Python]+StructuredPointProbe+ & \lstinline[language=Python]+nx,+ \lstinline[language=Python]+ny,+ \lstinline[language=Python]+nz,+ \lstinline[language=Python]+bounding_box=None+ & Returns a vtk structured object. \lstinline[language=Python]+nx, ny, nz+ are the number of points in the $x,\, y, \, z$ directions respectively. If \lstinline[language=Python]+bounding_box+ is not specified the bounding box of the domain is calculated automatically. If specified \lstinline[language=Python]+bounding_box = [xmin, xmax, ymin, ymax, zmin, zmax]+. \\ \hline
 %
-\lstinline[language=Python]+SubFieldFromField+ & \lstinline[language=Python]+fieldName,+ \lstinline[language=Python]+array,+ \lstinline[language=Python]+newFieldName=None+ & Subtracts \lstinline[language=Python]+array+ from the field called \lstinline[language=Python]+fieldName+. If \lstinline[language=Python]+newFieldName+ is specified  then a new field with that name is created that takes the values of the product, otherwise the original field is replaced. \\ \hline
-%
 \lstinline[language=Python]+Write+ & \lstinline[language=Python]+filename = []+ & Writes the data to a vtu file. If \lstinline[language=Python]+filename+ is not specified the name of the file originally read in will be used and therefore the input file will be overwritten. \\ \hline
 %
 % 
 
=== modified file 'python/vtktools.py'
--- python/vtktools.py	2012-03-19 21:00:50 +0000
+++ python/vtktools.py	2014-02-19 10:25:38 +0000
@@ -262,64 +262,8 @@
 
   def ProbeData(self, coordinates, name):
     """Interpolate field values at these coordinates."""
-
-    # Initialise locator
-    locator = vtk.vtkPointLocator()
-    locator.SetDataSet(self.ugrid)
-    locator.SetTolerance(10.0)
-    locator.Update()
-
-    # Initialise probe
-    points = vtk.vtkPoints()
-    points.SetDataTypeToDouble()
-    ilen, jlen = coordinates.shape
-    for i in range(ilen):
-      points.InsertNextPoint(coordinates[i][0], coordinates[i][1], coordinates[i][2])
-    polydata = vtk.vtkPolyData()
-    polydata.SetPoints(points)
-    probe = vtk.vtkProbeFilter()
-    probe.SetInput(polydata)
-    probe.SetSource(self.ugrid)
-    probe.Update()
-
-    # Generate a list invalidNodes, containing a map from invalid nodes in the
-    # result to their closest nodes in the input
-    valid_ids = probe.GetValidPoints()
-    valid_loc = 0
-    invalidNodes = []
-    for i in range(ilen):
-      if valid_ids.GetTuple1(valid_loc) == i:
-        valid_loc += 1
-      else:
-        nearest = locator.FindClosestPoint([coordinates[i][0], coordinates[i][1], coordinates[i][2]])
-        invalidNodes.append((i, nearest))
-
-    # Get final updated values
-    pointdata=probe.GetOutput().GetPointData()
-    vtkdata=pointdata.GetArray(name)
-    nc=vtkdata.GetNumberOfComponents()
-    nt=vtkdata.GetNumberOfTuples()
-    array = arr([vtkdata.GetValue(i) for i in range(nt * nc)])
-    
-    # Fix the point data at invalid nodes
-    if len(invalidNodes) > 0:
-      try:
-        oldField = self.ugrid.GetPointData().GetArray(name)
-        components = oldField.GetNumberOfComponents()
-      except:
-        try:
-          oldField = self.ugrid.GetCellData().GetArray(name)
-          components = oldField.GetNumberOfComponents()
-        except:
-          raise Exception("ERROR: couldn't find point or cell field data with name "+name+" in file "+self.filename+".")
-      for invalidNode, nearest in invalidNodes:
-        for comp in range(nc):
-          array[invalidNode * nc + comp] = oldField.GetValue(nearest * nc + comp)
-          
-    valShape = self.GetField(name)[0].shape
-    array.shape = tuple([nt] + list(valShape))
-          
-    return array
+    probe = VTU_Probe(self.ugrid, coordinates)
+    return probe.GetField(name)
 
   def RemoveField(self, name):
     """Removes said field from the unstructured grid."""
@@ -489,93 +433,6 @@
     probe.Update ()
 
     return probe.GetOutput()
-
-  ### Field manipulation methods ###
-
-  def ManipulateField(self, fieldName, manipFunc, newFieldName = None):
-    """
-    Generic field manipulation method. Applies the supplied manipulation function
-    manipFunc to the field fieldName. manipFunc must be a function of the form:
-
-      def manipFunc(field, index):
-        # ...
-        return fieldValAtIndex
-    """
-
-    field = self.GetField(fieldName)
-    if newFieldName is None or fieldName == newFieldName:
-      self.RemoveField(fieldName)
-      newFieldName = fieldName
-
-    field = arr([manipFunc(field, i) for i in range(len(field))])
-    self.AddField(newFieldName, field)
-
-    return
-
-  def AddFieldToField(self, fieldName, array, newFieldName = None):
-    def ManipFunc(field, index):
-      return field[index] + array[index]
-
-    self.ManipulateField(fieldName, ManipFunc, newFieldName)
-
-    return
-
-  def SubFieldFromField(self, fieldName, array, newFieldName = None):
-    def ManipFunc(field, index):
-      return field[index] - array[index]
-
-    self.ManipulateField(fieldName, ManipFunc, newFieldName)
-
-    return
-
-  def DotFieldWithField(self, fieldName, array, newFieldName = None):
-    """
-    Dot product
-    """
-
-    def ManipFunc(field, index):
-      sum = 0.0
-      for i, val in enumerate(field[index]):
-        sum += val * array[index][i]
-        
-      return sum
-
-    self.ManipulateField(fieldName, ManipFunc, newFieldName)
-
-    return
-
-  def CrossFieldWithField(self, fieldName, array, newFieldName = None, postMultiply = True):
-    """
-    Cross product
-    """
-
-    def ManipFunc(field, index):
-      if postMultiply:
-        return numpy.cross(field[index], array[index])
-      else:
-        return numpy.cross(array[index], field[index])
-
-    self.ManipulateField(fieldName, ManipFunc, newFieldName)
-
-    return
-
-  def MatMulFieldWithField(self, fieldName, array, newFieldName = None, postMultiply = True):
-    """
-    Matrix multiplication
-    """
-
-    def ManipFunc(field, index):
-      if postMultiply:
-        return numpy.matrix(field[i]) * numpy.matix(array[i])
-      else:
-        return numpy.matix(array[i]) * numpy.matrix(field[i])
-
-    self.ManipulateField(fieldName, ManipFunc, newFieldName)
-
-    return
-
-  # Default multiplication is dot product
-  MulFieldByField = DotFieldWithField
   
   def GetDerivative(self, name):
     """
@@ -629,6 +486,73 @@
     cdtpd.PassCellDataOn()
     cdtpd.Update()
     self.ugrid=cdtpd.GetUnstructuredGridOutput()
+
+class VTU_Probe(object):
+  """A class that combines a vtkProbeFilter with a list of invalid points (points that it failed to probe
+  where we take the value of the nearest point)"""
+
+  def __init__(self, ugrid, coordinates):
+    # Initialise locator
+    locator = vtk.vtkPointLocator()
+    locator.SetDataSet(ugrid)
+    locator.SetTolerance(10.0)
+    locator.Update()
+
+    # Initialise probe
+    points = vtk.vtkPoints()
+    points.SetDataTypeToDouble()
+    ilen, jlen = coordinates.shape
+    for i in range(ilen):
+      points.InsertNextPoint(coordinates[i][0], coordinates[i][1], coordinates[i][2])
+    polydata = vtk.vtkPolyData()
+    polydata.SetPoints(points)
+    self.probe = vtk.vtkProbeFilter()
+    self.probe.SetInput(polydata)
+    self.probe.SetSource(ugrid)
+    self.probe.Update()
+
+    # Generate a list invalidNodes, containing a map from invalid nodes in the
+    # result to their closest nodes in the input
+    valid_ids = self.probe.GetValidPoints()
+    valid_loc = 0
+    self.invalidNodes = []
+    for i in range(ilen):
+      if valid_ids.GetTuple1(valid_loc) == i:
+        valid_loc += 1
+      else:
+        nearest = locator.FindClosestPoint([coordinates[i][0], coordinates[i][1], coordinates[i][2]])
+        self.invalidNodes.append((i, nearest))
+    self.ugrid = ugrid
+
+  def GetField(self, name):
+    # Get final updated values
+    pointdata = self.probe.GetOutput().GetPointData()
+    vtkdata=pointdata.GetArray(name)
+    nc=vtkdata.GetNumberOfComponents()
+    nt=vtkdata.GetNumberOfTuples()
+    array = arr([vtkdata.GetValue(i) for i in range(nt * nc)])
+    
+    # Fix the point data at invalid nodes
+    if len(self.invalidNodes) > 0:
+      oldField = self.ugrid.GetPointData().GetArray(name)
+      if oldField is None:
+        oldField = self.ugrid.GetCellData().GetArray(name)
+        if oldField is None:
+          raise Exception("ERROR: couldn't find point or cell field data with name "+name+".")
+      components = oldField.GetNumberOfComponents()
+      for invalidNode, nearest in self.invalidNodes:
+        for comp in range(nc):
+          array[invalidNode * nc + comp] = oldField.GetValue(nearest * nc + comp)
+          
+    # this is a copy and paster from vtu.GetField above:
+    if nc==9:
+      return array.reshape(nt,3,3)
+    elif nc==4:
+      return array.reshape(nt,2,2)
+    else:
+      return array.reshape(nt,nc)
+          
+    return array
     
 def VtuMatchLocations(vtu1, vtu2, tolerance = 1.0e-6):
   """
@@ -698,6 +622,8 @@
 
   # If the input vtu point locations match, do not use probe
   useProbe  = not VtuMatchLocations(vtu1, vtu2)
+  if useProbe:
+    probe = VTU_Probe(vtu2.ugrid, vtu1.GetLocations())
 
   # Copy the grid from the first input vtu into the output vtu
   resultVtu.ugrid.DeepCopy(vtu1.ugrid)
@@ -707,65 +633,40 @@
   fieldNames1 = vtu1.GetFieldNames()
   fieldNames2 = vtu2.GetFieldNames()
   for fieldName in fieldNames1:
+    field1 = vtu1.GetField(fieldName)
     if fieldName in fieldNames2:
       if useProbe:
-        field2 = vtu2.ProbeData(vtu1.GetLocations(), fieldName)
+        field2 = probe.GetField(fieldName)
       else:
         field2 = vtu2.GetField(fieldName)
-      resultVtu.SubFieldFromField(fieldName, field2)
+      resultVtu.AddField(fieldName, field1-field2)
     else:
       resultVtu.RemoveField(fieldName)
 
+  # Also look for cell-based fields. This only works if we don't have
+  # to interpolate (both meshes are the same)
+  vtkdata=vtu1.ugrid.GetCellData()
+  fieldNames1 = [vtkdata.GetArrayName(i) for i in range(vtkdata.GetNumberOfArrays())]
+  vtkdata=vtu2.ugrid.GetCellData()
+  fieldNames2 = [vtkdata.GetArrayName(i) for i in range(vtkdata.GetNumberOfArrays())]
+  if useProbe:
+    # meshes are different - we can't interpolate cell-based fields so let's just remove them from the output
+    for fieldName in fieldNames1:
+      if fieldName=='vtkGhostLevels':
+        # this field should just be passed on unchanged
+        continue
+      resultVtu.RemoveField(fieldName)
+  else:
+    # meshes are the same - we can simply subtract
+    for fieldName in fieldNames1:
+      if fieldName=='vtkGhostLevels':
+        # this field should just be passed on unchanged
+        continue
+      elif fieldName in fieldNames2:
+        field1 = vtu1.GetField(fieldName)
+        field2 = vtu2.GetField(fieldName)
+        resultVtu.AddField(fieldName, field1-field2)
+      else:
+        resultVtu.RemoveField(fieldName)
+
   return resultVtu  
-
-def usage():
-  print 'Usage:'
-  print 'COMMAND LINE: vtktools.py [-h] [-p] [-e var1,var2, ...] INPUT_FILENAME'
-  print ''
-  print 'INPUT_FILENAME:'
-  print '           The input file name.'
-  print ''
-  print 'OPTIONS:'
-  print '   -h      Prints this usage message.'
-  print '   -p      Converts the coordinates from xyz to latitude and longitude.'
-  print '   -e      Extracts the data point from the variables provided.'
-
-if __name__ == "__main__":
-  import vtktools
-  import math
-  import getopt
-  
-  optlist, args = getopt.getopt(sys.argv[1:], 'hpe:')
-   
-  v = vtktools.vtu(args[0])
-  
-  # Parse arguments
-  LongLat = False
-  for o,a in optlist:
-    if o == '-h':
-      usage()
-    elif o == '-p':
-      LongLat = True
-    elif o == '-e':
-      scalars = a.strip().split(",") 
-
-  # Project domain if necessary
-  if(LongLat):
-    v.ApplyEarthProjection()
- 
-  # Extract variables
-  if(scalars):
-    npoints = v.ugrid.GetNumberOfPoints ()
-    nvar = len(scalars)
-    for i in range (npoints):
-      (x,y,z) = v.ugrid.GetPoint (i)
-      line = "%lf "%x+"%lf "%y
-      for scalar in scalars:
-        line = line+" %lf"%v.ugrid.GetPointData().GetArray(scalar).GetTuple1(i)
-      print line
-
-
-
-
-
-