Creating OpenGeoSys Mesh with Inclined Borehole Heat Exchangers

This page is based on a Jupyter notebook.

This tutorial is made to illustrate the procedure of creating an OGS mesh file with inclined Borehole Heat Exchangers (BHEs) in it. Such mesh uses prism, tetrahedron and pyramid elements for the soil part, and line elements for the BHEs. In this tutorial of inclined BHEs, a layer of hexagonal shape prism elements are created around each BHE for optimal accuracy (Diersch et al. 2011) and all other parts of the geometry consist of tetrahedron and pyramid elements. The produced mesh file is made explicitly for the HEAT_TRANSPORT_BHE module in OGS and will NOT work with other modules. For better understanding, an image of 1D inclined BHEs is presented.

First, external packages have been imported and Gmsh is initialized.

import os… (click to toggle)

import os
import sys
from pathlib import Path

import gmsh

gmsh.initialize()

The geometry is a 3D structure that has 3 boreholes (2 inclined and 1 vertical) in it. Similar to the previous tutorial of BHEs, the first step is to create the surface 1 with all the necessary points which regulate the borehole locations, as well as the mesh size. Now, we define the basic geometry of the BHEs, as well as the element sizes around them.

# environment variable for output path… (click to toggle)

# environment variable for output path
out_dir = os.environ.get("OGS_TESTRUNNER_OUT_DIR", "_out")
if not Path(out_dir).exists():
    Path(out_dir).mkdir(parents=True)

exe_dir = os.environ.get("OGS_BINARY_DIR")

# mesh file name
bhe_mesh_file_name = "inclined_bhe_mesh_file"

# geometry parameters
width = 20
length = 20
depth = 12
bhe_depth = depth - 2

# element sizes
bhe_radius = 0.07
alpha = 6.134  # see Diersch et al. 2011 Part 2
delta = (
    alpha * bhe_radius
)  # meshsize at BHE and distance of the surrounding optimal mesh points
elem_size = 0.5

In this step, we are going to create the top surface using the python interface of Gmsh. To create a point with the built-in CAD kernel, the Python API function gmsh.model.geo.addPoint() is used. The first 3 arguments are the point coordinates (x, y, z). The next (optional) argument is the target mesh size close to the point. The last (optional) argument is the point tag (a strictly positive integer that uniquely identifies the point). Here, we have assigned 4 boundary points.

# cube surface… (click to toggle)

# cube surface
gmsh.model.geo.addPoint(-width / 2.0, 0.0, 0.0, elem_size, 1)
gmsh.model.geo.addPoint(width / 2.0, 0.0, 0.0, elem_size, 2)
gmsh.model.geo.addPoint(width / 2.0, length, 0.0, elem_size, 3)
gmsh.model.geo.addPoint(-width / 2.0, length, 0.0, elem_size, 4)

Next, the API function gmsh.model.geo.addLine is used to create straight-line segments with the built-in kernel follows the same conventions: the first 2 arguments are point tags (the start and end points of the line), and the last (optional) is the line tag. Note that curve tags are separate from point tags. Hence we can reuse tag ‘1’ for our first curve. And as a general rule, elementary entity tags in Gmsh have to be unique per geometrical dimension.

gmsh.model.geo.addLine(1, 2, 1)… (click to toggle)

gmsh.model.geo.addLine(1, 2, 1)
gmsh.model.geo.addLine(2, 3, 2)
gmsh.model.geo.addLine(3, 4, 3)
gmsh.model.geo.addLine(4, 1, 4)

Next, a curve loop has been defined by using the ordered list of 4 connected lines. The API function to create curve loops takes a list of integers as the first argument, and the curve loop tag (which must be unique amongst curve loops) as the second (optional) argument.

gmsh.model.geo.addCurveLoop([1, 2, 3, 4], 1)

In this step, the structure of the geometry is added. In order to create the structure, bottom surface and other surrounding surfaces have been created using similar method to the top surface. Hence, four points of the bottom surface have been added. Later, all necessary lines have been added to build the 3D geometry. Next, 5 curve loops (4 surroundings and 1 bottom) have been defined and corresponding plane surface is added using the API function gmsh.model.geo.addPlaneSurface.

gmsh.model.geo.addPoint(-width / 2.0, 0.0, -depth, elem_size, 101)… (click to toggle)

gmsh.model.geo.addPoint(-width / 2.0, 0.0, -depth, elem_size, 101)
gmsh.model.geo.addPoint(width / 2.0, 0.0, -depth, elem_size, 102)
gmsh.model.geo.addPoint(width / 2.0, length, -depth, elem_size, 103)
gmsh.model.geo.addPoint(-width / 2.0, length, -depth, elem_size, 104)

