Crack lips closed in 3D
Hi Christophe,
I was trying to generate a mesh with an edge crack in 3D following the 2D example. Somehow even after marking the OpenBoundaryPhysicalGroup
's the crack lips are still coming out as sealed. The approach I follow is to generate a box, generate a rectangle to fragment the box (after rotating it) and then generate the crack seam.
Here is the relevant code in 2d:
import gmsh
import sys
gmsh.initialize(sys.argv)
gmsh.model.add("square with edge crack")
gmsh.option.set_number("Mesh.CharacteristicLengthMin", 0.01)
gmsh.option.set_number("Mesh.CharacteristicLengthMax", 10)
gmsh.option.setNumber('General.NumThreads', 5)
surf1 = 1
L = 1
H = 1
c = 0.5
gmsh.model.occ.addRectangle(0, 0, 0, L, H, surf1)
pt1 = gmsh.model.occ.addPoint(0.0, H/2., 0)
pt2 = gmsh.model.occ.addPoint(c, H/2., 0)
line1 = gmsh.model.occ.addLine(pt1, pt2)
o, m = gmsh.model.occ.fragment([(2, surf1)], [(1, line1)])
gmsh.model.occ.synchronize()
eps=1.e-3
bdry = gmsh.model.get_entities_in_bounding_box(-eps, -eps, -eps, eps, 1+eps, eps, dim=1)
bdry_tags = [bdry[i][1] for i in range(len(bdry))] # should have 2 because of the cut line
bdry_points = gmsh.model.get_entities_in_bounding_box(-eps, H/2.-eps, -eps, eps, H/2.+eps, eps, dim=0)
# m contains, for each input entity (surf1, line1), the child entities
# (if any) after the fragmentation, as lists of tuples. To apply the crack
# plugin we group all the intersecting lines in a physical group
new_surf = m[0][0][1]
new_lines = [item[1] for sublist in m[1:] for item in sublist]
gmsh.model.addPhysicalGroup(2, [new_surf], 100)
gmsh.model.addPhysicalGroup(1, new_lines, 101)
# gmsh.model.addPhysicalGroup(1, bdry_tags, 102)
gmsh.model.addPhysicalGroup(0, [bdry_points[0][1]], 102)
gmsh.model.mesh.generate(dim=2)
gmsh.plugin.setNumber("Crack", "PhysicalGroup", 101)
gmsh.plugin.setNumber("Crack", "OpenBoundaryPhysicalGroup", 102)
gmsh.plugin.run("Crack")
name = 'square_with_edge_crack.msh'
gmsh.write(name)
gmsh.finalize()
# read and export to xdmf
import meshio
mesh = meshio.read(name)
points, cells = mesh.points, mesh.get_cells_type("triangle")
line_cells = mesh.get_cells_type("line")
cell_data_tri = mesh.get_cell_data("gmsh:physical", "triangle")
cell_data_line = mesh.get_cell_data("gmsh:physical", "line")
tri_mesh = meshio.Mesh(points[:, :2], cells={"triangle": cells}, cell_data={"domain": [cell_data_tri]})
line_mesh = meshio.Mesh(points[:, :2], cells={"line": line_cells}, cell_data={"domain_bdry": [cell_data_line]})
meshio.xdmf.write("square_with_edge_crack.xdmf", tri_mesh)
meshio.xdmf.write("square_with_edge_crack_line.xdmf", line_mesh)
which worked perfectly
whereas for the 3D mesh
from numpy import deg2rad
import gmsh
import sys
gmsh.initialize(sys.argv)
gmsh.model.add("cuboid with edge crack")
gmsh.option.set_number("Mesh.CharacteristicLengthMin", 0.01)
gmsh.option.set_number("Mesh.CharacteristicLengthMax", 10)
gmsh.option.setNumber('General.NumThreads', 5)
surf1 = 1
L = 1
H = 1
thick = 0.05 # dim/20 just for check
c = 0.5
gmsh.model.occ.add_box(0, 0, 0, L, H, thick, tag=1)
eps_length = 1.e-3
# in 3D you should add a rectangular facet for the crack, right?
rectangle1 = gmsh.model.occ.addRectangle(0, H/2. - thick/2., thick/2., c, thick)
gmsh.model.occ.rotate([(2, rectangle1)], 0., H/2., thick/2., 1.0, 0.0, 0.0, deg2rad(90)) # this should work right?
o, m = gmsh.model.occ.fragment([(3, surf1)], [(2, rectangle1)])
gmsh.model.occ.synchronize()
eps=1.e-3
bdry = gmsh.model.get_entities_in_bounding_box(-eps, -eps, -eps, eps, H + eps, thick + eps, dim=3)
bdry_tags = [bdry[i][1] for i in range(len(bdry))] # should have 2 because of the cut rectangle
bdry_points = gmsh.model.get_entities_in_bounding_box(-eps, H/2.-eps, -eps, c, H/2.+eps, thick+eps, dim=0)
bdry_lines = gmsh.model.get_entities_in_bounding_box(-eps, H/2.-eps, -eps, c, H/2.+eps, thick+eps, dim=1)
# m contains, for each input entity (surf1, line1), the child entities
# (if any) after the fragmentation, as lists of tuples. To apply the crack
# plugin we group all the intersecting lines in a physical group
new_vol = m[0][0][1]
new_surfs = [item[1] for sublist in m[1:] for item in sublist]
gmsh.model.addPhysicalGroup(3, [new_vol], 100)
gmsh.model.addPhysicalGroup(2, new_surfs, 101)
# gmsh.model.addPhysicalGroup(1, bdry_tags, 102)
gmsh.model.addPhysicalGroup(1, [bdry_lines[0][1]], 102)
gmsh.model.addPhysicalGroup(0, [bdry_points[0][1]], 103)
gmsh.model.mesh.generate(3)
gmsh.plugin.setNumber("Crack", "Dimension", 2)
gmsh.plugin.setNumber("Crack", "PhysicalGroup", 101)
gmsh.plugin.setNumber("Crack", "OpenBoundaryPhysicalGroup", 102)
# gmsh.plugin.setNumber("Crack", "OpenBoundaryPhysicalGroup", 103)
gmsh.plugin.run("Crack")
name = 'cuboid_with_edge_crack.msh'
gmsh.write(name)
gmsh.finalize()
it seems that the crack lips are closed (the image is basically a slice through the thickness of the mesh).
Could you help me pin down the issue? I also took a look at the crack3d.py
tutorial but can't seem to figure out how to mark the surfaces/lines properly to open the crack.
Thanks