Kernel crash during RegionDifference of complex 3D meshes + BoundaryMeshRegion::bsuncl error

Question

I'm modeling a screw rotor with cavities (for CAD/CFD analysis). The workflow involves:

Generating the rotor body from 5 parametric surfaces.
Creating cavity profiles from 4 extruded segments.
Subtracting cavities from the rotor using RegionDifference.

Problem:

With high PlotPoints: Kernel crashes silently during RegionDifference[bmr, bmr2].
With reduced PlotPoints (e.g., 30): BoundaryMeshRegion::bsuncl error (boundary not closed).

Error Messages:

BoundaryMeshRegion::bsuncl: The boundary surface is not closed because the edges [...] only come from a single face.

Code Summary:

Segment A


{X[1], Y[1]} = {R1 Cos [t], R1 Sin[t]};(*Inner Radius*)

dd = -ArcCos[Rp/R2] // N;
{X[2], Y[2]} = {-(2 Rp) Cos[dd*t] + R2 Cos [2 t*dd], -(2 Rp) Sin[t*dd] + 
    R2 Sin [2 t*dd]};(*Epicycloid Curve*)

{X[3], Y[3]} = {R2 Cos [t], R2 Sin[t]};(*Outer Rad*)

{r1} = {R2 - (R2 - 
        Rp)/(360 Degree - (Pi + γ - 10 Degree)) (t - (Pi + γ - 
         10 Degree))};(*Transition Curve A*)

{X[4], Y[4]} = {r1 Cos[t], r1 Sin[t]};

{r2} = {Rp - ((Rp - R1)/(β)) t};(*Transition Curve B*)

{X[5], Y[5]} = {r2 Cos[t], r2 Sin[t]};

{min1[1], max1[1]} = {β, Pi};
{min1[2], max1[2]} = {1, 0};
{min1[3], max1[3]} = {Pi, Pi + γ - 10 Degree};
{min1[4], max1[4]} = {Pi + γ - 10 Degree, 2 Pi};
{min1[5], max1[5]} = {0, β};

For[i = 1, i ⩽ 5, i++,
  {x1[i], y1[i], z1[i], dum} = M . {X[i], Y[i], 0, 1}; 
  ];

