Adding a crack seems to mess up the node ordering. Consider the indices array. The values there range range from 0 to 50, while you only have 43 points in points.shape.
You would need to renumber the topology.
I’ve done this below.
import gmsh
import numpy as np
from mpi4py import MPI
import dolfinx
import dolfinx.io.gmshio as gmshio
from dolfinx.io.utils import distribute_entity_data
L = 1
H = 0.1
l = 0.05 #length of the crack
lc = 0.01 # mesh size around the crack
tdim = 2 #domain dimension
gdim = 2
# Initialise gmsh and set options
gmsh.initialize()
model = gmsh.model()
# surface without crack
s = model.occ.addRectangle(0.0, 0.0, 0.0, H, L)
# crack
p4 = model.occ.addPoint(H / 2, L / 2 - l / 2, 0, lc, tag=10)
p5 = model.occ.addPoint(H / 2, L / 2 + l / 2, 0, lc, tag=11)
crack = model.occ.addLine(p4, p5, tag=12)
o, m = gmsh.model.occ.fragment([(tdim, s)], [(tdim - 1, crack)])
model.occ.synchronize()
new_surf = m[0][0][1]
model.addPhysicalGroup(tdim, [new_surf], 100)
model.setPhysicalName(tdim, new_surf, "Rectangular surface")
model.addPhysicalGroup(tdim - 1, [12], tag=12)
model.setPhysicalName(tdim - 1, 12, "crack")
gmsh.model.addPhysicalGroup(0, [p4], 13)
model.mesh.generate(tdim)
indices1, points1, _ = model.mesh.getNodes()
gmsh.plugin.setNumber("Crack", "Dimension", tdim - 1)
gmsh.plugin.setNumber("Crack", "PhysicalGroup", 12)
gmsh.plugin.setNumber("Crack", "DebugView", tdim - 1)
gmsh.plugin.setNumber("Crack", "OpenBoundaryPhysicalGroup", 13)
gmsh.plugin.run("Crack")
indices, points, _ = model.mesh.getNodes()
points = points.reshape(-1, 3)
# Gmsh indices starts at 1. We therefore subtract one to use
# zero-based numbering
indices -= 1
# In some cases, Gmsh does not return the points in the same
# order as their unique node index. We therefore sort nodes in
# geometry according to the unique index
perm_sort = np.argsort(indices)
points = points[perm_sort]
# Remap missing points to a contiguous index
index_to_renumbered = np.full(int(indices[perm_sort[-1]]+1), -1, dtype=np.int64)
index_to_renumbered[indices[perm_sort]] = np.arange(len(indices))
topologies = gmshio.extract_topology_and_markers(model)
comm = MPI.COMM_WORLD
rank = 0
for key, val in topologies.items():
val["topology"] = index_to_renumbered[val["topology"]]
# Extract Gmsh cell id, dimension of cell and number of nodes to
# cell for each
num_cell_types = len(topologies.keys())
cell_information = dict()
cell_dimensions = np.zeros(num_cell_types, dtype=np.int32)
for i, element in enumerate(topologies.keys()):
_, dim, _, num_nodes, _, _ = model.mesh.getElementProperties(element)
cell_information[i] = {"id": element, "dim": dim, "num_nodes": num_nodes}
cell_dimensions[i] = dim
# Sort elements by ascending dimension
perm_sort = np.argsort(cell_dimensions)
# Broadcast cell type data and geometric dimension
cell_id = cell_information[perm_sort[-1]]["id"]
tdim = cell_information[perm_sort[-1]]["dim"]
num_nodes = cell_information[perm_sort[-1]]["num_nodes"]
cell_id, num_nodes = comm.bcast([cell_id, num_nodes], root=rank)
# Check for facet data and broadcast relevant info if True
has_facet_data = False
if tdim - 1 in cell_dimensions:
has_facet_data = True
dtype = dolfinx.default_real_type
has_facet_data = comm.bcast(has_facet_data, root=rank)
if has_facet_data:
num_facet_nodes = comm.bcast(cell_information[perm_sort[-2]]["num_nodes"], root=rank)
gmsh_facet_id = cell_information[perm_sort[-2]]["id"]
marked_facets = np.asarray(topologies[gmsh_facet_id]["topology"], dtype=np.int64)
facet_values = np.asarray(topologies[gmsh_facet_id]["cell_data"], dtype=np.int32)
cells = np.asarray(topologies[cell_id]["topology"], dtype=np.int64)
cell_values = np.asarray(topologies[cell_id]["cell_data"], dtype=np.int32)
else:
cell_id, num_nodes = comm.bcast([None, None], root=rank)
cells, x = np.empty([0, num_nodes], dtype=np.int32), np.empty([0, gdim], dtype=dtype)
cell_values = np.empty((0,), dtype=np.int32)
has_facet_data = comm.bcast(None, root=rank)
if has_facet_data:
num_facet_nodes = comm.bcast(None, root=rank)
marked_facets = np.empty((0, num_facet_nodes), dtype=np.int32)
facet_values = np.empty((0,), dtype=np.int32)
# Create distributed mesh
ufl_domain = gmshio.ufl_mesh(cell_id, gdim, dtype=dtype)
gmsh_cell_perm = gmshio.cell_perm_array(dolfinx.cpp.mesh.to_type(str(ufl_domain.ufl_cell())), num_nodes)
cells = cells[:, gmsh_cell_perm].copy()
mesh = dolfinx.mesh.create_mesh(comm, cells, points[:, :gdim].astype(dtype, copy=False), ufl_domain)
# Create MeshTags for cells
local_entities, local_values = distribute_entity_data(
mesh._cpp_object, mesh.topology.dim, cells, cell_values
)
mesh.topology.create_connectivity(mesh.topology.dim, 0)
adj = dolfinx.cpp.graph.AdjacencyList_int32(local_entities)
ct = dolfinx.mesh.meshtags_from_entities(
mesh, mesh.topology.dim, adj, local_values.astype(np.int32, copy=False)
)
ct.name = "Cell tags"
# Create MeshTags for facets
topology = mesh.topology
tdim = topology.dim
if has_facet_data:
# Permute facets from MSH to DOLFINx ordering
# FIXME: This does not work for prism meshes
if topology.cell_type == dolfinx.mesh.CellType.prism or topology.cell_type == dolfinx.mesh.CellType.pyramid:
raise RuntimeError(f"Unsupported cell type {topology.cell_type}")
facet_type = dolfinx.cpp.mesh.cell_entity_type(
dolfinx.cpp.mesh.to_type(str(ufl_domain.ufl_cell())), tdim - 1, 0
)
gmsh_facet_perm = gmshio.cell_perm_array(facet_type, num_facet_nodes)
marked_facets = marked_facets[:, gmsh_facet_perm]
local_entities, local_values = distribute_entity_data(
mesh._cpp_object, tdim - 1, marked_facets, facet_values
)
mesh.topology.create_connectivity(topology.dim - 1, tdim)
adj = dolfinx.cpp.graph.AdjacencyList_int32(local_entities)
ft = dolfinx.mesh.meshtags_from_entities(mesh, tdim - 1, adj, local_values.astype(np.int32, copy=False))
ft.name = "Facet tags"
else:
ft = dolfinx.mesh.meshtags(mesh, tdim - 1, np.empty(0, dtype=np.int32), np.empty(0, dtype=np.int32))
import ufl
print("Crack surface", dolfinx.fem.assemble_scalar(dolfinx.fem.form(1*ufl.ds(domain=mesh, subdomain_data=ft, subdomain_id=12))),2*l)
print("vol", dolfinx.fem.assemble_scalar(dolfinx.fem.form(1*ufl.dx(domain=mesh))),L*H)
print(dolfinx.fem.assemble_scalar(dolfinx.fem.form(1*ufl.ds(domain=mesh))),2*(L+H+l))
# b_indices = dolfinx.mesh.exterior_facet_indices(mesh.topology)
# ft = dolfinx.mesh.meshtags(mesh, mesh.topology.dim-1, b_indices, np.arange(len(b_indices), dtype=np.int32))
with dolfinx.io.XDMFFile(MPI.COMM_WORLD, "mesh.xdmf", "w") as xdmf:
xdmf.write_mesh(mesh)
xdmf.write_meshtags(ct, mesh.geometry)
xdmf.write_meshtags(ft, mesh.geometry)
It seems like the crack in the mesh is slightly bigger than what has been marked with 12.
I do not have time to debug this further today.