diff --git a/docs/examples/howto_meshlib/plot_bhe_mesh.py b/docs/examples/howto_meshlib/plot_bhe_mesh.py
new file mode 100644
index 0000000000000000000000000000000000000000..a956ba13ec012511776fa237ddfe4cfc026c57d1
--- /dev/null
+++ b/docs/examples/howto_meshlib/plot_bhe_mesh.py
@@ -0,0 +1,79 @@
+"""
+Creating a BHE mesh (Borehole Heat Exchanger)
+=============================================
+
+This example demonstrates how to create a Borehole Heat Exchanger (BHE) mesh.
+"""
+
+# %%
+from pathlib import Path
+from tempfile import mkdtemp
+
+import pyvista as pv
+from pyvista.plotting import Plotter
+
+from ogstools.meshlib.gmsh_meshing import bhe_mesh
+from ogstools.msh2vtu import msh2vtu
+
+# %% [markdown]
+# Generate a customizable BHE mesh (using gmsh):
+
+# %%
+tmp_dir = Path(mkdtemp())
+msh_file = tmp_dir / "bhe.msh"
+bhe_mesh(
+ width=20,
+ length=30,
+ depth=40,
+ x_BHE=10,
+ y_BHE=10,
+ bhe_depth=25,
+ out_name=msh_file,
+)
+
+# %% [markdown]
+# Now we convert the gmsh mesh to the VTU format with msh2vtu.
+# Passing the list of dimensions [1, 3] to msh2vtu ensures, that the line
+# elements will also be part of the domain mesh.
+
+# %%
+msh2vtu(
+ msh_file,
+ output_path=tmp_dir,
+ rdcd=True,
+ ogs=True,
+ dim=[1, 3],
+ log_level="ERROR",
+)
+
+# %% [markdown]
+# Load the domain mesh and extract BHE line:
+
+# %%
+mesh = pv.read(tmp_dir / "bhe_domain.vtu")
+bhe_line = mesh.extract_cells_by_type(pv.CellType.LINE)
+
+# %% [markdown]
+# Visualize the mesh:
+
+# %%
+p = Plotter()
+p.add_mesh(mesh, style="wireframe", color="grey")
+p.add_mesh(
+ mesh.clip("x", bhe_line.center, crinkle=True),
+ show_edges=True,
+ scalars="MaterialIDs",
+ cmap="Accent",
+ categories=True,
+ scalar_bar_args={"vertical": True, "n_labels": 2, "fmt": "%.0f"},
+)
+p.add_mesh(bhe_line, color="r", line_width=3)
+p.show()
+
+# sphinx_gallery_start_ignore
+
+from shutil import rmtree # noqa: E402 # pylint: disable=C0413, C0411
+
+rmtree(tmp_dir)
+
+# sphinx_gallery_end_ignore
diff --git a/ogstools/meshlib/gmsh_meshing.py b/ogstools/meshlib/gmsh_meshing.py
index 542a1e76a720a7eb26d853ad7dc446ed10d55ea8..51862b6852226dccfb7f7d29687a80c1e8160ab2 100644
--- a/ogstools/meshlib/gmsh_meshing.py
+++ b/ogstools/meshlib/gmsh_meshing.py
@@ -107,3 +107,68 @@ def cuboid(
gmsh.model.mesh.setOrder(order)
gmsh.write(str(out_name))
gmsh.finalize()
+
+
+def bhe_mesh(
+ width: float = 10.0,
+ length: float = 10.0,
+ depth: float = 30.0,
+ x_BHE: float = 5,
+ y_BHE: float = 5,
+ bhe_depth: float = 20,
+ order: int = 1,
+ out_name: Path = Path("bhe_mesh.msh"),
+):
+ gmsh.initialize()
+ model, geo = (gmsh.model, gmsh.model.geo)
+ model.add("Testmodel")
+
+ geo.addPoint(0.0, 0.0, 0.0)
+ geo.addPoint(width, 0.0, 0.0)
+ geo.addPoint(width, length, 0.0)
+ geo.addPoint(0.0, length, 0.0)
+
+ geo.addLine(1, 2)
+ geo.addLine(2, 3)
+ geo.addLine(3, 4)
+ geo.addLine(4, 1)
+
+ geo.addCurveLoop([1, 2, 3, 4], 1)
+ geo.addPlaneSurface([1], tag=1)
+
+ alpha = 6.134 # valid only for 6 surronding points
+ bhe_radius = 0.076
+ delta = alpha * bhe_radius
+ delta_2 = 0.866 * delta
+ delta_3 = 0.5 * delta
+
+ geo.addPoint(x_BHE, y_BHE, 0, delta)
+ d1 = geo.addPoint(x_BHE, y_BHE - delta, 0, delta)
+ d2 = geo.addPoint(x_BHE, y_BHE + delta, 0, delta)
+ d3 = geo.addPoint(x_BHE + delta_2, y_BHE + delta_3, 0, delta)
+ d4 = geo.addPoint(x_BHE - delta_2, y_BHE + delta_3, 0, delta)
+ d5 = geo.addPoint(x_BHE + delta_2, y_BHE - delta_3, 0, delta)
+ d6 = geo.addPoint(x_BHE - delta_2, y_BHE - delta_3, 0, delta)
+ geo.synchronize()
+ model.mesh.embed(0, [d1, d2, d3, d4, d5, d6], 2, 1)
+
+ soil_1 = geo.extrude([(2, 1)], 0, 0, -depth / 2, [6], [1], True)
+ soil_2 = geo.extrude([soil_1[0]], 0, 0, -depth / 2, [6], [1], True)
+ bhe = geo.extrude([(0, 5)], 0, 0, -bhe_depth, [7], [1], True)
+
+ top_soil_tag = model.addPhysicalGroup(dim=3, tags=[soil_1[1][1]])
+ bottom_soil_tag = model.addPhysicalGroup(dim=3, tags=[soil_2[1][1]])
+ bhe_tag = model.addPhysicalGroup(dim=1, tags=[bhe[1][1]])
+
+ model.setPhysicalName(dim=3, tag=top_soil_tag, name="top")
+ model.setPhysicalName(dim=3, tag=bottom_soil_tag, name="bottom")
+ model.setPhysicalName(dim=1, tag=bhe_tag, name="bhe")
+
+ geo.synchronize()
+ model.mesh.generate(3)
+ gmsh.option.setNumber("Mesh.SecondOrderIncomplete", 1)
+ model.mesh.setOrder(order)
+ model.mesh.removeDuplicateNodes()
+
+ gmsh.write(str(out_name))
+ gmsh.finalize()
diff --git a/ogstools/msh2vtu/__init__.py b/ogstools/msh2vtu/__init__.py
index db5df401181a35a6d0f42e44d774173f056d7019..aaeb297945220c98a4ebf19cbe55b8b231c3b430 100644
--- a/ogstools/msh2vtu/__init__.py
+++ b/ogstools/msh2vtu/__init__.py
@@ -119,11 +119,8 @@ def find_connected_domain_cells(
# there should be one domain cell for each boundary cell, however cells
# of boundary dimension may be in the domain (e.g. as sources)
if number_of_connected_domain_cells == 1:
- domain_cells_array[
- cell_index
- ] = (
- common_domain_cells.pop()
- ) # assign only one (unique) connected dmain cell
+ # assign only one (unique) connected dmain cell
+ domain_cells_array[cell_index] = common_domain_cells.pop()
elif number_of_connected_domain_cells < 1 and not warned_lt1:
logging.warning(
"find_connected_domain_cells: cell %s"
@@ -156,7 +153,7 @@ def msh2vtu(
input_filename: Path,
output_path: Path = Path(),
output_prefix: str = "",
- dim: int = 0,
+ dim: Union[int, list[int]] = 0,
delz: bool = False,
swapxy: bool = False,
rdcd: bool = True,
@@ -253,6 +250,7 @@ def msh2vtu(
# set spatial dimension of mesh
if dim == 0:
+ assert isinstance(dim, int)
# automatically detect spatial dimension of mesh
_dim = dim0 # initial value
for test_dim, test_cell_types in available_cell_types.items():
@@ -265,7 +263,8 @@ def msh2vtu(
logging.info("Detected mesh dimension: %s", str(_dim))
logging.info("##")
else:
- _dim = dim # trust the user
+ # trust the user
+ _dim = max(dim) if isinstance(dim, list) else dim
# delete third dimension if wanted by user
if delz:
@@ -292,6 +291,8 @@ def msh2vtu(
)
domain_cell_types = existing_cell_types.intersection(
available_cell_types[domain_dim]
+ if isinstance(dim, int)
+ else set().union(*[available_cell_types[d] for d in dim])
)
else:
logging.warning("Error, invalid dimension dim=%s!", str(_dim))
@@ -639,7 +640,13 @@ def msh2vtu(
for cell_type in subdomain_cell_types:
# cells
- selection_index = mesh.cell_sets_dict[name][cell_type]
+ all_false = np.full(
+ cell_data_dict[gmsh_physical_cell_data_key][cell_type].shape,
+ False,
+ )
+ selection_index = mesh.cell_sets_dict[name].get(
+ cell_type, all_false
+ )
selection_cells_values = cells_dict[cell_type][selection_index]
if len(selection_cells_values): # if there are some data
selection_cells_block = (cell_type, selection_cells_values)
diff --git a/ogstools/msh2vtu/_cli.py b/ogstools/msh2vtu/_cli.py
index 36917f79875fa76208a963c904f7602b818b96e4..495b1a0b361a45c281b37c773d8e1b09c87a87bf 100644
--- a/ogstools/msh2vtu/_cli.py
+++ b/ogstools/msh2vtu/_cli.py
@@ -43,10 +43,12 @@ def argparser():
"-d",
"--dim",
type=int,
+ nargs="*",
default=0,
help=(
"spatial dimension (1, 2 or 3), trying automatic detection, "
- "if not given"
+ "if not given. If multiple dimensions are provided, all elements of"
+ "these dimensions are embedded in the resulting domain mesh."
),
)
parser.add_argument(
diff --git a/tests/test_msh2vtu.py b/tests/test_msh2vtu.py
index a69c0f6c62e242642b6f4e9721a3b2d3fc9431e4..f35deb74920bd9e1e7b95db75177d68e60da09af 100644
--- a/tests/test_msh2vtu.py
+++ b/tests/test_msh2vtu.py
@@ -1,6 +1,7 @@
"""
Tests (pytest) for msh2vtu
"""
+
import os
import runpy
import subprocess
@@ -115,6 +116,28 @@ def test_cuboid(tmp_path: Path):
assert msh2vtu(msh_file, tmp_path, output_prefix="cuboid") == 0
+def test_bhe_mesh(tmp_path: Path):
+ """Create bhe gmsh mesh and convert with msh2vtu."""
+ msh_file = Path(tmp_path, "bhe.msh")
+ permutations = product(
+ [10.0, 20.0], [15.0, 30.0], [40.0, 80.0], [20.0, 30.0]
+ )
+ for width, length, depth, bhe_depth in permutations:
+ gmsh_meshing.bhe_mesh(
+ width=width,
+ length=length,
+ depth=depth,
+ x_BHE=5.0,
+ y_BHE=5.0,
+ bhe_depth=bhe_depth,
+ order=1,
+ out_name=msh_file,
+ )
+ assert msh2vtu(msh_file, rdcd=True, ogs=True, dim=[1, 3]) == 0
+ mesh = pv.read("bhe_domain.vtu")
+ assert set(mesh.celltypes) == {3, 13} # wedges (3D) and lines (1D)
+
+
def test_gmsh(tmp_path: Path):
os.chdir(tmp_path)
for script in [