gmsh.model.geo.addLine(1, 101, 101)
gmsh.model.geo.addLine(2, 102, 102)
gmsh.model.geo.addLine(3, 103, 103)
gmsh.model.geo.addLine(4, 104, 104)
gmsh.model.geo.addLine(101, 102, 105)
gmsh.model.geo.addLine(102, 103, 106)
gmsh.model.geo.addLine(103, 104, 107)
gmsh.model.geo.addLine(104, 101, 108)


gmsh.model.geo.addCurveLoop([105, 106, 107, 108], 5)
gmsh.model.geo.addCurveLoop([2, 103, -106, -102], 6)
gmsh.model.geo.addCurveLoop([3, 104, -107, -103], 7)
gmsh.model.geo.addCurveLoop([4, 101, -108, -104], 8)
gmsh.model.geo.addCurveLoop([1, 102, -105, -101], 9)

gmsh.model.geo.addPlaneSurface([5], 5)
gmsh.model.geo.addPlaneSurface([6], 6)
gmsh.model.geo.addPlaneSurface([7], 7)
gmsh.model.geo.addPlaneSurface([8], 8)
gmsh.model.geo.addPlaneSurface([9], 9)

Now, 3 BHE surfaces have been created using parameters such as BHE coordinates, delta. The BHE node is surrounded by the 6 additional nodes in hexagonal shape which will create better mesh for optimal accuracy (Diersch et al. 2011 Part 2 DOI:10.1016/j.cageo.2010.08.002).

For each BHE node, 6 additional nodes and 6 corresponding line are added. After that, a curve loop is created using 6 lines. Then, the BHE surface is defined using that curve loop. Following the same method, BHE surface 2 and 3 have been created.

# BHE surfaces 1… (click to toggle)

# BHE surfaces 1
gmsh.model.geo.addPoint(-6, 10, 0, delta, 5)
gmsh.model.geo.addPoint(-6, 10 - delta, 0, delta, 6)
gmsh.model.geo.addPoint(-6, 10 + delta, 0, delta, 7)
gmsh.model.geo.addPoint(-6 + 0.866 * delta, 10 + 0.5 * delta, 0, delta, 8)
gmsh.model.geo.addPoint(-6 - 0.866 * delta, 10 + 0.5 * delta, 0, delta, 9)
gmsh.model.geo.addPoint(-6 + 0.866 * delta, 10 - 0.5 * delta, 0, delta, 10)
gmsh.model.geo.addPoint(-6 - 0.866 * delta, 10 - 0.5 * delta, 0, delta, 11)

gmsh.model.geo.addLine(7, 8, 5)
gmsh.model.geo.addLine(8, 10, 6)
gmsh.model.geo.addLine(10, 6, 7)
gmsh.model.geo.addLine(6, 11, 8)
gmsh.model.geo.addLine(11, 9, 9)
gmsh.model.geo.addLine(9, 7, 10)

gmsh.model.geo.addCurveLoop([5, 6, 7, 8, 9, 10], 2)
gmsh.model.geo.addPlaneSurface([2], 2)

# BHE surfaces 2
gmsh.model.geo.addPoint(6, 10, 0, delta, 12)
gmsh.model.geo.addPoint(6, 10 - delta, 0, delta, 13)
gmsh.model.geo.addPoint(6, 10 + delta, 0, delta, 14)
gmsh.model.geo.addPoint(6 + 0.866 * delta, 10 + 0.5 * delta, 0, delta, 15)
gmsh.model.geo.addPoint(6 - 0.866 * delta, 10 + 0.5 * delta, 0, delta, 16)
gmsh.model.geo.addPoint(6 + 0.866 * delta, 10 - 0.5 * delta, 0, delta, 17)
gmsh.model.geo.addPoint(6 - 0.866 * delta, 10 - 0.5 * delta, 0, delta, 18)

gmsh.model.geo.addLine(14, 15, 11)
gmsh.model.geo.addLine(15, 17, 12)
gmsh.model.geo.addLine(17, 13, 13)
gmsh.model.geo.addLine(13, 18, 14)
gmsh.model.geo.addLine(18, 16, 15)
gmsh.model.geo.addLine(16, 14, 16)

gmsh.model.geo.addCurveLoop([11, 12, 13, 14, 15, 16], 3)
gmsh.model.geo.addPlaneSurface([3], 3)

# BHE surfaces 3
gmsh.model.geo.addPoint(0, 10, 0, delta, 19)
gmsh.model.geo.addPoint(0, 10 - delta, 0, delta, 20)
gmsh.model.geo.addPoint(0, 10 + delta, 0, delta, 21)
gmsh.model.geo.addPoint(0 + 0.866 * delta, 10 + 0.5 * delta, 0, delta, 22)
gmsh.model.geo.addPoint(0 - 0.866 * delta, 10 + 0.5 * delta, 0, delta, 23)
gmsh.model.geo.addPoint(0 + 0.866 * delta, 10 - 0.5 * delta, 0, delta, 24)
gmsh.model.geo.addPoint(0 - 0.866 * delta, 10 - 0.5 * delta, 0, delta, 25)

gmsh.model.geo.addLine(21, 22, 17)
gmsh.model.geo.addLine(22, 24, 18)
gmsh.model.geo.addLine(24, 20, 19)
gmsh.model.geo.addLine(20, 25, 20)
gmsh.model.geo.addLine(25, 23, 21)
gmsh.model.geo.addLine(23, 21, 22)

gmsh.model.geo.addCurveLoop([17, 18, 19, 20, 21, 22], 4)
gmsh.model.geo.addPlaneSurface([4], 4)
gmsh.model.geo.synchronize()

The gmsh.model.geo.extrude command extrudes BHE surface 1, 2 and 3 along the z axis and automatically creates a new volume (as well as all the needed points, curves and surfaces). The function takes a vector of (dim, tag) pairs as input as well as the translation vector, and returns a vector of (dim, tag) pairs as output. For the BHE 1 and 2 in translational vector, x-coordinate is set to 2 and -2 respectively in order to create inclined BHE along x-direction. As a result, the end point of BHEs will be deflected by 2 units along (+) and (-) x-direction.

The 2D mesh extrusion is done with the same extrude() function, but by specifying element ‘Layers’ (Here, one layer each with 10 subdivisions). The number of elements for each layer and the (end) height of each layer are specified in two vectors. The last (optional) argument for the extrude() function specifies whether the extruded mesh should be recombined or not. In this case, it is True since we want to recombine and produce prism mesh elements.

Later, BHE points 5, 12 and 19 have been extruded in the same way to create line elements (Lr, Ll and Lm) which represent BHEs.

P = gmsh.model.geo.extrude([(2, 2)], 2, 0, -10, [10], [1], True) # BHE surface 1… (click to toggle)

P = gmsh.model.geo.extrude([(2, 2)], 2, 0, -10, [10], [1], True)  # BHE surface 1
Q = gmsh.model.geo.extrude([(2, 3)], -2, 0, -10, [10], [1], True)  # BHE surface 2
R = gmsh.model.geo.extrude([(2, 4)], 0, 0, -10, [10], [1], True)  # BHE surface 3

Lr = gmsh.model.geo.extrude([(0, 5)], 2, 0, -10, [10], [1], True)  # BHE 1
Ll = gmsh.model.geo.extrude([(0, 12)], -2, 0, -10, [10], [1], True)  # BHE 2
Lm = gmsh.model.geo.extrude([(0, 19)], 0, 0, -10, [10], [1], True)  # BHE 3

gmsh.model.geo.synchronize()

Finally, geometry surface 1 is created by combining curve loop 1, 2, 3 and 4. Here, 2, 3, 4 are BHE curve loops and 1 is geometry curve loop. After that, 6 surfaces of geometry and 21 surfaces from 3 BHE hexagonal objects are combined to create a surface loop. Later, volume is added to the surface loop using API function gmsh.model.geo.addVolume. Before creating mesh, gmsh.model.geo.synchronize is used to synchronize the CAD entities with the Gmsh model, which will create the relevant Gmsh data structures.

gmsh.model.geo.addPlaneSurface([1, 2, 3, 4], 1)… (click to toggle)

gmsh.model.geo.addPlaneSurface([1, 2, 3, 4], 1)

sl = gmsh.model.geo.addSurfaceLoop(
    [
        1,
        5,
        6,
        7,
        8,
        9,
        204,
        183,
        187,
        191,
        195,
        199,
        203,
        140,
        119,
        123,
        127,
        131,
        135,
        139,
        172,
        151,
        155,
        159,
        163,
        167,
        171,
    ]
)
v = gmsh.model.geo.addVolume([sl])

gmsh.model.geo.synchronize()

Later gmsh.model.addPhysicalGroup command used to group elementary geometrical entities to define material properties. Gmsh will export in output files only mesh elements that belong to at least one physical group. Physical groups are also identified by tags, i.e. strictly positive integers, that should be unique per dimension (0D, 1D, 2D or 3D). Here, BHE 1, 2 and 3 are tagged by 2, 3, 4 physical group respectively.

gmsh.model.addPhysicalGroup(3, [P[1][1], Q[1][1], R[1][1], v], 1)… (click to toggle)

gmsh.model.addPhysicalGroup(3, [P[1][1], Q[1][1], R[1][1], v], 1)
gmsh.model.addPhysicalGroup(1, [Lr[1][1]], 2)
gmsh.model.addPhysicalGroup(1, [Ll[1][1]], 3)
gmsh.model.addPhysicalGroup(1, [Lm[1][1]], 4)

