#!/usr/bin/env python
# coding: utf-8
# # 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)}$.
# In[1]:
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.
# In[2]:
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.
# In[3]:
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`).
# In[4]:
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.
# In[5]:
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.
# In[6]:
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)
# 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.
# In[7]:
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](https://jsdokken.com/converted_files/tutorial_pygmsh.html) to split the `meshio` mesh into triangle and line meshes containing the corresponding markers.
# In[8]:
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)