def Creating_z(Plane,Slope,Sediment,Start_time,End_time,Time_step):

	vtuObject = vtktools.vtu(Plane)
	vtuObject.GetFieldNames()
	gradient = vtuObject.GetScalarField('u')
	ugrid = vtk.vtkUnstructuredGrid()
	gridreader=vtk.vtkXMLUnstructuredGridReader()
	gridreader.SetFileName(Plane)
	gridreader.Update()
	ugrid = gridreader.GetOutput()
	points = ugrid.GetPoints()

	nPoints = ugrid.GetNumberOfPoints()
	for p in range(0,nPoints):
		x = (points.GetPoint(p)[:2] + (gradient[p],))
		points.SetPoint(p,x)
    
	ugrid.Update()
###################################################################################################################
	t = Start_time
	dt = Time_step
	et = End_time
	while t <= et:

		Import = Sediment + str(t) +'000000.vtu'
		NewSave = Sediment + str(t) + '_sed_slope.pvd'
		vtuObjectSed = vtktools.vtu(Import)
		vtuObjectSed.GetFieldNames()
		gradientSed = vtuObjectSed.GetScalarField('u')
		sedgrid = vtk.vtkUnstructuredGrid()
		sedgridreader=vtk.vtkXMLUnstructuredGridReader()
		sedgridreader.SetFileName(Import)
		sedgridreader.Update()
		sedgrid = sedgridreader.GetOutput()
		s = sedgrid.GetPoints()
	
		for p in range(0,nPoints):
			x = ((s.GetPoint(p)[0],) + (s.GetPoint(p)[1],) + ((gradientSed[p]+gradient[p]),))
			s.SetPoint(p,x)

		writer = vtk.vtkUnstructuredGridWriter()
		writer.SetFileName(NewSave)
		writer.SetInput(sedgrid)
		writer.Update()
		writer.Write()
		t += dt
	writer = vtk.vtkUnstructuredGridWriter()
	writer.SetFileName(Slope)
	writer.SetInput(ugrid)
	writer.Update()
	writer.Write()

def read(filenames, timestep=None):
    '''Reads an unstructured mesh with added data.

    :param filenames: The files to read from.
    :type filenames: str
    :param timestep: Time step to read from, in case of an Exodus input mesh.
    :type timestep: int, optional
    :returns mesh{2,3}d: The mesh data.
    :returns point_data: Point data read from file.
    :type point_data: dict
    :returns field_data: Field data read from file.
    :type field_data: dict
    '''
    if isinstance(filenames, (list, tuple)) and len(filenames)==1:
        filenames = filenames[0]

    if isinstance(filenames, basestring):
        filename = filenames
        # serial files
        extension = os.path.splitext(filename)[1]

        import re
        # setup the reader
        # TODO Most readers have CanReadFile() -- use that.
        if extension == '.vtu':
            from vtk import vtkXMLUnstructuredGridReader
            reader = vtkXMLUnstructuredGridReader()
            vtk_mesh = _read_vtk_mesh(reader, filename)
        elif extension == '.vtk':
            from vtk import vtkUnstructuredGridReader
            reader = vtkUnstructuredGridReader()
            vtk_mesh = _read_vtk_mesh(reader, filename)
        elif extension == '.xmf':
            from vtk import vtkXdmfReader
            reader = vtkXdmfReader()
            vtk_mesh = _read_vtk_mesh(reader, filename)
        elif extension in [ '.ex2', '.exo', '.e' ]:
            from vtk import vtkExodusIIReader
            reader = vtkExodusIIReader()
            reader.SetFileName( filename )
            vtk_mesh = _read_exodusii_mesh(reader, timestep=timestep)
        elif re.match('[^\.]*\.e\.\d+\.\d+', filename):
            # Parallel Exodus files.
            # TODO handle with vtkPExodusIIReader
            from vtk import vtkExodusIIReader
            reader = vtkExodusIIReader()
            reader.SetFileName( filenames[0] )
            vtk_mesh = _read_exodusii_mesh(reader, timestep=timestep)
        else:
            raise RuntimeError( 'Unknown file type \'%s\'.' % filename )
    else:
        # Parallel files.
        # Assume Exodus format as we don't know anything else yet.
        from vtk import vtkPExodusIIReader
        # TODO Guess the file pattern or whatever.
        reader = vtkPExodusIIReader()
        reader.SetFileNames( filenames )
        vtk_mesh = _read_exodusii_mesh(reader, filename, timestep=timestep)

    return vtk_mesh

