Shapes

sies.shapes.boundary

Discretized \(C^2\)-smooth closed boundaries.

A C2Boundary stores the sampled boundary points of a simply connected planar domain together with the tangent, acceleration and outward normal vectors. It supports rigid motions and scaling, which return new objects (instances are immutable by convention).

C2Boundary

C2Boundary(points, tvec, avec, normal, center_of_mass=None, name='')

A discretized \(C^2\)-smooth closed curve in the plane.

Parameters:

Name	Type	Description	Default
points	(ndarray, shape(2, n))	Coordinates of the boundary points. The curve must not be tied off: the first and last points are distinct.	required
tvec	(ndarray, shape(2, n))	Tangent vectors at the boundary points.	required
avec	(ndarray, shape(2, n))	Acceleration (second derivative) vectors.	required
normal	(ndarray, shape(2, n))	Outward unit normal vectors.	required
center_of_mass	(ndarray, shape(2))	Center of mass of the enclosed domain. Computed from the Stokes formula if not provided.	None
name	str	Human-readable name of the shape.	""

Attributes:

Name	Type	Description
points	(ndarray, shape(2, n))	Boundary points.
tvec	(ndarray, shape(2, n))	Tangent vectors.
avec	(ndarray, shape(2, n))	Acceleration vectors.
normal	(ndarray, shape(2, n))	Outward unit normal vectors.
name	str	Name of the shape.

Source code in src/sies/shapes/boundary.py

def __init__(
    self,
    points: NDArray,
    tvec: NDArray,
    avec: NDArray,
    normal: NDArray,
    center_of_mass: NDArray | None = None,
    name: str = "",
):
    for arr in (points, tvec, avec, normal):
        if arr.ndim != 2 or arr.shape[0] != 2:
            raise ValueError("Boundary arrays must have shape (2, n).")
    if not points.shape == tvec.shape == avec.shape == normal.shape:
        raise ValueError("All boundary arrays must have the same shape.")

    self.points = points
    self.tvec = tvec
    self.avec = avec
    self.normal = normal
    self.name = name

    if center_of_mass is None:
        center_of_mass = self.stokes_center_of_mass(points, tvec, normal)
    self.center_of_mass = np.asarray(center_of_mass, dtype=float).reshape(2)

nb_points property

nb_points

int: Number of discretization points.

theta property

theta

ndarray: Uniform parameterization in [0, 2 pi), not tied off.

cpoints property

cpoints

ndarray: Boundary points as complex numbers x + i y.

tvec_norm property

tvec_norm

ndarray: Euclidean norm of the tangent vectors.

sigma property

sigma

ndarray: Curve integration elements.

The boundary integral of a function f is approximated by sum(f(points) * sigma).

perimeter property

perimeter

float: Perimeter of the boundary.

area property

area

float: Area of the enclosed domain, from the Stokes formula.

diameter property

diameter

float: Upper bound of the shape diameter, from the center of mass.

box property

box

Tuple of float: Width and height of the minimal bounding box.

principal_direction property

principal_direction

ndarray: Principal direction of the shape (unit vector).

The sign is normalized so that the first nonzero component is positive (the direction is defined up to a sign).

__add__

__add__(z0)

Translate the boundary by the vector z0.

Source code in src/sies/shapes/boundary.py

def __add__(self, z0: NDArray) -> "C2Boundary":
    """Translate the boundary by the vector `z0`."""
    z0 = np.asarray(z0, dtype=float).reshape(2)
    out = copy.copy(self)
    out.points = self.points + z0[:, np.newaxis]
    out.center_of_mass = self.center_of_mass + z0
    return out

__sub__

__sub__(z0)

Translate the boundary by the vector -z0.

Source code in src/sies/shapes/boundary.py

def __sub__(self, z0: NDArray) -> "C2Boundary":
    """Translate the boundary by the vector ``-z0``."""
    return self + (-np.asarray(z0, dtype=float))

__mul__

__mul__(s)

Scale the boundary by the positive factor s.

Source code in src/sies/shapes/boundary.py

