meshconvert.py

利用C

        if sect == "node":
            # Read node definition
            m = re_node.match(l)
            if m is None:
                handler.warn("Node on line %d is on unsupported format" % (lineno,))
                continue
            idx, c0, c1, c2 = m.groups()
            try: coords = [float(c) for c in (c0, c1, c2)]
            except ValueError:
                handler.warn("Node on line %d contains non-numeric coordinates"
                        % (lineno,))
                continue
            nodes[int(idx)] = coords
        elif sect == "element":
            if not supported_elem:
                continue
            m = re_tetra.match(l)
            if m is None:
                handler.error("Node on line %d badly specified (expected 3 "
                        "coordinates)" % (lineno,))
            idx, n0, n1, n2, n3 = [int(x) for x in m.groups()]
            elems[idx] = (tp, n0, n1, n2, n3)
            eid2elset.setdefault(elset, set()).add(idx)

    ifile.close()

    # Note that vertices/cells must be consecutively numbered, which isn't
    # necessarily the case in Abaqus. Therefore we enumerate and translate
    # original IDs to sequence indexes.

    handler.start_vertices(len(nodes))
    nodeids = nodes.keys()
    nodeids.sort()
    for idx, nid in enumerate(nodeids):
        handler.add_vertex(idx, nodes[nid])
    handler.end_vertices()

    handler.start_cells(len(elems))
    elemids = elems.keys()
    elemids.sort()
    for idx, eid in enumerate(elemids):
        elem = elems[eid]
        tp = elem[0]
        elemnodes = []
        for nid in elem[1:]:
            try: elemnodes.append(nodeids.index(nid))
            except ValueError:
                handler.error("Element %s references non-existent node %s" % (eid, nid))
        handler.add_cell(idx, elemnodes)
    handler.end_cells()

    # Define the material function for the cells

    num_entities = 0
    for matname, elsetids in material2elsetids.items():
        if matname not in materials:
            handler.error("Unknown material %s referred to for element sets %s" %
                    (matname, ", ".join(elsetids)))
        num_entities += len(elsetids)
    handler.start_meshfunction("material", 3, num_entities)
    # Each material is associated with a number of element sets
    for i, matname in enumerate(materials):
        try: elsetids = material2elsetids[matname]
        except KeyError:
            # No elements for this material
            continue
        # For each element set associated with this material
        elsets = []
        for eid in elsetids:
            try: elsets.append(eid2elset[eid])
            except KeyError:
                handler.error("Material '%s' is assigned to undefined element "
                        "set '%s'" % (matname, eid))
        for elset in elsets:
            for elemid in elset:
                handler.add_entity_meshfunction(elemids.index(elemid), i)
    handler.end_meshfunction()

def netcdf2xml(ifilename,ofilename):
    "Convert from NetCDF format to DOLFIN XML."

    print "Converting from NetCDF format (.ncdf) to DOLFIN XML format"

    # Open files
    ifile = open(ifilename, "r")
    ofile = open(ofilename, "w")


    cell_type = None
    dim = 0

    # Scan file for dimension, number of nodes, number of elements
    while 1:
        line = ifile.readline()
        if not line:
            error("Empty file")
        if re.search(r"num_dim.*=", line):
            dim = int(re.match(".*\s=\s(\d+)\s;",line).group(1))
        if re.search(r"num_nodes.*=", line):
            num_vertices = int(re.match(".*\s=\s(\d+)\s;",line).group(1))
        if re.search(r"num_elem.*=", line):
            num_cells = int(re.match(".*\s=\s(\d+)\s;",line).group(1))
        if re.search(r"connect1 =",line):
            break            

    num_dims=dim
    
    # Set cell type
    if dim == 2:
        cell_type ="triangle"
    if dim == 3:
        cell_type ="tetrahedron"
        
    # Check that we got the cell type
    if cell_type == None:
        error("Unable to find cell type.")
            
    # Write header
    write_header_mesh(ofile, cell_type, dim)
    write_header_cells(ofile, num_cells)
    num_cells_read = 0

    # Read and write cells
    while 1:
        # Read next line
        line = ifile.readline()
        if not line: 
            break
        connect=re.split("[,;]",line)
        if num_dims == 2:
            n0 = int(connect[0])-1
            n1 = int(connect[1])-1
            n2 = int(connect[2])-1
            write_cell_triangle(ofile, num_cells_read, n0, n1, n2)
        elif num_dims == 3:
            n0 = int(connect[0])-1
            n1 = int(connect[1])-1
            n2 = int(connect[2])-1
            n3 = int(connect[3])-1
            write_cell_tetrahedron(ofile, num_cells_read, n0, n1, n2, n3)
        num_cells_read +=1
        if num_cells == num_cells_read:
           write_footer_cells(ofile)
           write_header_vertices(ofile, num_vertices)
           break

    num_vertices_read = 0
    coords = [[],[],[]]
    coord = -1

    while 1:
        line = ifile.readline()
        if not line:
            error("Missing data")
        if re.search(r"coord =",line):
            break            

    # Read vertices
    while 1:
        line = ifile.readline()
