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 findPointsInCell(points,
                     cell,
                     verbose=1):

    ugrid_cell = vtk.vtkUnstructuredGrid()
    ugrid_cell.SetPoints(cell.GetPoints())
    cell = vtk.vtkHexahedron()
    for k_point in xrange(8): cell.GetPointIds().SetId(k_point, k_point)
    cell_array_cell = vtk.vtkCellArray()
    cell_array_cell.InsertNextCell(cell)
    ugrid_cell.SetCells(vtk.VTK_HEXAHEDRON, cell_array_cell)

    geometry_filter = vtk.vtkGeometryFilter()
    geometry_filter.SetInputData(ugrid_cell)
    geometry_filter.Update()
    cell_boundary = geometry_filter.GetOutput()

    pdata_points = vtk.vtkPolyData()
    pdata_points.SetPoints(points)

    enclosed_points_filter = vtk.vtkSelectEnclosedPoints()
    enclosed_points_filter.SetSurfaceData(cell_boundary)
    enclosed_points_filter.SetInputData(pdata_points)
    enclosed_points_filter.Update()

    points_in_cell = [k_point for k_point in xrange(points.GetNumberOfPoints()) if enclosed_points_filter.GetOutput().GetPointData().GetArray('SelectedPoints').GetTuple(k_point)[0]]
    return points_in_cell

    def _set_input (self):        
        """ This function tries its best to generate an appropriate
        input for the Normals.  If one has an input StructuredGrid or
        StructuredPoints or even a RectilinearGrid the PolyDataNormals
        will not work.  In order for it to work an appropriate
        intermediate filter is used to create the correct output."""        
        debug ("In PolyDataNormals::_set_input ()")
        out = self.prev_fil.GetOutput ()
        f = None
        if out.IsA ('vtkStructuredGrid'):
            f = vtk.vtkStructuredGridGeometryFilter ()
        elif out.IsA ('vtkRectilinearGrid'):
            f = vtk.vtkRectilinearGridGeometryFilter ()
        elif out.IsA ('vtkStructuredPoints') or out.IsA('vtkImageData'):
            if hasattr (vtk, 'vtkImageDataGeometryFilter'):
                f = vtk.vtkImageDataGeometryFilter ()
            else:
                f = vtk.vtkStructuredPointsGeometryFilter ()
        elif out.IsA('vtkUnstructuredGrid'):
            f = vtk.vtkGeometryFilter()
        elif out.IsA('vtkPolyData'):
            f = None
        else:
            msg = "This module does not support the given "\
                  "output - %s "%(out.GetClassName ())
            raise Base.Objects.ModuleException, msg

        if f:
            f.SetInput (out)
            self.fil.SetInput (f.GetOutput ())
        else:
            self.fil.SetInput(out)

    def clip(inp, origin, normal, inside_out=False, take_cell=False):
        """
        VTK operation: clip.  A vtkGeometryFilter is used to convert the
        resulted vtkUnstructuredMesh object into a vtkPolyData object.

        @param inp: input VTK object.
        @type inp: vtk.vtkobject
        @param origin: a 3-tuple for cut origin.
        @type origin: tuple
        @param normal: a 3-tuple for cut normal.
        @type normal: tuple
        @keyword inside_out: make inside out.  Default False.
        @type inside_out: bool
        @keyword take_cell: treat the input VTK object with values on cells.
            Default False.
        @type: take_cell: bool
        @return: output VTK object.
        @rtype: vtk.vtkobject
        """
        import vtk

        pne = vtk.vtkPlane()
        pne.SetOrigin(origin)
        pne.SetNormal(normal)
        clip = vtk.vtkClipDataSet()
        if take_cell:
            clip.SetInput(inp)
        else:
            clip.SetInputConnection(inp.GetOutputPort())
        clip.SetClipFunction(pne)
        if inside_out:
            clip.InsideOutOn()
        parison = vtk.vtkGeometryFilter()
        parison.SetInputConnection(clip.GetOutputPort())
        return parison

