def testLinear(self):
        pts = vtk.vtkPoints()
        pts.SetNumberOfPoints(4)
        pts.InsertPoint(0, (0, 0, 0))
        pts.InsertPoint(1, (1, 0, 0))
        pts.InsertPoint(2, (0.5, 1, 0))
        pts.InsertPoint(3, (0.5, 0.5, 1))

        te = vtk.vtkTetra()
        ptIds = te.GetPointIds()
        for i in range(4):
            ptIds.SetId(i, i)

        ghosts = vtk.vtkUnsignedCharArray()
        ghosts.SetName("vtkGhostLevels")
        ghosts.SetNumberOfTuples(4)
        ghosts.SetValue(0, 1)
        ghosts.SetValue(1, 1)
        ghosts.SetValue(2, 1)
        ghosts.SetValue(3, 0)

        grid = vtk.vtkUnstructuredGrid()
        grid.Allocate(1, 1)
        grid.InsertNextCell(te.GetCellType(), te.GetPointIds())
        grid.SetPoints(pts)
        grid.GetPointData().AddArray(ghosts)

        dss = vtk.vtkDataSetSurfaceFilter()
        dss.SetInputData(grid)
        dss.Update()
        self.assertEqual(dss.GetOutput().GetNumberOfCells(), 3)

def makeTetraGrid(nCubes, flip=0):
    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 = 5*(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):
                ulf = (i+1)*(max_y)*(max_z)+j*(max_z)+k        # upper left front
                urf = (i+1)*(max_y)*(max_z)+(j+1)*(max_z)+k    # upper right front
                lrf = i*(max_y)*(max_z)+(j+1)*(max_z)+k        # lower right front
                llf = i*(max_y)*(max_z)+j*(max_z)+k            # lower left front 
                ulb = (i+1)*(max_y)*(max_z)+j*(max_z)+k+1      # upper left back
                urb = (i+1)*(max_y)*(max_z)+(j+1)*(max_z)+k+1  # upper right back 
                lrb = i*(max_y)*(max_z)+(j+1)*(max_z)+k+1      # lower right back
                llb = i*(max_y)*(max_z)+j*(max_z)+k+1          # lower left back
                
                point_order = [  # not flip
                    [[llf,urf,lrf,lrb],
                     [llf,ulf,urf,ulb],
                     [lrb,urf,urb,ulb],
                     [llf,lrb,llb,ulb],
                     [llf,ulb,urf,lrb]],
                    # flip
                    [[llf,ulf,lrf,llb],
                     [ulf,urf,lrf,urb],
                     [ulf,ulb,urb,llb],
                     [lrf,urb,lrb,llb],
                     [ulf,urb,lrf,llb]
                     ]]
                
                for o in point_order[flip]:
                    cell = vtk.vtkTetra()
                    id=0
                    for p in o:
                        cell.GetPointIds().SetId(id,p)
                        id+=1
                    meshGrid.InsertNextCell(cell.GetCellType(),cell.GetPointIds())

                flip = not flip

            if (max_z-1)%2==0:
                flip = not flip
        if (max_y-1)%2==0:
            flip = not flip

    return meshGrid

 def create_cell(elem):
     tetra = vtk.vtkTetra()
     ids = tetra.GetPointIds()
     ids.SetId(0, elem[0])
     ids.SetId(1, elem[1])
     ids.SetId(2, elem[2])
     ids.SetId(3, elem[3])
     return tetra

def meshToUnstructeredGrid(mesh):
	
	"""Converts a FiPy mesh structure to a vtkUnstructuredGrid.
	
	Works for 2D and 3D meshes.
	
	Args:
		mesh (fipy.GmshImporter3D): Some Fipy mesh object.
		
	Returns:
		vtk.vtkUnstructuredGrid	
	"""
	
	
	# Get vertex coordinates
	coords=mesh.vertexCoords
	
	if len(coords)==2:
		x,y=coords
		dim=2
	else:
		x,y,z=coords
		dim=3
		
	# Insert them as points
	points = vtk.vtkPoints()
	for i in range(len(x)):
		if dim==2:
			points.InsertNextPoint(x[i], y[i],0)
		else:	
			points.InsertNextPoint(x[i], y[i],z[i])

	# Insert tetrahedrons
	verts=mesh._getOrderedCellVertexIDs().T
		
	cellArray = vtk.vtkCellArray()
	for j,vert in enumerate(verts):
		
		if dim==3:
			tetra = vtk.vtkTetra()
		else:
			tetra = vtk.vtkTriangle()
			
		for i,v in enumerate(vert):
			tetra.GetPointIds().SetId(i, v)
		cellArray.InsertNextCell(tetra)

	# Grid
	grid = vtk.vtkUnstructuredGrid()
	grid.SetPoints(points)
	
	if dim==3:
		grid.SetCells(vtk.VTK_TETRA, cellArray)
	else:
		grid.SetCells(vtk.VTK_TRIANGLE, cellArray)
	
	return grid

  def unstructuredgrid( self, points, npars=None ):
    """add unstructured grid"""

    points = _nansplit( points )
    #assert isinstance( points, (list,tuple,numpy.ndarray) ), 'Expected list of point arrays'

    import vtk

    vtkPoints = vtk.vtkPoints()
    vtkPoints.SetNumberOfPoints( sum(pts.shape[0] for pts in points) )

    cnt = 0
    for pts in points:

      np, ndims = pts.shape
      if not npars:
        npars = ndims

      vtkelem   = None

      if np == 2:
        vtkelem = vtk.vtkLine()
      elif np == 3:
        vtkelem = vtk.vtkTriangle()
      elif np == 4:  
        if npars == 2:
          vtkelem = vtk.vtkQuad()
        elif npars == 3:
          vtkelem = vtk.vtkTetra()
      elif np == 8:
        vtkelem = vtk.vtkVoxel() # TODO hexahedron for not rectilinear NOTE ordering changes!

      if not vtkelem:
        raise Exception( 'not sure what to do with cells with ndims=%d and npoints=%d' % (ndims,np) )

      if ndims < 3:
        pts = numpy.concatenate([pts,numpy.zeros(shape=(pts.shape[0],3-ndims))],axis=1)

      cellpoints = vtkelem.GetPointIds()

      for i,point in enumerate(pts):
        vtkPoints .SetPoint( cnt, point )
        cellpoints.SetId( i, cnt )
        cnt +=1
    
      self.vtkMesh.InsertNextCell( vtkelem.GetCellType(), cellpoints )

    self.vtkMesh.SetPoints( vtkPoints )

	def initTetras(self):
		
		"""Sets up tetrahedras describing mesh."""
		
		verts=self.embryo.simulation.mesh.mesh._getOrderedCellVertexIDs().T
		
		cellArray = vtk.vtkCellArray()
		for j,vert in enumerate(verts):

			tetra = vtk.vtkTetra()
			
			for i,v in enumerate(vert):
				tetra.GetPointIds().SetId(i, v)
			cellArray.InsertNextCell(tetra)
			
		self.grid.SetCells(vtk.VTK_TETRA, cellArray)

    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 get_linear_cell(cell):
    """ Get equivalent linear cell to vtkCell cell"""
    if cell.GetCellType() in (vtk.VTK_POLY_LINE,):
        linear_cell = vtk.vtkLine()
        linear_cell.GetPoints().SetPoint(0, cell.GetPoints().GetPoint(0))
        linear_cell.GetPoints().SetPoint(1, cell.GetPoints().GetPoint(1))
    elif cell.GetCellType() in (vtk.VTK_QUADRATIC_TRIANGLE,):
        linear_cell = vtk.vtkTriangle()
        linear_cell.GetPoints().SetPoint(0, cell.GetPoints().GetPoint(0))
        linear_cell.GetPoints().SetPoint(1, cell.GetPoints().GetPoint(1))
        linear_cell.GetPoints().SetPoint(2, cell.GetPoints().GetPoint(2))
    elif cell.GetCellType() in (vtk.VTK_QUADRATIC_TETRA,):
        linear_cell = vtk.vtkTetra()
        linear_cell.GetPoints().SetPoint(0, cell.GetPoints().GetPoint(0))
        linear_cell.GetPoints().SetPoint(1, cell.GetPoints().GetPoint(1))
        linear_cell.GetPoints().SetPoint(2, cell.GetPoints().GetPoint(2))
        linear_cell.GetPoints().SetPoint(3, cell.GetPoints().GetPoint(3))
    else:
        linear_cell = cell

    return linear_cell

    def doLinear(self, ghosts, ncells):
        pts = vtk.vtkPoints()
        pts.SetNumberOfPoints(4)
        pts.InsertPoint(0, (0, 0, 0))
        pts.InsertPoint(1, (1, 0, 0))
        pts.InsertPoint(2, (0.5, 1, 0))
        pts.InsertPoint(3, (0.5, 0.5, 1))

        te = vtk.vtkTetra()
        ptIds = te.GetPointIds()
        for i in range(4):
            ptIds.SetId(i, i)

        grid = vtk.vtkUnstructuredGrid()
        grid.Allocate(1, 1)
        grid.InsertNextCell(te.GetCellType(), te.GetPointIds())
        grid.SetPoints(pts)
        grid.GetPointData().AddArray(ghosts)

        dss = vtk.vtkDataSetSurfaceFilter()
        dss.SetInputData(grid)
        dss.Update()
        self.assertEqual(dss.GetOutput().GetNumberOfCells(), ncells)

    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 MakeTetrahedron():
    '''
      Make a tetrahedron.
    '''
    numberOfVertices = 4
 
    points = vtk.vtkPoints()
    points.InsertNextPoint(0, 0, 0)
    points.InsertNextPoint(1, 0, 0)
    points.InsertNextPoint(1, 1, 0)
    points.InsertNextPoint(0, 1, 1)
 
    tetra = vtk.vtkTetra()
    for i in range(0, numberOfVertices):
        tetra.GetPointIds().SetId(i, i)
 
    cellArray = vtk.vtkCellArray()
    cellArray.InsertNextCell(tetra)
 
    unstructuredGrid = vtk.vtkUnstructuredGrid()
    unstructuredGrid.SetPoints(points)
    unstructuredGrid.SetCells(vtk.VTK_TETRA, cellArray)
 
    return unstructuredGrid

    def load_fast_geometry(self, fgrid_filename, dirname, name='main', plot=True):
        skip_reading = self._remove_old_geometry(fgrid_filename)
        if skip_reading:
            return

        model = FGridReader(log=self.log, debug=False)

        base_filename, ext = os.path.splitext(fgrid_filename)
        if '.fgrid' == ext:
            dimension_flag = 3
        #elif '.ele' == ext:
            #dimension_flag = 3
        else:
            raise RuntimeError('unsupported extension.  Use "cogsg" or "front".')

        read_loads = True
        model.read_fgrid(fgrid_filename, dimension_flag)

        dimension_flag = 3
        nodes = model.nodes
        tris = model.tris - 1
        tets = model.tets - 1

        nnodes = nodes.shape[0]
        ntris = tris.shape[0]
        ntets = tets.shape[0]

        #print('node0 = %s' % str(nodes[0, :]))
        #print('node%i = %s' % (1, str(nodes[1, :])))
        #print('node%i = %s' % (2, str(nodes[2, :])))
        #print('node%i = %s' % (nnodes, str(nodes[-1, :])))
        #print('tris.max/min = ', tris.max(), tris.min())
        #print('tets.max/min = ', tets.max(), tets.min())
        #bcs = model.bcs
        #mapbc = model.mapbc
        #loads = model.loads

        self.nNodes = nnodes
        self.nElements = ntris + ntets

        #print("nNodes = %i" % self.nNodes)
        #print("nElements = %i" % self.nElements)

        self.grid.Allocate(self.nElements, 1000)
        #self.gridResult.SetNumberOfComponents(self.nElements)

        points = vtk.vtkPoints()
        points.SetNumberOfPoints(self.nNodes)
        self.nid_map = {}
        if 0:
            fraction = 1. / self.nNodes  # so you can color the nodes by ID
            for nid, node in sorted(iteritems(nodes)):
                points.InsertPoint(nid - 1, *node)
                self.gridResult.InsertNextValue(nid * fraction)

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

        nid = 0
        for i in range(nnodes):
            points.InsertPoint(nid, nodes[i, :])
            nid += 1

        if dimension_flag == 2:
            for (n0, n1, n2) in tris:
                elem = vtkTriangle()
                #node_ids = elements[eid, :]
                elem.GetPointIds().SetId(0, n0)
                elem.GetPointIds().SetId(1, n1)
                elem.GetPointIds().SetId(2, n2)
                self.grid.InsertNextCell(5, elem.GetPointIds())  #elem.GetCellType() = 5  # vtkTriangle
        elif dimension_flag == 3:
            if ntets:
                for (n0, n1, n2, n3) in tets:
                    elem = vtkTetra()
                    elem.GetPointIds().SetId(0, n0)
                    elem.GetPointIds().SetId(1, n1)
                    elem.GetPointIds().SetId(2, n2)
                    elem.GetPointIds().SetId(3, n3)
                    self.grid.InsertNextCell(10, elem.GetPointIds())  #elem.GetCellType() = 5  # vtkTriangle
        else:
            raise RuntimeError('dimension_flag=%r' % dimension_flag)

        self.grid.SetPoints(points)
        self.grid.Modified()
        if hasattr(self.grid, 'Update'):
            self.grid.Update()

        # regions/loads
        self. turn_text_on()
        self.scalarBar.Modified()

        cases = {}
        #cases = self.result_cases
        self._fill_fast_results(cases, model, results=False)
        self._finish_results_io(cases)

  if type.GetElementTypeId() == elements.ELEMENT_QUAD:
    newNodes = copy.deepcopy(nodes)
    newNodes[2] = nodes[3]
    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):
    """

    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])
#.........这里部分代码省略.........

aHexahedronGrid.InsertNextCell(aHexahedron.GetCellType(),aHexahedron.GetPointIds())
aHexahedronGrid.SetPoints(hexahedronPoints)
aHexahedronMapper = vtk.vtkDataSetMapper()
aHexahedronMapper.SetInputData(aHexahedronGrid)
aHexahedronActor = vtk.vtkActor()
aHexahedronActor.SetMapper(aHexahedronMapper)
aHexahedronActor.AddPosition(2,0,0)
aHexahedronActor.GetProperty().BackfaceCullingOn()
# Tetra
tetraPoints = vtk.vtkPoints()
tetraPoints.SetNumberOfPoints(4)
tetraPoints.InsertPoint(0,0,0,0)
tetraPoints.InsertPoint(1,1,0,0)
tetraPoints.InsertPoint(2,0,1,0)
tetraPoints.InsertPoint(3,1,1,1)
aTetra = vtk.vtkTetra()
aTetra.GetPointIds().SetId(0,0)
aTetra.GetPointIds().SetId(1,1)
aTetra.GetPointIds().SetId(2,2)
aTetra.GetPointIds().SetId(3,3)
aTetraGrid = vtk.vtkUnstructuredGrid()
aTetraGrid.Allocate(1,1)
aTetraGrid.InsertNextCell(aTetra.GetCellType(),aTetra.GetPointIds())
aTetraGrid.SetPoints(tetraPoints)
aTetraMapper = vtk.vtkDataSetMapper()
aTetraMapper.SetInputData(aTetraGrid)
aTetraActor = vtk.vtkActor()
aTetraActor.SetMapper(aTetraMapper)
aTetraActor.AddPosition(4,0,0)
aTetraActor.GetProperty().BackfaceCullingOn()
# Wedge

    def load_tetgen_geometry(self, smesh_filename, dirname, plot=True):
        print("load_tetgen_geometry...")
        skipReading = self.removeOldGeometry(smesh_filename)
        if skipReading:
            return

        model = TetgenReader(log=self.log, debug=False)

        base_filename, ext = os.path.splitext(smesh_filename)
        node_filename = base_filename + '.node'
        ele_filename = base_filename + '.ele'
        if '.smesh' == ext:
            dimension_flag = 2
        elif '.ele' == ext:
            dimension_flag = 3
        else:
            raise RuntimeError('unsupported extension.  Use "smesh" or "ele".')
        model.read_tetgen(node_filename, smesh_filename, ele_filename, dimension_flag)
        nodes = model.nodes
        tris = model.tri
        tets = model.tet

        self.nNodes, three = nodes.shape
        ntris = 0
        ntets = 0
        if dimension_flag == 2:
            ntris, three = tris.shape
        elif dimension_flag == 3:
            ntets, four = tets.shape
        else:
            raise RuntimeError()
        self.nElements = ntris + ntets

        print("nNodes = ",self.nNodes)
        print("nElements = ", self.nElements)

        self.grid.Allocate(self.nElements, 1000)
        #self.gridResult.SetNumberOfComponents(self.nElements)
        self.grid2.Allocate(1, 1000)

        points = vtk.vtkPoints()
        points.SetNumberOfPoints(self.nNodes)
        #self.gridResult.Allocate(self.nNodes, 1000)
        #vectorReselt.SetNumberOfComponents(3)
        self.nidMap = {}
        #elem.SetNumberOfPoints(nNodes)
        if 0:
            fraction = 1. / self.nNodes  # so you can color the nodes by ID
            for nid, node in sorted(iteritems(nodes)):
                points.InsertPoint(nid - 1, *node)
                self.gridResult.InsertNextValue(nid * fraction)
                #print(str(element))

                #elem = vtk.vtkVertex()
                #elem.GetPointIds().SetId(0, i)
                #self.aQuadGrid.InsertNextCell(elem.GetCellType(), elem.GetPointIds())
                #vectorResult.InsertTuple3(0, 0.0, 0.0, 1.0)

        assert nodes is not None
        nnodes, three = nodes.shape

        nid = 0
        print("nnodes=%s" % nnodes)
        for i in range(nnodes):
            points.InsertPoint(nid, nodes[i, :])
            nid += 1

        #elements -= 1
        if dimension_flag == 2:
            for (n0, n1, n2) in tris:
                elem = vtkTriangle()
                #node_ids = elements[eid, :]
                elem.GetPointIds().SetId(0, n0)
                elem.GetPointIds().SetId(1, n1)
                elem.GetPointIds().SetId(2, n2)
                self.grid.InsertNextCell(5, elem.GetPointIds())  #elem.GetCellType() = 5  # vtkTriangle
        elif dimension_flag == 3:
            for (n0, n1, n2, n3) in tets:
                elem = vtkTetra()
                assert elem.GetCellType() == 10, elem.GetCellType()
                elem.GetPointIds().SetId(0, n0)
                elem.GetPointIds().SetId(1, n1)
                elem.GetPointIds().SetId(2, n2)
                elem.GetPointIds().SetId(3, n3)
                self.grid.InsertNextCell(10, elem.GetPointIds())  #elem.GetCellType() = 5  # vtkTriangle
        else:
            raise RuntimeError()

        self.grid.SetPoints(points)
        self.grid.Modified()
        self.grid.Update()
        print("updated grid")

        # loadSTLResults - regions/loads
        self.TurnTextOn()
        self.scalarBar.VisibilityOff()
        self.scalarBar.Modified()

        cases = {}
        ID = 1
#.........这里部分代码省略.........

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 mapElements(self, points, points2, nidMap, model, j):
        self.eidMap = {}
        i = 0
        for (eid, element) in sorted(model.elements.iteritems()):
            self.eidMap[eid] = i
            #print element.type
            if isinstance(element, CTRIA3) or isinstance(element, CTRIAR):
                #print "ctria3"
                elem = vtkTriangle()
                nodeIDs = element.nodeIDs()
                elem.GetPointIds().SetId(0, nidMap[nodeIDs[0]])
                elem.GetPointIds().SetId(1, nidMap[nodeIDs[1]])
                elem.GetPointIds().SetId(2, nidMap[nodeIDs[2]])
                self.grid.InsertNextCell(elem.GetCellType(),
                                         elem.GetPointIds())
            elif isinstance(element, CTRIA6):
                nodeIDs = element.nodeIDs()
                if None not in nodeIDs:
                    elem = vtkQuadraticTriangle()
                    elem.GetPointIds().SetId(3, nidMap[nodeIDs[3]])
                    elem.GetPointIds().SetId(4, nidMap[nodeIDs[4]])
                    elem.GetPointIds().SetId(5, nidMap[nodeIDs[5]])
                else:
                    elem = vtkTriangle()
                elem.GetPointIds().SetId(0, nidMap[nodeIDs[0]])
                elem.GetPointIds().SetId(1, nidMap[nodeIDs[1]])
                elem.GetPointIds().SetId(2, nidMap[nodeIDs[2]])
                self.grid.InsertNextCell(elem.GetCellType(),
                                         elem.GetPointIds())
            elif isinstance(element, CTRIAX6):
                # midside nodes are required, nodes out of order
                nodeIDs = element.nodeIDs()
                if None not in nodeIDs:
                    elem = vtkQuadraticTriangle()
                    elem.GetPointIds().SetId(3, nidMap[nodeIDs[1]])
                    elem.GetPointIds().SetId(4, nidMap[nodeIDs[3]])
                    elem.GetPointIds().SetId(5, nidMap[nodeIDs[5]])
                else:
                    elem = vtkTriangle()
                elem.GetPointIds().SetId(0, nidMap[nodeIDs[0]])
                elem.GetPointIds().SetId(1, nidMap[nodeIDs[2]])
                elem.GetPointIds().SetId(2, nidMap[nodeIDs[4]])
                #a = [0,2,4]
                #msg = "CTRIAX6 %i %i %i" %(nidMap[nodeIDs[a[0]]],
                #                           nidMap[nodeIDs[a[1]]],
                #                           nidMap[nodeIDs[a[2]]] )
                #raise RuntimeError(msg)
                #sys.stdout.flush()

                #elem.GetPointIds().SetId(0, nidMap[nodeIDs[0]])
                #elem.GetPointIds().SetId(1, nidMap[nodeIDs[1]])
                #elem.GetPointIds().SetId(2, nidMap[nodeIDs[2]])
                self.grid.InsertNextCell(elem.GetCellType(),
                                         elem.GetPointIds())

            elif (isinstance(element, CQUAD4) or isinstance(element, CSHEAR) or
                  isinstance(element, CQUADR)):
                nodeIDs = element.nodeIDs()
                elem = vtkQuad()
                elem.GetPointIds().SetId(0, nidMap[nodeIDs[0]])
                elem.GetPointIds().SetId(1, nidMap[nodeIDs[1]])
                elem.GetPointIds().SetId(2, nidMap[nodeIDs[2]])
                elem.GetPointIds().SetId(3, nidMap[nodeIDs[3]])
                self.grid.InsertNextCell(elem.GetCellType(),
                                         elem.GetPointIds())
            elif isinstance(element, CQUAD8):
                nodeIDs = element.nodeIDs()
                if None not in nodeIDs:
                    elem = vtkQuadraticQuad()
                    elem.GetPointIds().SetId(4, nidMap[nodeIDs[4]])
                    elem.GetPointIds().SetId(5, nidMap[nodeIDs[5]])
                    elem.GetPointIds().SetId(6, nidMap[nodeIDs[6]])
                    elem.GetPointIds().SetId(7, nidMap[nodeIDs[7]])
                else:
                    elem = vtkQuad()
                elem.GetPointIds().SetId(0, nidMap[nodeIDs[0]])
                elem.GetPointIds().SetId(1, nidMap[nodeIDs[1]])
                elem.GetPointIds().SetId(2, nidMap[nodeIDs[2]])
                elem.GetPointIds().SetId(3, nidMap[nodeIDs[3]])
                self.grid.InsertNextCell(elem.GetCellType(),
                                         elem.GetPointIds())
            elif isinstance(element, CTETRA4):
                elem = vtkTetra()
                nodeIDs = element.nodeIDs()
                elem.GetPointIds().SetId(0, nidMap[nodeIDs[0]])
                elem.GetPointIds().SetId(1, nidMap[nodeIDs[1]])
                elem.GetPointIds().SetId(2, nidMap[nodeIDs[2]])
                elem.GetPointIds().SetId(3, nidMap[nodeIDs[3]])
                self.grid.InsertNextCell(elem.GetCellType(),
                                         elem.GetPointIds())
            elif isinstance(element, CTETRA10):
                nodeIDs = element.nodeIDs()
                if None not in nodeIDs:
                    elem = vtkQuadraticTetra()
                    elem.GetPointIds().SetId(4, nidMap[nodeIDs[4]])
                    elem.GetPointIds().SetId(5, nidMap[nodeIDs[5]])
                    elem.GetPointIds().SetId(6, nidMap[nodeIDs[6]])
                    elem.GetPointIds().SetId(7, nidMap[nodeIDs[7]])
                    elem.GetPointIds().SetId(8, nidMap[nodeIDs[8]])
                    elem.GetPointIds().SetId(9, nidMap[nodeIDs[9]])
#.........这里部分代码省略.........

#.........这里部分代码省略.........
        elif dim == 2:
            for i in xrange(n):
                points.SetPoint(i, x[i], y[i], 0.0)
        else:
            for i in xrange(n):
                points.SetPoint(i, x[i], y[i], z[i])
        vtu.SetPoints(points)

        id_list = vtk.vtkIdList()
        if dim == 1:
            if degree in [0, 1]:
                cell_type = vtk.vtkLine().GetCellType()
                id_list.SetNumberOfIds(2)
                cell_map = None
            elif degree == 2:
                cell_type = vtk.vtkQuadraticEdge().GetCellType()
                id_list.SetNumberOfIds(3)
                cell_map = None
            else:
                cell_type = vtk.vtkCubicLine().GetCellType()
                id_list.SetNumberOfIds(4)
                cell_map = None
        elif dim == 2:
            if degree in [0, 1]:
                cell_type = vtk.vtkTriangle().GetCellType()
                id_list.SetNumberOfIds(3)
                cell_map = None
            else:
                cell_type = vtk.vtkQuadraticTriangle().GetCellType()
                id_list.SetNumberOfIds(6)
                cell_map = {0:0, 1:1, 2:2, 3:5, 4:3, 5:4}
        else:
            if degree in [0, 1]:
                cell_type = vtk.vtkTetra().GetCellType()
                id_list.SetNumberOfIds(4)
                cell_map = None
            else:
                cell_type = vtk.vtkQuadraticTetra().GetCellType()
                id_list.SetNumberOfIds(10)
                cell_map = {0:0, 1:1, 2:2, 3:3, 4:9, 5:6, 6:8, 7:7, 8:5, 9:4}
        for i in xrange(mesh.num_cells()):
            cell = xdof.cell_dofs(i)
            assert(len(cell) == id_list.GetNumberOfIds())
            if not cell_map is None:
                cell = [cell[cell_map[j]] for j in xrange(len(cell))]
            for j in xrange(len(cell)):
                id_list.SetId(j, xnode_map[cell[j]])
            vtu.InsertNextCell(cell_type, id_list)

        if degree == 0:
            for fn in fns[e]:
                if not fn.value_rank() == 0:
                    raise NotImplementException("Function rank %i not supported by write_vtu" % fn.value_rank())
                data = fn.vector().gather(nodes)
                cell_data = vtk.vtkDoubleArray()
                cell_data.SetNumberOfComponents(1)
                cell_data.SetNumberOfValues(mesh.num_cells())
                cell_data.SetName(fn.name())
                for i, datum in enumerate(data):
                    cell_data.SetValue(i, datum)
                vtu.GetCellData().AddArray(cell_data)
            vtu.GetCellData().SetActiveScalars(names[e][0])
        else:
            for fn in fns[e]:
                if not fn.value_rank() == 0:
                    raise NotImplementException("Function rank %i not supported by write_vtu" % fn.value_rank())

    def mapElements(self, points, points2, nidMap, model, j):
        #self.eidMap = {}
        i = 0
        for (eid, element) in sorted(model.elements.iteritems()):
            self.eidMap[eid] = i
            #print element.type
            if isinstance(element, CTRIA3) or isinstance(element, CTRIAR):
                elem = vtkTriangle()
                nod