def initVTK(self, key):
        print "initializing VTK objects"
        self.__setMinMaxRange(key)
        self.__normalize(key)

        # Create our Data object
        vtk.vtkStructuredGrid()
        self.vtkStructuredGrid = vtk.vtkStructuredGrid()
        self.vtkStructuredGrid.SetDimensions([128]*3)

        # Set up Points
        self.vtkPoints = vtk.vtkPoints()
        self.vtkPoints.Allocate(128*3)

        # Set up Cells
        # self.vtkCells = vtk.vtkCellArray()

        # Setup the velocity vectors
        self.vtkVectors = vtk.vtkDoubleArray()
        self.vtkVectors.SetNumberOfComponents(3)
        self.vtkVectors.SetNumberOfTuples(self.sdfIdent.size)

        # Setup the Scalars
        self.vtkScalars = vtk.vtkDoubleArray()
        self.vtkScalars.SetName(key)

        # Allocate points, cells, scalars, and vector fields
        self.AllocateData(self.sdfIdent.size, key)

        # Now attach newly allocated objects to the grid
        self.vtkStructuredGrid.SetPoints(self.vtkPoints)
        # self.vtkStructuredGrid.SetVerts(self.vtkCells)
        self.vtkStructuredGrid.GetPointData().SetVectors(self.vtkVectors)
        self.vtkStructuredGrid.GetPointData().SetScalars(self.vtkScalars)
        self.vtkStructuredGrid.GetPointData().SetActiveScalars(key)

def _write_vtk_box(box_points, filename, dimensions):
	"""
	Private method that writes a vtk file containing FFD control points.
	
	:param numpy.ndarray box_points: coordinates of the FFD control points.
	:param string filename: name of the output file.
	:param list dimensions: dimension of the lattice in (x, y, z) directions.
	"""
	# setup points and vertices
	points = vtk.vtkPoints()
	
	for index in range(0, box_points.shape[1]):
		ind = points.InsertNextPoint(box_points[0, index], box_points[1, index], box_points[2, index])
		
	grid = vtk.vtkStructuredGrid()
	
	grid.SetPoints(points)
	grid.SetDimensions(dimensions)
	grid.Modified()
	
	writer = vtk.vtkStructuredGridWriter()
	writer.SetFileName(filename)
	
	if vtk.VTK_MAJOR_VERSION <= 5:
		grid.Update()
		writer.SetInput(grid)
	else:
		writer.SetInputData(grid)
		
	writer.Write()

    def testStructured(self):
        rt = vtk.vtkRTAnalyticSource()
        rt.SetWholeExtent(-5, 5, -5, 5, -5, 5)
        rt.Update()
        i = rt.GetOutput()

        st = vtk.vtkStructuredGrid()
        st.SetDimensions(i.GetDimensions())

        nps = i.GetNumberOfPoints()
        ps = vtk.vtkPoints()
        ps.SetNumberOfPoints(nps)
        for idx in xrange(nps):
            ps.SetPoint(idx, i.GetPoint(idx))

        st.SetPoints(ps)

        s = vtk.vtkSphere()
        s.SetRadius(2)
        s.SetCenter(0,0,0)

        c = vtk.vtkTableBasedClipDataSet()
        c.SetInputData(st)
        c.SetClipFunction(s)
        c.SetInsideOut(1)

        c.Update()

        self.assertEqual(c.GetOutput().GetNumberOfCells(), 64)

 def export2vts(self, filename):
     import vtk
     grid = self.info['grid']
     # Set grid points
     points = vtk.vtkPoints()
     for z in range (grid[2]):
       for y in range (grid[1]):
         for x in range (grid[0]):
           points.InsertNextPoint( [x , y , z ] )
     struGrid = vtk.vtkStructuredGrid()
     struGrid.SetDimensions( grid )
     struGrid.SetPoints( points )
     struGrid.Update()
     # Set point data
     grid_pt_num = struGrid.GetNumberOfPoints()
     array = vtk.vtkIntArray()
     array.SetNumberOfComponents(1) # this is 3 for a vector
     array.SetNumberOfTuples(grid_pt_num)
     array.SetName('MicrostructureID')
     matPoint = self.microstructure.reshape( [grid_pt_num] )
     for i in range(grid_pt_num):
         array.SetValue(i, matPoint[i])
     struGrid.GetPointData().AddArray(array)
     struGrid.Update()
     # Write output file
     outWriter = vtk.vtkXMLStructuredGridWriter()
     outWriter.SetDataModeToBinary()
     outWriter.SetCompressorTypeToZLib()
     outWriter.SetFileName( filename + '.vts' )
     outWriter.SetInput(struGrid)
     outWriter.Update()
     outWriter.Write()

def _write_vtk_box(box_points, filename, dimensions):
	"""
	Private method that writes a vtk file containing FFD control points.
	
	:param numpy.ndarray box_points: coordinates of the FFD control points.
	:param string filename: name of the output file.
	:param list dimensions: dimension of the lattice in (x, y, z) directions.
	
	.. warning::
			If you want to visualize in paraview the inner points, 
			you have to slice the lattice because paraview does not visualize them automatically
			even in the wireframe visualization.
	"""
	# setup points and vertices
	points = vtk.vtkPoints()
	
	for index in range(0, box_points.shape[1]):
		points.InsertNextPoint(box_points[0, index], box_points[1, index], box_points[2, index])
		
	grid = vtk.vtkStructuredGrid()
	
	grid.SetPoints(points)
	grid.SetDimensions(dimensions)
	grid.Modified()
	
	writer = vtk.vtkStructuredGridWriter()
	writer.SetFileName(filename)
	
	if vtk.VTK_MAJOR_VERSION <= 5:
		grid.Update()
		writer.SetInput(grid)
	else:
		writer.SetInputData(grid)
		
	writer.Write()

def putMaskOnVTKGrid(data,grid,actorColor=None,deep=True):
  #Ok now looking 
  msk = data.mask
  imsk =  VN.numpy_to_vtk(msk.astype(numpy.int).flat,deep=deep)
  mapper = None
  if msk is not numpy.ma.nomask:
      msk =  VN.numpy_to_vtk(numpy.logical_not(msk).astype(numpy.uint8).flat,deep=deep)
      if actorColor is not None:
          grid2 = vtk.vtkStructuredGrid()
          grid2.CopyStructure(grid)
          geoFilter = vtk.vtkDataSetSurfaceFilter()
          if grid.IsA("vtkStructuredGrid"):
              grid2.GetPointData().SetScalars(imsk)
              #grid2.SetCellVisibilityArray(imsk)
              p2c = vtk.vtkPointDataToCellData()
              p2c.SetInputData(grid2)
              geoFilter.SetInputConnection(p2c.GetOutputPort())
          else:
              grid2.GetCellData().SetScalars(imsk)
              geoFilter.SetInputData(grid)
          geoFilter.Update()
          mapper = vtk.vtkPolyDataMapper()
          mapper.SetInputData(geoFilter.GetOutput())
          lut = vtk.vtkLookupTable()
          lut.SetNumberOfTableValues(1)
          r,g,b = actorColor
          lut.SetNumberOfTableValues(2)
          lut.SetTableValue(0,r/100.,g/100.,b/100.)
          lut.SetTableValue(1,r/100.,g/100.,b/100.)
          mapper.SetLookupTable(lut)
          mapper.SetScalarRange(1,1)
      if grid.IsA("vtkStructuredGrid"):
          grid.SetPointVisibilityArray(msk)
          #grid.SetCellVisibilityArray(msk)
  return mapper

