# Difficulty in specifying mesh refinement

I am trying to give a region within a 3D volume a finer mesh than the rest of the volume. My problem is more complicated but here is a minimum working example.

I define a 3D cuboid and try and have a fine mesh in a region between an inner and outer hemisphere.

Needs["NDSolveFEM"];
Len = 150/1000;  (*length of plate  *)
ht = 15/1000; (* height *)
wd = 80/1000; (* width *)
reg = BoundaryDiscretizeRegion[
Cuboid[{-Len/2, -wd/2, 0}, {Len/2, wd/2, ht}]]


   rcp = 0.02;
mrf = Compile[{{c, _Real, 2}, {a, _Real, 0}},
Block[{d, com},
com = Total[c]/3;
d = Norm[com];
If[0.8 rcp < d < 1.2 rcp && a > 10^-10, True, False]
]
];
mesh = ToElementMesh[reg, MeshRefinementFunction -> mrf];
Show[mesh["Wireframe"],
PlotRange -> {{-Len/2, Len/2}, {-wd/2, wd/2}, {0, ht/4}}]


This clearly shows an annulus on the boundary but it is not the correct size. If I extract the boundary points and plot them together with the circles defining the region of refinement we see that the annulus is too small.

cc = mesh["Coordinates"];
c1 = Select[cc, #[[3]] == 0 &];
ListPlot[c1[[All, {1, 2}]],
PlotRange -> {{-0.03, 0.03}, {-0.03, 0.03}},
AspectRatio -> Automatic,
Epilog -> {Circle[{0, 0}, 0.8 rcp], Circle[{0, 0}, 1.2 rcp],
Red, PointSize[0.02], Point[{0, 0}]}]


Have I done something silly? What is wrong? I know I am taking my reference from the centroid of each element and this will make a slight difference but if I increase the "MaxCellMeasure" it does not change the position.

Edit

Due to Henrik Schumacher and user21, who have very sharp eyes, they show an elementary mistake. I averaged the 4 points defining the element by dividing by 3. I probably picked this up by copying the MeshRefinementFunction from a 2D problem. Here is the corrected version. I use Mean as a variant on the two other methods given in their answer.

rcp = 0.02;
mrf = Compile[{{c, _Real, 2}, {a, _Real, 0}},
Block[{d, com},
com = Mean[c];
d = Norm[com];
If[0.8 rcp < d < 1.2 rcp && a > 10^-10, True, False]
]
];
mesh = ToElementMesh[reg, MeshRefinementFunction -> mrf];
Show[mesh["Wireframe"],
PlotRange -> {{-Len/2, Len/2}, {-wd/2, wd/2}, {0, ht/4}}]

cc = mesh["Coordinates"];
c1 = Select[cc, #[[3]] == 0 &];
ListPlot[c1[[All, {1, 2}]],
PlotRange -> {{-0.03, 0.03}, {-0.03, 0.03}},
AspectRatio -> Automatic,
Epilog -> {Circle[{0, 0}, 0.8 rcp], Circle[{0, 0}, 1.2 rcp],
Red, PointSize[0.02], Point[{0, 0}]}]


All working properly now.

We want to generate a tet mesh, so

com = Total[c]/4;


or better

com = Total[c]/Length[c];


would be more appropriate.

• We are dealing with tets here. It's a coincidence that a 2D quad has 4 vertices and a 3D tet has 4 too. – user21 Jun 6 '19 at 5:32
• But it's a good idea to use Length[c] in stead of a fixed number and I updated the ToElementMesh documentation to use that. – user21 Jun 6 '19 at 5:38
• @user21 Aaah, that makes sense. – Henrik Schumacher Jun 6 '19 at 5:44
• I made an edit. I hope that if fine? reg is not relevant here. It's only the boundary that is given to TetGen and then, during the mesh generation of the tets the relevant tets are refined. – user21 Jun 6 '19 at 5:48
• @user21 It's all fine. Thank you! My quick fix was a bit too quick... =) – Henrik Schumacher Jun 6 '19 at 5:50