.. include:: ../../macros.hrst
.. include:: ../../abbreviations.hrst

.. _chapter:PythonAndBatchPythonCalculatorProgrammableSourceAndFilter:

Python & Batch: Python Calculator, Programmable Source & Filter
###############################################################

Introduction
============

|ParaView| has three filters that give a user access to python math
functions as well as the underlying VTK library. These are the Python
Calculator, the Programmable Source and Filter. This tutorial explains both.

Python Calculator
=================

The Python Calculator allows a user to apply calculations that are
available in Python. These include such functions as get volume of
cells, get area of cells, get the cross product, dot product, curl, etc.
The whole formula must fit on one line. Lets go over the Python
Calculator with an example.

-  Lets create a new point variable, loaded with 5.

   -  Open can.ex2.
   -  Turn all variables on.
   -  **Apply**.
   -  **Filters → Alphabetical → Python Calculator**.
   -  Change Expression to **5**.
   -  Set **Array Association** to **Point Data**
   -  Change **Array Name** to **Calculated Variable**.
   -  **Apply**.
   -  Paint by **Calculated Variable**.

-  Lets create a variable equal to two times displacement.

   -  Change Expression to **DISPL*2**.
   -  **Apply**.
   -  Go to last timestep.
   -  Rescale to Data Range.
   -  Note: A more complete way to access displacement data. We explicitly pull
      **DISPL** from the first input to the filter. Don't forget that the first
      input is [0], the second is [1], etc.

   -  Change Expression to **(inputs[0].PointData['DISPL']) * 2**.
   -  **Apply**.
   -  Rescale to Data Range.

-  Here is how to multiply a vector by a global variable. Lets multiply
   DISPL by the timestep (TMSTEP)

   -  Change expression to **inputs[0].FieldData[‘NSTEPS’][time_index] * DISPL**.
   -  **Apply**.
   -  Rescale to Data Range.
   -  To recap, If a function requires a variable input, use the string above. If the
      function needs the input mesh, use inputs[].

   .. admonition:: **Did you know?**
      :class: tip

      Python uses square brackets for arrays, and parenthesis
      to group parts of your formula. Also, to designate a string (i.e.,
      variable name), use either single or double quotation marks.

-  Lets create a variable to contains the cell volume.

   -  Change expression to **volume(inputs[0])**.
   -  **Apply**.
   -  Rescale to Data Range.

-  Numerous functions can be combined in one expression. For instance, a
   nonsensical expression would be to take a sin of a divergence of a
   curl. Here is the expression:

   -  Change expression to **sin(divergence(curl('ACCL')))**. Apply.
      Rescale to Data Range.

Interesting are functions available through the Python Calculator:

-  area(dataset)
-  aspect(dataset)
-  cos(array)
-  cross(X,Y) where X and Y are two 3D vector arrays
-  curl(array)
-  divergence(array)
-  dot(a1,a2)
-  eigenvalue and eigenvector(array)
-  gradient(array)
-  max(array)
-  mean(array)
-  min(array)
-  norm(array)
-  sin(array)
-  strain(array)
-  volume(array)
-  vorticity(array)

The complete list can be found in :numref:`sec:PythonCalculator` of the User's Guide.


Programmable Filter
===================

The Programmable Filter is used to modify data in the pipeline using
Python. An example of a Programmable Filter dividing a variable by 2:

-  Divide **ACCL** by 2

   -  Open can.ex2.
   -  Turn all variables on.
   -  **Apply**.
   -  Select **Filters → Alphabetical → Programmable Filter**.
   -  Leave Output Data
   -  Set Type as **Same as Input**.
   -  Enter the following into the Script window:

      .. code-block:: python

         input0 = inputs[0]
         dataArray = input0.PointData["ACCL"] / 2.0
         output.PointData.append(dataArray, "ACCL_half")

   -  Set **Coloring** by **ACCL_half**

