Issue with Importing Physical Groups into FEniCS

[Sander_Buyse]

Hi,

I’m working with FEniCS and meshio to load a mesh from a GMSH file. I have both physical subdomains and boundaries defined in the mesh, but I’m encountering an issue with getting these physical groups into FEniCS properly. I need to use these physical groups for applying boundary conditions and material properties, but they aren’t showing up as expected in the simulation.

Here’s my approach:
import meshio
from fenics import *
import p1_simulation_config

# ────────── Mesh Loading ──────────
def read_and_load_mesh(mesh_filename):
    # Read and Convert Mesh
    msh = meshio.read(mesh_filename)

    # Extract elements and their physical regions
    triangle_cells = msh.cells_dict["triangle"]
    physical_regions = msh.cell_data_dict["gmsh:physical"]["triangle"]

    line_cells = msh.cells_dict["line"]
    physical_boundaries = msh.cell_data_dict["gmsh:physical"]["line"]

    # Convert to XDMF for FEniCS compatibility
    meshio.write("mesh.xdmf", meshio.Mesh(points=msh.points, cells={"triangle": triangle_cells}))
    meshio.write("boundaries.xdmf", meshio.Mesh(points=msh.points, cells={"line": line_cells}, 
                                              cell_data={"line": [physical_boundaries]}))
    meshio.write("subdomains.xdmf", meshio.Mesh(points=msh.points, cells={"triangle": triangle_cells}, 
                                              cell_data={"name": [physical_regions]}))

    # Load Mesh into FEniCS
    mesh = Mesh()
    with XDMFFile("mesh.xdmf") as infile:
        infile.read(mesh)

    boundaries = MeshFunction("size_t", mesh, mesh.topology().dim() - 1, mesh.domains())
    with XDMFFile("boundaries.xdmf") as infile:
        infile.read(boundaries)

    subdomains = MeshFunction("size_t", mesh, mesh.topology().dim())
    with XDMFFile("subdomains.xdmf") as infile:
        infile.read(subdomains)

    return mesh, boundaries, subdomains

The issue is that when I try to load the physical groups (boundaries and subdomains), they don’t seem to be correctly read by FEniCS. Could someone help me understand why the physical groups aren’t being imported correctly and how to resolve this issue? Specifically, I need guidance on:

Ensuring that the physical boundaries and subdomains are properly recognized by FEniCS.
Correctly referencing these physical groups when applying boundary conditions and material properties in the simulation.

Any help would be greatly appreciated!

[dokken]

As you have not provided the mesh, it is hard to give any further guidance.
You have also not provided any error message.

See for instance the recent post:

Fail to read tagged facets from mesh generated with GMSH

You are here trying to integrate over the vertices of your mesh? Then you are not integrating over a facet (“ds”) but a vertex (“dp”). However, lets take 10 steps backwards, and consider: import meshio import gmsh gmsh.initialize() gmsh.model.add("my model") c_r = [0, 0, 0] c_R = [0, 0, 0] r = 1 R = 2 resolution = 0.1 print(f"Mesh resolution = {resolution}") disk_r = gmsh.model.occ.addDisk(c_r[0], c_r[1], c_r[2], r, r) disk_R = gmsh.model.occ.addDisk(c_R[0], c_R[1], c_R[2], R, R) ring = gm…

for some guidance.

[Sander_Buyse]

Sorry, I shal give both. Thank you for the help!

Code mesh_generating:

import gmsh
import p1_simulation_config




