{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Generating a shell model for a shallow arch I-beam with the `Gmsh` Python API\n", "\n", "This notebook illustrates the use of the `Gmsh` Python API for generating a shell model of a curved I-beam model.\n", "\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "## Geometric parameters\n", "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)}$." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import gmsh\n", "import numpy as np\n", "\n", "filename = \"I_beam\"\n", "# I-beam profile \n", "wb = 0.2 # bottom flange width\n", "wt = 0.3 # top flange width\n", "h = 0.5 # web height\n", "# Arch geometry\n", "theta = np.pi/6 # arch opening half-angle\n", "L = 10. # chord length\n", "R = L/2/np.sin(theta) # arch radius\n", "f = R-L/2/np.tan(theta) # arch rise" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Defining points and lines for the left-hand side cross-section\n", "\n", "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." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "gmsh.initialize() \n", "geom = gmsh.model.geo" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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. \n", "\n", "> 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." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "lcar = 0.1 # characteristic mesh size density (will not be used)\n", "bottom_points = [geom.addPoint(0, -wb/2, -h/2, lcar),\n", " geom.addPoint(0, 0, -h/2, lcar),\n", " geom.addPoint(0, wb/2, -h/2, lcar)]\n", "top_points = [geom.addPoint(0, -wt/2, h/2, lcar),\n", " geom.addPoint(0, 0, h/2, lcar),\n", " geom.addPoint(0, wt/2, h/2, lcar)]\n", "bottom_flange = [geom.addLine(bottom_points[0], bottom_points[1]),\n", " geom.addLine(bottom_points[1], bottom_points[2])]\n", "web = [geom.addLine(bottom_points[1], top_points[1])]\n", "top_flange = [geom.addLine(top_points[0], top_points[1]),\n", " geom.addLine(top_points[1], top_points[2])]\n", "start_section = bottom_flange + web + top_flange" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Cross-section rotation and extrusion\n", "\n", "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`)." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "dimTags = [(1, l) for l in start_section]\n", "geom.rotate(dimTags, 0, 0, 0, 0, 1, 0, -theta)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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.\n", "\n", "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:\n", "$[\\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." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "end_bottom_flange = []\n", "end_top_flange = []\n", "end_web = []\n", "surfaces = []\n", "for ini, end in zip([bottom_flange, web, top_flange],\n", " [end_bottom_flange, end_web, end_top_flange]):\n", " for l in ini:\n", " outDimTags = geom.revolve([(1, l)], L/2, 0, -(R-f), 0, 1, 0, 2*theta,\n", " numElements=[50], heights=[1.0])\n", " end.append(outDimTags[0][1])\n", " surfaces.append(outDimTags[1][1])\n", " geom.synchronize()\n", "end_section = end_bottom_flange + end_web + end_top_flange" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Tagging and mesh generation\n", "\n", "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." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "2" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "for f in bottom_flange + top_flange + end_bottom_flange + end_top_flange:\n", " geom.mesh.setTransfiniteCurve(f, 6)\n", "for w in web + end_web:\n", " geom.mesh.setTransfiniteCurve(w, 11)\n", "\n", "gmsh.model.addPhysicalGroup(2, surfaces, 1)\n", "gmsh.model.addPhysicalGroup(1, start_section, 1)\n", "gmsh.model.addPhysicalGroup(1, end_section, 2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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.\n", "\n", "> Note that you can export the corresponding `geo` file by writing to a file with a `.geo_unrolled` extension." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "gmsh.model.mesh.generate(dim=2)\n", "gmsh.write(filename + \".msh\")\n", "\n", "gmsh.finalize()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Reading with `meshio` and export to XDMF format\n", "\n", "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." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "import meshio\n", "\n", "# from J. Dokken website\n", "def create_mesh(mesh, cell_type):\n", " cells = np.vstack([cell.data for cell in mesh.cells if cell.type==cell_type])\n", " data = np.hstack([mesh.cell_data_dict[\"gmsh:physical\"][key]\n", " for key in mesh.cell_data_dict[\"gmsh:physical\"].keys() if key==cell_type])\n", " mesh = meshio.Mesh(points=mesh.points, cells={cell_type: cells}, \n", " cell_data={\"name_to_read\": [data]})\n", " return mesh\n", "\n", "mesh = meshio.read(filename + \".msh\")\n", "\n", "shell_mesh = create_mesh(mesh, \"triangle\")\n", "line_mesh = create_mesh(mesh, \"line\")\n", "meshio.write(filename + \".xdmf\", shell_mesh)\n", "meshio.write(filename + \"_facet_region.xdmf\", line_mesh)" ] } ], "metadata": { "celltoolbar": "Format de la Cellule Texte Brut", "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.8.10" }, "toc": { "base_numbering": 1, "nav_menu": {}, "number_sections": false, "sideBar": true, "skip_h1_title": false, "title_cell": "Table of Contents", "title_sidebar": "Contents", "toc_cell": false, "toc_position": {}, "toc_section_display": true, "toc_window_display": false } }, "nbformat": 4, "nbformat_minor": 4 }