3D FEM Mesh problem - either bad mesh or inaccurate geometry

Question

I am having trouble obtaining a good 3D mesh. I have tried lots of stuff, either resulting in a bad mesh or inaccurate geometry.

My geometry consists of a cuboid with a quarter cylinder cut out on a side. I am running Mathematica 12.1 and can not update as I am adding to some existing code.

I am building a half cylinder using BooleanRegion which I then use to cut it out of a cuboid using RegionDifference.

What am I missing here? This is quite an easy geometry and every proper meshing tool would handle this with ease. How can I generate a proper mesh that is sufficiently representing my desired geometry?

Bad mesh:

I was able to achieve an adequate representation of the geometry by using BoundaryDiscretizeRegion on the tool-bodies that build the quarter cylinder. This results in lots of useless, small clusters of cells in the later mesh. The ToElementMesh seems to be forced to follow the meshregion from the BoundaryDiscreticeRegion.

Bad Geometry:

Just using the boolean operations on the primitives results in really poor representation of the geometry. I was able to achieve a somewhat decent, yet flawed geometry building it up in a complicated way.

Note the corner being cut here!

ClearAll["Global`*"]
Needs["NDSolve`FEM`"]
(**Parameters*)
sChannel = 0.8 ;(distance channels)
dChannel = 0.4 ;(diameter channel)
dPlatte = 0.5; (thickness plate*)
L = 3;
(Flow in x, height z, width y)
top = dPlatte + (dPlatte - dChannel/2)/2
bot = - dPlatte - dChannel/2 - (dPlatte - dChannel/2)
left = 0;
right = sChannel/2;
(******BAD MESH STUFF****)
(Geometry**)
maxCell = 0.001;
(block)
block = Cuboid[{0, left, bot},
   {L, right, top}];
Print["block"]
RegionPlot3D[block, Mesh -> All, MeshFunctions -> Automatic]
(channelMid)
channelMid = BooleanRegion[[Not] #2 [And] #1 &, 
   {BoundaryDiscretizeRegion[
     Cylinder [
      {{0 - 0.1, 0, 0},
       {L + 0.1, 0, 0}},
      dChannel/2],
     MaxCellMeasure -> maxCell],
    BoundaryDiscretizeRegion[
     Cuboid[
      {-0.2, -dChannel/2 - 0.1, 0},
      {L + 0.2, dChannel + 0.1, 0 + dChannel/2 + 0.1}],
     MaxCellMeasure -> maxCell]
    }
   ];
Print["channelMid"]
RegionPlot3D[channelMid, Mesh -> All, MeshFunctions -> Automatic]
(cut channel out of block)
domain = RegionDifference[block, channelMid];
Print["domain"]
RegionPlot3D[domain, Mesh -> All, MeshFunctions -> Automatic]
(**mesh**)
Print["Mesh0"]
mesh0 = ToElementMesh[domain]
mesh0["Wireframe"[
"MeshElementStyle" -> {FaceForm[Green], FaceForm[Red]}]]
(******BAD GEOMETRY STUFF*****)
(Geometry*)
blockBot = Cuboid[{0, left, -0.5},
   {L, right, 0}];
RegionPlot3D[blockBot, Mesh -> All, MeshFunctions -> Automatic]
cyli = Cylinder[{{0 - 0.1, 0, 0},
    {L + 0.1, 0, 0}},
   dChannel/2];
RegionPlot3D[cyli, Mesh -> All, MeshFunctions -> Automatic]
blockBot = BooleanRegion[[Not] #2 [And] #1 &, {blockBot, cyli}];
RegionPlot3D[blockBot, Mesh -> All, MeshFunctions -> Automatic]
blockTop1 = Cuboid[{0, left + dChannel/2, 0},
   {L, right, 0.5}];
domain1 = RegionUnion[blockBot, blockTop1];
RegionPlot3D[domain1, Mesh -> All, MeshFunctions -> Automatic]
blockTop2 = Cuboid[{0, left, 0},
   {L, +dChannel/2, 0.5}];
domain1 = RegionUnion[domain1, blockTop2];
RegionPlot3D[domain, Mesh -> All, MeshFunctions -> Automatic]
(**mesh**)
mesh1 = ToElementMesh[domain1]
mesh1["Wireframe"[
  "MeshElementStyle" -> {FaceForm[Green], FaceForm[Red]}]]

user21 · Accepted Answer · 2021-08-06T10:19:50.947

10

In 12.1 you will have to choose a different boundary mesh generator (this is no longer necessary in 12.3)

Set up the geometry:

block = Cuboid[{0, left, bot}, {L, right, top}];
channelMid = 
  RegionDifference[
   Cylinder[{{0 - 0.1, 0, 0}, {L + 0.1, 0, 0}}, dChannel/2], 
   Cuboid[{-2/10, -dChannel/2 - 1/10, 0}, {L + 2/10, dChannel + 0.1, 
     0 + dChannel/2 + 1/10}]];
region = RegionDifference[block, channelMid];

Mesh it:

mesh = ToElementMesh[region, 
   "BoundaryMeshGenerator" -> {"OpenCascade"}, 
   "MaxCellMeasure" -> Infinity];

Visualize:

groups = mesh["BoundaryElementMarkerUnion"];
temp = Most[Range[0, 1, 1/(Length[groups])]];
colors = ColorData["BrightBands"][#] & /@ temp;
mesh["Wireframe"["MeshElementStyle" -> FaceForm /@ colors, 
  ViewPoint -> {2, 2, 1}]]

Spend some quality time with the documentation, it saves much frustration.

edited Aug 06 '21 at 10:19

answered Aug 06 '21 at 05:27

user21

39,710
8
110
167

Thank you very much, this works perfectly! I was on the wrong track, as I thought the problem was not on the meshing side but occuring with the "channelMid" part. When I visualize the half cylinder with RegionPlot3D it still has it's sharp edges cut. – Tobias Aug 06 '21 at 07:19
@Tobias, RegionPlot uses a different meshing algorithm then the actual mesh generators. Use the mesh visualization to check the quality of the mesh. – user21 Aug 06 '21 at 11:00

3D FEM Mesh problem - either bad mesh or inaccurate geometry

1 Answers1