OpenGeoSys

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.

inclined_bhe_1D.png

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)
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)
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)
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)
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).

Inclined_bhe_top_view.png

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)
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.00134748s, CPU 0.002237s)
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.147829s, CPU 0.1477s)
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 8828 nodes...

Info    : Done tetrahedrizing 8836 nodes (Wall 0.0757951s, CPU 0.070336s)
Info    : Reconstructing mesh...
Info    :  - Creating surface mesh
Info    :  - Identifying boundary edges
Info    :  - Recovering boundary
Info    :  - Added 1 Steiner point
Info    : Done reconstructing mesh (Wall 0.208402s, CPU 0.176202s)

Info    : Found volume 4
Info    : It. 0 - 0 nodes created - worst tet radius 11.0009 (nodes removed 0 0)
Info    : It. 500 - 492 nodes created - worst tet radius 3.3295 (nodes removed 0 8)
Info    : It. 1000 - 992 nodes created - worst tet radius 2.72916 (nodes removed 0 8)
Info    : It. 1500 - 1492 nodes created - worst tet radius 2.4216 (nodes removed 0 8)
Info    : It. 2000 - 1990 nodes created - worst tet radius 2.23048 (nodes removed 0 10)
Info    : It. 2500 - 2490 nodes created - worst tet radius 2.08119 (nodes removed 0 10)

Info    : It. 3000 - 2988 nodes created - worst tet radius 1.9662 (nodes removed 0 12)
Info    : It. 3500 - 3484 nodes created - worst tet radius 1.87854 (nodes removed 0 16)
Info    : It. 4000 - 3978 nodes created - worst tet radius 1.80735 (nodes removed 0 22)
Info    : It. 4500 - 4454 nodes created - worst tet radius 1.74434 (nodes removed 0 46)
Info    : It. 5000 - 4924 nodes created - worst tet radius 1.69214 (nodes removed 0 76)
Info    : It. 5500 - 5384 nodes created - worst tet radius 1.64791 (nodes removed 0 116)
Info    : It. 6000 - 5843 nodes created - worst tet radius 1.60634 (nodes removed 1 156)
Info    : It. 6500 - 6271 nodes created - worst tet radius 1.57395 (nodes removed 7 222)
Info    : It. 7000 - 6689 nodes created - worst tet radius 1.54215 (nodes removed 9 302)
Info    : It. 7500 - 7134 nodes created - worst tet radius 1.51052 (nodes removed 10 356)
Info    : It. 8000 - 7583 nodes created - worst tet radius 1.47978 (nodes removed 13 404)

Info    : It. 8500 - 8045 nodes created - worst tet radius 1.45186 (nodes removed 13 442)
Info    : It. 9000 - 8536 nodes created - worst tet radius 1.42602 (nodes removed 17 447)
Info    : It. 9500 - 9026 nodes created - worst tet radius 1.39904 (nodes removed 18 456)
Info    : It. 10000 - 9517 nodes created - worst tet radius 1.37516 (nodes removed 18 465)
Info    : It. 10500 - 10006 nodes created - worst tet radius 1.35324 (nodes removed 18 476)
Info    : It. 11000 - 10498 nodes created - worst tet radius 1.33302 (nodes removed 22 480)
Info    : It. 11500 - 10991 nodes created - worst tet radius 1.31569 (nodes removed 24 485)
Info    : It. 12000 - 11481 nodes created - worst tet radius 1.29775 (nodes removed 24 495)
Info    : It. 12500 - 11970 nodes created - worst tet radius 1.28033 (nodes removed 26 504)
Info    : It. 13000 - 12455 nodes created - worst tet radius 1.2657 (nodes removed 28 517)

Info    : It. 13500 - 12945 nodes created - worst tet radius 1.25057 (nodes removed 29 526)
Info    : It. 14000 - 13422 nodes created - worst tet radius 1.23703 (nodes removed 34 544)
Info    : It. 14500 - 13894 nodes created - worst tet radius 1.22327 (nodes removed 35 571)
Info    : It. 15000 - 14369 nodes created - worst tet radius 1.21004 (nodes removed 37 594)
Info    : It. 15500 - 14857 nodes created - worst tet radius 1.19728 (nodes removed 37 606)
Info    : It. 16000 - 15343 nodes created - worst tet radius 1.18476 (nodes removed 38 619)
Info    : It. 16500 - 15826 nodes created - worst tet radius 1.17404 (nodes removed 41 633)
Info    : It. 17000 - 16317 nodes created - worst tet radius 1.16245 (nodes removed 41 642)
Info    : It. 17500 - 16801 nodes created - worst tet radius 1.76256 (nodes removed 45 654)
Info    : It. 18000 - 17263 nodes created - worst tet radius 1.14079 (nodes removed 47 690)
Info    : It. 18500 - 17754 nodes created - worst tet radius 1.13005 (nodes removed 48 698)