#        print line
        if not line: 
            break
        if re.search(r"\A\s\s\S+,",line):
#            print line
            coord+=1
            print "Found x_"+str(coord)+" coordinates"
        coords[coord] += line.split()
        if re.search(r";",line):
            break

    # Write vertices
    for i in range(num_vertices):
        if num_dims == 2:
            x = float(re.split(",",coords[0].pop(0))[0])
            y = float(re.split(",",coords[1].pop(0))[0])
            z = 0
        if num_dims == 3:
            x = float(re.split(",",coords[0].pop(0))[0])
            y = float(re.split(",",coords[1].pop(0))[0])
            z = float(re.split(",",coords[2].pop(0))[0])
        write_vertex(ofile, i, x, y, z)


    # Write footer
    write_footer_vertices(ofile)
    write_footer_mesh(ofile)

    # Close files
    ifile.close()
    ofile.close()

def exodus2xml(ifilename,ofilename):
    "Convert from Exodus II format to DOLFIN XML."

    print "Converting from Exodus II format to NetCDF format"

    name = ifilename.split(".")[0]
    netcdffilename = name +".ncdf"
    os.system('ncdump '+ifilename + ' > '+netcdffilename)
    netcdf2xml(netcdffilename,ofilename)

# Write mesh header
def write_header_mesh(ofile, cell_type, dim):
    ofile.write("""\
<?xml version=\"1.0\" encoding=\"UTF-8\"?>

<dolfin xmlns:dolfin=\"http://www.fenics.org/dolfin/\">
  <mesh celltype="%s" dim="%d">
""" % (cell_type, dim))

# Write graph header
def write_header_graph(ofile, graph_type):
    ofile.write("""\
<?xml version=\"1.0\" encoding=\"UTF-8\"?>

<dolfin xmlns:dolfin=\"http://www.fenics.org/dolfin/\">
  <graph type="%s">
""" % (graph_type))

# Write mesh footer
def write_footer_mesh(ofile):
    ofile.write("""\
  </mesh>
</dolfin>
""")

# Write graph footer
def write_footer_graph(ofile):
    ofile.write("""\
  </graph>
</dolfin>
""")

def write_header_vertices(ofile, num_vertices):
    "Write vertices header"
    print "Expecting %d vertices" % num_vertices
    ofile.write("    <vertices size=\"%d\">\n" % num_vertices)

def write_footer_vertices(ofile):
    "Write vertices footer"
    ofile.write("    </vertices>\n")
    print "Found all vertices"

def write_header_edges(ofile, num_edges):
    "Write edges header"
    print "Expecting %d edges" % num_edges
    ofile.write("    <edges size=\"%d\">\n" % num_edges)

def write_footer_edges(ofile):
    "Write edges footer"
    ofile.write("    </edges>\n")
    print "Found all edges"

def write_vertex(ofile, vertex, x, y, z):
    "Write vertex"
    ofile.write("      <vertex index=\"%d\" x=\"%s\" y=\"%s\" z=\"%s\"/>\n" % \
        (vertex, x, y, z))

def write_graph_vertex(ofile, vertex, num_edges, weight = 1):
    "Write graph vertex"
    ofile.write("      <vertex index=\"%d\" num_edges=\"%d\" weight=\"%d\"/>\n" % \
        (vertex, num_edges, weight))

def write_graph_edge(ofile, v1, v2, weight = 1):
	 "Write graph edge"
	 ofile.write("      <edge v1=\"%d\" v2=\"%d\" weight=\"%d\"/>\n" % \
        (v1, v2, weight))

def write_header_cells(ofile, num_cells):
    "Write cells header"
    ofile.write("    <cells size=\"%d\">\n" % num_cells)
    print "Expecting %d cells" % num_cells

def write_footer_cells(ofile):
    "Write cells footer"
    ofile.write("    </cells>\n")
    print "Found all cells"

def write_cell_triangle(ofile, cell, n0, n1, n2):
    "Write cell (triangle)"
    ofile.write("      <triangle index=\"%d\" v0=\"%d\" v1=\"%d\" v2=\"%d\"/>\n" % \
        (cell, n0, n1, n2))

def write_cell_tetrahedron(ofile, cell, n0, n1, n2, n3):
    "Write cell (tetrahedron)"
    ofile.write("      <tetrahedron index=\"%d\" v0=\"%d\" v1=\"%d\" v2=\"%d\" v3=\"%d\"/>\n" % \
        (cell, n0, n1, n2, n3))

def _error(message):
    "Write an error message"
    for line in message.split("\n"):
        print "*** %s" % line
    sys.exit(2)

def convert2xml(ifilename, ofilename, iformat=None):
    """ Convert a file to the DOLFIN XML format.
    """
    convert(ifilename, XmlHandler(ofilename), iformat=iformat)

