# Visualising Finite Element Mesh in 3D with MeshRefinementFunction

How do you check that a MeshRefinementFunction is doing what you want in a 3D mesh where you only see the boundary mesh?

I started by making a block with a hole as follows

Needs["NDSolveFEM"]
Lx = 0.1; Ly = 0.05;
h = 0.04;
r = 0.01;
ir = ImplicitRegion[x^2 + y^2 >= r^2, {x, y, z}];
mesh = ToElementMesh[ir, {{-Lx, Lx}, {-Ly, Ly}, {0, h}}];
mesh["Wireframe"] I am interested in the region around the top of the hole so I decided to refine this using the following MeshRefinmentFunction

mrf = Compile[{{c, _Real, 2}, {a, _Real, 0}},
Block[{d, com, c0 = {0, 0, h}},
com = Total[c]/3;
d = Norm[com - c0];
If[d < 3 r && a > 10^-9, True, False]
]
];
ir = ImplicitRegion[x^2 + y^2 >= r^2, {x, y, z}];
mesh = ToElementMesh[ir, {{-Lx, Lx}, {-Ly, Ly}, {0, h}},
MeshRefinementFunction -> mrf];
mesh["Wireframe"] The problem here is that you only see the boundary mesh and the top of the hole looks like it does not have a refined mesh. I thus looked a bit closer with this

Show[mesh["Wireframe"], PlotRange -> {{-Lx/4, Lx/4}, {0, Ly/2}, All}] This does look like the top of the inside of the hole is not refined. However, this could be just some illusion so my attempt to check further was to get the coordinates and plot those in slices going from the bottom of the hole to the top

coords = mesh["Coordinates"];
inc = h/20;
Table[Graphics3D[{Point[coords]},
PlotRange -> {{-Lx/4, Lx/4}, {-Ly/2, Ly/2}, {z, z + inc}},
ViewPoint -> {0, 0, 1}], {z, 0, h - inc, inc}] I noted that I have lots of black points where the refinement starts but as I work up the density of black points reduces around the edges of the hole. This seems to confirm that the mesh is not refined near the top of the hole. However, this could be an illusion. So the first question is how best to check to see if the mesh has been properly refined?

I am also aware that my MeshRefinementFunction could be wrong or that for some reason the refinement process can't see the top of the hole. So the second question is: Is my MeshRefinementFunction wrong?

Thanks for any help.

Edit 1

A suggestion by Henrik Schumacher is to plot the MeshElements directly. Using this we can again plot in slices (fewer this time).

inc = h/20;
Table[Show[
mesh["Wireframe"["MeshElement" -> "MeshElements",
"MeshElementStyle" -> EdgeForm[{Black, Thin, Opacity[0.3]}]]],
PlotRange -> {{-Lx/4, Lx/4}, {-Ly/2, Ly/2}, {z, z + inc}},
ViewPoint -> {0, 0, 1}], {z, 0, h - inc, inc}] This again suggests that the refinement is not working near the top of the hole.

This is a short coming in TetGen and there is not much Wolfram Research can do about that. I think it is similar if not the same issue reported here.

How can you see that? If you sow and reap all tet coordinates that are marked for refinement you will find:

Needs["NDSolveFEM"]
Lx = 0.1; Ly = 0.05;
h = 0.04;
r = 0.01;
mrf = With[{r = r},
Compile[{{c, _Real, 2}, {a, _Real, 0}},
Block[{d, com, c0 = {0, 0, h}}, com = Total[c]/3;
d = Norm[com - c0];
If[d < 3 r && a > 10^-9, Sow[c]; True, False]]]];
ir = ImplicitRegion[x^2 + y^2 >= r^2, {x, y, z}];
re = Reap[
mesh = ToElementMesh[ir, {{-Lx, Lx}, {-Ly, Ly}, {0, h}},
MeshRefinementFunction -> mrf]];
Graphics3D[Tetrahedron[re[[2, 1]]], Boxed -> False] All those tets are marked for refinement. And if you look at

mesh["Wireframe"["MeshElement" -> "MeshElements",
PlotRange -> {{-Lx/4, Lx/4}, {0, Ly/2}, All}]] You see that those tets are refined. It's the tets at the boundary TetGen does not refine. But that can be remedied by refining the boundary.

We create the region in a different way. We discretize the primitives before we call RegionDifference. This allows us to refine the boundary of the cylinder.

br = RegionDifference[
BoundaryDiscretizeRegion[Cuboid[{-Lx, -Ly, 0}, {Lx, Ly, h}]],
BoundaryDiscretizeRegion[Cylinder[{{0, 0, -h}, {0, 0, 2 h}}, r],
"MaxCellMeasure" -> {"Area" -> 0.00000075}]];
mrf = With[{r = r},
Compile[{{c, _Real, 2}, {a, _Real, 0}},
Block[{d, com, c0 = {0, 0, h}}, com = Total[c]/3;
d = Norm[com - c0];
If[d < 3 r && a > 10^-9, True, False]]]];
mesh = ToElementMesh[br, MeshRefinementFunction -> mrf];
mesh["Wireframe"[PlotRange -> {{-Lx/4, Lx/4}, {0, Ly/2}, All}]] The mesh elements then look like this:

