def makeVoxelGrid(nCubes):
    max_x=nCubes+1
    max_y=nCubes+1
    max_z=nCubes+1
    scale=20./nCubes #*(19./20)
    meshPoints = vtk.vtkPoints()
    meshPoints.SetNumberOfPoints(max_x*max_y*max_z)
    #i*(max_y)*(max_z)+j*(max_z)+k
    for i in xrange(max_x):
        for j in xrange(max_y):
            for k in xrange(max_z):
                meshPoints.InsertPoint(i*(max_y)*(max_z)+j*(max_z)+k,scale*i,scale*j,scale*k)

    nelements = (max_x-1)*(max_y-1)*(max_z-1)
    meshGrid = vtk.vtkUnstructuredGrid()
    meshGrid.Allocate(nelements, nelements)
    meshGrid.SetPoints(meshPoints)

    for i in range(max_x-1):                  
        for j in range(max_y-1):              
            for k in range(max_z-1):          
                cell = vtk.vtkHexahedron()
                # This is the order we need such that the normals point out.
                cell.GetPointIds().SetId(0, (i+1)*(max_y)*(max_z)+j*(max_z)+k)       # upper left front
                cell.GetPointIds().SetId(1, (i+1)*(max_y)*(max_z)+(j+1)*(max_z)+k)   # upper right front
                cell.GetPointIds().SetId(2, i*(max_y)*(max_z)+(j+1)*(max_z)+k)       # lower right front
                cell.GetPointIds().SetId(3, i*(max_y)*(max_z)+j*(max_z)+k)           # lower left front node index
                cell.GetPointIds().SetId(4, (i+1)*(max_y)*(max_z)+j*(max_z)+k+1)     # upper left back
                cell.GetPointIds().SetId(5, (i+1)*(max_y)*(max_z)+(j+1)*(max_z)+k+1) # upper right back
                cell.GetPointIds().SetId(6, i*(max_y)*(max_z)+(j+1)*(max_z)+k+1)     # lower right back
                cell.GetPointIds().SetId(7, i*(max_y)*(max_z)+j*(max_z)+k+1)         # lower left back
                meshGrid.InsertNextCell(cell.GetCellType(), cell.GetPointIds())

    return meshGrid

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 makeVoxelGrid(cubesize, scale):
    max_x=cubesize
    max_y=cubesize
    max_z=cubesize
    meshPoints = vtk.vtkPoints()
    meshPoints.SetNumberOfPoints(max_x*max_y*max_z)
    #i*(max_y)*(max_z)+j*(max_z)+k
    for i in xrange(max_x):
        for j in xrange(max_y):
            for k in xrange(max_z):
                meshPoints.InsertPoint(i*(max_y)*(max_z)+j*(max_z)+k,scale*i,scale*j,scale*k)

    nelements = (max_x-1)*(max_y-1)*(max_z-1)
    meshGrid = vtk.vtkUnstructuredGrid()
    meshGrid.Allocate(nelements, nelements)
    meshGrid.SetPoints(meshPoints)

    for i in range(max_x-1):                  
        for j in range(max_y-1):              
            for k in range(max_z-1):          
                cell = vtk.vtkHexahedron()
                cell.GetPointIds().SetId(0, i*(max_y)*(max_z)+j*(max_z)+k)           # lower left front node index
                cell.GetPointIds().SetId(1, i*(max_y)*(max_z)+(j+1)*(max_z)+k)       # lower right front
                cell.GetPointIds().SetId(2, (i+1)*(max_y)*(max_z)+(j+1)*(max_z)+k)   # upper right front
                cell.GetPointIds().SetId(3, (i+1)*(max_y)*(max_z)+j*(max_z)+k)       # upper left front
                cell.GetPointIds().SetId(4, i*(max_y)*(max_z)+j*(max_z)+k+1)         # lower left back
                cell.GetPointIds().SetId(5, i*(max_y)*(max_z)+(j+1)*(max_z)+k+1)     # lower right back
                cell.GetPointIds().SetId(6, (i+1)*(max_y)*(max_z)+(j+1)*(max_z)+k+1) # upper right back
                cell.GetPointIds().SetId(7, (i+1)*(max_y)*(max_z)+j*(max_z)+k+1)     # upper left back
                meshGrid.InsertNextCell(cell.GetCellType(), cell.GetPointIds())

    return meshGrid

 def getHex(el):
     hex = vtk.vtkHexahedron()
     hex.GetPointIds().SetId(0, el[0])
     hex.GetPointIds().SetId(1, el[1])
     hex.GetPointIds().SetId(2, el[2])
     hex.GetPointIds().SetId(3, el[3])
     hex.GetPointIds().SetId(4, el[4])
     hex.GetPointIds().SetId(5, el[5])
     hex.GetPointIds().SetId(6, el[6])
     hex.GetPointIds().SetId(7, el[7])
     return hex

 def defineMesh(self,nodes,elements):
     self._points    = vtk.vtkPoints()
     self._hexs = vtk.vtkUnstructuredGrid() #vtk.vtkCellArray()
     for i in range(len(nodes)):
          self._points.InsertNextPoint(nodes[i].x(),nodes[i].y(),nodes[i].z())
     for el  in elements:
          hexa  = vtk.vtkHexahedron()
          conn = el.connectivite()
          for ii,n in enumerate(conn):
              hexa.GetPointIds().SetId(ii,n)
          self._hexs.InsertNextCell(hexa.GetCellType(), hexa.GetPointIds())
     self._hexs.SetPoints(self._points)

