v0.14.0
sdf_wavy_2d.py
Go to the documentation of this file.
1 import numpy as np
2 from scipy.optimize import brentq
3
4 A = 0.0002
5 wave_len = 2
6 w = 2. * np.pi / wave_len
7 ind = 0.00159624
8
9 class WavySurface:
10  def F(p, x, y, t):
11  q = A * (1. - np.cos(w * p)) - ind * t
12  lag = 2. * (y - q)
13  return 2. * (p - x) - lag * A * w * np.sin(w * p)
14
15  def dist(x, y, a, b, t):
16  try:
17  p0 = brentq(WavySurface.F, a, b, args=(x, y, t))
18  q0 = A * (1. - np.cos(w * p0)) - ind * t
19  return np.sqrt((p0 - x)**2 + (q0 - y)**2)
20  except ValueError:
21  return 1e12
22
23  def sDF(x, y,z, t):
24  y0 = A * (1. - np.cos(w * x)) - ind * t
25  if np.any(np.abs(y0 - y)) < 1e-12:
26  return 0.0
27
28  x = np.abs(x) % wave_len
29  x = np.where(x > wave_len / 2, wave_len - x, x)
30
31  eps = 1e-12
32
33  d1 = np.vectorize(WavySurface.dist)(x, y, eps, wave_len / 2 - eps, t)
34  wave_y = A * (1 - np.cos(w * x)) - ind * t
35  sgn_d1 = np.sign(wave_y - y)
36
37  wave_y_d2 = -ind * t
38  sgn_d2 = np.sign(wave_y_d2 - y)
39  d2 = np.sqrt(x**2 + (wave_y_d2 - y)**2) # when xp to be 0
40
41  wave_y_d3 = A * (1 - np.cos(w * (wave_len / 2))) - ind * t
42  sgn_d3 = np.sign(wave_y_d3 - y)
43  d3 = np.sqrt((wave_len / 2 - x)**2 + (A * (1 - np.cos(w * (wave_len / 2))) - ind * t - y)**2) # when xp to be wave_len/2
44
45  d_min = np.minimum(d2, d3)
46  d_final = np.where(d1 < d_min, d1, np.where(d2 < d3, d2, d3))
47  sgn_final = np.where(d1 < d_min, sgn_d1, np.where(d2 < d3, sgn_d2, sgn_d3))
48  sdf = sgn_final * d_final
49  return sdf
50
51
52  def gradSDF(x, y, z, t):
53  delta = 1e-6
54  df_dx = ((WavySurface.sDF(x + delta, y, z, t) - WavySurface.sDF(x - delta, y, z, t)) / (2 * delta)).reshape((-1, 1))
55  df_dy = ((WavySurface.sDF(x, y + delta,z, t) - WavySurface.sDF(x, y - delta,z, t)) / (2 * delta)).reshape((-1, 1))
56  df_dz = np.zeros_like(df_dy).reshape((-1, 1))
57  return np.hstack([df_dx, df_dy, df_dz])
58
59  def hessSDF(x, y, z, t):
60  delta = 1e-6
61  df2_dx2 = (WavySurface.sDF(x + delta, y,z,t) - 2 * WavySurface.sDF(x, y,z,t) + WavySurface.sDF(x - delta, y,z, t)) / (delta ** 2)
62  df2_dy2 = (WavySurface.sDF(x, y + delta,z,t) - 2 * WavySurface.sDF(x, y,z,t) + WavySurface.sDF(x, y - delta,z, t)) / (delta ** 2)
63  df2_dxdy = (WavySurface.sDF(x + delta, y + delta,z,t) - WavySurface.sDF(x + delta, y - delta, z, t) - WavySurface.sDF(x - delta, y + delta,z, t) + WavySurface.sDF(x - delta, y - delta,z,t)) / (4 * delta ** 2)
64
65  return np.hstack([df2_dx2.reshape((-1,1)), df2_dxdy.reshape((-1,1)), np.zeros_like(df2_dy2).reshape((-1,1)),df2_dy2.reshape((-1,1)), np.zeros_like(df2_dy2).reshape((-1,1)), np.zeros_like(df2_dy2).reshape((-1,1))])
66
67 def sdf(delta_t, t, x, y, z, tx, ty, tz, block_id):
68  return WavySurface.sDF(x, y, z, t)
69
70 def grad_sdf(delta_t, t, x, y, z, tx, ty, tz, block_id):
72
73 def hess_sdf(delta_t, t, x, y, z, tx, ty, tz, block_id):
74  return WavySurface.hessSDF(x, y,z, t)
sdf_wavy_2d.hess_sdf
def hess_sdf(delta_t, t, x, y, z, tx, ty, tz, block_id)
Definition: sdf_wavy_2d.py:73
sdf_wavy_2d.sdf
def sdf(delta_t, t, x, y, z, tx, ty, tz, block_id)
Definition: sdf_wavy_2d.py:67
Definition: sdf_wavy_2d.py:52
sdf_wavy_2d.WavySurface.sDF
def sDF(x, y, z, t)
Definition: sdf_wavy_2d.py:23
sdf_wavy_2d.WavySurface.F
def F(p, x, y, t)
Definition: sdf_wavy_2d.py:10
sdf_wavy_2d.WavySurface
Definition: sdf_wavy_2d.py:9
sdf_wavy_2d.WavySurface.dist
def dist(x, y, a, b, t)
Definition: sdf_wavy_2d.py:15