meshconvert.py

利用C

        # Modify line
        if "xmlns" in line:
            line = "<dolfin xmlns:dolfin=\"http://www.fenics.org/dolfin/\">\n"
        if "<mesh>" in line:
            line = "  <mesh celltype=\"%s\" dim=\"%d\">\n" % (cell_type, dim)
        if dim == 2 and " z=\"0.0\"" in line:
            line = line.replace(" z=\"0.0\"", "")
        if " name=" in line:
            line = line.replace(" name=", " index=")
        if " name =" in line:
            line = line.replace(" name =", " index=")
        if "n0" in line:
            line = line.replace("n0", "v0")
        if "n1" in line:
            line = line.replace("n1", "v1")
        if "n2" in line:
            line = line.replace("n2", "v2")
        if "n3" in line:
            line = line.replace("n3", "v3")

        # Write line
        ofile.write(line)

    # Close files
    ifile.close();
    ofile.close();
    print "Conversion done"

def metis_graph2graph_xml(ifilename, ofilename):
    "Convert from Metis graph format to DOLFIN Graph XML."

    print "Converting from Metis graph format to DOLFIN Graph XML."
    
    # Open files
    ifile = open(ifilename, "r")
    ofile = open(ofilename, "w")

    # Read number of vertices and edges
    line = ifile.readline()
    if not line:
       _error("Empty file")

    (num_vertices, num_edges) = line.split()

    write_header_graph(ofile, "undirected")
    write_header_vertices(ofile, int(num_vertices))

    for i in range(int(num_vertices)):
        line = ifile.readline()
        edges = line.split()
        write_graph_vertex(ofile, i, len(edges))

    write_footer_vertices(ofile)
    write_header_edges(ofile, int(num_edges))

    # Step to beginning of file and skip header info
    ifile.seek(0)
    ifile.readline()
    for i in range(int(num_vertices)):
        line = ifile.readline()
        edges = line.split()
        for e in edges:
            if i < int(e):
                write_graph_edge(ofile, i, int(e))

    write_footer_edges(ofile)
    write_footer_graph(ofile)

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

def scotch_graph2graph_xml(ifilename, ofilename):
    "Convert from Scotch graph format to DOLFIN Graph XML."

    print "Converting from Scotch graph format to DOLFIN Graph XML."
    
    # Open files
    ifile = open(ifilename, "r")
    ofile = open(ofilename, "w")

    # Skip graph file version number
    ifile.readline()

    # Read number of vertices and edges
    line = ifile.readline()
    if not line:
       _error("Empty file")

    (num_vertices, num_edges) = line.split()

    # Read start index and numeric flag
    # Start index is 0 or 1 (C/Fortran)
    # Numeric flag is 3 bits where bit 1 enables vertex labels
    # bit 2 enables edge weights and bit 3 enables vertex weights

    line = ifile.readline()
    (start_index, numeric_flag) = line.split()

    # Handling not implented
    if not numeric_flag == "000":
       _error("Handling of scotch vertex labels, edge- and vertex weights not implemented")

    write_header_graph(ofile, "undirected")
    write_header_vertices(ofile, int(num_vertices))

    # Read vertices and edges, first number gives number of edges from this vertex (not used)
    for i in range(int(num_vertices)):
        line = ifile.readline()
        edges = line.split()
        write_graph_vertex(ofile, i, len(edges)-1)

    write_footer_vertices(ofile)
    write_header_edges(ofile, int(num_edges)/2)


    # Step to beginning of file and skip header info
    ifile.seek(0)
    ifile.readline()
    ifile.readline()
    ifile.readline()
    for i in range(int(num_vertices)):
        line = ifile.readline()

        edges = line.split()
        for j in range(1, len(edges)):
            if i < int(edges[j]):
                write_graph_edge(ofile, i, int(edges[j]))

    write_footer_edges(ofile)
    write_footer_graph(ofile)

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

def write_header_meshfunction(ofile, dimensions, size):
    header = """<?xml version="1.0" encoding="UTF-8"?>
<dolfin xmlns:dolfin="http://www.fenics.org/dolfin/">
  <meshfunction type="uint" dim="%d" size="%d">
""" % (dimensions, size)
    ofile.write(header)

def write_entity_meshfunction(ofile, index, value):
    ofile.write("""    <entity index=\"%d\" value=\"%d\"/>
""" % (index, value))

def write_footer_meshfunction(ofile):
    ofile.write("""  </meshfunction>
</dolfin>""")