def boxesActor(boxes, color = (0, 0, 0, 100), wireframe=False):
    """Create a vtkActor representing a list of hexahedron.

    The hexahedrons are assumed to be aligned with the coordinate axis, and are
    given using their extremal points.


    Args:
      box: [[[xmin, ymin, zmin], [xmax, ymax, zmax]], ...]
      color: RGBA color
      wireframe: if True, the boxes are drawn in wireframe mode
    """
    grid = vtk.vtkUnstructuredGrid()
    points = vtk.vtkPoints()
    colors = vtk.vtkUnsignedCharArray()
    colors.SetNumberOfComponents(4)
    colors.SetName("Colors")
    for (index, box) in enumerate(boxes):
        colors.InsertNextTuple4(color[0], color[1], color[2], 100)
        P0 = [box[0][0], box[0][1], box[0][2]]
        P1 = [box[1][0], box[0][1], box[0][2]]
        P2 = [box[1][0], box[1][1], box[0][2]]
        P3 = [box[0][0], box[1][1], box[0][2]]
        P4 = [box[0][0], box[0][1], box[1][2]]
        P5 = [box[1][0], box[0][1], box[1][2]]
        P6 = [box[1][0], box[1][1], box[1][2]]
        P7 = [box[0][0], box[1][1], box[1][2]]
        points.InsertNextPoint(P0)
        points.InsertNextPoint(P1)
        points.InsertNextPoint(P2)
        points.InsertNextPoint(P3)
        points.InsertNextPoint(P4)
        points.InsertNextPoint(P5)
        points.InsertNextPoint(P6)
        points.InsertNextPoint(P7)
        hexa = vtk.vtkHexahedron()
        hexa.GetPointIds().SetNumberOfIds(8)
        for i in range(8):
            hexa.GetPointIds().SetId(i, 8 * index + i)
        grid.InsertNextCell(hexa.GetCellType(), hexa.GetPointIds())
    grid.SetPoints(points)
    grid.GetCellData().SetScalars(colors)
    mapper = vtk.vtkDataSetMapper()
    set_mapper_input(mapper, grid)
    actor = vtk.vtkActor()
    actor.SetMapper(mapper)
    if wireframe:
        actor.GetProperty().SetRepresentationToWireframe()
        actor.GetProperty().SetLineWidth(2.)
    return actor

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

    def buildPartialVTKGrid(self, factag):
        #get points required for factag
        pointsList = self.getPointsWithFactag(factag)

        #create a lookup table so we can map the
        #cells from the global list to a local list
        points = vtk.vtkPoints()
        localIdx = 0
        ptMap = {}
        for pt in pointsList:
            ptMap[int(pt)] = localIdx
            localIdx = localIdx + 1
            p = self.mesh.coords[(pt*3):(pt*3+3)]
            points.InsertNextPoint(p)
        
        vtkgrid = vtk.vtkUnstructuredGrid()
        vtkgrid.SetPoints(points)

        #get elements that have desired factag
        felements = self.getElementsWithFactag(factag)

        #build the vtk elements
        for element in felements:
            type = element.getType()
            nodes = element.nodes
            if type == eTypes.TRI:
                cell = vtk.vtkTriangle()
            elif type == eTypes.QUAD:
                cell = vtk.vtkQuad()
            elif type == eTypes.TET:
                cell = vtk.vtkTetra()
            elif type == eTypes.PYRAMID:
                cell = vtk.vtkPyramid()
            elif type == eTypes.PRISM:
                cell = vtk.vtkWedge()  #prism
            elif type == eTypes.HEX:
                cell = vtk.vtkHexahedron()
            else:
                raise # throw an exception
            j = 0
            for n in nodes:
                localId = ptMap[int(n)]
                cell.GetPointIds().SetId(j,localId)
                j = j+1
            vtkgrid.InsertNextCell(cell.GetCellType(), cell.GetPointIds())
        return vtkgrid