Meshes generated with Gmsh must be converted to VTU file format later. Currently, the only supported Gmsh format is 2.2

gmsh.option.setNumber("Mesh.MshFileVersion", 2.2)

Then We can then generate a 3D mesh and save it to disk

gmsh.model.mesh.generate(3)… (click to toggle)

gmsh.model.mesh.generate(3)

gmsh.model.mesh.removeDuplicateNodes()

gmsh.write(f"{out_dir}/{bhe_mesh_file_name}.msh")

Info    : Meshing 1D...
Info    : [  0%] Meshing curve 1 (Line)
Info    : [ 10%] Meshing curve 2 (Line)
Info    : [ 10%] Meshing curve 3 (Line)
Info    : [ 10%] Meshing curve 4 (Line)
Info    : [ 10%] Meshing curve 5 (Line)
Info    : [ 10%] Meshing curve 6 (Line)
Info    : [ 10%] Meshing curve 7 (Line)
Info    : [ 20%] Meshing curve 8 (Line)
Info    : [ 20%] Meshing curve 9 (Line)
Info    : [ 20%] Meshing curve 10 (Line)
Info    : [ 20%] Meshing curve 11 (Line)
Info    : [ 20%] Meshing curve 12 (Line)
Info    : [ 20%] Meshing curve 13 (Line)
Info    : [ 20%] Meshing curve 14 (Line)
Info    : [ 30%] Meshing curve 15 (Line)
Info    : [ 30%] Meshing curve 16 (Line)
Info    : [ 30%] Meshing curve 17 (Line)
Info    : [ 30%] Meshing curve 18 (Line)
Info    : [ 30%] Meshing curve 19 (Line)
Info    : [ 30%] Meshing curve 20 (Line)
Info    : [ 30%] Meshing curve 21 (Line)
Info    : [ 40%] Meshing curve 22 (Line)
Info    : [ 40%] Meshing curve 101 (Line)
Info    : [ 40%] Meshing curve 102 (Line)
Info    : [ 40%] Meshing curve 103 (Line)
Info    : [ 40%] Meshing curve 104 (Line)
Info    : [ 40%] Meshing curve 105 (Line)
Info    : [ 40%] Meshing curve 106 (Line)
Info    : [ 50%] Meshing curve 107 (Line)
Info    : [ 50%] Meshing curve 108 (Line)
Info    : [ 50%] Meshing curve 110 (Extruded)
Info    : [ 50%] Meshing curve 111 (Extruded)
Info    : [ 50%] Meshing curve 112 (Extruded)
Info    : [ 50%] Meshing curve 113 (Extruded)
Info    : [ 50%] Meshing curve 114 (Extruded)
Info    : [ 60%] Meshing curve 115 (Extruded)
Info    : [ 60%] Meshing curve 117 (Extruded)
Info    : [ 60%] Meshing curve 118 (Extruded)
Info    : [ 60%] Meshing curve 122 (Extruded)
Info    : [ 60%] Meshing curve 126 (Extruded)
Info    : [ 60%] Meshing curve 130 (Extruded)
Info    : [ 60%] Meshing curve 134 (Extruded)
Info    : [ 70%] Meshing curve 142 (Extruded)
Info    : [ 70%] Meshing curve 143 (Extruded)
Info    : [ 70%] Meshing curve 144 (Extruded)
Info    : [ 70%] Meshing curve 145 (Extruded)
Info    : [ 70%] Meshing curve 146 (Extruded)
Info    : [ 70%] Meshing curve 147 (Extruded)
Info    : [ 70%] Meshing curve 149 (Extruded)
Info    : [ 80%] Meshing curve 150 (Extruded)
Info    : [ 80%] Meshing curve 154 (Extruded)
Info    : [ 80%] Meshing curve 158 (Extruded)
Info    : [ 80%] Meshing curve 162 (Extruded)
Info    : [ 80%] Meshing curve 166 (Extruded)
Info    : [ 80%] Meshing curve 174 (Extruded)
Info    : [ 80%] Meshing curve 175 (Extruded)
Info    : [ 90%] Meshing curve 176 (Extruded)
Info    : [ 90%] Meshing curve 177 (Extruded)
Info    : [ 90%] Meshing curve 178 (Extruded)
Info    : [ 90%] Meshing curve 179 (Extruded)
Info    : [ 90%] Meshing curve 181 (Extruded)
Info    : [ 90%] Meshing curve 182 (Extruded)
Info    : [ 90%] Meshing curve 186 (Extruded)
Info    : [100%] Meshing curve 190 (Extruded)
Info    : [100%] Meshing curve 194 (Extruded)
Info    : [100%] Meshing curve 198 (Extruded)
Info    : [100%] Meshing curve 205 (Extruded)
Info    : [100%] Meshing curve 206 (Extruded)
Info    : [100%] Meshing curve 207 (Extruded)
Info    : Done meshing 1D (Wall 0.00144701s, CPU 0.00247s)
Info    : Meshing 2D...
Info    : [  0%] Meshing surface 1 (Plane, Frontal-Delaunay)
Info    : [ 10%] Meshing surface 2 (Plane, Frontal-Delaunay)
Info    : [ 10%] Meshing surface 3 (Plane, Frontal-Delaunay)
Info    : [ 20%] Meshing surface 4 (Plane, Frontal-Delaunay)
Info    : [ 20%] Meshing surface 5 (Plane, Frontal-Delaunay)
Info    : [ 20%] Meshing surface 6 (Plane, Frontal-Delaunay)
Info    : [ 30%] Meshing surface 7 (Plane, Frontal-Delaunay)
Info    : [ 30%] Meshing surface 8 (Plane, Frontal-Delaunay)
Info    : [ 30%] Meshing surface 9 (Plane, Frontal-Delaunay)
Info    : [ 40%] Meshing surface 119 (Extruded)
Info    : [ 40%] Meshing surface 123 (Extruded)
Info    : [ 40%] Meshing surface 127 (Extruded)
Info    : [ 50%] Meshing surface 131 (Extruded)
Info    : [ 50%] Meshing surface 135 (Extruded)
Info    : [ 50%] Meshing surface 139 (Extruded)
Info    : [ 60%] Meshing surface 140 (Extruded)
Info    : [ 60%] Meshing surface 151 (Extruded)
Info    : [ 60%] Meshing surface 155 (Extruded)
Info    : [ 70%] Meshing surface 159 (Extruded)
Info    : [ 70%] Meshing surface 163 (Extruded)
Info    : [ 70%] Meshing surface 167 (Extruded)
Info    : [ 80%] Meshing surface 171 (Extruded)
Info    : [ 80%] Meshing surface 172 (Extruded)
Info    : [ 80%] Meshing surface 183 (Extruded)
Info    : [ 90%] Meshing surface 187 (Extruded)
Info    : [ 90%] Meshing surface 191 (Extruded)
Info    : [ 90%] Meshing surface 195 (Extruded)
Info    : [100%] Meshing surface 199 (Extruded)
Info    : [100%] Meshing surface 203 (Extruded)
Info    : [100%] Meshing surface 204 (Extruded)
Info    : Done meshing 2D (Wall 0.121486s, CPU 0.118755s)
Info    : Meshing 3D...
Info    : Meshing volume 1 (Extruded)
Info    : Meshing volume 2 (Extruded)
Info    : Meshing volume 3 (Extruded)
Info    : 3D Meshing 1 volume with 1 connected component
Info    : Tetrahedrizing 8829 nodes...
Info    : Done tetrahedrizing 8837 nodes (Wall 0.0708499s, CPU 0.069992s)

Info    : Reconstructing mesh...
Info    :  - Creating surface mesh
Info    :  - Identifying boundary edges
Info    :  - Recovering boundary
Info    : Done reconstructing mesh (Wall 0.1974s, CPU 0.191186s)
Info    : Found volume 4

Info    : It. 0 - 0 nodes created - worst tet radius 11.0009 (nodes removed 0 0)
Info    : It. 500 - 500 nodes created - worst tet radius 3.32294 (nodes removed 0 0)
Info    : It. 1000 - 1000 nodes created - worst tet radius 2.71805 (nodes removed 0 0)
Info    : It. 1500 - 1500 nodes created - worst tet radius 2.42141 (nodes removed 0 0)
Info    : It. 2000 - 1998 nodes created - worst tet radius 2.22495 (nodes removed 0 2)
Info    : It. 2500 - 2498 nodes created - worst tet radius 2.07679 (nodes removed 0 2)
Info    : It. 3000 - 2995 nodes created - worst tet radius 1.97088 (nodes removed 0 5)

Info    : It. 3500 - 3491 nodes created - worst tet radius 1.87597 (nodes removed 0 9)
Info    : It. 4000 - 3982 nodes created - worst tet radius 1.80491 (nodes removed 0 18)
Info    : It. 4500 - 4461 nodes created - worst tet radius 1.74426 (nodes removed 0 39)
Info    : It. 5000 - 4929 nodes created - worst tet radius 1.69201 (nodes removed 1 70)
Info    : It. 5500 - 5391 nodes created - worst tet radius 1.64551 (nodes removed 1 108)
Info    : It. 6000 - 5852 nodes created - worst tet radius 1.60605 (nodes removed 2 146)
Info    : It. 6500 - 6267 nodes created - worst tet radius 1.57436 (nodes removed 11 222)
Info    : It. 7000 - 6683 nodes created - worst tet radius 1.54205 (nodes removed 17 300)
Info    : It. 7500 - 7131 nodes created - worst tet radius 1.5109 (nodes removed 18 351)
Info    : It. 8000 - 7582 nodes created - worst tet radius 1.48016 (nodes removed 20 398)

