# Mesh Markers in GMSH **URL:** \ **Category:** variational formulation\ **Tags:** mesh, boundary-conditions, dolfinx, mpi, gmsh\ **Created:** [March 2, 2025, 7:58pm UTC](https://fenicsproject.discourse.group/t/mesh-markers-in-gmsh/17124 "2025-03-02T19:58:57Z")\ **Posts on this page:** 3\ **Page:** 1
**Author:** ![umarkhayaz](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/umarkhayaz/32/6239_2.png) [@umarkhayaz](https://fenicsproject.discourse.group/u/umarkhayaz)\ **Post date:** [March 2, 2025, 7:58pm UTC](https://fenicsproject.discourse.group/t/mesh-markers-in-gmsh/17124/1 "2025-03-02T19:58:57Z")
I am currently working on solving a single edge notch problem implementing phase-field formulation. I have successfully written the code to mark boundary subdomains in FEniCS. However, I now want to mark these boundary subdomains in GMSH. Part of my .gmsh script is structured as follows: ```auto // ------------------------------ // Geometry: Define Points // ------------------------------ Point(1) = {0.0, L/2, 0.0, h}; Point(2) = {W, L/2, 0, h}; Point(3) = {W, -L/2, 0, h}; Point(4) = {0.0, -L/2, 0, h}; Point(5) = {0.0, 0.00, 0, h}; Point(6) = {L/2, 0.00, 0, h}; // ------------------------------ // Geometry: Define Outer Boundary and Surface // ------------------------------ Line(7) = {1, 2}; Line(8) = {2, 3}; Line(9) = {3, 4}; Line(10) = {4, 5}; Line(11) = {5, 1}; Line Loop(12) = {7, 8, 9, 10, 11}; Plane Surface(13) = {12}; Physical Surface(14) = {13}; // ------------------------------ // Geometry: Define Internal Crack // ------------------------------ Line(15) = {5, 6}; Curve{15} In Surface{13}; Physical Curve(16) = {15}; Physical Point(17) = {5}; // ------------------------------ // Mark Boundary Subdomains (Physical Lines) // ------------------------------ Physical Point("RightTop") = {2}; // RighTop boundary point (x, y = W, L/2) Physical Line("Top") = {7}; // Top boundary (y = L/2) Physical Line("Bottom") = {9}; // Bottom boundary (y = -L/2) Physical Line("Left") = {8}; // Left boundary (x = 0) Physical Line("Right") = {10}; // Right boundary (x = W) Physical Line("Crack") = {15}; // Crack (y = 0) ``` Additionally, the FEniCS code I wrote to define the Dirichlet boundary conditions is as follows: ```auto # Create mesh and define function space mesh=Mesh("crack.xml") h=FacetArea(mesh) #area/length of a cell facet on a given mesh h_avg = (h('+') + h('-'))/2 n=FacetNormal(mesh) # Assume the mesh conversion produces a MeshFunction of type "size_t" for boundaries. boundary_markers = MeshFunction("size_t", mesh, mesh.topology().dim() - 1) # Define ds measure using the loaded markers ds = ds(subdomain_data=boundary_markers) physical_tags = { "RightTop": 2, # Example tag number for the RightTop point "Top": 7, # Example tag number for Top boundary "Bottom": 9, # Example tag number for Bottom boundary "Left": 8, # Example tag number for Left boundary "Right": 10, # Example tag number for Right boundary "Crack": 15 # Example tag number for Crack curve } # Define Dirichlet boundary conditions c = Expression(...) r = Expression(...) # Here, for example, we set the displacement in the x-direction to zero on the RightTop boundary. bcl = DirichletBC(V.sub(0), Constant(0.0), boundary_markers, physical_tags["RightTop"]) # Prescribe the y-displacement r on both Top and Bottom boundaries. bcb2 = DirichletBC(V.sub(1), r, boundary_markers, physical_tags["Bottom"]) bct2 = DirichletBC(V.sub(1), r, boundary_markers, physical_tags["Top"]) bcs = [bcl, bcb2, bct2] ``` Despite this, unlike the code I wrote for marking boundary subdomains directly in FEniCS, which gives me correct results, the current approach using GMSH to mark boundary subdomains is not converging and is producing incorrect results. Any guidance or suggestions on this issue would be greatly appreciated. ---
**Author:** ![Kei\_Yamamoto](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/kei_yamamoto/32/2904_2.png) [@Kei\_Yamamoto](https://fenicsproject.discourse.group/u/Kei_Yamamoto)\ **Post date:** [March 4, 2025, 9:22pm UTC](https://fenicsproject.discourse.group/t/mesh-markers-in-gmsh/17124/2 "2025-03-04T21:22:36Z")
You posted this question with `dolfinx` tag, but the code seems to be using legacy-dolfin? Either way, I think you need to clarify how you ended up obtaining `crack.xml` in your code. It’s probably better to convert to XDMF file, as explained in the link below. > [@Transitioning from mesh.xml to mesh.xdmf, from dolfin-convert to meshio](https://fenicsproject.discourse.group/t/transitioning-from-mesh-xml-to-mesh-xdmf-from-dolfin-convert-to-meshio/412/3): > > This can be done using meshio and writing the explicit data to xdmf file. import meshio msh = meshio.read("test.msh") meshio.write("mesh.xdmf", meshio.Mesh(points=msh.points, cells={"tetra": msh.cells["tetra"]})) meshio.write("mf.xdmf", meshio.Mesh(points=msh.points, cells={"triangle": msh.cells["triangle"]}, cell\_data={"triangle": {"name\_to\_read": msh.cell\_data["triangle"]["gmsh:physical"]}})) from dolfin import \* mesh = Mesh() with XDMFFile("mesh.xdmf") … ---
**Author:** ![umarkhayaz](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/umarkhayaz/32/6239_2.png) [@umarkhayaz](https://fenicsproject.discourse.group/u/umarkhayaz)\ **Post date:** [March 5, 2025, 7:55pm UTC](https://fenicsproject.discourse.group/t/mesh-markers-in-gmsh/17124/3 "2025-03-05T19:55:37Z")
Yes, I think I accidentally used the wrong tag. The code is actually utilizing legacy-dolfin. In the scenario mentioned, I was creating the `crack.xml` file by running `dolfin-convert crack.msh crack.xml`. Later, I made a few updates to the FEniCS code. Despite these changes, I still wasn’t getting the correct results. As a result, I had to switch to a different approach to generate `crack.xdmf` and `facet_crack.xdmf` using the following code: ```auto import sys filename = "crack" filenamemsh = (filename + ".msh") filenamexdmf = (filename + ".xdmf") filenamexdmf_facet = ("facet_" + filename + ".xdmf") import meshio mesh_from_file = meshio.read(filenamemsh) import numpy def create_mesh(mesh, cell_type, prune_z=False): cells = mesh.get_cells_type(cell_type) cell_data = mesh.get_cell_data("gmsh:physical", cell_type) points = mesh.points[:,:2] if prune_z else mesh.points out_mesh = meshio.Mesh(points=points, cells={cell_type: cells},\ cell_data={"name_to_read":[cell_data]}) return out_mesh line_mesh = create_mesh(mesh_from_file, "line", prune_z=True) meshio.write(filenamexdmf_facet, line_mesh) triangle_mesh = create_mesh(mesh_from_file, "triangle", prune_z=True) meshio.write(filenamexdmf, triangle_mesh) ``` Once the `crack.xdmf` and `facet_crack.xdmf` files were generated, I imported them into the FEniCS code using the following: ```auto # Initialize an empty mesh object mesh = Mesh() meshname = "crack.xdmf" facet_meshname = "facet_crack.xdmf" # Read the .xdmf file data into mesh object with XDMFFile(meshname) as infile: infile.read(mesh) # Extract initial mesh coords x = SpatialCoordinate(mesh) # Read the subdomain data stored in the *.xdmf file mvc = MeshValueCollection("size_t", mesh, 2) with XDMFFile(meshname) as infile: infile.read(mvc, "name_to_read") mf = cpp.mesh.MeshFunctionSizet(mesh, mvc) # Read mesh boundary facets. # Specify parameters of 1-D mesh function object mvc_1d = MeshValueCollection("size_t", mesh, 1) # Read the *.xdmf data into mesh function object with XDMFFile(facet_meshname) as infile: infile.read(mvc_1d, "name_to_read") # Store boundary facets data facets = cpp.mesh.MeshFunctionSizet(mesh, mvc_1d) ds = Measure('ds', domain=mesh, subdomain_data=facets) dx = Measure('dx', domain=mesh, subdomain_data=mf) n=FacetNormal(mesh) ################################################################################## # Define Dirichlet boundary conditions ################################################################################## c=Expression(...) r=Expression(...) bcl= DirichletBC(V.sub(0), Constant(0.0), facets, 20) bcb2 = DirichletBC(V.sub(1), r, facets, 22) bct2 = DirichletBC(V.sub(1), r, facets, 21) bcs = [bcl, bcb2, bct2] ``` In the code, I am computing the J-integral. When I execute the code, the plate doesn’t open along the geometric crack of some length (in a single edge notch specimen) even though the displacement is applied to the top and bottom boundaries. I can observe the deformed shape when visualizing the results in Paraview, but neither the geometric crack open, nor the crack propagates as expected. Additionally, the J-integral values I obtain are significantly larger than the critical energy release rate (Gc), and these values remain constant across all time steps.