Info    : It. 19000 - 18243 nodes created - worst tet radius 1.12024 (nodes removed 49 708)
Info    : It. 19500 - 18732 nodes created - worst tet radius 1.1107 (nodes removed 50 718)
Info    : It. 20000 - 19215 nodes created - worst tet radius 1.10161 (nodes removed 55 730)
Info    : It. 20500 - 19695 nodes created - worst tet radius 1.09284 (nodes removed 55 750)
Info    : It. 21000 - 20185 nodes created - worst tet radius 1.0845 (nodes removed 55 760)
Info    : It. 21500 - 20671 nodes created - worst tet radius 1.07675 (nodes removed 55 774)
Info    : It. 22000 - 21159 nodes created - worst tet radius 1.06917 (nodes removed 55 786)
Info    : It. 22500 - 21638 nodes created - worst tet radius 1.06147 (nodes removed 55 807)
Info    : It. 23000 - 22115 nodes created - worst tet radius 1.05439 (nodes removed 55 830)
Info    : It. 23500 - 22606 nodes created - worst tet radius 1.04733 (nodes removed 56 838)

Info    : It. 24000 - 23095 nodes created - worst tet radius 1.03991 (nodes removed 56 849)
Info    : It. 24500 - 23582 nodes created - worst tet radius 1.03306 (nodes removed 59 859)
Info    : It. 25000 - 24052 nodes created - worst tet radius 1.02655 (nodes removed 59 889)
Info    : It. 25500 - 24542 nodes created - worst tet radius 1.0201 (nodes removed 59 899)
Info    : It. 26000 - 25036 nodes created - worst tet radius 1.01311 (nodes removed 59 905)
Info    : It. 26500 - 25528 nodes created - worst tet radius 1.00714 (nodes removed 61 911)
Info    : It. 27000 - 26011 nodes created - worst tet radius 1.26426 (nodes removed 65 924)
Info    : 3D refinement terminated (34966 nodes total):
Info    :  - 293 Delaunay cavities modified for star shapeness
Info    :  - 991 nodes could not be inserted
Info    :  - 194580 tetrahedra created in 1.19667 sec. (162601 tets/s)
Info    : 270 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

Info    :  detmin =  1.218811746621236E-02 eps=  1.000000000000000E-06 65 iter term 1
Info    :  detmin =  1.218817896100575E-02 eps=  1.000000000000000E-06 1 iter term 1
Info    : Done optimizing pyramids for hybrid mesh

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

Info    : Optimization starts (volume = 4779.95) with worst = 0.00180778 / average = 0.781798:
Info    : 0.00 < quality < 0.10 :       403 elements
Info    : 0.10 < quality < 0.20 :      1284 elements
Info    : 0.20 < quality < 0.30 :      2154 elements
Info    : 0.30 < quality < 0.40 :      3489 elements
Info    : 0.40 < quality < 0.50 :      5543 elements
Info    : 0.50 < quality < 0.60 :      9026 elements
Info    : 0.60 < quality < 0.70 :     18281 elements
Info    : 0.70 < quality < 0.80 :     42696 elements
Info    : 0.80 < quality < 0.90 :     73446 elements
Info    : 0.90 < quality < 1.00 :     38253 elements
Info    : 3753 edge swaps, 0 node relocations (volume = 4779.95): worst = 0.0795196 / average = 0.793812 (Wall 0.104553s, CPU 0.104654s)
Info    : 3827 edge swaps, 0 node relocations (volume = 4779.95): worst = 0.0795196 / average = 0.793967 (Wall 0.132366s, CPU 0.132431s)
Info    : No ill-shaped tets in the mesh :-)
Info    : 0.00 < quality < 0.10 :         2 elements
Info    : 0.10 < quality < 0.20 :        12 elements
Info    : 0.20 < quality < 0.30 :        46 elements
Info    : 0.30 < quality < 0.40 :      3410 elements
Info    : 0.40 < quality < 0.50 :      5326 elements
Info    : 0.50 < quality < 0.60 :      8754 elements
Info    : 0.60 < quality < 0.70 :     18185 elements
Info    : 0.70 < quality < 0.80 :     43382 elements
Info    : 0.80 < quality < 0.90 :     74073 elements
Info    : 0.90 < quality < 1.00 :     37955 elements
Info    : Done optimizing mesh (Wall 0.456397s, CPU 0.450133s)
Info    : 35025 nodes 209389 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/vZ6vnZiU/0/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/vZ6vnZiU/0/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.

