[arvindpj]

Dear All,

Based on the feature metiond in the paraview blog, I am trying to use the matplotlib library to plot the data filtered through "PlotOverLine" filter. I believe this is of type vtkMultiBlockDataSet with 1 child (vtkPolyData)

I am having issues with SetAttributeArrayStatus(…) for an array in vtkMultiBlockDataSet.

Here is the script, I used:

def setup_data(view):
  from paraview.numpy_support import vtk_to_numpy
  for i in xrange(view.GetNumberOfVisibleDataObjects()):
    dataObject = view.GetVisibleDataObjectForSetup(i).GetBlock(0)   #Accessing the point data
    pressure = dataObject.GetPointData().GetArray("p")
    print vtk_to_numpy(pressure)  # prints the pressure array
    view.DisableAllAttributeArrays()
    view.SetAttributeArrayStatus(i, vtkDataObject.POINT, "p", 1)  # code breaks here! Could not pass the vtkDataObject.POINT array

def render(view, width, height):
  from paraview.numpy_support import vtk_to_numpy
  from paraview import python_view
  figure = python_view.matplotlib_figure(width, height)
  for i in xrange(view.GetNumberOfVisibleDataObjects()):
    dataObject = view.GetVisibleDataObjectForRendering(i) .GetBlock(0)   
    pressure = dataObject.GetPointData().GetArray("p")   
    print vtk_to_numpy(pressure)

Any help is much appreciated.
Thanks,
Jay

[arvindpj]

import vtk
def setup_data(view):
  from paraview.numpy_support import vtk_to_numpy
  for i in xrange(view.GetNumberOfVisibleDataObjects()):
    dataObject = view.GetVisibleDataObjectForSetup(i).GetBlock(0)   #Accessing the point data
    pressure = dataObject.GetPointData().GetArray("P")
    print vtk_to_numpy(pressure)  # prints the pressure array
    view.DisableAllAttributeArrays()
    view.SetAttributeArrayStatus(i, vtk.vtkDataObject.POINT, "P", 1)  
    
def render(view, width, height):
  from paraview.numpy_support import vtk_to_numpy
  from paraview import python_view
  figure = python_view.matplotlib_figure(width, height)
  for i in xrange(view.GetNumberOfVisibleDataObjects()):
    dataObject = view.GetVisibleDataObjectForRendering(i) .GetBlock(0)   
    pressure = dataObject.GetPointData().GetArray("P")   
    print vtk_to_numpy(pressure)

And this works!

[francois]

Nice !
Please, could you detail your workflow to explain this very interesting feature to the community.
Many thanks

[arvindpj]

One could incorporate matplotlib plots using the new type of view in paraview 4.1 called the "python view".

More details are available here.

Here is a typical script to use inside python view to plot line plots. (bar charts and cluster plots could also be made):

import vtk
def setup_data(view):
  from paraview.numpy_support import vtk_to_numpy
  for i in xrange(view.GetNumberOfVisibleDataObjects()):
    dataObject = view.GetVisibleDataObjectForSetup(i)
    view.DisableAllAttributeArrays()
    #view.EnableAllAttributeArrays()
    view.SetAttributeArrayStatus(i, vtk.vtkDataObject.POINT, "mag(wallShearStress)", 1)
    print "Done.. setup_data!"

def render(view, width, height):
  from paraview.numpy_support import vtk_to_numpy
  from paraview import python_view
  #from matplotlib import style
  #style.use('ggplot') 

  figure = python_view.matplotlib_figure(width, height)
  plt = figure.add_subplot(1,1,1)
  plt.set_xlabel('$Length \quad [m]$', fontsize=24)
  plt.set_ylabel(r'$\tau_{wall}\quad   [m/s^{2}]$  ', fontsize=24)

  for tick in plt.xaxis.get_major_ticks():
    tick.label.set_fontsize(16)
  for tick in plt.yaxis.get_major_ticks():
    tick.label.set_fontsize(16)

  for i in xrange(view.GetNumberOfVisibleDataObjects()):
    dataObject = view.GetVisibleDataObjectForRendering(i)

    tauwall = vtk_to_numpy(dataObject.GetPointData().GetArray("mag(wallShearStress)"))
    
    coord = vtk_to_numpy(dataObject.GetPoints().GetData())
    x,y,z = coord[:,0],coord[:,1],coord[:,2]
      
    plt.plot(y,tauwall,'r-', linewidth=2)
    figure.savefig('wallShear.png')
    print "Done.. render_data!"
  return python_view.figure_to_image(figure)

