[pythag0ra5]

Hi everyone,

I would like to get your input on how to mesh a geometry that consists of a large surrounding domain containing a high number of small cylindrical holes. The simplified example below shows 19 holes, but in my actual case the number will easily exceed 100. The geometry and all dimensions are fully defined, and the holes follow a regular pattern.



For several reasons I cannot apply periodic boundary conditions, so I need to simulate the full perforated region with all holes explicitly included.

I have tried using snappyHexMesh, which works in principle, but it produces extremely large meshes. In addition, I cannot enforce clean O-grids or purely hexahedral meshes, which I would strongly prefer.

I am therefore looking for alternative strategies for meshing this type of geometry. One idea I considered is to create a blockMesh region around a single cylinder and then replicate or "tile" this block multiple times. I am not sure, however, whether this is a practical approach.

During my research I found an example on the CFMesh website that is very close to my case, but unfortunately it refers to the commercial version, which I don’t have access to.

I would be very interested in hearing about your experiences or suggestions for similar geometries. In particular, I am curious whether this could be accomplished using pure blockMesh, possibly supported by Python scripting. I also came across the classy-blocks toolbox, which looks promising, but I have no hands-on experience with it yet.

Any suggestions or insights would be greatly appreciated!

[Yann]

Hi!

Originally Posted by pythag0ra5

I am therefore looking for alternative strategies for meshing this type of geometry. One idea I considered is to create a blockMesh region around a single cylinder and then replicate or "tile" this block multiple times. I am not sure, however, whether this is a practical approach.

Never tried this myself, but maybe you could give a try to mirrorMesh?

[pythag0ra5]

Hi Yann,

thanks a lot for your reply. Yes, using mirrorMesh could be an approach. However, as I understand it, this would mainly help with mirroring the left part of the domain onto the right. Definitely a simplification, but it still doesn’t address the core question.

What I’m still unsure about is how to mesh one section of the block and then replicate this meshed pattern many times before applying the mirroring. Do you have any idea how this could be done?

[kandelabr]

Hello,


yes, as per your inquiry, this is doable with classy_blocks; a tileable 2D sketch, connecting a cylinder and a ring to a square or a hexagon would be the base; then a simple grid would cover the most of the domain.


I can prepare something quick this week and let you know so you can tell how close it is to what you need.

[GerhardHolzinger]

Some general points about meshing with blockMesh:

• Can you draw it using only rectangles?
• Can you identify a repeating pattern in your geometry?

The first point is an exercise in creativity. This makes or brakes your attempt of using blockMesh.


The second point can save you tons of work. Re-use is the name of the game. You brought up class_blocks in your question: this is good thinking. Automation-by-scripting, which implies clever use of data structures, turns a potential slog of tedious work into something practical.


Even when I don't use Python in conjunction with blockMesh, for very simple geometries, I use OpenFOAM's inline calculations.


A final note: don't fixate on blockMesh only. Using blockMesh in combination with subsetMesh and mirrorMesh may unlock otherwise impossible or harder-to-achieve features. Throw in extrudeMesh as well as combineMesh/stitchMesh for good measure. The attached mesh is such a case. blockMesh-only, this would be quite hard.


Another final note: if you geometry turns out to be too much do deal with in a blockMesh-based meshing strategy; you could carry a lot of the considerations over to a Salome-based meshing strategy. In my experience, following a bottom-up meshing approach in Salome is also quite good.

[kandelabr]

This is what I can whip up with ~130 lines and a couple of hours.
It's fully parametric but it will still take some extra elbow grease to fill the gaps to get a rectangular domain. Also, patches and cell zones are missing (but it's pretty easy to add those).

See the attached script. To keep things as simple as possible, I allowed myself to utilize some greater classy kung fu but I think it should be pretty self-explanatory.

[GerhardHolzinger]

I have no experience with classy_blocks, but if you can do two more additional types blocks, then you can make the transition to a rectangular overall domain.


Orange is basically a half-hexagon. Grey is the other block needed to fill in the gaps.

[kandelabr]

Sure, it's just a matter of sacrificing more of my time to the gods of meshing. This is what I could do before lunch

[GerhardHolzinger]

I was thinking more about the person who initially posted the question.


You, kandelabr, do sterling service, when you provide random strangers on the internet, vital building blocks upon which they can build their work.

[pythag0ra5]

Dear kandelabr,
dear GerhardHolzinger,

