Python vtk.vtkUnstructuredGrid() Examples

def createDEM_v2():

    ds = xr.open_dataset('../data/output/Peru_20160601-20180530_comp4.nc')
    
    points = vtk.vtkPoints()
    quad = vtk.vtkQuad()
    cells = vtk.vtkCellArray()
    
    numPoints = ds.south_north.size*ds.west_east.size
   
    print('Write points \n')
    for i,j in product(ds.south_north.values,ds.west_east.values):
            points.InsertNextPoint(ds.lat.isel(south_north=i,west_east=j), ds.lon.isel(south_north=i,west_east=j), ds.HGT.sel(south_north=i,west_east=j).values/6370000.0)
    
    print('Write cells \n')
    for idx in range(points.GetNumberOfPoints()-ds.west_east.size):
        if (idx%ds.west_east.size != 0):
            quad.GetPointIds().SetId(0,idx)
            quad.GetPointIds().SetId(1,idx+1)
            quad.GetPointIds().SetId(2,idx+ds.west_east.size+1)
            quad.GetPointIds().SetId(3,idx+ds.west_east.size)
            cells.InsertNextCell(quad)

    print('Create unstructured grid \n') 
    grid = vtk.vtkUnstructuredGrid()
    grid.SetPoints(points)
    grid.SetCells(vtk.VTK_QUAD, cells)
    
    writer = vtk.vtkXMLUnstructuredGridWriter()
    writer.SetFileName('cosipy.vtu')
    writer.SetInputData(grid)
    writer.Write() 

def create_voxel(self):
        numberOfVertices = 8

        points = vtk.vtkPoints()
        points.InsertNextPoint(0, 0, 0)
        points.InsertNextPoint(1, 0, 0)
        points.InsertNextPoint(0, 1, 0)
        points.InsertNextPoint(1, 1, 0)
        points.InsertNextPoint(0, 0, 1)
        points.InsertNextPoint(1, 0, 1)
        points.InsertNextPoint(0, 1, 1)
        points.InsertNextPoint(1, 1, 1)

        voxel = vtk.vtkVoxel()
        for i in range(0, numberOfVertices):
            voxel.GetPointIds().SetId(i, i)

        ugrid = vtk.vtkUnstructuredGrid()
        ugrid.SetPoints(points)
        ugrid.InsertNextCell(voxel.GetCellType(), voxel.GetPointIds())

        gfilter = vtk.vtkGeometryFilter()
        gfilter.SetInput(ugrid)
        gfilter.Update()
        return gfilter 

def create_voxel(self):
        numberOfVertices = 8

        points = vtk.vtkPoints()
        points.InsertNextPoint(0, 0, 0)
        points.InsertNextPoint(1, 0, 0)
        points.InsertNextPoint(0, 1, 0)
        points.InsertNextPoint(1, 1, 0)
        points.InsertNextPoint(0, 0, 1)
        points.InsertNextPoint(1, 0, 1)
        points.InsertNextPoint(0, 1, 1)
        points.InsertNextPoint(1, 1, 1)

        voxel = vtk.vtkVoxel()
        for i in range(0, numberOfVertices):
            voxel.GetPointIds().SetId(i, i)

        ugrid = vtk.vtkUnstructuredGrid()
        ugrid.SetPoints(points)
        ugrid.InsertNextCell(voxel.GetCellType(), voxel.GetPointIds())

        gfilter = vtk.vtkGeometryFilter()
        gfilter.SetInput(ugrid)
        gfilter.Update()
        return gfilter 

def surface_geom_to_vtk(surfgeom):
    """Convert the triangulated surface to a :class:`pyvista.UnstructuredGrid`
    object

    Args:
        surfgeom (:class:`omf.surface.SurfaceGeometry`): the surface geomotry to
            convert
    """

    output = vtk.vtkUnstructuredGrid()
    pts = vtk.vtkPoints()
    cells = vtk.vtkCellArray()

    # Generate the points
    pts.SetNumberOfPoints(surfgeom.num_nodes)
    pts.SetData(nps.numpy_to_vtk(surfgeom.vertices))

    # Generate the triangle cells
    cellConn = surfgeom.triangles.array
    cellsMat = np.concatenate(
        (np.ones((cellConn.shape[0], 1), dtype=np.int64)*cellConn.shape[1], cellConn),
        axis=1).ravel()
    cells = vtk.vtkCellArray()
    cells.SetNumberOfCells(cellConn.shape[0])
    cells.SetCells(cellConn.shape[0],
            nps.numpy_to_vtk(cellsMat, deep=True, array_type=vtk.VTK_ID_TYPE))

    # Add to output
    output.SetPoints(pts)
    output.SetCells(vtk.VTK_TRIANGLE, cells)
    return pyvista.wrap(output) 

def test(self):
        """`PercentThreshold`: make sure no errors arise"""
        data = PVGeo.model_build.CreateTensorMesh().apply()
        thresh = PercentThreshold(percent=75).apply(data, 'Random Data')
        self.assertTrue(isinstance(thresh, vtk.vtkUnstructuredGrid))
        return True


############################################################################### 

def __init__(self, **kwargs):
        FilterBase.__init__(self,
                            nInputPorts=1, inputType='vtkPointSet',
                            nOutputPorts=1, outputType='vtkUnstructuredGrid')
        self.__dx = kwargs.get('dx', None)
        self.__dy = kwargs.get('dy', None)
        self.__dz = kwargs.get('dz', None)
        self.__estimate_grid = kwargs.get('estimate', True)
        self.__safe = kwargs.get('safe', 10.0)
        self.__unique = kwargs.get('unique', True)
        self.__tolerance = kwargs.get('tolerance', None)
        self.__angle = kwargs.get('angle', 0.0) 

def clear_mesh(self):
        """Use to clean/rebuild the mesh.
        """
        self.__mesh = vtk.vtkUnstructuredGrid()
        ubcMeshReaderBase.clear_models(self) 

def __ubc_octree(self, filename_mesh, filename_models, output):
        """Wrapper to Read UBC GIF OcTree mesh and model file pairs. UBC OcTree
        models are defined using a 2-file format. The "mesh" file describes how
        the data is descritized. The "model" file lists the physical property
        values for all cells in a mesh. A model file is meaningless without an
        associated mesh file. This only handles OcTree formats

        Args:
            filename_mesh (str): The OcTree Mesh filename as an absolute path
                for the input mesh file in UBC OcTree Mesh Format
        filename_models (list(str)): The model filenames as absolute paths for
            the input model timesteps in UBC OcTree Model Format.
            output (vtkUnstructuredGrid): The output data object

        Return:
            vtkUnstructuredGrid:
                A ``vtkUnstructuredGrid`` generated from the UBC 2D/3D Mesh grid.
                Mesh is defined by the input mesh file. Cell data is defined by
                the input model file.
        """
        if self.need_to_readMesh():
            # Construct/read the mesh
            self.ubc_octree_mesh(filename_mesh, pdo=output)
            self.need_to_readMesh(flag=False)
        output.DeepCopy(self.__mesh.to_vtk())
        if self.need_to_readModels() and self.this_has_models():
            # Read the model data
            self.__models = []
            for f in filename_models:
                # Read the model data
                self.__models.append(ubcMeshReaderBase.ubc_model_3d(f))
            self.need_to_readModels(flag=False)
        return output 

def place_model_on_octree_mesh(mesh, model, data_name='Data'):
        """Places model data onto a mesh. This is for the UBC Grid data reaers
        to associate model data with the mesh grid.

        Args:
            mesh (vtkUnstructuredGrid): The ``vtkUnstructuredGrid`` that is the
                mesh to place the model data upon. Needs to have been read in by ubcOcTree
            model (np.ndarray): A NumPy float array that holds all of the data
                to place inside of the mesh's cells.
            data_name (str): The name of the model data array once placed on the
                ``vtkUnstructuredGrid``.

        Return:
            vtkUnstructuredGrid:
                The input ``vtkUnstructuredGrid`` with model data appended.
        """
        if isinstance(model, dict):
            for key in model.keys():
                mesh = OcTreeReader.place_model_on_octree_mesh(mesh, model[key], data_name=key)
            return mesh
        # Make sure this model file fits the dimensions of the mesh
        numCells = mesh.GetNumberOfCells()
        if (numCells < len(model)):
            raise _helpers.PVGeoError('This model file has more data than the given mesh has cells to hold.')
        elif (numCells > len(model)):
            raise _helpers.PVGeoError('This model file does not have enough data to fill the given mesh\'s cells.')

        # This is absolutely crucial!
        # Do not play with unless you know what you are doing!
        # Also note that this assumes ``discretize`` handles addin this array
        ind_reorder = nps.vtk_to_numpy(
            mesh.GetCellData().GetArray('index_cell_corner'))

        model = model[ind_reorder]

        # Convert data to VTK data structure and append to output
        c = interface.convert_array(model, name=data_name, deep=True)
        # THIS IS CELL DATA! Add the model data to CELL data:
        mesh.GetCellData().AddArray(c)
        return mesh 

def ubc_octree_mesh(self, FileName, pdo=None):
        """This method reads a UBC OcTree Mesh file and builds a
        ``vtkUnstructuredGrid`` of the data in the file. This method generates
        the ``vtkUnstructuredGrid`` without any data attributes.

        Args:
            FileName (str): The mesh filename as an absolute path for the input
                mesh file in UBC OcTree format.
            pdo (vtkUnstructuredGrid): A pointer to the output data object.

        Return:
            vtkUnstructuredGrid:
                a ``vtkUnstructuredGrid`` generated from the UBCMesh grid.
                Mesh is defined by the input mesh file.
                No data attributes here, simply an empty mesh. Use the
                ``place_model_on_octree_mesh()`` method to associate with model data.
        """
        try:
            self.__mesh = discretize.TreeMesh.readUBC(FileName)
        except (IOError, OSError) as fe:
            raise _helpers.PVGeoError(str(fe))
        if pdo is None:
            pdo = self.__mesh.to_vtk()
        else:
            pdo.DeepCopy(self.__mesh.to_vtk())
        return pdo 

def __init__(self, nOutputPorts=1, outputType='vtkUnstructuredGrid', **kwargs):
        ubcMeshReaderBase.__init__(self,
                                   nOutputPorts=nOutputPorts, outputType=outputType,
                                   **kwargs)

        self.__mesh = None
        self.__models = [] 

def writeVTK(mesh, filename, models=None, directory=''):
        """Makes and saves a VTK object from this mesh and given models

        Parameters
        ----------

        filename : str
            path to the output vtk file or just its name if directory is specified

        models : dict
            dictionary of numpy.array - Name('s) and array('s). Match number of cells

        directory : str
            directory where the UBC GIF file lives


        """
        vtkObj = InterfaceVTK.to_vtk(mesh, models=models)
        writers = {
            'vtkUnstructuredGrid' : InterfaceVTK._save_unstructured_grid,
            'vtkRectilinearGrid' : InterfaceVTK._save_rectilinear_grid,
            'vtkStructuredGrid' : InterfaceVTK._save_structured_grid,
            }
        key = vtkObj.GetClassName()
        try:
            write = writers[key]
        except:
            raise RuntimeError('VTK data type `%s` is not currently supported.' % key)
        return write(filename, vtkObj, directory=directory) 

def _save_unstructured_grid(filename, vtkUnstructGrid, directory=''):
        """Saves a VTK unstructured grid file (vtu) for an already generated
        :class:`pyvista.UnstructuredGrid` object.

        Parameters
        ----------

        filename : str
            path to the output vtk file or just its name if directory is specified

        directory : str
            directory where the UBC GIF file lives

        """
        if not isinstance(vtkUnstructGrid, _vtk.vtkUnstructuredGrid):
            raise RuntimeError('`_save_unstructured_grid` can only handle `vtkUnstructuredGrid` objects. `%s` is not supported.' % vtkUnstructGrid.__class__)
        # Check the extension of the filename
        fname = os.path.join(directory, filename)
        ext = os.path.splitext(fname)[1]
        if ext is '':
            fname = fname + '.vtu'
        elif ext not in '.vtu':
            raise IOError('{:s} is an incorrect extension, has to be .vtu'.format(ext))
        # Make the writer
        vtuWriteFilter = _vtkUnstWriter()
        if float(_vtk_version.split('.')[0]) >= 6:
            vtuWriteFilter.SetInputDataObject(vtkUnstructGrid)
        else:
            vtuWriteFilter.SetInput(vtkUnstructGrid)
        vtuWriteFilter.SetFileName(fname)
        # Write the file
        vtuWriteFilter.Update() 

def selectionPointId(self, grid, ptid):
        if self.actorpt is not None:
            self._vtkren.RemoveActor(self.actorpt)
        if grid is None:
            return
        ids = vtk.vtkIdTypeArray()
        ids.SetNumberOfComponents(1)
        ids.InsertNextValue(ptid)
        selectionNode = vtk.vtkSelectionNode()
        selectionNode.SetFieldType(1)
        selectionNode.SetContentType(4)
        selectionNode.SetSelectionList(ids)
        selection = vtk.vtkSelection()
        selection.AddNode(selectionNode)
        extractSelection = vtk.vtkExtractSelection()
        if VTK_VERSION_MAJOR < 8:
            raise RuntimeError("VTK version is too old, please upgrade")
        extractSelection.SetInputData(0, grid)
        extractSelection.SetInputData(1, selection)
        extractSelection.Update()
        selected = vtk.vtkUnstructuredGrid()
        selected.ShallowCopy(extractSelection.GetOutput())
        selectedMapper = vtk.vtkDataSetMapper()
        if VTK_VERSION_MAJOR < 8:
            raise RuntimeError("VTK version is too old, please upgrade")
        selectedMapper.SetInputData(selected)
        self.actorpt = vtk.vtkActor()
        self.actorpt.SetMapper(selectedMapper)
        self.actorpt.GetProperty().SetColor(0, 1, 0)
        self.actorpt.GetProperty().SetPointSize(6)
        self.actorpt.PickableOff()
        self.actorpt.DragableOff()
        self._vtkren.AddActor(self.actorpt)
        self._iren.Render() 

def selectionCellId(self, grid, ptid):
        if self.actorpt is not None:
            self._vtkren.RemoveActor(self.actorpt)
        if grid is None:
            return
        ids = vtk.vtkIdTypeArray()
        ids.SetNumberOfComponents(1)
        ids.InsertNextValue(ptid)
        selectionNode = vtk.vtkSelectionNode()
        selectionNode.SetFieldType(0)
        selectionNode.SetContentType(4)
        selectionNode.SetSelectionList(ids)
        selection = vtk.vtkSelection()
        selection.AddNode(selectionNode)
        extractSelection = vtk.vtkExtractSelection()
        if VTK_VERSION_MAJOR < 8:
            raise RuntimeError("VTK version is too old, please upgrade")
        extractSelection.SetInputData(0, grid)
        extractSelection.SetInputData(1, selection)
        extractSelection.Update()
        selected = vtk.vtkUnstructuredGrid()
        selected.ShallowCopy(extractSelection.GetOutput())
        selectedMapper = vtk.vtkDataSetMapper()
        if VTK_VERSION_MAJOR < 8:
            raise RuntimeError("VTK version is too old, please upgrade")
        selectedMapper.SetInputData(selected)
        self.actorpt = vtk.vtkActor()
        self.actorpt.SetMapper(selectedMapper)
        self.actorpt.GetProperty().SetColor(0, 1, 0)
        self.actorpt.GetProperty().SetPointSize(6)
        self.actorpt.PickableOff()
        self.actorpt.DragableOff()
        self._vtkren.AddActor(self.actorpt)
        self._iren.Render() 