def Slice_VTK_data_to_VTK(inputFileName,
                             outputFileName,
                                 point, normal,
                                 resolution
                                           ):
    reader = vtk.vtkXMLUnstructuredGridReader()
    reader.SetFileName(inputFileName)
    reader.Update()
    
    plane = vtk.vtkPlane()
    plane.SetOrigin(point)
    plane.SetNormal(normal)

    cutter = vtk.vtkFiltersCorePython.vtkCutter()
    cutter.SetCutFunction(plane)
    cutter.SetInputConnection(reader.GetOutputPort())
    cutter.Update()

    #Convert tht polydata structure générated by cutter into unstructured grid by triangulation
    triFilter = vtk.vtkDataSetTriangleFilter()
    triFilter.SetInputConnection(cutter.GetOutputPort())
    triFilter.Update()
    
    writer = vtk.vtkXMLUnstructuredGridWriter()
    writer.SetInputData(triFilter.GetOutput())
    writer.SetFileName(outputFileName)
    writer.Write()

    def _test(self, fname):
        reader = vtk.vtkXMLUnstructuredGridReader()
        reader.SetFileName(VTK_DATA_ROOT + fname)

        elev = vtk.vtkElevationFilter()
        elev.SetInputConnection(reader.GetOutputPort())
        elev.SetLowPoint(-0.05, 0.05, 0)
        elev.SetHighPoint(0.05, 0.05, 0)

        grad = vtk.vtkGradientFilter()
        grad.SetInputConnection(elev.GetOutputPort())
        grad.SetInputArrayToProcess(0, 0, 0, vtk.vtkDataObject.FIELD_ASSOCIATION_POINTS, "Elevation")

        grad.Update()

        vals = (10, 0, 0)

        for i in range(3):
            r = grad.GetOutput().GetPointData().GetArray("Gradients").GetRange(i)

            self.assertTrue(abs(r[0] - vals[i]) < 1E-4)
            self.assertTrue(abs(r[1] - vals[i]) < 1E-4)

        elev.SetLowPoint(0.05, -0.05, 0)
        elev.SetHighPoint(0.05, 0.05, 0)
        grad.Update()

        vals = (0, 10, 0)

        for i in range(3):
            r = grad.GetOutput().GetPointData().GetArray("Gradients").GetRange(i)

            self.assertTrue(abs(r[0] - vals[i]) < 1E-4)
            self.assertTrue(abs(r[1] - vals[i]) < 1E-4)

    def __init__(self,offset,get=1,file='data',type='vtu',mirrorPlane=None,silent=0):
        self.offset=offset
        self.filenameout = file
        self.type = type
        self.isLoaded = False
        self.mirrorPlane=mirrorPlane
        self.silent = silent

        if type == 'vtu':
            self.filename=''.join([self.filenameout,repr(offset).zfill(4),'.vtu'])
            self.datareader = v.vtkXMLUnstructuredGridReader()
        elif type == 'pvtu':
            self.filename=''.join([self.filenameout,repr(offset).zfill(4),'.pvtu'])
            self.datareader = v.vtkXMLPUnstructuredGridReader()
        elif type == 'vti':
            self.filename=''.join([self.filenameout,repr(offset).zfill(4),'.vti'])
            self.datareader = v.vtkXMLImageDataReader()
        else:
            print 'Unknown filetype'
            
        if (self.silent == 0): print '========================================'
        if (self.silent == 0): print  'loading file %s' % (self.filename)

        if get != None:
            self.getAll()