# ────────── Mesh Generation ──────────
def create_mesh(visualize = False, mesh_filename = "simulation_mesh.msh"):
    gmsh.initialize()
    gmsh.model.add("simulation_mesh")

    # Calculate the total height of the domain (glass + air + glass)
    domain_size = p1_simulation_config.h_glass_1 + p1_simulation_config.h_air + p1_simulation_config.h_glass_2

    # Define characteristic length based on the domain size and the refinement factor
    l_c = domain_size / p1_simulation_config.N_refinement       # Characteristic length for mesh resolution


    # ────────── Define Geometry Points ──────────
    # Define the points for the simulation domain (glass layers, air gap, and conductors)
    p_1 = gmsh.model.geo.addPoint(0, 0, 0, l_c)
    p_2 = gmsh.model.geo.addPoint(p1_simulation_config.w_simulation, 0, 0, l_c)

    p_3 = gmsh.model.geo.addPoint(0, p1_simulation_config.h_glass_1, 0, l_c)
    p_4 = gmsh.model.geo.addPoint(p1_simulation_config.w_simulation, p1_simulation_config.h_glass_1, 0, l_c)
    
    p_5 = gmsh.model.geo.addPoint(0, p1_simulation_config.h_glass_1 + p1_simulation_config.h_air, 0, l_c)
    p_6 = gmsh.model.geo.addPoint(p1_simulation_config.w_simulation, p1_simulation_config.h_glass_1 + p1_simulation_config.h_air, 0, l_c)
    
    p_7 = gmsh.model.geo.addPoint(0, p1_simulation_config.h_glass_1 + p1_simulation_config.h_air + p1_simulation_config.h_glass_2, 0, l_c)
    p_8 = gmsh.model.geo.addPoint(p1_simulation_config.w_simulation, p1_simulation_config.h_glass_1 + p1_simulation_config.h_air + p1_simulation_config.h_glass_2, 0, l_c)
    
    p_9 = gmsh.model.geo.addPoint(p1_simulation_config.x_conductor_1, p1_simulation_config.h_glass_1, 0, l_c)
    p_10 = gmsh.model.geo.addPoint(p1_simulation_config.x_conductor_1 + p1_simulation_config.w_conductor_1, p1_simulation_config.h_glass_1, 0, l_c)
    
    p_11 = gmsh.model.geo.addPoint(p1_simulation_config.x_conductor_2, p1_simulation_config.h_glass_1, 0, l_c)
    p_12 = gmsh.model.geo.addPoint(p1_simulation_config.x_conductor_2 + p1_simulation_config.w_conductor_2, p1_simulation_config.h_glass_1, 0, l_c)

    p_13 = gmsh.model.geo.addPoint(p1_simulation_config.x_conductor_1, p1_simulation_config.h_glass_1 + p1_simulation_config.h_conductor, 0, l_c)
    p_14 = gmsh.model.geo.addPoint(p1_simulation_config.x_conductor_1 + p1_simulation_config.w_conductor_1, p1_simulation_config.h_glass_1 + p1_simulation_config.h_conductor, 0, l_c)
    
    p_15 = gmsh.model.geo.addPoint(p1_simulation_config.x_conductor_2, p1_simulation_config.h_glass_1 + p1_simulation_config.h_conductor, 0, l_c)
    p_16 = gmsh.model.geo.addPoint(p1_simulation_config.x_conductor_2 + p1_simulation_config.w_conductor_2, p1_simulation_config.h_glass_1 + p1_simulation_config.h_conductor, 0, l_c)
    

    # ────────── Define Geometry Lines ──────────
    # Define the lines that form the boundaries and connections between the domain points
    l_1 = gmsh.model.geo.addLine(p_1, p_3)
    l_2 = gmsh.model.geo.addLine(p_3, p_5)
    l_3 = gmsh.model.geo.addLine(p_5, p_7)
    l_4 = gmsh.model.geo.addLine(p_7, p_8)
    l_5 = gmsh.model.geo.addLine(p_6, p_8)
    l_6 = gmsh.model.geo.addLine(p_4, p_6)
    l_7 = gmsh.model.geo.addLine(p_2, p_4)
    l_8 = gmsh.model.geo.addLine(p_1, p_2)

    l_9 = gmsh.model.geo.addLine(p_3, p_9)
    l_10 = gmsh.model.geo.addLine(p_9, p_10)
    l_11 = gmsh.model.geo.addLine(p_10, p_11)
    l_12 = gmsh.model.geo.addLine(p_11, p_12)
    l_13 = gmsh.model.geo.addLine(p_4, p_12)
    l_14 = gmsh.model.geo.addLine(p_5, p_6)

    l_15 = gmsh.model.geo.addLine(p_9, p_13)
    l_16 = gmsh.model.geo.addLine(p_13, p_14)
    l_17 = gmsh.model.geo.addLine(p_10, p_14)
    l_18 = gmsh.model.geo.addLine(p_11, p_15)
    l_19 = gmsh.model.geo.addLine(p_15, p_16)
    l_20 = gmsh.model.geo.addLine(p_12, p_16)


    # ────────── Surface Definitions ──────────
    # Define curve loops and corresponding surfaces for each region (glass layers and air gap)
    loop_glass_1 = gmsh.model.geo.addCurveLoop([l_1, l_9, l_10, l_11, l_12, -l_13, -l_7, -l_8])
    surface_glass_1 = gmsh.model.geo.addPlaneSurface([loop_glass_1])

    loop_air = gmsh.model.geo.addCurveLoop([l_2, l_14, -l_6, l_13, l_20, -l_19, -l_18, -l_11, l_17, -l_16, -l_15, -l_9])
    surface_air = gmsh.model.geo.addPlaneSurface([loop_air])

    loop_glass_2 = gmsh.model.geo.addCurveLoop([l_3, l_4, -l_5, -l_14])
    surface_glass_2 = gmsh.model.geo.addPlaneSurface([loop_glass_2])


    # ────────── Assign Physical Groups ──────────
    # Define physical groups for material regions (surfaces)
    glass_1_physical = gmsh.model.addPhysicalGroup(2, [surface_glass_1])
    gmsh.model.setPhysicalName(2, glass_1_physical, "Glass_layer_1")

    air_physical = gmsh.model.addPhysicalGroup(2, [surface_air])
    gmsh.model.setPhysicalName(2, air_physical, "Air_layer")

    glass_2_physical = gmsh.model.addPhysicalGroup(2, [surface_glass_2])
    gmsh.model.setPhysicalName(2, glass_2_physical, "Glass_layer_2")

    # Define physical groups for conductor boundaries
    conductor_1_physical = gmsh.model.addPhysicalGroup(1, [l_10, l_15, l_16, l_17])
    gmsh.model.setPhysicalName(1, conductor_1_physical, "Conductor_1")

    conductor_2_physical = gmsh.model.addPhysicalGroup(1, [l_12, l_18, l_19, l_20])
    gmsh.model.setPhysicalName(1, conductor_2_physical, "Conductor_2")

    # Define physical groups for domain boundaries
    lower_boundary_physical = gmsh.model.addPhysicalGroup(1, [l_8])
    gmsh.model.setPhysicalName(1, lower_boundary_physical, "Lower_boundary")

    upper_boundary_physical = gmsh.model.addPhysicalGroup(1, [l_4])
    gmsh.model.setPhysicalName(1, upper_boundary_physical, "Upper_boundary")

    left_boundary_physical = gmsh.model.addPhysicalGroup(1, [l_1, l_2, l_3])
    gmsh.model.setPhysicalName(1, left_boundary_physical, "Left_boundary")

    right_boundary_physical = gmsh.model.addPhysicalGroup(1, [l_5, l_6, l_7])
    gmsh.model.setPhysicalName(1, right_boundary_physical, "Right_boundary")

    # Define physical groups for glass-air interface boundaries
    interface_boundary_physical = gmsh.model.addPhysicalGroup(1, [l_9, l_11, l_13, l_14])
    gmsh.model.setPhysicalName(1, interface_boundary_physical, "Glass-air interface boundaries")


    # ────────── Finalize Geometry and Generate Mesh ──────────
    gmsh.model.geo.synchronize()

    # ── TODO: Improve mesh refinement by adjusting characteristic length (l_c) for different regions. 
    # Specifically, refine mesh around conductors or other regions of interest for better accuracy.
    gmsh.model.mesh.generate(2)


    # ────────── Visualization ──────────
    if visualize:
        gmsh.fltk.run()


    # ────────── Finalize and Save Mesh ──────────
    # Save the mesh to a file in GMSH format
    gmsh.write(mesh_filename)
    gmsh.finalize()

    # Return the filename for further use
    return mesh_filename




