def quadsActor(bounds, color):
    """Create solid, axis-aligned quads at 0 in Z.

    Args:
      bounds: [[[xmin, ymin], [xmax, ymax]], ...]
      color: [R, G, B, A]
    """
    points = vtk.vtkPoints()
    quads = vtk.vtkCellArray()
    colors = vtk.vtkUnsignedCharArray()
    colors.SetNumberOfComponents(4)
    for (index, bound) in enumerate(bounds):
        colors.InsertNextTuple4(*color)
        (low, high) = bound
        points.InsertNextPoint(low[0], low[1], 0)
        points.InsertNextPoint(high[0], low[1], 0)
        points.InsertNextPoint(high[0], high[1], 0)
        points.InsertNextPoint(low[0], high[1], 0)
        quad = vtk.vtkQuad()
        for i in range(4):
            quad.GetPointIds().SetId(i, 4 * index + i)
        quads.InsertNextCell(quad)
    poly_data = vtk.vtkPolyData()
    poly_data.SetPoints(points)
    poly_data.SetPolys(quads)
    poly_data.GetCellData().SetScalars(colors)
    mapper = vtk.vtkPolyDataMapper()
    set_mapper_input(mapper, poly_data)
    actor = vtk.vtkActor()
    actor.SetMapper(mapper)
    actor.RotateZ(180)
    actor.RotateY(180)
    return actor

  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 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 _create_tecplot_shells(self, is_quads, quads, is_tris, tris):
     if is_quads:
         elements = quads
         for iface, face in enumerate(quads):
             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
             #break
     if is_tris:
         elements = tris
         for iface, face in enumerate(tris):
             elem = vtkTriangle()
             epoints = elem.GetPointIds()
             epoints.SetId(0, face[0])
             epoints.SetId(1, face[1])
             epoints.SetId(2, face[2])
             self.grid.InsertNextCell(5, elem.GetPointIds())  #elem.GetCellType() = 5  # vtkTriangle

  def load_point_cloud(self, points, colors=None, mesh=False):
    # first, generate quad mesh, if desired
    if mesh:
      mesh = vtk.vtkCellArray()
      if len(points.shape) != 3 or points.shape[2] != 3:
        raise Exception('For the meshing to work, the points array must have ' +
          'the form MxNx3')

      for i in xrange(points.shape[0] - 1):
        offset = i * points.shape[1]
        for j in xrange(points.shape[1] - 1):
          quad = vtk.vtkQuad()
          pids = quad.GetPointIds()
          pids.SetNumberOfIds(4)
          idx = offset + j
          pids.SetId(0, idx)
          pids.SetId(1, idx + 1)
          pids.SetId(2, idx + points.shape[1] + 1)
          pids.SetId(3, idx + points.shape[1])
          mesh.InsertNextCell(quad)

    poly_data = vtk.vtkPolyData()
    poly_data.SetPoints(self._convert_points_array(points))

    if mesh is not False:
      poly_data.SetPolys(mesh)
    else: # no mesh data; only show vertices
      vertices = vtk.vtkCellArray()
      vertices.InsertNextCell(points.shape[0], np.arange(points.shape[0]))
      poly_data.SetVerts(vertices)

    # Visualize
    actor = self._actor_from_poly_data(poly_data)
    actor.GetProperty().SetPointSize(self.point_size)

    if colors is not None:
      self._add_color_to_actor(actor, colors)

    return self._add_mesh_actor(actor)

    def set_quad_grid(self, name, nodes, elements, color, line_width=5, opacity=1.):
        """
        Makes a CQUAD4 grid
        """
        self.create_alternate_vtk_grid(name, color=color, line_width=line_width,
                                       opacity=opacity, representation='wire')

        nnodes = nodes.shape[0]
        nquads = elements.shape[0]
        #print(nodes)
        if nnodes == 0:
            return
        if nquads == 0:
            return

        #print('adding quad_grid %s; nnodes=%s nquads=%s' % (name, nnodes, nquads))
        points = vtk.vtkPoints()
        points.SetNumberOfPoints(nnodes)
        for nid, node in enumerate(nodes):
            #print(nid, node)
            points.InsertPoint(nid, *list(node))

        #assert vtkQuad().GetCellType() == 9, elem.GetCellType()
        self.alt_grids[name].Allocate(nquads, 1000)
        for element in elements:
            elem = vtk.vtkQuad()
            point_ids = elem.GetPointIds()
            point_ids.SetId(0, element[0])
            point_ids.SetId(1, element[1])
            point_ids.SetId(2, element[2])
            point_ids.SetId(3, element[3])
            self.alt_grids[name].InsertNextCell(9, elem.GetPointIds())
        self.alt_grids[name].SetPoints(points)

        self._add_alt_actors({name : self.alt_grids[name]})

        #if name in self.geometry_actors:
        self.geometry_actors[name].Modified()