def save_lesion(output_file, input_file, field, limits):
    reader = v.vtkXMLUnstructuredGridReader()
    reader.SetFileName(input_file)
    reader.Update()

    threshold = v.vtkThreshold()
    threshold.SetInput(reader.GetOutput())

    if limits[0] is None:
        threshold.ThresholdByLower(limits[1])
    elif limits[1] is None:
        threshold.ThresholdByUpper(limits[0])
    else:
        threshold.ThresholdBetween(*limits)

    threshold.SetInputArrayToProcess(0, 0, 0, v.vtkDataObject.FIELD_ASSOCIATION_CELLS, field)
    threshold.Update()

    extract_surface = v.vtkDataSetSurfaceFilter()
    extract_surface.SetInput(threshold.GetOutput())
    extract_surface.Update()

    writer = v.vtkXMLPolyDataWriter()
    writer.SetFileName(output_file)
    writer.SetInput(extract_surface.GetOutput())
    writer.Write()

def vtu_readvf3d(name):
        
    #generate reader and read file
    reader = vtk.vtkXMLUnstructuredGridReader()
    reader.SetFileName(name)
    reader.Update()
    vtkdata = reader.GetOutput()
    
    #Obtain vectors and points from data in file
    nPoints = vtkdata.GetNumberOfPoints()
    pointData = vtkdata.GetPointData()
    vectorData = pointData.GetVectors()
    scalarData = pointData.GetScalars()    
    try:
        vect_list = vtk2np(vectorData)
    except AttributeError:
        print 'there is no vector data'
        vect_list = None
    try:
        scal_list = vtk2np(scalarData)
    except AttributeError:
        #print 'there is no scalar data'
        scal_list = None
    pnts_list = N.array([vtkdata.GetPoint(i) for i in range(nPoints)])

    #rearange data to 3d matices    
    bounds = vtkdata.GetBounds()#(Xmin,Xmax,Ymin,Ymax,Zmin,Zmax) 
    [X,Y,Z,U,V,W] = ijk(pnts_list,bounds,vect_list) 
    
    return [X,Y,Z,U,V,W]

def Create_Topo(Slope,Plane):
	ugrid = vtk.vtkUnstructuredGrid()
	gridreader=vtk.vtkUnstructuredGridReader()
	gridreader.SetFileName(Slope)
	gridreader.Update()
	ugrid = gridreader.GetOutput()
	GeomFilt1 = vtk.vtkGeometryFilter()
	GeomFilt1.SetInput(ugrid)
	GeomFilt1.Update()
	x = GeomFilt1.GetOutput()

	u = vtk.vtkUnstructuredGrid()
	bgridreader=vtk.vtkXMLUnstructuredGridReader()
	bgridreader.SetFileName(Plane)
	bgridreader.Update()
	u = bgridreader.GetOutput() 
	GeomFilt2 = vtk.vtkGeometryFilter()
	GeomFilt2.SetInput(u)
	GeomFilt2.Update()
	y = GeomFilt2.GetOutput()

	append = vtk.vtkAppendPolyData()
	append.AddInput(x)
	append.AddInput(y)
	data = append.GetOutput()

	d = vtk.vtkDelaunay3D()
	d.SetInput(data)
	d.Update
	d.BoundingTriangulationOff()  
	d.SetTolerance(0.00001)
	d.SetAlpha(0)
	d.Update()
	z = d.GetOutput()
	return z

def readUGrid(
        filename,
        verbose=1):

    myVTK.myPrint(verbose, "*** readUGrid: " + filename + " ***")

    if ('vtk' in filename):
        ugrid_reader = vtk.vtkUnstructuredGridReader()
    elif ('vtu' in filename):
        ugrid_reader = vtk.vtkXMLUnstructuredGridReader()
    else:
        assert 0, "File must be .vtk or .vtu. Aborting."

    assert (os.path.isfile(filename)), "Wrong filename. Aborting."

    ugrid_reader.SetFileName(filename)
    ugrid_reader.Update()
    ugrid = ugrid_reader.GetOutput()

    if (verbose):
        n_points = ugrid.GetNumberOfPoints()
        print 'n_points =', n_points

        n_cells = ugrid.GetNumberOfCells()
        print 'n_cells =', n_cells

    return ugrid

