In [30]:
%load_ext autoreload
import numpy as np
import functools
%aimport FE
%aimport common_meshes
import matplotlib.pyplot as plt
import triangle
import pandas as pd
import matplotlib.cm as cm
%matplotlib inline
print(plt.style.available)
plt.style.use('seaborn-paper')
In [2]:
def cartesian_product(arrays):
broadcastable = np.ix_(*arrays)
broadcasted = np.broadcast_arrays(*broadcastable)
rows, cols = functools.reduce(np.multiply, broadcasted[0].shape), len(broadcasted)
out = np.empty(rows * cols, dtype=broadcasted[0].dtype)
start, end = 0, rows
for a in broadcasted:
out[start:end] = a.reshape(-1)
start, end = end, end + rows
return out.reshape(cols, rows).T
In [5]:
def is_in_ellipse(L, W, lam, pt):
x,y = pt
return int((x*np.pi/L)**2 + (y*np.pi/W)**2 < lam)
def count(L, W, lam, dirichlet):
"""Compute the eigenvalue counting function of a rectangle"""
# if dirichlet, start all counting ranges at 1. otherwise, start at 0.
start = int(dirichlet)
max_x = np.floor(L*np.sqrt(lam)/np.pi)
max_y = np.floor(W*np.sqrt(lam)/np.pi)
x = np.arange(start, max_x+1, 1)
y = np.arange(start, max_y+1, 1)
grid = cartesian_product([x,y])
def test(pt):
return is_in_ellipse(L, W, lam, pt)
if not list(grid):
return 0
else:
return np.sum(np.apply_along_axis(test, 1, grid))
count(5, 3, 10, True)
Out[5]:
In [65]:
first_eigenvalues = np.array(
sorted([(np.pi*m/1.)**2 + (np.pi*n/1.)**2 for m in range(10) for n in range(25)])[:25]
)
In [43]:
print(first_eigenvalues)
In [78]:
skeleton = {"vertices": np.array([[0., 0.], [1., 0.], [1., 1.], [0., 1.]]),
"segments": np.array([[0, 1], [1, 2], [2, 3], [3, 0]])}
meshes = {}
fineness = ["0.00005", "0.00007", "0.0001", "0.0005"]
for a in fineness:
meshes[str(a)] = triangle.triangulate(skeleton, 'pqa'+str(a))
In [79]:
%%time
eigvals = {}
for a in fineness:
eigvals[str(a)] = FE.findEigs(meshes[str(a)], 25)[0]
In [80]:
eigval_array = np.array([first_eigenvalues] + [eigvals[str(a)] for a in fineness]).T
eigval_df = pd.DataFrame(eigval_array, columns=["ground_truth"] + [str(a) for a in fineness])
In [81]:
eigval_df
Out[81]:
In [82]:
colors = cm.viridis(np.linspace(0, 1, len(fineness)+1))
for y, c in zip(eigval_df.columns, colors):
print(c)
plt.scatter(eigval_df.index, eigval_df[y], label=y, color=c)
plt.legend(loc='upper left')
plt.show()
In [83]:
colors = cm.viridis(np.linspace(0, 1, len(fineness)+1))
for y, c in zip(eigval_df.columns, colors):
print(c)
plt.scatter(eigval_df.index, (eigval_df[y] - eigval_df["ground_truth"])/eigval_df["ground_truth"],
label=y, color=c)
plt.legend(loc='upper left')
plt.show()
In [ ]: