{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Checkout www.pygimli.org for more examples\n%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\nCAD to mesh tutorial\n====================\nIn this example you will learn how to create a geometry in\n`FreeCAD `_ and then export and mesh it using\n`Gmsh `_.\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Gmsh comes with a build-in CAD engine for defining a geometry, as shown\nin the `flexible mesh generation example `_,\nbut using a parametric CAD program such as FreeCAD is much more intuitive and flexible.\n\nFor this tutorial you will need Gmsh and its Python API (application\nprogramming interface). These can be installed by the command below\ninside your (new) conda environment.\n\n.. code-block:: bash\n\n conda install -c conda-forge gmsh python-gmsh\n\nIn case you also want to try out FreeCAD, installing it from their\n`website `_ will give you the\nmost up to date version.\n\nThis example is based on an ERT modeling and inversion experiment on a\nsmall dike. However, this FreeCAD \u2192 Gmsh workflow can easily be\ntranslated to other geophysical methods. The geometry and acquisition\ndesign come from the IDEA League `master\nthesis `_\nof Joost Gevaert. The target in this example is to find the geometry\nof a sand channel underneath the dike.\n\nFreeCAD: create the geometry\n----------------------------\n\nTwo geometries have to be created. One for modeling and one for\ninversion. When the same meshes are used for modeling and inversion,\nthe geometry of the sand channel is alreadyincluded in the structure\nof the mesh. Therefore, the mesh itself would act as prior information\nto the inversion. The modeling geometry consists of three regions: the\nouter region; the inner region (same as inversion region in this\nexample) and the sand channel. The inversion geometry consists of two\nregions: the outer region and the inversion region.\n\nThe geometries are defined in three steps:\n\n1. Each region of the geometry designed separately in the Part\n workbench, or in the Part Design workbench for more complicated\n geometries. To get familiar with the part design workbench, this\n `FreeCAD-tutorial `_\n with some videos is great.\n2. Merge all regions into one single \u201dcompsolid\u201d, i.e.composite\n solid. Meaning one object that consists of multiple solids that\n share the interfaces between the solids.\n3. Export the geometry in ``.brep`` *\n\n(1) The outer and inversion regions of this dike example were created\nin the Part Design workbench, by making a sketch and then extruding\nit with the Pad option. See the Inversion-Region in the object tree\nin the figure below. You can also have a look at how these geometries\nwere created by `downloading `_\nthe .FCStd FreeCAD files and playing around with them. The sand\nchannel is a simple cube, created in the Part workbench. Dimensions:\nL = 8.0 m ; W = 15.0 m ; H = 2.0 m. Position: x = 7.5 m ; y = -1.5 m\n; z = -2.3 m.\n\n(2) The trick then lies in merging these shapes into a single\ncompsolid. This is done in the following steps:\n\n1. Open a new project and merge all objects, i.e. regions (File \u2192\n Merge project...) into this project\n2. In the Part workbench, select all objects and create Boolean\n Fragments (Part \u2192 Split \u2192 Boolean Fragments)\n3. Select the newly created BooleanFragments in the object tree and\n change its Mode property to CompSolid, see the figure below.\n4. Keep BooleanFragments selected and then apply a Compound Filter to\n it (Part \u2192 Compound \u2192 Compound Filter)\n5. Quality check the obtained geometry.\n Select the newly created CompoundFilter from the object tree and\n click Check Geometry (Part \u2192 Check Geometry).\n SOLID: in the Shape Content, should match the number of objects\n merged when creating the Boolean Fragments, 3 in this example.\n COMPSOLID: should be 1. Always, also for other geometries.\n COMPOUND: should be 0. Always.\n COMPSOLID: 1 and COMPOUND: 0 indicates that the objects were\n indeed merged correctly to one single compsolid, see the figure\n below.\n\n(3) Select the CompounSolid from the object tree and export (File \u2192\nExport...) as .brep.\n\n.. figure:: ../../_static/cad_tutorial/freecad_fig.png\n :align: center\n\n FreeCAD important dialogs for making a correct compsolid.\n\n* It must be ``.brep``. This is the native format of the OpenCascade\nCAD engine on which both FreeCAD and Gmsh run. ``.step`` (also\n``.stp``) is the standardized CAD exchange format, for some reason\nthis format does not export the shape as a compound solid. Gmsh can\nalso read ``.stl`` and ``.iges`` files. ``.stl`` files only contain\nsurface information and cannot easily be reedited. ``.iges`` is an\nold format for which development stopped after 1996 and geometries\nare not always imported correctly.\n\nGmsh: mesh the geometry\n-----------------------\n\nMeshing with Gmsh is incredibly versatile, but has a very steep\nlearning curve. Here we use the Python Application Programming\nInterface (API). To get familiar with the Python API, the Gmsh\n`tutorials `_\n(`overview `_) were\nconverted to `Python\nscripts `_\nand additional\n`demos `_\nare also provided. I will mention or provide links to relevant\ntutorials and demos, have a look at these for extra context.\n\nLet's start by importing our geometry into Gmsh:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import numpy as np\nimport pygimli as pg\ngmsh = pg.optImport(\"gmsh\", \"do this tutorial. Install by running: pip install gmsh\")\n\n# Download all nessesary files\ngeom_filename = pg.getExampleFile(\"cad/dike_mod.brep\")\nelec_pos_filename = pg.getExampleFile(\"cad/elec_pos.csv\")\n\nif gmsh:\n # Starting it up (tutorial t1.py)\n gmsh.initialize()\n gmsh.option.setNumber(\"General.Terminal\", 1)\n gmsh.model.add(\"dike_mod\")\n # Load a BREP file (t20.py & demo step_assembly.py)\n # .brep files don't contain info about units, so scaling has to be applied\n gmsh.option.setNumber(\"Geometry.OCCScaling\", 0.001)\n volumes = gmsh.model.occ.importShapes(geom_filename)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Before diving into local mesh refinement, putting the electrodes in\nthe mesh and assigning region, boundary and electrode markers, the\n.brep geometry file should be checked. Especially check whether the\nmeshes of two adjacent regions share nodes on their interfaces. The\nmesh can be viewed by running the following lines of code:\n\n.. code-block:: python\n\n # Run this code after every change in the mesh to see what changed.\n gmsh.model.occ.synchronize()\n gmsh.model.mesh.generate(3)\n gmsh.fltk.run()\n\nTips for viewing the mesh:\n\n1. Double left clicking opens a menu in where you can set geometry\n and mesh visibility.\n2. Tools \u2192 Visibility opens a window in which you can select parts of\n the mesh and geometry. Here you can find the tags of the\n elementary entities of the geometry. It is also handy later to QC\n whether physical groups were set correctly.\n3. Clip the mesh and geometry with Tools \u2192 Clipping.\n4. The number of elements ect. can be found in the Tools \u2192 Statistics\n window.\n\nMake sure to quickly write down the Gmsh volume tags of the outer\nregion, dike and channel and the surface tags of the free surface and\nthe underground boundary of the box. You will need this in the next\nstep.\n\nIf importing and meshing the .brep geometry went correctly, great!\nNext we include the electrodes (from Excel file) into the geometry\nand define the Characteristic Length (CL) for each region and the\nelectrodes. The CL is defined at each Point and dictates the mesh\nsize at that point. The mesh size between points is interpolated\nlinearly, by default.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "cl_elec = 0.1\ncl_dike = 0.6\ncl_outer = 30\n# Gmsh geometry tags of relevant parts. Find the tags in the Gmsh interface.\ntags = {\"outer region\": 2,\n \"dike\": 3,\n \"channel\": 1,\n \"surface\": [7, 11, 12, 13, 21, 23, 24,\n 25, 27, 29, 30, 31],\n \"boundary\": [8, 14, 15, 16, 20], # \"Underground Box Boundary\"\n \"electrodes\": []}\n\nif gmsh:\n # Syncronize CAD representation with the Gmsh model (t1.py)\n # Otherwise gmsh.model.get* methods don't work.\n gmsh.model.occ.synchronize()\n # Set mesh sizes for the dike and outer region.\n # The order, in which mesh sizes are set, matters. Big -> Small\n gmsh.model.mesh.setSize( # Especially t16.py, also t2; 15; 18; 21\n gmsh.model.getBoundary( # get dimTags of boundary elements of\n (3, tags[\"outer region\"]), # dimTag: (dim, tag)\n recursive=True), # recursive -> dimTags of points\n cl_outer)\n gmsh.model.mesh.setSize(\n gmsh.model.getBoundary((3, tags[\"dike\"]),recursive=True),\n cl_dike)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now reload the script, mesh the geometry again and have a look how\nthe mesh changed. The next step is adding the electrodes to the mesh.\nThe grid on the dike has 152 electrodes. These points are added in\nGmsh as points 201-352, to prevent clashing with points already\ndefined in the geometry.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "if gmsh:\n # positions: np.array([elec#, x, y, z, y \"over ground\"])\n pos = np.genfromtxt(elec_pos_filename, delimiter=\",\", skip_header=1)\n # Electrodes are put at 2 cm depth, such that they can be embeded in the volume of the dike.\n # Embeding the electrodes into the surface elements complicates meshing.\n elec_depth = 0.02 # elec depth [m]\n pos[:, 3] = pos[:, 3] - elec_depth\n # Add the electrodes to the Gmsh model and put the tags into the Dict\n for xyz in pos:\n tag = int(200 + xyz[0])\n gmsh.model.occ.addPoint(xyz[1], xyz[2], xyz[3], cl_elec, tag)\n tags[\"electrodes\"].append(tag)\n # Embed electrodes in dike volume. (t15.py)\n gmsh.model.occ.synchronize()\n gmsh.model.mesh.embed(0, tags[\"electrodes\"], 3, tags[\"dike\"])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Reload the Gmsh script and mesh it again to see the result. Further\nmesh refinement is then possible with so-called background fields.\nTaking a quick look at Gmsh tutorial\n`t10.geo `_\nis highly recommended. It shows a wide range of possible background\nfields. In this example a Distance field is defined from the\nelectrodes and then a Threshold field is applied as the background\nfield:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# LcMax - /------------------\n# /\n# /\n# /\n# LcMin -o----------------/\n# | | |\n# Point DistMin DistMax\n# Field 1: Distance to electrodes\n\nif gmsh:\n gmsh.model.mesh.field.add(\"Distance\", 1)\n gmsh.model.mesh.field.setNumbers(1, \"NodesList\", tags[\"electrodes\"])\n # Field 2: Threshold that dictates the mesh size of the background field\n gmsh.model.mesh.field.add(\"Threshold\", 2)\n gmsh.model.mesh.field.setNumber(2, \"IField\", 1)\n gmsh.model.mesh.field.setNumber(2, \"LcMin\", cl_elec)\n gmsh.model.mesh.field.setNumber(2, \"LcMax\", cl_dike)\n gmsh.model.mesh.field.setNumber(2, \"DistMin\", 0.2)\n gmsh.model.mesh.field.setNumber(2, \"DistMax\", 1.5)\n gmsh.model.mesh.field.setNumber(2, \"StopAtDistMax\", 1)\n gmsh.model.mesh.field.setAsBackgroundMesh(2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Again reload the Gmsh script and mesh it, to see the result. As the\nlast step in creating the mesh, the physical groups have to be\ndefined, such PyGIMLi recognize regions, boundaries and the\nelectrodes, see ``help(pg.meshtools.readGmsh)``. Make sure to follow\nthe same Physical Group tag number conventions for marking the\nregions, surfaces and points as used in PyGIMLi.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "if gmsh:\n # Physical Volumes, \"Regions\" in pyGIMLi\n pgrp = gmsh.model.addPhysicalGroup(3, [tags[\"outer region\"]], 1) #(dim, tag, pgrp tag)\n gmsh.model.setPhysicalName(3, pgrp, \"Outer Region\") # Physical group name in Gmsh\n pgrp = gmsh.model.addPhysicalGroup(3, [tags[\"dike\"]], 2)\n gmsh.model.setPhysicalName(3, pgrp, \"Dike\")\n pgrp = gmsh.model.addPhysicalGroup(3, [tags[\"channel\"]], 3)\n gmsh.model.setPhysicalName(3, pgrp, \"Channel\")\n # Physical Surfaces, \"Boundaries\" in pyGIMLi,\n # pgrp tag = 1 --> Free Surface | pgrp tag > 1 --> Mixed BC\n pgrp = gmsh.model.addPhysicalGroup(2, tags[\"surface\"], 1)\n gmsh.model.setPhysicalName(2, pgrp, \"Surface\")\n pgrp = gmsh.model.addPhysicalGroup(2, tags[\"boundary\"], 2)\n gmsh.model.setPhysicalName(2, pgrp, \"Underground Boundary\")\n # Physical Points, \"Electrodes / Sensors\" in pyGIMLi, pgrp tag 99\n pgrp = gmsh.model.addPhysicalGroup(0, tags[\"electrodes\"], 99)\n gmsh.model.setPhysicalName(0, pgrp, \"Electrodes\")\n\n # Generate the mesh and write the mesh file\n gmsh.model.occ.synchronize()\n gmsh.model.mesh.generate(3)\n gmsh.write(\"dike_mod.msh\")\n gmsh.finalize()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The final mesh should look something like the figure below. Check\nwhether the Physical Groups are defined correctly using the\nVisibility window as shown in the figure. Finally make the inversion\nmesh in the same way as the modeling mesh. The differences being that\n(a) there should be no sand channel in the geometry of the inversion\nmesh. Meaning that there are also only 2 volumes: the outer region\nand the dike, i.e. inversion region. And (b) that the mesh does not\nhave to be as fine. cl elec = 0.25 and cl dike = 1.2 were used for\nthe inversion mesh in the attached .msh file. Besides changing the\nmesh size by playing around with the CL, the general mesh size can\nalso be adapted in Gmsh by changing the Global mesh size factor\n(double left click).\n\n.. figure:: ../../_static/cad_tutorial/gmsh_fig.png\n :align: center\n\n Gmsh visualizatin of the modeling mesh of the dike, with visibility dialog.\n\nAdditional very useful material\n-------------------------------\n\n- `Meshing terrain from a .STL file with\n Gmsh `_\n- `Meshing with Gmsh from\n QGIS `_\n- `FreeCAD GeoMatics\n workbench `_\n (`replaces GeoData\n workbench `_)\n allows for GPS, LiDAR and GIS data to be imported to FreeCAD\n\n" ] } ], "metadata": { "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.7.3" } }, "nbformat": 4, "nbformat_minor": 0 }