Error in using meshio.write() and create_mesh() for a mesh imported from gmsh

1

As an MWE, I came up with this simple problem:

I created a simple geometry in gmsh, imported it to FEniCS. The geometry looks like below (followed by the gmsh script). I extracted the annulus part and tried to create xdmf files for the cells and the boundaries of the annulus using meshio.write().

Screenshot (118)
Screenshot (118)1257Ă—707 41.3 KB

gmsh script:

SetFactory("OpenCASCADE");
//+
Circle(1) = {0, 0, 0, 1, 0, 2*Pi};
//+
Circle(2) = {0, 0, 0, 2, 0, 2*Pi};
//+
Curve Loop(1) = {1};
//+
Plane Surface(1) = {1};
//+
Curve Loop(2) = {2};
//+
Curve Loop(3) = {1};
//+
Plane Surface(2) = {2, 3};
//+
Physical Curve("C1", 1) = {1};
//+
Physical Curve("C2", 2) = {2};
//+
Physical Surface("S1", 3) = {1};
//+
Physical Surface("S2", 4) = {2};

Python script:

from dolfin import *     
import meshio
from dolfin import Mesh, XDMFFile, File, MeshValueCollection, cpp
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)
    out_mesh = meshio.Mesh(points=mesh.points, cells={
                           cell_type: cells}, cell_data={"name_to_read": [cell_data]})
    if prune_z:
        out_mesh.prune_z_0()
    return out_mesh
#--------------------------------------------------------------IMPORTING GEOMETRY FROM GMSH-----------------------------------------------------------------------------------------------

msh = meshio.read("annulus_mesh.msh")
triangle_mesh = create_mesh(msh, "triangle", True)
line_mesh = create_mesh(msh, "line", True)
meshio.write("mesh.xdmf", triangle_mesh)
meshio.write("mf.xdmf", line_mesh) 
from dolfin import *
mesh = Mesh()
mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim())
with XDMFFile("mesh.xdmf") as infile:
   infile.read(mesh)
   infile.read(mvc, "name_to_read")
cf = cpp.mesh.MeshFunctionSizet(mesh, mvc)

mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim()-1)
with XDMFFile("mf.xdmf") as infile:
    infile.read(mvc, "name_to_read")
mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)


# Extracting the annulus
annulus = SubMesh(mesh,cf,4)


xdmf = XDMFFile("annulus_mesh.xdmf")
xdmf.write(annulus)
xdmf.close()

annulus_mesh = meshio.read("annulus_mesh.xdmf")
meshio.write("mesh_annulus.msh", annulus_mesh)


mesh_from_file = meshio.read("mesh_annulus.msh")
triangle_mesh = create_mesh(mesh_from_file, "triangle", prune_z=True)
line_mesh = create_mesh(mesh_from_file, "line", prune_z=True)
meshio.write("mesh_2d.xdmf", triangle_mesh)
meshio.write("mesh_1d.xdmf", line_mesh)

The error reads:

triangle_mesh = create_mesh(mesh_from_file, "triangle", prune_z=True)
TypeError: create_mesh(): incompatible function arguments. The following argument types are supported:
    1. (arg0: dolfin.cpp.function.Function) -> dolfin.cpp.mesh.Mesh
    2. (arg0: object) -> dolfin.cpp.mesh.Mesh

Invoked with: <meshio mesh object>
  Number of points: 983
  Number of cells:
    triangle: 1798
  Point data: gmsh:dim_tags
  Cell data: gmsh:geometrical, 'triangle'; kwargs: prune_z=True

Any advice would be extremely helpful.

2

The dolfin package contains a function named create_mesh. The second time you use from dolfin import *, you overwrite the function create_mesh in your Python namespace. If you remove the second import, this error should be eliminated.

3

@conpierce8 , Now, I get the error:

triangle_mesh = create_mesh(mesh_from_file, "triangle", prune_z=True)
  File "fenicsgmshfenics.py", line 17, in create_mesh
    cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
  File "/home/WIN/avishmj/.local/lib/python3.8/site-packages/meshio/_mesh.py", line 187, in get_cell_data
    [d for c, d in zip(self.cells, self.cell_data[name]) if c.type == cell_type]
KeyError: 'gmsh:physical'
4

It appears that this new error occurs because you only write the mesh to the file “annulus_mesh.xdmf”, and not the MeshFunctions that mark the domains and boundaries. Therefore, when you you create “mesh_annulus.msh”, there are no physical groups stored in the file, and therefore you cannot access physical group information with “gmsh:physical”.

It doesn’t look like there’s a straightforward way to transfer your MeshFunctions mf and cf (defined on mesh) to the SubMesh annulus. But if you have created your mesh in Gmsh, maybe you could extract the submesh by first meshing the whole domain and exporting to a file, e.g. annulus_mesh.msh, then deleting the elements from the inner disk and exporting the remaining elements to a different file, e.g. submesh.msh?

5

That’s the problem. I need the mesh there to solve another variational form on the disc. If I delete it, I have to remesh it again. Then the problem will be that the boundary nodes of the disc and the inner boundary nodes of the annulus may not align. And I need them to be aligned.

6

If you use the Gmsh GUI, you can:

  1. Mesh the entire domain (disk + annulus)
  2. Export the mesh (e.g. “mesh.msh”)
  3. Click “Mesh > Delete > Surfaces” and select the disk. This will remove the mesh elements from the disk, but it doesn’t remove the underlying geometric surface. Furthermore, it will leave the mesh over the annulus intact.
  4. Export the mesh (e.g. “annulus_mesh.msh”).

The meshed annulus produced in this way will be conforming with the disk.

7

@conpierce8 ,

My actual problem consists of a for loop of length 100 and I need to do what you suggested for each iteration. I need FEniCS to communicate to gmsh inside the loop. Right now, I don’t know how to make FEniCS talk to gmsh.

8

You can use the following Python code, which creates the geometry in your .geo file above, and exports the entire mesh and only the annulus to two files. If you have your geometry defined in a .geo file already, you can simply use gmsh.open("file_name.geo") to load your geometry.

import gmsh
import numpy as np

# Always call gmsh.initialize() first
gmsh.initialize()

# Create geometry
gmsh.model.occ.addCircle(x=0, y=0, z=0, r=1, tag=1)
gmsh.model.occ.addCircle(x=0, y=0, z=0, r=2, tag=2)

gmsh.model.occ.addCurveLoop([1], tag=1)
gmsh.model.occ.addPlaneSurface([1], tag=1)

gmsh.model.occ.addCurveLoop([2], tag=2)
gmsh.model.occ.addCurveLoop([1], tag=3)
gmsh.model.occ.addPlaneSurface([2, 3], tag=2)

gmsh.model.occ.synchronize()

# Create physical groups
gmsh.model.addPhysicalGroup(dim=1, tags=[1], tag=1)
gmsh.model.addPhysicalGroup(dim=1, tags=[2], tag=2)
gmsh.model.addPhysicalGroup(dim=2, tags=[1], tag=3)
gmsh.model.addPhysicalGroup(dim=2, tags=[2], tag=4)

# Generate mesh
gmsh.model.mesh.generate(dim=2)

# Export full mesh
gmsh.write("entire_domain.msh")

# Clear mesh from the center disk and export annulus mesh
gmsh.model.mesh.clear(dimTags=[(2, 1)])
gmsh.write("annulus.msh")

# Always call gmsh.finalize() at the end
gmsh.finalize()
1 Like