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,
)
|