Its
"BoundaryFunction"property withFindRootis useful to "improve"RegionPlotoutput. I don't know if we can easily build this compiled function directly from a genericRegion.For some misterious (for me) reason, the output of
RegionPlot[numericalRegion["Predicates"], ...]appears more accurate thanRegionPlot[symbolicRegion, ...]. This turns out to be very important to get a valid MeshRegionMeshRegionwith the actual implementation because...
In the RegionPlot output improvement phase can happens thatboundary vertices are moved, and a boundary vertex moves andcan cross a mesh line. In this caseSo, at some times we need to reorder the vertices on a boundary, otherwise we can get an invalid / self-crossing boundary. This becomes a problem for fine meshes (in the following samples for MeshMesh option >= 50. A more accurate RegionPlot output is helpful to reduce the need of reordering.
I'd appreciate any suggestion to detect and fix this problem and to improve the Wholewhole implementation.
<< NDSolve`FEM`
Options[ConstrainedRegionPlotMeshGenerator] = {Mesh -> Automatic,
PlotPoints -> Automatic, MaxRecursion -> Automatic};
ConstrainedRegionPlotMeshGenerator[region_NumericalRegion,
opts : OptionsPattern[]] :=
ConstrainedRegionPlotMeshGeneratorCore[region, opts] // ToBoundaryMesh
ConstrainedRegionPlotMeshGenerator[region_ /; ConstantRegionQ[region],
opts : OptionsPattern[]] :=
ConstrainedRegionPlotMeshGeneratorCore[ToNumericalRegion@region, opts]
ConstrainedRegionPlotMeshGeneratorCore[region_NumericalRegion,
opts : OptionsPattern[]] :=
Module[{g, mptsx, mptsy, mpts, blns, mptq, grids, xg, yg, \[Delta],
pts, ptsx, ptsy, x, y},
(* draw region with RegionPlot and extract the GraphicsComplex *)
g = First@RegionPlot[
(*region["SymbolicRegion"],*)
region["Predicates"],
Evaluate[
Sequence @@
MapThread[
Prepend, {region["Bounds"], region["PredicateVariables"]}]],
Frame -> False, PlotStyle -> None, BoundaryStyle -> Red,
MeshStyle -> Blue,
Evaluate@FilterRules[{opts}, {Mesh, PlotPoints, MaxRecursion}]];
(* get mesh lines endpoints and boundary lines *)
{mptsx, mptsy} =
Union @@@
Cases[g, {Blue, lines___Line} :> {lines}, \[Infinity]][[All, All,
1, {1, -1}]];
If[Intersection[mptsx, mptsy] =!= {}, Return[$Failed]]; (*
Unsupported at present *)
mpts = DeleteDuplicates@Join[mptsx, mptsy];
blns =
Cases[g, {Directive[___, Red, ___], lines___Line} :>
lines, \[Infinity]];
(* build closed boundary lines selecting boundary lines vertices \
that are also on mesh lines *)
mptq = AssociationThread[mpts, Range@Length@mpts];
blns = Map[DeleteMissing[mptq /@ #] &, blns, {2}];
blns =
Map[If[Last@# == First@#, #, Append[#, First@#]] &, blns, {2}];
(* improve boundary vertices position *)
\[Delta] = region["BoundaryFunction"];
ptsx = {#[[1]],
y /. FindRoot[\[Delta][{#[[1]], y}], {y, #[[2]]}]} & /@
g[[1, mptsx]];
ptsy = {x /.
FindRoot[\[Delta][{x, #[[2]]}], {x, #[[1]]}], #[[2]]} & /@
g[[1, mptsy]];
pts = Join[ptsx, ptsy];
(* build result*)
BoundaryMeshRegion[pts, Sequence @@ blns]
];
For example, with a uniform mesh, and a MeshRegionMeshRegion output:
With an explicit, non-uniform grid, and a MeshRegionMeshRegion output:
With an explicit, non-uniform grid, and an ElementMeshElementMesh output:
