Mesh Based Regions Workshop

Rob Knapp

Wolfram Research Numerical Computation Development Team

SetDirectory[NotebookDirectory[]]

Hand Built Boundary Mesh

To get started, let construct a boundary mesh that approximates a unit disk.

The set of all points inside a closed curve (in 2D) or surface (in 3D) is known as a boundary representation for the point set.  
A BoundaryMeshRegion is a boundary representation for a region where the closed curve or surface is given as a mesh.
In 2D, the curve mesh is a collection of Line segments.

circlePoints = Compile[{{n, _Integer}},Module[{θ=0.,Δθ=2. Pi/n},Table[θ=Δθ i;{Cos[θ],Sin[θ]},{i,n}]]];

bdisk[n_] := BoundaryMeshRegion[circlePoints[n],Line[Append[Range[n],1]]]

bdisk[10]

Experiment with this and see if you can come up with an approximate model for how fast the area converges to π as a function of n.
Hint: RegionMeasure might be useful

Convergence

data = Table[{2^k,π-RegionMeasure[bdisk[2^k]]},{k,2,10}]

ListPlot[Log[data]]

Clear[n];
model[n_] =  Exp[Fit[Log[data],{1,Log[n]},Log[n]]]

Triangulating

An important operation on a boundary mesh is to subdivide it into triangles.  This can be done with TriangulateMesh

tdisk = TriangulateMesh[bdisk[5]]

Experiment with triangulating boundary disks with different numbers of points.  You might want to modify the size of the triangles using the TriangulateMesh option MaxCellMeasure.

Manipulate[TriangulateMesh[bdisk[Ceiling[10^ln]],MaxCellMeasure->10.^lm],{{ln,1,"(Subscript[log, 10]) n"},1,4},{{lm,-2, "(Subscript[log, 10])m"},-4,0}]

Mesh Representations

To better understand how to work with meshes, it helps to understand the representation used.

A MeshRegion consists of a collection of cells.  A cell can be any of the following

The indexes (Subscript[i, k])refer to positions in the coordinates.

Note that there is a heirarchy of cells -- e.g. a Tetrahedron, has faces that are Polygon, and each of those Polygon faces has edges that are Line and each of those Line edges has vertices that are Point.

A BoundaryMeshRegion in dimension d consists only of cells of dimension d-1 that together make a closed surface (or set of closed surfaces).

Exercise

Construct a mesh region that is a unit square with 2 triangle cells

coordinates = {{0,0},{1,0},{0,1},{1,1}};

Solution

mr = MeshRegion[coordinates,Polygon[{{1,2,4},{1,3,4}}]]

Orientation

When possible the input is corrected so that the orientation of cells is consistent.
In 2D orientation of polygons is chosen to be counterclockwise.

MeshCells[mr,2]

MeshRegions are Regions

A MeshRegion or BoundaryMeshRegion is a region

RegionQ[mr]

so functions that work on regions generally work with meshes too.  These include:

RegionDimension

RegionEmbeddingDimension

RegionMeasure

RegionCentroid

RegionMember

RegionNearest

RegionDistance

SignedRegionDistance

Manipulate[ContourPlot[Evaluate[SignedRegionDistance[bdisk[n]][{x,y}]],{x,-1,1},{y,-1,1}],{{n,10},3,100,1}]

For all but the simplest symbolic regions, these commands can be much faster with meshes than the corresponding symbolic region

test = RandomReal[{-1,1},{10^6,2}];
Map[First,{AbsoluteTiming[SignedRegionDistance[bdisk[10],test];],AbsoluteTiming[smf = SignedRegionDistance[Disk[],test];]}]

Exercise

As an exercise, write an acceptance-rejection method to find random points that lie in the space between the discretized disk with n sides and the actual disk.
Hint: RegionMember

A simple solution

This is a text cell.

randomBetween[n_, npts_]:= Module[{df = RegionMember[bdisk[n]],sf =RegionMember[Disk[]],data,i=0,p},
data = Reap[
While[i < npts,
p=RandomReal[{-1,1},2];
If[UnsameQ[df[p],sf[p]],i++;Sow[p]]
]];
data[[-1,1]]
]

ListPlot[randomBetween[4,1000],AspectRatio->1]

Mesh Cell Indexing

In a MeshRegion and cells are indexed by {d, i} where d is the dimension of the cell and i is the index number for the cell in that dimension.

Sets of cells can be referred to by {d,ispec} where ispec is a specification for cell indexes, such as “Boundary” or a pattern.

For a simple example, consider

mesh = TriangulateMesh[bdisk[5],MaxCellMeasure->.25]

Table[HighlightMesh[mesh,Labeled[{d,All},Placed["Index","Centroid"]]],{d,0,2}]

bcells = MeshCells[mesh,{1,"Boundary"}]

MeshCellIndex[mesh,{1,"Boundary"}]

