网格划分方法

• Larger desired lengths than actual lengths gives “internal pressure” • Find force equilibrium using Forward Euler pn+1 = pn + ∆tF (pn )
8
Boundary Constraints
23
Gradient Limiting
• The Gradient Limiting Equation: ∂h + |∇h| = min(|∇h|, Gc ), h(0) = h0 ∂t • Analytical solution (convex domains, use Lax’ theorem): h(x) = min(h0 (y ) + Gc x − y ).
1. Distribute points using size function h(x, y ), reject points outside 2. Obtain topology by Delaunay triangulation 3. Find force equilibrium
7
Force Equilibrium
Cartesian
Octree
Unstructured
16
Curvature Adaptation
• On the boundaries h(x) ≤ 1/K |κ(x)|, resolution K • Curvature κ computed explicitly or using distance function
Contours of φ(x, y ) and shock y
j+1 j i−2 i−1 i i+1 i+2 i+3 x
φ
φ(x, j ), p1 (x), p2 (x) p2 p1
φmax
i−2 i−1
i
i+1 i+2 i+3
x
20
Feature Size Adaptation - Example
1
Topics
1. Introduction 2. The New Meshing Algorithm 3. Automatic Mesh Size Control 4. Applications
2
Mesh Generation
• Given a geometry, determine node points and element connectivity • Applications:
1. Identify medial axis by shock detection in φ(x) (skeletonization) 2. Solve ∇φMA
· ∇lfs = 0 with lfs(x) = 2φ(x) on medial axis
19
Shock Detection in Distance Function φ(x)
y
• PDE-based formulation gives straight-forward extension to:
– Higher order – Unstructured meshes – Space and solution dependent Gc (x, h)
• Solve in O(n log n) time using modified Fast Marching Method
5
Topics
1. Introduction 2. The New Meshing Algorithm 3. Automatic Mesh Size Control 4. Applications
6
The New Meshing Algorithm
Distribute points Triangulate Force equilibrium
11
Mesh Quality
Delaunay Refinement with Laplacian Smoothing
150
# Elements
100
50
0 0.7 150
0.8
0.9
1
Force Equilibrium by distmesh2d
# Elements
100
50
0 0.7
0.8 0.9 Element Quality
• Piecewise linear force function (current length , desired length k ( − ) if < , 0 0 f ( , 0) = 0 if ≥ 0 . •
0 0 ):
= C1 h(x, y ) for non-uniform element sizes
• Use distance function to project points to boundary • Points can move tangentially along curve
Updated Node Location = (x,y) − ∇d(x,y)⋅d(x,y) (x,y)
x,y
– Better treatment of boundaries, graded meshes more efficient – FEM/FVM/BEM methods better suited
• Distance functions sometimes more natural
– Image processing, MRI scans, biomechanics ´ • Easier to implement and manipulate than B-rep and Bezier patches
4
Why Create Meshes on Level Sets?
• Level Set Methods superior for interface propagation with large
deformations and topology changes
• But: Unstructured grids better for the physical problem
1
12
Future Work
• Performance
– Local Delaunay updates, 2-D easy, 3-D studied by [Devillers] – Hierarchical solver – Implicit instead of Forward Euler
• Robustness
17
Feature Size Adaptation
• h(x) ≤ lfs(x)/R • The local feature size lfs is “width” of geometry • Ruppert: lfs(x) is “the larger distance from x to the closest two
d( x,y )
−∇
d(
)
9
MATLAB Demo
10
4-D Hypersphere
• d = r − 1 with r =
4 i=1ቤተ መጻሕፍቲ ባይዱ
x2 i
• h0 = 0.2 gives 3, 458 nodes and 59, 222 elements • Mesh volume V4 = 4.74 (expected π 2 /2 ≈ 4.93) • Hyper-surface area S4 = 19.2 (expected 2π 2 ≈ 19.7) • Poisson’s equation −∇2 u = 1, bnd cond’s u = 0|r=1 . Analytical solution u = (1 − r 2 )/8, error e ∞ = 7.9 · 10−4 .
non-adjacent polytopes [of the boundary]”
18
Feature Size Adaptation
• Alternative definition: lfs(x) is equal to the smallest distance between the boundary and the point on the medial axis closest to x • Algorithm:
– Numerical Solution of PDEs (FEM, FVM, BEM), Interpolation – Computer Graphics, Visualization
• Delaunay refinement, Advancing Front, Octree
– Explicit boundary representations
13
Topics
1. Introduction 2. The New Meshing Algorithm 3. Automatic Mesh Size Control 4. Applications
14
Automatic Mesh Size Control
• Important preprocessing step for many mesh generators • Find mesh size function h(x) considering 1. Curvature, 2. Local feature
3
The Level Set Method [Osher, Sethian]
• Represent geometry by signed distance function φ(x) • Discretize φ on Cartesian, unstructured, or Octree grid • Propagate boundary with PDEs, e.g. φt + F |∇φ| = 0, φt + v · ∇φ = 0 • Handles very large deformations • Higher stability than explicit methods • Fast Marching Method: Solve |∇φ| = 1 in O(n log n) time
size, 3. Non-geometric data (user, numerical), and 4. Grading
15
Background Grids
• Discretize mesh size function h(x) on a coarse background grid • Mesh from previous iteration for adaptive FEM solvers
– Repeated discrete refinements [Borouchaki et al] – Balanced Octrees [Blacker, Stephenson]
• New continuous formulation:
– Steady-state solution to Hamilton-Jacobi equation
– Respect boundaries – Control logic and termination criterion involving element qualities
• Other ideas
– Quad and Block meshing – Force-based sliver exudation
21
Curvature + Feature Size Adaptation
22
Gradient Limiting
• Max size ratio G between elements (hi ≤ Ghj for neighboring i, j ) • Continuous analogue: |∇h(x)| ≤ G − 1 ≡ Gc • Traditional methods:
Mesh Generation using Level Set Methods and PDE-Based Size Control
Per-Olof Persson (persson@) Department of Mathematics, MIT
MIT EECS/ME/AA, March 15, 2004