.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "examples/co2_leakage_along_a_fault/1_generate_mesh_in_python_with_gmsh.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code .. rst-class:: sphx-glr-example-title .. _sphx_glr_examples_co2_leakage_along_a_fault_1_generate_mesh_in_python_with_gmsh.py: Generate mesh in Python with Gmsh ================================= This examples describes how to generate a mesh using Gmsh in Python. It requires :mod:`pygmsh` v7 to be installed: .. code-block:: bash pip install pygmsh --user The mesh can also be generated from scratch directly in Gmsh either using its GUI and/or its internal scripting language. Note that :mod:`toughio` is not required at this preliminary stage of the pre-processing. This example only intends to show how the mesh used in this sample problem has been generated. The user is free to use any meshing software as long as the mesh format is supported by :mod:`toughio` (through :mod:`meshio`). .. GENERATED FROM PYTHON SOURCE LINES 18-19 First, we import :mod:`numpy` and :mod:`pygmsh`. .. GENERATED FROM PYTHON SOURCE LINES 19-23 .. code-block:: Python import numpy as np import pygmsh .. GENERATED FROM PYTHON SOURCE LINES 26-27 We can define a bunch of useful variables such as the characteristic length or some parameters to characterize the model. .. GENERATED FROM PYTHON SOURCE LINES 27-40 .. code-block:: Python lc = 100.0 # Characteristic length of mesh xmin, xmax = 0.0, 2000.0 # X axis boundaries zmin, zmax = -500.0, -2500.0 # Z axis boundaries inj_z = -1500.0 # Depth of injection flt_offset = 50.0 # Offset of fault flt_thick = 10.0 # Thickness of fault tana = np.tan(np.deg2rad(80.0)) # Tangeant of dipping angle of fault dist = 500.0 - 0.5 * flt_thick # Distance from injection point (0.0, -1500.0) to left wall of fault bnd_thick = 10.0 # Thickness of boundary elements .. GENERATED FROM PYTHON SOURCE LINES 43-47 We start by defining the geometrical entity representing the fault. The fault is represented as a finite-thickness element that intersects all the layers of the model. To ensure conformity of the final mesh, each wall of the fault is represented by a segmented line where the positions of the nodes correspond to the intersections of the hanging (left) and foot (right) walls with the different layers. Note that the foot wall is inverted so that the fault entity forms a closed rectangular loop. .. GENERATED FROM PYTHON SOURCE LINES 47-56 .. code-block:: Python depths = [zmin, -1300.0, -1450.0, -1550.0, -1700.0, zmax] fault_left = [[dist + (z - inj_z) / tana, z, 0.0] for z in depths] depths = [zmin, -1300.0 + flt_offset, -1450.0 + flt_offset, -1550.0 + flt_offset, -1700.0 + flt_offset, zmax] fault_right = [[dist + (z - inj_z) / tana + flt_thick, z, 0.0] for z in depths] fault_pts = fault_left + fault_right[::-1] .. GENERATED FROM PYTHON SOURCE LINES 59-60 Now, we define the aquifer located at the left side of the fault. .. GENERATED FROM PYTHON SOURCE LINES 60-92 .. code-block:: Python cenaq_left_pts = [ [xmin, -1450.0, 0.0], [dist + (-1450.0 - inj_z) / tana, -1450.0, 0.0], [dist + (-1550.0 - inj_z) / tana, -1550.0, 0.0], [xmin, -1550.0, 0.0], ] capro_top_left_pts = [ [xmin, -1300.0, 0.0], [dist + (-1300.0 - inj_z) / tana, -1300.0, 0.0], [dist + (-1450.0 - inj_z) / tana, -1450.0, 0.0], [xmin, -1450.0, 0.0], ] capro_bot_left_pts = [ [xmin, -1550.0, 0.0], [dist + (-1550.0 - inj_z) / tana, -1550.0, 0.0], [dist + (-1700.0 - inj_z) / tana, -1700.0, 0.0], [xmin, -1700.0, 0.0], ] uppaq_left_pts = [ [xmin, zmin, 0.0], [dist + (zmin - inj_z) / tana, zmin, 0.0], [dist + (-1300.0 - inj_z) / tana, -1300.0, 0.0], [xmin, -1300.0, 0.0], ] basaq_left_pts = [ [xmin, -1700.0, 0.0], [dist + (-1700.0 - inj_z) / tana, -1700.0, 0.0], [dist + (zmax - inj_z) / tana, zmax, 0.0], [xmin, zmax, 0.0], ] .. GENERATED FROM PYTHON SOURCE LINES 95-96 Likewise, we also define the aquifer located at the right side of the fault. .. GENERATED FROM PYTHON SOURCE LINES 96-128 .. code-block:: Python cenaq_right_pts = [ [dist + (-1450.0 - inj_z + flt_offset) / tana + flt_thick, -1450.0 + flt_offset, 0.0], [xmax, -1450.0 + flt_offset, 0.0], [xmax, -1550.0 + flt_offset, 0.0], [dist + (-1550.0 - inj_z + flt_offset) / tana + flt_thick, -1550.0 + flt_offset, 0.0], ] capro_top_right_pts = [ [dist + (-1300.0 - inj_z + flt_offset) / tana + flt_thick, -1300.0 + flt_offset, 0.0], [xmax, -1300.0 + flt_offset, 0.0], [xmax, -1450.0 + flt_offset, 0.0], [dist + (-1450.0 - inj_z + flt_offset) / tana + flt_thick, -1450.0 + flt_offset, 0.0], ] capro_bot_right_pts = [ [dist + (-1550.0 - inj_z + flt_offset) / tana + flt_thick, -1550.0 + flt_offset, 0.0], [xmax, -1550.0 + flt_offset, 0.0], [xmax, -1700.0 + flt_offset, 0.0], [dist + (-1700.0 - inj_z + flt_offset) / tana + flt_thick, -1700.0 + flt_offset, 0.0], ] uppaq_right_pts = [ [dist + (zmin - inj_z) / tana + flt_thick, zmin, 0.0], [xmax, zmin, 0.0], [xmax, -1300.0 + flt_offset, 0.0], [dist + (-1300.0 - inj_z + flt_offset) / tana + flt_thick, -1300.0 + flt_offset, 0.0], ] basaq_right_pts = [ [dist + (-1700.0 - inj_z + flt_offset) / tana + flt_thick, -1700.0 + flt_offset, 0.0], [xmax, -1700.0 + flt_offset, 0.0], [xmax, zmax, 0.0], [dist + (zmax-inj_z) / tana + flt_thick, zmax, 0.0], ] .. GENERATED FROM PYTHON SOURCE LINES 131-135 Then, we define the boundary elements. In this sample problem, a no-flow boundary condition is imposed on the left side of the model (default in TOUGH), and Dirichlet boundary conditions are imposed elsewhere. Thus, physical boundary elements must be defined at the top, right and bottom sides of the model. Similarly to the fault entity, boundary entities are segmented to ensure conformity of the final mesh. .. GENERATED FROM PYTHON SOURCE LINES 135-169 .. code-block:: Python bound_right_pts = [ [xmax, zmin, 0.0], [xmax + bnd_thick, zmin, 0.0], [xmax + bnd_thick, zmax, 0.0], [xmax, zmax, 0.0], [xmax, -1700.0 + flt_offset, 0.0], [xmax, -1550.0 + flt_offset, 0.0], [xmax, -1450.0 + flt_offset, 0.0], [xmax, -1300.0 + flt_offset, 0.0], ] bound_top_pts = [ [xmin, zmin, 0.0], [dist + (zmin - inj_z) / tana, zmin, 0.0], [dist + (zmin - inj_z) / tana + flt_thick, zmin, 0.0], [xmax, zmin, 0.0], [xmax + bnd_thick, zmin, 0.0], [xmax + bnd_thick, zmin + bnd_thick, 0.0], [dist + (zmin - inj_z + bnd_thick) / tana + flt_thick, zmin + bnd_thick, 0.0], [dist + (zmin - inj_z + bnd_thick) / tana, zmin + bnd_thick, 0.0], [xmin, zmin + bnd_thick, 0.0], ] bound_bot_pts = [ [xmin, zmax, 0.0], [dist + (zmax - inj_z) / tana, zmax, 0.0], [dist + (zmax - inj_z) / tana + flt_thick, zmax, 0.0], [xmax, zmax, 0.0], [xmax + bnd_thick, zmax, 0.0], [xmax + bnd_thick, zmax - bnd_thick, 0.0], [dist + (zmax - inj_z - bnd_thick) / tana + flt_thick, zmax - bnd_thick, 0.0], [dist + (zmax - inj_z - bnd_thick) / tana, zmax - bnd_thick, 0.0], [xmin, zmax - bnd_thick, 0.0], ] .. GENERATED FROM PYTHON SOURCE LINES 172-175 Once all the points have been created, we can now generate the geometry, assign rock types/materials as Gmsh physical properties, and generate the mesh. To refine the mesh in the injection zone, the characteristic length of each layer entity is increased the farther we get from the injection point. Note that layers are defined such that their characteristic lengths are increasing. This is because Gmsh keeps the first node defined in the geometry in case it detects duplicated nodes. .. GENERATED FROM PYTHON SOURCE LINES 175-225 .. code-block:: Python with pygmsh.geo.Geometry() as geo: # Define polygons fault = geo.add_polygon(fault_pts, mesh_size=0.1 * lc) cenaq_left = geo.add_polygon(cenaq_left_pts, mesh_size=0.1 * lc) capro_top_left = geo.add_polygon(capro_top_left_pts, mesh_size=0.2 * lc) capro_bot_left = geo.add_polygon(capro_bot_left_pts, mesh_size=0.2 * lc) uppaq_left = geo.add_polygon(uppaq_left_pts, mesh_size=2.0 * lc) basaq_left = geo.add_polygon(basaq_left_pts, mesh_size=2.0 * lc) cenaq_right = geo.add_polygon(cenaq_right_pts, mesh_size=0.75 * lc) capro_top_right = geo.add_polygon(capro_top_right_pts, mesh_size=0.75 * lc) capro_bot_right = geo.add_polygon(capro_bot_right_pts, mesh_size=0.75 * lc) uppaq_right = geo.add_polygon(uppaq_right_pts, mesh_size=2.0 * lc) basaq_right = geo.add_polygon(basaq_right_pts, mesh_size=2.0 * lc) bound_right = geo.add_polygon(bound_right_pts, mesh_size=lc) bound_top = geo.add_polygon(bound_top_pts, mesh_size=lc) bound_bot = geo.add_polygon(bound_bot_pts, mesh_size=lc) # Define materials geo.add_physical([uppaq_left, uppaq_right], "UPPAQ") geo.add_physical([capro_top_left, capro_bot_left, capro_top_right, capro_bot_right], "CAPRO") geo.add_physical([cenaq_left, cenaq_right], "CENAQ") geo.add_physical([basaq_left, basaq_right], "BASAQ") geo.add_physical(fault, "FAULT") geo.add_physical([bound_right, bound_top, bound_bot], "BOUND") # Remove duplicate entities geo.env.removeAllDuplicates() mesh = geo.generate_mesh(dim=2, algorithm=6) # Convert cell sets to material cell_data = [np.empty(len(c.data), dtype=int) for c in mesh.cells] field_data = {} for i, (k, v) in enumerate(mesh.cell_sets.items()): if k: field_data[k] = np.array([i + 1, 3]) for ii, vv in enumerate(v): cell_data[ii][vv] = i + 1 mesh.cell_data["material"] = cell_data mesh.field_data.update(field_data) mesh.cell_sets = {} # Remove lower dimension entities idx = [i for i, cell in enumerate(mesh.cells) if cell.type == "triangle"] mesh.cells = [mesh.cells[i] for i in idx] mesh.cell_data = {k: [v[i] for i in idx] for k, v in mesh.cell_data.items()} # Export the mesh for post-processing mesh.write("mesh.vtu") .. GENERATED FROM PYTHON SOURCE LINES 228-229 The generated mesh can be visualized in Python with :mod:`pyvista`. .. GENERATED FROM PYTHON SOURCE LINES 229-243 .. code-block:: Python import pyvista pyvista.set_plot_theme("document") p = pyvista.Plotter(window_size=(800, 800)) p.add_mesh( mesh=pyvista.from_meshio(mesh), scalar_bar_args={"title": "Materials"}, show_scalar_bar=True, show_edges=True, ) p.view_xy() p.show() .. image-sg:: /examples/co2_leakage_along_a_fault/images/sphx_glr_1_generate_mesh_in_python_with_gmsh_001.png :alt: 1 generate mesh in python with gmsh :srcset: /examples/co2_leakage_along_a_fault/images/sphx_glr_1_generate_mesh_in_python_with_gmsh_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 2.090 seconds) .. _sphx_glr_download_examples_co2_leakage_along_a_fault_1_generate_mesh_in_python_with_gmsh.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: 1_generate_mesh_in_python_with_gmsh.ipynb <1_generate_mesh_in_python_with_gmsh.ipynb>` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: 1_generate_mesh_in_python_with_gmsh.py <1_generate_mesh_in_python_with_gmsh.py>` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_