mesh["Wireframe"["MeshElement" -> "MeshElements",
PlotRange -> {{-Lx/4, Lx/4}, {0, Ly/2}, All}]] • Very impressive. You always seem to have the right answer at your fingertips. Thanks – Hugh May 23 '18 at 5:35
• @Hugh, I am glad I could help. – user21 May 23 '18 at 5:51

Edit:

All the functions from this answer are also included in a convenient package MeshTools.

I know that my answer does not directly address OP's question about visualization of 3D mesh, but I would like to show alternative approach to generate 3D meshes of requested shape. It relies on the fact that this cuboid geometry with a cylindrical hole can be thought of as an extrusion of 2D mesh. The workflow is a bit hacky and relies on external packages and FEMAddOns.

We use dimensions from OP's question and generate 2D mesh of triangles.

Needs["NDSolveFEM"]
Lx = 0.1;
Ly = 0.05;
r = 0.01;
ir = ImplicitRegion[x^2 + y^2 >= r^2, {x, y}];

(* Parameters in the functions are adjusted for 2D mesh. *)
mrf = Compile[{{c, _Real, 2}, {a, _Real, 0}},
Block[{d, com, c0 = {0, 0}},
com = Total[c]/3;
d = Norm[com - c0];
If[d < 2 r && a > 2*10^-5, True, False]]
];

mesh = ToElementMesh[ir, {{-Lx, Lx}, {-Ly, Ly}},
MeshRefinementFunction -> mrf,
"MeshQualityGoal" -> 1,
"MeshOrder" -> 1,
MaxCellMeasure -> 1
];
mesh["Wireframe"] Then we load the other packages (that we have already installed) and convert 2D mesh of triangles to quadrilaterals and apply smoothing to increase the quality of the mesh.

<< AceFEM Then we extrude the quadrilateral mesh for chosen thickness and specified number of elements through thickness. Definition of function is given bellow. I think this way we have more control over mesh density and often for FEM analyses hexahedral mesh is even more desired. Otherwise one could still write a function to split hexahedra to tetrahedra if that is necessary.

mesh3D = ExtrudeMesh[quadMesh, 0.04, 10]
mesh3D["Wireframe"["MeshElementStyle" -> FaceForm[LightBlue]]] (* Basics of this function are taken from AceFEM documentation. *)
ExtrudeMesh//ClearAll
ExtrudeMesh::badType="Only first order 2D quadrilateral mesh is supported.";
ExtrudeMesh[mesh_ElementMesh,thickness_/;thickness>0,layers_Integer?Positive]:=Module[{
n2D,nodes2D,nodes3D,elements2D,elements3D
},
If[
];

nodes2D=mesh["Coordinates"];
elements2D=Join@@ElementIncidents[mesh["MeshElements"]];
n2D=Length@nodes2D;
nodes3D=With[{dz=thickness/layers},
Flatten[#,1]&@Table[
Map[Join[#,{(l-1)dz}]&,nodes2D],
{l,layers+1}
]
];

elements3D=Flatten[#,1]&@Table[
Map[Join[n2D*(l-1)+#,n2D*l+#]&,elements2D],
{l,layers}
];

ToElementMesh[
"Coordinates"->nodes3D,
"MeshElements"->{HexahedronElement[elements3D]}
]
]

• Thanks for this. Good to see an alternative. How would you do the mesh refinement near the top of the hole? – Hugh May 23 '18 at 13:04
• A wild idea: one could try 1/4 symmetry of the 2D domain and use a StructuredMesh for that. Then project the coordinates to the other 3 quadrants and increase the incidents for the other 3 quadrants. Put all of that into ToElementMesh - the hope would be that it removes the duplicate coordinates. Then do the extrusion to 3D. Just a thought not sure if it will work. (I wish I had another +1 for FEMAddOns of course:-) – user21 May 23 '18 at 13:11
• @Hugh Do you mean refinement in thickness direction? I guess you can fiddle with definition of dz in ExtrudeMesh function and make it non constant. Or you could make 2 meshes, one with thinner elements through thickness, and stack one on top of another. – Pinti May 23 '18 at 18:05
• Yes in thickness direction. Those ideas might work well. Thanks – Hugh May 23 '18 at 18:06
• @user21 I am quite confident that would work. Although then extra refinement around the hole is not possible any more, because you are bound to regular spacing of structured mesh. – Pinti May 23 '18 at 18:13

I found the following here.

You can try out things like

mesh["Wireframe"[
"MeshElement" -> "MeshElements",
"MeshElementStyle" -> EdgeForm[{Black, Thin, Opacity[0.3]}]
]
] • Thanks this is very helpful. I have added this visualisation to my question. It suggests that the refinement function is not working proper. Would you agree? – Hugh May 22 '18 at 15:13
• Yes I do. Fortunately, user21 already resolved the problem. Please do not forget to accept his/her answer. – Henrik Schumacher May 23 '18 at 7:28