def diffpack2xml(ifilename, ofilename):
    "Convert from Diffpack tetrahedral grid format to DOLFIN XML."

    print diffpack2xml.__doc__

    # Format strings for MeshFunction XML files
    meshfunction_header = """\
<?xml version="1.0" encoding="UTF-8"?>\n
<dolfin xmlns:dolfin="http://www.fenics.org/dolfin/">
  <meshfunction type="uint" dim="%d" size="%d">\n"""
    meshfunction_entity = "    <entity index=\"%d\" value=\"%d\"/>\n"
    meshfunction_footer = "  </meshfunction>\n</dolfin>"

    # Open files
    ifile = open(ifilename, "r")
    ofile = open(ofilename, "w")
    ofile_mat = open(ofilename.split(".")[0]+"_mat.xml", "w")
    ofile_bi = open(ofilename.split(".")[0]+"_bi.xml", "w")
   
    # Read and analyze header
    while 1:
        line = ifile.readline()
        if not line:
           _error("Empty file")
        if line[0] == "#":
            break
        if re.search(r"Number of elements", line):
	    num_cells = int(re.match(r".*\s(\d+).*", line).group(1))
        if re.search(r"Number of nodes", line):
	    num_vertices = int(re.match(r".*\s(\d+).*", line).group(1))

    write_header_mesh(ofile, "tetrahedron", 3)
    write_header_vertices(ofile, num_vertices)
    ofile_bi.write(meshfunction_header % (0, num_vertices))    
    ofile_mat.write(meshfunction_header % (3, num_cells))    

    # Read & write vertices
    # Note that only first boundary indicator is rewriten into XML
    for i in range(num_vertices):
        line = ifile.readline()
        m = re.match(r"^.*\(\s*(.*)\s*\).*\](.*)$", line)
        x = re.split("[\s,]+", m.group(1))
        write_vertex(ofile, i, x[0], x[1], x[2])
        tmp = m.group(2).split()
        if len(tmp) > 0:
            bi = int(tmp[0])     
        else:
            bi = 0
        ofile_bi.write(meshfunction_entity % (i, bi))

    write_footer_vertices(ofile)
    write_header_cells(ofile, num_cells)

    # Ignore comment lines
    while 1:
        line = ifile.readline()
        if not line:
           _error("Empty file")
        if line[0] == "#":
            break

    # Read & write cells
    for i in range(int(num_cells)):
        line = ifile.readline()
        v = line.split();
        if v[1] != "ElmT4n3D":
           _error("Only tetrahedral elements (ElmT4n3D) are implemented.") 
        write_cell_tetrahedron(ofile, i, int(v[3])-1, int(v[4])-1, int(v[5])-1, int(v[6])-1)
        ofile_mat.write(meshfunction_entity % (i, int(v[2])))
        
    write_footer_cells(ofile)
    write_footer_mesh(ofile)
    ofile_bi.write(meshfunction_footer)    
    ofile_mat.write(meshfunction_footer)    

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

class ParseError(Exception):
    """ Error encountered in source file.
    """

class DataHandler(object):
    """ Baseclass for handlers of mesh data.

    The actual handling of mesh data encountered in the source file is
    delegated to a polymorfic object. Typically, the delegate will write the
    data to XML.
    @ivar _state: the state which the handler is in, one of State_*.
    @ivar _cell_type: cell type in mesh. One of CellType_*.
    @ivar _dim: mesh dimensions.
    """
    State_Invalid, State_Init, State_Vertices, State_Cells, State_MeshFunction = range(5)
    CellType_Tetrahedron, CellType_Triangle = range(2)

    def __init__(self):
        self._state = self.State_Invalid

    def set_mesh_type(self, cell_type, dim):
        assert self._state == self.State_Invalid
        self._state = self.State_Init
        if cell_type == "tetrahedron":
            self._cell_type = self.CellType_Tetrahedron
        elif cell_type == "triangle":
            self._cell_type = self.CellType_Triangle
        self._dim = dim

    def start_vertices(self, num_vertices):
        assert self._state == self.State_Init
        self._state = self.State_Vertices

    def add_vertex(self, vertex, coords):
        assert self._state == self.State_Vertices

    def end_vertices(self):
        assert self._state == self.State_Vertices
        self._state = self.State_Init

    def start_cells(self, num_cells):
        assert self._state == self.State_Init
        self._state = self.State_Cells

    def add_cell(self, cell, nodes):
        assert self._state == self.State_Cells

    def end_cells(self):
        assert self._state == self.State_Cells
        self._state = self.State_Init

    def start_meshfunction(self, name, dim, size):
        assert self._state == self.State_Init
        self._state = self.State_MeshFunction

    def add_entity_meshfunction(self, index, value):
        assert self._state == self.State_MeshFunction

    def end_meshfunction(self):
        assert self._state == self.State_MeshFunction
        self._state = self.State_Init

    def warn(self, msg):
        """ Issue warning during parse.
        """
        warnings.warn(msg)

    def error(self, msg):
        """ Raise error during parse.

        This method is expected to raise ParseError.
        """
        raise ParseError(msg)

    def close(self):
        self._state = self.State_Invalid