def __mul__(self, s: float) -> "C2Boundary":
    """Scale the boundary by the positive factor `s`."""
    if s <= 0:
        raise ValueError("Scaling factor must be positive.")
    out = copy.copy(self)
    out.points = self.points * s
    out.tvec = self.tvec * s
    out.avec = self.avec * s
    out.center_of_mass = self.center_of_mass * s
    for attr in self._scale_attrs:
        setattr(out, attr, getattr(self, attr) * s)
    return out

rotate

rotate(phi)

Rotate the boundary around the origin.

Parameters:

Name	Type	Description	Default
phi	float	Rotation angle in radians (counterclockwise).	required

Returns:

Type	Description
C2Boundary	The rotated boundary.

Source code in src/sies/shapes/boundary.py

def rotate(self, phi: float) -> "C2Boundary":
    """Rotate the boundary around the origin.

    Parameters
    ----------
    phi : float
        Rotation angle in radians (counterclockwise).

    Returns
    -------
    C2Boundary
        The rotated boundary.
    """
    rot = np.array([[np.cos(phi), -np.sin(phi)], [np.sin(phi), np.cos(phi)]])
    out = copy.copy(self)
    out.points = rot @ self.points
    out.tvec = rot @ self.tvec
    out.avec = rot @ self.avec
    out.normal = rot @ self.normal
    out.center_of_mass = rot @ self.center_of_mass
    return out

is_inside

is_inside(x)

Check whether a point lies in the ball circumscribing the shape.

Parameters:

Name	Type	Description	Default
x	(ndarray, shape(2))	Query point.	required

Returns:

Type	Description
bool	True if x is inside the ball centered at the center of mass with radius diameter / 2.

Source code in src/sies/shapes/boundary.py

def is_inside(self, x: NDArray) -> bool:
    """Check whether a point lies in the ball circumscribing the shape.

    Parameters
    ----------
    x : ndarray, shape (2,)
        Query point.

    Returns
    -------
    bool
        True if `x` is inside the ball centered at the center of
        mass with radius ``diameter / 2``.
    """
    return bool(
        np.linalg.norm(np.asarray(x).reshape(2) - self.center_of_mass) < self.diameter / 2
    )

is_disjoint

is_disjoint(other)

Check whether two shapes have disjoint circumscribed balls.

Parameters:

Name	Type	Description	Default
other	C2Boundary	The other shape.	required

Returns:

Type	Description
bool	True if the balls B(com, diameter / 2) of both shapes do not intersect.

Source code in src/sies/shapes/boundary.py

def is_disjoint(self, other: "C2Boundary") -> bool:
    """Check whether two shapes have disjoint circumscribed balls.

    Parameters
    ----------
    other : C2Boundary
        The other shape.

    Returns
    -------
    bool
        True if the balls ``B(com, diameter / 2)`` of both shapes do
        not intersect.
    """
    dist = np.linalg.norm(self.center_of_mass - other.center_of_mass)
    return bool(dist > (self.diameter + other.diameter) / 2)

plot

plot(ax=None, **kwargs)

Plot the boundary curve.

Parameters:

Name	Type	Description	Default
ax	Axes	Axes to draw on. A new figure is created if omitted.	None
**kwargs		Forwarded to plot.	{}

Returns:

Type	Description
Axes	The axes containing the plot.

Source code in src/sies/shapes/boundary.py

def plot(self, ax=None, **kwargs):
    """Plot the boundary curve.

    Parameters
    ----------
    ax : matplotlib.axes.Axes, optional
        Axes to draw on. A new figure is created if omitted.
    **kwargs
        Forwarded to `plot`.

    Returns
    -------
    matplotlib.axes.Axes
        The axes containing the plot.
    """
    import matplotlib.pyplot as plt

    if ax is None:
        _, ax = plt.subplots()
    closed = np.column_stack([self.points, self.points[:, :1]])
    ax.plot(closed[0], closed[1], **kwargs)
    ax.set_aspect("equal")
    return ax

stokes_center_of_mass staticmethod

stokes_center_of_mass(points, tvec, normal)

Compute the center of mass of a domain by the Stokes formula.

