Boundary Conforming Unstructured Mesh Generation using Advancing Fronts

1
Boundary Conforming Unstructured Mesh Generation
using Advancing Fronts
By Ryan Chilton Advisor Robert Lee
The ElectroScience Laboratory Department of
Electrical Engineering The Ohio State University
2
Mesh Generation in Computational EM
Both IE and PDE methods require geometry
discretizations (meshes).
3D isoparametric surface mesh. Suitable for
general SIEs.
3D linear tetrahedra. Suitable for VIE, FEM.
2D linear triangle mesh. Suitable for 2D FEM
We focus here on Boundary Conforming meshes,
leaning towards FE
1.
A pre-discretized boundary definition is given.
Respect it.
Tile and cover the domain with elements of
prescribed size.
2.
3
Relevant Element Quality Metrics
A comprehensive metric of element shape is given
by the ratio
2D
3D
Poorer elements have more extreme (smaller)
ratios
Dont Forget!
Element size is also important.
Governs modeling accuracy for all.
Governs dispersion relation in PDE.
Typically quoted in elements per
Very small ratio, poor quality.
Large ratio, high quality.
4
Delaunay Triangulation1
A mesh is Delaunay iff it has the empty
circle(sphere) property.
Definition
For each triangle T(A,B,C), the circle
circumscribing T contains no nodes other than A,
B, and C. 2
A (the) delaunay triangulation of a set of points.
Alternate, non-Delaunay mesh.
Delaunay triangulations can be excellent for
numerical methods
For a given point set, Delaunay maxes the min
angle. (!) 3
Qualifier Delaunay tetrahedralization does not
have this property.
For n points, finding the delaunay triangulation
is O(n log n) 4
5
Constrained Delaunay Triangulation
So, that means Delaunay is awesome right?
Fast.
Does not conform naturally to boundaries.
Correct.
Presumes Nodes is known a priori!
Provably optimum
Incremental insertion prohibits QHull.
Delaunay triangulation does not solve the stated
problem.
A relaxed version of the definition of inside
the circle
Definition
New Steiner Points
Node D is only inside the circumcircle of
T(A,B,C) if D is visible from T.
Boris Nikolaevich Delone (1890-1980)
and allowing Steiner Points
permits efficient boundary Constrained Delaunay
Triangulation 56.
6
Greedy Strategies
Advancing Fronts (AF) is a kind of Greedy
Algorithm 7. Greedy algorithms exploit an
optimal substructure property
Lemma
A mesh M is high quality iff every submesh of M
is high quality. (Applies all the way down to
single element meshes!)
Greedy algorithm correctness requires the
greedy-choice property
Definition
This criterion is not met.
A globally optimal solution can be arrived at by
making a locally optimal (greedy) choice.
(More on this later!)
7
Advancing Fronts
Definitions Galore!
So, what is a front?
A front is a set of non-crossing orientable
edges where each node has is same in-degree as
out-degree.
Orientable Edges have direction. The edge from
I to J is not the same as the edge from J to I.
Node A point location, (x,y) ordered pair.
Edge A connection between two nodes.
In-degree the in-degree of a node is the of
incoming edges.
Out-degree the out-degree of a node is the of
outgoing edges.
Here are some fronts
These are not fronts
Bodini MT
CopperPlate
SHOWCARD
Theyre not Fronts, theyre Fonts! (Ha, Ha)
8
Advancing Fronts
The Greedy Strategy While the front edge set is
nonempty,
Pick a front edge E. Create equilateral
triangle T, growing towards Remove E from the
front, add the two new edges of T.
1. 2. 3.
this is why edges are orientable.
Demonstration
Inward normal.
New front edges.
Selected edge.
New triangle.
After more iterations.
Red edges still satisfy front definition.
Loop Invariant
The front is the boundary between meshed and
unmeshed geometry.
9
Pathological Outcomes of Greedy Strategy
Greedy-choice property is not satisfied. An
optimal triangle T might
Opt
A bad triangle lurking beyond the horizon.
Crossing front edges, violates front property.
Opt
A
A
B
Selected edge.
B
Violate the front definition.
Leave a bad choice for later.
Small zone around optimum position.
Snap to existing node.
C
A
B
Avoid these by implementing penalties and a
smidge of look-ahead.
10
Revised Greedy Strategy
Algorithm
while (Front is nonempty) Select the
shortest Front Edge (A,B). Compute the optimum
node placement, Opt, for Edge(A,B) delta
length(Edge(A,B)) Define the set of candidate
nodes as empty set. (C ) while (C is
empty) C All nodes which are within a
ball of size delta to Opt. For each candidate
node Ci If Edge (A,Ci) or Edge (B,Ci)
cross a Front Edges, remove Ci from C If
Ci is A, remove Ci from C If Ci is B,
remove Ci from C. If Triangle (A,B,Ci)
has negative or zero area, remove Ci from C
If (C is empty) If Edge(A,Opt) and
Edge(B,Opt) cross no Front Edges, insert Opt into
C Else delta delta2 For each
candidate node Ci Compute quality scores for
Triangle (A,B,Ci) Create and store the
Triangle, Ti, which has the maximum quality
score. Remove Edge (A,B) from the Front. If
Edge(A,Ci) is in the Front, remove Edge(A,Ci)
from the Front. Else add Edge (A,Ci) to the
Front. If Edge(B,Ci) is in the Front,
remove Edge(B,Ci) from the Front. Else add
Edge (B,Ci) to the Front.
Greed
Lookahead
Penalties
Greed
Advance Maintain Front Invariant
11
Advancing Fronts in Action
Meshing a simple closed domain.
Meshing a domain with a hole.
The 2D algorithm is correct, but not optimal.
Domain is partitioned and covered. Each element
valid.
Not every element is high quality. Must cleanup.
12
Mesh Improvement Topology Changes
Local correction strategies for element
improvement exist. 8
1.
3.
Algorithms
2.
Offending Edge
Scope of operation.
13
Mesh Improvement Smoothing
Smoothing is moving the nodes without changing
their connectivity.
Algorithm
For each movable node, move it to the centroid
of its incidence polygon.
Note how similar this is to solving Laplaces
equation using Jacobis or Gauss/Seidel Method.
More sophisticated update algorithms exist as
well. 9
Massive Nodes.
Springy Edges.
1.
2.
Short edges try to expand.
Beware, some schemes exceed the difficulty of our
PDEs!
Long edges in tension.
14
Sample meshes after improvement (1/4)
272 Triangles 447 Edges 175 Nodes 1.59 sec CPU
Square with circle target.
15
Sample meshes after improvement (2/4)
1956 Triangles 3030 Edges 1075 Nodes 11.83 sec CPU
Rectangle with ogival target.
16
Sample meshes after improvement (3/4)
6041 Triangles 9201 Edges 3160 Nodes 57.7 sec CPU
Two nearby Nodes.
Degenerate Triangle
GIGO
A Tiny Edge
A Tiny Angle
Rectangle with NACA airfoil target.
17
Sample meshes after improvement (4/4)
Mom
4743 Triangles 7230 Edges 2487 Nodes 50.38 sec
CPU
Scooter
My old Kentucky home.
University of Kentucky
18
Computational Complexity (1/3)
Presuming the algorithm terminates, well count
operations.
(In 2D, this assumption is provable. Not so easy
in 3D)
Definition
Number of Front Edges when the tth triangle is
being created.
F(t)
Termination Assumption
while (Front is nonempty) Select the
shortest Front Edge (A,B). Compute the optimum
node placement, Opt, for Edge(A,B). Candidacy
set C Opt. while (No valid candidate has
been found) Increase candidacy ball.
Compute the candidacy node set. For each
candidate node Ci If Triangle (A,B,Ci)
cross a Front Edge, remove Ci. For each
candidate node Ci Compute quality scores
for Triangle (A,B,Ci) Create and store the
triangle which has the maximum quality score.
Update the Front.
O(T) O(log F(t)) O(1) O(1) O(1)
O(1) O(F(t)) O(1) O(F(t))
O(1) O(1) O(1) O(log F(t))
Leading Order Term
(Enforcing front-crossing penalty)
Total running time is O(TF(t)). But F(t) is
tricky to bound tightly.
19
Computational Complexity (2/3)
Intuition Eulers formula suggest O(T) O(N)
O(E).
N Nodes T Triangles E Edges
Leonhard Euler 1707 - 1783
But is F(t) O(E), thus O(T)? Empirical evidence
suggests its smaller
Square E-Shape Kentucky
Square Meshes (various sizes)
20
Computational Complexity (3/3)
How does Favg(T) compare to T?
Square Meshes (various sizes)
Square Meshes (various sizes)
Favg(T) is probably O(T0.5). Think about
dimensional analysis
T is an area metric. AreaO(T) F is a perimeter
metric. PerimeterO(F)
The total running time for AF would be
O(T1.5)O(N1.5).
21
Acceleration via a DC implementation?
Largest complexity term is wasteful front
management
Divide!
Conquer!
Sun Tzu famous proponent of divide and conquer.
Adding T(a,b,c) requires intersection test
against entire front, O(E)
Cannot check only nearby edges determining them
is O(E)
Partition front set about a cutting line?
Almost think big.
Partition about a cutting
line!
entire problem
22
Divide and conquer strategy. (1/3)
1.
2.
1.
Original problem description, consists of all
front edges.
2.
Cut the domain, introducing non-crossable,
non-meshable boundaries. Splits the front edges
into two disjoint sets.
23
Divide and conquer strategy. (2/3)
Cut edge.
Exposed edges.
New front.
3.
4.
3.
Mesh each sub domain independently.
4.
Sew up by restarting the algorithm with new
fronts, which are composed of cut exposed edges
Overall, this is the same basic concept as
merge-sort.
24
Divide and conquer strategy. (3/3)
Why stop at one level partitioning? For huge
domains, keep cutting!
The tail of the recursion a domain falls below
an area threshold.
Partitioning order.
Grow out until box encountered.
Seed triangle. Introduces 3 new front edges.
Nice work!
Domains without any front edges? Seed a new
triangle grow out.
25
Divide conquer results (1/4)
Test platform 2.4GHz P4, 2GB g v3.4.4 g
Iterative
1707 Triangles 2.89 sec
Divide Conquer
1583 Triangles 1.08 sec
Color denotes membership in a common partition.
26
Divide conquer results. (2/4)
Inset (next pg)
A few slivers (lt10)
Comparable quality, superior speed.
27
Divide conquer results. (3/4)
Empty zone.
Merge areas.
Quality in merge areas improves if a small empty
zone is left between leaf domains.
Zone size controls tradeoff between quality and
merge complexity.
28
Divide conquer results. (4/4)
x
DC is meshing much faster than the iterative
implementation
Each data point is a mesh.
106.37 triangles (2.4M) 103.33 sec
(3553) (Scalability, wow!)
1.5 decades
1 decade
1 decade
Disparity in slopes suggests improvement in
asymptotic complexity.
Probably running in O(n log2n), not O(n)
Consistent with previous O(n1.5) conjecture.
29
Parallelization strategy. (1/2)
Note that each of the leaf meshes can be solved
simultaneously
Task assignment.
Domain partitioning sequence.
However, the sewing operations introduce
workflow dependencies

