# Mesh 2D gmsh to fenics legacy (2019.2)

**URL:** <https://fenicsproject.discourse.group/t/mesh-2d-gmsh-to-fenics-legacy-2019-2/15863>\
**Category:** mesh\
**Created:** [September 22, 2024, 6:56pm UTC](https://fenicsproject.discourse.group/t/mesh-2d-gmsh-to-fenics-legacy-2019-2/15863 "2024-09-22T18:56:01Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![Ceesar](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/ceesar/32/7506_2.png) [@Ceesar](https://fenicsproject.discourse.group/u/Ceesar)\
**Post date:** [September 22, 2024, 6:56pm UTC](https://fenicsproject.discourse.group/t/mesh-2d-gmsh-to-fenics-legacy-2019-2/15863/1 "2024-09-22T18:56:01Z")

</div>

Hi, I am generating my mesh and boundary conditions in GMSH. I use meshio to generate an XDMF file for use in legacy FEniCS code:

```auto
import gmsh
import meshio

# Função para inicializar o GMSH e configurar opções
def initialize_gmsh():
    gmsh.initialize()
    gmsh.option.setNumber("General.Terminal", 1) # Ativa saída no terminal
    gmsh.model.add("modelo_4_retangulos") # Adiciona o modelo

# Função para criar pontos e linhas para um retângulo
def create_rectangle(x_min, y_min, x_max, y_max):
    geom = gmsh.model.geo
    p1 = geom.addPoint(x_min, y_min, 0)
    p2 = geom.addPoint(x_max, y_min, 0)
    p3 = geom.addPoint(x_max, y_max, 0)
    p4 = geom.addPoint(x_min, y_max, 0)

    l1 = geom.addLine(p1, p2) # Linha inferior
    l2 = geom.addLine(p2, p3) # Linha direita
    l3 = geom.addLine(p3, p4) # Linha superior
    l4 = geom.addLine(p4, p1) # Linha esquerda

    loop = geom.addCurveLoop([l1, l2, l3, l4])
    return loop, [l1, l2, l3, l4] # Retorna também as linhas

# Função para configurar a malha transfinita
def set_transfinite_mesh(curves, nx, ny):
    for i, tag in enumerate(curves):
        gmsh.model.geo.mesh.setTransfiniteCurve(tag, (nx if i % 2 == 0 else ny) + 1)

# Função principal para gerar a malha com condições de contorno
def generate_mesh():
    geom = gmsh.model.geo

    # Dimensões dos retângulos
    small_width, small_height = 1, 1 # Dimensões dos retângulos pequenos
    large_width, large_height = 3, 2 # Dimensões do retângulo maior (domínio de projeto)

    # Cria os quatro retângulos
    loop1, curves1 = create_rectangle(0, 0, small_width, small_height) # Retângulo pequeno inferior esquerdo
    loop2, curves2 = create_rectangle(small_width, 0, small_width + large_width, large_height / 2) # Retângulo grande inferior
    loop3, curves3 = create_rectangle(small_width, large_height / 2, small_width + large_width, large_height) # Retângulo grande superior
    loop4, curves4 = create_rectangle(small_width + large_width, 0, small_width + large_width + small_width, small_height) # Retângulo pequeno inferior direito

    # Define a malha transfinita para cada conjunto de curvas
    set_transfinite_mesh(curves1, 10, 10) # Retângulo pequeno inferior esquerdo
    set_transfinite_mesh(curves2, 30, 10) # Retângulo grande inferior
    set_transfinite_mesh(curves3, 30, 10) # Retângulo grande superior
    set_transfinite_mesh(curves4, 10, 10) # Retângulo pequeno inferior direito

    # Adiciona condições de contorno (grupos físicos)
    gmsh.model.addPhysicalGroup(1, [curves1[3]], 201) # Inlet (vertical esquerda do primeiro retângulo)
    gmsh.model.setPhysicalName(1, 201, "Inlet")

    gmsh.model.addPhysicalGroup(1, [curves4[1]], 202) # Outlet (vertical direita do último retângulo)
    gmsh.model.setPhysicalName(1, 202, "Outlet")

    # Paredes: todas as demais arestas externas
    walls = [
        curves1[0], # Linha inferior do 1º retângulo 
        curves1[2], # Linha superior do 1º retângulo 

        curves2[0], # Linha superior do 2º retângulo

        curves3[1], # Linha direita do 3º retângulo
        curves3[2], # Linha superior do 3º retângulo
        curves3[3], # Linha esquerda do 3º retângulo

        curves4[0], # Linha inferior do 4º retângulo 
        curves4[2], # Linha superior do 4º retângulo 
    ]
    
    # Define o grupo de paredes
    gmsh.model.addPhysicalGroup(1, walls, 203) # Paredes externas
    gmsh.model.setPhysicalName(1, 203, "Walls")

    # Cria superfícies
    s1 = geom.addPlaneSurface([loop1])
    s2 = geom.addPlaneSurface([loop2])
    s3 = geom.addPlaneSurface([loop3])
    s4 = geom.addPlaneSurface([loop4])

    # Define o grupo físico para as superfícies (domínio)
    gmsh.model.addPhysicalGroup(2, [s1, s2, s3, s4], 101)
    gmsh.model.setPhysicalName(2, 101, "Domain")

    # Define malha estruturada nas superfícies
    geom.mesh.setTransfiniteSurface(s1)
    geom.mesh.setTransfiniteSurface(s2)
    geom.mesh.setTransfiniteSurface(s3)
    geom.mesh.setTransfiniteSurface(s4)

    # Sincroniza a geometria
    geom.synchronize()

    # Gera a malha
    gmsh.model.mesh.generate(2)

    # Salva a malha no formato .msh
    gmsh.write("heatsink.msh")

# Converte msh em um único xdmf
def convert_msh_to_xdmf(msh_filename, xdmf_filename):
    # Lê o arquivo .msh usando meshio
    mesh = meshio.read(msh_filename)

    # Filtra as células de domínio e fronteiras
    cells = [("triangle", mesh.cells_dict["triangle"]),
             ("line", mesh.cells_dict["line"])]

    # Armazena os dados do domínio e das fronteiras
    cell_data = {"gmsh:physical": [mesh.cell_data_dict["gmsh:physical"]["triangle"], 
                                   mesh.cell_data_dict["gmsh:physical"]["line"]]}

    # Cria a malha combinada
    combined_mesh = meshio.Mesh(
        points=mesh.points,
        cells=cells,
        cell_data=cell_data
    )

    # Salva a malha combinada no formato .xdmf
    meshio.write(xdmf_filename, combined_mesh)

# Função para finalizar o GMSH
def finalize_gmsh():
    gmsh.finalize()

# Executa o fluxo completo de GMSH e Meshio
def main():
    initialize_gmsh() # Inicializa GMSH
    generate_mesh() # Gera a malha com condições de contorno
    gmsh.fltk.run() # Abre a interface gráfica (opcional)
    convert_msh_to_xdmf("heatsink.msh", "heatsink.xdmf") # Converte a malha em um único arquivo .xdmf
    finalize_gmsh() # Finaliza GMSH

if __name__ == " __main__":
    main()

```

I have already tried several ways to import this mesh and boundaries, following responses in older threads, but I haven’t been successful, as some problem always comes up. Could someone help me by explaining the proper way to do it?

---

<div class="post-metadata">

**Author:** ![Ceesar](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/ceesar/32/7506_2.png) [@Ceesar](https://fenicsproject.discourse.group/u/Ceesar)\
**Post date:** [September 23, 2024, 12:42am UTC](https://fenicsproject.discourse.group/t/mesh-2d-gmsh-to-fenics-legacy-2019-2/15863/2 "2024-09-23T00:42:35Z")

</div>

I also tried using meshx to convert in a more practical way, following an old tutorial on another site, but again I had problems loading the file. Well, it’s a different case from the issue above, but either way, I’m informing you of my attempts ( [How to use MeshX to convert the mesh file from Gmsh to XDMF file in FEniCS – Computational Mechanics](https://computationalmechanics.in/how-to-use-meshx-to-convert-the-mesh-file-from-gmsh-to-xdmf-file-in-fenics/))

**GMSH:**

```auto
Point(1) = {0, 0, 0, 1};
Point(2) = {1, 0, 0, 1};
Point(3) = {1, 1, 0, 1};
Point(4) = {0, 1, 0, 1};
Point(5) = {1, 0.75, 0, 1};
Point(6) = {1, 0.5, 0, 1};
Point(7) = {0.5, 0.5, 0, 1};
Point(8) = {0.5, 0.75, 0, 1};
​
Line(1) = {1, 2};
Line(2) = {2, 6};
Line(3) = {6, 7};
Line(4) = {7, 8};
Line(5) = {8, 5};
Line(6) = {5, 6};
Line(7) = {5, 3};
Line(8) = {3, 4};
Line(9) = {4, 1};
​
Curve Loop(1) = {8, 9, 1, 2, 3, 4, 5, 7};
Plane Surface(1) = {1};
Curve Loop(2) = {5, 6, 3, 4};
Plane Surface(2) = {2};
​
Physical Curve("top", 10) = {8};
Physical Curve("left", 11) = {9};
Physical Curve("right", 12) = {7, 6, 2};
Physical Curve("bottom", 13) = {1};
​
Physical Surface("Domain", 14) = {1};
Physical Surface("Obstacle", 15) = {2};

```

**Meshx:**

`meshx file.msh`

**FEniCS:**

```auto
# Define mesh
mesh = Mesh()
with XDMFFile("mesh/triangle.xdmf") as infile:
    infile.read(mesh)
​
# mesh value collection.
mvc_1 = MeshValueCollection("size_t", mesh, 1) #for curve with dim 1
mvc_2 = MeshValueCollection("size_t", mesh, 2) #for surface with dim 2
​
#import the json file so that tag names can be used
f = open('mesh/tags.json') 
tags = json.load(f)
​
with XDMFFile("mesh/line.xdmf") as infile:
    print("Reading 1d line data into dolfin mvc")
    infile.read(mvc_1, "tag")
​
print("Constructing MeshFunction from MeshValueCollection")
boundaries = cpp.mesh.MeshFunctionSizet(mesh, mvc_1) # creating mesh function for curve 
​
with XDMFFile("mesh/triangle.xdmf") as infile:
    print("Reading 2d surface data into dolfin mvc")
    infile.read(mvc_2, "tag")
​
print("Constructing MeshFunction from MeshValueCollection")
domains = cpp.mesh.MeshFunctionSizet(mesh, mvc_2) # creating mesh function for surface 
​
#--------------------------------------------------------------------------------------------------------------------------------
​
# Define Dirichlet boundary conditions at top and bottom boundaries
bcs = [DirichletBC(V, 5.0, boundaries, tags['Top']),
       DirichletBC(V, 0.0, boundaries, tags['Bottom'])]
​
#--------------------------------------------------------------------------------------------------------------------------------
​
# Define new measures associated with the interior domains and
# exterior boundaries
dx = Measure("dx")(subdomain_data=domains)
ds = Measure("ds")(subdomain_data=boundaries)
​
#--------------------------------------------------------------------------------------------------------------------------------
​
# Define variational form
F = (inner(a0*grad(u), grad(v))*dx(tags['Domain']) + inner(a1*grad(u), grad(v))*dx(tags['Obstacle'])
     - g_L*v*ds(tags['Left']) - g_R*v*ds(tags['Right'])
     - f*v*dx(tags['Domain']) - f*v*dx(tags['Obstacle']))

```

**Error:**

```auto
String tag Num Tag Dim
------------ --------- -----
top 10 1
left 11 1
right 12 1
bottom 13 1
Domain 14 2
Obstacle 15 2
Creating gmsh:physical mesh
There was some error 😱
Traceback (most recent call last):
  File "/usr/local/bin/meshx", line 8, in <module>
    sys.exit(main())
  File "/usr/local/lib/python3.10/dist-packages/meshx/_cli.py", line 13, in main
    write(args)
  File "/usr/local/lib/python3.10/dist-packages/meshx/_helpers.py", line 13, in write
    converter.write_sub_domains()
  File "/usr/local/lib/python3.10/dist-packages/meshx/_helpers.py", line 86, in write_sub_domains
    infile.read(mesh)
RuntimeError: 

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
*** https://fenicsproject.discourse.group/
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error: Unable to open XDMF file.
*** Reason: XDMF file "mesh/line.xdmf" does not exist.
*** Where: This error was encountered inside XDMFFile.cpp.
*** Process: 0
*** 
*** DOLFIN version: 2019.2.0.13.dev0
*** Git changeset: ubuntu
*** -------------------------------------------------------------------------

```

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [September 23, 2024, 6:18am UTC](https://fenicsproject.discourse.group/t/mesh-2d-gmsh-to-fenics-legacy-2019-2/15863/3 "2024-09-23T06:18:49Z")

</div>

A general advice is to set physical group after all synchronisation is done.

Please also note that you could have presented a way smaller example.

Also note that you cannot write both triangles and lines to the same xdmf file and be compatible with legacy dolfin, this is for instance shown here:

> [@Error defining surface and measures for a 3D mesh imported from Gmsh](https://fenicsproject.discourse.group/t/error-defining-surface-and-measures-for-a-3d-mesh-imported-from-gmsh/14056/4):
>
> Edge integrals are not supported in 3D meshes without using MeshView. You have corrected the integral measure for the surfaces (ds, not dS) Consider the following: from dolfin import \* import meshio 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]})…

And

> **[Mesh generation and conversion with GMSH and PYGMSH](https://jsdokken.com/src/pygmsh_tutorial.html#second)**
>
> Personal webpage

---

<div class="post-metadata">

**Author:** ![Ceesar](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/ceesar/32/7506_2.png) [@Ceesar](https://fenicsproject.discourse.group/u/Ceesar)\
**Post date:** [September 27, 2024, 7:01pm UTC](https://fenicsproject.discourse.group/t/mesh-2d-gmsh-to-fenics-legacy-2019-2/15863/4 "2024-09-27T19:01:56Z")

</div>

Thank you for the response, I am now able to import XDMF meshes into FEniCS Legacy. I am basing my work on [Mesh generation and conversion with GMSH and PYGMSH | Jørgen S. Dokken](https://jsdokken.com/src/pygmsh_tutorial.html#second) (DFG-2D 2 Turek benchmark).

I made modifications for my problem, essentially a T-shape geometry (the central rectangle is a design subdomain for topological optimization). However, there seems to be some issue with the boundary conditions, so I returned to the DFG-2D 2 Turek benchmark problem. The mesh below is the same as in the tutorial (just to provide an example with fewer lines of code):

```auto
import meshio
import gmsh
import pygmsh

resolution = 0.01
# Channel parameters
L = 2.2
H = 0.41
c = [0.2, 0.2, 0]
r = 0.05

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

# Add points with finer resolution on left side
points = [
    model.add_point((0, 0, 0), mesh_size=resolution),
    model.add_point((L, 0, 0), mesh_size=5 * resolution),
    model.add_point((L, H, 0), mesh_size=5 * resolution),
    model.add_point((0, H, 0), mesh_size=resolution),
]

# Add lines between all points creating the rectangle
channel_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
channel_loop = model.add_curve_loop(channel_lines)
plane_surface = model.add_plane_surface(channel_loop, holes=[circle.curve_loop])

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

volume_marker = 6
model.add_physical([plane_surface], "Volume")
model.add_physical([channel_lines[0]], "Inflow")
model.add_physical([channel_lines[2]], "Outflow")
model.add_physical([channel_lines[1], channel_lines[3]], "Walls")
model.add_physical(circle.curve_loop.curves, "Obstacle")

geometry.generate_mesh(dim=2)
gmsh.write("mesh.msh")
gmsh.clear()
geometry. __exit__ ()

mesh_from_file = meshio.read("mesh.msh")

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

line_mesh = create_mesh(mesh_from_file, "line", prune_z=True)
meshio.write("facet_mesh.xdmf", line_mesh)

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

```

As I noticed issues with the results of my problem, I created a simple Stokes case just to evaluate the behavior of the Inlet. It’s as if the condition were applied only to the extreme points of the Inlet line. I have already checked other codes, but I couldn’t identify what the problem might be:

```auto
from dolfin import *

# Function to read the generated mesh and boundary data
def ler_malha():
    mesh = Mesh()
    xdmf = XDMFFile(mesh.mpi_comm(), "mesh.xdmf")
    xdmf.read(mesh)

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

    print("Boundary markers present:", set(cf.array()))

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

    print("Boundary markers present:", set(sf.array()))

    return mesh, cf, sf

def main():
    mesh, cf, sf = ler_malha() # Read the mesh and mesh functions

    V_h = VectorElement("CG", mesh.ufl_cell(), 2)
    Q_h = FiniteElement("CG", mesh.ufl_cell(), 1)
    W = FunctionSpace(mesh, V_h * Q_h)
    V, Q = W.split()

    v, q = TestFunctions(W)
    x = TrialFunction(W)
    u, p = split(x)
    s = Function(W, name="State")
    V_collapse = V.collapse()
    g = Function(V_collapse, name="Control")

    nu = Constant(1)

    # Define the expression for parabolic velocity
    inflow_velocity = Expression(("1", "0"), degree=1)

    # Boundary conditions
    bc_inlet = DirichletBC(W.sub(0), inflow_velocity, sf, 2) # Inlet marker
    bc_walls = DirichletBC(W.sub(0), Constant((0, 0)), sf, 4) # Walls marker
    bc_outlet = DirichletBC(W.sub(1), Constant(0), sf, 3) # Outlet marker (pressure)

    bcs = [bc_inlet, bc_walls, bc_outlet]

    # Weak formulation of Stokes
    a = (nu * inner(grad(u), grad(v)) * dx
        - inner(p, div(v)) * dx
        - inner(q, div(u)) * dx)
    L = inner(Constant((0, 0)), v) * dx

    # Solve the problem
    A, b = assemble_system(a, L, bcs)
    solve(A, s.vector(), b)

    # Extract the solution
    u, p = s.split()

    # Save results to files
    file_velocity = File("resultado_velocidade.pvd")
    file_pressure = File("resultado_pressao.pvd")
    file_velocity << u
    file_pressure << p

if __name__ == " __main__":
    main()

```

 ![image](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/f/fec03d921dbc69bfb247dc2f0f4d4eb60af360b7.png)