volumential.singular_integral_2d#

The 2D singular integrals are computed using the transform described in http://link.springer.com/10.1007/BF00370482.

This module owns the host-side (numpy/scipy/mpmath) quadrature rules used to build near-field interaction tables:

volumential.singular_integral_2d.adaptive_quadrature(func, a, b, args=(), tol=1.49e-08, rtol=1.49e-08, maxiter=50, vec_func=False, miniter=1)[source]#

Approximate the removed scipy.integrate.quadrature API.

The legacy code relied on quadrature order refinement semantics from older SciPy. Recreate the essential behavior by increasing the fixed Gauss order until successive iterates stabilize.

volumential.singular_integral_2d.quad(func, a, b, args=(), tol=1.49e-08, rtol=1.49e-08, maxiter=50, vec_func=False, miniter=1)#

Approximate the removed scipy.integrate.quadrature API.

The legacy code relied on quadrature order refinement semantics from older SciPy. Recreate the essential behavior by increasing the fixed Gauss order until successive iterates stabilize.

volumential.singular_integral_2d.update_qquad_leggauss_formula(deg1, deg2) → None[source]#

Refresh the module-level tensor-product Gauss-Legendre weights.

Warning

Only quad_weights is updated; the node arrays quad_points_x and quad_points_y are shadowed by locals here and therefore left untouched. None of the module’s quadrature routines read these globals – they build their own rules – so this helper is kept only for backwards compatibility.

volumential.singular_integral_2d.qquad(func, a, b, c, d, args=(), tol=1.49e-08, rtol=1.49e-08, maxitero=50, maxiteri=50, vec_func=False, minitero=1, miniteri=1, method='Adaptive')[source]#

Computes a (tensor product) double integral.

Integrate func on [a, b]X[c, d] using Gaussian quadrature with absolute tolerance tol.

Parameters:

  • func (collections.abc.Callable) – A double variable Python function or method to integrate.

  • a (float) – Lower-left corner of integration region.

  • b (float) – Lower-right corner of integration region.

  • c (float) – Upper-left corner of integration region.

  • d (float) – Upper-right corner of integration region.

  • args (tuple) – Extra arguments to pass to function.

  • tol (float) – rtol Iteration stops when error between last two iterates is less than tol OR the relative change is less than rtol.

  • rtol (float) – Iteration stops when error between last two iterates is less than tol OR the relative change is less than rtol.

  • maxitero (int) – Maximum order of outer Gaussian quadrature.

  • maxiteri (int) – Maximum order of inner Gaussian quadrature.

  • vec_func (bool) – True if func handles arrays as arguments (is a “vector” function). Default is True.

  • minitero (int) – Minimum order of outer Gaussian quadrature.

  • miniteri (int) – Minimum order of inner Gaussian quadrature.

Returns:

  • val: Gaussian quadrature approximation (within tolerance) to integral.

  • err: Difference between last two estimates of the integral.

Return type:

tuple[float, float]

volumential.singular_integral_2d.solve_affine_map_2d(source_tria, target_tria)[source]#

Computes the affine map and its inverse that maps the source_tria to target_tria.

Parameters:
Returns:

  • mapping: the forward map.

  • J: the Jacobian.

  • invmap: the inverse map.

  • invJ: the Jacobian of inverse map.

Return type:

tuple[collections.abc.Callable, float, collections.abc.Callable, float]

volumential.singular_integral_2d.tria2rect_map_2d()[source]#

Returns the mapping and its inverse that maps a template triangle to a template rectangle.

  • Template triangle [T]: (0,0)–(1,0)–(0,1)–(0,0)

  • Template rectangle [R]: (0,0)–(1,0)–(1,pi/2)–(0,pi/2)–(0,0)

Returns:

The mapping, its Jacobian, its inverse, and the Jacobian of its inverse. Note that the Jacobians are returned as lambdas since they are not constants.

Return type:

tuple[collections.abc.Callable, collections.abc.Callable, collections.abc.Callable, collections.abc.Callable]

volumential.singular_integral_2d.is_in_t(pt)[source]#

Checks if a point is in the template triangle T.

Parameters:

pt (tuple[float, float]) – The point to be checked.

Returns:

True if pt is in T.

Return type:

bool

volumential.singular_integral_2d.is_in_r(pt, a=0, b=1, c=0, d=1.5707963267948966)[source]#

Checks if a point is in the (template) rectangle R.

Parameters:

pt (tuple[float, float]) – The point to be checked.

Returns:

True if pt is in R.

Return type:

bool

volumential.singular_integral_2d.is_collinear(p0, p1, p2) → bool[source]#

Return whether the three 2D points are (numerically) collinear.

volumential.singular_integral_2d.is_positive_triangle(tria) → bool[source]#

Return whether the triangle’s vertices are counter-clockwise ordered.

volumential.singular_integral_2d.tria_quad(func, tria, args=(), tol=1.49e-08, rtol=1.49e-08, maxiter=50, vec_func=True, miniter=1)[source]#

Computes a double integral on a general triangular region.

Integrate func on tria by transforming the region into a rectangle and using Gaussian quadrature with absolute tolerance tol.

The integrand, func, is allowed to have singularity at most $O(r)$ at the first vertex of the triangle. It is okay if func does not evaluate at the singular point. This function handles that automatically.

Parameters:

  • func (collections.abc.Callable) – A double variable Python function or method to integrate.

  • tria (tuple[tuple[float, float], tuple[float, float], tuple[float, float]]) – The triangular region to do quadrature.

  • args (tuple) – Extra arguments to pass to function.

  • tol (float) – rtol Iteration stops when error between last two iterates is less than tol OR the relative change is less than rtol.

  • rtol (float) – Iteration stops when error between last two iterates is less than tol OR the relative change is less than rtol.

  • maxiter (int) – Maximum order of Gaussian quadrature.

  • vec_func (bool) – True if func handles arrays as arguments (is a “vector” function). Default is True.

  • miniter (int) – Minimum order of Gaussian quadrature.

Returns:

  • val: Gaussian quadrature approximation (within tolerance)

    to integral.

  • err: Difference between last two estimates of the integral.

Return type:

tuple[float, float]

volumential.singular_integral_2d.tria_quad_duffy_radial(func, tria, args=(), radial_rule='tanh-sinh', deg_theta=20, radial_quad_order=61, mp_dps=50)[source]#

Integrate func over a triangle with a Duffy-type radial rule.

The triangle is mapped to the template rectangle by solve_affine_map_2d() followed by tria2rect_map_2d(), so an integrable singularity at tria[0] is absorbed into the Jacobian.

Parameters:

radial_rule – one of "tanh-sinh" (mpmath, most accurate), "tanh-sinh-fast" (precomputed double-precision nodes) or "adaptive" (Gauss order refinement).

Returns:

(value, error_estimate); the error estimate is always zero.

volumential.singular_integral_2d.quadri_quad_duffy_radial(func, quadrilateral, singular_point, args=(), radial_rule='tanh-sinh', deg_theta=20, radial_quad_order=61, mp_dps=50)[source]#

Duffy-type quadrature over a quadrilateral split at singular_point.

Returns:

(value, error_estimate).

volumential.singular_integral_2d.box_quad_duffy_radial(func, a, b, c, d, singular_point, args=(), radial_rule='tanh-sinh', deg_theta=20, radial_quad_order=61, mp_dps=50)[source]#

Duffy-type quadrature over [a, b] x [c, d].

A singular_point outside the box is projected onto its boundary.

Returns:

(value, error_estimate).

volumential.singular_integral_2d.tria_quad_tanh_sinh_radial(func, tria, args=(), deg_theta=20, mp_dps=50)[source]#

tria_quad_duffy_radial() with the "tanh-sinh" radial rule.

volumential.singular_integral_2d.quadri_quad_tanh_sinh_radial(func, quadrilateral, singular_point, args=(), deg_theta=20, mp_dps=50)[source]#

quadri_quad_duffy_radial() with the "tanh-sinh" radial rule.

volumential.singular_integral_2d.box_quad_tanh_sinh_radial(func, a, b, c, d, singular_point, args=(), deg_theta=20, mp_dps=50)[source]#

box_quad_duffy_radial() with the "tanh-sinh" radial rule.

volumential.singular_integral_2d.box_quad_duffy_radial_nd(func, bounds, singular_point, args=(), radial_rule='tanh-sinh-fast', deg_regular=20, radial_quad_order=61, mp_dps=50)[source]#

Duffy-type quadrature over an axis-aligned box in any dimension.

The box is split into orthants around singular_point and each orthant is integrated once per axis ordering, so that the radial coordinate always absorbs the singularity.

Parameters:

  • bounds – per-axis (lower, upper) pairs; their number sets the dimension.

  • singular_point – singular point, projected into the box.

  • radial_rule – see tria_quad_duffy_radial(); additionally accepts "adaptive".

Returns:

(value, error_estimate); the error estimate is always zero.

volumential.singular_integral_2d.box_quad(func, a, b, c, d, singular_point, args=(), tol=1.49e-08, rtol=1.49e-08, maxiter=50, vec_func=True, miniter=1)[source]#

Compute a singular 2D integral over an axis-aligned box.

Integrate func on [a, b] x [c, d] using transformed Gaussian quadrature around singular_point.

Parameters:

  • func – callable integrand of two variables.

  • a – lower x bound.

  • b – upper x bound.

  • c – lower y bound.

  • d – upper y bound.

  • singular_point – singular point as (x, y).

  • args – extra positional arguments passed to func.

  • tol – absolute tolerance for adaptive quadrature.

  • rtol – relative tolerance for adaptive quadrature.

  • maxiter – maximum adaptive quadrature order.

  • vec_func – whether func accepts vectorized array inputs.

  • miniter – minimum adaptive quadrature order.

Returns:

(value, error_estimate).

volumential.singular_integral_2d.quadri_quad(func, quadrilateral, singular_point, args=(), tol=1.49e-08, rtol=1.49e-08, maxiter=50, vec_func=True, miniter=1)[source]#

Compute a singular 2D integral over a quadrilateral.

Parameters:

  • func – callable integrand of two variables.

  • quadrilateral – vertices ((x1, y1), ..., (x4, y4)).

  • singular_point – singular point as (x, y).

  • args – extra positional arguments passed to func.

  • tol – absolute tolerance for adaptive quadrature.

  • rtol – relative tolerance for adaptive quadrature.

  • maxiter – maximum adaptive quadrature order.

  • vec_func – whether func accepts vectorized array inputs.

  • miniter – minimum adaptive quadrature order.

Returns:

(value, error_estimate).