if __name__ == "__main__":
    create_mesh()

Error message:
File “/home/sande/EP&I - Liquid Crystal Simulation/Simulation - basic air slab/p3_electrostatic_solver.py”, line 39, in read_and_load_mesh
infile.read(boundaries, “line”)
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

---

*** fenics-support@googlegroups.com

---

*** Remember to include the error message listed below and, if possible,
*** include a minimal running example to reproduce the error.

---

*** -------------------------------------------------------------------------
*** Error: Unable to open MeshFunction for reading.
*** Reason: Mesh Grid with data Attribute not found in XDMF.
*** Where: This error was encountered inside XDMFFile.cpp.
*** Process: 0

---

*** DOLFIN version: 2019.1.0
*** Git changeset: 997d6afd946dfa52f86ee2737314d99cdae5383f
*** -------------------------------------------------------------------------

[dokken]

@Sander_Buyse

    boundaries = MeshFunction("size_t", mesh, mesh.topology().dim() - 1, mesh.domains())
    with XDMFFile("boundaries.xdmf") as infile:
        infile.read(boundaries)

    subdomains = MeshFunction("size_t", mesh, mesh.topology().dim())
    with XDMFFile("subdomains.xdmf") as infile:
        infile.read(subdomains)

This is not the right way of reading in data from the meshio files. As shown in:

@dokken

See for instance the recent post:

You are here trying to integrate over the vertices of your mesh? Then you are not integrating over a facet (“ds”) but a vertex (“dp”). However, lets take 10 steps backwards, and consider: import meshio import gmsh gmsh.initialize() gmsh.model.add(“my model”) c_r = [0, 0, 0] c_R = [0, 0, 0] r = 1 R = 2 resolution = 0.1 print(f"Mesh resolution = {resolution}") disk_r = gmsh.model.occ.addDisk(c_r[0], c_r[1], c_r[2], r, r) disk_R = gmsh.model.occ.addDisk(c_R[0], c_R[1], c_R[2], R, R) ring = gm

One read in such data with

dim()}")

# read the triangles
mvc = MeshValueCollection( "size_t", mesh, mesh.topology().dim() )
with XDMFFile( "triangle_mesh.xdmf" ) as infile:
    infile.read( mvc, "name_to_read" )
vf = cpp.mesh.MeshFunctionSizet( mesh, mvc )
xdmf.close()

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

#read the vertices
mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim()-2)
with XDMFFile("vertex_mesh.xdmf") as infile:
    infile.read(mvc, "name_to_read")
pf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
xdmf.close()