def check_min_max_radii():
  # Loop over vtu's and check that each core does not extend from inner radius to outer radius:
  filelist = glob.glob("Stokes_sphere_0/Stokes_sphere_0_*.vtu")
  filelist.sort()

  smallest_radius_max = constants.outer_radius
  largest_radius_min = constants.inner_radius

  for filename in filelist:
    print('Working on file:', filename)

    # Set up reader (note that vtktools does this automatically, but vtktools isn't being used here):
    if filename[-4:] == ".vtu":
      reader=vtk.vtkXMLUnstructuredGridReader()
    else:
      raise Exception("ERROR: don't recognise file extension" + filename)
    reader.SetFileName(filename)
    reader.Update()

    # Extract and store all field data:
    data = reader.GetOutput()

    # Extract radius information:
    radius_data = data.GetPointData().GetScalars("Radius")

    # Convert radii into a nice numpy array:
    radius = numpy.array([radius_data.GetValue(p) for p in xrange(radius_data.GetNumberOfTuples())])
    
    # Check for smallest maximum radius and largest minimum radius on each core:
    smallest_radius_max = min(smallest_radius_max, numpy.max(radius))
    largest_radius_min = max(largest_radius_min, numpy.min(radius))

  return smallest_radius_max, largest_radius_min

def read(filetype, filename):
    if filetype == 'vtk':
        reader = vtk.vtkUnstructuredGridReader()
        reader.SetFileName(filename)
        reader.Update()
        vtk_mesh = reader.GetOutput()
    elif filetype == 'vtu':
        reader = vtk.vtkXMLUnstructuredGridReader()
        reader.SetFileName(filename)
        reader.Update()
        vtk_mesh = reader.GetOutput()
    elif filetype == 'exodus':
        reader = vtk.vtkExodusIIReader()
        reader.SetFileName(filename)
        vtk_mesh = _read_exodusii_mesh(reader)
    else:
        raise RuntimeError('Unknown file type \'%s\'.' % filename)

    # Explicitly extract points, cells, point data, field data
    points = vtk.util.numpy_support.vtk_to_numpy(
        vtk_mesh.GetPoints().GetData())
    cells = _read_cells(vtk_mesh)
    point_data = _read_data(vtk_mesh.GetPointData())
    cell_data = _read_data(vtk_mesh.GetCellData())
    field_data = _read_data(vtk_mesh.GetFieldData())

    return points, cells, point_data, cell_data, field_data

def get_pointdata(offset,filenameout='data',type='pvtu',attribute_mode='cell'):

    if type == 'vtu':
        filename=''.join([filenameout,repr(offset).zfill(4),'.vtu'])
        datareader = v.vtkXMLUnstructuredGridReader()
    elif type == 'pvtu':
        filename=''.join([filenameout,repr(offset).zfill(4),'.pvtu'])
        datareader = v.vtkXMLPUnstructuredGridReader()
        
    print '=== Reading ',filename,' ==='
        
    datareader.SetFileName(filename)
    datareader.Update()
    data = datareader.GetOutput()

    if attribute_mode == 'cell':
        c2p = v.vtkCellDataToPointData()
        c2p.SetInput(data)
        pointdata=c2p.GetOutput()
        pointdata.Update()
    else:
        pointdata = data 


    vtk_points=data.GetPoints().GetData()
    points=ah.vtk2array(vtk_points)
 
    print '=== Done with reading data! ==='
    return {'pointdata': pointdata, 'points': points}

