Subject: [Viper-dev] nice 1d plotting

From: Dag Lindbo <[email-masked]>
Date: Tue, 15 Apr 2008 10:03:41 +0200

Dear Viper developers!

I have written a simple python function that produce a nice x-versus-y 
plot for (numpy) vectors x and y (using VTK XYPlotActor). This is then 
wrapped in a function that extracts these vectors from a DOLFIN 1d 
function. The plot(x,y) function by itself is useful as a development 
tool (imho).

Pleas take a look at the code attached and give me feedback. I hope to 
insert this into Viper somehow.

/Dag
<?xml version="1.0" encoding="UTF-8"?>

<dolfin xmlns:dolfin="http://www.fenics.org/dolfin/">
  <function> 
  <mesh celltype="interval" dim="1">
    <vertices size="51">
      <vertex index="0" x="0"/>
      <vertex index="1" x="0.02"/>
      <vertex index="2" x="0.04"/>
      <vertex index="3" x="0.06"/>
      <vertex index="4" x="0.08"/>
      <vertex index="5" x="0.1"/>
      <vertex index="6" x="0.12"/>
      <vertex index="7" x="0.14"/>
      <vertex index="8" x="0.16"/>
      <vertex index="9" x="0.18"/>
      <vertex index="10" x="0.2"/>
      <vertex index="11" x="0.22"/>
      <vertex index="12" x="0.24"/>
      <vertex index="13" x="0.26"/>
      <vertex index="14" x="0.28"/>
      <vertex index="15" x="0.3"/>
      <vertex index="16" x="0.32"/>
      <vertex index="17" x="0.34"/>
      <vertex index="18" x="0.36"/>
      <vertex index="19" x="0.38"/>
      <vertex index="20" x="0.4"/>
      <vertex index="21" x="0.42"/>
      <vertex index="22" x="0.44"/>
      <vertex index="23" x="0.46"/>
      <vertex index="24" x="0.48"/>
      <vertex index="25" x="0.5"/>
      <vertex index="26" x="0.52"/>
      <vertex index="27" x="0.54"/>
      <vertex index="28" x="0.56"/>
      <vertex index="29" x="0.58"/>
      <vertex index="30" x="0.6"/>
      <vertex index="31" x="0.62"/>
      <vertex index="32" x="0.64"/>
      <vertex index="33" x="0.66"/>
      <vertex index="34" x="0.68"/>
      <vertex index="35" x="0.7"/>
      <vertex index="36" x="0.72"/>
      <vertex index="37" x="0.74"/>
      <vertex index="38" x="0.76"/>
      <vertex index="39" x="0.78"/>
      <vertex index="40" x="0.8"/>
      <vertex index="41" x="0.82"/>
      <vertex index="42" x="0.84"/>
      <vertex index="43" x="0.86"/>
      <vertex index="44" x="0.88"/>
      <vertex index="45" x="0.9"/>
      <vertex index="46" x="0.92"/>
      <vertex index="47" x="0.94"/>
      <vertex index="48" x="0.96"/>
      <vertex index="49" x="0.98"/>
      <vertex index="50" x="1"/>
    </vertices>
    <cells size="50">
      <interval index="0" v0="0" v1="1"/>
      <interval index="1" v0="1" v1="2"/>
      <interval index="2" v0="2" v1="3"/>
      <interval index="3" v0="3" v1="4"/>
      <interval index="4" v0="4" v1="5"/>
      <interval index="5" v0="5" v1="6"/>
      <interval index="6" v0="6" v1="7"/>
      <interval index="7" v0="7" v1="8"/>
      <interval index="8" v0="8" v1="9"/>
      <interval index="9" v0="9" v1="10"/>
      <interval index="10" v0="10" v1="11"/>
      <interval index="11" v0="11" v1="12"/>
      <interval index="12" v0="12" v1="13"/>
      <interval index="13" v0="13" v1="14"/>
      <interval index="14" v0="14" v1="15"/>
      <interval index="15" v0="15" v1="16"/>
      <interval index="16" v0="16" v1="17"/>
      <interval index="17" v0="17" v1="18"/>
      <interval index="18" v0="18" v1="19"/>
      <interval index="19" v0="19" v1="20"/>
      <interval index="20" v0="20" v1="21"/>
      <interval index="21" v0="21" v1="22"/>
      <interval index="22" v0="22" v1="23"/>
      <interval index="23" v0="23" v1="24"/>
      <interval index="24" v0="24" v1="25"/>
      <interval index="25" v0="25" v1="26"/>
      <interval index="26" v0="26" v1="27"/>
      <interval index="27" v0="27" v1="28"/>
      <interval index="28" v0="28" v1="29"/>
      <interval index="29" v0="29" v1="30"/>
      <interval index="30" v0="30" v1="31"/>
      <interval index="31" v0="31" v1="32"/>
      <interval index="32" v0="32" v1="33"/>
      <interval index="33" v0="33" v1="34"/>
      <interval index="34" v0="34" v1="35"/>
      <interval index="35" v0="35" v1="36"/>
      <interval index="36" v0="36" v1="37"/>
      <interval index="37" v0="37" v1="38"/>
      <interval index="38" v0="38" v1="39"/>
      <interval index="39" v0="39" v1="40"/>
      <interval index="40" v0="40" v1="41"/>
      <interval index="41" v0="41" v1="42"/>
      <interval index="42" v0="42" v1="43"/>
      <interval index="43" v0="43" v1="44"/>
      <interval index="44" v0="44" v1="45"/>
      <interval index="45" v0="45" v1="46"/>
      <interval index="46" v0="46" v1="47"/>
      <interval index="47" v0="47" v1="48"/>
      <interval index="48" v0="48" v1="49"/>
      <interval index="49" v0="49" v1="50"/>
    </cells>
  </mesh>
  <vector size="51"> 
    <entry row="0" value="0"/>
    <entry row="1" value="0.374764603239811"/>
    <entry row="2" value="0.742910749583345"/>
    <entry row="3" value="1.09805444426946"/>
    <entry row="4" value="1.43427231086925"/>
    <entry row="5" value="1.74631142984911"/>
    <entry row="6" value="2.02977542586399"/>
    <entry row="7" value="2.2812802115458"/>
    <entry row="8" value="2.49857387048554"/>
    <entry row="9" value="2.68061643249474"/>
    <entry row="10" value="2.82761671506785"/>
    <entry row="11" value="2.94102492591846"/>
    <entry row="12" value="3.02348128864796"/>
    <entry row="13" value="3.07872251150828"/>
    <entry row="14" value="3.1114494126491"/>
    <entry row="15" value="3.12716039129139"/>
    <entry row="16" value="3.13195664419376"/>
    <entry row="17" value="3.13232602771564"/>
    <entry row="18" value="3.13491322227123"/>
    <entry row="19" value="3.1462843412132"/>
    <entry row="20" value="3.17269432299366"/>
    <entry row="21" value="3.21986534685115"/>
    <entry row="22" value="3.29278412175838"/>
    <entry row="23" value="3.39552522977078"/>
    <entry row="24" value="3.53110678192579"/>
    <entry row="25" value="3.70138350015476"/>
    <entry row="26" value="3.90698101283387"/>
    <entry row="27" value="4.14727369158694"/>
    <entry row="28" value="4.42040681448263"/>
    <entry row="29" value="4.72336227048347"/>
    <entry row="30" value="5.05206547753406"/>
    <entry row="31" value="5.40152972666169"/>
    <entry row="32" value="5.7660328386278"/>
    <entry row="33" value="6.13931987498029"/>
    <entry row="34" value="6.5148247223665"/>
    <entry row="35" value="6.8859027003722"/>
    <entry row="36" value="7.246065952638"/>
    <entry row="37" value="7.58921328240525"/>
    <entry row="38" value="7.90984629045302"/>
    <entry row="39" value="8.2032641586316"/>
    <entry row="40" value="8.46573017868908"/>
    <entry row="41" value="8.69460412702404"/>
    <entry row="42" value="8.88843579592293"/>
    <entry row="43" value="9.04701636789127"/>
    <entry row="44" value="9.17138581311754"/>
    <entry row="45" value="9.26379604801074"/>
    <entry row="46" value="9.32763115993897"/>
    <entry row="47" value="9.36728752424725"/>
    <entry row="48" value="9.38801806046922"/>
    <entry row="49" value="9.39574614503377"/>
    <entry row="50" value="9.39685577270204"/>
  </vector>
  <finiteelement signature="Lagrange finite element of degree 1 on a interval"/>
  <dofmap signature="FFC dof map for Lagrange finite element of degree 1 on a interval"/>
  </function> 
</dolfin>
import vtk
from dolfin import *

# Simple x-y plotting using VTK. For use with DOLFIN
# 
# Dag Lindbo, [email-masked], 2008-04-15

def plot_func_1D(filename,linespec):

    # extract mesh vertices and function value at vertices
    u = Function(filename)
    mesh = u.mesh()
    x = mesh.coordinates()
    y = x.copy()
    u.interpolate(y)
    
    # plot
    plot_xy(x,y,filename,linespec)

def plot_xy(x,y,title,linespec):
    
    xVal = vtk.vtkFloatArray()
    yVal = vtk.vtkFloatArray()
    
    for i in range(len(x)):
        xVal.InsertNextTuple1(x[i])
        yVal.InsertNextTuple1(y[i])
        
    # Field data is the appropriate aggregate
    # (thanks to forum posting by Sander Niemeijer)
    xy_field = vtk.vtkFieldData()
    xy_field.AllocateArrays(2)
    xy_field.AddArray(xVal)
    xy_field.AddArray(yVal)
        
    xy_fieldObject = vtk.vtkDataObject()
    xy_fieldObject.SetFieldData(xy_field)
    
    # Set up the xyplot actor
    xyplot = vtk.vtkXYPlotActor()

    text_prop = xyplot.GetTitleTextProperty()
    text_prop.SetColor(0.0, 0.0, 0.0)
    text_prop.SetFontFamilyToArial()
    
    xyplot.AddDataObjectInput(xy_fieldObject)
    xyplot.GetProperty().SetColor(0,0,0)
    xyplot.SetBorder(10);
    xyplot.GetPositionCoordinate().SetValue(0.0, 0.0, 0);
    xyplot.GetPosition2Coordinate().SetValue(1.0, 1.0, 0);
    xyplot.GetProperty().SetLineWidth(1);
    xyplot.GetProperty().SetPointSize(5);
    xyplot.SetPlotColor(0,0,0,1)
    xyplot.PlotPointsOff()
    xyplot.PlotLinesOff()
    xyplot.SetTitle(title)
    xyplot.SetXTitle("x")
    xyplot.SetYTitle("u(x)")
    xyplot.SetXValuesToValue()
    xyplot.SetAxisTitleTextProperty(text_prop)
    xyplot.SetAxisLabelTextProperty(text_prop)
    xyplot.SetTitleTextProperty(text_prop)
    xyplot.SetDataObjectXComponent(0,0)
    xyplot.SetDataObjectYComponent(0,1)
    
    # determine line style
    if linespec == ".":
        xyplot.PlotCurvePointsOn()
        xyplot.SetPlotPoints(0,1)
    elif linespec == "-":
        xyplot.PlotCurveLinesOn()
    elif linespec == ".-" or linespec == "-.":
        xyplot.PlotCurvePointsOn()
        xyplot.PlotCurveLinesOn()
        xyplot.SetPlotPoints(0,1)
    else:
        print "Invalid line spec. Exiting"
        exit(0)

    # set up the renderer
    ren = vtk.vtkRenderer()
    ren.SetBackground(1,1,1)
    
    renWin = vtk.vtkRenderWindow()
    renWin.AddRenderer(ren)
    renWin.SetSize(800,600)

    ren.AddActor2D(xyplot)
    
    # render
    iren = vtk.vtkRenderWindowInteractor()
    iren.SetRenderWindow(renWin)
    iren.Initialize()
    renWin.Render()
    iren.Start()
    
file = "poisson.xml"
plot_func_1D(file,"-.")