Info    : It. 8500 - 8039 nodes created - worst tet radius 1.45303 (nodes removed 20 441)
Info    : It. 9000 - 8526 nodes created - worst tet radius 1.42546 (nodes removed 25 449)
Info    : It. 9500 - 9020 nodes created - worst tet radius 1.39807 (nodes removed 26 454)
Info    : It. 10000 - 9508 nodes created - worst tet radius 1.37355 (nodes removed 26 466)
Info    : It. 10500 - 9999 nodes created - worst tet radius 1.35375 (nodes removed 26 475)
Info    : It. 11000 - 10475 nodes created - worst tet radius 1.33413 (nodes removed 30 495)
Info    : It. 11500 - 10971 nodes created - worst tet radius 1.31699 (nodes removed 32 497)
Info    : It. 12000 - 11450 nodes created - worst tet radius 1.29948 (nodes removed 32 518)
Info    : It. 12500 - 11947 nodes created - worst tet radius 1.28227 (nodes removed 34 519)

Info    : It. 13000 - 12441 nodes created - worst tet radius 1.26672 (nodes removed 36 523)
Info    : It. 13500 - 12915 nodes created - worst tet radius 1.59762 (nodes removed 38 547)
Info    : It. 14000 - 13390 nodes created - worst tet radius 1.23693 (nodes removed 42 568)
Info    : It. 14500 - 13867 nodes created - worst tet radius 1.22324 (nodes removed 42 591)
Info    : It. 15000 - 14338 nodes created - worst tet radius 1.20993 (nodes removed 44 618)
Info    : It. 15500 - 14822 nodes created - worst tet radius 1.19817 (nodes removed 44 634)
Info    : It. 16000 - 15307 nodes created - worst tet radius 1.18664 (nodes removed 47 646)
Info    : It. 16500 - 15803 nodes created - worst tet radius 1.17376 (nodes removed 47 650)
Info    : It. 17000 - 16289 nodes created - worst tet radius 1.16171 (nodes removed 47 664)
Info    : It. 17500 - 16762 nodes created - worst tet radius 1.1504 (nodes removed 51 687)
Info    : It. 18000 - 17235 nodes created - worst tet radius 1.14108 (nodes removed 51 714)

Info    : It. 18500 - 17710 nodes created - worst tet radius 1.13065 (nodes removed 53 737)
Info    : It. 19000 - 18200 nodes created - worst tet radius 1.12139 (nodes removed 54 746)
Info    : It. 19500 - 18695 nodes created - worst tet radius 1.11225 (nodes removed 55 750)
Info    : It. 20000 - 19183 nodes created - worst tet radius 1.10326 (nodes removed 56 761)
Info    : It. 20500 - 19662 nodes created - worst tet radius 1.09408 (nodes removed 57 781)
Info    : It. 21000 - 20143 nodes created - worst tet radius 1.08515 (nodes removed 58 799)
Info    : It. 21500 - 20622 nodes created - worst tet radius 1.07752 (nodes removed 61 817)
Info    : It. 22000 - 21111 nodes created - worst tet radius 1.06929 (nodes removed 61 828)
Info    : It. 22500 - 21598 nodes created - worst tet radius 1.0618 (nodes removed 61 841)
Info    : It. 23000 - 22084 nodes created - worst tet radius 1.05421 (nodes removed 61 855)

Info    : It. 23500 - 22563 nodes created - worst tet radius 1.04737 (nodes removed 64 873)
Info    : It. 24000 - 23057 nodes created - worst tet radius 1.03968 (nodes removed 64 879)
Info    : It. 24500 - 23547 nodes created - worst tet radius 1.67283 (nodes removed 67 886)
Info    : It. 25000 - 24037 nodes created - worst tet radius 1.02645 (nodes removed 67 896)
Info    : It. 25500 - 24531 nodes created - worst tet radius 1.01986 (nodes removed 67 902)
Info    : It. 26000 - 25012 nodes created - worst tet radius 1.01376 (nodes removed 69 919)

Info    : It. 26500 - 25502 nodes created - worst tet radius 1.00753 (nodes removed 71 927)
Info    : It. 27000 - 25990 nodes created - worst tet radius 1.00137 (nodes removed 72 938)
Info    : 3D refinement terminated (34932 nodes total):
Info    :  - 299 Delaunay cavities modified for star shapeness
Info    :  - 1026 nodes could not be inserted
Info    :  - 194333 tetrahedra created in 1.32214 sec. (146983 tets/s)
Info    : 287 node relocations