def makeHexahedra(xSideLength=1.0, ySideLength=1.0, zSideLength=1.0):
    """
    This function makes a generic hexahedra object with the given side lengths.
    The angles between every edge will be 90 degrees.
    The sides will be parallel x,y,z axes of the global coordinate system.

    .. NOTE:: The side lengths have default values of 1.0. If no input is given, then the variables will be defined
    with the given default value. In this case, the value is 1.0

    :param xSideLength: float, the length of the side that runs parallel with the x-axis
    :param ySideLength: float, the length of the side that runs parallel with the y-axis
    :param zSideLength: float, the length of the side that runs parallel with the z-axis
    :return: vtkUnstructuredGrid, vtkPoints, An object from the vtk module that is used to define shapes. Also return the original points that were used to define the hexahedra.
    """
    coords = np.array([[0.0, 0.0, 0.0], # Define the points that make up the hexahedra
                       [1.0, 0.0, 0.0],
                       [1.0, 1.0, 0.0],
                       [0.0, 1.0, 0.0],
                       [0.0, 0.0, 1.0],
                       [1.0, 0.0, 1.0],
                       [1.0, 1.0, 1.0],
                       [0.0, 1.0, 1.0]])

    # multiply the coordinates by the side length. .. NOTE:: This is not matrix multiplication.
    coords = coords*np.array([xSideLength, ySideLength, zSideLength])

    hex = vtk.vtkHexahedron() # Initialize the hexahedron class instance.
    points = vtk.vtkPoints() # Initialize a vtkPoints class instance. This is specific to the VTK module that is used for visualization.

    for i in range(8): # Iterate over a list with 8 indicies (.. NOTE:: 'range(8)' is a bulit in function of python).
        points.InsertNextPoint(coords[i]) # Insert the ith coord into the points object.

        hex.GetPointIds().SetId(i, i) # While we are iteration from 0 to 7, insert the point ids

    uGrid = vtk.vtkUnstructuredGrid() # Initialize the vtk unstructured grid object.

    uGrid.SetPoints(points) # Set the points in the uGrid as the points that we previously defined.
    uGrid.InsertNextCell(hex.GetCellType(), hex.GetPointIds()) # Set the hexahedra that we previously defined as the only cell in the grid.
    return uGrid, points

def MakeHexahedron():
    '''
      A regular hexagon (cube) with all faces square and three squares around
       each vertex is created below.
 
      Setup the coordinates of eight points
       (the two faces must be in counter clockwise
       order as viewed from the outside).
 
      As an exercise you can modify the coordinates of the points to create
       seven topologically distinct convex hexahedras.
    '''
    numberOfVertices = 8
 
    # Create the points
    points = vtk.vtkPoints()
    points.InsertNextPoint(0.0, 0.0, 0.0)
    points.InsertNextPoint(1.0, 0.0, 0.0)
    points.InsertNextPoint(1.0, 1.0, 0.0)
    points.InsertNextPoint(0.0, 1.0, 0.0)
    points.InsertNextPoint(0.0, 0.0, 1.0)
    points.InsertNextPoint(1.0, 0.0, 1.0)
    points.InsertNextPoint(1.0, 1.0, 1.0)
    points.InsertNextPoint(0.0, 1.0, 1.0)
 
    # Create a hexahedron from the points
    hex_ = vtk.vtkHexahedron()
    for i in range(0, numberOfVertices):
        hex_.GetPointIds().SetId(i, i)
 
    # Add the points and hexahedron to an unstructured grid
    uGrid = vtk.vtkUnstructuredGrid()
    uGrid.SetPoints(points)
    uGrid.InsertNextCell(hex_.GetCellType(), hex_.GetPointIds())
 
    return uGrid

def getPointsInCell(
        points,
        cell,
        verbose=0):

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

    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()
    if (vtk.vtkVersion.GetVTKMajorVersion() >= 6):
        geometry_filter.SetInputData(ugrid_cell)
    else:
        geometry_filter.SetInput(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)
    if (vtk.vtkVersion.GetVTKMajorVersion() >= 6):
        enclosed_points_filter.SetInputData(pdata_points)
    else:
        enclosed_points_filter.SetInput(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').GetTuple1(k_point)]
    return points_in_cell

    def buildFullVTKGrid(self):
        # Create the points for VTK
        points = vtk.vtkPoints()
        for i in range(0, len(self.mesh.coords)/3):
            p = self.mesh.coords[(i*3):(i*3+3)]
            points.InsertNextPoint(p)

        #add the points and cells to unstructured grid
        vtkgrid = vtk.vtkUnstructuredGrid()
        vtkgrid.SetPoints(points)

        #add the VTK elements to the mesh
        for element in self.mesh.elements:
            type = element.getType()
            nodes = element.nodes
            cell = vtk.vtkTriangle()
            if type == eTypes.TRI:
                cell = vtk.vtkTriangle()
            elif type == eTypes.QUAD:
                cell = vtk.vtkQuad()
            elif type == eTypes.TET:
                cell = vtk.vtkTetra()
            elif type == eTypes.PYRAMID:
                cell = vtk.vtkPyramid()
            elif type == eTypes.PRISM:
                cell = vtk.vtkWedge()  #prism
            elif type == eTypes.HEX:
                cell = vtk.vtkHexahedron()
            else:
                raise # throw an exception
            j = 0
            for n in nodes:
                cell.GetPointIds().SetId(j,n)
                j = j+1
            vtkgrid.InsertNextCell(cell.GetCellType(), cell.GetPointIds())
        return vtkgrid

    def _make_tecplot_geometry(self, model, quads_only=False):
        nodes = model.xyz
        nnodes = self.nNodes

        points = vtk.vtkPoints()
        points.SetNumberOfPoints(nnodes)
        #self.gridResult.Allocate(self.nNodes, 1000)
        #vectorReselt.SetNumberOfComponents(3)
        #self.nidMap = {}
        #elem.SetNumberOfPoints(nNodes)

        #assert nodes is not None
        #nnodes = nodes.shape[0]

        mmax = amax(nodes, axis=0)
        mmin = amin(nodes, axis=0)
        dim_max = (mmax - mmin).max()
        self.update_axes_length(dim_max)
        for i in range(nnodes):
            points.InsertPoint(i, nodes[i, :])

        #elements = model.elements
        quads = model.quad_elements
        hexas = model.hexa_elements
        tets = model.tet_elements
        tris = model.tri_elements

        is_quads = len(quads)
        is_tris = len(tris)
        is_hexas = len(hexas)
        is_tets = len(tets)

        is_shells = is_quads + is_tris
        is_solids = is_tets + is_hexas
        if is_shells:
            is_surface = True
            self.self._create_tecplot_shells(is_quads, quads, is_tris, tris)

        elif is_solids:
            if 0:
                tris, quads = model.skin_elements()
                is_tris = bool(len(tris))
                is_quads = bool(len(quads))
                self.self._create_tecplot_shells(is_quads, quads, is_tris, tris)
            else:
                if is_tets:
                    elements = tets
                    is_surface = False
                    self.nElements = model.nelements
                    self.grid.Allocate(self.nElements, 1000)

                    nelements = elements.shape[0]
                    for eid in range(nelements):
                        elem = vtkTetra()
                        node_ids = elements[eid, :]
                        epoints = elem.GetPointIds()
                        epoints.SetId(0, node_ids[0])
                        epoints.SetId(1, node_ids[1])
                        epoints.SetId(2, node_ids[2])
                        epoints.SetId(3, node_ids[3])
                        self.grid.InsertNextCell(elem.GetCellType(), elem.GetPointIds())  #elem.GetCellType() = 5  # vtkTriangle


                if is_hexas:
                    elements = hexas
                    is_surface = True
                    # is_surface = False
                    is_volume = not is_surface

                    if is_surface:
                        self.nElements = model.nelements
                        free_faces = array(model.get_free_faces(), dtype='int32')# + 1
                        nfaces = len(free_faces)
                        elements = free_faces
                        self.grid.Allocate(nfaces, 1000)

                        for face in free_faces:
                            elem = vtkQuad()
                            epoints = elem.GetPointIds()
                            epoints.SetId(0, face[0])
                            epoints.SetId(1, face[1])
                            epoints.SetId(2, face[2])
                            epoints.SetId(3, face[3])
                            self.grid.InsertNextCell(elem.GetCellType(), elem.GetPointIds())  #elem.GetCellType() = 5  # vtkTriangle

                    elif is_volume:
                        self.nElements = model.nelements
                        self.grid.Allocate(self.nElements, 1000)

                        nelements = elements.shape[0]
                        for eid in range(nelements):
                            elem = vtkHexahedron()
                            node_ids = elements[eid, :]
                            epoints = elem.GetPointIds()
                            epoints.SetId(0, node_ids[0])
                            epoints.SetId(1, node_ids[1])
                            epoints.SetId(2, node_ids[2])
                            epoints.SetId(3, node_ids[3])
                            epoints.SetId(4, node_ids[4])
                            epoints.SetId(5, node_ids[5])
#.........这里部分代码省略.........

    def __init__(self, parent = None):
        super(VTKFrame, self).__init__(parent)

        self.vtkWidget = QVTKRenderWindowInteractor(self)
        vl = QtGui.QVBoxLayout(self)
        vl.addWidget(self.vtkWidget)
        vl.setContentsMargins(0, 0, 0, 0)
 
        self.ren = vtk.vtkRenderer()
        self.ren.SetBackground(0.1, 0.2, 0.4)
        self.vtkWidget.GetRenderWindow().AddRenderer(self.ren)
        self.iren = self.vtkWidget.GetRenderWindow().GetInteractor()
 
        # Create source
        #Setup the coordinates of eight points
        #(the two faces must be in counter clockwise order as viewd from the outside)
        P0 = [0.0, 0.0, 0.0];
        P1 = [1.0, 0.0, 0.0];
        P2 = [1.0, 1.0, 0.0];
        P3 = [0.0, 1.0, 0.0];
        P4 = [0.0, 0.0, 1.0];
        P5 = [1.0, 0.0, 1.0];
        P6 = [1.0, 1.0, 1.0];
        P7 = [0.0, 1.0, 1.0];
         
         
        #Create the points
        points = vtk.vtkPoints();
        points.InsertNextPoint(P0);
        points.InsertNextPoint(P1);
        points.InsertNextPoint(P2);
        points.InsertNextPoint(P3);
        points.InsertNextPoint(P4);
        points.InsertNextPoint(P5);
        points.InsertNextPoint(P6);
        points.InsertNextPoint(P7);
         
        #Create a hexahedron from the points
        hexa = vtk.vtkHexahedron();
        hexa.GetPointIds().SetId(0,0);
        hexa.GetPointIds().SetId(1,1);
        hexa.GetPointIds().SetId(2,2);
        hexa.GetPointIds().SetId(3,3);
        hexa.GetPointIds().SetId(4,4);
        hexa.GetPointIds().SetId(5,5);
        hexa.GetPointIds().SetId(6,6);
        hexa.GetPointIds().SetId(7,7);
         
        #Add the hexahedron to a cell array
        hexs = vtk.vtkCellArray();
        hexs.InsertNextCell(hexa);
         
        #Add the points and hexahedron to an unstructured grid
        uGrid =vtk.vtkUnstructuredGrid();
        uGrid.SetPoints(points);
        uGrid.InsertNextCell(hexa.GetCellType(), hexa.GetPointIds());
         
        surface=vtk.vtkDataSetSurfaceFilter()
        surface.SetInput(uGrid)
        surface.Update()
         
        aBeamMapper = vtk.vtkDataSetMapper()
        aBeamMapper.SetInput(surface.GetOutput())
        #aBeamMapper.SetInput(uGrid)
        aBeamActor = vtk.vtkActor()
        aBeamActor.SetMapper(aBeamMapper)
        aBeamActor.AddPosition(0,0,0)
        aBeamActor.GetProperty().SetColor(1,1,0)
        aBeamActor.GetProperty().SetOpacity(0.60)
        aBeamActor.GetProperty().EdgeVisibilityOn()
        aBeamActor.GetProperty().SetEdgeColor(1,1,1)
        aBeamActor.GetProperty().SetLineWidth(1.5)
         
        #create a plane to cut,here it cuts in the XZ direction (xz normal=(1,0,0);XY =(0,0,1),YZ =(0,1,0)
        plane=vtk.vtkPlane()
        plane.SetOrigin(0.5,0,0)
        plane.SetNormal(1,0,0)
         
        #create cutter
        cutter=vtk.vtkCutter()
        cutter.SetCutFunction(plane)
        cutter.SetInput(aBeamActor.GetMapper().GetInput())
        cutter.Update()
        cutterMapper=vtk.vtkDataSetMapper()
        cutterMapper.SetInputConnection( cutter.GetOutputPort())
         
        #create plane actor
        planeActor=vtk.vtkActor()
        planeActor.GetProperty().SetColor(1,0.5,0.5)
        planeActor.GetProperty().SetLineWidth(2)
        planeActor.SetMapper(cutterMapper)
         
        #Add the actor to the scene
        self.ren.AddActor(aBeamActor)
        self.ren.AddActor(planeActor)
        self.ren.ResetCamera()

        self._initialized = False

def readAbaqusMeshFromINP(
        mesh_filename,
        elem_types="all",
        verbose=1):

    myVTK.myPrint(verbose, "*** readAbaqusMeshFromINP: " + mesh_filename + " ***")

    points = vtk.vtkPoints()
    cell_array = vtk.vtkCellArray()

    mesh_file = open(mesh_filename, "r")

    context = ""
    for line in mesh_file:
        if (line[-1:] == "\n"): line = line[:-1]
        #myVTK.myPrint(verbose, "line =", line)

        if line.startswith("**"): continue

        if (context == "reading nodes"):
            if line.startswith("*"):
                context = ""
            else:
                splitted_line = line.split(",")
                points.InsertNextPoint([float(coord) for coord in splitted_line[1:4]])

        if (context == "reading elems"):
            if line.startswith("*"):
                context = ""
            else:
                splitted_line = line.split(",")
                assert (len(splitted_line) == 1+cell_n_points), "Wrong number of elements in line. Aborting."
                for k_point in xrange(cell_n_points): cell.GetPointIds().SetId(k_point, int(splitted_line[1+k_point])-1)
                cell_array.InsertNextCell(cell)

        if line.startswith("*NODE"):
            context = "reading nodes"
        if line.startswith("*ELEMENT"):
            if ("TYPE=F3D4" in line) and (("quad" in elem_types) or ("all" in elem_types)):
                context = "reading elems"
                cell_vtk_type = vtk.VTK_QUAD
                cell_n_points = 4
                cell = vtk.vtkQuad()
            elif ("TYPE=C3D4" in line) and (("tet" in elem_types) or ("all" in elem_types)):
                context = "reading elems"
                cell_vtk_type = vtk.VTK_TETRA
                cell_n_points = 4
                cell = vtk.vtkTetra()
            elif ("TYPE=C3D8" in line) and (("hex" in elem_types) or ("all" in elem_types)):
                context = "reading elems"
                cell_vtk_type = vtk.VTK_HEXAHEDRON
                cell_n_points = 8
                cell = vtk.vtkHexahedron()
            else:
                print "Warning: elements not read: " + line + "."

    mesh_file.close()

    ugrid = vtk.vtkUnstructuredGrid()
    ugrid.SetPoints(points)
    ugrid.SetCells(cell_vtk_type, cell_array)

    myVTK.myPrint(verbose, "n_cells = " + str(ugrid.GetNumberOfCells()))

    return ugrid

def addVoxels(
        ugrid,
        verbose=1):

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

    n_points = ugrid.GetPoints().GetNumberOfPoints()
    seed_point = int(n_points/2)

    point_locator = myVTK.getPointLocator(ugrid)

    n_closest_points = 100
    closest_points = vtk.vtkIdList()
    point_locator.FindClosestNPoints(
        n_closest_points,
        ugrid.GetPoint(seed_point),
        closest_points)

    P0 = numpy.array(ugrid.GetPoint(closest_points.GetId(0)))
    #print "P0 = " + str(P0)

    P1 = numpy.array(ugrid.GetPoint(closest_points.GetId(1)))
    #print "P1 = " + str(P1)

    X = P1-P0
    dX = numpy.linalg.norm(X)
    X /= dX
    #print "X = " + str(X)
    #print "dX = " + str(dX)

    for k_point in xrange(2, n_closest_points):
        #print "k_point = " + str(k_point)
        P2 = numpy.array(ugrid.GetPoint(closest_points.GetId(k_point)))

        Y = P2-P0
        dY = numpy.linalg.norm(Y)
        Y /= dY
        #print "P2 = " + str(P2)
        #print "Y = " + str(Y)
        #print "dY = " + str(dY)
        #print "numpy.dot(Y, X) = " + str(numpy.dot(Y, X))
        if (abs(numpy.dot(Y, X)) < 0.001):
            Y = P2-P0
            dY = numpy.linalg.norm(Y)
            Y /= dY
            break
    #print "Y = " + str(Y)
    #print "dY = " + str(dY)

    Z = numpy.cross(X, Y)
    dZ_list = []

    for k_point in xrange(2, n_closest_points):
        #print "k_point = " + str(k_point)
        P3 = numpy.array(ugrid.GetPoint(closest_points.GetId(k_point)))

        ZZ = P3-P0
        dZ = numpy.linalg.norm(ZZ)
        ZZ /= dZ
        #print "P3 = " + str(P3)
        #print "ZZ = " + str(ZZ)
        #print "dZ = " + str(dZ)
        #print "numpy.dot(ZZ, Z) = " + str(numpy.dot(ZZ, Z))
        if (abs(numpy.dot(ZZ, Z)) > 0.999):
            dZ_list.append(dZ)
    dZ = min(dZ_list)
    #print "Z = " + str(Z)
    #print "dZ = " + str(dZ)

    #print "numpy.dot(Y, X) = " + str(numpy.dot(Y, X))
    #print "numpy.dot(Z, X) = " + str(numpy.dot(Z, X))
    #print "numpy.dot(Z, Y) = " + str(numpy.dot(Z, Y))

    points = vtk.vtkPoints()

    point_locator = vtk.vtkPointLocator()
    point_locator.InitPointInsertion(points, ugrid.GetBounds())
    radius = min(dX, dY, dZ)/2

    cell = vtk.vtkHexahedron()
    cell_array = vtk.vtkCellArray()

    for k_point in xrange(n_points):
        P = numpy.array(ugrid.GetPoint(k_point))

        point_ids = []
        for pm_Z in [-1,+1]:
            for pm_Y in [-1,+1]:
                if (pm_Y == -1):   pm_X_list = [-1,+1]
                elif (pm_Y == +1): pm_X_list = [+1,-1]
                for pm_X in pm_X_list:
                    PP = P + pm_Z * (dZ/2) * Z + pm_Y * (dY/2) * Y + pm_X * (dX/2) * X

                    if (points.GetNumberOfPoints() == 0):
                        point_id = point_locator.InsertNextPoint(PP)
                    else:
                        point_id = point_locator.FindClosestInsertedPoint(PP)
                        dist = numpy.linalg.norm(points.GetPoint(point_id)-PP)
                        if (dist > radius):
                            point_id = point_locator.InsertNextPoint(PP)
#.........这里部分代码省略.........

P7 = [0.0, 1.0, 1.0]
 
 
# Create the points
points = vtk.vtkPoints()
points.InsertNextPoint(P0)
points.InsertNextPoint(P1)
points.InsertNextPoint(P2)
points.InsertNextPoint(P3)
points.InsertNextPoint(P4)
points.InsertNextPoint(P5)
points.InsertNextPoint(P6)
points.InsertNextPoint(P7)
 
# Create a hexahedron from the points
hex = vtk.vtkHexahedron() 
hex.GetPointIds().SetId(0,0)
hex.GetPointIds().SetId(1,1)
hex.GetPointIds().SetId(2,2)
hex.GetPointIds().SetId(3,3)
hex.GetPointIds().SetId(4,4)
hex.GetPointIds().SetId(5,5)
hex.GetPointIds().SetId(6,6)
hex.GetPointIds().SetId(7,7)
 
# Add the hexahedron to a cell array
hexs = vtk.vtkCellArray()
hexs.InsertNextCell(hex)
 
# Add the points and hexahedron to an unstructured grid
uGrid = vtk.vtkUnstructuredGrid()

    def read(self):
        self.num_times = 1
        self.times = [0.0]

        filename = self.config['filename']
        # check if we are reading a ground or air conc
        if filename.find('Gconc') != -1:
            dim = 2
            cell_type = vtk.VTK_QUAD
        else:
            dim = 3
            cell_type = vtk.VTK_HEXAHEDRON

        print 'reading data...'
        data = np.genfromtxt(filename, delimiter=',')
        data = reorder.reorder(data, 0, dim)
        print 'done'

        num_pts = data.shape[0]

        var = vtk.vtkFloatArray()
        varname = os.path.basename(filename).split('_')[0]
        var.SetName(varname)
        var.SetNumberOfValues(num_pts)

        print 'creating points...'
        pts = vtk.vtkPoints()
        pts.SetNumberOfPoints(num_pts)
        if dim == 2:
            for i in xrange(0, num_pts):
                pts.SetPoint(i, data[i,0], data[i,1], 0.0)
                var.SetValue(i, data[i,2])
        else:
            for i in xrange(0, num_pts):
                pts.SetPoint(i, data[i,0], data[i,1], data[i,2])
                var.SetValue(i, data[i,3])
        print 'done'

        print 'creating cells...'
        cell_array = vtk.vtkCellArray()
        coord_x, indices = np.unique(data[:,0], return_index=True)
        indices_x = np.append(indices, num_pts)
        if dim == 2:
            for i in xrange(0, coord_x.shape[0]-1):
                for j in xrange(0, indices_x[i+1]-indices_x[i]-1):
                    quad = vtk.vtkQuad()
                    quad.GetPointIds().SetId(0, indices_x[i]  +j)
                    quad.GetPointIds().SetId(1, indices_x[i+1]+j)
                    quad.GetPointIds().SetId(2, indices_x[i+1]+j+1)
                    quad.GetPointIds().SetId(3, indices_x[i]  +j+1)
                    cell_array.InsertNextCell(quad)

        else:
            for i in xrange(0, coord_x.shape[0]-1):
                plane_b = data[indices_x[i]:indices_x[i+1]]
                coord_y_b, indices_y_b = np.unique(plane_b[:,1], return_index=True)
                indices_y_b = np.append(indices_y_b, plane_b.shape[0])

                plane_f = data[indices_x[i+1]:indices_x[i+2]]
                coord_y_f, indices_y_f = np.unique(plane_f[:,1], return_index=True)
                indices_y_f = np.append(indices_y_f, plane_f.shape[0])

                for j in xrange(0, coord_y_b.shape[0]-1):
                    for k in xrange(0, indices_y_b[j+2]-indices_y_b[j+1]-1):
                        hexa = vtk.vtkHexahedron()
                        p = []
                        p.append(indices_x[i]  +indices_y_b[j]  +k  )
                        p.append(indices_x[i]  +indices_y_b[j+1]+k  )
                        p.append(indices_x[i]  +indices_y_b[j+1]+k+1)
                        p.append(indices_x[i]  +indices_y_b[j]  +k+1)
                        p.append(indices_x[i+1]+indices_y_f[j]  +k  )
                        p.append(indices_x[i+1]+indices_y_f[j+1]+k  )
                        p.append(indices_x[i+1]+indices_y_f[j+1]+k+1)
                        p.append(indices_x[i+1]+indices_y_f[j]  +k+1)
                        hexa = vtk.vtkHexahedron()
                        for l in range(0,8):
                            hexa.GetPointIds().SetId(l, p[l])
                        cell_array.InsertNextCell(hexa)
        print 'done'

        self.grid.append(vtk.vtkUnstructuredGrid())
        self.grid[0].SetPoints(pts)
        self.grid[0].SetCells(cell_type, cell_array)
        self.grid[0].GetPointData().SetScalars(var)

    newNodes[3] = nodes[2]
      
  return newNodes

if VtkSupport():
  VTK_UNKNOWN = None
  VTK_EMPTY_CELL = vtk.vtkEmptyCell().GetCellType()
  VTK_VERTEX = vtk.vtkVertex().GetCellType()
  VTK_LINE = vtk.vtkLine().GetCellType()
  VTK_QUADRATIC_LINE = vtk.vtkQuadraticEdge().GetCellType()
  VTK_TRIANGLE = vtk.vtkTriangle().GetCellType()
  VTK_QUADRATIC_TRIANGLE = vtk.vtkQuadraticTriangle().GetCellType()
  VTK_TETRAHEDRON = vtk.vtkTetra().GetCellType()
  VTK_QUADRATIC_TETRAHEDRON = vtk.vtkQuadraticTetra().GetCellType()
  VTK_QUAD = vtk.vtkQuad().GetCellType()
  VTK_HEXAHEDRON = vtk.vtkHexahedron().GetCellType()
   
  vtkTypeIds = ( \
      VTK_UNKNOWN, \
      VTK_EMPTY_CELL, \
      VTK_VERTEX, \
      VTK_LINE, VTK_QUADRATIC_LINE, \
      VTK_TRIANGLE, VTK_QUADRATIC_TRIANGLE, VTK_QUAD, \
      VTK_TETRAHEDRON, VTK_QUADRATIC_TETRAHEDRON, VTK_HEXAHEDRON \
    )
   
  class VtkType(elements.ElementType):
    """
    Class defining a VTK element type
    """

def readDynaMesh(lsdyna_mesh_filename,
                 cell_type='hexahedron',
                 verbose=0):

    myVTK.myPrint(verbose, "*** readDynaMesh: "+lsdyna_mesh_filename+" ***")

    points = vtk.vtkPoints()

    if (cell_type == 'vertex'):
        cell_vtk_type = vtk.VTK_VERTEX
        cell = vtk.vtkVertex()
        cell_array = vtk.vtkCellArray()
    elif (cell_type == 'hexahedron'):
        cell_vtk_type = vtk.VTK_HEXAHEDRON
        cell = vtk.vtkHexahedron()
        cell_array = vtk.vtkCellArray()
    else:
        print 'Wrong cell type. Aborting.'
        exit()

    myVTK.myPrint(verbose-1, "Reading Dyna mesh file...")

    mesh_file = open(lsdyna_mesh_filename, 'r')

    context = ''
    for li