Mshr is not working

General
1

Hello All,

It seems that mshr is no longer supported. I am interested in a simple 2D mesh of a rectangle - circle. After some search, I found that one can use pygmsh and meshio. I have installed both of them in my env, but I am not sure how to proceed further. Any guidelines would be highly appreciated. I believe this can be done in a few lines, but I am totally new to pygmsh and meshio, and how can I generate something compatible with FENICS 2019

2

@dokken has some wonderful tutorials on this topic: Tutorials | Jørgen S. Dokken

1 Like
3

Thanks, @nate!

I am trying to reproduce the following

domain = Rectangle(Point(0., 0.), Point(Length, Length)) - Circle(Point(0., 0.), Radius, 100)
mesh   = generate_mesh(domain, mesh_density)
dimens = mesh.geometry().dim()

This was done in mshr, but I am still unable to reproduce it using pygmsh and meshio.

4

Without posting the code that doesn’t work as expected, you are unlikely to get any help.
Please post your attempted solution.

1 Like
5

@dokken, thanks for the reply and tutorials!

I am trying to generate a mesh for the top-right quarter of “a plate with a hole in the center.” My attempt was based on your 2D tutorial.

This was my attempt:

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

def generate_geom():
    resolution = 0.01
    # Channel parameters
    L = 1.0
    H = 1.0
    c = [0.0, 0.0, 0]
    r = 0.1

    # Initialize empty geometry using the build in kernel in GMSH
    geometry = pg.geo.Geometry()
    # Fetch model we would like to add data to
    model = geometry.__enter__()
    # Add circle
    circle = model.add_circle(c, r, mesh_size=3*resolution)

    points = [model.add_point((0, 0, 0), mesh_size=resolution),
              model.add_point((L, 0, 0), mesh_size=resolution),
              model.add_point((L, H, 0), mesh_size=resolution),
              model.add_point((0, H, 0), mesh_size=resolution)]

    # Add lines between all points creating the rectangle
    rectangle_lines = [model.add_line(points[i], points[i+1])
                    for i in range(-1, len(points)-1)]

    # Create a line loop and plane surface for meshing
    rectangle_loop = model.add_curve_loop(rectangle_lines)
    plane_surface  = model.add_plane_surface(rectangle_loop, holes=[circle.curve_loop])

    # Call gmsh kernel before add physical entities
    model.synchronize()

    model.add_physical([plane_surface], "Plate")
    model.add_physical([rectangle_lines[0]], "left")
    model.add_physical([rectangle_lines[1]], "bottom")
    model.add_physical([rectangle_lines[2]], "right")
    model.add_physical([rectangle_lines[3]], "top")

    return model


geometry = generate_geom()
geometry.generate_mesh(dim=2)
gmsh.write("mesh.msh")
gmsh.clear()
geometry.__exit__()
mesh_from_file = meshio.read("mesh.msh")

triangle_mesh = create_mesh(mesh_from_file, "triangle", prune_z=True)
meshio.write("mesh.xdmf", triangle_mesh)

I got the following error:

Exception: Unable to recover the edge 43704 (100/148) on curve 4 (on surface 2)

Also, how can I visualize it in order to see if it is doing the right thing? And before that I think I am missing a command to define it as a mesh compatible with FEniCS, i.e., what would be the FEniCS command to define the mesh for the FE analysis?

Best

6

@dokken,

In the tutorial, you have two meshes: line_mesh and triangle_mesh. Why is that? Do you need both of them for the FE analysis?

7

First of all, you can open an msh file with gmsh, and look at the mesh.
Next, you can look at the mesh.xdmf file with Paraview, to see if the conversion has worked.

I have these two meshes as I want to use tagged facets from the mesh in my analysis (for boundary conditions).

These can be read into DOLFIN as a MeshValueCollection or `MeshFunction’, see various other posts on the forum regarding this.

As you can see from the error message

  File "/root/shared/mwe.py", line 57, in <module>
    geometry.generate_mesh(dim=2)
  File "/usr/local/lib/python3.10/dist-packages/pygmsh/common/geometry.py", line 374, in generate_mesh
    gmsh.model.mesh.generate(dim)
  File "/usr/local/lib/gmsh.py", line 1973, in generate
    raise Exception(logger.getLastError())
Exception: Unable to recover the edge 43704 (100/148) on curve 4 (on surface 2)

this is a problem when generating the mesh with gmsh.

You can add gmsh.fltk.run() before geometry.generate_mesh(dim=2) to inspect the objects you have created.

In general, I no longer use pygmsh, and prefer to stick to the GMSH Python API, see for instance: Using the GMSH Python API to generate complex meshes | Jørgen S. Dokken

For a simple 2D geometry with a hole in it, see: https://github.com/jorgensd/dolfinx_ipcs/blob/main/create_and_convert_2D_mesh.py

1 Like
8

Many thanks, @dokken

I should explore gmsh python api then. For the time being, do you think there is a quick fix for this? I added the command you suggested and got the following:

image
image924Ă—918 3.24 KB

Obviously, it needs to be trimmed.

Thanks again

9

As I do not use pygmsh much I cannot help you much with that. That is why i linked you to cases that uses the GMSH Python API directly, which also introduces holes in the domain.

1 Like