def setVisibilityGrids(self, n):
        a = self._vtkren.GetActors()
        a.InitTraversal()
        s = a.GetNextItem()
        while s:
            grid = s.GetMapper().GetInput()
            if vtk.vtkUnstructuredGrid.SafeDownCast(grid):
                s.SetVisibility(n)
            s = a.GetNextItem() 

def test_multi_block_append(ant, sphere, uniform, airplane, rectilinear):
    """This puts all of the example data objects into a a MultiBlock container"""
    multi = MultiBlock()
    # Add and test examples
    datasets = (ant, sphere, uniform, airplane, rectilinear)
    for i, dataset in enumerate(datasets):
        multi.append(dataset)
        assert multi.n_blocks == i + 1
        assert isinstance(multi[i], type(dataset))
    assert multi.bounds is not None
    # Now overwrite a block
    multi[4] = pyvista.Sphere()
    assert isinstance(multi[4], PolyData)
    multi[4] = vtk.vtkUnstructuredGrid()
    assert isinstance(multi[4], pyvista.UnstructuredGrid) 

def test_multi_block_init_vtk():
    multi = vtk.vtkMultiBlockDataSet()
    multi.SetBlock(0, vtk.vtkRectilinearGrid())
    multi.SetBlock(1, vtk.vtkStructuredGrid())
    multi = MultiBlock(multi)
    assert isinstance(multi, MultiBlock)
    assert multi.n_blocks == 2
    assert isinstance(multi.GetBlock(0), RectilinearGrid)
    assert isinstance(multi.GetBlock(1), StructuredGrid)
    multi = vtk.vtkMultiBlockDataSet()
    multi.SetBlock(0, vtk.vtkRectilinearGrid())
    multi.SetBlock(1, vtk.vtkStructuredGrid())
    multi = MultiBlock(multi, deep=True)
    assert isinstance(multi, MultiBlock)
    assert multi.n_blocks == 2
    assert isinstance(multi.GetBlock(0), RectilinearGrid)
    assert isinstance(multi.GetBlock(1), StructuredGrid)
    # Test nested structure
    multi = vtk.vtkMultiBlockDataSet()
    multi.SetBlock(0, vtk.vtkRectilinearGrid())
    multi.SetBlock(1, vtk.vtkImageData())
    nested = vtk.vtkMultiBlockDataSet()
    nested.SetBlock(0, vtk.vtkUnstructuredGrid())
    nested.SetBlock(1, vtk.vtkStructuredGrid())
    multi.SetBlock(2, nested)
    # Wrap the nested structure
    multi = MultiBlock(multi)
    assert isinstance(multi, MultiBlock)
    assert multi.n_blocks == 3
    assert isinstance(multi.GetBlock(0), RectilinearGrid)
    assert isinstance(multi.GetBlock(1), UniformGrid)
    assert isinstance(multi.GetBlock(2), MultiBlock) 

