"""Kernel functions for boundary integral operators."""
from __future__ import annotations
import numpy as np
from numpy.typing import NDArray
TWOPI = 2.0 * np.pi
[docs]
def laplace_single_layer_kernel(
x: NDArray[np.floating],
y: NDArray[np.floating],
) -> float:
"""Single-layer kernel G(x,y) = (1/2π) log|x-y| for 2D Laplace.
Parameters
----------
x : array, shape (2,)
Target point
y : array, shape (2,)
Source point
Returns
-------
float
Kernel value (undefined for x == y)
"""
r = np.linalg.norm(x - y)
if r < 1e-15:
return 0.0 # Caller must handle diagonal separately
return np.log(r) / TWOPI
[docs]
def laplace_double_layer_kernel(
x: NDArray[np.floating],
y: NDArray[np.floating],
n_y: NDArray[np.floating],
) -> float:
"""Double-layer kernel ∂G/∂n_y for 2D Laplace.
∂G/∂n_y = (1/2π) * (x - y) · n_y / |x - y|²
Parameters
----------
x : array, shape (2,)
Target point
y : array, shape (2,)
Source point
n_y : array, shape (2,)
Outward unit normal at y
Returns
-------
float
Kernel value (undefined for x == y)
"""
diff = x - y
r2 = np.dot(diff, diff)
if r2 < 1e-30:
return 0.0 # Caller must handle diagonal separately
return np.dot(diff, n_y) / (TWOPI * r2)
[docs]
def laplace_double_layer_diagonal(kappa: float) -> float:
"""Diagonal limit of double-layer kernel for smooth curve.
lim_{y→x} ∂G(x,y)/∂n_y = κ(x) / 4π
where κ is the signed curvature.
"""
return kappa / (4.0 * np.pi)
# =============================================================================
# Vectorized versions for efficiency
# =============================================================================
[docs]
def laplace_single_layer_matrix(
targets: NDArray[np.floating],
sources: NDArray[np.floating],
weights: NDArray[np.floating],
) -> NDArray[np.floating]:
"""Assemble single-layer operator matrix.
A[i,j] = G(x_i, y_j) * w_j
Parameters
----------
targets : array, shape (M, 2)
Target points
sources : array, shape (N, 2)
Source points
weights : array, shape (N,)
Quadrature weights
Returns
-------
A : array, shape (M, N)
Operator matrix
"""
M = targets.shape[0]
N = sources.shape[0]
# Compute all pairwise distances
diff = targets[:, np.newaxis, :] - sources[np.newaxis, :, :] # (M, N, 2)
r2 = np.sum(diff ** 2, axis=2) # (M, N)
r = np.sqrt(r2)
# Avoid log(0)
r = np.maximum(r, 1e-15)
A = np.log(r) / TWOPI * weights[np.newaxis, :]
return A
[docs]
def laplace_double_layer_matrix(
targets: NDArray[np.floating],
sources: NDArray[np.floating],
normals: NDArray[np.floating],
weights: NDArray[np.floating],
curvature: NDArray[np.floating] | None = None,
same_points: bool = False,
) -> NDArray[np.floating]:
"""Assemble double-layer operator matrix.
A[i,j] = ∂G(x_i, y_j)/∂n_j * w_j (off-diagonal)
A[i,i] = κ_i / 4π * w_i (diagonal, if same_points)
Parameters
----------
targets : array, shape (M, 2)
Target points
sources : array, shape (N, 2)
Source points (boundary nodes)
normals : array, shape (N, 2)
Outward unit normals at sources
weights : array, shape (N,)
Quadrature weights
curvature : array, shape (N,), optional
Signed curvature at sources (required if same_points=True)
same_points : bool
If True, targets == sources and diagonal correction is applied
Returns
-------
A : array, shape (M, N)
Operator matrix
"""
M = targets.shape[0]
N = sources.shape[0]
# Compute all pairwise differences
diff = targets[:, np.newaxis, :] - sources[np.newaxis, :, :] # (M, N, 2)
r2 = np.sum(diff ** 2, axis=2) # (M, N)
# Avoid division by zero
r2_safe = np.maximum(r2, 1e-30)
# dot(diff, normal) for each pair
dot_n = np.sum(diff * normals[np.newaxis, :, :], axis=2) # (M, N)
A = dot_n / (TWOPI * r2_safe) * weights[np.newaxis, :]
# Apply diagonal correction if same_points
if same_points:
if curvature is None:
raise ValueError("curvature required when same_points=True")
if M != N:
raise ValueError("M must equal N when same_points=True")
np.fill_diagonal(A, curvature / (4.0 * np.pi) * weights)
return A