-  Subtract two datasets from each other. Use two instances of
   disk_out_ref.ex2 as two datasets, subtract **GaMe3** from **AsH3**.

   -  Open disk_out_ref.ex2.
   -  Turn all variables on.
   -  **Apply**.
   -  Open disk_out_ref.ex2.
   -  Turn all variables on.
   -  **Apply**.
   -  Highlight both datasets. Use the key with the mouse to select more
      than one input.
   -  Select **Filters → Alphabetical → Programmable Filter**.
   -  Leave Output Data
   -  Set Type as **Same as Input**.
   -  Enter the following into the Script window:

      .. code-block:: python

         v_0 = inputs[0].PointData['AsH3']
         v_1 = inputs[1].PointData['GaMe3']
         output.PointData.append(v_1 - v_0, 'difference')

   -  Set **Coloring** to **difference**

-  Create a tensor

   .. code-block:: python

      from paraview.vtk.numpy_interface import dataset_adapter as dsa
      import numpy
      def make_tensor(xx,yy,zz, xy, yz, xz):

            t = numpy.vstack([xx,yy,zz,xy, yz,
            xz]).transpose().view(dsa.VTKArray)
            t.DataSet = xx.DataSet
            t.Association = xx.Association
            return t

      xx = inputs[0].PointData["sigma_xx"]
      yy = inputs[0].PointData["sigma_yy"]
      zz = inputs[0].PointData["sigma_zz"]
      xy = inputs[0].PointData["sigma_xy"]
      yz = inputs[0].PointData["sigma_yz"]
      xz = inputs[0].PointData["sigma_xz"]
      output.PointData.append(make_tensor(xx,yy,zz,xy,yz,xz), "tensor")


-  Subtract two timesteps from each other

   -  Load can.ex2
   -  Select **Filters → Alphabetical → Force Time**.
   -  Set the timestep to the one whose displacement you want to be zero.
   -  The output of this filter is the can.ex2 dataset "frozen" at the chosen timestep. It will not change as
      you advance through time.
   -  Select both can.ex2 (using the **CTRL** key) and the **ForceTime1** sources in
      the Pipeline Browser. Note that the original can.ex2 source will
      update as the timestep changes, but the **ForceTime1** source will not change.
   -  Select **Filters → Alphabetical → Python Calculator**.
   -  Substract the displacement in the "frozen" dataset from the current timestep in can.ex2.
   -  Set the expression to **inputs[0].PointData['DISPL'] - inputs[1].PointData['DISPL']**
   -  The order of inputs into the **Python Calculator** is not well
      defined, so you may need to swap the indices in inputs[0] and
      inputs[1] to get the correct sign on the result, but this should work.
   -  (If desired,) Finally, add a Plot Selection over Time filter. This
      filter will run over all time steps, subtracting the data from the
      current timestep from the data in the "frozen" timestep produced by the Force Time
      filter, and plot the result in a graph. The first timestep should have a value of 0.

