Unable to use gmsh generated mesh

teja781 · December 4, 2023, 11:28pm

Hello,

I am trying to work with a 3D mesh generated by gmsh with refinement in the area of interest. gmsh is able to create the mesh successfully but I get an error message while trying to get it to work with fenicsx with the model_to_mesh command of gmshio. Here is the python code and the generated output.

import numpy as np
import gmsh
from dolfinx.io.gmshio import model_to_mesh
import sys
from mpi4py import MPI

##Start geometry and mesh ##
L=1.2 #length
H=1.0 #horizontal dim
a0 = 0.25
h0 = (0.063**2 - (0.063/2)**2)**0.5
pc = 0.5-a0-h0  # pre-crack length
at = pc+a0+h0   # total initial crack length
thickness = 0.5
mrf = 5
lc = 15e-3*mrf
top_marker = 1
bottom_marker = 2
right_marker = 3
left_top_marker = 4
left_bottom_marker = 5
plate_marker = 6

proc = MPI.COMM_WORLD.rank
gmsh.initialize()
if proc == 0:
	gmsh.model.add('fracture_geometry')

	gmsh.model.geo.addPoint(-H/2,-L/2,0,lc,1)
	gmsh.model.geo.addPoint(H/2,-L/2,0,lc,2)
	gmsh.model.geo.addPoint(H/2,L/2,0,lc,3)
	gmsh.model.geo.addPoint(-H/2,L/2,0,lc,4)
	gmsh.model.geo.addPoint(-H/2,0.063/2,0,lc,5)
	gmsh.model.geo.addPoint(-H/2+a0,0.063/2,0,lc,6)
	gmsh.model.geo.addPoint(-H/2+a0+h0,0.0025,0,lc,7)
	gmsh.model.geo.addPoint(-H/2+a0+h0+pc,0.0,0,lc,8)
	gmsh.model.geo.addPoint(-H/2+a0+h0,-0.0025,0,lc,9)
	gmsh.model.geo.addPoint(-H/2+a0,-0.063/2,0,lc,10)
	gmsh.model.geo.addPoint(-H/2,-0.063/2,0,lc,11)

	gmsh.model.geo.addLine(1,2,1)
	gmsh.model.geo.addLine(2,3,2)
	gmsh.model.geo.addLine(3,4,3)
	gmsh.model.geo.addLine(4,5,4)
	gmsh.model.geo.addLine(5,6,5)
	gmsh.model.geo.addLine(6,7,6)
	gmsh.model.geo.addLine(7,8,7)
	gmsh.model.geo.addLine(8,9,8)
	gmsh.model.geo.addLine(9,10,9)
	gmsh.model.geo.addLine(10,11,10)
	gmsh.model.geo.addLine(11,1,11)

	gmsh.model.geo.addCurveLoop([1,2,3,4,5,6,7,8,9,10,11],1)
	pl = gmsh.model.geo.addPlaneSurface([1],1)
	vl = gmsh.model.geo.extrude([(2,1)],0,0,thickness)

	gmsh.model.geo.synchronize()
	# center plane for refinement
	pointEntities = gmsh.model.getEntities(dim=0)
	for p in pointEntities:
		tag = p[1]
		pointCoords = gmsh.model.getValue(0, tag, [])
		if np.isclose(pointCoords,[-H/2+a0+h0+pc,0,thickness]).all():rsp = [tag]
	gmsh.model.geo.addPoint(H/2,0,0,lc,102)
	gmsh.model.geo.addPoint(H/2,0,thickness,lc,103)
	gmsh.model.geo.addLine(8,102,102)
	gmsh.model.geo.addLine(102,103,103)
	gmsh.model.geo.addLine(103,rsp[0],104)
	gmsh.model.geo.addLine(rsp[0],8,105)
	gmsh.model.geo.addCurveLoop([102,103,104,105],200)
	gmsh.model.geo.addPlaneSurface([200],200)

	gmsh.model.geo.synchronize()
	
	gmsh.model.mesh.field.add("Distance", 1)
	gmsh.model.mesh.field.setNumbers(1, "SurfacesList", [200])
	gmsh.model.mesh.field.setNumber(1, "Sampling", 100)
	
	# We then define a `Threshold' field, which uses the return value of the
	# `Distance' field 1 in order to define a simple change in element size
	# depending on the computed distances
	#
	# SizeMax -                     /------------------
	#                              /
	#                             /
	#                            /
	# SizeMin -o----------------/
	#              |                |    |
	#        Point         DistMin  DistMax
	gmsh.model.mesh.field.add("Threshold", 2)
	gmsh.model.mesh.field.setNumber(2, "InField", 1)
	gmsh.model.mesh.field.setNumber(2, "SizeMin", lc/mrf)
	gmsh.model.mesh.field.setNumber(2, "SizeMax", lc)
	gmsh.model.mesh.field.setNumber(2, "DistMin", 0.15)
	gmsh.model.mesh.field.setNumber(2, "DistMax", 0.5)

	gmsh.model.mesh.field.setAsBackgroundMesh(2)
	gmsh.option.setNumber("Mesh.Algorithm", 5)
	gmsh.model.mesh.generate(3)
	surfaceEntities = gmsh.model.getEntities(dim=2)
	for e in surfaceEntities:
	    tag = e[1]
	    nodeTags, nodeCoords, nodeParams = gmsh.model.mesh.getNodes(2, tag)
	    if np.isclose(nodeCoords[0],-0.5) and nodeCoords[1]>0:gmsh.model.addPhysicalGroup(2,[tag],left_top_marker)
	    elif np.isclose(nodeCoords[0],-0.5) and nodeCoords[1]<0:gmsh.model.addPhysicalGroup(2,[tag],left_bottom_marker)
	    else:continue
	gmsh.model.addPhysicalGroup(3,[1],1)