def create_mesh(M,P,T):
	n_naca_pts = 50
	write_test_file(M,P,T,-5.,5.,nsections=5)
	wing=create_wing('current_wing','output')
	n_sections=len(wing)
	n_section_pts = 2*n_naca_pts-1

	# build mesh
	vtk_model = vtk.vtkStructuredGrid()
	vtk_model.SetDimensions(n_section_pts,n_sections,1)
	# build points
	vtk_points = vtk.vtkPoints()

	for j in xrange(n_sections):
	    upper_pts = numpy.array([wing[j][1],wing[j][3]]).T
	    lower_pts = numpy.array([wing[j][2],wing[j][4]]).T
	    section_pts = numpy.concatenate((lower_pts[::-1],upper_pts[1:]))
	    for i in xrange(n_section_pts):
		vtk_points.InsertNextPoint(section_pts[i,0],wing[j][0],section_pts[i,1])
	# set points
	vtk_model.SetPoints(vtk_points)

	# convert to poly data    
	pdata_filter = vtk.vtkGeometryFilter()
	if vtk.VTK_MAJOR_VERSION <= 5:
	    pdata_filter.SetInput(vtk_model)
	else:
	    pdata_filter.SetInputData(vtk_model)
	pdata_filter.Update()
	poly_data = pdata_filter.GetOutput()

	# compute normals
	norms = vtk.vtkPolyDataNormals()
	if vtk.VTK_MAJOR_VERSION <= 5:
	    norms.SetInput(poly_data)
	else:
	    norms.SetInputData(poly_data)
	norms.ComputePointNormalsOff()
	norms.ComputeCellNormalsOn()
	norms.ConsistencyOn()
	norms.Update()

	# clean poly data
	clean_poly = vtk.vtkCleanPolyData()
	clean_poly.ToleranceIsAbsoluteOn()
	clean_poly.SetAbsoluteTolerance(1.e-6)
	if vtk.VTK_MAJOR_VERSION <= 5:
	    clean_poly.SetInput(norms.GetOutput())
	else:
	    clean_poly.SetInputData(norms.GetOutput())
	clean_poly.Update()

	# write output mesh
	writer = vtk.vtkXMLPolyDataWriter()
	if vtk.VTK_MAJOR_VERSION <= 5:
	    writer.SetInput(clean_poly.GetOutput())
	else:
	    writer.SetInputData(clean_poly.GetOutput())
	writer.SetFileName('output.vtp')
	writer.Write()

def genGridOnPoints(data1,data2,gm,deep=True):
  continents = False
  xm,xM,ym,yM = None, None, None, None
  useStructuredGrid = True
  try:
    g=data1.getGrid()
    x = g.getLongitude()[:]
    y = g.getLatitude()[:]
    continents=True
    wrap=[0,360]
    if isinstance(g,cdms2.gengrid.AbstractGenericGrid): # Ok need unstrctured grid
      useStructuredGrid = False
  except:
    #hum no grid that's much easier
    x=data1.getAxis(-1)[:]
    y=data1.getAxis(-2)[:]
    wrap=None

  if x.ndim==1:
    y = y[:,numpy.newaxis]*numpy.ones(x.shape)[numpy.newaxis,:]
    x = x[numpy.newaxis,:]*numpy.ones(y.shape)
  x=x.flatten()
  y=y.flatten()
  sh =list(x.shape)
  sh.append(1)
  x=numpy.reshape(x,sh)
  y=numpy.reshape(y,sh)
  #Ok we have our points in 2D let's create unstructured points grid
  xm=x.min()
  xM=x.max()
  ym=y.min()
  yM=y.max()
  z = numpy.zeros(x.shape)
  m3 = numpy.concatenate((x,y),axis=1)
  m3 = numpy.concatenate((m3,z),axis=1)
  deep = True
  pts = vtk.vtkPoints()
  ## Convert nupmy array to vtk ones
  ppV = numpy_to_vtk_wrapper(m3,deep=deep)
  pts.SetData(ppV)
  projection = vcs.elements["projection"][gm.projection]
  xm,xM,ym,yM = getRange(gm,xm,xM,ym,yM)
  geo, geopts = project(pts,projection,[xm,xM,ym,yM])
  ## Sets the vertics into the grid
  if useStructuredGrid:
    vg = vtk.vtkStructuredGrid()
    vg.SetDimensions(y.shape[1],y.shape[0],1)
  else:
    vg = vtk.vtkUnstructuredGrid()
  vg.SetPoints(geopts)
  return vg,xm,xM,ym,yM,continents,wrap,geo

    def addGrid(self, grid):
        nx, ny, nz = grid.shape[1:]

        self.display.append(True)

        self.grids.append(vtk.vtkStructuredGrid())

        self.grids[-1].SetExtent(0, nz-1, 0, ny-1, 0, nx-1)
        p = vtk.vtkPoints()

        shp = grid.shape
        grid.shape = (3, nx*ny*nz)
        p.SetData(vtknp.numpy_to_vtk(np.ascontiguousarray(grid.T), deep=True, array_type=vtknp.get_vtk_array_type(grid.dtype)))
        grid.shape = shp
        self.grids[-1].SetPoints(p)

        #Couleur
        color = np.random.rand(3)
        #Create a vtkOutlineFilter to draw the bounding box of the data set.
        ol = vtk.vtkOutlineFilter()
        if (vtk.vtkVersion().GetVTKMajorVersion()>=6):
          ol.SetInputData(self.grids[-1])
        else:
          ol.SetInput(self.grids[-1])      
        olm = vtk.vtkPolyDataMapper()
        olm.SetInputConnection(ol.GetOutputPort())
        ola = vtk.vtkActor()
        ola.SetMapper(olm)
        ola.GetProperty().SetColor(color)

        s=vtk.vtkShrinkFilter()
        if (vtk.vtkVersion().GetVTKMajorVersion()>=6):
          s.SetInputData(self.grids[-1])
        else:
          s.SetInput(self.grids[-1])      
        s.SetShrinkFactor(0.8)
        #
        mapper = vtk.vtkDataSetMapper()
        #map.SetInputData(data)
        mapper.SetInputConnection(s.GetOutputPort())
        act = vtk.vtkLODActor()
        act.SetMapper(mapper)
        #act.GetProperty().SetRepresentationToWireframe()
        #act.GetProperty().SetRepresentationToPoints()	
        act.GetProperty().SetColor(color)
        act.GetProperty().SetEdgeColor(color)
        act.GetProperty().EdgeVisibilityOff()	
        self.actors.append(act)
        self.setBounds()
        self.ren.SetActiveCamera(self.cam)