thank you for your replies. Please excuse my late answer. I was ill this week and only now had the time to look at your comments in detail.

• Can you draw it using only rectangles?
• Can you identify a repeating pattern in your geometry?

Yes. In this case both conditions are met. I guess this setup is quite typical for a heat-exchanger-type arrangement.

A final note: don't fixate on blockMesh only. Using blockMesh in combination with subsetMesh and mirrorMesh may unlock otherwise impossible or harder-to-achieve features.

Could you give me some guidance on how subsetMesh could be applied here in a useful way? I have not used it in practice...

I agree with your general point. I also built my own workflow for generating blockMeshDicts and meshes in general. By now I construct many meshes directly in Python and do not rely on classyBlocks (altough this seems like a really interessting approach to me!).

For your reference, here is my current approach. It generates a rectangle of width W and height H with a circular cutout of diameter D. The center coordinates CX0 and CY0 control the circle’s position. This works reliably and is fully parameterized.

#!/usr/bin/env python3
# -*- coding: utf-8 -*-

"""
blockMeshDict generator
Rectangular 2D domain (W x H x T) with a circular cut-out (diameter D),
centered at (CX0, CY0). Cross topology (4 hex blocks). Circle rotated by 45°
so that each block owns exactly one quarter-arc as a true block edge.

Patches:
- cylinder              : inner wall (4 curved quads)
- frontAndBack          : empty (enforces 2D; requires NZ = 1)
- left|right|bottom|top : explicit outer walls

If you need real 3D change 'frontAndBack' to 'wall'.
"""

import math
import os

# ----------------------------
# User parameters
# ----------------------------
W   = 0.05      # domain width  [m]
H   = 0.05      # domain height [m]
T   = 0.01      # thickness (z) [m]

D   = 0.02      # hole diameter [m]
CX0 = 2.0       # hole center x [m]
CY0 = 1.0       # hole center y [m]

# base cell counts per block (x, y, z)
NX, NY, NZ = 20, 20, 2

# radial grading factor (>1 = finer at circle, coarser to the outside)
G_RAD = 5.0

# ----------------------------
# Derived geometry
# ----------------------------
R   = D / 2.0
X0  = CX0 - W/2.0          # domain min x
X1  = CX0 + W/2.0          # domain max x
Y0  = CY0 - H/2.0          # domain min y
Y1  = CY0 + H/2.0          # domain max y
Z1  = -T/2.0               # bottom z
Z2  =  T/2.0               # top z

# 45° offset for circle vertices (endpoints at 45°, 135°, 225°, 315°)
A = R / math.sqrt(2.0)

def P(x: float, y: float, z: float) -> str:
    """Format a point for OpenFOAM."""
    return f"({x:.6f} {y:.6f} {z:.6f})"

# ----------------------------
# Vertices (fixed indices)
# 00..03 domain bottom: 0=LL, 1=LR, 2=UR, 3=UL
# 04..07 domain top   : 4=LL, 5=LR, 6=UR, 7=UL
# 08..11 circle bottom  (45°, 135°, 225°, 315°)
# 12..15 circle top     (45°, 135°, 225°, 315°)
# ----------------------------
verts = [
    # domain bottom
    P(X0, Y0, Z1),  # 0
    P(X1, Y0, Z1),  # 1
    P(X1, Y1, Z1),  # 2
    P(X0, Y1, Z1),  # 3
    # domain top
    P(X0, Y0, Z2),  # 4
    P(X1, Y0, Z2),  # 5
    P(X1, Y1, Z2),  # 6
    P(X0, Y1, Z2),  # 7
    # circle bottom (rotated +45°)
    P(CX0 + A, CY0 + A, Z1),  # 8
    P(CX0 - A, CY0 + A, Z1),  # 9
    P(CX0 - A, CY0 - A, Z1),  # 10
    P(CX0 + A, CY0 - A, Z1),  # 11
    # circle top
    P(CX0 + A, CY0 + A, Z2),  # 12
    P(CX0 - A, CY0 + A, Z2),  # 13
    P(CX0 - A, CY0 - A, Z2),  # 14
    P(CX0 + A, CY0 - A, Z2),  # 15
]