# Launch the GUI to see the results:
if '-nopopup' not in sys.argv:
	gmsh.fltk.run()
fracture_mesh, cell_tags, facet_tags = model_to_mesh(gmsh.model, MPI.COMM_WORLD, 0, 3)
gmsh.finalize()

Info    : Meshing 1D...
Info    : [  0%] Meshing curve 1 (Line)
Info    : [ 10%] Meshing curve 2 (Line)
Info    : [ 10%] Meshing curve 3 (Line)
Info    : [ 10%] Meshing curve 4 (Line)
Info    : [ 20%] Meshing curve 5 (Line)
Info    : [ 20%] Meshing curve 6 (Line)
Info    : [ 20%] Meshing curve 7 (Line)
Info    : [ 20%] Meshing curve 8 (Line)
Info    : [ 30%] Meshing curve 9 (Line)
Info    : [ 30%] Meshing curve 10 (Line)
Info    : [ 30%] Meshing curve 11 (Line)
Info    : [ 30%] Meshing curve 13 (Line)
Info    : [ 40%] Meshing curve 14 (Line)
Info    : [ 40%] Meshing curve 15 (Line)
Info    : [ 40%] Meshing curve 16 (Line)
Info    : [ 50%] Meshing curve 17 (Line)
Info    : [ 50%] Meshing curve 18 (Line)
Info    : [ 50%] Meshing curve 19 (Line)
Info    : [ 50%] Meshing curve 20 (Line)
Info    : [ 60%] Meshing curve 21 (Line)
Info    : [ 60%] Meshing curve 22 (Line)
Info    : [ 60%] Meshing curve 23 (Line)
Info    : [ 60%] Meshing curve 25 (Line)
Info    : [ 70%] Meshing curve 26 (Line)
Info    : [ 70%] Meshing curve 30 (Line)
Info    : [ 70%] Meshing curve 34 (Line)
Info    : [ 80%] Meshing curve 38 (Line)
Info    : [ 80%] Meshing curve 42 (Line)
Info    : [ 80%] Meshing curve 46 (Line)
Info    : [ 80%] Meshing curve 50 (Line)
Info    : [ 90%] Meshing curve 54 (Line)
Info    : [ 90%] Meshing curve 58 (Line)
Info    : [ 90%] Meshing curve 62 (Line)
Info    : [ 90%] Meshing curve 102 (Line)
Info    : [100%] Meshing curve 103 (Line)
Info    : [100%] Meshing curve 104 (Line)
Info    : [100%] Meshing curve 105 (Line)
Info    : Done meshing 1D (Wall 0.0317913s, CPU 0.030543s)
Info    : Meshing 2D...
Info    : [  0%] Meshing surface 1 (Plane, Delaunay)
Info    : [ 10%] Meshing surface 27 (Surface, Delaunay)
Info    : [ 20%] Meshing surface 31 (Surface, Delaunay)
Info    : [ 30%] Meshing surface 35 (Surface, Delaunay)
Info    : [ 30%] Meshing surface 39 (Surface, Delaunay)
Info    : [ 40%] Meshing surface 43 (Surface, Delaunay)
Info    : [ 50%] Meshing surface 47 (Surface, Delaunay)
Info    : [ 50%] Meshing surface 51 (Surface, Delaunay)
Info    : [ 60%] Meshing surface 55 (Surface, Delaunay)
Info    : [ 70%] Meshing surface 59 (Surface, Delaunay)
Info    : [ 80%] Meshing surface 63 (Surface, Delaunay)
Info    : [ 80%] Meshing surface 67 (Surface, Delaunay)
Info    : [ 90%] Meshing surface 68 (Plane, Delaunay)
Info    : [100%] Meshing surface 200 (Plane, Delaunay)
Info    : Done meshing 2D (Wall 0.175926s, CPU 0.169952s)
Info    : Meshing 3D...
Info    : 3D Meshing 1 volume with 1 connected component
Info    : Tetrahedrizing 6883 nodes...
Info    : Done tetrahedrizing 6891 nodes (Wall 0.0669155s, CPU 0.063153s)
Info    : Reconstructing mesh...
Info    :  - Creating surface mesh
Info    :  - Identifying boundary edges
Info    :  - Recovering boundary
Info    : Done reconstructing mesh (Wall 0.163732s, CPU 0.155422s)
Info    : Found volume 1
Info    : It. 0 - 0 nodes created - worst tet radius 18.6038 (nodes removed 0 0)
Info    : It. 500 - 500 nodes created - worst tet radius 3.49053 (nodes removed 0 0)
Info    : It. 1000 - 1000 nodes created - worst tet radius 2.78519 (nodes removed 0 0)
Info    : It. 1500 - 1500 nodes created - worst tet radius 2.47897 (nodes removed 0 0)
Info    : It. 2000 - 2000 nodes created - worst tet radius 2.2595 (nodes removed 0 0)
Info    : It. 2500 - 2500 nodes created - worst tet radius 2.09824 (nodes removed 0 0)
Info    : It. 3000 - 3000 nodes created - worst tet radius 1.98587 (nodes removed 0 0)
Info    : It. 3500 - 3500 nodes created - worst tet radius 1.88906 (nodes removed 0 0)
Info    : It. 4000 - 4000 nodes created - worst tet radius 1.80282 (nodes removed 0 0)
Info    : It. 4500 - 4500 nodes created - worst tet radius 1.73497 (nodes removed 0 0)
Info    : It. 5000 - 5000 nodes created - worst tet radius 1.67743 (nodes removed 0 0)
Info    : It. 5500 - 5500 nodes created - worst tet radius 1.6294 (nodes removed 0 0)
Info    : It. 6000 - 6000 nodes created - worst tet radius 1.58251 (nodes removed 0 0)
Info    : It. 6500 - 6500 nodes created - worst tet radius 1.54318 (nodes removed 0 0)
Info    : It. 7000 - 7000 nodes created - worst tet radius 1.50811 (nodes removed 0 0)
Info    : It. 7500 - 7500 nodes created - worst tet radius 1.47679 (nodes removed 0 0)
Info    : It. 8000 - 8000 nodes created - worst tet radius 1.44479 (nodes removed 0 0)
Info    : It. 8500 - 8500 nodes created - worst tet radius 1.41765 (nodes removed 0 0)
Info    : It. 9000 - 9000 nodes created - worst tet radius 1.391 (nodes removed 0 0)
Info    : It. 9500 - 9500 nodes created - worst tet radius 1.36782 (nodes removed 0 0)
Info    : It. 10000 - 10000 nodes created - worst tet radius 1.34567 (nodes removed 0 0)
Info    : It. 10500 - 10500 nodes created - worst tet radius 1.32442 (nodes removed 0 0)
Info    : It. 11000 - 11000 nodes created - worst tet radius 1.30449 (nodes removed 0 0)
Info    : It. 11500 - 11500 nodes created - worst tet radius 1.28522 (nodes removed 0 0)
Info    : It. 12000 - 12000 nodes created - worst tet radius 1.26796 (nodes removed 0 0)
Info    : It. 12500 - 12500 nodes created - worst tet radius 1.2512 (nodes removed 0 0)
Info    : It. 13000 - 13000 nodes created - worst tet radius 1.23498 (nodes removed 0 0)
Info    : It. 13500 - 13500 nodes created - worst tet radius 1.21898 (nodes removed 0 0)
Info    : It. 14000 - 14000 nodes created - worst tet radius 1.20545 (nodes removed 0 0)
Info    : It. 14500 - 14500 nodes created - worst tet radius 1.19249 (nodes removed 0 0)
Info    : It. 15000 - 15000 nodes created - worst tet radius 1.17991 (nodes removed 0 0)
Info    : It. 15500 - 15500 nodes created - worst tet radius 1.16669 (nodes removed 0 0)
Info    : It. 16000 - 16000 nodes created - worst tet radius 1.15447 (nodes removed 0 0)
Info    : It. 16500 - 16500 nodes created - worst tet radius 1.14309 (nodes removed 0 0)
Info    : It. 17000 - 17000 nodes created - worst tet radius 1.13216 (nodes removed 0 0)
Info    : It. 17500 - 17500 nodes created - worst tet radius 1.12188 (nodes removed 0 0)
Info    : It. 18000 - 18000 nodes created - worst tet radius 1.11118 (nodes removed 0 0)
Info    : It. 18500 - 18500 nodes created - worst tet radius 1.10158 (nodes removed 0 0)
Info    : It. 19000 - 19000 nodes created - worst tet radius 1.09207 (nodes removed 0 0)
Info    : It. 19500 - 19500 nodes created - worst tet radius 1.0831 (nodes removed 0 0)
Info    : It. 20000 - 20000 nodes created - worst tet radius 1.07419 (nodes removed 0 0)
Info    : It. 20500 - 20500 nodes created - worst tet radius 1.06572 (nodes removed 0 0)
Info    : It. 21000 - 21000 nodes created - worst tet radius 1.05682 (nodes removed 0 0)
Info    : It. 21500 - 21500 nodes created - worst tet radius 1.04838 (nodes removed 0 0)
Info    : It. 22000 - 22000 nodes created - worst tet radius 1.04086 (nodes removed 0 0)
Info    : It. 22500 - 22500 nodes created - worst tet radius 1.03249 (nodes removed 0 0)
Info    : It. 23000 - 23000 nodes created - worst tet radius 1.02571 (nodes removed 0 0)
Info    : It. 23500 - 23500 nodes created - worst tet radius 1.01846 (nodes removed 0 0)
Info    : It. 24000 - 24000 nodes created - worst tet radius 1.01119 (nodes removed 0 0)
Info    : It. 24500 - 24500 nodes created - worst tet radius 1.00424 (nodes removed 0 0)
Info    : 3D refinement terminated (33258 nodes total):
Info    :  - 2 Delaunay cavities modified for star shapeness
Info    :  - 0 nodes could not be inserted
Info    :  - 179422 tetrahedra created in 1.20757 sec. (148581 tets/s)
Info    : Done meshing 3D (Wall 1.74878s, CPU 1.72775s)
Info    : Optimizing mesh...
Info    : Optimizing volume 1
Info    : Optimization starts (volume = 0.590953) with worst = 0.00355257 / average = 0.777654:
Info    : 0.00 < quality < 0.10 :       355 elements
Info    : 0.10 < quality < 0.20 :      1100 elements
Info    : 0.20 < quality < 0.30 :      1979 elements
Info    : 0.30 < quality < 0.40 :      3308 elements
Info    : 0.40 < quality < 0.50 :      5090 elements
Info    : 0.50 < quality < 0.60 :      8939 elements
Info    : 0.60 < quality < 0.70 :     18753 elements
Info    : 0.70 < quality < 0.80 :     40008 elements
Info    : 0.80 < quality < 0.90 :     66006 elements
Info    : 0.90 < quality < 1.00 :     33884 elements
Info    : 3394 edge swaps, 111 node relocations (volume = 0.590953): worst = 0.157839 / average = 0.789587 (Wall 0.0837764s, CPU 0.083692s)
Info    : 3413 edge swaps, 111 node relocations (volume = 0.590953): worst = 0.300068 / average = 0.789657 (Wall 0.111156s, CPU 0.111071s)
Info    : No ill-shaped tets in the mesh :-)
Info    : 0.00 < quality < 0.10 :         0 elements
Info    : 0.10 < quality < 0.20 :         0 elements
Info    : 0.20 < quality < 0.30 :         0 elements
Info    : 0.30 < quality < 0.40 :      3240 elements
Info    : 0.40 < quality < 0.50 :      4925 elements
Info    : 0.50 < quality < 0.60 :      8691 elements
Info    : 0.60 < quality < 0.70 :     18549 elements
Info    : 0.70 < quality < 0.80 :     40681 elements
Info    : 0.80 < quality < 0.90 :     66635 elements
Info    : 0.90 < quality < 1.00 :     33626 elements
Info    : Done optimizing mesh (Wall 0.35338s, CPU 0.345206s)
Info    : 33258 nodes 193665 elements
X_ChangeProperty: BadValue (integer parameter out of range for operation) 0x0
-------------------------------------------------------
Version       : 4.11.1-git-4608e03
License       : GNU General Public License
Build OS      : Linux64-sdk
Build date    : 20231005
Build host    : buildkitsandbox
Build options : 64Bit ALGLIB[contrib] ANN[contrib] Bamg Blossom DIntegration Dlopen DomHex Eigen[contrib] Fltk Gmm[contrib] Hxt Kbipack LinuxJoystick MathEx[contrib] Mesh Metis[contrib] Mpeg Netgen ONELAB ONELABMetamodel OpenCASCADE OpenCASCADE-CAF OpenGL OpenMP OptHom Parser Plugins Png Post QuadMeshingTools QuadTri Solver TetGen/BR Voro++[contrib] WinslowUntangler Zlib
FLTK version  : 1.3.8
OCC version   : 7.5.1
Packaged by   : root
Web site      : https://gmsh.info
Issue tracker : https://gitlab.onelab.info/gmsh/gmsh/issues
-------------------------------------------------------
python3: /src/dolfinx/cpp/dolfinx/common/MPI.h:117: constexpr int dolfinx::MPI::index_owner(int, std::size_t, std::size_t): Assertion `index < N' failed.

Loguru caught a signal: SIGABRT
Stack trace:
25      0x55b2a2d12ba5 _start + 37
24      0x7fd3a763ce40 __libc_start_main + 128
23      0x7fd3a763cd90 /lib/x86_64-linux-gnu/libc.so.6(+0x29d90) [0x7fd3a763cd90]
22      0x55b2a2d12cad Py_BytesMain + 45
21      0x55b2a2d3c41e Py_RunMain + 702
20      0x55b2a2d49653 _PyRun_AnyFileObject + 67
19      0x55b2a2d49a08 _PyRun_SimpleFileObject + 424
18      0x55b2a2d4a525 /usr/bin/python3(+0x264525) [0x55b2a2d4a525]
17      0x55b2a2d440bb /usr/bin/python3(+0x25e0bb) [0x55b2a2d440bb]
16      0x55b2a2d4a7d8 /usr/bin/python3(+0x2647d8) [0x55b2a2d4a7d8]
15      0x55b2a2d1fcf6 PyEval_EvalCode + 134
14      0x55b2a2d1fe56 /usr/bin/python3(+0x239e56) [0x55b2a2d1fe56]
13      0x55b2a2c2ee0d _PyEval_EvalFrameDefault + 1725
12      0x55b2a2c4670c _PyFunction_Vectorcall + 124
11      0x55b2a2c351f1 _PyEval_EvalFrameDefault + 27297
10      0x55b2a2c3c5eb _PyObject_MakeTpCall + 603
9       0x55b2a2c45e0e /usr/bin/python3(+0x15fe0e) [0x55b2a2c45e0e]
8       0x7fd39369cae8 /usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/cpp.cpython-310-x86_64-linux-gnu.so(+0xb2ae8) [0x7fd39369cae8]
7       0x7fd3937f26c8 /usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/cpp.cpython-310-x86_64-linux-gnu.so(+0x2086c8) [0x7fd3937f26c8]
6       0x7fd39352079f dolfinx::io::xdmf_utils::distribute_entity_data(dolfinx::mesh::Topology const&, std::vector<long> const&, long, dolfinx::fem::ElementDofLayout const&, std::mdspan<int const, std::extents<unsigned long, 18446744073709551615ul, 18446744073709551615ul>, std::layout_right, std::default_accessor<int const> >, int, std::span<long const, 18446744073709551615ul>, std::span<int const, 18446744073709551615ul>) + 11151
5       0x7fd3a764ce96 /lib/x86_64-linux-gnu/libc.so.6(+0x39e96) [0x7fd3a764ce96]
4       0x7fd3a763b71b /lib/x86_64-linux-gnu/libc.so.6(+0x2871b) [0x7fd3a763b71b]
3       0x7fd3a763b7f3 abort + 211
2       0x7fd3a7655476 raise + 22
1       0x7fd3a76a99fc pthread_kill + 300
0       0x7fd3a7655520 /lib/x86_64-linux-gnu/libc.so.6(+0x42520) [0x7fd3a7655520]
2023-12-04 18:12:44.038 (   6.791s) [main            ]                       :0     FATL| Signal: SIGABRT
Aborted (core dumped)

Could I please get some help in understanding what the issue is?

Thank you!

francesco-ballarin · December 5, 2023, 6:10am

I think that (at least one of) the issue is calling gmsh.model.addPhysicalGroup after gmsh.model.mesh.generate.

If I change the final part of the if statement to be:

    ....

    gmsh.model.mesh.field.setAsBackgroundMesh(2)
    gmsh.option.setNumber("Mesh.Algorithm", 5)

    """
    surfaceEntities = gmsh.model.getEntities(dim=2)
    for e in surfaceEntities:
        tag = e[1]
        nodeTags, nodeCoords, nodeParams = gmsh.model.mesh.getNodes(2, tag)
        if np.isclose(nodeCoords[0],-0.5) and nodeCoords[1]>0:gmsh.model.addPhysicalGroup(2,[tag],left_top_marker)
        elif np.isclose(nodeCoords[0],-0.5) and nodeCoords[1]<0:gmsh.model.addPhysicalGroup(2,[tag],left_bottom_marker)
        else:continue
    """

    gmsh.model.addPhysicalGroup(3,[1],1)
    gmsh.model.mesh.generate(3)

I do not get any error, while I do get an error if swapping the last two lines in my snippet.

Unfortunately, requiring gmsh.model.mesh.generate to be the last call means that I had to comment the part of your code related to gmsh.model.addPhysicalGroup(2, because one of gmsh.model.getEntities or gmsh.model.mesh.getNodes seems to require the mesh to be available. I don’t understand what is supposed to happen there, but I would suggest to have a look if it was possible to rewrite this in an equivalent way that makes use of the ids you get from gmsh.model.geo.add* calls.

teja781 · December 5, 2023, 5:51pm

Thank you for your reply! In the portion that you commented out I am adding physicalGroups using the node coordinates on the surface which is why I put them after generating the mesh. It seems like a strange limitation that I won’t be able to do that for the mesh to be compatible with dolfinx.

I’ve changed the code to add all the physical groups using bounding_box of surfaces, before mesh generation and it works now.

Here is the working code.

import numpy as np
import gmsh
from dolfinx.io.gmshio import model_to_mesh
from dolfinx import io
import sys
from mpi4py import MPI

##Start geometry and mesh ##
L=1.2 #length
H=1.0 #horizontal dim
a0 = 0.25
h0 = (0.063**2 - (0.063/2)**2)**0.5
pc = 0.5-a0-h0  # pre-crack length
at = pc+a0+h0   # total initial crack length
thickness = 0.5
mrf = 5
lc = 15e-3*mrf
top_marker = 1
bottom_marker = 2
right_marker = 3
left_top_marker = 4
left_bottom_marker = 5
plate_marker = 6

proc = MPI.COMM_WORLD.rank
gmsh.initialize()
if proc == 0:
    gmsh.model.add('fracture_geometry')

    gmsh.model.geo.addPoint(-H/2,-L/2,0,lc,1)
    gmsh.model.geo.addPoint(H/2,-L/2,0,lc,2)
    gmsh.model.geo.addPoint(H/2,L/2,0,lc,3)
    gmsh.model.geo.addPoint(-H/2,L/2,0,lc,4)
    gmsh.model.geo.addPoint(-H/2,0.063/2,0,lc,5)
    gmsh.model.geo.addPoint(-H/2+a0,0.063/2,0,lc,6)
    gmsh.model.geo.addPoint(-H/2+a0+h0,0.0025,0,lc,7)
    gmsh.model.geo.addPoint(-H/2+a0+h0+pc,0.0,0,lc,8)
    gmsh.model.geo.addPoint(-H/2+a0+h0,-0.0025,0,lc,9)
    gmsh.model.geo.addPoint(-H/2+a0,-0.063/2,0,lc,10)
    gmsh.model.geo.addPoint(-H/2,-0.063/2,0,lc,11)

    gmsh.model.geo.addLine(1,2,1)
    gmsh.model.geo.addLine(2,3,2)
    gmsh.model.geo.addLine(3,4,3)
    gmsh.model.geo.addLine(4,5,4)
    gmsh.model.geo.addLine(5,6,5)
    gmsh.model.geo.addLine(6,7,6)
    gmsh.model.geo.addLine(7,8,7)
    gmsh.model.geo.addLine(8,9,8)
    gmsh.model.geo.addLine(9,10,9)
    gmsh.model.geo.addLine(10,11,10)
    gmsh.model.geo.addLine(11,1,11)

    gmsh.model.geo.addCurveLoop([1,2,3,4,5,6,7,8,9,10,11],1)
    pl = gmsh.model.geo.addPlaneSurface([1],1)
    vl = gmsh.model.geo.extrude([(2,1)],0,0,thickness)

    gmsh.model.geo.synchronize()

    surfaceEntities = gmsh.model.getEntities(dim=2)
    for e in surfaceEntities:
        boxCoords = gmsh.model.getBoundingBox(e[0], e[1])
        if np.isclose(boxCoords[::3],-0.5).all() and (np.array(boxCoords[1::3])>0).all():
            gmsh.model.addPhysicalGroup(2,[e[1]],left_top_marker)
        elif np.isclose(boxCoords[::3],-0.5).all() and (np.array(boxCoords[1::3])<0).all():
            gmsh.model.addPhysicalGroup(2,[e[1]],left_bottom_marker)
        else:continue
    gmsh.model.addPhysicalGroup(3,[1],1)

    # center plane for refinement
    pointEntities = gmsh.model.getEntities(dim=0)
    for p in pointEntities:
        tag = p[1]
        pointCoords = gmsh.model.getValue(0, tag, [])
        if np.isclose(pointCoords,[-H/2+a0+h0+pc,0,thickness]).all():rsp = [tag]
    gmsh.model.geo.addPoint(H/2,0,0,lc,102)
    gmsh.model.geo.addPoint(H/2,0,thickness,lc,103)
    gmsh.model.geo.addLine(8,102,102)
    gmsh.model.geo.addLine(102,103,103)
    gmsh.model.geo.addLine(103,rsp[0],104)
    gmsh.model.geo.addLine(rsp[0],8,105)
    gmsh.model.geo.addCurveLoop([102,103,104,105],200)
    gmsh.model.geo.addPlaneSurface([200],200)

    gmsh.model.geo.synchronize()
    
    gmsh.model.mesh.field.add("Distance", 1)
    gmsh.model.mesh.field.setNumbers(1, "SurfacesList", [200])
    gmsh.model.mesh.field.setNumber(1, "Sampling", 100)
    
    # We then define a `Threshold' field, which uses the return value of the
    # `Distance' field 1 in order to define a simple change in element size
    # depending on the computed distances
    #
    # SizeMax -                     /------------------
    #                              /
    #                             /
    #                            /
    # SizeMin -o----------------/
    #              |                |    |
    #        Point         DistMin  DistMax
    gmsh.model.mesh.field.add("Threshold", 2)
    gmsh.model.mesh.field.setNumber(2, "InField", 1)
    gmsh.model.mesh.field.setNumber(2, "SizeMin", lc/mrf)
    gmsh.model.mesh.field.setNumber(2, "SizeMax", lc)
    gmsh.model.mesh.field.setNumber(2, "DistMin", 0.15)
    gmsh.model.mesh.field.setNumber(2, "DistMax", 0.5)


    gmsh.model.mesh.field.setAsBackgroundMesh(2)
    gmsh.option.setNumber("Mesh.Algorithm", 5)
    gmsh.model.mesh.generate(3)

    gmsh.write(f'fracture_geometry.msh')

fracture_mesh, cell_tags, facet_tags = model_to_mesh(gmsh.model, MPI.COMM_WORLD, 0, 3)
gmsh.finalize()

dokken · December 5, 2023, 6:17pm

The dolfinx.gmshio uses Gmsh Python API to extract the mesh and corresponding markers.
This is then a limitation in how Gmsh treats physical groups.

Since you are doing this directly with mesh nodes, why not do it once the mesh has been read into dolfinx instead?

teja781 · December 5, 2023, 6:36pm

Yes, that would be a more general solution for something like this. Since I’ve already managed to solve this issue for now, I’ll keep that in mind for the future.

Thank you!

Topic		Replies	Views
Create a mesh with Gmsh General	2	287	September 20, 2023
Problem with `model_to_mesh` in `dolfinx.io.gmshio` mesh	2	252	March 25, 2024
Problem reading mesh using meshio or interfacing msh with dolfinx dolfinx mesh	2	344	December 15, 2023
Mesh 2D gmsh to fenics legacy (2019.2) mesh	3	32	September 27, 2024
Dolfin:Problem importing a 3D mesh generated with gmsh mesh gmsh	9	521	January 2, 2024

Unable to use gmsh generated mesh

Related topics