Hello everyone,
I’m now trying to convert my .msh file to .xml files, all works fine except for the facet_region.xml. It seems that it transforms only the lines defined as physical lines in .geo file and eliminate all edges of mesh.
this is my inital .geo file:
Blockquote
d = 0.015;
R = 0.25;
Point(1) = {0., 0., 0., d};
Point(2) = {0.4-R, 0., 0., d};
Point(3) = {0.4, 0., 0., d};
Point(4) = {0.4+R, 0., 0., d};
Point(5) = {1., 0., 0., d};
Point(6) = {1., 0.5, 0., d};
Point(7) = {0.6+R, 0.5, 0., d};
Point(8) = {0.6, 0.5, 0., d};
Point(9) = {0.6-R, 0.5, 0., d};
Point(10) = {0., 0.5, 0., d};
Line(1) = {1, 2};
Line(2) = {2, 3};
Line(3) = {3, 4};
Line(4) = {4, 5};
Line(5) = {5, 6};
Line(6) = {6, 7};
Line(7) = {7, 8};
Line(8) = {8, 9};
Line(9) = {9, 10};
Line(10) = {10, 1};
Circle(11) = {4,3,2};
Circle(12) = {9,8,7};
//+
Line Loop(1) = {4, 5, 6, -12, 9, 10, 1, -11};
//+
Plane Surface(1) = {1};
//+
Line Loop(2) = {11, 2, 3};
//+
Plane Surface(2) = {2};
//+
Line Loop(3) = {12, 7, 8};
//+
Plane Surface(3) = {3};
Physical Line(1) = {10};
Physical Line(2) = {5};
Physical Line(3) = {11, 12};
//+
Physical Surface(1) = {1};
Physical Surface(2) = {2,3};
my purpose is to give labels for the facets of the mesh, either exterior ones (labelled 1 and 2 for the left and right boundary respectively) or labelled 3 for the interface between the two materials, and all the others are labelled by default value 0.
but by runding this code:
def collect_segments(mesh):
segment_data = {}
segment_types = {"line", "line3"} # Assuming 'line' and 'line3' are your segment types
if 'gmsh:physical' in mesh.cell_data:
# Iterate over each cell block and the associated physical tags
for cell_block, physical_tags in zip(mesh.cells, mesh.cell_data['gmsh:physical']):
if cell_block.type in segment_types:
segment_data[cell_block.type] = dict()
for tag, cell in zip(physical_tags, cell_block.data):
segment_data[cell_block.type][tuple(sorted(cell))] = tag
return segment_data
i’ve got 174 facets (segments) in total, while the real number should be 9156. (it ignores all lines with label 0 )
To solve this problem, i’ve tried a new approach:
def collect_segments(mesh):
segments = set()
for cell_block in mesh.cells:
for cell in cell_block.data:
# Check all two-node combinations within each cell
if len(cell) >= 2:
for i in range(len(cell)):
for j in range(i + 1, len(cell)):
# Create a sorted tuple to avoid duplicate edges being counted
edge = tuple(sorted([cell[i], cell[j]]))
# print(edge)
segments.add(edge)
return segments
Fortunately, it gives out all 9156 elements, but since this approach is not Mesh Dependency, it seems to have a problem with the ordering of the tuple of the segments. and the resulting facet_region.xml file can not be correctly read by fenics.
Is there any efficient approach to solve this problem?