Info    : Generating pyramids for hybrid mesh...
Info    : Done generating 180 pyramids for hybrid mesh
Info    : Optimizing pyramids for hybrid mesh...
Info    : Using new pyramid optimization

 detmin =  5.272180366111974E-03 eps=  1.000000000000000E-06 59 iter term 1
 detmin =  5.272179788285682E-03 eps=  1.000000000000000E-06 1 iter term 1
Info    : Done optimizing pyramids for hybrid mesh

Info    : Done meshing 3D (Wall 3.45751s, CPU 3.42692s)
Info    : Optimizing mesh...
Info    : Optimizing volume 4

Info    : Optimization starts (volume = 4780.02) with worst = 0.000318721 / average = 0.781741:
Info    : 0.00 < quality < 0.10 :       439 elements
Info    : 0.10 < quality < 0.20 :      1261 elements
Info    : 0.20 < quality < 0.30 :      2177 elements
Info    : 0.30 < quality < 0.40 :      3526 elements
Info    : 0.40 < quality < 0.50 :      5356 elements
Info    : 0.50 < quality < 0.60 :      8829 elements
Info    : 0.60 < quality < 0.70 :     18524 elements
Info    : 0.70 < quality < 0.80 :     42819 elements
Info    : 0.80 < quality < 0.90 :     73139 elements
Info    : 0.90 < quality < 1.00 :     38259 elements
Info    : 3776 edge swaps, 0 node relocations (volume = 4780.02): worst = 0.096823 / average = 0.793933 (Wall 0.10177s, CPU 0.101843s)
Info    : 3834 edge swaps, 0 node relocations (volume = 4780.02): worst = 0.096823 / average = 0.794038 (Wall 0.140097s, CPU 0.140115s)
Info    : No ill-shaped tets in the mesh :-)
Info    : 0.00 < quality < 0.10 :         2 elements
Info    : 0.10 < quality < 0.20 :        18 elements
Info    : 0.20 < quality < 0.30 :        63 elements
Info    : 0.30 < quality < 0.40 :      3433 elements
Info    : 0.40 < quality < 0.50 :      5145 elements
Info    : 0.50 < quality < 0.60 :      8535 elements
Info    : 0.60 < quality < 0.70 :     18471 elements
Info    : 0.70 < quality < 0.80 :     43472 elements
Info    : 0.80 < quality < 0.90 :     73781 elements
Info    : 0.90 < quality < 1.00 :     37969 elements
Info    : Done optimizing mesh (Wall 0.507414s, CPU 0.508011s)
Info    : 34992 nodes 209120 elements

Info    : Removing duplicate mesh nodes...
Info    : Found 33 duplicate nodes 
Info    : Removed 33 duplicate mesh nodes
Info    : Done removing duplicate mesh nodes
Info    : Writing '/var/lib/gitlab-runner/builds/F1XUyv4cx/1/ogs/build/release-all/web/content/docs/tutorials/Inclined_bhe_meshing/notebook-inclined_bhe_meshing/inclined_bhe_mesh_file.msh'...

Info    : Done writing '/var/lib/gitlab-runner/builds/F1XUyv4cx/1/ogs/build/release-all/web/content/docs/tutorials/Inclined_bhe_meshing/notebook-inclined_bhe_meshing/inclined_bhe_mesh_file.msh'

Launch the GUI to see the results. Later gmsh.finalize() will be called when done using Gmsh Python API

if "CI" not in os.environ:… (click to toggle)

if "CI" not in os.environ:
    gmsh.fltk.run()

gmsh.finalize()

If everything runs well, you will see the following mesh with incline BHEs.

2D version will look like this.

Now checking whether the Gmsh format mesh file is properly created. If not give an error message.

check_file = Path(f"{out_dir}/{bhe_mesh_file_name}.msh").is_file()… (click to toggle)

check_file = Path(f"{out_dir}/{bhe_mesh_file_name}.msh").is_file()
if check_file:
    print("Creation of BHE mesh in Gmsh format was successful.")
else:
    print("Error! Gmsh file is not properly created in the BHE meshing tutorial.")
    raise SystemExit()

Creation of BHE mesh in Gmsh format was successful.

Finally, the mesh file which has been created using the python interface of Gmsh, will be converted to OGS mesh, in particular to VTU file format. Please, add ogs.exe to the directory of this example file to run GMSH2OGS or give the full path of your ogs executable. Here, option -v (–validation) validate the mesh and shows crucial information of mesh. Option -i takes Gmsh input file name as string and -o is output file name as string as well.

!GMSH2OGS -i {out_dir}/{bhe_mesh_file_name}.msh -o {out_dir}/{bhe_mesh_file_name}.vtu -v