def plotWithVTK(n, x, u):
    # Stockage des donnees au format VTK
    pts = vtk.vtkPoints()
    pts.SetData(numpy_support.numpy_to_vtk(x))

    uvtk = numpy_support.numpy_to_vtk(u)
    uvtk.SetName("u")

    data = vtk.vtkStructuredGrid()
    nn = n + [1]*(3 - len(n))
    data.SetDimensions(nn)
    data.SetPoints(pts)
    data.GetPointData().SetScalars(uvtk)
    data.Update()

    # Sortie fichier
    w = vtk.vtkXMLStructuredGridWriter()
    w.SetFileName("poisson.vts")
    w.SetInput(data)
    w.Write()

    def testStructured2D(self):
        planes = ['XY', 'XZ', 'YZ']
        expectedNCells = [42, 34, 68]
        for plane, nCells in zip(planes,expectedNCells):
            rt = vtk.vtkRTAnalyticSource()
            if plane == 'XY':
                rt.SetWholeExtent(-5, 5, -5, 5, 0, 0)
            elif plane == 'XZ':
                rt.SetWholeExtent(-5, 5, 0, 0, -5, 5)
            else:
                rt.SetWholeExtent(0, 0, -5, 5, -5, 5)
            rt.Update()
            i = rt.GetOutput()

            st = vtk.vtkStructuredGrid()
            st.SetDimensions(i.GetDimensions())

            nps = i.GetNumberOfPoints()
            ps = vtk.vtkPoints()
            ps.SetNumberOfPoints(nps)
            for idx in range(nps):
                ps.SetPoint(idx, i.GetPoint(idx))

            st.SetPoints(ps)

            cyl = vtk.vtkCylinder()
            cyl.SetRadius(2)
            cyl.SetCenter(0,0,0)
            transform = vtk.vtkTransform()
            transform.RotateWXYZ(45,20,1,10)
            cyl.SetTransform(transform)

            c = vtk.vtkTableBasedClipDataSet()
            c.SetInputData(st)
            c.SetClipFunction(cyl)
            c.SetInsideOut(1)

            c.Update()

            self.assertEqual(c.GetOutput().GetNumberOfCells(), nCells)

ExtractGrid1.SetIncludeBoundary(0)
ExtractGrid2 = vtk.vtkExtractGrid()
ExtractGrid2.SetInputData(output)
ExtractGrid2.SetVOI(29,56,0,32,0,24)
ExtractGrid2.SetSampleRate(1,1,1)
ExtractGrid2.SetIncludeBoundary(0)
LineSourceWidget0 = vtk.vtkLineSource()
LineSourceWidget0.SetPoint1(3.05638,-3.00497,28.2211)
LineSourceWidget0.SetPoint2(3.05638,3.95916,28.2211)
LineSourceWidget0.SetResolution(20)
mbds = vtk.vtkMultiBlockDataSet()
mbds.SetNumberOfBlocks(3)
i = 0
while i < 3:
    locals()[get_variable_name("ExtractGrid", i, "")].Update()
    locals()[get_variable_name("sg", i, "")] = vtk.vtkStructuredGrid()
    locals()[get_variable_name("sg", i, "")].ShallowCopy(locals()[get_variable_name("ExtractGrid", i, "")].GetOutput())
    mbds.SetBlock(i,locals()[get_variable_name("sg", i, "")])
    del locals()[get_variable_name("sg", i, "")]
    i = i + 1

Stream0 = vtk.vtkStreamTracer()
Stream0.SetInputData(mbds)
Stream0.SetSourceConnection(LineSourceWidget0.GetOutputPort())
Stream0.SetIntegrationStepUnit(2)
Stream0.SetMaximumPropagation(20)
Stream0.SetInitialIntegrationStep(0.5)
Stream0.SetIntegrationDirection(0)
Stream0.SetIntegratorType(0)
Stream0.SetMaximumNumberOfSteps(2000)
Stream0.SetTerminalSpeed(1e-12)

points = vtk.vtkPoints()
points.InsertNextPoint(-1,1,0)
points.InsertNextPoint(0,1,0)
points.InsertNextPoint(1,1,0)
points.InsertNextPoint(-1,0,0)
points.InsertNextPoint(0,0,0)
points.InsertNextPoint(1,0,0)
points.InsertNextPoint(-1,-1,0)
points.InsertNextPoint(0,-1,0)
points.InsertNextPoint(1,-1,0)
faceColors = vtk.vtkFloatArray()
faceColors.InsertNextValue(0)
faceColors.InsertNextValue(1)
faceColors.InsertNextValue(1)
faceColors.InsertNextValue(2)
sgrid = vtk.vtkStructuredGrid()
sgrid.SetDimensions(3,3,1)
sgrid.SetPoints(points)
sgrid.GetCellData().SetScalars(faceColors)
Cell2Point = vtk.vtkCellDataToPointData()
Cell2Point.SetInputData(sgrid)
Cell2Point.PassCellDataOn()
Cell2Point.Update()
mapper = vtk.vtkDataSetMapper()
mapper.SetInputData(Cell2Point.GetStructuredGridOutput())
mapper.SetScalarModeToUsePointData()
mapper.SetScalarRange(0,2)
actor = vtk.vtkActor()
actor.SetMapper(mapper)
# Add the actors to the renderer, set the background and size
ren1.AddActor(actor)

blankedPlane = vtk.vtkStructuredGridGeometryFilter()
blankedPlane.SetInputConnection(blankIt.GetOutputPort())
blankedPlane.SetExtent(0, 100, 0, 100, 0, 0)

planeMapper = vtk.vtkPolyDataMapper()
planeMapper.SetInputConnection(blankedPlane.GetOutputPort())
planeMapper.SetScalarRange(0.197813, 0.710419)

planeActor = vtk.vtkActor()
planeActor.SetMapper(planeMapper)

# The second blanking technique uses grid data values to create the blanking.
# Here we borrow the image data and threshold on that.
#
anotherGrid = vtk.vtkStructuredGrid()
anotherGrid.CopyStructure(plane.GetOutput())
anotherGrid.GetPointData().SetScalars(blankImage.GetPointData().GetScalars())

blankGrid = vtk.vtkBlankStructuredGrid()
blankGrid.SetInputData(anotherGrid)
blankGrid.SetArrayName("blankScalars")
blankGrid.SetMinBlankingValue(-0.5)
blankGrid.SetMaxBlankingValue(0.5)

