[msh2vtu] fix for BHE meshes

- allow to embed multiple dims in domain (3 and 1 for BHE meshes)
- added a test for this
- added a docu example

docs/examples/howto_meshlib/plot_bhe_mesh.py

+"""
+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

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

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)

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(

tests/test_msh2vtu.py

 """
 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 [