def ugrid2pdata(
        ugrid,
        only_trianlges=False,
        verbose=0):

    mypy.my_print(verbose, "*** ugrid2pdata ***")

    filter_geometry = vtk.vtkGeometryFilter()
    if (vtk.vtkVersion.GetVTKMajorVersion() >= 6):
        filter_geometry.SetInputData(ugrid)
    else:
        filter_geometry.SetInput(ugrid)
    filter_geometry.Update()
    pdata = filter_geometry.GetOutput()

    if (only_trianlges):
        filter_triangle = vtk.vtkTriangleFilter()
        if (vtk.vtkVersion.GetVTKMajorVersion() >= 6):
            filter_triangle.SetInputData(pdata)
        else:
            filter_triangle.SetInput(pdata)
        filter_triangle.Update()
        pdata = filter_triangle.GetOutput()

    return pdata

    def findLocalCsys(self, filename):
        poly = vtk.vtkGeometryFilter()
        poly.SetInputData(self.vtkMeshes[filename])
        poly.Update()

        distanceFilter = vtk.vtkImplicitPolyDataDistance()
        distanceFilter.SetInput(poly.GetOutput())

        gradients = vtk.vtkFloatArray()
        gradients.SetNumberOfComponents(3)
        gradients.SetName("LevelSet Normals")

        distances = vtk.vtkFloatArray()
        distances.SetNumberOfComponents(1)
        distances.SetName("Signed Distances")

        N = self.vtkMeshes[filename].GetCellData().GetArray(
            "Centroids").GetNumberOfTuples()
        for i in range(N):
            g = np.zeros(3, np.float32)
            p = self.vtkMeshes[filename].GetCellData().GetArray(
                "Centroids").GetTuple3(i)
            d = distanceFilter.EvaluateFunction(p)
            distanceFilter.EvaluateGradient(p, g)
            g = old_div(np.array(g), np.linalg.norm(np.array(g)))
            gradients.InsertNextTuple(g)
            distances.InsertNextValue(d)
        self.vtkMeshes[filename].GetCellData().AddArray(gradients)
        self.vtkMeshes[filename].GetCellData().AddArray(distances)

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 setVtkGeometryFilter(self):
     """
     Set the vtk_geometry_filter.  We coerce STRUCTURED_GRID and
     UNSTRUCTURED_GRID vtk_dataset_types to be POLYDATA, so we only
     need a vtk.vtkGeometryFilter().
     """
     self.vtk_geometry_filter = vtk.vtkGeometryFilter()

 def CreateSurfaceCells(self,inMesh):
     #Remove the surface cells from the mesh
     cellDimFilter = vtkvmtkcontrib.vtkvmtkCellDimensionFilter()
     cellDimFilter.SetInput(inMesh)
     cellDimFilter.ThresholdByUpper(3)
     cellDimFilter.Update()
     volumetricMesh = cellDimFilter.GetOutput()
     
     #Get new surface cells
     geomFilter = vtk.vtkGeometryFilter()
     geomFilter.SetInput(cellDimFilter.GetOutput())
     geomFilter.Update()
     newSurfaceCells = geomFilter.GetOutput()
     
     #If the celEntityIdArray exist, project the original entity ids
     cellEntityIdsArray = newSurfaceCells.GetCellData().GetArray(self.CellEntityIdsArrayName)
     if (cellEntityIdsArray != None):
         #Convert the surface cells to poly data
         surfaceCellsToSurface = vmtkscripts.vmtkMeshToSurface()
         surfaceCellsToSurface.Mesh = newSurfaceCells
         surfaceCellsToSurface.Execute()
     
         #Get the original surface cells
         meshThreshold = vtk.vtkThreshold()
         meshThreshold.SetInput(self.Mesh)
         meshThreshold.ThresholdByUpper(self.WallCellEntityId+0.5)
         meshThreshold.SetInputArrayToProcess(0,0,0,1,self.CellEntityIdsArrayName)
         meshThreshold.Update()
         
         meshToSurface = vmtkscripts.vmtkMeshToSurface()
         meshToSurface.Mesh = meshThreshold.GetOutput()
         meshToSurface.Execute()
         
         #Project the entity ids form the old surface cells to the new surface cells
         #TODO: This is hackish(need for a tolerance), find a beeter way
         projector = vtkvmtkcontrib.vtkvmtkSurfaceProjectCellArray()
         projector.SetInput(surfaceCellsToSurface.Surface)
         projector.SetReferenceSurface(meshToSurface.Surface)
         projector.SetProjectedArrayName(self.CellEntityIdsArrayName)
         projector.SetDefaultValue(self.WallCellEntityId)
         projector.SetDistanceTolerance(self.Tolerance)
         projector.Update()
         
         #Convert the surface cells back to unstructured grid
         surfaceToMesh = vmtkscripts.vmtkSurfaceToMesh()
         surfaceToMesh.Surface = projector.GetOutput()
         surfaceToMesh.Execute()
         
         newSurfaceCells = surfaceToMesh.Mesh
 
 
     #append the new surface cells to the volumetric elements
     appendFilter = vtkvmtk.vtkvmtkAppendFilter()
     appendFilter.AddInput(volumetricMesh)
     appendFilter.AddInput(newSurfaceCells)
     appendFilter.Update()
     
     return appendFilter.GetOutput()

