Mesh 2D gmsh to fenics legacy (2019.2)

Ceesar · September 22, 2024, 6:56pm

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:

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?

Ceesar · September 23, 2024, 12:42am

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)

GMSH:

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:

# 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:

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
*** -------------------------------------------------------------------------

dokken · September 23, 2024, 6:18am

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:

And

Ceesar · September 27, 2024, 7:01pm

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 (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):

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:

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()

Topic		Replies	Views
Create a mesh with Gmsh General	2	297	September 20, 2023
Problem converting gmsh to xdmf for fenics mesh	1	514	December 1, 2022
Gmsh and adaptation	4	618	April 27, 2022
What is the right way to create 3D meshes for FEniCS with Gmsh? mesh	1	1493	October 19, 2021
Import from gmsh and define boundary conditions	1	1574	March 28, 2020

Mesh 2D gmsh to fenics legacy (2019.2)

Related topics