"""Parametrizations of quadrics.
A note on numerical accuracy: Due to the exponential growth of cosh and sinh,
the relation ``cosh(t) ** 2 - sinh(t) ** 2 = 1`` Can break earlier than perhaps
expected. Here is the maximum error of this relation for some randomly sampled
parameter values, which should be taken into account when choosing tolerances:
+-----------------+---------------+-------------------------+
| Parameter value | Approx. error | Approx. size of cosh(t) |
+-----------------+---------------+-------------------------+
| 0 < t < 1 | 1e-15 | 1e2 |
+-----------------+---------------+-------------------------+
| 1 < t < 2 | 1e-15 | 1e2 |
+-----------------+---------------+-------------------------+
| 2 < t < 3 | 1e-14 | 1e3 |
+-----------------+---------------+-------------------------+
| 3 < t < 4 | 1e-13 | 1e3 |
+-----------------+---------------+-------------------------+
| 4 < t < 5 | 1e-12 | 1e4 |
+-----------------+---------------+-------------------------+
| 5 < t < 6 | 1e-11 | 1e4 |
+-----------------+---------------+-------------------------+
| 6 < t < 7 | 1e-10 | 1e4 |
+-----------------+---------------+-------------------------+
| 7 < t < 8 | 1e-9 | 1e5 |
+-----------------+---------------+-------------------------+
| 8 < t < 9 | 1e-9 | 1e5 |
+-----------------+---------------+-------------------------+
| 9 < t < 10 | 1e-8 | 1e6 |
+-----------------+---------------+-------------------------+
| 10 < t < 11 | 1e-7 | 1e6 |
+-----------------+---------------+-------------------------+
| 11 < t < 12 | 1e-6 | 1e7 |
+-----------------+---------------+-------------------------+
| 12 < t < 13 | 1e-5 | 1e7 |
+-----------------+---------------+-------------------------+
| 13 < t < 14 | 1e-4 | 1e8 |
+-----------------+---------------+-------------------------+
| 14 < t < 15 | 1e-3 | 1e8 |
+-----------------+---------------+-------------------------+
| 15 < t < 16 | 1e-3 | 1e8 |
+-----------------+---------------+-------------------------+
| 16 < t < 17 | 1e-2 | 1e9 |
+-----------------+---------------+-------------------------+
| 17 < t < 18 | 1e-1 | 1e9 |
+-----------------+---------------+-------------------------+
| 18 < t < 19 | 1e0 | 1e10 |
+-----------------+---------------+-------------------------+
| 19 < t < 20 | 1e1 | 1e10 |
+-----------------+---------------+-------------------------+
"""
import numpy as np
import ddg.datastructures.nets.net as net
import ddg.datastructures.nets.utils as nutils
# ===========
# 0D quadrics
# ===========
[docs]def two_points(name="TwoPoints"):
"""x**2 == 1"""
point1 = net.PointNet([-1.0])
point2 = net.PointNet([1.0])
return net.NetCollection([point1, point2], name=name)
# ======
# Conics
# ======
[docs]def circle(name="Ellipse"):
"""x**2 + y**2 == 1"""
f_circle = lambda u: np.array([np.sin(u), np.cos(u)])
circle = net.SmoothNet(f_circle, [[0, 2 * np.pi, True]], name=name)
return circle
[docs]def hyperbola(rot=False, name="Hyperbola"):
"""x**2 - y**2 == -1"""
if not rot:
f_hyperbola1 = lambda u: np.array([np.sinh(u), np.cosh(u)])
f_hyperbola2 = lambda u: np.array([np.sinh(u), -np.cosh(u)])
else:
f_hyperbola1 = lambda u: np.array([np.cosh(u), np.sinh(u)])
f_hyperbola2 = lambda u: np.array([-np.cosh(u), np.sinh(u)])
hyperbola1 = net.SmoothNet(f_hyperbola1, [[-np.inf, np.inf]])
hyperbola2 = net.SmoothNet(f_hyperbola2, [[-np.inf, np.inf]])
return net.NetCollection([hyperbola1, hyperbola2], name=name)
[docs]def linepair(name="ParallelLines"):
"""x**2 == 1"""
f_line1 = lambda u: np.array([1, u])
f_line2 = lambda u: np.array([-1, u])
line1 = net.SmoothNet(f_line1, [[-np.inf, np.inf]])
line2 = net.SmoothNet(f_line2, [[-np.inf, np.inf]])
return net.NetCollection([line1, line2], name=name)
[docs]def interline(name="IntersectingLines"):
"""x**2 - y**2 == 0"""
f_line1 = lambda u: np.array([u, u])
f_line2 = lambda u: np.array([u, -u])
line1 = net.SmoothNet(f_line1, [[-np.inf, np.inf]])
line2 = net.SmoothNet(f_line2, [[-np.inf, np.inf]])
return net.NetCollection([line1, line2], name=name)
[docs]def cline(rot=False, name="Line"):
"""x**2 == 0"""
if not rot:
f_line = lambda u: np.array([0, u])
else:
f_line = lambda u: np.array([u, 0])
line = net.SmoothNet(f_line, [[-np.inf, np.inf]], name=name)
return line
[docs]def parabola(name="Parabola"):
"""x**2 + 2*y == 0"""
f_parabola = lambda u: np.array([u, -(u**2) / 2])
parabola = net.SmoothNet(f_parabola, [[-np.inf, np.inf]], name=name)
return parabola
# ========
# Quadrics
# ========
# -- non-degenerate --
[docs]def sphere(axis=2, name="Ellipsoid"):
"""x**2 + y**2 + z**2 == 1"""
f_circle = lambda u: np.array([np.sin(u), np.cos(u)])
half_circle = net.SmoothNet(f_circle, [[0, np.pi]])
return nutils.surface_of_revolution(half_circle, axis=axis, name=name)
[docs]def hyperboloid_one_sheeted(axis=2, name="HyperboloidOneSheeted"):
"""x**2 + y**2 - z**2 == 1"""
f_hyperbola = lambda u: np.array([np.cosh(u), np.sinh(u)])
hyperbola = net.SmoothNet(f_hyperbola, [[-np.inf, np.inf]])
return nutils.surface_of_revolution(hyperbola, axis=axis, name=name)
[docs]def hyperboloid_two_sheeted(axis=2, name="HyperboloidTwoSheeted"):
"""x**2 + y**2 - z**2 == -1"""
f_hyperbola1 = lambda u: np.array([np.sinh(u), np.cosh(u)])
f_hyperbola2 = lambda u: np.array([np.sinh(u), -np.cosh(u)])
hyperbola1 = net.SmoothNet(f_hyperbola1, [[0, np.inf]])
hyperbola2 = net.SmoothNet(f_hyperbola2, [[0, np.inf]])
hyperbola1 = nutils.surface_of_revolution(hyperbola1, axis=axis, name=name)
hyperbola2 = nutils.surface_of_revolution(hyperbola2, axis=axis, name=name)
return net.NetCollection([hyperbola1, hyperbola2], name=name)
[docs]def paraboloid_elliptic(axis=2, name="ParaboloidElliptic"):
"""x**2 + y**2 + 2*z == 0"""
f_parabola = lambda u: np.array([u, -(u**2) / 2])
parabola = net.SmoothNet(f_parabola, [[0, np.inf]])
return nutils.surface_of_revolution(parabola, axis=axis, name=name)
[docs]def paraboloid_hyperbolic(axis=2, name="ParaboloidHyperbolic"):
"""x**2 - y**2 + 2*z == 0"""
i, j, k = 0, 1, 2
if axis == 0:
i, j, k = 1, 2, 0
elif axis == 1:
i, j, k = 0, 2, 1
def f(s, t):
x = np.zeros(3)
x[i] = (s + t) / 2
x[j] = (s - t) / 2
x[k] = -s * t / 2
return x
return net.SmoothNet(f, ([-np.inf, np.inf], [-np.inf, np.inf]), name=name)
# -- degenerate --
[docs]def cylinder_elliptic(axis=2, name="CylinderElliptic"):
"""x**2 + y**2 == 1"""
f_circle = lambda u: np.array([np.cos(u), np.sin(u)])
circle = net.SmoothNet(f_circle, [[0.0, 2 * np.pi, True]])
return nutils.cylinder(circle, axis=axis, name=name)
[docs]def cylinder_hyperbolic(axis=2, rot=False, name="CylinderHyperbolic"):
"""x**2 - y**2 == 1"""
if not rot:
f_hyperbola1 = lambda u: np.array([np.cosh(u), np.sinh(u)])
f_hyperbola2 = lambda u: np.array([-np.cosh(u), np.sinh(u)])
else:
f_hyperbola1 = lambda u: np.array([np.sinh(u), np.cosh(u)])
f_hyperbola2 = lambda u: np.array([np.sinh(u), -np.cosh(u)])
hyperbola1 = net.SmoothNet(f_hyperbola1, [[-np.inf, np.inf]])
hyperbola2 = net.SmoothNet(f_hyperbola2, [[-np.inf, np.inf]])
hyperbola1 = nutils.cylinder(hyperbola1, axis=axis, name=name)
hyperbola2 = nutils.cylinder(hyperbola2, axis=axis, name=name)
return net.NetCollection([hyperbola1, hyperbola2], name=name)
[docs]def cylinder_parabolic(axis=1, name="CylinderParabolic"):
"""x**2 + 2*z == 0"""
f_parabola = lambda u: np.array([u, -(u**2) / 2])
parabola = net.SmoothNet(f_parabola, [[-np.inf, np.inf]])
return nutils.cylinder(parabola, axis=axis, name=name)
[docs]def cone(axis=2, name="Cone"):
"""x**2 + y**2 - z**2 == 0"""
def f_circle(u):
x = np.array([np.cos(u), np.sin(u), 1.0])
x = np.insert(x, axis, 1.0)
x = x[:-1]
return x
circle = net.SmoothNet(f_circle, [[0.0, 2 * np.pi, True]])
return nutils.cone(circle, [0.0, 0.0, 0.0], name=name)
[docs]def line(axis=2, name="Line"):
"""x**2 + y**2 == 0"""
def f_line(u):
x = np.zeros(3)
x[axis] = u
return x
return net.SmoothNet(f_line, [[-np.inf, np.inf]], name=name)
[docs]def planes_intersecting(axis=2, name="Planes"):
"""x**2 - y**2 == 0"""
def f_plane1(u, v):
x = np.array([u, u, 1.0])
x = np.insert(x, axis, v)
x = x[:-1]
return x
def f_plane2(u, v):
x = np.array([u, -u, 1.0])
x = np.insert(x, axis, v)
x = x[:-1]
return x
s1 = net.SmoothNet(f_plane1, ([-np.inf, np.inf], [-np.inf, np.inf]))
s2 = net.SmoothNet(f_plane2, ([-np.inf, np.inf], [-np.inf, np.inf]))
return net.NetCollection([s1, s2], name=name)
[docs]def planes_parallel(axis=0, name="PlanesParallel"):
"""x**2 = 1"""
def f_plane1(u, v):
x = np.array([u, v, 1.0])
x = np.insert(x, axis, 1)
x = x[:-1]
return x
def f_plane2(u, v):
x = np.array([u, v, 1.0])
x = np.insert(x, axis, -1)
x = x[:-1]
return x
s1 = net.SmoothNet(f_plane1, ([-np.inf, np.inf], [-np.inf, np.inf]))
s2 = net.SmoothNet(f_plane2, ([-np.inf, np.inf], [-np.inf, np.inf]))
return net.NetCollection([s1, s2], name=name)
[docs]def plane(axis=0, name="Plane"):
"""double x**2 = 0, or single 2*x == 0"""
def f_plane(u, v):
x = np.array([u, v, 1.0])
x = np.insert(x, axis, 0.0)
x = x[:-1]
return x
s = net.SmoothNet(f_plane, ([-np.inf, np.inf], [-np.inf, np.inf]), name=name)
return s