1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
|
"""
"""
import random
import numpy as np
import perfplot
from scipy.spatial import Delaunay
import meshplex
def setup(n):
radius = 1.0
k = np.arange(n)
boundary_pts = radius * np.column_stack(
[np.cos(2 * np.pi * k / n), np.sin(2 * np.pi * k / n)]
)
# Compute the number of interior points such that all triangles can be somewhat
# equilateral.
edge_length = 2 * np.pi * radius / n
domain_area = np.pi - n * (radius**2 / 2 * (edge_length - np.sin(edge_length)))
cell_area = np.sqrt(3) / 4 * edge_length**2
target_num_cells = domain_area / cell_area
# Euler:
# 2 * num_points - num_boundary_edges - 2 = num_cells
# <=>
# num_interior_points ~= 0.5 * (num_cells + num_boundary_edges) + 1 - num_boundary_points
m = int(0.5 * (target_num_cells + n) + 1 - n)
# generate random points in circle;
# <https://mathworld.wolfram.com/DiskPointPicking.html>
for seed in range(0, 255):
np.random.seed(seed)
r = np.random.rand(m)
alpha = 2 * np.pi * np.random.rand(m)
interior_pts = np.column_stack(
[np.sqrt(r) * np.cos(alpha), np.sqrt(r) * np.sin(alpha)]
)
pts = np.concatenate([boundary_pts, interior_pts])
tri = Delaunay(pts)
# Make sure there are exactly `n` boundary points
mesh0 = meshplex.MeshTri(pts, tri.simplices)
mesh1 = meshplex.MeshTri(pts, tri.simplices)
if np.sum(mesh0.is_boundary_point) == n:
break
mesh0.create_edges()
mesh1.create_edges()
num_interior_edges = np.sum(mesh0.is_interior_edge)
idx = random.sample(range(num_interior_edges), n // 10)
print(num_interior_edges, len(idx), len(idx) / num_interior_edges)
# # move interior points a little bit such that we have edges to flip
# max_step = np.min(mesh.cell_inradius) / 2
# mesh.points = mesh.points + max_step * np.random.rand(*mesh.points.shape)
# print(mesh.num_delaunay_violations)
return mesh0, mesh1, idx
def flip_old(data):
mesh0, mesh1, idx = data
mesh0.flip_interior_edges_old(idx)
def flip_new(data):
mesh0, mesh1, idx = data
mesh1.flip_interior_edges(idx)
perfplot.show(
setup=setup,
kernels=[flip_old, flip_new],
n_range=[2**k for k in range(5, 13)],
equality_check=None,
# set target time to 0 to avoid more than one repetition
target_time_per_measurement=0.0,
)
|