blankedPlane2 = vtk.vtkStructuredGridGeometryFilter()
blankedPlane2.SetInputConnection(blankGrid.GetOutputPort())
blankedPlane2.SetExtent(0, 100, 0, 100, 0, 0)

planeMapper2 = vtk.vtkPolyDataMapper()
planeMapper2.SetInputConnection(blankedPlane2.GetOutputPort())

def GetSource(dataType):
    s = vtk.vtkRTAnalyticSource()
    # Fake serial source
    if rank == 0:
        s.Update()

    if dataType == 'ImageData':
        return s.GetOutput()

    elif dataType == 'UnstructuredGrid':
        dst = vtk.vtkDataSetTriangleFilter()
        dst.SetInputData(s.GetOutput())
        dst.Update()
        return dst.GetOutput()

    elif dataType == 'RectilinearGrid':

        input = s.GetOutput()

        rg = vtk.vtkRectilinearGrid()
        rg.SetExtent(input.GetExtent())
        dims = input.GetDimensions()
        spacing = input.GetSpacing()

        x = vtk.vtkFloatArray()
        x.SetNumberOfTuples(dims[0])
        for i in range(dims[0]):
            x.SetValue(i, spacing[0]*i)

        y = vtk.vtkFloatArray()
        y.SetNumberOfTuples(dims[1])
        for i in range(dims[1]):
            y.SetValue(i, spacing[1]*i)

        z = vtk.vtkFloatArray()
        z.SetNumberOfTuples(dims[2])
        for i in range(dims[2]):
            z.SetValue(i, spacing[2]*i)

        rg.SetXCoordinates(x)
        rg.SetYCoordinates(y)
        rg.SetZCoordinates(z)

        rg.GetPointData().ShallowCopy(input.GetPointData())

        return rg

    elif dataType == 'StructuredGrid':

        input = s.GetOutput()

        sg = vtk.vtkStructuredGrid()
        sg.SetExtent(input.GetExtent())
        pts = vtk.vtkPoints()
        sg.SetPoints(pts)
        npts = input.GetNumberOfPoints()
        for i in xrange(npts):
            pts.InsertNextPoint(input.GetPoint(i))
        sg.GetPointData().ShallowCopy(input.GetPointData())

        return sg

VTK_DATA_ROOT = vtkGetDataRoot()

# Remove cullers so single vertex will render
ren1 = vtk.vtkRenderer()
ren1.GetCullers().RemoveAllItems()
renWin = vtk.vtkRenderWindow()
renWin.AddRenderer(ren1)
iren = vtk.vtkRenderWindowInteractor()
iren.SetRenderWindow(renWin)
cell = vtk.vtkGenericCell()
ptIds = vtk.vtkIdList()
# 0D
ZeroDPts = vtk.vtkPoints()
ZeroDPts.SetNumberOfPoints(1)
ZeroDPts.SetPoint(0,0,0,0)
ZeroDGrid = vtk.vtkStructuredGrid()
ZeroDGrid.SetDimensions(1,1,1)
ZeroDGrid.SetPoints(ZeroDPts)
ZeroDGrid.GetCell(0)
ZeroDGrid.GetCell(0,cell)
ZeroDGrid.GetCellPoints(0,ptIds)
ZeroDGeom = vtk.vtkStructuredGridGeometryFilter()
ZeroDGeom.SetInputData(ZeroDGrid)
ZeroDGeom.SetExtent(0,2,0,2,0,2)
ZeroDMapper = vtk.vtkPolyDataMapper()
ZeroDMapper.SetInputConnection(ZeroDGeom.GetOutputPort())
ZeroDActor = vtk.vtkActor()
ZeroDActor.SetMapper(ZeroDMapper)
ZeroDActor.SetPosition(0,0,0)
ren1.AddActor(ZeroDActor)
# 1D - X

from NACA import create_wing
import numpy
import vtk
import dakota_utils

n_naca_pts = 50
wing=create_wing("/home/david/Documents/Projet/Optim_INSA_5GMM/examples/Dakota/test_wing/TEST",n_naca_pts) # desole pour le chemin
n_sections=len(wing)
n_section_pts = 2*n_naca_pts-1

