Mesh Markers in GMSH

variational formulation
, , , ,
1

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:



// ------------------------------
// 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:



# 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.

2

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.

3

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:

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:

# 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.