# ----------------------------
# Blocks (lower 4 CCW, then upper 4 CCW)
# Quadrant mapping (outer corners in brackets):
#   NE (right/top)   [1,2]  carries arc 11→8
#   NW (top/left)    [2,3]  carries arc 8→9
#   SW (left/bottom) [3,0]  carries arc 9→10
#   SE (bottom/right)[0,1]  carries arc 10→11
# ----------------------------
block_NE = (11, 1, 2, 8,   15, 5, 6, 12)
block_NW = (8,  2, 3, 9,   12, 6, 7, 13)
block_SW = (9,  3, 0, 10,  13, 7, 4, 14)
block_SE = (10, 0, 1, 11,  14, 4, 5, 15)

# ----------------------------
# Arc edges: quarter circles (midpoints on cardinal axes)
# ----------------------------
arcs_bottom = [
    (8,  9,  P(CX0,     CY0 + R, Z1)),  # 45°→135°
    (9,  10, P(CX0 - R, CY0,     Z1)),  # 135°→225°
    (10, 11, P(CX0,     CY0 - R, Z1)),  # 225°→315°
    (11, 8,  P(CX0 + R, CY0,     Z1)),  # 315°→45°
]
arcs_top = [
    (12, 13, P(CX0,     CY0 + R, Z2)),
    (13, 14, P(CX0 - R, CY0,     Z2)),
    (14, 15, P(CX0,     CY0 - R, Z2)),
    (15, 12, P(CX0 + R, CY0,     Z2)),
]

# ----------------------------
# Front/back faces
# ----------------------------
front_back = [
    (11, 1, 2, 8),   (15, 5, 6, 12),   # NE
    (8,  2, 3, 9),   (12, 6, 7, 13),   # NW
    (9,  3, 0, 10),  (13, 7, 4, 14),   # SW
    (10, 0, 1, 11),  (14, 4, 5, 15),   # SE
]

# ----------------------------
# Cylinder faces (inner wall)
# ----------------------------
cyl_faces = [
    (11, 8,  12, 15),  # NE
    (8,  9,  13, 12),  # NW
    (9,  10, 14, 13),  # SW
    (10, 11, 15, 14),  # SE
]

# ----------------------------
# Outer boundary faces
# ----------------------------
left   = (0, 4, 7, 3)  # x = X0
right  = (1, 2, 6, 5)  # x = X1
bottom = (0, 1, 5, 4)  # y = Y0
top    = (3, 7, 6, 2)  # y = Y1

# ----------------------------
# Emit blockMeshDict
# ----------------------------
out = []

# Header
out.append(
    "FoamFile\n"
    "{\n"
    "    version     2.0;\n"
    "    format      ascii;\n"
    "    class       dictionary;\n"
    "    object      blockMeshDict;\n"
    "}\n\n"
    "scale 1;\n\n"
)

# Vertices
out.append("vertices\n(\n")
for v in verts:
    out.append(f"    {v}\n")
out.append(");\n\n")

# Blocks
out.append("blocks\n(\n")
out.append(f"    hex ({' '.join(map(str, block_NE))}) ({NX} {NY} {NZ}) simpleGrading ({G_RAD} 1 1)\n")
out.append(f"    hex ({' '.join(map(str, block_NW))}) ({NX} {NY} {NZ}) simpleGrading ({G_RAD} 1 1)\n")
out.append(f"    hex ({' '.join(map(str, block_SW))}) ({NX} {NY} {NZ}) simpleGrading ({G_RAD} 1 1)\n")
out.append(f"    hex ({' '.join(map(str, block_SE))}) ({NX} {NY} {NZ}) simpleGrading ({G_RAD} 1 1)\n")
out.append(");\n\n")

# Edges
out.append("edges\n(\n")
for i, j, mid in arcs_bottom:
    out.append(f"    arc {i} {j} {mid}\n")
for i, j, mid in arcs_top:
    out.append(f"    arc {i} {j} {mid}\n")
out.append(");\n\n")

# Boundary
out.append("boundary\n(\n")

# inner wall
out.append("    cylinder\n    {\n        type wall;\n        faces (\n")
for f4 in cyl_faces:
    out.append(f"            ({' '.join(map(str, f4))})\n")
out.append("        );\n    }\n\n")

# front/back = empty (2D)
out.append("    frontAndBack\n    {\n        type empty;\n        faces (\n")
for f4 in front_back:
    out.append(f"            ({' '.join(map(str, f4))})\n")
out.append("        );\n    }\n\n")

# explicit outer patches
out.append("    left\n    {\n        type wall;\n        faces (\n")
out.append(f"            ({' '.join(map(str, left))})\n")
out.append("        );\n    }\n\n")