MeshCellIndex[mesh,bcells]

Experiment with different sets of cells and see what they are with the labels above

Properties

Each mesh cell has a set of properties, some of that are intrinsic and some that may be assigned.
They may be accessed via PropertyValue

Show[mesh,Graphics[{Blue,Point[PropertyValue[{mesh,2},MeshCellCentroid]]}]]

PropertyValue[{mesh, {1,"Boundary"}},MeshCellMeasure]

Total[%]

Properties that are not intrinsic (such as display properties) can also be set.

PropertyValue[{mesh,{1,"Boundary"}},MeshCellStyle] = Red;

mesh

areas = PropertyValue[{mesh,{2,All}},MeshCellMeasure];
max = Max[areas];
PropertyValue[{mesh,{2,All}},MeshCellStyle]=Map[Hue[#/max]&,areas];
mesh

Meshes from Points

Here is a function that generates points in 3D that are uniformly distributed on the surface of the unit sphere.

sup[n_]:=Module[{θ,cosϕ,sinϕ},
θ=RandomReal[2Pi,n];
cosϕ=RandomReal[{-1,1},n];
sinϕ=Sin[ArcCos[cosϕ]];
Transpose[{Cos[θ]sinϕ,Sin[θ]sinϕ,cosϕ}]];

A set of 47 such points -- you may want to experiment with the number.

points = sup[47];

Construct a MeshRegion from the points with no other cells

pm=MeshRegion[points,Point[Range[Length[points]]]]

HighlightMesh[pm,Labeled[0,"Index"]]

DelaunayMesh

DelaunayMesh does a Delaunay triangulation of a set of points.
VoronoiMesh constructs the dual of the Delaunay mesh.

dm=DelaunayMesh[points]

HighlightMesh[dm,{Style[3,Opacity[.1]],Labeled[0,"Index"]}]

How would you verify the Delaunay property?
Hint: Circumsphere.

Verifying the Delaunay property

testCell[Tetrahedron[vertices_],data_] := Not[Apply[Or,MemberQ[Circumsphere[vertices],points]]]

MeshPrimitives[dm,{3,1}]

testCell[%,points]

Map[testCell[#,points]&,MeshPrimitives[dm,{3,All}]]

ConvexHullMesh

This gives the boundary mesh for the convex hull of the given points.

ch = ConvexHullMesh[points]

See if you can get a good approximation of a sphere, a ball?

ConvexHullMesh[sup[1000]]

Meshes from Regions

There are 4 commands that do this:
DiscretizeRegion, DiscretizeGraphics
BoundaryDiscretizeRegion,BoundaryDiscretizeGraphics

The Graphics converters handle “special” regions that are effectively graphics primitives, like Disk, Ball, Cone, etc. (and also graphics in general)
and are typically faster when they work.

The Region converters will try harder to convert any region (RegionQ) into a mesh.

BoundaryDiscretizeGraphics[Ball[]]

iball = ImplicitRegion[x^2+y^2+z^2≤1,{x,y,z}];

BoundaryDiscretizeGraphics[iball]

BoundaryDiscretizeRegion[iball]

Note that the mesh generated is not nearly as good, but ...

BoundaryDiscretizeRegion[Ball[]]

All of the commands take bounds, e.g. the part of the ball in the first octant

BoundaryDiscretizeRegion[Ball[],{0,1}]

DIscretizeGraphics can handle graphics more generally too:

g = Plot3D[Sin[x y],{x,0,4},{y,0,4}]

DiscretizeGraphics[Normal[g]]

Some regions to try

Spend some time getting meshes from regions -- try your own favorites, or here are some to try

tr[1]= ImplicitRegion[Mod[x,7]==1,{{x,0,47}}];
tr[2]=ImplicitRegion[Sin[x^2]>=0,{{x,0,10}}];
tr[3]=Disk[{0,0},{1,3}];
tr[4]=Disk[{0,0},{3,1},{0,π}];
tr[5] = RegionUnion[tr[3],tr[4]];
tr[6] = RegionIntersection[tr[3],tr[4]];
tr[7]=ImplicitRegion[1≤x^2+y^2≤2,{x,y}];
tr[8]=ImplicitRegion[x^2+2y^2≥1 || x^2-y^2==.1,{x,y}];
tr[9,a_:0]:=
ImplicitRegion[(-1+x) (-3+2 x) (-x^2+y^2)-4 (-2 x+x^2+y^2)^2≥a,{x,y}]; 
tr[10,a_:0] := ImplicitRegion[-(1/2)+(((-3 x+4 (x^3)))^2)-(Sqrt[2]) (-3 x+4 (x^3)) (-3 y+4 (y^3))+(((-3 y+4 (y^3)))^2)<= a,{x,y}];tr[11] = RegionProduct[tr[2],tr[3]];
tr[12]=Cone[{{0,0,0},{1,1,1}}];

Some of these give different looking results from DIscretizeRegion and BoundaryDiscretizeRegion --why?
For the RegionProduct, try the product of the discretized arguments (same for RegionUnion, RegionIntersection)

A word of caution regarding meshes

Now that you know how to discretize regions, its tempting to go off the deep end and discretize everything you see, but, for many regions, it is better to keep the region in symbolic form up until the point that you need to do something explicit with it.  There are two important reasons for this:
1) Usually the symbolic form is a compact representation of the region -- the corresponding mesh may use much more memory
2) Once a region has been discretized, things done with the mesh can only be (at best) as accurate as the discretization.

Here is an example of number 2 above:

NIntegrate[1,{x,y}∈Disk[]]-π

NIntegrate[1,{x,y}∈DiscretizeRegion[Disk[],MaxCellMeasure->0.00001]]-π

Working with Meshes

As an example of working with meshes, lets write a simple (not very accurate) numerical integrator for 2D regions.

Here is a simple integration rule for a triangle given the values of the function at its 3 vertices {(Subscript[f, 1]),(Subscript[f, 2]),(Subscript[f, 3])}

(∫_T)f≃(1/3)(Subscript[A, T])((Subscript[f, 1])+(Subscript[f, 2])+(Subscript[f, 3]))
where (Subscript[A, T]) is the area of the triangle.

Write a function that will use this to approximate the integral of a function f over a triangulated mesh m.

tintegrate[f_,m_MeshRegion] := ...

For simplicity, you may assume that f is a pure function.

Try your function with

f[x_,y_]:= Sin[2π x^2]Cos[3π (y^2-x)];
Plot3D[f[x,y],{x,-1,1},{y,-1,1}]

On discretizations of the regions

Ω1=Rectangle[];

Ω2=Disk[];

See how quickly it converges to carefully computed values for the integrals:

i1=-0.07873298731206846`;
i2 = -0.018888805449916633`;

One solution

Clear[tintegrate];
tintegrate[f_,m_MeshRegion]:=Module[{v, fv, tv,ta},
v = MeshCoordinates[m];
fv = Apply[f,v,{1}];
tv = Part[MeshCells[m,2],All,1];
ta = PropertyValue[{m,2},MeshCellMeasure];
Total[MapThread[(Total[fv[[#1]]]#2/3)&,{tv,ta}]]
]

edata1 = Table[
{h = 2.^-k,
tintegrate[f,DiscretizeRegion[Ω1,MaxCellMeasure->h]]-i1},{k,3,16}]

ListPlot[RealExponent[edata1]]

Fit[RealExponent[edata1],{1,x},x]

edata2 = Table[
{h = 2.^-k,
tintegrate[f,DiscretizeRegion[Ω2,MaxCellMeasure->h]]-i2},{k,3,16}]

ListPlot[RealExponent[edata2]]

Fit[RealExponent[edata2],{1,x},x]

Importing Meshes

This is a text cell.

Ω=BoundaryDiscretizeGraphics[Import["gear.stl"]]

RegionBounds[Ω]

RegionMeasure[Ω]

NIntegrate[1,{x,y,z}∈Ω]

Solving a PDE

Solve the heat equation where the inner (axle hole) of the gear is given a particular temperature.

First get the inner circle radius:

vertices = MeshCoordinates[Ω];
radii = Map[Norm[#[[{1,2}]]-{100,100}]&,vertices];
ListPlot[radii]

{Min[radii],Max[radii]}

r1 = 25;
r2=100;usol = NDSolveValue[{Laplacian[u[x,y,z],{x,y,z}]==0,DirichletCondition[u[x,y,z]==1,(x-100)^2+(y-100)^2≤ r1^2],DirichletCondition[u[x,y,z]==0,(x-100)^2+(y-100)^2≥ r2^2]},u,{x,y,z}∈Ω,Method->{"FiniteElement","MeshOptions"->{"MeshOrder"->1}}]

NDSolve actually creates a special form of mesh, an ElementMesh that is very similar to a MeshRegion, but is based on elements that work with the finite element method directly.

emesh = usol["ElementMesh"]

It can be converted to a MeshRegion quite easily:

mr = MeshRegion[emesh]

Can you use the mesh properties to make a visualization of solution on the boundary?

bcent = PropertyValue[{mr,{2,"Boundary"}},MeshCellCentroid];

bval = Apply[usol,bcent,{1}];
{Min[bval],Max[bval]}

cfun=ColorData["TemperatureMap"]

PropertyValue[{mr,{2,"Boundary"}},MeshCellStyle]=Map[cfun,bval];

mr

Cheating a bit (functions like this are described in the finite element documentation)

NDSolve`FEM`ElementMeshSurfacePlot3D[usol]

Summary

In this workshop, we’ve just scratched the boundary of what it possible with mesh based regions and region discretization.
There are lots of examples in the documentation.