gmshtoucd.py

Go to the documentation of this file.
00001 #!/usr/bin/env python
00002 
00003 """
00004 gmshtoucd - converts Gmsh mesh files into unstructured cell data format
00005 
00006 Copyright (C) 2002-2010 Werner Scholz, Richard Boardman, Hans Fangohr
00007 
00008 gmshtoucd is free software; you can redistribute it and/or modify it
00009 under the terms of the GNU General Public License as published by the
00010 Free Software Foundation; either version 2
00011 
00012 gmshtoucd is distributed in the hope that it will be useful, but WITHOUT
00013 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00014 FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
00015 for more details.
00016 
00017 You should have received a copy of the GNU General Public License
00018 along with magpar; if not, write to the Free Software Foundation,
00019 Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
00020 
00021 This program was designed for magpar, so results may
00022 vary when used with other programs accepting AVS/UCD
00023 as an input file format
00024 
00025 This program is based on ngtoucd by Richard Boardman, Hans Fangohr.
00026 
00027 """
00028 
00029 VERSION="0.1"
00030 
00031 import os, string, sys, time
00032 
00033 def readmesh(filename):
00034     """This will read the Gmsh mesh 'filename' and return
00035     points, tetrahedra and triangles"""
00036 
00037     # expected Gmsh msh file format:
00038     # http://www.geuz.org/gmsh/doc/texinfo/gmsh_10.html#SEC62
00039     
00040     f = open(filename, 'r')
00041 
00042     #################################################################
00043     #
00044     # read the first line to get the number of points
00045     #
00046     #################################################################
00047 
00048     thisline = f.readline()
00049     if thisline != "$MeshFormat\n" :
00050         bail("Unexpected input " + thisline)
00051 
00052     thisline = f.readline()
00053     segments  = string.split(thisline[:-1])
00054     # `.msh' file format, version 2.0
00055     if segments[0] != "2" and segments[0] != "2.1":
00056       bail("Unexpected input " + thisline)
00057     # file-type: ASCII format
00058     if int(segments[1])!=0 :
00059       bail("Unexpected input " + thisline)
00060     # data-size (ignored)
00061 
00062     thisline = f.readline()
00063     if thisline != "$EndMeshFormat\n" :
00064         bail("Unexpected input " + thisline)
00065 
00066     #################################################################
00067     #
00068     # read the points
00069     #
00070     #################################################################
00071 
00072     thisline = f.readline()
00073     if thisline != "$Nodes\n" :
00074         bail("Unexpected input " + thisline)
00075       
00076     thisline = f.readline()
00077     number_of_points = int(thisline)
00078 
00079     points = []
00080 
00081     for i in range(number_of_points):
00082         thisline  = f.readline()
00083         segments  = string.split(thisline[:-1])
00084         thispoint = list(map(lambda p: float(p), segments[1:4]))
00085         points.append(thispoint)
00086 
00087     thisline = f.readline()
00088     if thisline != "$EndNodes\n" :
00089         bail("Unexpected input " + thisline)
00090       
00091     #################################################################
00092     #
00093     # now read all the volume elements
00094     # first, ascertain the number of volume elements
00095     #
00096     #################################################################
00097 
00098     thisline = f.readline()
00099     if thisline != "$Elements\n" :
00100         bail("Unexpected input " + thisline)
00101       
00102     thisline             = f.readline()
00103     number_of_elements   = int(thisline)
00104 
00105     #################################################################
00106     #
00107     # then read and populate, keeping track
00108     # of the tetrahedra subdomain
00109     #
00110     #################################################################
00111     
00112     triangles           = []
00113     triangles_subdomain = []
00114 
00115     tetrahedra           = []
00116     tetrahedra_subdomain = []
00117 
00118     triangles_id=0
00119     tetrahedra_id=0
00120 
00121     for i in range(number_of_elements):
00122         thisline = f.readline()
00123         segments = string.split(thisline[:-1])
00124         tags = int(segments[2])
00125 
00126         # element type 2: 3 node triangle
00127         if int(segments[1]) == 2 :
00128             segments = string.split(thisline[:-1])
00129             if len(segments) != 6+tags:
00130                 print "Found", len(segments), "elements, but expected ", 6+tags
00131                 print segments
00132                 raise "Oops", "File format error"
00133 
00134             triangles_subdomain.append(triangles_id)
00135             # similar adjustment for materials as with tetrahedra
00136             triangle = [list(map(lambda t : int(t)-1, segments[3+tags:5+tags])), int(segments[3])]
00137             triangles.append(triangle)
00138             triangles_id=triangles_id+1
00139 
00140         # element type 4: 4 node tetrahedron
00141         if int(segments[1]) == 4 :
00142             if len(segments) != 7+tags:
00143                 print "Found", len(segments), "elements, but expected ", 7+tags
00144                 print segments
00145                 raise "Oops", "File format error"
00146 
00147             tetrahedra_subdomain.append(tetrahedra_id)
00148 
00149             # rpb, 17th March 2006 - updated to handle materials (note:
00150             # in Gmsh these are top-level objects [tlos])
00151             # original mapping was [a, b, c, d] (where a-d are points)
00152             # new mapping is [[a,b,c,d],m] (where m is a material ID)
00153             tetrahedron = [list(map(lambda t : int(t)-1, segments[3+tags:7+tags])), int(segments[3])]
00154 
00155             # Gmsh counts from 1 rather than from 0, hence int(t)-1
00156 
00157             tetrahedra.append(tetrahedron)
00158             tetrahedra_id=tetrahedra_id+1
00159 
00160     thisline = f.readline()
00161     if thisline != "$EndElements\n" :
00162         bail("Unexpected input " + thisline)
00163       
00164     #################################################################
00165     #
00166     # send back the read mesh
00167     #
00168     #################################################################
00169     
00170     return (points, tetrahedra, triangles, tetrahedra_subdomain, triangles_subdomain)
00171 
00172 def bail(message="Sorry, I don't understand. Bailing out..."):
00173     """A 'get-out' clause"""
00174     print message
00175     sys.exit(1)
00176 
00177 def initialise_file(file):
00178     """Attempt to remove the file first to avoid append issues"""
00179     try:
00180         os.remove(file)
00181     except:
00182         # no problem
00183         pass
00184 
00185 def writeln(line, file):
00186     """One-shot file line append; this keeps the code terse"""
00187     f = open(file, 'a')
00188     f.write(line + '\n')
00189     f.close()
00190 
00191 def stamp():
00192     """Format a time and a date for stdout feedback"""
00193     thistime = time.localtime()
00194     return "[%04d/%02d/%02d %02d:%02d:%02d] " % (thistime[0],
00195                                                  thistime[1],
00196                                                  thistime[2],
00197                                                  thistime[3],
00198                                                  thistime[4],
00199                                                  thistime[5])
00200 
00201 
00202 def info():
00203     """Output GPL details"""
00204     print """gmshtoucd, Copyright (C) 2002-2010
00205              Werner Scholz, Richard Boardman, Hans Fangohr
00206              gmshtoucd comes with ABSOLUTELY NO WARRANTY;
00207              for details see the file COPYING"""
00208 
00209 def writeucd():
00210     """Call readmesh to read the Gmsh mesh, then convert and write the UCD mesh"""
00211     al = len(sys.argv)
00212     if al < 2 or al > 3 or sys.argv[1] == '--help' or sys.argv[1] == '-h':
00213         bail("Usage: " + sys.argv[0] + " gmsh.msh [ucdmesh.inp]")
00214 
00215     inputfile = sys.argv[1]
00216 
00217     if len(sys.argv) == 3:
00218         outputfile = sys.argv[2]
00219     else:
00220         outputfile = os.path.splitext(inputfile)[0] + ".inp"
00221 
00222     print
00223     print "About to read, check and process the gmsh file ", inputfile
00224     print "and then create the AVS/UCD mesh", outputfile
00225     print
00226 
00227     print stamp() + "Initialising file" + outputfile
00228     initialise_file(outputfile)
00229     print stamp() + "File " + outputfile + " initialised"
00230 
00231     print stamp() + "Attemping read of input " + inputfile
00232     (points, tetrahedra, triangles, tetrasub, trisub) = readmesh(inputfile)
00233     print stamp() + "Input file successfully read"
00234 
00235     print stamp() + "Points:     %d" % len(points)
00236     print stamp() + "Tetrahedra: %d [%d in subdomain]" % (len(tetrahedra), len(tetrasub))
00237     print stamp() + "Triangles:  %d [%d in subdomain]" % (len(triangles),  len(trisub))
00238 
00239     #################################################################
00240     #
00241     # the first line of a UCD file reads:
00242     # a b c d e
00243     #
00244     # where a is the number of nodes
00245     #       b is the number of cells
00246     #       c is the length of vector data associated with the nodes
00247     #       d is the length of vector data associated with the cells
00248     #       e is the length of vector data associated with the model
00249     #
00250     # example: 12 2 1 0 0
00251     #
00252     #################################################################
00253 
00254     # n.b. here the third integer indicates the number of placeholders
00255     
00256     print stamp() + "Creating descriptor [UCD]"
00257     writeln(str(len(points)) + " " + str(len(tetrahedra)) + " 3 0 0", outputfile)
00258     print stamp() + "UCD descriptor created"
00259     
00260     #################################################################
00261     #
00262     # then we have nodes in threespace, one line per node
00263     # n x y z
00264     # where n is the node ID -- integer (not necessarily sequential)
00265     #       x,y,z are the x, y and z coordinates
00266     #
00267     #################################################################
00268 
00269     print stamp() + "Now converting nodes"
00270     for i in range(len(points)):
00271         x, y, z = points[i][0], points[i][1], points[i][2]
00272         writeln(str(i+1) + " " + str(x) + " " + str(y) + " " + str(z), outputfile)
00273     print stamp() + "Nodes converted"
00274 
00275     #################################################################
00276     #
00277     # now the cells, one line/cell:
00278     #
00279     # c m t n1 n2 ... nn
00280     #
00281     # where c is the cell ID
00282     #       m is the material type (int, leave as 1 if we don't care)
00283     #       t is the cell type (prism|hex|pyr|tet|quad|tri|line|pt)
00284     #       n1...nn is a node ID list corresponding to cell vertices
00285     #
00286     #################################################################
00287 
00288     cellctr = 0
00289 
00290     print stamp() + "Converting tetrahedra"
00291 
00292     for tetra in tetrahedra:
00293         tet = tetra[0] # just the tetrahedron; tetra[1] is the material
00294         tetorder  = [0, 1, 2, 3]
00295         tetstring =  " " + str(tet[tetorder[0]]+1)
00296         tetstring += " " + str(tet[tetorder[1]]+1)
00297         tetstring += " " + str(tet[tetorder[2]]+1)
00298         tetstring += " " + str(tet[tetorder[3]]+1)
00299         writeln(str(cellctr+1) + " " + str(tetra[1]) + " tet" + tetstring, outputfile)
00300         cellctr   += 1
00301 
00302     print stamp() + "Tetrahedra converted"
00303 
00304     #################################################################
00305     #
00306     # uncomment the following section to convert triangles as well
00307     # - since magpar will ignore these, leave away for the time being
00308     #
00309     # print stamp() + "Converting triangles"
00310     #
00311     # cellctr = 0
00312     # for triang in triangles:
00313     #     tri = triang[0] # just the triangle; triang[1] is the material
00314     #     tristring = ""
00315     #     for node in tri:
00316     #         tristring += " " + str(node+1)
00317     #     writeln(str(cellctr+1) + " " + str(triang[1]) + " tri" + ts, outputfile)
00318     #
00319     # print stamp() + "Triangles converted"
00320     #
00321     #################################################################
00322 
00323     #################################################################
00324     #
00325     # for the data vector associated with the nodes:
00326     #
00327     # first line tells us into which components the vector is divided
00328     #
00329     #   example: vector of 5 floats could be 3 3 1 1
00330     #            node scalar could be 1 1
00331     #
00332     # next lines, for each data component, use a cs label/unit pair
00333     #
00334     #   example: temperature, kelvin
00335     #
00336     # subsequent lines, for each node, the vector of associated data
00337     # in this order
00338     #
00339     #   example: 1 10\n2 15\n3 12.4\n4 9
00340     #
00341     #################################################################
00342 
00343     print stamp() + "Initialising placeholders"
00344 
00345     writeln("3 1 1 1", outputfile)
00346     writeln("M_x, none", outputfile)
00347     writeln("M_y, none", outputfile)
00348     writeln("M_z, none", outputfile)
00349     
00350     #################################################################
00351     #
00352     # create some initial scalar values; here, generally
00353     # +x=1.0 +y=+z=0.01*x (everything else)="\epsilon"
00354     #
00355     #################################################################
00356 
00357     scalarstring =  " 1.0 0.0 0.0"
00358 
00359     for i in range(len(points)):
00360         writeln(str(i+1) + scalarstring, outputfile)
00361 
00362     print stamp() + "Placeholders inserted"
00363 
00364     # all done
00365 
00366     print stamp() + "All finished"
00367 
00368 
00369     
00370 if __name__=="__main__":
00371     writeucd()

magpar - Parallel Finite Element Micromagnetics Package
Copyright (C) 2002-2010 Werner Scholz