def quad_poly_data(width, height):
    from vtk import vtkPoints, vtkQuad, vtkCellArray, vtkPolyData

    points = vtkPoints()
    points.InsertNextPoint((-1 * int(width / 2.0), int(height / 2.0), 0))
    points.InsertNextPoint((int(width / 2.0), int(height / 2.0), 0))
    points.InsertNextPoint((int(width / 2.0), -1 * int(height / 2.0), 0))
    points.InsertNextPoint((-1 * int(width / 2.0), -1 * int(height / 2.0), 0))

    quad = vtkQuad()
    quad.GetPointIds().SetId(0, 0)
    quad.GetPointIds().SetId(1, 1)
    quad.GetPointIds().SetId(2, 2)
    quad.GetPointIds().SetId(3, 3)

    arr = vtkCellArray()
    arr.InsertNextCell(quad)

    pd = vtkPolyData()
    pd.SetPoints(points)
    pd.SetPolys(arr)

    return pd

    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 load_degen_geom_geometry(self, csv_filename, dirname, name='main', plot=True):
        #key = self.case_keys[self.icase]
        #case = self.result_cases[key]

        skip_reading = self._remove_old_adb_geometry(csv_filename)
        if skip_reading:
            return

        model = DegenGeom(log=self.log, debug=False)
        self.model_type = 'vspaero'
        #self.model_type = model.model_type
        model.read_degen_geom(csv_filename)
        for name, comps in sorted(model.components.items()):
            print('name = %r' % name)
            #print(comp)
            print('------------')
            for comp in comps:
                nodes = comp.xyz
                elements = comp.elements
                nnodes = nodes.shape[0]
                nelements = elements.shape[0]

        self.nNodes = nnodes
        self.nElements = nelements

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

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

        assert nodes is not None

        nid = 0
        #print("nxyz_nodes=%s" % nxyz_nodes)
        mmax = amax(nodes, axis=0)
        mmin = amin(nodes, axis=0)
        dim_max = (mmax - mmin).max()
        self.create_global_axes(dim_max)

        for i in range(nnodes):
            points.InsertPoint(nid, nodes[i, :])
            nid += 1
        #self.log.info('nxyz_nodes=%s nwake_nodes=%s total=%s' % (
            #nnodes, nwake_nodes, nxyz_nodes + nwake_nodes))
        #self.log.info('nxyz_elements=%s nwake_elements=%s total=%s' % (
            #nxyz_elements, nwake_elements, nxyz_elements + nwake_elements))

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

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

        # load results - regions/loads
        self. turn_text_on()
        self.scalarBar.VisibilityOn()
        self.scalarBar.Modified()

        #mach = model.machs[0]
        #alpha = model.alphas[0]
        #beta = model.betas[0]
        #note = ':  Mach=%.2f, alpha=%.1f, beta=%.1f' % (mach, alpha, beta)
        note = 'name=%s' % name
        self.iSubcaseNameMap = {1: ['OpenVSP%s' % note, '']}
        cases = {}
        ID = 1

        form, cases = self._fill_degen_geom_case(cases, ID, model, nnodes, nelements)
        self._finish_results_io2(form, cases)

    def load_surf_geometry(self, surf_filename, dirname, plot=True):
        #skip_reading = self.remove_old_openfoam_geometry(openfoam_filename)
        #if skip_reading:
        #    return

        model = SurfReader()

        self.model_type = 'surf'
        print('surf_filename = %s' % surf_filename)

        model.read_surf(surf_filename)
        nnodes = model.nodes.shape[0]
        ntris = model.tris.shape[0]
        nquads = model.quads.shape[0]
        nelements = ntris + nquads

        nodes = model.nodes
        self.nElements = nelements
        self.nNodes = nnodes

        print("nNodes = %s" % self.nNodes)
        print("nElements = %s" % self.nElements)
        assert nelements > 0, nelements

        self.grid.Allocate(self.nElements, 1000)

        points = vtk.vtkPoints()
        points.SetNumberOfPoints(self.nNodes)

        mmax = amax(nodes, axis=0)
        mmin = amin(nodes, axis=0)
        dim_max = (mmax - mmin).max()
        self.create_global_axes(dim_max)
        self.log.info('max = %s' % mmax)
        self.log.info('min = %s' % mmin)

        for inode, node in enumerate(nodes):
            points.InsertPoint(inode, node)

        tris = model.tris - 1
        quads = model.quads - 1


        if ntris:
            for eid, element in enumerate(tris):
                elem = vtkTriangle()
                elem.GetPointIds().SetId(0, element[0])
                elem.GetPointIds().SetId(1, element[1])
                elem.GetPointIds().SetId(2, element[2])
                self.grid.InsertNextCell(elem.GetCellType(),
                                         elem.GetPointIds())
        if nquads:
            for eid, element in enumerate(quads):
                elem = vtkQuad()
                elem.GetPointIds().SetId(0, element[0])
                elem.GetPointIds().SetId(1, element[1])
                elem.GetPointIds().SetId(2, element[2])
                elem.GetPointIds().SetId(3, element[3])
                self.grid.InsertNextCell(elem.GetCellType(),
                                         elem.GetPointIds())

        model.read_surf_failnode(surf_filename)
        if len(model.nodes_failed):
            if 'failed_nodes' not in self.alt_grids:
                yellow = (1., 1., 0.)
                self.create_alternate_vtk_grid('failed_nodes', color=yellow, line_width=3, opacity=1.0)

            ifailed = where(model.nodes_failed == 1)[0]
            nfailed = len(ifailed)
            self.alt_grids['failed_nodes'].Allocate(nfailed, 1000)
            grid2 = self.alt_grids['failed_nodes']
            points2 = vtk.vtkPoints()
            points2.SetNumberOfPoints(nfailed)

            for j, nid in enumerate(model.nodes_failed):
                elem = vtk.vtkVertex()
                c = nodes[nid - 1, :]
                print(nid, c)
                points2.InsertPoint(j, *c)
                elem.GetPointIds().SetId(0, j)
                self.alt_grids['failed_nodes'].InsertNextCell(elem.GetCellType(), elem.GetPointIds())
            self.alt_grids['failed_nodes'].SetPoints(points2)
            self._add_alt_actors(self.alt_grids)

            actor = self.geometry_actors['failed_nodes']
            actor.Modified()
            prop = actor.GetProperty()
            prop.SetRepresentationToPoints()
            prop.SetPointSize(10)


            # self.

        self.nElements = nelements
        self.grid.SetPoints(points)
        self.grid.Modified()
        if hasattr(self.grid, 'Update'):
            self.grid.Update()
        #print("updated grid")

