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

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?

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.

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["NDSolveFEM"]

(****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;

(****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]}]]

(*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]}]]

• Have you seen this? Aug 5, 2021 at 16:25

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,
"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.

• 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. Aug 6, 2021 at 7:19
• @Tobias, RegionPlot` uses a different meshing algorithm then the actual mesh generators. Use the mesh visualization to check the quality of the mesh. Aug 6, 2021 at 11:00