# build mesh
vtk_model = vtk.vtkStructuredGrid()
vtk_model.SetDimensions(n_section_pts,n_sections,1)
# build points
vtk_points = vtk.vtkPoints()

for j in xrange(n_sections):
    upper_pts = numpy.array([wing[j][1],wing[j][3]]).T
    lower_pts = numpy.array([wing[j][2],wing[j][4]]).T
    section_pts = numpy.concatenate((lower_pts[::-1],upper_pts[1:]))
    for i in xrange(n_section_pts):
        vtk_points.InsertNextPoint(section_pts[i,0],wing[j][0],section_pts[i,1])
# set points
vtk_model.SetPoints(vtk_points)

# convert to poly data    
pdata_filter = vtk.vtkGeometryFilter()
if vtk.VTK_MAJOR_VERSION <= 5:
    pdata_filter.SetInput(vtk_model)
else:
    pdata_filter.SetInputData(vtk_model)

def GetSource(dataType):
    s = vtk.vtkRTAnalyticSource()

    if dataType == 'ImageData':
        return s

    elif dataType == 'UnstructuredGrid':
        dst = vtk.vtkDataSetTriangleFilter()
        dst.SetInputConnection(s.GetOutputPort())
        return dst

    elif dataType == 'RectilinearGrid':
        s.Update()

        input = s.GetOutput()

        rg = vtk.vtkRectilinearGrid()
        rg.SetExtent(input.GetExtent())
        dims = input.GetDimensions()
        spacing = input.GetSpacing()

        x = vtk.vtkFloatArray()
        x.SetNumberOfTuples(dims[0])
        for i in range(dims[0]):
            x.SetValue(i, spacing[0]*i)

        y = vtk.vtkFloatArray()
        y.SetNumberOfTuples(dims[1])
        for i in range(dims[1]):
            y.SetValue(i, spacing[1]*i)

        z = vtk.vtkFloatArray()
        z.SetNumberOfTuples(dims[2])
        for i in range(dims[2]):
            z.SetValue(i, spacing[2]*i)

        rg.SetXCoordinates(x)
        rg.SetYCoordinates(y)
        rg.SetZCoordinates(z)

        rg.GetPointData().ShallowCopy(input.GetPointData())

        pf.SetInputData(rg)
        return pf

    elif dataType == 'StructuredGrid':
        s.Update()

        input = s.GetOutput()

        sg = vtk.vtkStructuredGrid()
        sg.SetExtent(input.GetExtent())
        pts = vtk.vtkPoints()
        sg.SetPoints(pts)
        npts = input.GetNumberOfPoints()
        for i in range(npts):
            pts.InsertNextPoint(input.GetPoint(i))
        sg.GetPointData().ShallowCopy(input.GetPointData())

        pf.SetInputData(sg)
        return pf

    elif dataType == 'Table':
        s.Update()
        input = s.GetOutput()

        table = vtk.vtkTable()
        RTData = input.GetPointData().GetArray(0)
        nbTuples = RTData.GetNumberOfTuples()

        array = vtk.vtkFloatArray()
        array.SetName("RTData")
        array.SetNumberOfTuples(nbTuples)

        for i in range(0, nbTuples):
            array.SetTuple1(i, float(RTData.GetTuple1(i)))

        table.AddColumn(array)

        pf.SetInputData(table)
        return pf

    def render_geometry(self, vx, vy, vz, mv, **kwargs):
        """ create geometry """

        self.hasData = False
        gridData = kwargs.get('gridData', True)

        if gridData:
            self.gridfunc = vtk.vtkStructuredGrid()
        else:
            self.gridfunc = vtk.vtkRectilinearGrid()

        # check data structure
        if type(vx) is not numpy.ndarray or type(vy) is not numpy.ndarray:
            logging.error('X,Y vectors must be numpy arrays')
            return False

        if type(vz) is not numpy.ndarray:
            logging.error('Z vector must be numpy array')
            return False

        if isinstance(mv, (list, tuple)) and len(mv) > 0:
            logging.error('V scalars must be numpy array')
            return False

        if len(vx) == 0 or len(vy) == 0 or len(vz) == 0:
            logging.error('Zero size data')
            return False

        mva = isinstance(mv, numpy.ndarray) and mv.any()

        if vz.ndim != 2:
            logging.error('Z must be 2 dimensional numpy matrix')
            return False

        if vz[0].size != vy.size:
            logging.error('Y dimension not match')
            return False

        if vz.transpose()[0].size != vx.size:
            logging.error('X dimension not match')
            return False

        if mva:
            if mv.size != vz.size or vz[0].size != mv[0].size or vz[1].size != mv[1].size:
                logging.error('Z and V dimension not match')
                return False

        Nx = vx.size
        Ny = vy.size
        Nz = vz.size

        # put data, z, into a 2D structured grid
        self.Points = vtk.vtkPoints()
        [self.Points.InsertNextPoint(vx[i], vy[j], vz[i, j]) for j in xrange(Ny) for i in xrange(Nx)]

        if gridData:
            self.gridfunc.SetDimensions(Nx, Ny, 1)
            self.gridfunc.SetPoints(self.Points)
        else:
            xCoords = vtk.vtkFloatArray()
            [xCoords.InsertNextValue(vx[i]) for i in xrange(Nx)]

            yCoords = vtk.vtkFloatArray()
            [yCoords.InsertNextValue(vy[j]) for j in xrange(Ny)]

            s = list(vz.flatten())
            vz = numpy.array(s)
            # vz.sort()
            Nz = vz.size

            zCoords = vtk.vtkFloatArray()
            [zCoords.InsertNextValue(vz[k]) for k in xrange(Nz)]

            vCoords = vtk.vtkFloatArray()
            [vCoords.InsertNextValue(n) for n in xrange(1)]

            self.gridfunc.SetDimensions(Nx, Ny, 1)
            self.gridfunc.SetXCoordinates(xCoords)
            self.gridfunc.SetYCoordinates(yCoords)
            self.gridfunc.SetZCoordinates(vCoords)

        # get scalar field from z/v-values
        self.Colors = vtk.vtkFloatArray()
        self.Colors.SetNumberOfComponents(1)
        self.Colors.SetNumberOfTuples(Nx * Ny)

        pt = mv if mva else vz
        [self.Colors.InsertComponent(i + j * Nx, 0, pt[i, j]) for j in xrange(Ny) for i in xrange(Nx)]

        self.gridfunc.GetPointData().SetScalars(self.Colors)
        self.hasData = True

        # scaling
        Xrange = vx.max() - vx.min()
        Yrange = vy.max() - vy.min()

        # filter none
        Zrange = numpy.nanmax(vz) - numpy.nanmin(vz)

        # must have x or y ranges