def create_probe(filename):
    if(verbose):
        print "Creating probe from ", filename

    pd = vtk.vtkStructuredPoints()

    reader = vtk.vtkXMLUnstructuredGridReader()
    reader.SetFileName(filename)
    ugrid = reader.GetOutput()
    ugrid.Update()

    b = ugrid.GetBounds()
    pd.SetOrigin(b[0], b[2], b[4])
    l = [b[1] - b[0], b[3] - b[2], b[5] - b[4]]
    tot_len = float(l[0] + l[1] + l[2])
    dims = intervals
    pd.SetExtent(0, dims[0]-1, 0, dims[1] -1, 0, dims[2] -1)
    pd.SetUpdateExtent(0, dims[0]-1, 0, dims[1] -1, 0, dims[2] -1)
    pd.SetWholeExtent(0, dims[0]-1, 0, dims[1] -1, 0, dims[2] -1)
    pd.SetDimensions(dims)
    dims = [max(1, x-1) for x in dims]
    l = [max(1e-3, x) for x in l]
    sp = [l[0]/dims[0], l[1]/dims[1], l[2]/dims[2]]
    pd.SetSpacing(sp)

    return pd

 def __init__(self, module_manager):
     SimpleVTKClassModuleBase.__init__(
         self, module_manager,
         vtk.vtkXMLUnstructuredGridReader(), 'Reading vtkXMLUnstructuredGrid.',
         (), ('vtkXMLUnstructuredGrid',),
         replaceDoc=True,
         inputFunctions=None, outputFunctions=None)

 def ReadVTKXMLMeshFile(self):
     if (self.InputFileName == ''):
         self.PrintError('Error: no InputFileName.')
     self.PrintLog('Reading VTK XML mesh file.')
     reader = vtk.vtkXMLUnstructuredGridReader()
     reader.SetFileName(self.InputFileName)
     reader.Update()
     self.Mesh = reader.GetOutput()

def test_clean_unstructured_grid():
    """ Test the clean_unstructured_grid routine"""
    reader = vtk.vtkXMLUnstructuredGridReader()
    del reader
    ugrid = vtk.vtkUnstructuredGrid()
    clean_grid = IO.clean_unstructured_grid(ugrid)

    assert clean_grid

def parseDump (filename):
    reader=vtk.vtkXMLUnstructuredGridReader()
    reader.SetFileName(filename)
    reader.Update()
    output=reader.GetOutput()
    celldata=output.GetCellData()
    cells=output.GetCells()
    points=output.GetPoints()
    return cells,celldata,points

    def __init__(self, filename, mesh=None, deepcopy=False):
        """
        Initialize a the reader with a pvd or a vtu filename
        """
        if not os.path.isfile(filename):
            raise IOError("File '%s' does not excist"%filename)
        filetype = filename.split(".")[-1]
        self._name = ".".join(filename.split(".")[0:-1])
        if filetype not in ["pvd", "vtu"]:
            raise TypeError("Expected a 'pvd' or a 'vtu' file")

        # Get dirname
        dirname = os.path.dirname(filename)
        
        # Check mesh argument
        if mesh is not None and not isinstance(mesh, Mesh):
            raise TypeError, "Expected a 'Mesh' for the mesh arguments"

        # Store deepcopy argument
        self._deepcopy = deepcopy
        
        # Store mesh
        self._mesh = mesh
        
        # Initialize the filename cache
        self._filenames = []
        if filetype == "vtu":
            self._filenames.append(filename)
            self._times = np.array([])
        else:
            # Parse pvd file
            tree = ElementTree(file=filename)
            times = []
            for item in tree.iter():
                if item.tag == "DataSet":
                    self._filenames.append(os.path.join(\
                        dirname,item.attrib["file"]))
                    times.append(float(item.attrib["timestep"]))
            
            times = np.array(times, dtype='d')

            # If there are no time data stored in the file use an empty array
            if np.all(np.diff(times)==1):
                times = np.array([])

            # Store time data
            self._times = times

        # Construct file reader
        self.reader = vtk.vtkXMLUnstructuredGridReader()
        
        # Read in data from file
        self._update_vtk_data()

        # Init dolfin structures (Function, FunctionSpace)
        self._init_dolfin_data()

def vtu_points(offset,filenameout='data',):

    filename=''.join([filenameout,repr(offset).zfill(4),'.vtu'])

    datareader = v.vtkXMLUnstructuredGridReader()
    datareader.SetFileName(filename)
    datareader.Update()
    data = datareader.GetOutput()
    
    vtk_points=data.GetPoints().GetData()
    points=ah.vtk2array(vtk_points)
    return points

    def open(self, k):
        """ Open a file for reading."""
        rdr = vtk.vtkXMLUnstructuredGridReader()

        print 'loading %s'%self.data[k][1]
        rdr.SetFileName(self.data[k][1])
        rdr.Update()

        self.data[k][2] = rdr.GetOutput()
        cloc = vtk.vtkCellLocator()
        cloc.SetDataSet(self.data[k][2])
        cloc.BuildLocator()
        self.data[k][3] = cloc