I have encountered a meshing problem and I am a bit stuck with it. I tried to replicate the problem while removing most of the original code to be more concise.

I am using OpenCascadeLink to perform 3D meshing as it is very convenient and works well. The design consists of a 3D serpentine shaped that is in contact with a very large cuboid placed below. Below is an image of the top view of the simplified design so it is easier to understand.

The serpentine shape has some thickness in the z direction tt and the cuboid that is placed below as well (its thickness is tm). The serpentine is build by stitching together several cuboids. The small dashed square of width wl is only represented to build the serpentine shape, but is not part of the mesh.

The code below builds the shape and shows it roughly using Graphics3D:

Needs["NDSolveFEM"]
Needs["OpenCascadeLink"]

(*Chip thickness*)
tc = 250. 10^-6;

(*Large square width*)
wm = 0.5 10^-2;
(*Large square length*)
lm = wm;
(*Large square thickness*)
tm = 10. 10^-6 ;

(*Square width (x)*)
wl = wm/2;
(*Square length (y)*)
ll = wl;

(*Serpentine width*)
wt = 12. 10^-6;
(*Serpentine thickness*)
tt = 10. 10^-6;
(*Serpentine shape parameters*)
a = wl/6.;
b = 0.9*(ll);

(*large square design*)

memb = Cuboid[{-wm/2, -wm/2, tc}, {wm/2, wm/2, tc + tm}];

(*Serpentine design*)

serpentineLinesVertical = Table[{
Cuboid[{(-wl/2 + a/2) - wt/2, -b/2 - wt/2,
tc + tm}, {(-wl/2 + a/2) + wt/2, wt/2, tc + tm + tt}] /. {x_,
y_, z_} -> {x + n a, y, z}}, {n, {0, 5}}];
serpentineLinesVertical2 = Table[{
Cuboid[{(-wl/2 + a/2) + a - wt/2, -b/2 - wt/2,
tc + tm}, {(-wl/2 + a/2) + a + wt/2, b/2 + wt/2,
tc + tm + tt}] /. {x_, y_, z_} -> {x + n a, y, z}}, {n, 0, 3,
1}];
serpentineLinesHorizontal =
Table[{Cuboid[{-wl/2, wt/2, tc + tm}, {-wl/2 + a/2 + wt/2, -wt/2,
tc + tm + tt}] /. {x_, y_, z_} -> {x + n a, y, z}}, {n, {0}}];
serpentineLinesHorizontal2 =
Table[{Cuboid[{-wl/2 - wt/2, wt/2, tc + tm}, {-wl/2 + a/2, -wt/2,
tc + tm + tt}] /. {x_, y_, z_} -> {x + n a, y, z}}, {n, {11/
2}}];
serpentineLinesHorizontal3 =
Table[{Cuboid[{-wl/2 + 3 a/2 - wt/2, b/2 + wt/2,
tc + tm}, {-wl/2 + 5 a/2 + wt/2, b/2 - wt/2,
tc + tm + tt}] /. {x_, y_, z_} -> {x + n a, y, z}}, {n, {0,
2}}];
serpentineLinesHorizontal4 =
Table[{Cuboid[{-wl/2 + a/2 - wt/2, -b/2 + wt/2,
tc + tm}, {-wl/2 + 3 a/2 + wt/2, -b/2 - wt/2,
tc + tm + tt}] /. {x_, y_, z_} -> {x + n a, y, z}}, {n, {0, 2,
4}}];
partsSerpentine =
Flatten[Flatten[#] & /@ {serpentineLinesVertical,
serpentineLinesVertical2, serpentineLinesHorizontal,
serpentineLinesHorizontal2, serpentineLinesHorizontal3,
serpentineLinesHorizontal4}];

Graphics3D[{memb, partsSerpentine}, Boxed -> False]


This looks OK, and the different parts of the serpentine overlap to later form an union.

Now I try to build it using OpenCascadeLink. I want to remove the internal boundaries within the serpentine shape, but I want to keep them between the serpentine and the cuboid placed below. Therefore, following the tutorial, I split the faces and then put them back together to keep internal boundaries between the cuboid and the serpentine.

First dealing with the cuboid placed below the serpentine:

sMemb = OpenCascadeShape[memb];


Then dealing with the serpentine :

sSerpentine = OpenCascadeShape /@ partsSerpentine;
(*I first Union the shapes, and then create the faces because I want to get rid of internal boundaries*)

(*Union the faces of the two objects to preserve internal boundaries*)