Parameters:

Name	Type	Description	Default
points	(ndarray, shape(2, n))	Boundary points.	required
tvec	(ndarray, shape(2, n))	Tangent vectors.	required
normal	(ndarray, shape(2, n))	Outward unit normal vectors.	required

Returns:

Type	Description
(ndarray, shape(2))	Center of mass of the enclosed domain.

Source code in src/sies/shapes/boundary.py

@staticmethod
def stokes_center_of_mass(points: NDArray, tvec: NDArray, normal: NDArray) -> NDArray:
    """Compute the center of mass of a domain by the Stokes formula.

    Parameters
    ----------
    points : ndarray, shape (2, n)
        Boundary points.
    tvec : ndarray, shape (2, n)
        Tangent vectors.
    normal : ndarray, shape (2, n)
        Outward unit normal vectors.

    Returns
    -------
    ndarray, shape (2,)
        Center of mass of the enclosed domain.
    """
    nb_points = points.shape[1]
    sigma = 2 * np.pi / nb_points * np.linalg.norm(tvec, axis=0)
    mass = np.sum(points * normal * sigma) / 2
    first_moment = np.sum(points**2 * normal * sigma, axis=1) / 2
    return first_moment / mass

check_sampling staticmethod

check_sampling(points)

Check that the sampled curve is locally consistent.

By a Taylor expansion, a \(C^1\) simple curve sampled finely enough satisfies \(\langle f(t_{n+1}) - f(t_n), f(t_n) - f(t_{n-1}) \rangle > 0\).

Parameters:

Name	Type	Description	Default
points	(ndarray, shape(2, n))	Boundary points.	required

Returns:

Type	Description
bool	True if the condition holds at every point.

Source code in src/sies/shapes/boundary.py

@staticmethod
def check_sampling(points: NDArray) -> bool:
    r"""Check that the sampled curve is locally consistent.

    By a Taylor expansion, a $C^1$ simple curve sampled finely
    enough satisfies
    $\langle f(t_{n+1}) - f(t_n), f(t_n) - f(t_{n-1}) \rangle > 0$.

    Parameters
    ----------
    points : ndarray, shape (2, n)
        Boundary points.

    Returns
    -------
    bool
        True if the condition holds at every point.
    """
    prev = np.roll(points, 1, axis=1)
    nxt = np.roll(points, -1, axis=1)
    forward = nxt - points
    backward = points - prev
    inner = np.sum(forward * backward, axis=0)
    return bool(np.all(inner > 0))

sies.shapes.standard

Catalog of standard parametric shapes.

Every class derives from C2Boundary and provides the analytic tangent, acceleration and normal vectors of its boundary whenever possible. Shapes with corners (triangle, rectangle) are smoothed by spline resampling.

Ellipse

Ellipse(axis_a, axis_b, nb_points)

Bases: C2Boundary

An ellipse centered at the origin.

Parameters:

Name	Type	Description	Default
axis_a	float	Length of the semi-major axis.	required
axis_b	float	Length of the semi-minor axis; must satisfy axis_b <= axis_a.	required
nb_points	int	Number of boundary discretization points.	required

Attributes:

Name	Type	Description
axis_a, axis_b	float	Semi-axis lengths.
phi	float	Cumulated rotation angle of the shape.

Source code in src/sies/shapes/standard.py

def __init__(self, axis_a: float, axis_b: float, nb_points: int):
    if axis_a < axis_b:
        raise ValueError("The semi-major axis must be longer than the semi-minor one.")

    theta = 2 * np.pi * np.arange(nb_points) / nb_points
    points = np.vstack([axis_a * np.cos(theta), axis_b * np.sin(theta)])
    tvec = np.vstack([-axis_a * np.sin(theta), axis_b * np.cos(theta)])
    avec = -points
    name = "Circle" if axis_a == axis_b else "Ellipse"
    super().__init__(points, tvec, avec, _outward_normal(tvec), np.zeros(2), name)

    self.axis_a = axis_a
    self.axis_b = axis_b
    self.phi = 0.0

rotate

rotate(phi)