[2024-10-17 11:58:40.782] [ogs] [[32minfo[m] Reading /var/lib/gitlab-runner/builds/F1XUyv4cx/1/ogs/build/release-all/web/content/docs/tutorials/Inclined_bhe_meshing/notebook-inclined_bhe_meshing/inclined_bhe_mesh_file.msh.
[2024-10-17 11:58:40.782] [ogs] [[33m[1mwarning[m] If the mesh is generated with Gmsh version 2 and it is saved (or exported) as "Version 2 ASCII" format, the flag --gmsh2_physical_id must be used for a correct conversion from physical id to MaterialIDs.

[2024-10-17 11:58:42.824] [ogs] [[32minfo[m] 	... finished.
[2024-10-17 11:58:42.824] [ogs] [[32minfo[m] Nr. Nodes: 34959.
[2024-10-17 11:58:42.824] [ogs] [[32minfo[m] Nr. Elements: 191306.
[2024-10-17 11:58:42.825] [ogs] [[32minfo[m] Mem for mesh: 38 MiB
[2024-10-17 11:58:42.825] [ogs] [[32minfo[m] Time for reading: 2.042905 seconds.
[2024-10-17 11:58:42.825] [ogs] [[32minfo[m] Read 34959 nodes and 191306 elements.
[2024-10-17 11:58:42.825] [ogs] [[32minfo[m] Please check your mesh carefully!
[2024-10-17 11:58:42.825] [ogs] [[32minfo[m] Degenerated or redundant mesh elements can cause OGS to stop or misbehave.
[2024-10-17 11:58:42.825] [ogs] [[32minfo[m] Use the -e option to delete redundant line elements.
[2024-10-17 11:58:42.825] [ogs] [[32minfo[m] Axis aligned bounding box: 	x [-10, 10) (extent 20)
	y [0, 20) (extent 20)
	z [-12, 4.94066e-324) (extent 12)

[2024-10-17 11:58:42.833] [ogs] [[32minfo[m] Edge length: [0.0379059, 1.05802]
[2024-10-17 11:58:42.834] [ogs] [[32minfo[m] Number of elements in the mesh:
[2024-10-17 11:58:42.834] [ogs] [[32minfo[m] 	Lines: 30
[2024-10-17 11:58:42.834] [ogs] [[32minfo[m] 	Tetrahedrons: 190916
[2024-10-17 11:58:42.834] [ogs] [[32minfo[m] 	Pyramids: 180
[2024-10-17 11:58:42.834] [ogs] [[32minfo[m] 	Prisms: 180
[2024-10-17 11:58:42.834] [ogs] [[32minfo[m] There are 1 properties in the mesh:
[2024-10-17 11:58:42.834] [ogs] [[32minfo[m] 	MaterialIDs: (191306 values) [0, 3]
[2024-10-17 11:58:42.834] [ogs] [[32minfo[m] Mesh Quality Control:
[2024-10-17 11:58:42.834] [ogs] [[32minfo[m] Looking for unused nodes...
[2024-10-17 11:58:42.846] [ogs] [[32minfo[m] Found 0 potentially collapsible nodes.
[2024-10-17 11:58:42.857] [ogs] [[32minfo[m] Testing mesh element geometry:

[2024-10-17 11:58:42.896] [ogs] [[32minfo[m] No errors found.
[2024-10-17 11:58:42.896] [ogs] [[32minfo[m] No elements found with zero volume.

[2024-10-17 11:58:42.896] [ogs] [[32minfo[m] No elements found with non coplanar nodes.

[2024-10-17 11:58:42.896] [ogs] [[32minfo[m] No elements found with non-convex geometry.

[2024-10-17 11:58:42.896] [ogs] [[32minfo[m] No elements found with wrong node order.

[2024-10-17 11:58:42.914] [ogs] [[32minfo[m] 1 property vectors copied, 0 vectors skipped.
[2024-10-17 11:58:42.915] [ogs] [[32minfo[m] No holes found within the mesh.

The above conversion tool also shows that there exist line, tetrahedrons, pyramids and prism elements.The number of lines, tetrahedrons, pyramids and prism element is respectively 30, 190918, 180 and 180 with no error.

check_file = Path(f"{out_dir}/{bhe_mesh_file_name}.vtu").is_file()… (click to toggle)

check_file = Path(f"{out_dir}/{bhe_mesh_file_name}.vtu").is_file()
if check_file:
    print("Conversion of mesh file from Gmsh to VTU format was successful.")
else:
    print(
        "Error! Gmsh file is not properly converted to VTU format in the BHE meshing tutorial."
    )
    raise SystemExit()

Conversion of mesh file from Gmsh to VTU format was successful.