Meshing problem using OpenCascadeLink

Question

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["NDSolve`FEM`"]
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];
facesMemb = OpenCascadeShapeFaces[sMemb];

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*)
serpentineUnion = OpenCascadeShapeUnion[Flatten@sSerpentine];
facesSerpentine = OpenCascadeShapeFaces[serpentineUnion];

(*Union the faces of the two objects to preserve internal boundaries*)
union = OpenCascadeShapeUnion[Flatten[{facesSerpentine, facesMemb}]];

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.

Answer 1

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["NDSolve`FEM`"]
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.