dolfin team mailing list archive
-
dolfin team
-
Mailing list archive
-
Message #23928
Patch adding support for gmsh physical markers on facets
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: mail@xxxxxxxxxxxxx-20110624060943-tennso1l72yixuo2
# target_branch: lp:dolfin
# testament_sha1: 6655ebd42c9a5f775e41a997e1ff9fae9fc69ca8
# timestamp: 2011-06-24 15:30:00 +0900
# source_branch: .
# base_revision_id: johannr@xxxxxxxxx-20110623175610-6o0pvc93qejak5fx
#
# 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=
Follow ups