-  Read eigenvector, calculate an eigenvalue, and place it into three variables

   .. code-block:: python

      # ParaView Programmable Filter script.  ParaView 5.0.1.
      #
      # This code will read in an eigenvector, calculate an
      #    eigenvalue, and place it into three variables.
      #
      # Be sure to correct the input variables below.  Also, note that
      #   the code uses ZX, not XZ.
      #
      # This code only works with any multiblock or vtk
      #   datasets (including ones with only one block - i.e.,
      #   Exodus datasets as input).
      #
      # Usage:  Run the Programmable Filter.
      #   Cut and paste this file in the section named "Script".
      #   Leave "Output Data Set Type" as "Same as Input".
      #   Click Apply button
      #
      # Written by Jeff Mauldin and Alan Scott
      #

      import numpy as np

      def process_composite_dataset(input0):
        # Pick up input arrays
        xxar = input0.CellData["EPSXX"]
        xyar = input0.CellData["EPSXY"]
        zxar = input0.CellData["EPSZX"]
        yyar = input0.CellData["EPSYY"]
        yzar = input0.CellData["EPSYZ"]
        zzar = input0.CellData["EPSZZ"]

        #print `xxar`
        #print len(xxar.Arrays)

        # Set output arrays to same type as input array.
        # Do a multiply to make sure we don't just have a
        # pointer to the original.
        outarray0 = xxar*0.5
        outarray1 = xxar*0.5
        outarray2 = xxar*0.5


        # Run a for loop over all blocks
        numsubarrays = len(xxar.Arrays)
        for ii in range(0, numsubarrays):
          # pick up input arrays for each block.
          xxarsub = xxar.Arrays[ii]
          xyarsub = xyar.Arrays[ii]
          zxarsub = zxar.Arrays[ii]
          yyarsub = yyar.Arrays[ii]
          yzarsub = yzar.Arrays[ii]
          zzarsub = zzar.Arrays[ii]

          #print `xxarsub`

          # Transpose and calculate the principle strain.
          strain = np.transpose(
              np.array(
                [ [xxarsub, xyarsub, zxarsub],
                  [xyarsub, yyarsub, yzarsub],
                  [zxarsub, yzarsub, zzarsub] ] ),
                    (2,0,1))

          principal_strain = np.linalg.eigvalsh(strain)

          # Move principle strain to temp output arrays for this block
          outarray0.Arrays[ii] = principal_strain[:,0]
          outarray1.Arrays[ii] = principal_strain[:,1]
          outarray2.Arrays[ii] = principal_strain[:,2]

        #ps0 = principal_strain[:,0]
        #print "ps0 len: " + str(len(ps0))

        # Finally, move the temp arrays to output arrays
        output.CellData.append(outarray0, "principal_strain_0")
        output.CellData.append(outarray1, "principal_strain_1")
        output.CellData.append(outarray2, "principal_strain_2")


      def process_unstructured_dataset(input0):
        # Pick up input arrays
        xxar = input0.CellData["EPSXX"]
        xyar = input0.CellData["EPSXY"]
        zxar = input0.CellData["EPSZX"]
        yyar = input0.CellData["EPSYY"]
        yzar = input0.CellData["EPSYZ"]
        zzar = input0.CellData["EPSZZ"]

        #print `xxar`
        #print len(xxar.Arrays)

        # Transpose and calculate the principle strain.
        strain = np.transpose(
              np.array(
                [ [xxar, xyar, zxar],
                  [xyar, yyar, yzar],
                  [zxar, yzar, zzar] ] ),
                    (2,0,1))

        principal_strain = np.linalg.eigvalsh(strain)

        #ps0 = principal_strain[:,0]
        #print "ps0 len: " + str(len(ps0))

        # Finally, move the temp arrays to output arrays
        output.CellData.append(principal_strain[:,0],
           "principal_strain_0")
        output.CellData.append(principal_strain[:,1],
           "principal_strain_1")
        output.CellData.append(principal_strain[:,2],
           "principal_strain_2")

      input0 = inputs[0]

      if input0.IsA("vtkCompositeDataSet"):
        process_composite_dataset(input0)
      elif input0.IsA("vtkUnstructuredGrid"):
        process_unstructured_dataset(input0)
      else:
        print "Bad dataset type for this script"


Important note: Since the **Programmable Source** and **Programmable Filter** work at the
server level, paraview.simple cannot be loaded or used.

Mode details about the **Programmable Filter** can be found in :numref:`chapter:PythonProgrammableFilter` of the
Reference Manual.

Additionally, this is list of interesting Blog Posts related to the **Programmable Filter**:

-  https://blog.kitware.com/improved-vtk-numpy-integration/
-  https://blog.kitware.com/improved-vtk-numpy-integration-part-2/
-  https://blog.kitware.com/improved-vtk-numpy-integration-part-3/
-  https://blog.kitware.com/mpi4py-and-vtk/
-  https://blog.kitware.com/improved-vtk-numpy-integration-part-4/
-  https://blog.kitware.com/improved-vtk-numpy-integration-part-5/