def extractcells(polydata, idlist=[0, 1, 2]):
    """Extract cells from polydata whose cellid is in idlist."""
    cellids = vtk.vtkIdList()  # specify cellids
    cellids.Initialize()
    for i in idlist:
        cellids.InsertNextId(i)
    extract = vtk.vtkExtractCells()  # extract cells with specified cellids
    extract.SetInput(polydata)
    extract.AddCellList(cellids)
    extraction = extract.GetOutput()
    geometry = vtk.vtkGeometryFilter()  # unstructured grid to polydata
    geometry.SetInput(extraction)
    return geometry.GetOutput()

def vtk2stl(fn_in, fn_out):

    reader = vtk.vtkDataSetReader()
    reader.SetFileName(fn_in)
    reader.Update()

    gfilter = vtk.vtkGeometryFilter()
    gfilter.SetInput(reader.GetOutput())

    writer = vtk.vtkSTLWriter()
    writer.SetFileName(fn_out)
    writer.SetInput(gfilter.GetOutput())
    writer.Write()

    def Execute(self):

        if self.Mesh == None:
            self.PrintError('Error: No input mesh.')

        if self.GenerateCaps == 1:
            if not ((self.Mesh.IsHomogeneous() == 1) & (self.Mesh.GetCellType(0) == 5)):
                self.PrintError('Error: In order to generate caps, all input mesh elements must be triangles.')
            meshToSurfaceFilter = vtk.vtkGeometryFilter()
            meshToSurfaceFilter.SetInputData(self.Mesh)
            meshToSurfaceFilter.Update()
            cap = vtkvmtk.vtkvmtkSimpleCapPolyData()
            cap.SetInputConnection(meshToSurfaceFilter.GetOutputPort())
            cap.SetCellMarkerArrayName(self.FacetMarkerArrayName)
            cap.Update()
            surfacetomesh = vtkvmtk.vtkvmtkPolyDataToUnstructuredGridFilter()
            surfacetomesh.SetInputConnection(cap.GetOutputPort())
            surfacetomesh.Update()
            self.Mesh = surfacetomesh.GetOutput()

        tetgen = vtkvmtk.vtkvmtkTetGenWrapper()
        tetgen.SetInputData(self.Mesh)
        tetgen.SetPLC(self.PLC)
        tetgen.SetRefine(self.Refine)
        tetgen.SetCoarsen(self.Coarsen)
        tetgen.SetNoBoundarySplit(self.NoBoundarySplit)
        tetgen.SetQuality(self.Quality)
        tetgen.SetMinRatio(self.MinRatio)
        tetgen.SetMinDihedral(self.MinDihedral)
        tetgen.SetMaxDihedral(self.MaxDihedral)
        tetgen.SetVarVolume(self.VarVolume)
        tetgen.SetFixedVolume(self.FixedVolume)
        tetgen.SetMaxVolume(self.MaxVolume)
        tetgen.SetRemoveSliver(self.RemoveSliver)
        tetgen.SetRegionAttrib(self.RegionAttrib)
        tetgen.SetEpsilon(self.Epsilon)
        tetgen.SetNoMerge(self.NoMerge)
        tetgen.SetDetectInter(self.DetectInter)
        tetgen.SetCheckClosure(self.CheckClosure)
        tetgen.SetOrder(self.Order)
        tetgen.SetDoCheck(self.DoCheck)
        tetgen.SetVerbose(self.Verbose)
        tetgen.SetUseSizingFunction(self.UseSizingFunction)
        tetgen.SetCellEntityIdsArrayName(self.CellEntityIdsArrayName)
        tetgen.SetTetrahedronVolumeArrayName(self.TetrahedronVolumeArrayName)
        tetgen.SetSizingFunctionArrayName(self.SizingFunctionArrayName)
        tetgen.SetOutputSurfaceElements(self.OutputSurfaceElements)
        tetgen.SetOutputVolumeElements(self.OutputVolumeElements)
        tetgen.Update()

        self.Mesh = tetgen.GetOutput()

    def update(self):
        delaunay = vtkDelaunay3D()
        delaunay.SetInput(self.input_)
        delaunay.SetTolerance(self.tolerance)
        delaunay.SetAlpha(self.alpha)
        delaunay.Update()

        geom = vtkGeometryFilter()
        geom.SetInputConnection(delaunay.GetOutputPort() )

        triangle = vtkTriangleFilter()
        triangle.SetInputConnection(geom.GetOutputPort())
        triangle.Update()
        self.output_ = triangle.GetOutput()

    def removeEndCaps(self):
        self.PrintLog("Using thresholding to remove endcaps.")

        th = vtk.vtkThreshold()
        th.SetInputData(self.Surface)
        th.SetInputArrayToProcess(0, 0, 0, 1, self.CellEntityIdsArrayName)
        th.ThresholdBetween(self.EndcapsThresholdLow, self.EndcapsThresholdHigh)
        th.Update()

        gf = vtk.vtkGeometryFilter()
        gf.SetInputConnection(th.GetOutputPort())
        gf.Update()

        self.DoubleSurface = gf.GetOutput()

def ExtractVtuGeometry(inputVtu):
  """
  Extract the geometry of a vtu. In 3D, this extracts the surface mesh.
  """
  
  filter = vtk.vtkGeometryFilter()
  filter.SetInput(inputVtu.ugrid)
  filter.Update()
  surfacePoly = filter.GetOutput()
  
  # Construct output
  result = vtu()
  result.ugrid = PolyDataToUnstructuredGrid(surfacePoly)
  
  return result

    def Execute(self):

        if self.Mesh == None:
            self.PrintError('Error: No input mesh.')

        meshToSurfaceFilter = vtk.vtkGeometryFilter()
        meshToSurfaceFilter.SetInputData(self.Mesh)
        meshToSurfaceFilter.Update()

        self.Surface = meshToSurfaceFilter.GetOutput()

        if self.CleanOutput == 1:
            cleaner = vtk.vtkCleanPolyData()
            cleaner.SetInputConnection(meshToSurfaceFilter.GetOutputPort())
            cleaner.Update()
            self.Surface = cleaner.GetOutput()

    def do_contour (self, event=None):
        debug ("In BandedSurfaceMap::do_contour ()")
        Common.state.busy ()
        if self.contour_on.get ():
            if not self.mod_m.get_scalar_data_name ():
                self.contour_on.set (0)
                msg = "Warning: No scalar data present to contour!"
                Common.print_err (msg)
                Common.state.idle ()
                return
            out = self.mod_m.GetOutput ()

            if out.IsA('vtkPolyData'):
                f = None
            elif out.IsA ('vtkStructuredGrid'):
                f = vtk.vtkStructuredGridGeometryFilter ()
            elif out.IsA ('vtkRectilinearGrid'):
                f = vtk.vtkRectilinearGridGeometryFilter ()
            elif out.IsA ('vtkStructuredPoints') or \
                 out.IsA('vtkImageData'):
                if hasattr (vtk, 'vtkImageDataGeometryFilter'):
                    f = vtk.vtkImageDataGeometryFilter ()
                else:
                    f = vtk.vtkStructuredPointsGeometryFilter ()
            elif out.IsA('vtkUnstructuredGrid'):
                f = vtk.vtkGeometryFilter()
            else:
                msg = "This module does not support the given "\
                      "output - %s "%(out.GetClassName ())
                raise Base.Objects.ModuleException, msg

            if f:
                f.SetInput (out)
                self.cont_fil.SetInput (f.GetOutput())
            else:
                self.cont_fil.SetInput (out)                

            self.map.SetInput (self.cont_fil.GetOutput ())
            self.map.SetScalarModeToUseCellData()
        else:
            self.map.SetInput (self.mod_m.GetOutput ())
            self.map.SetScalarModeToDefault()
        self.change_contour ()
        Common.state.idle ()

def filterUGridIntoPData(
        ugrid,
        only_trianlges=False,
        verbose=1):

    myVTK.myPrint(verbose, "*** filterUGridIntoPData ***")

    filter_geometry = vtk.vtkGeometryFilter()
    filter_geometry.SetInputData(ugrid)
    filter_geometry.Update()
    pdata = filter_geometry.GetOutput()

    if (only_trianlges):
        filter_triangle = vtk.vtkTriangleFilter()
        filter_triangle.SetInputData(pdata)
        filter_triangle.Update()
        pdata = filter_triangle.GetOutput()

    return pdata

    def Execute(self):

        if self.Mesh == None:
            self.PrintError('Error: No input mesh.')

        meshToSurfaceFilter = vtk.vtkGeometryFilter()
        meshToSurfaceFilter.SetInput(self.Mesh)
        meshToSurfaceFilter.Update()

        self.Surface = meshToSurfaceFilter.GetOutput()

        if self.CleanOutput == 1:
            cleaner = vtk.vtkCleanPolyData()
            cleaner.SetInput(meshToSurfaceFilter.GetOutput())
            cleaner.Update()
            self.Surface = cleaner.GetOutput()

        if self.Surface.GetSource():
            self.Surface.GetSource().UnRegisterAllOutputs()

def vtk2stl(fn_in, fn_out):

    reader = vtk.vtkDataSetReader()
    reader.SetFileName(fn_in)
    reader.Update()

    gfilter = vtk.vtkGeometryFilter()
    if vtk.VTK_MAJOR_VERSION <= 5:
        gfilter.SetInput(reader.GetOutput())
    else:
        gfilter.SetInputData(reader.GetOutput())
        gfilter.Update()

    writer = vtk.vtkSTLWriter()
    writer.SetFileName(fn_out)
    if vtk.VTK_MAJOR_VERSION <= 5:
        writer.SetInput(gfilter.GetOutput())
    else:
        writer.SetInputData(gfilter.GetOutput())
    writer.Write()