# How discretize a region placing vertices on a specific non-uniform grid

Given a generic region, for example:

Ω =
ImplicitRegion[
2 x^2 + 3 y^2 + 2 x y - 2 <= 0 ∧ x^2 + y^2 > .1, {x, y}];


and a non-uniform grid, for example:

Ωb = RegionBounds[Ω];
{xg, yg} =
N@Map[bound \[Function]
Range[bound[[1]],
bound[[2]], (bound[[2]] - bound[[1]])/20], Ωb];
{xg, yg} =
Map[g \[Function] {g[[1]],
Sequence @@ (g[[2 ;; -2]] RandomReal[{1 - .04, 1 + .04},
Length[g] - 2]), g[[-1]]}, {xg, yg}];


which, together, can be shown with:

Ωg =
RegionPlot[Ω, AspectRatio -> Automatic];

gg = Graphics[{LightGray,
Table[
Line[{{x, Ωb[[2, 1]]}, {x, Ωb[[2,
2]]}}], {x, xg}],
Table[
Line[{{Ωb[[1, 1]], y}, {Ωb[[1, 2]],
y}}], {y, yg}]
}];

Show[gg, Ωg]


I'm searching a way to discretize the region to a MeshRegion or ElementMesh in such a way that all vertices (MeshCoordinates or incidents) are placed on at least one gridline.

At present I'm working on a routine that discretize the region with BoundaryDiscretizeRegion with a reasonable MaxCellMeasure. Then I'm trying to split all Line mesh cells crossing any gridline into two or more Line such that at least one end is on a gridline. Then I plan to delete vertices not on a grid line and properly reconnect the sourronding vertices.

This is the only way I could think to, but it appear difficult, because there are many branches tho examine.

A more elegant and simple approach would be helpful. For example, there is a way to use ToElementMesh and ToBoundaryMesh to accomplish this task? In the end, the routine should work with thousand of vertices and hundred of gridlines in a reasonable time.

• (1) Do you want the internal mesh coordinates to lie on the gridlines or just the boundary coordinates? (2) How generic? The example has nice algebraic boundaries without singular points. – Michael E2 Dec 24 '14 at 2:21
• @MichaelE2 As of my self-answer, I finally found a way, not simple, not elegant, not fast, but it works, more or less. (1) I'm only interested to mesh the boundary in this way. (2) The region can have corners like a square, but it's not "pathological"; can be a derived region or a ParametricRegion; generally it's something for which RegionQ gives True; I cannot count on having a representation as inequalties. – unlikely Dec 24 '14 at 16:30

regplt =   RegionPlot[\[CapitalOmega], AspectRatio -> Automatic];

ContourPlot[{2 x^2 + 3 y^2 + 2 x y - 2, x^2 + y^2 - .1}, {x, -1.25, 1.25}, {y, -1.25, 1.25},
Contours -> {{0}}, BaseStyle -> Thick,
GridLines -> {xg, yg}, Method -> {"GridLinesInFront" -> True},
MeshFunctions -> {#1 &, #2 &}, Mesh -> {xg, yg},
MeshStyle -> {Directive[{Red, PointSize[.01]}],
Directive[{Green, PointSize[.01]}]},
ImageSize -> 400, Prolog -> regplt[[1]]]


Update:

... if I have a Region and not its implicit description

you can recover the implicit description of the region using:

RegionRegionProperty[\[CapitalOmega], {x, y},"ImplicitDescription"]


or

RegionRegionProperty[ \[CapitalOmega], {x, y},"FastDescription"][[1, 2]]


to get

-2+2 x^2+2 x y+3 y^2 <= 0 && x^2+y^2 > 0.1

• +1 for some interesting tricks I wasn't aware of. – unlikely Dec 24 '14 at 11:54
• However: 1) This can be turned into a BoundaryMeshRegion extracting from the Graphics the boundaries? maybe... 2) mesh points are enough accurate? I suppose they are only accurate for Graphics purposes 3) most important... I cannot apply this if I have a Region and not its implicit description... – unlikely Dec 24 '14 at 15:20
• another interesting undocumeted trick... does it works for any Region that can also be discretized? – unlikely Dec 26 '14 at 10:59
• @unlikely, please see the answers by Silvia, m_ goldberg,belisarius, and Simon to this question and the q/a linked there. – kglr Dec 26 '14 at 11:16
• I think the "ImplicitDescription" has floated up to the top level in RegionMember[reg, {x, y}]. (Or Simplify[RegionMember[reg, {x, y}], (x | y) \[Element] Reals] if you want exactly the same expression.) – Michael E2 Dec 26 '14 at 18:11

