Subject: [Viper-dev] visualizing a a 2D cut plane

From: Johan Hake <[email-masked]>
Date: Mon, 11 Aug 2008 11:25:04 +0200

Hello!

I am trying to visualize a 2D cut plane of my 3D scalar field. I suppose this 
is not included in the Viper features, yet? :)

How can I do this most easy? Write the data to vtk format, through viper, and 
then use for example paraview to do the cut plane?

Johan

---

From: "Martin Sandve Alnæs" <[email-masked]>
Date: Mon, 11 Aug 2008 11:32:51 +0200

What's wrong with

f = dolfin.Function(...)
v = File("foo.pvd")
v << f

?

2008/8/11 Johan Hake <[email-masked]>:
> _______________________________________________
> Viper-dev mailing list
> [email-masked]
> http://fenics.org/mailman/listinfo/viper-dev
>

-- 
Martin

---

From: Johan Hake <[email-masked]>
Date: Mon, 11 Aug 2008 11:43:28 +0200

> v << f

If I had a dolfin function this would be the prefered case. But in the pycc 
framwork this is not the case. It probably make sense to create a dolfin 
function from the numpy array and then save this to file using the above 
syntax.

Johan

> > I am trying to visualize a 2D cut plane of my 3D scalar field. I suppose
> > this is not included in the Viper features, yet? :)
> >
> > How can I do this most easy? Write the data to vtk format, through viper,
> > and then use for example paraview to do the cut plane?
> >
> >
> > Johan
> >

---

From: Ola Skavhaug <[email-masked]>
Date: Mon, 11 Aug 2008 13:03:58 +0200

> is not included in the Viper features, yet? :)

This is not supported yet, but I wouldn't rule out the possibility of it
getting supported in the future.

Ola
 
> then use for example paraview to do the cut plane?

You could use the write_vtk() method of the viper object. Or, use the DOLFIN
writer.

Ola

---

From: Johan Hake <[email-masked]>
Date: Mon, 11 Aug 2008 13:49:33 +0200

> getting supported in the future.

I took a breif look and the filter "should" not be too difficult to create.
  
   plane = vtk.vtkPlane()
   plane.SetOrigin(0,0,0)
   plane.SetNormal(0,1,0)
   cutter = vtk.vtkCutter()
   cutter.SetCutFunction(plane)
   cutter.SetInput(vtkUnstructuredGrid)

But I have not enough vtk experience to know if this is right or to know how 
to utelize the filter further. Maybe the cutter in my example can be used 
instead of the vtkUnstructuredGrid data set that is used in viper?

Johan

> You could use the write_vtk() method of the viper object. Or, use the
> DOLFIN writer.
>
> Ola
>
> > Johan
> >

---

From: Johan Hake <[email-masked]>
Date: Mon, 11 Aug 2008 23:35:20 +0200

Hello!

I implemented this feature. See the patch.

It is only possible to use it for a scalar field in 3 dimension, which make 
sense I suppose, please correct me! It is configurable by setting the cut 
plane by an origin point and a normal vector defining a cut plane. This is 
done by the kwargs [cut_plane_origin,cut_plane_normal].

The camera position is aligned with the normal vector. I couldn't find out how 
to set the view up vector. This is needed so the plot is positioned 
correctly. Any clever ideas? Otherwise we need to hand that vector from the 
interface too.

It does not interact correctly with other functionality, e.g., write_png. This 
should be easy to correct, by introducing a more intuitive filter hierarchy.

Johan
# HG changeset patch
# User "Johan Hake <[email-masked]>"
# Date 1218490464 -7200
# Node ID 4cd1f1094400857b2ae6ede871cddadd0afd7e55
# Parent  7c7412da6f10617d3d4d5d491f7cd5318d278a7a
Implemented an interface to plot cut plane with viper