Here is the script to automate the whole process by saving it a "plotOverIntersection.py" and then can be executed using pvpython.

 pvpython < plotOverIntersection.py

import os
try: paraview.simple
except: from paraview.simple import *
paraview.simple._DisableFirstRenderCameraReset()

script = """

import vtk
def setup_data(view):
  from paraview.numpy_support import vtk_to_numpy
  for i in xrange(view.GetNumberOfVisibleDataObjects()):
    dataObject = view.GetVisibleDataObjectForSetup(i)
    view.DisableAllAttributeArrays()
    #view.EnableAllAttributeArrays()
    view.SetAttributeArrayStatus(i, vtk.vtkDataObject.POINT, "mag(wallShearStress)", 1)
    print "Done.. setup_data!"

def render(view, width, height):
  from paraview.numpy_support import vtk_to_numpy
  from paraview import python_view
  #from matplotlib import style
  #style.use('ggplot') 

  figure = python_view.matplotlib_figure(width, height)
  plt = figure.add_subplot(1,1,1)
  plt.set_xlabel('$Length \quad [m]$', fontsize=24)
  plt.set_ylabel(r'$\tau_{wall}\quad   [m/s^{2}]$  ', fontsize=24)

  for tick in plt.xaxis.get_major_ticks():
    tick.label.set_fontsize(16)
  for tick in plt.yaxis.get_major_ticks():
    tick.label.set_fontsize(16)

  for i in xrange(view.GetNumberOfVisibleDataObjects()):
    dataObject = view.GetVisibleDataObjectForRendering(i)

    tauwall = vtk_to_numpy(dataObject.GetPointData().GetArray("mag(wallShearStress)"))
    
    coord = vtk_to_numpy(dataObject.GetPoints().GetData())
    x,y,z = coord[:,0],coord[:,1],coord[:,2]
      
    plt.plot(y,tauwall,'r-', linewidth=2)
    figure.savefig('wallShear.png')
    print "Done.. render_data!"
  return python_view.figure_to_image(figure)

"""




# create a new 'OpenFOAMReader'
case_foam = OpenFOAMReader(FileName='case.foam')   
case_foam.MeshRegions = ['wall_pipe']  #MESH REGION, USUALLY A SURFACE INCASE OF PLOTTING OVER INTERSECTION
case_foam.CellArrays = ['U', 'mag(U)', 'mag(wallShearStress)', 'p', 'wallShearStress'] # DATA TO LOAD


SetActiveSource(case_foam)


# get animation scene
animationScene1 = GetAnimationScene()

# update animation scene based on data timesteps
#animationScene1.UpdateAnimationUsingDataTimeSteps()

animationScene1.GoToLast()  #SKIP TO LAST TIME STEP


# get active view
renderView1 = GetActiveViewOrCreate('RenderView')

# Create the first P-O-L, set the extremes of the line and then commit to render
plotOnIntersectionCurves1 = PlotOnIntersectionCurves(Input=case_foam)
plotOnIntersectionCurves1.SliceType = 'Plane'

# Properties modified on plotOnIntersectionCurves1.SliceType
plotOnIntersectionCurves1.SliceType.Origin = [0, 0, 0]
plotOnIntersectionCurves1.SliceType.Normal = [1, 0, 0]

# create a new 'Merge Blocks'
mergeBlocks1 = MergeBlocks(Input=plotOnIntersectionCurves1)

view = CreateView("PythonView")
view.Script = script


Show(mergeBlocks1, view)

pythonView1 = GetActiveViewOrCreate('PythonView')
pythonView1.ViewSize = [ 1200, 1200 ]


source = mergeBlocks1
writer = CreateWriter("plotOverIntersection1.csv", source)
writer.FieldAssociation = "Points" # or "Cells"
writer.UpdatePipeline()

#SaveScreenshot('wallShear.png', magnification=2, quality=100, view=pythonView1)
#WriteImage("wallShear2.png",2)

RenderAllViews()
del writer
del plotOnIntersectionCurves1
del case_foam

Do let me know, if you need any additional details.
Cheers :-)
-Jay