class XmlHandler(DataHandler):
    """ Data handler class which writes to Dolfin XML.
    """
    def __init__(self, ofilename):
        DataHandler.__init__(self)
        self._ofilename = ofilename
        self.__ofile = file(ofilename, "wb")
        self.__ofile_meshfunc = None

    def set_mesh_type(self, cell_type, dim):
        DataHandler.set_mesh_type(self, cell_type, dim)
        write_header_mesh(self.__ofile, cell_type, dim)

    def start_vertices(self, num_vertices):
        DataHandler.start_vertices(self, num_vertices)
        write_header_vertices(self.__ofile, num_vertices)

    def add_vertex(self, vertex, coords):
        DataHandler.add_vertex(self, vertex, coords)
        write_vertex(self.__ofile, vertex, *coords)

    def end_vertices(self):
        DataHandler.end_vertices(self)
        write_footer_vertices(self.__ofile)

    def start_cells(self, num_cells):
        DataHandler.start_cells(self, num_cells)
        write_header_cells(self.__ofile, num_cells)

    def add_cell(self, cell, nodes):
        DataHandler.add_cell(self, cell, nodes)
        if self._cell_type == self.CellType_Tetrahedron:
            func = write_cell_tetrahedron
        func(self.__ofile, cell, *nodes)

    def end_cells(self):
        DataHandler.end_cells(self)
        write_footer_cells(self.__ofile)

    def start_meshfunction(self, name, dim, size):
        DataHandler.start_meshfunction(self, name, dim, size)
        fname = os.path.splitext(self.__ofile.name)[0]
        self.__ofile_meshfunc = file("%s_%s.xml" % (fname, name), "wb")
        write_header_meshfunction(self.__ofile_meshfunc, dim, size)

    def add_entity_meshfunction(self, index, value):
        DataHandler.add_entity_meshfunction(self, index, value)
        write_entity_meshfunction(self.__ofile_meshfunc, index, value)

    def end_meshfunction(self):
        DataHandler.end_meshfunction(self)
        write_footer_meshfunction(self.__ofile_meshfunc)
        self.__ofile_meshfunc.close()
        self.__ofile_meshfunc = None

    def close(self):
        DataHandler.close(self)
        if self.__ofile.closed:
            return
        write_footer_mesh(self.__ofile)
        self.__ofile.close()
        if self.__ofile_meshfunc is not None:
            self.__ofile_meshfunc.close()

def _abaqus(ifilename, handler):
    """ Convert from Abaqus.

    The Abaqus format first defines a node block, then there should be a number
    of elements containing these nodes.
    """
    params = False
    ifile = file(ifilename, "rb")
    handler.set_mesh_type("tetrahedron", 3)

    # Read definitions

    def read_params(params_spec, pnames, lineno):
        params = {}
        for p in params_spec:
            m = re.match(r"(.+)=(.+)", p)
            if m is not None:
                pname, val = m.groups()
            else:
                handler.warn("Invalid parameter syntax on line %d: %s" % (lineno, p))
                continue
            for pn in pnames:
                if pn == pname:
                    params[pn] = val
                    break

        return params
        
    nodes = {}
    elems = {}
    eid2elset = {}
    material2elsetids = {}
    materials = []
    re_sect = re.compile(r"\*([^,]+)(?:,(.*))?")
    re_node = re.compile(r"(\d+),\s*(.+),\s*(.+),\s*(.+)")
    re_tetra = re.compile(r"(\d+),\s*(\d+),\s*(\d+),\s*(\d+),\s*(\d+)")
    sect = None
    for lineno, l in enumerate(ifile):
        l = l.strip().lower()
        m = re_sect.match(l)
        if m is not None:
            sect, params_str = m.groups()
            params_spec = ([s.strip() for s in params_str.split(",")] if params_str
                    else [])

            if sect == "element":
                pnames = ("type", "elset")
                params = read_params(params_spec, pnames, lineno)
                if "type" not in params:
                   handler.error("Element on line %d doesn't declare TYPE" %
                            (lineno,))
                tp, elset = params["type"], params.get("elset")
                if tp not in ("c3d4", "dc3d4"):
                    handler.warn("Unsupported element type '%s' on line %d" % (tp, lineno))
                    supported_elem = False
                else:
                    supported_elem = True
            elif sect == "solid section":
                pnames = ("material", "elset")
                params = read_params(params_spec, pnames, lineno)
                for pname in pnames:
                    if pname not in params:
                       handler.error("Solid section on line %d doesn't "
                                "declare %s" % (lineno, pname.upper()))
                matname = params["material"]
                material2elsetids.setdefault(matname, []).append(params["elset"])
            elif sect == "material":
                name = read_params(params_spec, ["name"], lineno)["name"]
                materials.append(name)
            # We've read the section's heading, continue to next line
            continue
        
        # Read section entry