def __tree_mesh_to_vtk(mesh, models=None):
        """
        Constructs a :class:`pyvista.UnstructuredGrid` object of this tree mesh and
        the given models as ``cell_arrays`` of that ``pyvista`` dataset.

        Parameters
        ----------

        mesh : discretize.TreeMesh
            The tree mesh to convert to a :class:`pyvista.UnstructuredGrid`

        models : dict(numpy.ndarray)
            Name('s) and array('s). Match number of cells

        """
        # Make the data parts for the vtu object
        # Points
        ptsMat = np.vstack((mesh.gridN, mesh.gridhN))

        # Adjust if result was 2D (voxels are pixels in 2D):
        VTK_CELL_TYPE = _vtk.VTK_VOXEL
        if ptsMat.shape[1] == 2:
            # Add Z values of 0.0 if 2D
            ptsMat = np.c_[ptsMat, np.zeros(ptsMat.shape[0])]
            VTK_CELL_TYPE = _vtk.VTK_PIXEL
        if ptsMat.shape[1] != 3:
            raise RuntimeError('Points of the mesh are improperly defined.')
        # Rotate the points to the cartesian system
        ptsMat = np.dot(ptsMat, mesh.rotation_matrix)
        # Grab the points
        vtkPts = _vtk.vtkPoints()
        vtkPts.SetData(_nps.numpy_to_vtk(ptsMat, deep=True))
        # Cells
        cellArray = [c for c in mesh]
        cellConn = np.array([cell.nodes for cell in cellArray])
        cellsMat = np.concatenate((np.ones((cellConn.shape[0], 1), dtype=int)*cellConn.shape[1], cellConn), axis=1).ravel()
        cellsArr = _vtk.vtkCellArray()
        cellsArr.SetNumberOfCells(cellConn.shape[0])
        cellsArr.SetCells(cellConn.shape[0], _nps.numpy_to_vtk(cellsMat, deep=True, array_type=_vtk.VTK_ID_TYPE))
        # Make the object
        output = _vtk.vtkUnstructuredGrid()
        output.SetPoints(vtkPts)
        output.SetCells(VTK_CELL_TYPE, cellsArr)
        # Add the level of refinement as a cell array
        cell_levels = np.array([cell._level for cell in cellArray])
        refineLevelArr = _nps.numpy_to_vtk(cell_levels, deep=1)
        refineLevelArr.SetName('octreeLevel')
        output.GetCellData().AddArray(refineLevelArr)
        ubc_order = mesh._ubc_order
        # order_ubc will re-order from treemesh ordering to UBC ordering
        # need the opposite operation
        un_order = np.empty_like(ubc_order)
        un_order[ubc_order] = np.arange(len(ubc_order))
        order = _nps.numpy_to_vtk(un_order)
        order.SetName('index_cell_corner')
        output.GetCellData().AddArray(order)
        # Assign the model('s) to the object
        return assign_cell_data(output, models=models) 

def __init__(self, *args, **kwargs):
        """Initialize the unstructured grid."""
        super().__init__()
        deep = kwargs.pop('deep', False)

        if not len(args):
            return
        if len(args) == 1:
            if isinstance(args[0], vtk.vtkUnstructuredGrid):
                if deep:
                    self.deep_copy(args[0])
                else:
                    self.shallow_copy(args[0])

            elif isinstance(args[0], str):
                self._load_file(args[0])

            elif isinstance(args[0], vtk.vtkStructuredGrid):
                vtkappend = vtk.vtkAppendFilter()
                vtkappend.AddInputData(args[0])
                vtkappend.Update()
                self.shallow_copy(vtkappend.GetOutput())

            else:
                itype = type(args[0])
                raise TypeError('Cannot work with input type %s' % itype)

        elif len(args) == 3 and VTK9:
            arg0_is_arr = isinstance(args[0], np.ndarray)
            arg1_is_arr = isinstance(args[1], np.ndarray)
            arg2_is_arr = isinstance(args[2], np.ndarray)

            if all([arg0_is_arr, arg1_is_arr, arg2_is_arr]):
                self._from_arrays(None, args[0], args[1], args[2], deep)
            else:
                raise TypeError('All input types must be np.ndarray')

        elif len(args) == 4:
            arg0_is_arr = isinstance(args[0], np.ndarray)
            arg1_is_arr = isinstance(args[1], np.ndarray)
            arg2_is_arr = isinstance(args[2], np.ndarray)
            arg3_is_arr = isinstance(args[3], np.ndarray)

            if all([arg0_is_arr, arg1_is_arr, arg2_is_arr, arg3_is_arr]):
                self._from_arrays(args[0], args[1], args[2], args[3], deep)
            else:
                raise TypeError('All input types must be np.ndarray')

        else:
            err_msg = 'Invalid parameters.  Initialization with arrays ' +\
                      'requires the following arrays:\n'
            if VTK9:
                raise TypeError(err_msg + '`cells`, `cell_type`, `points`')
            else:
                raise TypeError(err_msg + '`offset`, `cells`, `cell_type`, `points`') 