diff -r 7c7412da6f10 -r 4cd1f1094400 src/viper/viper.py
--- a/src/viper/viper.py	Thu Jul 31 15:10:46 2008 +0200
+++ b/src/viper/viper.py	Mon Aug 11 23:34:24 2008 +0200
@@ -61,7 +61,10 @@
 __init__(x=x, mesh=mymesh, displacement=d):
     plot scalar of vector field x over mymesh and use d as a displacement vector field
 __init__(coordinates=coords, cells=elements):
-    plot mesh given in the two numpy arrays coords and elements 
+    plot mesh given in the two numpy arrays coords and elements
+__init__(..., cut_plane_origin = origin, cut_plane_normal = normal)
+    plot a cut plane of a 3D scalar plot. Set the origin and normal to define the cut plane.
+    The normal defaults to (0,0,1)
 
 Additional parameters: vmin, vmax, lutfile, wireframe, celltype, title
         """
@@ -111,10 +114,32 @@
         self.displacement = self._resolve_displacement(displacement)
 
         self.frame = kwargs.get("frame", None)
+
+        # Extract cut plane information
+        cut_plane_origin = kwargs.get("cut_plane_origin", None)
+        cut_plane_normal = kwargs.get("cut_plane_normal", (0,0,1))
+        
         if len(x.shape) == 1:
             if self.warpscalar and minmax > 0:
                 self.vtkgrid = self.warp_scalar(self.x, minmax)
-            (self.iren, self.renWin, self.ren) = self.simple_plotter(self.vtkgrid, vmin, vmax, wireframe)
+            
+            # If plot a cut plane
+            if not cut_plane_origin is None:
+                cut_plane_data = self._cut_plane_filter(cut_plane_origin, cut_plane_normal)
+                # Create the plotter using the cut plane data
+                (self.iren, self.renWin, self.ren) = self.simple_plotter(cut_plane_data, vmin, vmax, wireframe)
+                # Set the view to the cut_plane_origin and align the camera
+                # with the cut plane normal
+                self.ren.GetActiveCamera().SetFocalPoint(cut_plane_origin)
+                self.ren.GetActiveCamera().SetPosition(cut_plane_normal)
+                # FIXME! Should we compute this automaticly or include it
+                #        as another parameter?
+                #self.ren.GetActiveCamera().SetViewUp(NOT_SET) 
+                self.ren.ResetCamera()
+                self.ren.GetActiveCamera().Dolly(1.5)
+            else:
+                # Create the plotter using the original data
+                (self.iren, self.renWin, self.ren) = self.simple_plotter(self.vtkgrid, vmin, vmax, wireframe)
         else:
             self.mode="vector"
             (self.iren, self.renWin, self.ren) = self.vector_plotter(self.vtkgrid, self.mesh.coordinates(), x, vmin, vmax, wireframe)
@@ -868,6 +893,32 @@
             act.GetProperty().SetColor(1, 0, 0)
             self.ren.AddActor(act)
 
+    def _cut_plane_filter(self, origo, normal):
+        """ Return a cut plan data """
+        if not self.mesh.cells().shape[1] == 4:
+            raise RuntimeError, "Can only cut scalar 3 dimensional plots"
+        
+        if not (isinstance(origo, (numpy.ndarray, tuple, list)) and 
+                len(origo) == 3):
+            raise RuntimeError, "Provide a tuple, list or numpy array of length 3 for cut_plane_origo"
+        
+        if not (isinstance(normal, (numpy.ndarray, tuple, list)) and 
+                len(normal) ==3):
+            raise RuntimeError, "Provide a tuple, list or numpy array of length 3 for cut_plane_normal"
+
+        # Define the cut plane
+        plane = vtk.vtkPlane()
+        plane.SetOrigin(origo)
+        plane.SetNormal(normal)
+
+        # Initialize the Cutter filter
+        cutter = vtk.vtkCutter()
+        cutter.SetCutFunction(plane)
+        cutter.SetInput(self.vtkgrid)
+        cut_data = cutter.GetOutput()
+        return cut_data
+        
+
     def set_sphere_opacity(self, val):
         """Change the sphere opacity (range [0,1])."""
         if val < 0.0: val = 0.0