Directed acyclic graph (tree) of dependence
relationships.
For single CPU DC, the call stack of the
recursion handles dependency.
Multi-CPU algorithm must emulate this behavior
with a scheduler.
30
Parallelization strategy. (2/2)
Use a workpool model of job allocation, master
/ slave architecture.
Master
Partition geometry. Issue jobs from
workpool. Maintain dependency tree.
Communication
Job descriptions
Completion messages.
Slaves
Perform merge operations.
Essentially, these represent a mechanism for
remote function calls.
Mesh leaves.
31
Early parallel observations.
Test platform OSC Itanium2 mpiCC v1.2.6
IO is terribly expensive, approximately 50 of
total effort.
Linear speedup regime. (JgtgtP)
Serialized regime. (JP)
JJobs (mesh domains) PProcess (meshers)
Low overhead. Comm. is O(Perimeter) or O(1)! ,
work is ?(Area)
Like any respectable tyrant, the master does
little work!
This Master/Slave architecture limits speedup to
at most a factor of P-1. (as opposed to P)
32
Parallel scalability. (1/4)
Test platform OSC Itanium2 (IO suppressed)
Each o is a median of 3 trials (meshes).
Linear speedup.
Maximum expectable speedup.
Examine this trial, next page.
Increasing problem size.
Clear saturation effect, more processors ? faster
turnaround.
Larger problems dont fix saturation problem?
(Red flag!)
33
Parallel scalability - postmortem. (2/4)
So many processes that master cannot distribute a
job to each one before the first slave finishes?
Want to work, but seldom issued a
job. (Disastrous!)
First slave always working.
Master always communicating.
Appears that sub-problem granularity is set too
small.
34
Parallel scalability. (3/4)
Control granularity by adjusting minimum area
cutoff of DC tail.
Linear speedup regime. (JgtgtP)
Serialized regime. (JP)
Previous granularity Domains lt 16? by 16? are
dispatched
Current granularity Domains lt 64? by 64? are
dispatched
Each process now given work to do.
Dependency
Master no longer in a frenzy of comm.
Lower wall clock.
35
Parallel scalability. (4/4)
Repeat the previous measurements with larger
process granularity.
Speedup?
Much better!
35 minutes, 53 seconds on P4 2.4GHz g test
rig. 38 seconds on 16 CPU Itanium 2 cluster
with mpiCC Estimate 4 days, 1.8 hours for
iterative algorithm.
Linear speedup.
Maximum expectable speedup.
Previous measurements
Ideally, task granularity should be a function of
problem size and P.
36
Extending to 3D mesh generation
In some ways, Advancing Fronts extends naturally
to 3D
Orientable edges become orientable faces.
Front still the boundary of meshed unmeshed
geometry.
EMCC Almond Profile
EMCC Almond Target
Local, topology changing correction strategies
exist in 3D.
If the presented AF alg. can be ported to 3D, it
is O(N5/3)
37
Every silver lining has a cloud.
Suggested 3D algorithm not guaranteed to halt on
all inputs.
(A problem this big deserves a few extra Xs)
Evil takes many forms.
The Schonhardt Prism
1.
Equilateral prism mesh.
2.
Twist the top to deform.
3.
Call AF, loop ad infinitum.
The Schonhardt Prism is one such form.
Maybe you wont mesh one explicitly as an input
but it (or an isomorphism of it) might pop up
when fronts collide!
Workaround Backtracking. Goodbye efficiency,
hello poor quality.
Definition
FYI even deciding if AF can tetrahedralize a
given hull (ie, failing) is NP-hard.
NP-hard Provably Evil.
Reduction to 3D art gallery problem.
38
References
Jin-Fa Lee Dyczij-Edlinger, R. Automatic mesh
generation using a modified Delaunay tessellation
Antennas and Propagation Magazine, IEEE Volume
39, Issue 1, Feb. 1997 Page(s)34 - 45
Jonathan Richard Shewchuk, Triangle Engineering
a 2D Quality Mesh Generator and Delaunay
Triangulator.Applied Computational Geometry
Towards Geometric Engineering'', pages 203-222,
Springer-Verlag, Berlin, May 1996. Charles L.
Lawson. Software for C1 Surface Interpolation.
Mathematical Software III, pp 161-194. Academic
Press, New York, 1977. Joseph ORourke,
Computational Geometry in C, Cambridge University
Press, NY NY, 1994. Jonathan Richard Shewchuk,
Delaunay Refinement Algorithms for Triangular
Mesh Generation, Computational Geometry Theory
and Applications 22(1-3)21-74, May 2002 Chew,
L. P. 1987. Constrained Delaunay triangulations.
In Proceedings of the Third Annual Symposium on
Computational Geometry (Waterloo, Ontario,
Canada, June 08 - 10, 1987). D. Soule, Ed. SCG
'87. ACM Press, New York, NY, 215-222. T.
Cormen, C. Leiserson, R. Rivest, Introduction to
Algorithms, MIT Press, Cambridge Massachusetts,
1990. Per-Olaf Persson, Mesh Generation for
Implicit Geometries. MIT University Press,
Cambridge Massachusetts, 1997. Freitag L. A.,
Plassmann P. E., Local Optimization-based
simplicial mesh untangling and improvement, 1997.
1
2
3
4
5
6
7
8
9