def createDEM_v1():

    ds = xr.open_dataset('../data/output/Peru_20160601-20180530_comp4.nc')
    
    points = vtk.vtkPoints()
    
    numPoints = ds.south_north.size*ds.west_east.size
   
    print('Write points \n')
    for i,j in product(ds.south_north.values,ds.west_east.values):
            points.InsertNextPoint(ds.lat.isel(south_north=i,west_east=j), ds.lon.isel(south_north=i,west_east=j), ds.HGT.isel(south_north=i,west_east=j).values/6370000.0)
    
    print('Create unstructured grid \n') 
    polydata = vtk.vtkPolyData()
    polydata.SetPoints(points)

    delaunay = vtk.vtkDelaunay2D()
    delaunay.SetInputData(polydata)
    delaunay.Update()

#    subdivision = vtk.vtkButterflySubdivisionFilter()
#    subdivision.SetInputConnection(delaunay.GetOutputPort())
#    subdivision.Update()

    #smoother = vtk.vtkWindowedSincPolyDataFilter()
    #smoother.SetInputConnection(delaunay.GetOutputPort())
    #smoother.SetNumberOfIterations(5)
    #smoother.BoundarySmoothingOff()
    #smoother.FeatureEdgeSmoothingOff()
    #smoother.SetFeatureAngle(120.0)
    #smoother.SetPassBand(.001)
    #smoother.NonManifoldSmoothingOff()
    #smoother.NormalizeCoordinatesOff()
    #smoother.Update()

    appendFilter = vtk.vtkAppendFilter()
    appendFilter.AddInputData(delaunay.GetOutput())
    appendFilter.Update()

    unstructuredGrid = vtk.vtkUnstructuredGrid()
    unstructuredGrid.ShallowCopy(appendFilter.GetOutput())

    writer = vtk.vtkXMLUnstructuredGridWriter()
    writer.SetFileName('cosipy.vtu')
    writer.SetInputData(unstructuredGrid)
    writer.Write() 

def create_actor_hexahedron(grid, color, **kwargs):
    """ Creates a VTK actor for rendering voxels using hexahedron elements.

    :param grid: grid
    :type grid: ndarray
    :param color: actor color
    :type color: list
    :return: a VTK actor
    :rtype: vtkActor
    """
    # Keyword arguments
    array_name = kwargs.get('name', "")
    array_index = kwargs.get('index', 0)

    # Create hexahedron elements
    points = vtk.vtkPoints()
    hexarray = vtk.vtkCellArray()
    for j, pt in enumerate(grid):
        tmp = vtk.vtkHexahedron()
        fb = pt[0]
        for i, v in enumerate(fb):
            points.InsertNextPoint(v)
            tmp.GetPointIds().SetId(i, i + (j * 8))
        ft = pt[-1]
        for i, v in enumerate(ft):
            points.InsertNextPoint(v)
            tmp.GetPointIds().SetId(i + 4, i + 4 + (j * 8))
        hexarray.InsertNextCell(tmp)

    # Create an unstructured grid object and add points & hexahedron elements
    ugrid = vtk.vtkUnstructuredGrid()
    ugrid.SetPoints(points)
    ugrid.SetCells(tmp.GetCellType(), hexarray)
    # ugrid.InsertNextCell(tmp.GetCellType(), tmp.GetPointIds())

    # Map unstructured grid to the graphics primitives
    mapper = vtk.vtkDataSetMapper()
    mapper.SetInputDataObject(ugrid)
    mapper.SetArrayName(array_name)
    mapper.SetArrayId(array_index)

    # Create an actor and set its properties
    actor = vtk.vtkActor()
    actor.SetMapper(mapper)
    actor.GetProperty().SetColor(*color)

    # Return the actor
    return actor