This work is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License

The corresponding files can be obtained from:

• Jupyter Notebook: I_beam_gmsh.ipynb
• Python script: I_beam_gmsh.py

Generating a shell model for a shallow arch I-beam with the Gmsh Python API

This notebook illustrates the use of the Gmsh Python API for generating a shell model of a curved I-beam model.

Geometric parameters

We consider a I-shaped constant cross-section for the beam with bottom and top flange widths wb and wt and of height h. The arch will be a portion of a circular arc of total angle opening \(2\theta\) and chord length \(L\). The circular arc radius is therefore given by \(R=\dfrac{L}{2\sin(\theta)}\) and the arch rise by \(f=R-\dfrac{L}{2\tan(\theta)}\).

import gmsh
import numpy as np

filename = "I_beam"
# I-beam profile
wb = 0.2                # bottom flange width
wt = 0.3                # top flange width
h = 0.5                 # web height
# Arch geometry
theta = np.pi/6         # arch opening half-angle
L = 10.                  # chord length
R = L/2/np.sin(theta)   # arch radius
f = R-L/2/np.tan(theta) # arch rise

Defining points and lines for the left-hand side cross-section

We first begin by initializing Gmsh ang use the built-in geo CAD kernel since we won’t require any boolean operation which require the use of the Open Cascade kernel.

gmsh.initialize()
geom = gmsh.model.geo

Next, we define the points and lines of the left-hand side cross-section. To facilitate their definition, these objects are first defined for a vertical cross-section located in the \(X=0\) plane and centered on the origin \((0,0,0)\). We split them into a group of bottom_points for the bottom_flange and a group of top_points for the top_flange. The web then connects the middle points of each group. The start_section contains the web, bottom and top flange lines.

Note that Gmsh requires passing a mesh density for each point definition. However, this lcar value will not be used in the end since we will later define the number of points on each line using the TransfiniteCurve command.

lcar = 0.1 # characteristic mesh size density (will not be used)
bottom_points = [geom.addPoint(0, -wb/2, -h/2, lcar),
                 geom.addPoint(0, 0, -h/2, lcar),
                 geom.addPoint(0, wb/2, -h/2, lcar)]
top_points = [geom.addPoint(0, -wt/2, h/2, lcar),
              geom.addPoint(0, 0, h/2, lcar),
              geom.addPoint(0, wt/2, h/2, lcar)]
bottom_flange = [geom.addLine(bottom_points[0], bottom_points[1]),
                 geom.addLine(bottom_points[1], bottom_points[2])]
web = [geom.addLine(bottom_points[1], top_points[1])]
top_flange = [geom.addLine(top_points[0], top_points[1]),
              geom.addLine(top_points[1], top_points[2])]
start_section = bottom_flange + web + top_flange

Cross-section rotation and extrusion

Now, we rotate this section by an angle equal to \(-\theta\) around the \(Y\) axis to obtain our initial starting section for the arch. The rotate command requires a dimTags first argument i.e. a list of tuple of the form (dim, label) where dim is the entity dimension (dim=1 here for lines) and label is the label entity. The following arguments are then x, y, z, ax, ay, az, angle defining a rotation of angle angle around an axis of revolution passing through a point (x, y, z) along a direction (ax, ay, az).

dimTags = [(1, l) for l in start_section]
geom.rotate(dimTags, 0, 0, 0, 0, 1, 0, -theta)

We will now extrude the previously defined cross-section along a circular arc. The extrude command is used for extruding along a line, here we will need the ̀revolve command to extrude the cross-section by following a rotation of angle \(2\theta\) which is defined similarly to the rotate command. Here, the axis follows the direction \((0, 1, 0)\) and passes through the point \((L/2, 0, -(R-f))\). We also pass a list of layers and corresponding numbers of elements generated during extrusion. Here, we have only one layer (normalized height = 1.0) of 50 elements.

The revolve command returns a list of dimTags corresponding to the newly created entities. For instance, when extruding a line, this will create a surface and three new lines. The output stores these entities as follows: \([\text{end line}, \text{ surface}, \text{ lateral line }1,\text{ lateral line }2]\). We loop on the various entities composing the starting cross-section and append the newly created end-lines and surfaces in corresponding lists. Note that we need to call synchronize to update the corresonding Gmsh model with the newly created entities.

end_bottom_flange = []
end_top_flange = []
end_web = []
surfaces = []
for ini, end in zip([bottom_flange, web, top_flange],
                    [end_bottom_flange, end_web, end_top_flange]):
    for l in ini:
        outDimTags = geom.revolve([(1, l)], L/2, 0, -(R-f), 0, 1, 0, 2*theta,
                                 numElements=[50], heights=[1.0])
        end.append(outDimTags[0][1])
        surfaces.append(outDimTags[1][1])
        geom.synchronize()
end_section = end_bottom_flange + end_web + end_top_flange

Tagging and mesh generation

We finish by specifing the number of elements for the flange and web discretization using the setTransfiniteCurve command. We also affect the physical tag 1 for the left-hand side cross-section and 2 for the right-hand side cross section. We do not forget to also add physical tags to the surfaces otherwise they well later be ignored when generating the mesh.

for f in bottom_flange + top_flange + end_bottom_flange + end_top_flange:
    geom.mesh.setTransfiniteCurve(f, 6)
for w in web + end_web:
    geom.mesh.setTransfiniteCurve(w, 11)

gmsh.model.addPhysicalGroup(2, surfaces, 1)
gmsh.model.addPhysicalGroup(1, start_section, 1)
gmsh.model.addPhysicalGroup(1, end_section, 2)

2

We can now call the mesh generate function (up to dimension dim=2 since we don’t have any 3D solid element here). The generated mesh is then exported to a .msh which we will now read and export to XDMF format. We finish by the finalize command since we are finished with the Gmsh API.

Note that you can export the corresponding geo file by writing to a file with a .geo_unrolled extension.

gmsh.model.mesh.generate(dim=2)
gmsh.write(filename + ".msh")

gmsh.finalize()

Reading with meshio and export to XDMF format

Now that we have a .msh file, we can read that back in with the meshio package. We will obtain a mesh containing triangular cells as well as line cells which will contain the boundary markers. We use the create_mesh function from Jørgen S. Dokken pygmsh tutorial to split the meshio mesh into triangle and line meshes containing the corresponding markers.

import meshio

# from J. Dokken website
def create_mesh(mesh, cell_type):
    cells = np.vstack([cell.data for cell in mesh.cells if cell.type==cell_type])
    data = np.hstack([mesh.cell_data_dict["gmsh:physical"][key]
                           for key in mesh.cell_data_dict["gmsh:physical"].keys() if key==cell_type])
    mesh = meshio.Mesh(points=mesh.points, cells={cell_type: cells},
                       cell_data={"name_to_read": [data]})
    return mesh

mesh = meshio.read(filename + ".msh")

shell_mesh = create_mesh(mesh, "triangle")
line_mesh = create_mesh(mesh, "line")
meshio.write(filename + ".xdmf", shell_mesh)
meshio.write(filename + "_facet_region.xdmf", line_mesh)