#.........这里部分代码省略.........

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 buildMesh(xNumCells, yNumCells, filename):
    
    polydata = vtk.vtkPolyData()
    points = vtk.vtkPoints()
    quads = vtk.vtkCellArray()

    counter = 0 
    
    xBase = 0
    yBase = 0
    
    xStep = quadLength/(xNumCells * 1.0) # Focing float division (python 2).
    yStep = quadHeight/(yNumCells * 1.0)
    
    branchId = vtk.vtkDoubleArray()
    branchId.SetName("branchId")

    
    # Three rectangles in space for a bifurcation.
    if bifurcation == True:
        stopAt = 3
    else:
        stopAt = 1

    for k in range(0, stopAt):
        if k == 1:
            xBase = - quadLength * (xQuads / 2)
            yBase = quadHeight * yQuads
            
        elif k == 2:
            xBase = + quadLength * (xQuads / 2) 
            yBase = quadHeight * yQuads     
                       
        # The framework for each quad. The total number is set by globals.
        for i in range(0, yQuads):
            for j in range(0, xQuads):
                
                    # For how ever many cells are within each main quad.
                    # The total number is decided by the function parameters.
                    for x in range(0, yNumCells):
                        for y in range(0, xNumCells):
                            
                            quad = vtk.vtkQuad()
                                                        
                            p0 = [y * xStep + xBase + (quadLength * j), 
                                  x * yStep + yBase + (quadHeight * i), 0]
                                  
                            p1 = [(y + 1) * xStep + xBase + (quadLength * j),
                                  x * yStep + yBase + (quadHeight * i), 0]
                                  
                            p2 = [(y + 1) * xStep + xBase + (quadLength * j), 
                                  (x + 1) * yStep + yBase + (quadHeight * i), 0]
                                  
                            p3 = [y * xStep + xBase + (quadLength * j),
                                  (x + 1) * yStep + yBase + (quadHeight * i), 0]
                                  
                            
                            points.InsertNextPoint(p0)
                            points.InsertNextPoint(p1)
                            points.InsertNextPoint(p2)
                            points.InsertNextPoint(p3)
                            
                            quad.GetPointIds().InsertId(0, counter)
                            quad.GetPointIds().InsertId(1, counter + 1)
                            quad.GetPointIds().InsertId(2, counter + 2)
                            quad.GetPointIds().InsertId(3, counter + 3)
                            counter += 4
                            
                            quads.InsertNextCell(quad)
                            branchId.InsertNextValue(k)

    polydata.SetPoints(points)
    polydata.SetPolys(quads)
    polydata.GetCellData().SetScalars(branchId)
    
    polyDataWriter = vtk.vtkXMLPolyDataWriter()
    polyDataWriter.SetFileName(filename)
    polyDataWriter.SetInputData(polydata)    
    polyDataWriter.Write()

    return polydata    

