How to refine a boundary mesh with MeshRefinementFunction?

MeshRefinementFunction works well to refine elements in a defined spatial domain (here x>0):

f2d = Function[{vertices, area},
Block[{x, y},
{x, y} = Mean[vertices];
If[x > 0 && area > 0.001, True, False]]
];
m = ToElementMesh[Disk[], MeshRefinementFunction -> f2d];
MeshRegion[m]


resulting in:

It also works for 3d meshes:

f3d3 = Function[{vertices, volume},
Block[{x, y, z},
{x, y, z} = Mean[vertices];
If[x > 0 && volume > 0.0001, True, False]]
];
m = ToElementMesh[Sphere[], MeshRefinementFunction -> f3d3];
m["Wireframe"["MeshElement" -> "MeshElements"]]


However, when I try that for a 2d surface mesh embedded in 3 dimensions, I do not get any refinement:

f3d2 = Function[{vertices, area},
Block[{x, y, z},
{x, y, z} = Mean[vertices];
If[x > 0 && area > 0.0001, True, False]]
];
bm = ToBoundaryMesh[Sphere[], MeshRefinementFunction -> f3d2];
MeshRegion[bm]


What am I doing wrong or does Mathematica not support spatial refinement for boundary meshes (I am using Version 12.1)? Is there any workaround?

I have also tried DiscretizeRegion[RegionBoundary[Sphere[]], MeshRefinementFunction -> f3d2] and BoundaryDiscretizeRegion[Sphere[], MeshRefinementFunction -> f3d2] without success. BoundaryDiscretizeRegion directly tells me that the option MeshRefinementFunction is unknown.

Edit: Unfortunately taking the boundary mesh of refined 3d mesh does not work either. The boundary points do not seem to be affected by the refinement of the bulk elements.

f3d3 = Function[{vertices, volume},
Block[{x, y, z},
{x, y, z} = Mean[vertices];
If[x > 0 && volume > 0.0001, True, False]]
];
m = ToElementMesh[Ball[], MeshRefinementFunction -> f3d3];
m["Wireframe"["MeshElement" -> "MeshElements"]]
bmesh = ToBoundaryMesh[m];
bmesh["Wireframe"["MeshElement" -> "PointElements"]]


