Hi Please find a small patch with modifications to meshconvert.py that allows for the handling of facet physical regions using gmsh. This is done in much the same way as the existing physical region (for cells) support.
I am working on a few test cases, but at least the existing meshconvert test cases do not fail :) Please let me know if there are any questions or comments. Evan -- www.evanlezar.com GoogleTalk: evanlezar Skype: evanlezar
# Bazaar merge directive format 2 (Bazaar 0.90) # revision_id: [email protected] # target_branch: lp:dolfin # testament_sha1: 6655ebd42c9a5f775e41a997e1ff9fae9fc69ca8 # timestamp: 2011-06-24 15:30:00 +0900 # source_branch: . # base_revision_id: [email protected] # # Begin patch === modified file 'site-packages/dolfin_utils/meshconvert.py' --- site-packages/dolfin_utils/meshconvert.py 2011-06-02 19:26:59 +0000 +++ site-packages/dolfin_utils/meshconvert.py 2011-06-24 06:09:43 +0000 @@ -28,6 +28,7 @@ # Modified by Kent-Andre Mardal (Star-CD function) # Modified by Nuno Lopes (fix for emc2 mesh format (medit version 0)) # Modified by Neilen Marais (add gmsh support for reading physical region) +# Modified by Evan Lezar (add support for reading gmsh physical regions on facets) # NOTE: This module does not depend on (py)dolfin beeing installed. # NOTE: If future additions need that please import dolfin in a try: except: @@ -243,13 +244,19 @@ """ print "Converting from Gmsh format (.msh, .gmsh) to DOLFIN XML format" - + + # The dimension of the gmsh element types supported here as well as the dolfin cell types for each dimension + gmsh_dim = {1: 1, 2: 2, 4: 3} + cell_type_for_dim = {1: "interval", 2: "triangle", 3: "tetrahedron" } + # the gmsh element types supported for conversion + supported_gmsh_element_types = [1, 2, 4] + # Open files ifile = open(ifilename, "r") # Scan file for cell type cell_type = None - dim = 0 + highest_dim = 0 line = ifile.readline() while line: @@ -261,63 +268,48 @@ if line.find("$Elements") == 0: line = ifile.readline() - num_cells = int(line) - num_cells_counted = 0 - if num_cells == 0: - _error("No cells found in gmsh file.") + num_elements = int(line) + if num_elements == 0: + _error("No elements found in gmsh file.") line = ifile.readline() # Now iterate through elements to find largest dimension. Gmsh # format might include elements of lower dimensions in the element list. # We also need to count number of elements of correct dimensions. # Also determine which vertices are not used. - dim_2_count = 0 - dim_3_count = 0 - vertices_2_used = [] - # Array used to store gmsh tags for 2D (type 2/triangular) elements - tags_2 = [] - # Array used to store gmsh tags for 3D (type 4/tet) elements - tags_3 = [] - vertices_3_used = [] + dim_count = {1: 0, 2: 0, 3: 0} + vertices_used_for_dim = {1: [], 2: [], 3: []} + # Array used to store gmsh tags for 1D (type 1/line), 2D (type 2/triangular) elements and 3D (type 4/tet) elements + tags_for_dim = {1: [], 2: [], 3: []} + while line.find("$EndElements") == -1: element = line.split() elem_type = int(element[1]) num_tags = int(element[2]) - if elem_type == 2: - if dim < 2: - cell_type = "triangle" - dim = 2 - node_num_list = [int(node) for node in element[3 + num_tags:]] - vertices_2_used.extend(node_num_list) - if num_tags > 0: - tags_2.append(tuple(int(tag) for tag in element[3:3+num_tags])) - dim_2_count += 1 - elif elem_type == 4: - if dim < 3: - cell_type = "tetrahedron" - dim = 3 - vertices_2_used = None - node_num_list = [int(node) for node in element[3 + num_tags:]] - vertices_3_used.extend(node_num_list) - if num_tags > 0: - tags_3.append(tuple(int(tag) for tag in element[3:3+num_tags])) - dim_3_count += 1 + if elem_type in supported_gmsh_element_types: + dim = gmsh_dim[elem_type] + if highest_dim < dim: + highest_dim = dim + node_num_list = [int(node) for node in element[3 + num_tags:]] + vertices_used_for_dim[dim].extend(node_num_list) + if num_tags > 0: + tags_for_dim[dim].append(tuple(int(tag) for tag in element[3:3+num_tags])) + dim_count[dim] += 1 + else: + #TODO: output a warning here. "gmsh element type %d not supported" % elem_type + pass line = ifile.readline() else: # Read next line line = ifile.readline() # Check that we got the cell type and set num_cells_counted - if cell_type == None: - _error("Unable to find cell type.") - if dim == 3: - num_cells_counted = dim_3_count - vertex_set = set(vertices_3_used) - vertices_3_used = None - elif dim == 2: - num_cells_counted = dim_2_count - vertex_set = set(vertices_2_used) - vertices_2_used = None + if highest_dim == 0: + _error("Unable to find cells of supported type.") + + num_cells_counted = dim_count[highest_dim] + vertex_set = set(vertices_used_for_dim[highest_dim]) + vertices_used_for_dim[highest_dim] = None vertex_dict = {} for n,v in enumerate(vertex_set): @@ -327,7 +319,7 @@ ifile.seek(0) # Set mesh type - handler.set_mesh_type(cell_type, dim) + handler.set_mesh_type(cell_type_for_dim[highest_dim], highest_dim) # Initialise node list (gmsh does not export all vertexes in order) nodelist = {} @@ -339,8 +331,21 @@ num_vertices_read = 0 num_cells_read = 0 + # Only import the dolfin objects if facet markings exist + process_facets = False + if len(tags_for_dim[highest_dim-1]) > 0: + # first construct the mesh + from dolfin import MeshEditor, Mesh + mesh = Mesh() + mesh_editor = MeshEditor () + mesh_editor.open( mesh, highest_dim, highest_dim ) + process_facets = True + else: + # TODO: Output a warning or an error here + me = None + while state != 10: - + # Read next line line = ifile.readline() if not line: break @@ -368,6 +373,8 @@ elif state == 4: num_vertices = len(vertex_dict) handler.start_vertices(num_vertices) + if process_facets: + mesh_editor.init_vertices ( num_vertices ) state = 5 elif state == 5: (node_no, x, y, z) = line.split() @@ -379,6 +386,14 @@ continue nodelist[int(node_no)] = num_vertices_read handler.add_vertex(num_vertices_read, [x, y, z]) + if process_facets: + if highest_dim == 1: + mesh_editor.add_vertex( num_vertices_read, x) + elif highest_dim == 2: + mesh_editor.add_vertex( num_vertices_read, x, y) + elif highest_dim == 3: + mesh_editor.add_vertex( num_vertices_read, x, y, z) + num_vertices_read +=1 if num_vertices == num_vertices_read: @@ -392,54 +407,92 @@ state = 8 elif state == 8: handler.start_cells(num_cells_counted) + if process_facets: + mesh_editor.init_cells( num_cells_counted ) + state = 9 elif state == 9: element = line.split() elem_type = int(element[1]) num_tags = int(element[2]) - if elem_type == 2 and dim == 2: - node_num_list = [vertex_dict[int(node)] for node in element[3 + num_tags:]] - for node in node_num_list: - if not node in nodelist: - _error("Vertex %d of triangle %d not previously defined." % - (node, num_cells_read)) - cell_nodes = [nodelist[n] for n in node_num_list] - handler.add_cell(num_cells_read, cell_nodes) - num_cells_read +=1 - elif elem_type == 4 and dim == 3: - node_num_list = [vertex_dict[int(node)] for node in element[3 + num_tags:]] - for node in node_num_list: - if not node in nodelist: - _error("Vertex %d of tetrahedron %d not previously defined." % - (node, num_cells_read)) - # import pdb ; pdb.set_trace() - cell_nodes = [nodelist[n] for n in node_num_list] - handler.add_cell(num_cells_read, cell_nodes) + if elem_type in supported_gmsh_element_types: + dim = gmsh_dim[elem_type] + else: + dim = 0 + if dim == highest_dim: + node_num_list = [vertex_dict[int(node)] for node in element[3 + num_tags:]] + for node in node_num_list: + if not node in nodelist: + _error("Vertex %d of %s %d not previously defined." % + (node, cell_type_for_dim[dim], num_cells_read)) + cell_nodes = [nodelist[n] for n in node_num_list] + handler.add_cell(num_cells_read, cell_nodes) + + if process_facets: + mesh_editor.add_cell( num_cells_read, *cell_nodes ) + num_cells_read +=1 if num_cells_counted == num_cells_read: handler.end_cells() + if process_facets: + mesh_editor.close() state = 10 elif state == 10: break - + # Write mesh function based on the Physical Regions defined by # gmsh, but only if they are not all zero. All zero physical # regions indicate that no physical regions were defined. - if dim == 2: - tags = tags_2 - elif dim == 3: - tags = tags_3 - else: + if highest_dim not in [1,2,3]: _error("Gmsh tags not supported for dimension %i. Probably a bug" % dim) - + + tags = tags_for_dim[highest_dim] physical_regions = tuple(tag[0] for tag in tags) - if not all(tag == 0 for tag in tags): - handler.start_meshfunction("physical_region", dim, num_cells) + if not all(tag == 0 for tag in physical_regions): + handler.start_meshfunction("physical_region", dim, num_cells_counted) for i, physical_region in enumerate(physical_regions): handler.add_entity_meshfunction(i, physical_region) handler.end_meshfunction() + # Now process the facet markers + tags = tags_for_dim[highest_dim-1] + if (len(tags) > 0) and (mesh is not None): + physical_regions = tuple(tag[0] for tag in tags) + if not all(tag == 0 for tag in physical_regions): + mesh.init(highest_dim-1,0) + + # Get the facet-node connectivity information (reshape as a row of node indices per facet) + facets_as_nodes = mesh.topology()(highest_dim-1,0)().reshape ( mesh.num_facets(), highest_dim ) + + facets_to_check = range( mesh.num_facets() ) + + data = [int(0*k) for k in range(len(facets_to_check)) ] + for i, physical_region in enumerate(physical_regions): + nodes = [n-1 for n in vertices_used_for_dim[highest_dim-1][2*i:(2*i+highest_dim)]] + nodes.sort() + + if physical_region != 0: + found = False + for j in range(len(facets_to_check)): + index = facets_to_check[j] + if all ( facets_as_nodes[index,k] == nodes[k] for k in range(len(nodes)) ): + found = True; + facets_to_check.pop(j) + # set the value of the mesh function +# facet_mark_function[index] = physical_region + data[index] = physical_region + break; + + if not found: + raise Exception ( "The facet (%d) was not found to mark: %s" % (i, nodes) ) + +# # Create and initialise the mesh function + handler.start_meshfunction("facet_region", highest_dim-1, mesh.num_facets() ) + for index, physical_region in enumerate ( data ): + handler.add_entity_meshfunction(index, physical_region) + handler.end_meshfunction() + # Check that we got all data if state == 10: print "Conversion done" # Begin bundle IyBCYXphYXIgcmV2aXNpb24gYnVuZGxlIHY0CiMKQlpoOTFBWSZTWRcd5iwABUN/gERURFB6f/// f8ePCr////pgC1M93O+wdzN9673Not6vsmfer51z7p9t3ajfXdWnQN70JJEykyeTKY0aGppqeE0j 1PSaPKZDTQGmmRpoANBA0CAQVNjKPST1NqGmgDQAAaNAeoMgQk1T2U9QYDUymhmmSDQaAAANAGgk KSnlP0U08kIGjajajIaNBoAMgA0AAIoqegptU8TTSehPSaaZpGmJkGQ0wQANGgAkSBCaABDSj2qe RPaUGI2kNqAaAaGQDAEAPqx4eUP4ukf7xSc0rU0bp6Vn+t+wePyRd0l29K+Az3SDTR1Y+v1btmu0 k0L8jSioFO9UCNrpR6ogtwxXTTCyXYZF0QILqw00xhVAaAo3N1XzLCEfaI2uFQxVdD8Wh8J3abYJ Bz7zY1nieW15IJtNlqEMsSpoIJJrRQULwXEIJnhUpCAzDO5vaL4HNl/fyqtp1tHP4FV6z22ehXLO RzNB/Y/mrMt3cduNhMwA42Xydz0y+ogplyjy3ptHhz7ukpicmMLTxIBTSsIiVE81EfSQKYrrAdpg CKoHmFjEDjeUhdJFQyY8YHa9BZ6okEDZrkGa3me3LEJnxKEo+pQ4kPvJi6oekIyARwhuMf5j4+BP qJ7pmNIIg6WSVZChGGBU7Ga5zweqPG4wqEWUvZwrWXpXew0Nu6QeVixHvoWWLCIZH4Vm420nM6bD maaH1heO0YA6Cow04riM8HbhqCMnvm1ZVmJvDlQ1FvWahXMsozu8ra1Ei2GBy6pB/tIUWsTSycDP M7isAaj5ggEnn/sksmfXfB2GeIvgVA+riD3Ez+sTGfijrjlgKzqKgYIwqIwM6SSaSVZiSpZBs+lx 7FVyFaRb2cUZtGSbRxZ9Vu0XbjydRrnYezhwOcbDam9BuZOheiQM0IwVrwYV2HSZBpXre9BiTQ/F Kkw8Ozr+lKwiDyXyslJRJnEQXThr5GulGy+cCRolWBXizeu/0qtvE0mQIlktHd7d+/OE0acRF4VD sdEXROYmheftr5d8EuKOcoYYxrCso9Tp1g/epsRjfAqrpwNgWcR3gKCyL+7GChtu9SKY20zcDiWN 8V5kDSuNYOTsVtVP32IUYcDs+vn4HaY7EOJoKWSMPA0GyzDhj0QonbW2y+rYKN6njIsz0VN2axS8 TQYcxfS7ILFFwX3r21tLpV05IIjGNVL9ODrcX3RSadVTsd8wBdcy3te5ZALMRhFoYy2MNk5zxD7w 2dNijEDvtOLklYsvgcOVZXbqIqgWPjNz038gPUEEzCRmscQNGeYvyiGtsINgrwZNdry7YqkQKdWc HwaLSgW64I7SdDN1dJgcZCvvB6yhm0N1l0MWq4qy56WzctMjmuIQ6kKG2RVPd1v2NhqNB4PLNgNs nOczERgZK7kkhggYJDOYMAirFFbCcnaWrUArCgWWMIXYotVKsOp7XXFpFxVJ04Ji5fGIdxMF9x2O KIGlQ3WBkSBCVSAw06nYOOF9cybmVtZHbTDFxzBQ4hSi+1BulYn1ikhNNdzhd+8IKpN9IhYzpYVP N/KtAdyigKamMqXgmeOKgHPs4OXfsGsoutkOxWZ9i1b9rPSm7RaTpmZWVbfywYOG0tdV0ZXEefab XaQV2TIhOroXhazx09HeNfMHq8IT0gZhmBX4zU8rPaCDP1uWGtKsPDTCzQyY7nzU8PcLew85bBo9 rXraPZrgPyP6IM79ToxnhBy4KBOpAzLjr5YEPFB72I6MLr8lqLPbmgdJaJdfJDfm5weHDTuWRUVm cxWjxpDRwUTb2IJVhkhZSQaiz9fV+3TgVTRudoiGNsjZCw2TRVEgn4ci5uqGKBUbYRlzzLgmhqbK Kami+9SoDOcKhtuDMZrSz0nVdIdkKmLm1JHp6Rk5ee4DTsCReie2sbbbadUKcxquiUxl0lNCHKMh YUW8OLIZdWVb15FAdhKSDmP6QgVK3cZ43w1aOki6fV2TdCSOHKWNdS8W/iQugnyy3z8a3/VtLQS1 WkicCkEve2eSKJ5GVHpA5966quFldUqxVG2mnoihyqkomOlglyjNra+qhgpy5SrmszADKbxufXT5 bScYJhAqAigtD2usWW0Dcb5YBRg2k2JcTTWtZYr+9k2BWC56ND72GQz5eRLhsQDlZGZ8mkHmQMwD hGAo9JOfyLFVKSTYwYksD6dna3vDotvNcsPPlry0NMSGNwMyHDQFcyykpDX0VEN6COSFDopBS5kA 2/Dh7vA55XAATwXHcxNl73F9aBhcCZhjSSnnwIDl1QvJBwsip2xwXMhsZGZW9BIQUrismBX8k5Dq LrapfFYDwMiKrWFaNOMizODKAbub1MGxMYjIYFgDxQJzFh/ITBkZGIP2lcil7TQzO1sz+mKWHazi bSbBNdxxwywcQwF8nk5xSaOmH6Ckuz14mBmX6j6WLgNpc7kZTG5wtDIoQmNsYUSLqaDGVFjcr5sM Bg32jgGjUgUtABPfgLkzBXKjSHw61wYjiyGJDaNhtweuwUR99w2gw2YIjKu9A2JLd0yP1g2ThZhN 5HYFOeqrCamtrFdN1woyaWAfhf3402mN27jK6M0H5xyySEWg5MSrMgWLGlMjLRAkzE1AybE2m0DP KGeXNwUaSp0Sw2leBWarZXS9/dv18c7kIMccGDjZHEu5jhsZ7A45XG+WMinBX2j2Dgtux15Nt/GQ OY3S0zy5GM8z2X6HjB5io2wxDOEeApBySkqJDWVEs/R3EluuUNMNLUGrVXgSsr68OkFCJJCDYkXF Zlxl0DYqojXjQQWGkMTG0IuE0xKyGGQ4j5+OnOzeHLMZDIJZCtGan2sJckciMmfO6yDPib+zVooU OhCYTQOY0nEDLiAgaBhpaY2nuNDdP9qFuyAxljjuKBJEqCtIPmwjQxCJtxBgGJ6YUCYOsDnFKyMp yZFY7nA6FnZmBMQ1ksCCmVFihhBg1IDYGpwtV25ir0S8UBYWt8L5v/ffcIcj3xrG8UTz3ayFBHlm 8LxfcHJHdQ5AaUFmiTQs676SSaJ13RmNTvglCsN5CUYGAqccAXnBiQyVpfG6/fwV7nFVuNrqmeY7 Gm0FtFokKEuT+BnHgV3YXtgYUQWN2sRm0c45FgtCzYNsLWkRgNkGESkTFvc4SQrmD7/LiFjmZW+P EKlf+rMRNNVRl27NPmsZMxfeGAM9N4FOo7LT0Bh+GQXHEIWJCUIB5PSGY9mqPjBUIUTnGRlqCswd mTBHG8FqkpEwKSCJZThTTZMgcElZKSLElRUE0Uci4whdQtunQb7lxQ0xoGmgadOl1kI60L6w16jM FdC35qCy6siC5XtECFMxDOWaOtktdb4IrfryDKShUG23EdrqpPFhIq4c1IumisirQKVSozzklSy9 0xBmU0yJQvJ1zm8VZYYGINqNwtk+/tqK86NEJGGuVEjZtNp1WhLF3rYjfgkahVWBq35Z9S12DDCJ GLRRV2s0M0oooxNM4t/bqcx5G4GTNMITAVJHOUixF5BRQiOEc7dWeYosdZUlZcS33k2pTa4dCmbQ TXw6+MsSpIgghMOVnZhLR4K7rYM1OOUcRPpHWJgR1IUd6VBKWx+EkgLN3HN/xdyRThQkBcd5iwA=
_______________________________________________ Mailing list: https://launchpad.net/~dolfin Post to : [email protected] Unsubscribe : https://launchpad.net/~dolfin More help : https://help.launchpad.net/ListHelp