inclined3DBHE.png

2D version will look like this. inclined2DBHE.png

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-12-20 13:28:33.886] [ogs] [info] Reading /var/lib/gitlab-runner/builds/vZ6vnZiU/0/ogs/build/release-all/web/content/docs/tutorials/Inclined_bhe_meshing/notebook-inclined_bhe_meshing/inclined_bhe_mesh_file.msh.
[2024-12-20 13:28:33.886] [ogs] [warning] 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-12-20 13:28:35.918] [ogs] [info] 	... finished.
[2024-12-20 13:28:35.918] [ogs] [info] Nr. Nodes: 34992.
[2024-12-20 13:28:35.918] [ogs] [info] Nr. Elements: 191577.
[2024-12-20 13:28:35.919] [ogs] [info] Mem for mesh: 38 MiB
[2024-12-20 13:28:35.919] [ogs] [info] Time for reading: 2.033110 seconds.
[2024-12-20 13:28:35.919] [ogs] [info] Read 34992 nodes and 191577 elements.
[2024-12-20 13:28:35.919] [ogs] [info] Please check your mesh carefully!
[2024-12-20 13:28:35.919] [ogs] [info] Degenerated or redundant mesh elements can cause OGS to stop or misbehave.
[2024-12-20 13:28:35.919] [ogs] [info] Use the -e option to delete redundant line elements.
[2024-12-20 13:28:35.920] [ogs] [info] Axis aligned bounding box: 	x [-10, 10) (extent 20)
	y [0, 20) (extent 20)
	z [-12, 4.94066e-324) (extent 12)
[2024-12-20 13:28:35.929] [ogs] [info] Edge length: [0.211602, 1.0749]
[2024-12-20 13:28:35.931] [ogs] [info] Number of elements in the mesh:
[2024-12-20 13:28:35.931] [ogs] [info] 	Lines: 30
[2024-12-20 13:28:35.931] [ogs] [info] 	Tetrahedrons: 191187
[2024-12-20 13:28:35.931] [ogs] [info] 	Pyramids: 180
[2024-12-20 13:28:35.931] [ogs] [info] 	Prisms: 180
[2024-12-20 13:28:35.931] [ogs] [info] There are 1 properties in the mesh:
[2024-12-20 13:28:35.931] [ogs] [info] 	MaterialIDs: (191577 values) [0, 3]
[2024-12-20 13:28:35.931] [ogs] [info] Mesh Quality Control:
[2024-12-20 13:28:35.931] [ogs] [info] Looking for unused nodes...

[2024-12-20 13:28:35.946] [ogs] [info] Found 0 potentially collapsible nodes.
[2024-12-20 13:28:35.961] [ogs] [info] Testing mesh element geometry:

[2024-12-20 13:28:36.003] [ogs] [info] No errors found.
[2024-12-20 13:28:36.003] [ogs] [info] No elements found with zero volume.

[2024-12-20 13:28:36.003] [ogs] [info] No elements found with non coplanar nodes.

[2024-12-20 13:28:36.003] [ogs] [info] No elements found with non-convex geometry.

[2024-12-20 13:28:36.003] [ogs] [info] No elements found with wrong node order.

[2024-12-20 13:28:36.027] [ogs] [info] 1 property vectors copied, 0 vectors skipped.
[2024-12-20 13:28:36.028] [ogs] [info] 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.
Previous
Mesh Tutorial
Next
Creating Mesh with integrated Borehole Heat Exchangers

This article was written by Joy Brato Shil, Haibing Shao. If you are missing something or you find an error please let us know.
Generated with Hugo 0.122.0 in CI job 504124 | Last revision: November 27, 2023