• The refinement function never seems to be applied in the last case. For instance, Reap[ ToBoundaryMesh[Sphere[], MeshRefinementFunction -> (Sow[{##}] &)];] should show you the list of simplices passed to the refinement function, and their measure. That list is empty for this last case. – MarcoB Feb 2 at 2:46
• The "MeshRefinementFunction" is a true mesh refinement function. Unfortunately, there is no "BoundaryMeshRefinementFunction" – user21 Feb 2 at 7:21

In my extended comment to the question MeshRefinementFunction on 2D surfaces embedded in 3D, I showed that MaxCellMeasurecould be applied to 2D surfaces embedded in 3D, but that the MeshRefinementFunction seems to be ignored.

A potential workaround is to use the functionality in FEMAddOns to join two boundary meshes meshed at different resolutions. A sample workflow using BoundaryElementMeshJoin is shown below:

(*Uncommented the following function if FEMAddOns not \
installed*)
Needs["FEMAddOns"];
rl = ImplicitRegion[x^2 + y^2 + z^2 == 1 && x <= 0, {x, y, z}];
rr = ImplicitRegion[x^2 + y^2 + z^2 == 1 && x >= 0, {x, y, z}];
(bml = ToBoundaryMesh[rl]);
(bmr = ToBoundaryMesh[rr, MaxCellMeasure -> 0.00005]);
(bmj = BoundaryElementMeshJoin[bml, bmr])[
"Wireframe"["MeshElementStyle" -> FaceForm[Yellow]]]


Update 1: MeshRefinementFunction on a 2D mesh mapped to 3D

Here is another workaround that may allow you to use more complicated refinement strategies without having to break up your model into many regions. In this simple approach, we will map a rectangular region to a spherical region. We will create a structured Quad mesh and extract a boundary mesh. By setting "SteinerPoints"->False, the nodes should line up on the seam. Of course, a rectangular to spherical mapping will cause some distortions and you will have to make a judgment call if the final mesh is suited for purpose.

Helper functions

First, will define some helper functions to create a structured Quad mesh with refinement around the equator. Note that not all functions will be used in the workflow.

(*Define Some Helper Functions For Structured Meshes*)
pointsToMesh[data_] :=
MeshRegion[Transpose[{data}],
Line@Table[{i, i + 1}, {i, Length[data] - 1}]];
unitMeshGrowth[n_, r_] :=
Table[(r^(j/(-1 + n)) - 1.)/(r - 1.), {j, 0, n - 1}]
meshGrowth[x0_, xf_, n_, r_] := (xf - x0) unitMeshGrowth[n, r] + x0
firstElmHeight[x0_, xf_, n_, r_] :=
Abs@First@Differences@meshGrowth[x0, xf, n, r]
lastElmHeight[x0_, xf_, n_, r_] :=
Abs@Last@Differences@meshGrowth[x0, xf, n, r]
findGrowthRate[x0_, xf_, n_, fElm_] :=
Quiet@Abs@
FindRoot[firstElmHeight[x0, xf, n, r] - fElm, {r, 1.0001, 1/fElm},
Method -> "Brent"][[1, 2]]
meshGrowthByElm[x0_, xf_, n_, fElm_] :=
N@Sort@Chop@meshGrowth[x0, xf, n, findGrowthRate[x0, xf, n, fElm]]
meshGrowthByElm0[len_, n_, fElm_] := meshGrowthByElm[0, len, n, fElm]
flipSegment[l_] := (#1 - #2) & @@ {First[#], #} &@Reverse[l];
leftSegmentGrowth[len_, n_, fElm_] := meshGrowthByElm0[len, n, fElm]
rightSegmentGrowth[len_, n_, fElm_] :=
Module[{seg}, seg = leftSegmentGrowth[len, n, fElm];
flipSegment[seg]]
reflectRight[pts_] :=
With[{rt = ReflectionTransform[{1}, {Last@pts}]},
Union[pts, Flatten[rt /@ Partition[pts, 1]]]]
reflectLeft[pts_] :=
With[{rt = ReflectionTransform[{-1}, {First@pts}]},
Union[pts, Flatten[rt /@ Partition[pts, 1]]]]
extendMesh[mesh_, newmesh_] := Union[mesh, Max@mesh + newmesh]


Workflow to map a mesh refined 2D mesh to a sphere

The following workflow will accomplish:

• Create a structured Quad mesh with refinement around the equator.
• Extract the boundary mesh from the Quad mesh.
• Create a sinusoidal parametric region centered at the equator.
• Create a mesh refinement function based on parametric region.
• Create 2D mesh.
• Map 2D mesh to a sphere.

The code:

Print["Azimuthal mesh (horizontal)"]
rh = pointsToMesh@Subdivide[-π, π, 72]
Print["Inclination mesh (vertical)"]
rv = pointsToMesh[
reflectRight@rightSegmentGrowth[π/2, 25, π/2/40]]
Print["RegionProduct mesh"]
rp = RegionProduct[rh, rv]
(*Create boundary mesh from RegionProduct mesh*)
bmrect = ToBoundaryMesh@rp;
(*Create sinusoidal parametric region*)
pr = ParametricRegion[{ϕ, π/2 +
1/2 Sin[4 ϕ]}, {{ϕ, -π, π}}];
(* Create Mesh Refinement Function *)
mrf = With[
{rdf = RegionDistance[DiscretizeRegion@pr]},
Function[
{vertices, area},
Block[
{x, y}, {x, y} = Mean[vertices];
area > 0.00125 (1 + 100000 rdf[{x, y}]^8)
]
]
];
mrect = ToElementMesh[bmrect, "MeshElementType" -> TriangleElement,
MeshRefinementFunction -> mrf, "MeshOrder" -> 1,
"SteinerPoints" -> False, MaxCellMeasure -> 0.01];
Print["Sphere mesh with refinement mapped to 2D"]
mrect["Wireframe"]
(*Extract coordinates and incidents from mesh*)
crd = mrect["Coordinates"];
inc = ElementIncidents[mrect["MeshElements"]][[1]];
(*Map Spherical coordinates to 3D Cartesian*)
crd3d = crd /. {{ϕ_, θ_} -> { Cos[ϕ] Sin[θ],
Sin[θ] Sin[ϕ], Cos[θ]}};
mrkrs = ConstantArray[1, First@Dimensions@inc];
(*FEM Create BoundaryMesh*)
bm = ToBoundaryMesh["Coordinates" -> crd3d,
"BoundaryElements" -> {TriangleElement[inc, mrkrs]}];
Print["2D mesh mapped to a sphere"]
bm["Wireframe"["MeshElementStyle" -> FaceForm[Yellow]]]


With this approach, the mesh refinement about the equator looks okay, but you will get low quality and high aspect ratio triangles as you approach the poles.

Update 2: Capping the poles

@user21's recent answer to the question problem with DelaunayMesh 3D coordinates shows that you can create a DelaunayMesh simply by passing a coordinate list to ToElementMesh. We can mitigate some of the high aspect ratio elements at the pole by separately meshing the spherical end caps and joining the coordinates of all the meshes and using @user21's technique to create a DelaunayMesh. The approach is not completely seamless, but may be acceptable depending on your needs as shown in the following workflow:

angle = 30 °;
Print["Azimuthal mesh (horizontal)"]
rh = pointsToMesh@Subdivide[-π, π, 72]
Print["Inclination mesh (vertical)"]
rv = pointsToMesh[
reflectRight@(rightSegmentGrowth[π/2 - angle,
25, (π/2 - angle)/40] + angle)]
Print["RegionProduct mesh"]
rp = RegionProduct[rh, rv]
(*Create boundary mesh from RegionProduct mesh*)
bmrect = ToBoundaryMesh@rp;
(* Create Exponential Mesh Refinement Function *)
mrf = With[
{center = {0, First@Mean[MeshCoordinates[rv]]}},
Function[
{vertices, area},
Block[
{x, y}, {x, y} = Mean[vertices] - center;
area > 0.0001 (1 + 0.5 Exp[5* Norm[{x, y}]])
]
]
];
mrect = ToElementMesh[bmrect, "MeshElementType" -> TriangleElement,
MeshRefinementFunction -> mrf, "MeshOrder" -> 1,
"SteinerPoints" -> False, MaxCellMeasure -> 0.01];
Print["Uncapped sphere mesh with refinement mapped to 2D"]
mrect["Wireframe"]
(*Extract coordinates and incidents from mesh*)
crd = mrect["Coordinates"];
inc = ElementIncidents[mrect["MeshElements"]][[1]];
(*Map Spherical coordinates to 3D Cartesian*)
crd3d = crd /. {{ϕ_, θ_} -> { Cos[ϕ] Sin[θ],
Sin[θ] Sin[ϕ], Cos[θ]}};
mrkrs = ConstantArray[1, First@Dimensions@inc];
(*FEM Create BoundaryMesh*)
bm = ToBoundaryMesh["Coordinates" -> crd3d,
"BoundaryElements" -> {TriangleElement[inc, mrkrs]}];
Print["2D mesh mapped to an uncapped sphere"]
bm["Wireframe"["MeshElementStyle" -> FaceForm[Yellow]]]
(*Create spherical cap meshes*)
bmcaps = ToBoundaryMesh[
ImplicitRegion[
x^2 + y^2 + z^2 == 1 && (Cos[angle]^2 <= z^2), {x, y, z}],
MaxCellMeasure -> 0.001];
Print["Spherical caps"]
bmcaps["Wireframe"["MeshElementStyle" -> FaceForm[Green]]]
Print["Complete seamed surface mesh"]
crdsphere = bmcaps["Coordinates"];
ToBoundaryMesh[Join[crdsphere, bm["Coordinates"]]][
"Wireframe"["MeshElementStyle" -> FaceForm[Yellow]]]


• So, spatial refinement for boundary meshes is not implemented yet in Mathematica then? Too bad :( Your workaround is nice, if one would need a smoother transition in triangle sizes, then more subparts are necessary right? Calling distmesh from FEMAddOns migth be another option then I guess. – Oscillon Feb 2 at 5:01
• @Oscillon From @user21's comment, it appears that it is not implemented. You would need more subparts for a smoother transition. I tried "IncludePoints", but I could not get that to work on boundary meshes in 3D. I also looked at the details of DistMesh and it works for 2D only. – Tim Laska Feb 2 at 15:30
• Using the distmesh v1.1 from the authors webpage in matlab, I got a refined mesh with distmeshsurface(). If I find the time I will try to port it to mathematica. – Oscillon Feb 5 at 8:17
• Wow that is quite an inspiring workaround :) One might nitpick that there is still a circular edge around the pole due to the uncapping, but in lieu of a better answer I will accept yours. – Oscillon Feb 5 at 8:20

Here is an idea that works in 2D (not sure if it is going to work in 3D)

First we generate a mesh with a MeshRefinementFunction

f2d = Function[{vertices, area}, Block[{x, y}, {x, y} = Mean[vertices];
If[x > 0 && area > 0.001, True, False]]];
m = ToElementMesh[Disk[], MeshRefinementFunction -> f2d];


Then we extract the boundary mesh from that:

bmesh = ToBoundaryMesh[m];
bmesh["Wireframe"["MeshElement" -> "PointElements"]]


Next, we combine the symbolic region with the boundary mesh in a NumericalRegion and mesh that:

nr = ToNumericalRegion[Disk[]];
SetNumericalRegionElementMesh[nr, bmesh];
mesh = ToElementMesh[nr];
mesh["Wireframe"]
`

Give it a shot and see if this works in 3D too, it might not.

• I tried out your suggestion and posted the result in an edit to the question. No success unfortunately. – Oscillon Feb 2 at 22:36