#.........这里部分代码省略.........
            #print(nodes)
            for inode, node in enumerate(nodes):
                points.InsertPoint(inode, node)

        #elements -= 1
        normals = None
        if is_3d_blockmesh:
            nelements, three = hexas.shape
            for eid, element in enumerate(hexas):
                #print(element)
                elem = vtkHexahedron()
                elem.GetPointIds().SetId(0, element[0])
                elem.GetPointIds().SetId(1, element[1])
                elem.GetPointIds().SetId(2, element[2])
                elem.GetPointIds().SetId(3, element[3])
                elem.GetPointIds().SetId(4, element[4])
                elem.GetPointIds().SetId(5, element[5])
                elem.GetPointIds().SetId(6, element[6])
                elem.GetPointIds().SetId(7, element[7])
                self.grid.InsertNextCell(elem.GetCellType(),
                                         elem.GetPointIds())

                #elem = vtkTriangle()
                #node_ids = elements[eid, :]
                #elem.GetPointIds().SetId(0, node_ids[0])
                #elem.GetPointIds().SetId(1, node_ids[1])
                #elem.GetPointIds().SetId(2, node_ids[2])

                #elem.GetCellType() = 5  # vtkTriangle
                #self.grid.InsertNextCell(5, elem.GetPointIds())
        elif is_surface_blockmesh:
            nelements, four = quads.shape
            for eid, element in enumerate(quads):
                elem = vtkQuad()
                elem.GetPointIds().SetId(0, element[0])
                elem.GetPointIds().SetId(1, element[1])
                elem.GetPointIds().SetId(2, element[2])
                elem.GetPointIds().SetId(3, element[3])
                self.grid.InsertNextCell(elem.GetCellType(),
                                         elem.GetPointIds())
        elif is_face_mesh:
            elems = quads
            nelements = quads.shape[0]
            nnames = len(names)
            normals = zeros((nelements, 3), dtype='float32')
            if nnames != nelements:
                msg = 'nnames=%s nelements=%s names.max=%s names.min=%s' % (
                    nnames, nelements, names.max(), names.min())
                raise RuntimeError(msg)
            for eid, element in enumerate(elems):
                #print('element = %s' % element)
                ineg = where(element == -1)[0]

                nnodes = 4
                if ineg:
                    nnodes = ineg.max()

                #pid = 1
                pid = names[eid]
                if nnodes == 3: # triangle!
                    f.write('CTRIA3,%i,%i,%i,%i,%i\n' % (
                        eid+1, pid, element[0]+1, element[1]+1, element[2]+1))
                    elem = vtkTriangle()
                    a = nodes[element[1], :] - nodes[element[0], :]
                    b = nodes[element[2], :] - nodes[element[0], :]
                    n = cross(a, b)

    def load_ugrid_geometry(self, ugrid_filename, dirname, name='main', plot=True):
        #skip_reading = self.remove_old_openfoam_geometry(openfoam_filename)
        #if skip_reading:
        #    return
        if is_binary_file(ugrid_filename):
            model = UGRID(log=self.log, debug=True)
            base, fmt, ext = os.path.basename(ugrid_filename).split('.')
            is_2d = False
        else:
            base, ext = os.path.basename(ugrid_filename).split('.')
            model = UGRID2D_Reader(log=self.log, debug=True)
            is_2d = True

        self.model_type = 'ugrid'
        print('ugrid_filename = %s' % ugrid_filename)


        assert ext == 'ugrid', ugrid_filename
        model.read_ugrid(ugrid_filename)

        if is_2d:
            tris = model.tris
            quads = model.quads
        else:
            tris = model.tris - 1
            quads = model.quads - 1

        #self.nodes = nodes
        #self.tris  = tris
        #self.quads = quads
        #self.pids = pids

        #self.tets = tets
        #self.penta5s = penta5s
        #self.penta6s = penta6s
        #self.hexas = hexas

        nnodes = model.nodes.shape[0]
        ntris = model.tris.shape[0]
        nquads = model.quads.shape[0]
        nelements = ntris + nquads

        nodes = model.nodes
        self.nElements = nelements
        self.nNodes = nnodes

        print("nNodes = %s" % self.nNodes)
        print("nElements = %s" % self.nElements)
        assert nelements > 0, nelements

        self.grid.Allocate(self.nElements, 1000)

        points = vtk.vtkPoints()
        points.SetNumberOfPoints(self.nNodes)

        mmax = amax(nodes, axis=0)
        mmin = amin(nodes, axis=0)
        dim_max = (mmax - mmin).max()
        self.create_global_axes(dim_max)
        self.log.info('max = %s' % mmax)
        self.log.info('min = %s' % mmin)

        diff_node_ids = model.check_hanging_nodes(stop_on_diff=False)
        if len(diff_node_ids):
            red = (1., 0., 0.)
            self.create_alternate_vtk_grid('hanging_nodes', color=red, line_width=5, opacity=1., point_size=10, representation='point')
            self._add_ugrid_nodes_to_grid('hanging_nodes', diff_node_ids, nodes)
            self._add_alt_actors(self.alt_grids)


        for inode, node in enumerate(nodes):
            points.InsertPoint(inode, node)

        if ntris:
            for eid, element in enumerate(tris):
                elem = vtkTriangle()
                elem.GetPointIds().SetId(0, element[0])
                elem.GetPointIds().SetId(1, element[1])
                elem.GetPointIds().SetId(2, element[2])
                self.grid.InsertNextCell(elem.GetCellType(),
                                         elem.GetPointIds())
        if nquads:
            for eid, element in enumerate(quads):
                elem = vtkQuad()
                elem.GetPointIds().SetId(0, element[0])
                elem.GetPointIds().SetId(1, element[1])
                elem.GetPointIds().SetId(2, element[2])
                elem.GetPointIds().SetId(3, element[3])
                self.grid.InsertNextCell(elem.GetCellType(),
                                         elem.GetPointIds())

        self.nElements = nelements
        self.grid.SetPoints(points)
        self.grid.Modified()
        print('update...')
        if hasattr(self.grid, 'Update'):
            self.grid.Update()
        #print("updated grid")

        # loadCart3dResults - regions/loads
#.........这里部分代码省略.........

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

#.........这里部分代码省略.........
                n_cohax4 += part.cohax4.shape[0]
            if part.cax3 is not None:
                n_cax3 += part.cax3.shape[0]
            if part.cax4r is not None:
                n_cax4r += part.cax4r.shape[0]

            # solids
            if part.c3d10h is not None:
                n_c3d10h += part.c3d10h.shape[0]

            if n_r2d2:
                eids = part.r2d2[:, 0]
                node_ids = part.r2d2[:, 1:] + nid_offset
                for eid, node_ids in zip(eids, node_ids):
                    elem = vtkLine()
                    elem.GetPointIds().SetId(0, node_ids[0])
                    elem.GetPointIds().SetId(1, node_ids[1])
                    self.grid.InsertNextCell(elem.GetCellType(), elem.GetPointIds())

            if n_cpe3:
                eids = part.cpe3[:, 0]
                node_ids = part.cpe3[:, 1:] + nid_offset
                for eid, node_ids in zip(eids, node_ids):
                    elem = vtkTriangle()
                    elem.GetPointIds().SetId(0, node_ids[0])
                    elem.GetPointIds().SetId(1, node_ids[1])
                    elem.GetPointIds().SetId(2, node_ids[2])
                    self.grid.InsertNextCell(5, elem.GetPointIds())

            if n_cpe4:
                eids = part.cpe4[:, 0]
                node_ids = part.cpe4[:, 1:] + nid_offset
                for eid, node_ids in zip(eids, node_ids):
                    elem = vtkQuad()
                    elem.GetPointIds().SetId(0, node_ids[0])
                    elem.GetPointIds().SetId(1, node_ids[1])
                    elem.GetPointIds().SetId(2, node_ids[2])
                    elem.GetPointIds().SetId(3, node_ids[3])
                    self.grid.InsertNextCell(elem.GetCellType(), elem.GetPointIds())

            if n_cpe4r:
                eids = part.cpe4r[:, 0]
                node_ids = part.cpe4r[:, 1:] + nid_offset
                for eid, node_ids in zip(eids, node_ids):
                    elem = vtkQuad()
                    elem.GetPointIds().SetId(0, node_ids[0])
                    elem.GetPointIds().SetId(1, node_ids[1])
                    elem.GetPointIds().SetId(2, node_ids[2])
                    elem.GetPointIds().SetId(3, node_ids[3])
                    self.grid.InsertNextCell(elem.GetCellType(), elem.GetPointIds())

            if n_coh2d4:
                eids = part.coh2d4[:, 0]
                node_ids = part.coh2d4[:, 1:] + nid_offset
                for eid, node_ids in zip(eids, node_ids):
                    elem = vtkQuad()
                    elem.GetPointIds().SetId(0, node_ids[0])
                    elem.GetPointIds().SetId(1, node_ids[1])
                    elem.GetPointIds().SetId(2, node_ids[2])
                    elem.GetPointIds().SetId(3, node_ids[3])
                    self.grid.InsertNextCell(elem.GetCellType(), elem.GetPointIds())

            if n_cohax4:
                eids = part.cohax4[:, 0]
                node_ids = part.cohax4[:, 1:] + nid_offset
                for eid, node_ids in zip(eids, node_ids):

    def load_panair_geometry(self, panairFileName, dirname, plot=True):
        self.nidMap = {}

        #key = self.caseKeys[self.iCase]
        #case = self.resultCases[key]

        skipReading = self.removeOldGeometry(panairFileName)
        if skipReading:
            return

        model = PanairGrid(log=self.log, debug=self.debug)
        self.modelType = model.modelType
        model.read_panair(panairFileName)

        nodes, elements, regions = model.getPointsElementsRegions()
        #for nid,node in enumerate(nodes):
            #print "node[%s] = %s" %(nid,str(node))

        self.nNodes = len(nodes)
        self.nElements = len(elements)

        #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)
        #elem.SetNumberOfPoints(nNodes)
        if 0:
            fraction = 1. / nNodes  # so you can color the nodes by ID
            for nid, node in sorted(iteritems(nodes)):
                points.InsertPoint(nid - 1, *point)
                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 len(nodes) > 0
        mmax = amax(nodes, axis=0)
        mmin = amin(nodes, axis=0)
        dim_max = (mmax - mmin).max()
        self.update_axes_length(dim_max)
        for nid, node in enumerate(nodes):
            points.InsertPoint(nid, *node)

        assert len(elements) > 0
        for eid, element in enumerate(elements):
            (p1, p2, p3, p4) = element
            #print "element = ",element
            elem = vtkQuad()
            elem.GetPointIds().SetId(0, p1)
            elem.GetPointIds().SetId(1, p2)
            elem.GetPointIds().SetId(2, p3)
            elem.GetPointIds().SetId(3, p4)
            self.grid.InsertNextCell(elem.GetCellType(), elem.GetPointIds())

        #print("eid = ", eid)
        self.grid.SetPoints(points)
        #self.grid2.SetPoints(points2)
        #self.grid.GetPointData().SetScalars(self.gridResult)
        #print dir(self.grid) #.SetNumberOfComponents(0)
        #self.grid.GetCellData().SetNumberOfTuples(1);
        #self.grid.GetCellData().SetScalars(self.gridResult)
        self.grid.Modified()
        self.grid2.Modified()
        if hasattr(self.grid, 'Update'):
            self.grid.Update()
            self.grid2.Update()
            print("updated grid")

        #return

        # loadCart3dResults - regions/loads
        self.TurnTextOn()
        self.scalarBar.VisibilityOn()
        self.scalarBar.Modified()

        self.iSubcaseNameMap = {1: ['Panair', '']}
        cases = {}
        ID = 1

        #print "nElements = ",nElements
        loads = []
        cases = self.fillPanairGeometryCase(cases, ID, nodes, elements, regions, loads)
        self._finish_results_io(cases)