RotorL = {};
For[i = 1, i ⩽ 5, i++,
 {x1[i], y1[i], z1[i]};
 L1 = ParametricPlot3D[{x1[i], y1[i], z1[i]}, {t, min1[i], max1[i]}, {ϕr,
     0, wrap}, MaxRecursion -> 0, PlotPoints -> {100, 100}, 
   MeshFunctions -> {#4 &}, Mesh -> {{0, wrap}}, 
   MeshStyle -> Directive@{Thick, White}, Boxed -> False, Axes -> False, 
   PlotRange -> All, Method -> {"BoundaryOffset" -> False}, PlotRange -> All, 
   PlotStyle -> {White, Opacity[0.9]}, AxesLabel -> {"X", "Y", "Z"}];
 AppendTo[RotorL, L1];]

rb1 = 5;
db = 20;
hb = 3;
zx = 5;
ν = 2.25;
dc = 12.85;(*distance between cavity*)

(*cavity profile in the rotor body *)
ν0 = -ν - π/2;
zp = ν - π/2;(*cavity position in z axis*)
prof3Dr = {};
zpos = (L /. ϕr -> ϕrz) - 2;

φ = ((2 rb1 + dc)/Out1 )*3;
dot3r = ( {
     {Cos[ν0], -Sin[ν0], 0, 0},
     {Sin[ν0], Cos[ν0], 0, 0},
     {0, 0, 1, 0},
     {0, 0, 0, 1}
    } ) . ( {
     {Cos[ϕr], Sin[ϕr], 0, 0},
     {-Sin[ϕr], Cos[ϕr], 0, 0},
     {0, 0, 1, L},
     {0, 0, 0, 1}
    } ) . ( {
     {1, 0, 0, 0},
     {0, 1, 0, Out1},
     {0, 0, 1, zpos},
     {0, 0, 0, 1}
    } );

{x[1], z[1]} = {rb1*Cos[t], rb1*Sin[t]};
{x[2], z[2]} = {t Sin[π] - rb1, t Cos[π] + 1/2 hb};
{x[3], z[3]} = {rb1*Cos[t], rb1*Sin[t] + hb};
{x[4], z[4]} = {t Sin[π] + rb1, t Cos[π] + 1/2 hb};
{umin[1], umax[1]} = {2 π, π};
{umin[2], umax[2]} = {hb/2, -hb/2};
{umin[3], umax[3]} = {π, 0};
{umin[4], umax[4]} = {-hb/2, hb/2};

prof = {};
For[a = 1, a <= 4, a++,
 φ = ((2 rb1 + dc)/Out1 )*a;
 dot3r = ( {
     {Cos[ν0], -Sin[ν0], 0, 0},
     {Sin[ν0], Cos[ν0], 0, 0},
     {0, 0, 1, 0},
     {0, 0, 0, 1}
    } ) . ( {
     {Cos[ϕr], Sin[ϕr], 0, 0},
     {-Sin[ϕr], Cos[ϕr], 0, 0},
     {0, 0, 1, L},
     {0, 0, 0, 1}
    } ) . ( {
     {1, 0, 0, 0},
     {0, 1, 0, Out1},
     {0, 0, 1, zpos},
     {0, 0, 0, 1}
    } );
 For[i = 1, i ⩽ 4, i++,
  {X[i], Y[i], Z[i], dum} = 
    dot3r . {x[i], y, z[i], 1} /. {ϕr -> φ + 0, ϕrz -> 
       zp} ; 
  ];
 For[i = 1, i <= 4, i++,
  pro = ParametricPlot3D[{X[i], Y[i], Z[i]}, {t, umin[i], umax[i]}, {y, -20, 
     0}, MaxRecursion -> 0, PlotPoints -> {125, 125}, MeshFunctions -> {#4 &},
     Mesh -> {{0, wrap}}, MeshStyle -> Directive@{Thick, White}, 
    Boxed -> False, Axes -> False, PlotRange -> All, 
    Method -> {"BoundaryOffset" -> False}, PlotRange -> All, 
    PlotStyle -> {Red, Opacity[0.9]}, AxesLabel -> {"X", "Y", "Z"}];
  AppendTo[prof, pro];
  ];
 ]
Place = Show[prof, PlotRange -> All]

Show[RotorL, Place, PlotRange -> All] 


Solid Shape Region

reg1 = DiscretizeGraphics[Show[RotorL, PlotRange -> All]];
reg2 = DiscretizeGraphics[Place, PlotRange -> All];

merged = RegionUnion[reg1];
merged2 = RegionUnion[reg2];

findCaps[mesh_?MeshRegionQ] := 
 Block[{connectivity, neighborCounts, edgeCount, edgeIndices, edgeVertices, 
   holeGraphs, holePolygons, caps}, 
  connectivity = MeshConnectivityGraph[mesh, {1, 2}];
  neighborCounts = VertexDegree[connectivity];
  edgeCount = MeshCellCount[mesh, 1];
  edgeIndices = Select[Range[edgeCount], neighborCounts[[#]] === 1 &];
  edgeVertices = 
   MeshCells[mesh, 1, "Multicells" -> True][[1, 1]][[edgeIndices]];
  holeGraphs = 
   Select[Graph[edgeVertices] // ConnectedGraphComponents, 
    Length[VertexList[#]] > 10 &];
  holePolygons = Polygon[FindHamiltonianPath /@ holeGraphs];
  caps = MeshRegion[MeshCoordinates[mesh], holePolygons];
  caps /; MeshRegionQ[caps]]

caps = findCaps[merged];
closed = RegionUnion[merged, caps];
bmr = BoundaryMeshRegion[MeshCoordinates[closed], 
  MeshCells[closed, 2, "Multicells" -> True]]

caps2 = findCaps[merged2];
closed2 = RegionUnion[merged2, caps2];
bmr2 = BoundaryMeshRegion[MeshCoordinates[closed2], 
  MeshCells[closed2, 2, "Multicells" -> True]]

reg = RegionDifference[bmr, bmr2]

(* CRASH/ERROR HERE *)  
reg = RegionDifference[bmr, bmr2];

What I've Tried:

Reduced PlotPoints: Avoids crash but causes ::bsuncl.

RepairMesh/OpenCascadeLink:

Needs["NDSolve`FEM`"];  
reg1 = RepairMesh[reg1]; (* Helps but not fully *)  
Needs["OpenCascadeLink`"];  
diff = OpenCascadeShapeDifference[OpenCascadeShape[reg1], OpenCascadeShape[reg2]];

→ OpenCascade fails with non-watertight meshes.

Adjusted curve boundaries: Overlapped segments (e.g., tmin - 0.1°, tmax + 0.1°) → still has gaps.

Suspected Issues:

Gaps in transition curves between segments 4/5 and adjacent surfaces:

{min1[4], max1[4]} = {Pi + γ - 10°, 2π};  (* Might not connect to segment 3 *)

Cavity extrusion: The y parameter in cavities ({y, -20, 0}) may not seal properly.
Mesh complexity: 125×125 points = ~15k polygons per surface → too heavy.

Key Questions:

How can I ensure watertight meshes for parametric surfaces with complex transitions?
Is there a better way to close holes than findCaps? (I’m open to alternatives!)
How to simplify meshes without creating gaps?

System Info:

Mathematica 13.2 (or 12.3)
32 GB RAM / 8-core CPU

cvgmt · Accepted Answer · 2025-07-02 02:44:36Z

Edition

If we using the original codes( decreasing the points, change PlotPoints -> {100, 100} to PlotPoints -> {80, 80} and PlotPoints -> {125, 125} to PlotPoints -> {20, 20}) and by using OpenCascadeLink we can get the result although it is not so faster.

Needs["OpenCascadeLink`"];
Needs["NDSolve`FEM`"];
boundarymeshregion2polyhedron[bm_] := 
  Polyhedron[MeshCoordinates[bm], List @@@ MeshCells[bm, 2]];
shape1 = OpenCascadeShape[boundarymeshregion2polyhedron[bmr]];
shapes = OpenCascadeShape@ToBoundaryMesh[bmr2];
bm2 = OpenCascadeShapeSurfaceMeshToBoundaryMesh[
   OpenCascadeShapeDifference[shape1, shapes]];
BoundaryMeshRegion[bm2, ViewPoint -> Left]

Original

Not a fully answer, only a method similar with my previous answer or https://mathematica.stackexchange.com/a/308774/72111 which does not depend on the findCaps function when we build the RotorL.
Please note that we Reverse the {min1[2], max1[2]} to keep the orientation of the curves.
We consider the prof and Place later.

Clear["Global`*"];
wrap = 4.5 π;(*Max Wrap Angle*)
Ps = 70;(* 136/80 Start Lead Length in mm*)
Pe = 35;(* 24/80 End Lead Length in mm*)
P = Ps + ((Pe - Ps)/wrap)*t1;(*Pitch Equation*)
L = (1/(2 Pi))*\!\(
\*SubsuperscriptBox[\(∫\), \(0\), \(ϕr\)]\(P \
\[DifferentialD]t1\)\);(*Screw Rotor Length in mm*)
sega = Plot[P, {t1, 0, wrap}];

Out1 = 60;
R2 = (Out1 - (L Tan[α])); (*Outer Radius 53*)
Rp = 40;(*Pitch Radius*)
R1 = 2 Rp - R2;(*Inner Radius*)
α = 3*Degree;(*Conical Angle*)
γ = 160 Degree;(*160 Tooth Width Angle*)
β = Pi - γ;(*0.5 Transition Curve Angle*)
ρ = 0.0072;(*Iron cast density in g/mm^3*)

M = ( {
        {Cos[ϕr], Sin[ϕr], 0, 0},
        {-Sin[ϕr], Cos[ϕr], 0, 0},
        {0, 0, 1, +L},
        {0, 0, 0, 1}
       } );

{X[1], Y[1]} = {R1 Cos [t], R1 Sin[t]};(*Inner Radius*)

dd = -ArcCos[Rp/R2] // N;
{X[2], Y[2]} = {-(2 Rp) Cos[dd*t] + 
    R2 Cos [2 t*dd], -(2 Rp) Sin[t*dd] + 
        R2 Sin [2 t*dd]};(*Epicycloid Curve*)

{X[3], Y[3]} = {R2 Cos [t], R2 Sin[t]};(*Outer Rad*)

r1 = R2 - (R2 - 
               
       Rp)/(360 Degree - (Pi + γ - 
         10 Degree)) (t - (Pi + γ - 
                 10 Degree));(*Transition Curve A*)

{X[4], Y[4]} = {r1 Cos[t], r1 Sin[t]};

{r2} = {Rp - ((Rp - R1)/(β)) t};(*Transition Curve B*)

{X[5], Y[5]} = {r2 Cos[t], r2 Sin[t]};

{min1[1], max1[1]} =  {β, Pi};
{min1[2], max1[2]} = 
  Reverse@{1, 0}; (* <---------- Add Reverse here  ------------>  *)

{min1[3], max1[3]} = {Pi, Pi + γ - 10 Degree};
{min1[4], max1[4]} = {Pi + γ - 10 Degree, 2 Pi};
{min1[5], max1[5]} = {0, β};
(* tran=TransformationFunction[M];
tran@{X[i],Y[i],0}==Most[M.{X[i],Y[i],0,1}]//Simplify; *)

(* ------ the main code ------ *)

ptss = Table[
   Flatten[Table[
     M . {X[i], Y[i], 0, 1} // Most, {i, 1, 5}, {t, 
      Rest@Subdivide[min1[i], max1[i], 100]}], 1], {ϕr, 
    Rest@Subdivide[0, wrap, 120]}];
{m, n, p} = Dimensions[ptss];
indexes = 
  Most@Partition[
    ArrayReshape[Range[m*n], {m, n}], {2, 2}, {1, 1}, {1, 1}];
bm = BoundaryMeshRegion[
  Flatten[ptss, 
   1], {Map[Polygon[{#[[1, 1]], #[[1, 2]], #[[2, 2]]}] &, indexes, {2}],
    Map[Polygon[{#[[2, 1]], #[[2, 2]], #[[1, 1]]}] &, indexes, {2}], 
   Polygon[Range[n]], Polygon[Range[m*n - n + 1, m*n]]}]

Edit

Test the OpenCascadeLink. We try to make the RotorL as aPolyhedron. But since I do not know how to do with the prof and Place, we have to only test holes by using StadiumShape. We decreasing the points to make the result faster.

ptss = Table[
      Flatten[Table[
          M . {X[i], Y[i], 0, 1} // Most, {i, 1, 5}, {t, 
            Rest@Subdivide[min1[i], max1[i], 30]}], 1], {ϕr, 
        Rest@Subdivide[0, wrap, 80]}];
{m, n, p} = Dimensions[ptss];
indexes = 
    Most@Partition[
        ArrayReshape[Range[m*n], {m, n}], {2, 2}, {1, 1}, {1, 1}];
bm = BoundaryMeshRegion[
     Flatten[ptss, 
       1], {Map[Polygon[{#[[1, 1]], #[[1, 2]], #[[2, 2]]}] &, 
     indexes, {2}],
        Map[Polygon[{#[[2, 1]], #[[2, 2]], #[[1, 1]]}] &, 
     indexes, {2}], 
       Polygon[Range[n]], Polygon[Range[m*n - n + 1, m*n]]}];
Needs["OpenCascadeLink`"];
boundarymeshregion2polyhedron[bm_] := 
  Polyhedron[MeshCoordinates[bm], List @@@ MeshCells[bm, 2]];
polyhedron = boundarymeshregion2polyhedron[bm];
shape1 = OpenCascadeShape[polyhedron];
regs = Table[
   TransformedRegion[
    BoundaryDiscretizeRegion@
     RegionProduct[StadiumShape[{{0, -3}, {0, 3}}, 5], 
      Line[{{0.}, {50}}]], 
    RotationTransform[20 Degree + l*30 Degree, {0, 0, 1}]@*
     TranslationTransform[{0, 80, 10}]@*
     RotationTransform[π/2, {1, 0, 0}]], {l, 0, 4}];
shapes = 
 OpenCascadeShapeUnion[
  OpenCascadeShape /@ (boundarymeshregion2polyhedron /@ regs)];
bm2 = OpenCascadeShapeSurfaceMeshToBoundaryMesh[
   OpenCascadeShapeDifference[shape1, shapes]];
BoundaryMeshRegion[bm2, ViewPoint -> Left]

test another polygon.

Clear["Global`*"];
    poly = RandomPolygon[{"Simple", 20}];
    pts = PadRight[#, 3] & /@ poly[[1]];
    tran = RotationTransform[-4 θ, {0, 0, 1}]@*
   TranslationTransform[{-θ, -θ, -θ^2 - 
      20 θ + 1}]@*ScalingTransform[25 {1, 1, 1}];
    ptss = Table[tran /@ pts, {θ, 0, 2 π, .01}];
    {m, n, p} = Dimensions[ptss];
    indexes = 
      Most@Partition[
            ArrayReshape[Range[m*n], {m, n}], {2, 2}, {1, 1}, {1, 1}];
    bm = BoundaryMeshRegion[
      Flatten[ptss, 
       1], {Map[Polygon[{#[[1, 1]], #[[1, 2]], #[[2, 2]]}] &, 
    indexes, {2}],
        Map[Polygon[{#[[2, 1]], #[[2, 2]], #[[1, 1]]}] &, 
    indexes, {2}], 
       Polygon[Range[n]], Polygon[Range[m*n - n + 1, m*n]]}, 
  MeshCellHighlight -> {Polygon[Range[n]] -> Red, 
    Polygon[Range[m*n - n + 1, m*n]] -> Orange}]

do you think it is possible to make more than one hole in those 1st shape, let say we made it loops in the shape2? — kai29lol
– kai29lol, Commented Jun 25, 2025 at 3:07
This is great! I think the OpenCascade approach is the way to go. (+1). It might be easier to use a rectangle and sweep that along the path of the screw. One would then have the create a union of that shape with a cuboid to cut of the top and bottom of the swept shape to make it plane at the top and bottom. — user21
– user21, Commented Jun 25, 2025 at 3:35
@kai29lol No general way to determine the plot points and I think that Mathematica is not very effective for the region functions. — cvgmt
– cvgmt, Commented Jun 25, 2025 at 14:17

Stack Exchange Network

Kernel crash during RegionDifference of complex 3D meshes + BoundaryMeshRegion::bsuncl error

1 Answer 1

Linked

Hot Network Questions

Kernel crash during RegionDifference of complex 3D meshes + BoundaryMeshRegion::bsuncl error

1 Answer 1

Linked

Related

Hot Network Questions