This is the same idea as unlikely's, which has same following drawback. The discretization of the region makes the boundary piecewise linear and thereby reduces the accuracy. The points found on the grid are points on the linear approximations of the boundary. This causes the boundary to contract further (on the concave side), thereby further reducing the accuracy. Ideally, a root-finding procedure should be applied to adjust the points found in this method toward the boundary of the region. Given that the regions are generic, it seemed difficult to hit on a good generic algorithm. When the equation of the boundary is known, FindRoot is the obvious choice. In other cases, it or something else might be efficient. In any case, I left that adjustment to be done.

The idea is simple. For each boundary line/curve/path, find the grid lines crossed by each segment and solve for the intersections. Do this in order, so that the boundary element of the region can be constructed. The intersections are simple linear equations that can be conveniently solved with Rescale.

snaptogrid[reg_, xg_, yg_] :=
Module[{ptslists, coords, tpos,
xsgn, xgCrossings, xt,   ysgn, ygCrossings, yt},

ptslists = Function[{component},
coords = Part[
MeshCoordinates[component],
Append[#[[All, 1]], #[[-1, 2]]] &[MeshCells[component, 1] /. Line -> Identity]];
DeleteDuplicates @
Flatten[
Table[
tpos = t + {0, 1};
xsgn = Evaluate@UnitStep[-Apply[Times, # - coords[[tpos, 1]]]] &;
xgCrossings = Extract[xg, SparseArray[xsgn@xg]["NonzeroPositions"]];
xt = Rescale[xgCrossings, coords[[tpos, 1]], tpos];     (* x grid crossing times *)
ysgn = Evaluate@UnitStep[-Apply[Times, # - coords[[tpos, 2]]]] &;
ygCrossings = Extract[yg, SparseArray[ysgn@yg]["NonzeroPositions"]];
yt = Rescale[ygCrossings, coords[[tpos, 2]], tpos];     (* y grid crossing times *)

Function[{t}, Rescale[t, tpos, #] & /@ Transpose@coords[[tpos]]] /@
Sort@Flatten[{xt, yt}],                (* convert crossing times to coordinates *)

{t, Length[coords] - 1}],
1]
] /@ ConnectedMeshComponents[RegionBoundary[reg]];
BoundaryMeshRegion[
Flatten[ptslists, 1],
Sequence @@
(Line /@ Partition[#, 2, 1, 1] & /@
(InternalPartitionRagged[Range@Total[#], #] &[Length /@ ptslists]))
]
]


OP's example:

bΩ = snaptogrid[BoundaryDiscretizeRegion[Ω], xg, yg]; // AbsoluteTiming
(*  {0.237352, Null} *)

Show[
bΩ, gg,
Graphics[{Red, Point[MeshCoordinates[bΩ]]}]
]


Note: For alternatives to InternalParititionRagged, see Partitioning with varying partition size

• Interesting, "compact" routine, I had to study many tricks you used to have an idea of how it works. However, on my tests with MaxCellMeasure->0.2/200 and Length[xg] == Length[yg] == 200 I get an AbsoluteTiming of about 5 sec, while a sligtly revised version of mine (avoiding NSolve) requires about 0.3 sec.. Maybe it's because you check every boundary element against possible intersections with every gridline? – unlikely Dec 26 '14 at 17:35
• @unlikely That's a good guess. There would be a lot unnecessary testing. – Michael E2 Dec 26 '14 at 17:51

## MOST RECENT UPDATE

I sligtly revised the code, mainly to avoid NSolve to compute lines intersections and to directly get intersection in the proper order, if possible.

This first helper function, with a binary search algorithm, give the gridline index to wich all coordinates vals given belong, or a half-odd if not on gridline; sameTest allow to configure the accuracy of the check.

BinPositions[vals_, brakes_, sameTest_] :=
Map[val \[Function]
Catch@Module[{lo = 1, mid, hi = Length[brakes], el, res},
While[lo <= hi, Which[
TrueQ@sameTest[val, el = brakes[[mid = Floor[(lo + hi)/2]]]],
Throw[mid],
el > val, hi = mid - 1,
True, lo = mid + 1]];
lo - 1/2],
vals]


This second helper function, based on some output of the previous functions, gives the indices of the crossed gridlines. Maybe we can do better.

NearestIntegersBetween = {m, n} \[Function]
With[{\[Delta] = n - m, p = IntegerQ[m]}, Which[
\[Delta] > 1 || \[Delta] == 1 && p, {Ceiling[m], Floor[n]},
\[Delta] > 0, {Floor[m] + 1},
\[Delta] == 0 && p, {m, n},
-\[Delta] > 1 || -\[Delta] == 1 && p, {Floor[m], Ceiling[n]},
-\[Delta] > 0, {Floor[n] + 1},
True, {}
]];


This third helper function compute the intersections of a line segment with the crossed gridlines, in the order of the segment.

SegmentGridIntersections =
{x1, y1, x2, y2, xl, yl} \[Function]
Module[{m11 = y2 - y1, m12 = x1 - x2, v1 = x1 y2 - x2 y1},
Which[
Length@xl == 0,
{(v1 - m12*yl)/m11, yl}\[Transpose],

Length@yl == 0,
{xl, (v1 - m11*xl)/m12}\[Transpose],

True, Join[
{(v1 - m12*yl)/m11,
yl}\[Transpose], {xl, (v1 - m11*xl)/m12}\[Transpose]
] // Sort // If[x1 <= x2, #, Reverse@#] &
]
];


This last helper function process a single contour of a MeshRegion:

AdjustPolygonToGrid[vertices_, grids_, sameTest_] :=
Module[{positions, lines},
positions = MapThread[BinPositions[#1, #2, SameTest -> sameTest] &, {vertices\[Transpose], grids}]\[Transpose];
lines = Apply[NearestIntegersBetween,Transpose[Partition[positions, 2, 1], {1, 3, 2}], {2}];
lines = MapThread[Extract, {grids, Map[List, lines\[Transpose], {2}]}]\[Transpose];
MapThread[SegmentGridIntersections[Sequence @@ Flatten@#1, Sequence @@ #2] &, {Partition[vertices, 2, 1], lines}] // Flatten[#, 1] & // Append[#, #[[1]]] &
];


This main function finally process a whole MeshRegion.

AdjustMeshToGrid[meshRegion_, grids_, sameTest_] :=
Module[{polygons, vertices, map},
polygons = meshRegion["BoundaryPolygons"][[All, 1]];
polygons = AdjustPolygonToGrid[#, grids, sameTest] & /@ polygons;

vertices = DeleteDuplicates[Join @@ polygons];

(* Restituisce la mesh adattata *)
BoundaryMeshRegion[vertices,
Sequence @@ Line /@ Map[map, polygons, {2}]]
]


With thid code there are the results compared to the routine proposed by @Michael E2 on a uniform grid.

timeAvg =
Function[func,
Do[If[# > 0.3, Return[#/5^i]] & @@ Timing@Do[func, {5^i}], {i, 0,
15}], HoldFirst];

\[CapitalOmega] =
ImplicitRegion[
2 x^2 + 3 y^2 + 2 x y - 2 <= 0 \[And] x^2 + y^2 > .1, {x, y}];
\[CapitalOmega]b = RegionBounds[\[CapitalOmega]];

data = Table[Module[{grids, mesh, n},
grids =
N@Map[range \[Function]
Range @@
Append[range, -Subtract @@ range/20/2^k], \[CapitalOmega]b];
(*grids=Union/@Map[g\[Function]{g[[1]],
Sequence@@(g[[2;;-2]]RandomReal[{1-.02,1+.02},Length[g]-2]),
g[[-1]]},grids];*)
n = Length@First@grids;
mesh =
BoundaryDiscretizeRegion[\[CapitalOmega],
MaxCellMeasure -> Mean@Flatten[Differences /@ grids]];
<|
"Gridlines Count" -> n,
Abs[#1 - #2] <= 10.^-10 &] // timeAvg),
"snaptogrid" -> (snaptogrid[mesh, Sequence @@ grids] // timeAvg)
|>
], {k, 0, 10}] // Dataset


ListPlot[
Values@*Normal /@ {data[
data[All, {"Gridlines Count", "snaptogrid"}]},
Joined -> True, PlotLegends -> {"AdjustMeshToGrid", "snaptogrid"},
Frame -> True, Mesh -> Full, GridLines -> Automatic,
FrameTicks -> Automatic, PlotRange -> All]


I finally found a relatively short way to do. Not perfect, many special cases should be handled. Michael E2 routine permorms better on small uniform grids. I don't know why but on small random grids the rating is reversed.

I can do:

mesh = BoundaryDiscretizeRegion[\[CapitalOmega]];
polygons = mesh["BoundaryPolygons"][[All, 1]];
polygons =
AdjustPolygonToGrid[#, {xg, yg}, Abs[#1 - #2] <= 10^-6 &] & /@
polygons;
vertices = DeleteDuplicates[Join @@ polygons];
Sequence @@ Line /@ Map[verticesMap, polygons, {2}]];

Frame -> True, GridLines -> {xg, yg},
GridLinesStyle -> Darker[LightGray]]


with the main helper's function:

AdjustPolygonToGrid[vertices_, grids_, sameTest_] :=
Module[{positions, v1, p1},
positions =
BinPositions[#1, #2,
SameTest -> sameTest] &, {vertices\[Transpose],
grids}]\[Transpose];
{v1, p1} = {vertices[[1]], positions[[1]]};
Apply[{v2, p2} \[Function] Module[{lines, vlnew},
If[AnyTrue[p2, IntegerQ], Sow[v2],
If[lines =!= {{}, {}},
lines = MapThread[#1[[#2]] &, {grids, lines}];

lines =
RegionUnion @@
Flatten@{InfiniteLine[{#, 0}, {0, 1}] & /@ lines[[1]],
InfiniteLine[{0, #}, {1, 0}] & /@ lines[[2]]};

vlnew =
Block[{x, y}, {x, y} /.
NSolve[{x, y} \[Element]
Line[{v1, v2}] \[And] {x, y} \[Element] lines, {x,
y}, Reals]];

vlnew =
If[OrderedQ[{v1, v2}], Sort[vlnew], Reverse@Sort[vlnew]];
Sow /@ vlnew;
];
];
{v1, p1} = {v2, p2}
], RotateLeft@Transpose@{vertices, positions}, {1}] //
Reap // Last // Last // Append[#, #[[1]]] &
]


and two simple helper's function:

BinPositions[vals_, brakes_, sameTest_] :=
Map[val \[Function]
Catch@Module[{lo = 1, mid, hi = Length[brakes], el, res},
While[lo <= hi, Which[
TrueQ@sameTest[val, el = brakes[[mid = Floor[(lo + hi)/2]]]],
Throw[mid],
el > val, hi = mid - 1,
True, lo = mid + 1]];
lo - 1/2],
vals]

NearestIntegersBetween = {m, n} \[Function]
Which[#1 < #2, {##}, #1 == #2, {#1},
True, {}] & @@ ({Floor[#1] + 1, Ceiling[#2 - 1]} & @@
Sort[{m, n}]);


For example original and ajusted small polygon:

polygon = mesh["BoundaryPolygons"][[1, 1]];
AdjustPolygonToGrid[polygon, {xg, yg}, Abs[#1 - #2] < 10^-6 &];
Graphics[{
Opacity[.5], LightGray, Polygon@polygon, Opacity[1],
Orange, AbsolutePointSize[Large], Point@polygon,
Blue, Line@%,
Red, AbsolutePointSize[Medium], Point@%,
Green, Point@polygon[[1]],
Yellow, Point@polygon[[2]]},
PlotRange -> RegionBounds@Line@polygon,
PlotRangePadding -> Scaled[.05], Frame -> True,
GridLines -> {xg, yg}, GridLinesStyle -> Darker[LightGray]]


The overall procedure is a bit slow, also with this simple mesh. Some special cases are not handled (for example if two or more vertices of the original mesh are on the same cell edge).

• I think you may be interested to know that ToBoundaryMesh and ToElementMesh have an option "BoundaryMeshGenerator" that you can define yourself. So you could use your boundary mesh generator for FEM, for example. If you try that, I'd be very interested to hear about the results. This approach would also allow you to get a quadratic (boundary)mesh. If you try this let me know. – user21 Jan 5 '15 at 14:32
• @user21 This would be probably the best approach, but I suppose it's not easy. Can you give me just an idea of how this is done? Maybe using the idea of kguler we can implement a "ContourPlot" generator, more or less like the builtin "RegionPlot" generator? – unlikely Jan 13 '15 at 15:44
• This should not be too hard. In the ref page of ToBoundaryMesh in the options section there is a BoundaryMeshGenerator section. Have a look there. Your function will get a NumericalRegion (also documented) from that you can extract the "Predicates" compute your boundary mesh and package that with ToBoundaryMesh["Coordinates"->coords, "BoundaryElements" ->{...}] let me know if this works for you. – user21 Jan 14 '15 at 7:59
• @user21 You're right, it's not so hard. Indeed, I previously skipped some part of the documentation of the "BoundaryMeshGenerator" option. I used the simplest trick of kguler but with RegionPlot instead of ContourPlot. If someone interested, I can post another self-answer. I'm concerned about the accuracy of this trick. About builtin "RegionPlot" method for "BoundaryMeshGenerator" it's stated that "improve" the RegionPlot output. Can you give me an idea of what this means? – unlikely Jan 15 '15 at 9:26
• I'd very much appreciate if you could post that! The improvement is two fold: 1) the boundary points are moved back onto the symbolic boundary if possible and the mesh can be second order. – user21 Jan 15 '15 at 10:00

Building on the trick of kguler and on the @user21 recommendations I ends up with the following approach.

For many reason my actual interest is on first-order meshes and on a MeshRegion-type output. However, I used the some services of FEM toolbox and particularly of NumericalRegion.

• Its "BoundaryFunction" property with FindRoot is useful to "improve" RegionPlot output. I don't know if we can easily build this compiled function directly from a generic Region.

• For some misterious (for me) reason, the output of RegionPlot[numericalRegion["Predicates"], ...] appears more accurate than RegionPlot[symbolicRegion, ...]. This turns out to be very important to get a valid MeshRegion with the actual implementation because...

In the RegionPlot output improvement phase boundary vertices are moved, and a boundary vertex can cross a mesh line. So, 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 Mesh 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 whole implementation.

<< NDSolveFEM

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 @@
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 *)
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 MeshRegion output:

ConstrainedRegionPlotMeshGenerator[\[CapitalOmega], Mesh -> {5, 10}]


With an explicit, non-uniform grid, and a MeshRegion output:

Show[
ConstrainedRegionPlotMeshGenerator[\[CapitalOmega], Mesh -> grids],
GridLines -> grids, GridLinesStyle -> LightGray,
Method -> {"GridLinesInFront" -> True}
]


With an explicit, non-uniform grid, and an ElementMesh output:

mesh = ToBoundaryMesh[\[CapitalOmega],
"BoundaryMeshGenerator" -> {ConstrainedRegionPlotMeshGenerator,
Mesh -> grids}]

Show[
mesh["Wireframe"],
mesh["Wireframe"["MeshElement" -> "PointElements"]],
GridLines -> grids, GridLinesStyle -> LightGray
]


• Very cool! Thanks for sharing. That's very helpful for me to understand what and if the functionality in the FEM context get's used. If you have suggestions for improving functionality / documentation in the FEM context let me know – user21 Jan 16 '15 at 8:51