out.append("    right\n    {\n        type wall;\n        faces (\n")
out.append(f"            ({' '.join(map(str, right))})\n")
out.append("        );\n    }\n\n")

out.append("    bottom\n    {\n        type wall;\n        faces (\n")
out.append(f"            ({' '.join(map(str, bottom))})\n")
out.append("        );\n    }\n\n")

out.append("    top\n    {\n        type wall;\n        faces (\n")
out.append(f"            ({' '.join(map(str, top))})\n")
out.append("        );\n    }\n")

out.append(");\n\n")

# No defaultPatch: everything is explicitly assigned
out.append("mergePatchPairs();\n")

# Write file
os.makedirs("system", exist_ok=True)
with open("system/blockMeshDict", "w") as f:
    f.writelines(out)

print("✔ Generated: system/blockMeshDict")

What I would like to extend is a complete grid of these "rectangle with circular cutout" blocks. I could expand my Python-script, but I am specifically interested in doing more of this with OpenFOAM’s native tools. One idea I considered is creating repeated tiles with translateMesh and connecting them with stitchMesh. I have no experience with this workflow, though.

Even if I created a stitched tile pattern as a single mesh, I still do not see a clean way to embed that pattern into a larger 2D mesh. I could script all outer boundaries in Python as well, but this seems impractical to me. It produces a long script tailored to one specific case and gives me few reusable building blocks for other projects.

A more intuitive strategy would be to generate

1. a large background mesh with blockMesh,
2. generate an x-by-y tile pattern that contains the circular cutouts,
3. and then merge or embed both meshes. snappyHexMesh might also be an option for this.

Thank you again for your suggestions. I would be glad to continue discussing these approaches with you.

Have a nice weekend.

P.S.: kandelabr, thanks a bunch for the code. I will only be able to go through it properly in the next few days. One quick question for clarification: You created the hexagonal tiling pattern to achieve the sequence of offset holes in each row, correct?

[kandelabr]

Hello,


yes, I chose a hexagonal pattern because that's the one you shared from cfMesh. If you need an orthogonal grid-like-pattern you'd only have to change the sketch and the rest will be greatly simplified.


I have played with my script a bit more - added the missing pieces at sides and top/bottom surfaces. I'll finish it this week and add it as an example into the main classy_blocks repository.


Will let you know (and you let me know if you need a square grid)!

[pythag0ra5]

Hello kandelabr,

Yes, I need a pattern where the rows alternate. One arrangement in the first row, a different one in the next.

In my case, however, the number of hexagons per row is constant. So the number of cells per row does not change. This afternoon I experimented further with my plain Python based approach. I ended up using several blocks per row (however, the same amount of them) and shifting the holes within each block by a fixed offset. That keeps an overall rectangular structure, but it makes the cells much more distorted than in your hexagonal layout.

That is why your hexagonal approach is very interesting for me. I am looking forward to your reply!

P.S.: I ended up coding everything in a single Python script, i.e. I have not used any OpenFOAM utilities such as translateMesh yet.

[kandelabr]

Well I guess you just need a different blocking pattern! Something like that maybe:

[pythag0ra5]

Hi Kandelabr,

I just wanted to give you a quick update from my side. And by the way, thanks a lot for your sketch, it helped me a lot!

At first I was actually trying to model each "tile" with four trapezoidal blocks. Only much later did it really sink in that this would cause problems, because the two rows are shifted by half the tile width. This shift means that a tile of width W will coincide with a strip / offset piece of width W/2.

The approach you showed in your sketch was exactly the right one: building a tile from eight trapezoidal blocks, so that each block is always congruent with the strip above or below it.

So my actual problem is solved now. For the sake of completeness, I’d still be curious to know whether this could be done more quickly or easily with classyBocks. My vibe-coding approach produced a Python-script with close to 400 lines of code. The code is fully parametric, meaning you can completely change the geometry of a single tile and also the geometry of the heat exchanger itself, i.e. the assembly of all tiles. That all works, but I have to say I would never have been able to develop all of that completely on my own from scratch.

In any case, it works now, so I’m basically done with my question. I’m just still curious how you would do this with classyBocks.

[kandelabr]

Great job!


The classy blocks script would look almost completely the same as the one for hexagonal layout, just the sketch would be a little different - and even that, only the positioning of points and quadrangle numbering.


I'll still add my version to the examples so you'll be able to compare.