I can now create and visualize the boundary mesh:

bmesh = OpenCascadeShapeSurfaceMeshToBoundaryMesh[union];
Show[bmesh["Edgeframe"], ViewPoint -> Above]


showing something that looks good.

In particular in this region :

Making the mesh

mesh = ToElementMesh[bmesh,
"RegionMarker" -> {
{(* 1 : Membrane*){0, 0, tc + tm/2}, 1, (wm*wm*tm)/2},
{(* 2 : Serpentine*){0, -b/2, tc + tm + tt/2}, 2, (tt^3)/3}}
]


leads to a mesh without problem. If I now want to scale down the serpentine part, I change the parameter wl, and set it to wl=wm/12 (its previous value was wl=wm/2).

For the Boundarymesh Edgeframe representation, zoomed on the serpentine part, it leads to :

and only the serpentine shape is meshed (not the cuboid placed below). If I use the options "SimplifyResult"->True  during the Union process of the serpentine parts

  OpenCascadeShapeUnion[Flatten@sSerpentine,
"SimplifyResult" -> True];


I get a different result for the Boundarymesh:

but no mesh is generated. Oddly, if I want to scale down the serpentine even more, by setting wl=wm/15, then the BoundaryMesh looks correct, and a proper mesh is generated.

I have tried playing with some OpenCascade functions such as OpenCascadeShapeFix or OpenCascadeSimplify, but I still do not understand why the mesh works in some cases and not in other cases. I also tried to Union the serpentine parts sequentially, but that did not work either. I feel that the way the serpentine shape is put together might not be correct and cause troubles, but I do not know how to do it differently (I also tried using RegionUnion but that did not work better). Any help would be appreciated. I am using Mathematica V13.2 and OpenCascade v 7.6.3.

• I took a look at this but I can not say what causes this. I have used "ShapeSurfaceMeshOptions" -> {"LinearDeflection" -> 0.0000001} but that did not help. Another curious finding is that using and ignoring internal boundarries union = OpenCascadeShapeUnion[ Flatten[{serpentineUnion, OpenCascadeShape[memb]}]] has the same problem. Then I tried to Rationalize the numbers for the geometry. Nope. I suspect this is an issue with OpenCascade. As an experiment try to submerge the serpentine shape into the box and see what that does. Ping me if you find anything new. Commented Jun 2, 2023 at 8:27
• Thank you very much for your response. I have tried what you mentioned with no success. Without the membrane (cuboid below) there is no problem though. I tried to play a bit with the serpentine parameters, for instance with wl=wm/12 and a = wl/5 instead of a = wl/6 it works fine. However, reducing the thickness and width of the serpentine to 0.5*10^-6 fails completely (even without the cuboid below). Could it be related to a tolerance issue ? Commented Jun 8, 2023 at 18:30

It seems this is a tolerance issue in/with OpenCascade. Then here is an idea for a workaround. The idea is to scale this into something where we do not meet the OpenCascade tolerance setting. For example multiply everything by, say, 10^3. Then you should be able to get a decent geometry and a boundary mesh. Then you can attribute a length scale to the mesh. Let's look at a simple example:

Needs["NDSolveFEM"]
bmesh = ToBoundaryMesh[Annulus[{1, 0}, {1, 2}]];

bmesh["LengthUnit"]
(* None *)

SetLengthUnit[bmesh, "Meters"]
(* "Meters" *)


You can the use ElementMeshCoordinateRescale to rescale the coordinates to what you need.

bmesh2 = ElementMeshCoordinateRescale[bmesh, "Milimeters"]


Keep in mind though, that then it might be numerically better to scale the PDE coefficients into something of the same scale. Have a look at the this section in the documentation that hints at how this could be done.

• Thank you very much, I will try this! Commented Jun 9, 2023 at 16:31
• I applied your suggestion and it works well for building a small geometry -- it solved many of the problems I had trying to model thin films using OpenCascade. Thank you! Is theScaleUnits option only available when defining material properties using Entity[]? For instance, if I use a "HeatTransferCoefficient"` value in a boundary condition, will the unit be scaled accordingly? Commented Jun 10, 2023 at 0:13
• @Dopseu, it should also work when you specify values with Quantities. Commented Jun 12, 2023 at 6:21
• I also added a section on tolerance issues in the using OpenCascadeLink tutorial. Commented Jun 12, 2023 at 6:55
• Got it, thanks. For now I cannot see the section on tolerance but perhaps it will take some time to update? Commented Jun 12, 2023 at 16:59