Rotate the ellipse around the origin by the angle phi.

Source code in src/sies/shapes/standard.py

def rotate(self, phi: float) -> "Ellipse":
    """Rotate the ellipse around the origin by the angle `phi`."""
    out = super().rotate(phi)
    out.phi = self.phi + phi
    return out

Flower

Flower(axis_a, axis_b, nb_points, nb_petals=5, epsilon=0.3, pertb=1)

Bases: C2Boundary

A flower-like perturbation of an ellipse.

The boundary is the curve \(t \mapsto R_\phi A \, (1 + \epsilon \cos^k(nt)) (\cos t, \sin t)\) where A = diag(a, b) and n is the number of petals.

Parameters:

Name	Type	Description	Default
axis_a	float	Length of the semi-major axis of the underlying ellipse.	required
axis_b	float	Length of the semi-minor axis.	required
nb_points	int	Number of boundary discretization points.	required
nb_petals	int	Number of petals.	5
epsilon	float	Strength of the petal perturbation.	0.3
pertb	int	Exponent k of the perturbation; must be a positive integer.	1

Attributes:

Name	Type	Description
nb_petals	int	Number of petals.
epsilon	float	Perturbation strength.

Source code in src/sies/shapes/standard.py

def __init__(
    self,
    axis_a: float,
    axis_b: float,
    nb_points: int,
    nb_petals: int = 5,
    epsilon: float = 0.3,
    pertb: int = 1,
):
    if pertb < 1:
        raise ValueError("The perturbation exponent must be a positive integer.")

    n, k, eps = nb_petals, pertb, epsilon
    theta = 2 * np.pi * np.arange(nb_points) / nb_points
    scale = np.diag([axis_a, axis_b])

    radius = 1 + eps * np.cos(n * theta) ** k
    points = scale @ np.vstack([np.cos(theta) * radius, np.sin(theta) * radius])

    dradius = -k * eps * n * np.sin(n * theta) * np.cos(n * theta) ** (k - 1)
    tvec = scale @ np.vstack(
        [
            -np.sin(theta) * radius + dradius * np.cos(theta),
            np.cos(theta) * radius + dradius * np.sin(theta),
        ]
    )

    ddradius = (
        -k
        * eps
        * n**2
        * (
            np.cos(n * theta) ** k
            - (k - 1) * np.sin(n * theta) ** 2 * np.cos(n * theta) ** (k - 2)
        )
    )
    avec = scale @ np.vstack(
        [
            -np.cos(theta) * radius + 2 * dradius * -np.sin(theta) + ddradius * np.cos(theta),
            -np.sin(theta) * radius + 2 * dradius * np.cos(theta) + ddradius * np.sin(theta),
        ]
    )

    super().__init__(points, tvec, avec, _outward_normal(tvec), np.zeros(2), "Flower")

    self.nb_petals = nb_petals
    self.epsilon = epsilon
    self.axis_a = axis_a
    self.axis_b = axis_b

Triangle

Triangle(side, angle, nb_points, downsample=10)

Bases: C2Boundary

An isosceles triangle with smoothed corners.

Parameters:

Name	Type	Description	Default
side	float	Length of the two equal sides.	required
angle	float	Angle between the two equal sides, in radians.	required
nb_points	int	Number of boundary discretization points.	required
downsample	int	Down-sampling factor used to smooth out the corners; values larger than one produce a \(C^2\) boundary.	10

Attributes:

Name	Type	Description
side	float	Length of the equal sides.
angle	float	Apex angle.

Source code in src/sies/shapes/standard.py

def __init__(self, side: float, angle: float, nb_points: int, downsample: int = 10):
    height = side * np.cos(angle / 2)
    half_base = side * np.sin(angle / 2)

    vertices = np.array(
        [
            [0, 2 / 3 * height],
            [-half_base, -height / 3],
            [half_base, -height / 3],
        ]
    ).T
    fractions = np.array([side, 2 * half_base, side]) / (2 * (side + half_base))
    points0 = _polygon_points(vertices, fractions, nb_points)

    theta = 2 * np.pi * np.arange(nb_points) / nb_points
    points, tvec, avec, normal = resample_curve(
        points0, theta, nb_points, downsample=downsample
    )
    super().__init__(points, tvec, avec, normal, np.zeros(2), "Triangle")

    self.side = side
    self.angle = angle

Rectangle

Rectangle(width, height, nb_points, downsample=10)

Bases: C2Boundary

An axis-aligned rectangle with smoothed corners.

Parameters:

Name	Type	Description	Default
width	float	Size along the horizontal axis.	required
height	float	Size along the vertical axis.	required
nb_points	int	Number of boundary discretization points.	required
downsample	int	Down-sampling factor used to smooth out the corners.	10

Attributes:

Name	Type	Description
width, height	float	Rectangle dimensions.

Source code in src/sies/shapes/standard.py

def __init__(self, width: float, height: float, nb_points: int, downsample: int = 10):
    a, b = height, width
    vertices = np.array([[-b, a], [-b, -a], [b, -a], [b, a]]).T / 2
    fractions = np.array([a, b, a, b]) / (2 * (a + b))
    points0 = _polygon_points(vertices, fractions, nb_points)

    theta = 2 * np.pi * np.arange(nb_points) / nb_points
    points, tvec, avec, normal = resample_curve(
        points0, theta, nb_points, downsample=downsample
    )
    name = "Square" if width == height else "Rectangle"
    super().__init__(points, tvec, avec, normal, np.zeros(2), name)

    self.width = width
    self.height = height

Banana

Banana(axis_a, axis_b, center, curvature_center, nb_points)

Bases: C2Boundary

A banana-shaped object (an ellipse bent along a circular arc).

Parameters:

Name	Type	Description	Default
axis_a	float	Length of the semi-major axis of the underlying ellipse.	required
axis_b	float	Length of the semi-minor axis.	required
center	(ndarray, shape(2))	Center of the underlying ellipse.	required
curvature_center	(ndarray, shape(2))	Center of the bending circle.	required
nb_points	int	Number of boundary discretization points.	required

Source code in src/sies/shapes/standard.py

def __init__(
    self,
    axis_a: float,
    axis_b: float,
    center: NDArray,
    curvature_center: NDArray,
    nb_points: int,
):
    if axis_a < axis_b:
        raise ValueError("The semi-major axis must be longer than the semi-minor one.")

    x0, y0 = np.asarray(center, dtype=float)
    xc, yc = np.asarray(curvature_center, dtype=float)

    radius = np.hypot(xc - x0, yc - y0)
    theta0 = np.arctan2(y0 - yc, x0 - xc)
    alpha = axis_a / radius
    b = axis_b

    theta = 2 * np.pi * np.arange(nb_points) / nb_points
    t = theta0 + alpha * np.cos(theta)
    rho = radius + b * np.sin(theta)
    points = np.vstack([xc + rho * np.cos(t), yc + rho * np.sin(t)])

    tvec = np.vstack(
        [
            b * np.cos(theta) * np.cos(t) + alpha * np.sin(theta) * rho * np.sin(t),
            b * np.cos(theta) * np.sin(t) - alpha * np.sin(theta) * rho * np.cos(t),
        ]
    )

    avec = np.vstack(
        [
            -b * np.sin(theta) * np.cos(t)
            + alpha * b * np.cos(theta) * np.sin(t)
            + alpha * np.cos(theta) * rho * np.sin(t)
            + alpha * b * np.cos(theta) ** 2 * np.sin(t)
            - alpha**2 * np.sin(theta) ** 2 * rho * np.cos(t),
            -b * np.sin(theta) * np.cos(t)
            - alpha * b * np.sin(theta) * np.cos(t)
            - alpha * np.cos(theta) * rho * np.cos(t)
            - alpha * b * np.cos(theta) ** 2 * np.cos(t)
            - alpha**2 * np.sin(theta) ** 2 * rho * np.sin(t),
        ]
    )

    # The banana boundary is parameterized clockwise, so the outward
    # normal uses the rotation opposite to `_outward_normal` (which
    # assumes counterclockwise curves).
    normal = np.vstack([-tvec[1], tvec[0]])
    normal = normal / np.linalg.norm(normal, axis=0)

    super().__init__(points, tvec, avec, normal, None, "Banana")

    self.axis_a = axis_a
    self.axis_b = axis_b

sies.shapes.resampling

Resampling of closed curves by periodic cubic splines.

These utilities reconstruct the tangent, acceleration and outward normal vectors of a closed curve from sampled points. Corner singularities can be smoothed out by down-sampling the curve before the spline fit.

resample_curve

resample_curve(points, theta, nb_points, box=None, downsample=1)

Resample a closed curve and compute its differential quantities.

The curve is optionally rescaled to fit a bounding box, down-sampled (which smooths out corner singularities), then re-interpolated on a uniform parameter grid with a periodic cubic spline.

Parameters:

Name	Type	Description	Default
points	(ndarray, shape(2, n))	Coordinates of the curve samples, not tied off.	required
theta	(ndarray, shape(n))	Parameterization of the samples in [0, 2 pi).	required
nb_points	int	Number of points of the resampled curve.	required
box	tuple of float	If given, target size (width, height) of the bounding box of the curve.	None
downsample	int	Down-sampling factor applied before interpolation; values larger than one smooth out corners.	1

Returns:

Name	Type	Description
points	(ndarray, shape(2, nb_points))	Resampled curve.
tvec	(ndarray, shape(2, nb_points))	Tangent vectors (first derivative).
avec	(ndarray, shape(2, nb_points))	Acceleration vectors (second derivative).
normal	(ndarray, shape(2, nb_points))	Outward unit normal vectors.

Source code in src/sies/shapes/resampling.py

def resample_curve(
    points: NDArray,
    theta: NDArray,
    nb_points: int,
    box: tuple[float, float] | None = None,
    downsample: int = 1,
) -> tuple[NDArray, NDArray, NDArray, NDArray]:
    """Resample a closed curve and compute its differential quantities.

    The curve is optionally rescaled to fit a bounding box, down-sampled
    (which smooths out corner singularities), then re-interpolated on a
    uniform parameter grid with a periodic cubic spline.

    Parameters
    ----------
    points : ndarray, shape (2, n)
        Coordinates of the curve samples, not tied off.
    theta : ndarray, shape (n,)
        Parameterization of the samples in ``[0, 2 pi)``.
    nb_points : int
        Number of points of the resampled curve.
    box : tuple of float, optional
        If given, target size ``(width, height)`` of the bounding box of
        the curve.
    downsample : int, default 1
        Down-sampling factor applied before interpolation; values larger
        than one smooth out corners.

    Returns
    -------
    points : ndarray, shape (2, nb_points)
        Resampled curve.
    tvec : ndarray, shape (2, nb_points)
        Tangent vectors (first derivative).
    avec : ndarray, shape (2, nb_points)
        Acceleration vectors (second derivative).
    normal : ndarray, shape (2, nb_points)
        Outward unit normal vectors.
    """
    downsample = int(np.ceil(downsample))
    if downsample < 1:
        raise ValueError("Down-sampling factor must be positive.")

    if box is not None:
        xmin, xmax = points[0].min(), points[0].max()
        ymin, ymax = points[1].min(), points[1].max()
        if np.isclose(xmax, xmin) or np.isclose(ymax, ymin):
            raise ValueError("Curve must have nonzero extent in both axes for box rescaling.")
        center = np.array([(xmin + xmax) / 2, (ymin + ymax) / 2])
        scale = np.array([box[0] / (xmax - xmin), box[1] / (ymax - ymin)])
        points = (points - center[:, np.newaxis]) * scale[:, np.newaxis]

    spline = _periodic_spline(theta[::downsample], points[:, ::downsample])

    theta_new = 2 * np.pi * np.arange(nb_points) / nb_points
    new_points = spline(theta_new)
    tvec = spline(theta_new, 1)
    avec = spline(theta_new, 2)

    normal = np.vstack([tvec[1], -tvec[0]])
    normal = normal / np.linalg.norm(normal, axis=0)
    return new_points, tvec, avec, normal