def convert(ifilename, handler, iformat=None):
    """ Convert a file using a provided data handler.

    Note that handler.close is called when this function finishes.
    @param ifilename: Name of input file.
    @param handler: The data handler (instance of L{DataHandler}).
    @param iformat: Format of input file.
    """
    if iformat is None:
        iformat = format_from_suffix(os.path.splitext(ifilename)[1][1:])
    # XXX: Backwards-compat
    if hasattr(handler, "_ofilename"):
        ofilename = handler._ofilename
    # Choose conversion
    if iformat == "mesh":
        # Convert from mesh to xml format
        mesh2xml(ifilename, ofilename)
    elif iformat == "gmsh":
        # Convert from gmsh to xml format
        gmsh2xml(ifilename, ofilename)
    elif iformat == "xml-old":
        # Convert from old to new xml format
        xml_old2xml(ifilename, ofilename)
    elif iformat == "metis":
        # Convert from metis graph to dolfin graph xml format
        metis_graph2graph_xml(ifilename, ofilename)
    elif iformat == "scotch":
        # Convert from scotch graph to dolfin graph xml format
        scotch_graph2graph_xml(ifilename, ofilename)
    elif iformat == "diffpack":
        # Convert from Diffpack tetrahedral grid format to xml format
        diffpack2xml(ifilename, ofilename)
    elif iformat == "abaqus":
        # Convert from abaqus to xml format
        _abaqus(ifilename, handler)
    elif iformat == "NetCDF":
        # Convert from NetCDF generated from ExodusII format to xml format
        netcdf2xml(ifilename, ofilename)      
    elif iformat =="ExodusII":
        # Convert from ExodusII format to xml format via NetCDF
        exodus2xml(ifilename, ofilename)
    elif iformat == "StarCD": 
        # Convert from Star-CD tetrahedral grid format to xml format
        starcd2xml(ifilename, ofilename)
    else:
        _error("Sorry, cannot convert between %s and DOLFIN xml file formats." % iformat)

    # XXX: handler.close messes things for other input formats than abaqus
    if iformat == "abaqus":
        handler.close()

def starcd2xml(ifilename, ofilename):
    "Convert from Star-CD tetrahedral grid format to DOLFIN XML."

    print starcd2xml.__doc__

    if not os.path.isfile(ifilename[:-3] + "vrt") or not os.path.isfile(ifilename[:-3] + "cel"):
        print "StarCD format requires one .vrt file and one .cel file"
        sys.exit(2)
      

    # open output file 
    ofile = open(ofilename, "w")

    # Open file, the vertices are in a .vrt file
    ifile = open(ifilename[:-3] + "vrt", "r")

    write_header_mesh(ofile, "tetrahedron", 3)


    # Read & write vertices

    # first, read all lines (need to sweep to times through the file)
    lines = ifile.readlines()

    # second, find the number of vertices
    num_vertices = -1 
    counter = 0 
    # nodenr_map is needed because starcd support node numbering like 1,2,4 (ie 3 is missing)  
    nodenr_map = {} 
    for line in lines: 
        nodenr = int(line[0:15])
        nodenr_map[nodenr] = counter
        counter += 1 
    num_vertices = counter
    
    # third, run over all vertices
    write_header_vertices(ofile, num_vertices)
    for line in lines: 
        nodenr = int(line[0:15])
        vertex0 = float(line[15:31])
        vertex1 = float(line[31:47])
        vertex2 = float(line[47:63])
        write_vertex(ofile, nodenr_map[nodenr], float(vertex0), float(vertex1), float(vertex2))  
    write_footer_vertices(ofile)

    # Open file, the cells are in a .cel file
    ifile = open(ifilename[:-3] + "cel", "r")

    # Read & write cells 

    # first, read all lines (need to sweep to times through the file)
    lines = ifile.readlines()

    # second, find the number of cells 
    num_cells = -1 
    counter = 0 
    for line in lines: 
        l = [int(a) for a in line.split()]
        cellnr, node0, node1, node2, node3, node4, node5, node6, node7, tmp1, tmp2  = l 
	if node4 > 0:
        	if node2 == node3 and node4 == node5 and node5 == node6 and node6 == node7: # these nodes should be equal  
                	counter += 1
		else:            
			print "The file does contain cells that are not tetraheders. The cell number is ", cellnr, " the line read was ", line                      
        else:                              
            # triangles on the surface 
#            print "The file does contain cells that are not tetraheders node4==0. The cell number is ", cellnr, " the line read was ", line 
            #sys.exit(2)
            pass
 
    num_cells = counter

    # third, run over all cells 
    write_header_cells(ofile, num_cells)
    counter = 0 
    for line in lines: 
        l = [int(a) for a in line.split()]
        cellnr, node0, node1, node2, node3, node4, node5, node6, node7, tmp1, tmp2  = l 
        if (node4 > 0): 
	        if node2 == node3 and node4 == node5 and node5 == node6 and node6 == node7: # these nodes should be equal  
	
			write_cell_tetrahedron(ofile, counter, nodenr_map[node0], nodenr_map[node1], nodenr_map[node2], nodenr_map[node4])
          		counter += 1 

       
    write_footer_cells(ofile)
    write_footer_mesh(ofile)


    # Close files
    ifile.close()
    ofile.close()