#.........这里部分代码省略.........

mbw.PointData.append(na2, 'foo')
assert mbw.GetBlock(0).GetPointData().GetNumberOfArrays() == 1
assert mbw.GetBlock(1).GetPointData().GetNumberOfArrays() == 0
assert mbw.GetBlock(0).GetPointData().GetArray(0).GetName() == 'foo'

mbw.PointData.append(algs.max(na2), "maxfoo")
assert mbw.GetBlock(0).GetPointData().GetNumberOfArrays() == 2
assert mbw.GetBlock(1).GetPointData().GetNumberOfArrays() == 1
assert mbw.GetBlock(0).GetPointData().GetArray(1).GetName() == 'maxfoo'

# --------------------------------------

mb = vtk.vtkMultiBlockDataSet()
mb.SetBlock(0, vtk.vtkImageData())
mb.SetBlock(1, vtk.vtkImageData())
assert dsa.WrapDataObject(mb).Points is na

mb = vtk.vtkMultiBlockDataSet()
mb.SetBlock(0, vtk.vtkStructuredGrid())
mb.SetBlock(1, vtk.vtkImageData())
assert dsa.WrapDataObject(mb).Points is na

mb = vtk.vtkMultiBlockDataSet()
sg = vtk.vtkStructuredGrid()
sg.SetPoints(vtk.vtkPoints())
mb.SetBlock(0, sg)
mb.SetBlock(1, vtk.vtkImageData())
assert dsa.WrapDataObject(mb).Points.Arrays[0] is not na
assert dsa.WrapDataObject(mb).Points.Arrays[1] is na