<filename>volumential/singular_integral_2d.py
from __future__ import absolute_import, division, print_function

__copyright__ = "Copyright (C) 2017 - 2018 <NAME>"

__license__ = """
Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in
all copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
THE SOFTWARE.
"""

import logging

import numpy as np
import scipy as sp
import scipy.integrate as sint
from scipy.integrate import quadrature as quad  # noqa:F401

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

.. autofunction:: box_quad
"""

logger = logging.getLogger(__name__)

# {{{ quadrature on rectangular box with no singularity

quad_points_x = np.array([])
quad_points_y = np.array([])
quad_weights = np.array([])


def update_qquad_leggauss_formula(deg1, deg2):

    x1, w1 = np.polynomial.legendre.leggauss(deg1)
    x1 = (x1 + 1) / 2
    w1 = w1 / 2

    x2, w2 = np.polynomial.legendre.leggauss(deg2)
    x2 = (x2 + 1) / 2
    w2 = w2 / 2

    quad_points_x, quad_points_y = np.meshgrid(x1, x2)
    ww1, ww2 = np.meshgrid(w1, w2)
    global quad_weights
    quad_weights = ww1 * ww2


def 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",
):
    """Computes a (tensor product) double integral.

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

    :param func: A double variable Python function or method to integrate.
    :type func: function.
    :param a: Lower-left corner of integration region.
    :type a: float.
    :param b: Lower-right corner of integration region.
    :type b: float.
    :param c: Upper-left corner of integration region.
    :type c: float.
    :param d: Upper-right corner of integration region.
    :type d: float.
    :param args: Extra arguments to pass to function.
    :type args: tuple, optional.
    :param tol: rtol Iteration stops when error between last two iterates is
                less than tol OR the relative change is less than rtol.
    :type tol: float, optional.
    :param rtol: Iteration stops when error between last two iterates is less
                than tol OR the relative change is less than rtol.
    :type rtol: float, optional.
    :param maxitero: Maximum order of outer Gaussian quadrature.
    :type maxitero: int, optional.
    :param maxiteri: Maximum order of inner Gaussian quadrature.
    :type maxiteri: int, optional.
    :param vec_func: True if func handles arrays as arguments (is a "vector"
                function). Default is True.
    :type vec_func: bool, optional.
    :param minitero: Minimum order of outer Gaussian quadrature.
    :type minitero: int, optional.
    :param miniteri: Minimum order of inner Gaussian quadrature.
    :type miniteri: int, optional.

    :returns:
        - **val**: Gaussian quadrature approximation (within tolerance) to integral.
        - **err**: Difference between last two estimates of the integral.
    :rtype: tuple(float,float).
    """

    l1 = b - a
    l2 = d - c
    assert l1 > 0
    assert l2 > 0
    toli = tol / l2
    rtoli = rtol / l2

    if method == "Adaptive":

        # Using lambda for readability
        def outer_integrand(y):
            return sint.quadrature(  # NOQA
                lambda x: func(x, y, *args),
                a,
                b,
                (),
                toli,
                rtoli,
                maxiteri,
                vec_func,
                miniteri,
            )[
                0
            ]  # NOQA

        # Is there a simple way to retrieve err info from the inner quad calls?

        val, err = sint.quadrature(
            outer_integrand, c, d, (), tol, rtol, maxitero, vec_func, minitero
        )

    elif method == "Gauss":

        # Gauss quadrature with orders equal to maxiters

        # Using lambda for readability
        def outer_integrand(y):
            return sp.integrate.fixed_quad(  # NOQA
                np.vectorize(lambda x: func(x, y, *args)), a, b, (), maxiteri
            )[
                0
            ]  # NOQA

        # Is there a simple way to retrieve err info from the inner quad calls?

        val, err = sp.integrate.fixed_quad(
            np.vectorize(outer_integrand), c, d, (), maxitero
        )

        assert err is None

    else:

        raise NotImplementedError("Unsupported quad method: " + method)

    return (val, err)


# }}}

# {{{ affine mappings


def solve_affine_map_2d(source_tria, target_tria):
    """Computes the affine map and its inverse that maps the source_tria to
       target_tria.

       :param source_tria: The triangle to be mapped.
       :type source_tria:
            tuple(tuple(float,float),tuple(float,float),tuple(float,float)).
       :param target_tria: The triangle to map to.
       :type target_tria:
            tuple(tuple(float,float),tuple(float,float),tuple(float,float)).

       :returns:
        - **mapping**: the forward map.
        - **J**: the Jacobian.
        - **invmap**: the inverse map.
        - **invJ**: the Jacobian of inverse map.
       :rtype:
        tuple(lambda, float, lambda, float)
    """
    assert len(source_tria) == 3
    for p in source_tria:
        assert len(p) == 2

    assert len(target_tria) == 3
    for p in target_tria:
        assert len(p) == 2

    # DOFs: A11, A12, A21, A22, b1, b2
    rhs = np.array(
        [
            target_tria[0][0],
            target_tria[0][1],
            target_tria[1][0],
            target_tria[1][1],
            target_tria[2][0],
            target_tria[2][1],
        ]
    )

    coef = np.array(
        [
            [source_tria[0][0], source_tria[0][1], 0, 0, 1, 0],
            [0, 0, source_tria[0][0], source_tria[0][1], 0, 1],
            [source_tria[1][0], source_tria[1][1], 0, 0, 1, 0],
            [0, 0, source_tria[1][0], source_tria[1][1], 0, 1],
            [source_tria[2][0], source_tria[2][1], 0, 0, 1, 0],
            [0, 0, source_tria[2][0], source_tria[2][1], 0, 1],
        ]
    )

    # x, residuals, _, _ = np.linalg.lstsq(coef, rhs)
    # assert (np.allclose(residuals, 0))
    try:
        x = np.linalg.solve(coef, rhs)
    except np.linalg.linalg.LinAlgError:
        print("")
        print("source:", source_tria)
        print("target:", target_tria)
        raise SystemExit("Error: Singular source triangle encountered")
    assert len(x) == 6
    assert np.allclose(np.dot(coef, x), rhs)

    a = np.array([[x[0], x[1]], [x[2], x[3]]])
    b = np.array([x[4], x[5]])

    # Using default value is the idiomatic way to "capture by value"
    mapping = lambda x, a=a, b=b: a.dot(np.array(x)) + b  # NOQA
    jacob = np.linalg.det(a)

    inva = np.linalg.inv(a)
    invb = -inva.dot(b)
    invmap = lambda x, a=inva, b=invb: inva.dot(np.array(x)) + invb  # NOQA
    inv_jacob = np.linalg.det(inva)

    assert np.abs(jacob * inv_jacob - 1) < 1e-12

    return (mapping, jacob, invmap, inv_jacob)


# }}}

# {{{ standard-triangle-to-rectangle mappings


def tria2rect_map_2d():
    """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.
       :rtype: tuple(lambda, lambda, lambda, lambda)
    """

    # (x,y) --> (rho, theta): T --> R
    def mapping(x):
        return (x[0] + x[1], np.arctan2(np.sqrt(x[1]), np.sqrt(x[0])))

    def jacob(x):
        return 1 / (2 * np.sqrt(x[0] * x[1]))

    # x = rho * cos^2(theta), y = rho * sin^2(theta)
    # J = rho * sin(2*theta)
    def invmap(u):
        return (u[0] * (np.cos(u[1]) ** 2), u[0] * (np.sin(u[1]) ** 2))

    def inv_jacob(u):
        return u[0] * np.sin(2 * u[1])

    return (mapping, jacob, invmap, inv_jacob)


def is_in_t(pt):
    """Checks if a point is in the template triangle T.

       :param pt: The point to be checked.
       :type pt: tuple(float,float).

       :returns: True if pt is in T.
       :rtype: bool.
    """
    flag = True
    if pt[0] < 0 or pt[1] < 0:
        flag = False
    if pt[0] + pt[1] > 1:
        flag = False
    return flag


def is_in_r(pt, a=0, b=1, c=0, d=np.pi / 2):
    """Checks if a point is in the (template) rectangle R.

       :param pt: The point to be checked.
       :type pt: tuple(float,float).

       :returns: True if pt is in R.
       :rtype: bool.
    """
    flag = True
    if pt[0] < a or pt[1] < c:
        flag = False
    if pt[0] > b or pt[1] > d:
        flag = False
    return flag


# }}}

# {{{ quadrature on arbitrary triangle


def is_collinear(p0, p1, p2):
    # v1 = p0 --> p1
    x1, y1 = p1[0] - p0[0], p1[1] - p0[1]
    # v2 = p0 --> p2
    x2, y2 = p2[0] - p0[0], p2[1] - p0[1]
    # v1 cross v2 == 0 <==> collinearity
    return np.abs(x1 * y2 - x2 * y1) < 1e-16


def is_positive_triangle(tria):
    p0 = tria[0]
    p1 = tria[1]
    p2 = tria[2]
    # v1 = p0 --> p1
    x1, y1 = p1[0] - p0[0], p1[1] - p0[1]
    # v2 = p0 --> p2
    x2, y2 = p2[0] - p0[0], p2[1] - p0[1]
    # v1 cross v2 > 0 <==> is positive
    return (x1 * y2 - x2 * y1) > 0


def tria_quad(
    func,
    tria,
    args=(),
    tol=1.49e-08,
    rtol=1.49e-08,
    maxiter=50,
    vec_func=True,
    miniter=1,
):
    """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 virtex of the tiangle. It is okay if func does not evaluate at the
    singular point. This function handles that automatically.

    :param func: A double variable Python function or method to integrate.
    :type func: function.
    :param tria: The triangular region to do quadrature.
    :type tria:
        tuple(tuple(float,float), tuple(float,float), tuple(float,float)).
    :param args: Extra arguments to pass to function.
    :type args: tuple, optional.
    :param tol: rtol Iteration stops when error between last two iterates is
        less than tol OR the relative change is less than rtol.
    :type tol: float, optional.
    :param rtol: Iteration stops when error between last two iterates is less
        than tol OR the relative change is less than rtol.
    :type rtol: float, optional.
    :param maxiter: Maximum order of Gaussian quadrature.
    :type maxiter: int, optional.
    :param vec_func: True if func handles arrays as arguments
        (is a "vector" function). Default is True.
    :type vec_func: bool, optional.
    :param miniter: Minimum order of Gaussian quadrature.
    :type miniter: int, optional.

    :returns:
        - **val**: Gaussian quadrature approximation (within tolerance)
            to integral.
        - **err**: Difference between last two estimates of the integral.
    :rtype: tuple(float,float).
    """

    assert len(tria) == 3
    for p in tria:
        assert len(p) == 2

    # Handle degenerate triangles
    if is_collinear(*tria):
        return (0.0, 0.0)

    # The function must be regular at the last two vertices
    assert np.isfinite(func(tria[1][0], tria[1][1], *args))
    assert np.isfinite(func(tria[2][0], tria[2][1], *args))

    # Solve for transforms
    template_tria = ((0, 0), (1, 0), (0, 1))
    afmp, j_afmp, inv_afmp, j_inv_afmp = solve_affine_map_2d(tria, template_tria)
    nlmp, j_nlmp, inv_nlmp, j_inv_nlmp = tria2rect_map_2d()

    # tria --> rect
    def mapping(x, y):
        return nlmp(afmp((x, y)))

    def jacobian(x, y):
        return j_afmp * j_nlmp(afmp((x, y)))

    # rect --> tria
    def inv_mapping(rho, theta):
        return inv_afmp(inv_nlmp((rho, theta)))

    def inv_jacobian(rho, theta):
        return j_inv_afmp * j_inv_nlmp((rho, theta))

    # Transformed function is defined on [0,1]X[0,pi/2]
    def transformed_func(rho, theta):
        preimage = inv_mapping(rho, theta)
        return func(preimage[0], preimage[1], *args)

    # Transformed function, when multiplied by jacobian, should have no
    # singularity (numerically special treatment still needed)

    # integrand = func * jacobian

    def integrand(rho, theta):
        prior = transformed_func(rho, theta) * inv_jacobian(rho, theta)
        # If something blows up, it is near the singular point
        if ~np.isfinite(prior):
            assert rho < 1e-3
            assert inv_jacobian(rho, theta) < 1e-6
            prior = 0
        return prior

    return qquad(
        func=integrand,
        a=0,
        b=1,
        c=0,
        d=np.pi / 2,
        args=(),
        tol=tol,
        rtol=rtol,
        maxiteri=maxiter,
        maxitero=maxiter,
        vec_func=False,
        miniteri=miniter,
        minitero=miniter,
    )


# }}}

# {{{ quadrature on a 2d box with a singular point inside


def 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,
):
    """Computes a (tensor product) double integral, with the integrand being
        singular at some point inside the region.

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

    :param func: A double variable Python function or method to integrate.
    :type func: function.
    :param a: Lower-left corner of integration region.
    :type a: float.
    :param b: Lower-right corner of integration region.
    :type b: float.
    :param c: Upper-left corner of integration region.
    :type c: float.
    :param d: Upper-right corner of integration region.
    :type d: float.
    :param singular_point: The singular point of the integrand func.
    :type singular_point: tuple(float, float).
    :param args: Extra arguments to pass to function.
    :type args: tuple, optional.
    :param tol: rtol Iteration stops when error between last two iterates is
                less than tol OR the relative change is less than rtol.
    :type tol: float, optional.
    :param rtol: Iteration stops when error between last two iterates is less
                than tol OR the relative change is less than rtol.
    :type rtol: float, optional.
    :param maxiter: Maximum order of Gaussian quadrature.
    :type maxiter: int, optional.
    :param vec_func: True if func handles arrays as arguments (is a "vector"
                function). Default is True.
    :type vec_func: bool, optional.
    :param miniter: Minimum order of Gaussian quadrature.
    :type miniter: int, optional.

    :returns:
        - **val**: Gaussian quadrature approximation (within tolerance) to integral.
        - **err**: Difference between last two estimates of the integral.
    :rtype: tuple(float,float).
    """
    box = ((a, c), (b, c), (b, d), (a, d))

    if not isinstance(singular_point, tuple):
        singular_point = (singular_point[0], singular_point[1])

    # When singular point is outside, project it onto the box bounday
    # This can import speed by not integrating around the actual singularity
    # when not necessary. (The splitting is still needed since it can be quite
    # close to singular).
    singular_point = (max(singular_point[0], a), max(singular_point[1], c))
    singular_point = (min(singular_point[0], b), min(singular_point[1], d))

    return quadri_quad(
        func, box, singular_point, args, tol, rtol, maxiter, vec_func, miniter
    )


# quadrilateral
def quadri_quad(
    func,
    quadrilateral,
    singular_point,
    args=(),
    tol=1.49e-08,
    rtol=1.49e-08,
    maxiter=50,
    vec_func=True,
    miniter=1,
):
    """Computes a double integral over a (non-twisted) quadrilateral, with the
        integrand being singular at some point inside the region.

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

    :param func: A double variable Python function or method to integrate.
    :type func: function.
    :param quadrilateral: The integration region.
    :type quadrilateral: tuple(tuple(float,float),
        tuple(float,float), tuple(float,float), tuple(float,float)).
    :param singular_point: The singular point of the integrand func.
    :type singular_point: tuple(float, float).
    :param args: Extra arguments to pass to function.
    :type args: tuple, optional.
    :param tol: rtol Iteration stops when error between last two iterates is
                less than tol OR the relative change is less than rtol.
    :type tol: float, optional.
    :param rtol: Iteration stops when error between last two iterates is less
                than tol OR the relative change is less than rtol.
    :type rtol: float, optional.
    :param maxiter: Maximum order of Gaussian quadrature.
    :type maxiter: int, optional.
    :param vec_func: True if func handles arrays as arguments (is a "vector"
                function). Default is True.
    :type vec_func: bool, optional.
    :param miniter: Minimum order of Gaussian quadrature.
    :type miniter: int, optional.

    :returns:
        - **val**: Gaussian quadrature approximation (within tolerance) to integral.
        - **err**: Difference between last two estimates of the integral.
    :rtype: tuple(float,float).
    """
    assert len(quadrilateral) == 4
    for p in quadrilateral:
        assert len(p) == 2

    # split the quadrilateral into four triangles
    trias = [
        (singular_point, quadrilateral[0], quadrilateral[1]),
        (singular_point, quadrilateral[1], quadrilateral[2]),
        (singular_point, quadrilateral[2], quadrilateral[3]),
        (singular_point, quadrilateral[3], quadrilateral[0]),
    ]

    for tria in trias:
        if not is_positive_triangle(tria):
            assert is_collinear(*tria)

    val = np.zeros(4)
    err = np.zeros(4)

    for i in range(4):
        val[i], err[i] = tria_quad(
            func, trias[i], args, tol, rtol, maxiter, vec_func, miniter
        )

    integral = np.sum(val)
    error = np.linalg.norm(err)

    return (integral, error)


# }}}
'''
class DesingularizationMapping:
    def __init__(self, nquad_points_1d):

def build_singular_box_quadrature(
        kernel,
        desing_mapping,

        ):
    """
    :arg kernel: an instance of :class:`sumpy.kernel.Kernel`
    :arg desing_mapping: an instance of :class:`sumpy.kernel.Kernel`
    """
'''

# vim: filetype=pyopencl.python:fdm=marker

<filename>volumential/singular_integral_2d.py
from __future__ import absolute_import, division, print_function

__copyright__ = "Copyright (C) 2017 - 2018 <NAME>"

__license__ = """
Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in
all copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
THE SOFTWARE.
"""

import logging

import numpy as np
import scipy as sp
import scipy.integrate as sint
from scipy.integrate import quadrature as quad  # noqa:F401

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

.. autofunction:: box_quad
"""

logger = logging.getLogger(__name__)

# {{{ quadrature on rectangular box with no singularity

quad_points_x = np.array([])
quad_points_y = np.array([])
quad_weights = np.array([])


def update_qquad_leggauss_formula(deg1, deg2):

    x1, w1 = np.polynomial.legendre.leggauss(deg1)
    x1 = (x1 + 1) / 2
    w1 = w1 / 2

    x2, w2 = np.polynomial.legendre.leggauss(deg2)
    x2 = (x2 + 1) / 2
    w2 = w2 / 2

    quad_points_x, quad_points_y = np.meshgrid(x1, x2)
    ww1, ww2 = np.meshgrid(w1, w2)
    global quad_weights
    quad_weights = ww1 * ww2


def 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",
):
    """Computes a (tensor product) double integral.

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

    :param func: A double variable Python function or method to integrate.
    :type func: function.
    :param a: Lower-left corner of integration region.
    :type a: float.
    :param b: Lower-right corner of integration region.
    :type b: float.
    :param c: Upper-left corner of integration region.
    :type c: float.
    :param d: Upper-right corner of integration region.
    :type d: float.
    :param args: Extra arguments to pass to function.
    :type args: tuple, optional.
    :param tol: rtol Iteration stops when error between last two iterates is
                less than tol OR the relative change is less than rtol.
    :type tol: float, optional.
    :param rtol: Iteration stops when error between last two iterates is less
                than tol OR the relative change is less than rtol.
    :type rtol: float, optional.
    :param maxitero: Maximum order of outer Gaussian quadrature.
    :type maxitero: int, optional.
    :param maxiteri: Maximum order of inner Gaussian quadrature.
    :type maxiteri: int, optional.
    :param vec_func: True if func handles arrays as arguments (is a "vector"
                function). Default is True.
    :type vec_func: bool, optional.
    :param minitero: Minimum order of outer Gaussian quadrature.
    :type minitero: int, optional.
    :param miniteri: Minimum order of inner Gaussian quadrature.
    :type miniteri: int, optional.

    :returns:
        - **val**: Gaussian quadrature approximation (within tolerance) to integral.
        - **err**: Difference between last two estimates of the integral.
    :rtype: tuple(float,float).
    """

    l1 = b - a
    l2 = d - c
    assert l1 > 0
    assert l2 > 0
    toli = tol / l2
    rtoli = rtol / l2

    if method == "Adaptive":

        # Using lambda for readability
        def outer_integrand(y):
            return sint.quadrature(  # NOQA
                lambda x: func(x, y, *args),
                a,
                b,
                (),
                toli,
                rtoli,
                maxiteri,
                vec_func,
                miniteri,
            )[
                0
            ]  # NOQA

        # Is there a simple way to retrieve err info from the inner quad calls?

        val, err = sint.quadrature(
            outer_integrand, c, d, (), tol, rtol, maxitero, vec_func, minitero
        )

    elif method == "Gauss":

        # Gauss quadrature with orders equal to maxiters

        # Using lambda for readability
        def outer_integrand(y):
            return sp.integrate.fixed_quad(  # NOQA
                np.vectorize(lambda x: func(x, y, *args)), a, b, (), maxiteri
            )[
                0
            ]  # NOQA

        # Is there a simple way to retrieve err info from the inner quad calls?

        val, err = sp.integrate.fixed_quad(
            np.vectorize(outer_integrand), c, d, (), maxitero
        )

        assert err is None

    else:

        raise NotImplementedError("Unsupported quad method: " + method)

    return (val, err)


# }}}

# {{{ affine mappings


def solve_affine_map_2d(source_tria, target_tria):
    """Computes the affine map and its inverse that maps the source_tria to
       target_tria.

       :param source_tria: The triangle to be mapped.
       :type source_tria:
            tuple(tuple(float,float),tuple(float,float),tuple(float,float)).
       :param target_tria: The triangle to map to.
       :type target_tria:
            tuple(tuple(float,float),tuple(float,float),tuple(float,float)).

       :returns:
        - **mapping**: the forward map.
        - **J**: the Jacobian.
        - **invmap**: the inverse map.
        - **invJ**: the Jacobian of inverse map.
       :rtype:
        tuple(lambda, float, lambda, float)
    """
    assert len(source_tria) == 3
    for p in source_tria:
        assert len(p) == 2

    assert len(target_tria) == 3
    for p in target_tria:
        assert len(p) == 2

    # DOFs: A11, A12, A21, A22, b1, b2
    rhs = np.array(
        [
            target_tria[0][0],
            target_tria[0][1],
            target_tria[1][0],
            target_tria[1][1],
            target_tria[2][0],
            target_tria[2][1],
        ]
    )

    coef = np.array(
        [
            [source_tria[0][0], source_tria[0][1], 0, 0, 1, 0],
            [0, 0, source_tria[0][0], source_tria[0][1], 0, 1],
            [source_tria[1][0], source_tria[1][1], 0, 0, 1, 0],
            [0, 0, source_tria[1][0], source_tria[1][1], 0, 1],
            [source_tria[2][0], source_tria[2][1], 0, 0, 1, 0],
            [0, 0, source_tria[2][0], source_tria[2][1], 0, 1],
        ]
    )

    # x, residuals, _, _ = np.linalg.lstsq(coef, rhs)
    # assert (np.allclose(residuals, 0))
    try:
        x = np.linalg.solve(coef, rhs)
    except np.linalg.linalg.LinAlgError:
        print("")
        print("source:", source_tria)
        print("target:", target_tria)
        raise SystemExit("Error: Singular source triangle encountered")
    assert len(x) == 6
    assert np.allclose(np.dot(coef, x), rhs)

    a = np.array([[x[0], x[1]], [x[2], x[3]]])
    b = np.array([x[4], x[5]])

    # Using default value is the idiomatic way to "capture by value"
    mapping = lambda x, a=a, b=b: a.dot(np.array(x)) + b  # NOQA
    jacob = np.linalg.det(a)

    inva = np.linalg.inv(a)
    invb = -inva.dot(b)
    invmap = lambda x, a=inva, b=invb: inva.dot(np.array(x)) + invb  # NOQA
    inv_jacob = np.linalg.det(inva)

    assert np.abs(jacob * inv_jacob - 1) < 1e-12

    return (mapping, jacob, invmap, inv_jacob)


# }}}

# {{{ standard-triangle-to-rectangle mappings


def tria2rect_map_2d():
    """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.
       :rtype: tuple(lambda, lambda, lambda, lambda)
    """

    # (x,y) --> (rho, theta): T --> R
    def mapping(x):
        return (x[0] + x[1], np.arctan2(np.sqrt(x[1]), np.sqrt(x[0])))

    def jacob(x):
        return 1 / (2 * np.sqrt(x[0] * x[1]))

    # x = rho * cos^2(theta), y = rho * sin^2(theta)
    # J = rho * sin(2*theta)
    def invmap(u):
        return (u[0] * (np.cos(u[1]) ** 2), u[0] * (np.sin(u[1]) ** 2))

    def inv_jacob(u):
        return u[0] * np.sin(2 * u[1])

    return (mapping, jacob, invmap, inv_jacob)


def is_in_t(pt):
    """Checks if a point is in the template triangle T.

       :param pt: The point to be checked.
       :type pt: tuple(float,float).

       :returns: True if pt is in T.
       :rtype: bool.
    """
    flag = True
    if pt[0] < 0 or pt[1] < 0:
        flag = False
    if pt[0] + pt[1] > 1:
        flag = False
    return flag


def is_in_r(pt, a=0, b=1, c=0, d=np.pi / 2):
    """Checks if a point is in the (template) rectangle R.

       :param pt: The point to be checked.
       :type pt: tuple(float,float).

       :returns: True if pt is in R.
       :rtype: bool.
    """
    flag = True
    if pt[0] < a or pt[1] < c:
        flag = False
    if pt[0] > b or pt[1] > d:
        flag = False
    return flag


# }}}

# {{{ quadrature on arbitrary triangle


def is_collinear(p0, p1, p2):
    # v1 = p0 --> p1
    x1, y1 = p1[0] - p0[0], p1[1] - p0[1]
    # v2 = p0 --> p2
    x2, y2 = p2[0] - p0[0], p2[1] - p0[1]
    # v1 cross v2 == 0 <==> collinearity
    return np.abs(x1 * y2 - x2 * y1) < 1e-16


def is_positive_triangle(tria):
    p0 = tria[0]
    p1 = tria[1]
    p2 = tria[2]
    # v1 = p0 --> p1
    x1, y1 = p1[0] - p0[0], p1[1] - p0[1]
    # v2 = p0 --> p2
    x2, y2 = p2[0] - p0[0], p2[1] - p0[1]
    # v1 cross v2 > 0 <==> is positive
    return (x1 * y2 - x2 * y1) > 0


def tria_quad(
    func,
    tria,
    args=(),
    tol=1.49e-08,
    rtol=1.49e-08,
    maxiter=50,
    vec_func=True,
    miniter=1,
):
    """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 virtex of the tiangle. It is okay if func does not evaluate at the
    singular point. This function handles that automatically.

    :param func: A double variable Python function or method to integrate.
    :type func: function.
    :param tria: The triangular region to do quadrature.
    :type tria:
        tuple(tuple(float,float), tuple(float,float), tuple(float,float)).
    :param args: Extra arguments to pass to function.
    :type args: tuple, optional.
    :param tol: rtol Iteration stops when error between last two iterates is
        less than tol OR the relative change is less than rtol.
    :type tol: float, optional.
    :param rtol: Iteration stops when error between last two iterates is less
        than tol OR the relative change is less than rtol.
    :type rtol: float, optional.
    :param maxiter: Maximum order of Gaussian quadrature.
    :type maxiter: int, optional.
    :param vec_func: True if func handles arrays as arguments
        (is a "vector" function). Default is True.
    :type vec_func: bool, optional.
    :param miniter: Minimum order of Gaussian quadrature.
    :type miniter: int, optional.

    :returns:
        - **val**: Gaussian quadrature approximation (within tolerance)
            to integral.
        - **err**: Difference between last two estimates of the integral.
    :rtype: tuple(float,float).
    """

    assert len(tria) == 3
    for p in tria:
        assert len(p) == 2

    # Handle degenerate triangles
    if is_collinear(*tria):
        return (0.0, 0.0)

    # The function must be regular at the last two vertices
    assert np.isfinite(func(tria[1][0], tria[1][1], *args))
    assert np.isfinite(func(tria[2][0], tria[2][1], *args))

    # Solve for transforms
    template_tria = ((0, 0), (1, 0), (0, 1))
    afmp, j_afmp, inv_afmp, j_inv_afmp = solve_affine_map_2d(tria, template_tria)
    nlmp, j_nlmp, inv_nlmp, j_inv_nlmp = tria2rect_map_2d()

    # tria --> rect
    def mapping(x, y):
        return nlmp(afmp((x, y)))

    def jacobian(x, y):
        return j_afmp * j_nlmp(afmp((x, y)))

    # rect --> tria
    def inv_mapping(rho, theta):
        return inv_afmp(inv_nlmp((rho, theta)))

    def inv_jacobian(rho, theta):
        return j_inv_afmp * j_inv_nlmp((rho, theta))

    # Transformed function is defined on [0,1]X[0,pi/2]
    def transformed_func(rho, theta):
        preimage = inv_mapping(rho, theta)
        return func(preimage[0], preimage[1], *args)

    # Transformed function, when multiplied by jacobian, should have no
    # singularity (numerically special treatment still needed)

    # integrand = func * jacobian

    def integrand(rho, theta):
        prior = transformed_func(rho, theta) * inv_jacobian(rho, theta)
        # If something blows up, it is near the singular point
        if ~np.isfinite(prior):
            assert rho < 1e-3
            assert inv_jacobian(rho, theta) < 1e-6
            prior = 0
        return prior

    return qquad(
        func=integrand,
        a=0,
        b=1,
        c=0,
        d=np.pi / 2,
        args=(),
        tol=tol,
        rtol=rtol,
        maxiteri=maxiter,
        maxitero=maxiter,
        vec_func=False,
        miniteri=miniter,
        minitero=miniter,
    )


# }}}

# {{{ quadrature on a 2d box with a singular point inside


def 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,
):
    """Computes a (tensor product) double integral, with the integrand being
        singular at some point inside the region.

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

    :param func: A double variable Python function or method to integrate.
    :type func: function.
    :param a: Lower-left corner of integration region.
    :type a: float.
    :param b: Lower-right corner of integration region.
    :type b: float.
    :param c: Upper-left corner of integration region.
    :type c: float.
    :param d: Upper-right corner of integration region.
    :type d: float.
    :param singular_point: The singular point of the integrand func.
    :type singular_point: tuple(float, float).
    :param args: Extra arguments to pass to function.
    :type args: tuple, optional.
    :param tol: rtol Iteration stops when error between last two iterates is
                less than tol OR the relative change is less than rtol.
    :type tol: float, optional.
    :param rtol: Iteration stops when error between last two iterates is less
                than tol OR the relative change is less than rtol.
    :type rtol: float, optional.
    :param maxiter: Maximum order of Gaussian quadrature.
    :type maxiter: int, optional.
    :param vec_func: True if func handles arrays as arguments (is a "vector"
                function). Default is True.
    :type vec_func: bool, optional.
    :param miniter: Minimum order of Gaussian quadrature.
    :type miniter: int, optional.

    :returns:
        - **val**: Gaussian quadrature approximation (within tolerance) to integral.
        - **err**: Difference between last two estimates of the integral.
    :rtype: tuple(float,float).
    """
    box = ((a, c), (b, c), (b, d), (a, d))

    if not isinstance(singular_point, tuple):
        singular_point = (singular_point[0], singular_point[1])

    # When singular point is outside, project it onto the box bounday
    # This can import speed by not integrating around the actual singularity
    # when not necessary. (The splitting is still needed since it can be quite
    # close to singular).
    singular_point = (max(singular_point[0], a), max(singular_point[1], c))
    singular_point = (min(singular_point[0], b), min(singular_point[1], d))

    return quadri_quad(
        func, box, singular_point, args, tol, rtol, maxiter, vec_func, miniter
    )


# quadrilateral
def quadri_quad(
    func,
    quadrilateral,
    singular_point,
    args=(),
    tol=1.49e-08,
    rtol=1.49e-08,
    maxiter=50,
    vec_func=True,
    miniter=1,
):
    """Computes a double integral over a (non-twisted) quadrilateral, with the
        integrand being singular at some point inside the region.

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

    :param func: A double variable Python function or method to integrate.
    :type func: function.
    :param quadrilateral: The integration region.
    :type quadrilateral: tuple(tuple(float,float),
        tuple(float,float), tuple(float,float), tuple(float,float)).
    :param singular_point: The singular point of the integrand func.
    :type singular_point: tuple(float, float).
    :param args: Extra arguments to pass to function.
    :type args: tuple, optional.
    :param tol: rtol Iteration stops when error between last two iterates is
                less than tol OR the relative change is less than rtol.
    :type tol: float, optional.
    :param rtol: Iteration stops when error between last two iterates is less
                than tol OR the relative change is less than rtol.
    :type rtol: float, optional.
    :param maxiter: Maximum order of Gaussian quadrature.
    :type maxiter: int, optional.
    :param vec_func: True if func handles arrays as arguments (is a "vector"
                function). Default is True.
    :type vec_func: bool, optional.
    :param miniter: Minimum order of Gaussian quadrature.
    :type miniter: int, optional.

    :returns:
        - **val**: Gaussian quadrature approximation (within tolerance) to integral.
        - **err**: Difference between last two estimates of the integral.
    :rtype: tuple(float,float).
    """
    assert len(quadrilateral) == 4
    for p in quadrilateral:
        assert len(p) == 2

    # split the quadrilateral into four triangles
    trias = [
        (singular_point, quadrilateral[0], quadrilateral[1]),
        (singular_point, quadrilateral[1], quadrilateral[2]),
        (singular_point, quadrilateral[2], quadrilateral[3]),
        (singular_point, quadrilateral[3], quadrilateral[0]),
    ]

    for tria in trias:
        if not is_positive_triangle(tria):
            assert is_collinear(*tria)

    val = np.zeros(4)
    err = np.zeros(4)

    for i in range(4):
        val[i], err[i] = tria_quad(
            func, trias[i], args, tol, rtol, maxiter, vec_func, miniter
        )

    integral = np.sum(val)
    error = np.linalg.norm(err)

    return (integral, error)


# }}}
'''
class DesingularizationMapping:
    def __init__(self, nquad_points_1d):

def build_singular_box_quadrature(
        kernel,
        desing_mapping,

        ):
    """
    :arg kernel: an instance of :class:`sumpy.kernel.Kernel`
    :arg desing_mapping: an instance of :class:`sumpy.kernel.Kernel`
    """
'''

# vim: filetype=pyopencl.python:fdm=marker

Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:  The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. # noqa:F401 The 2D singular integrals are computed using the transform described in http://link.springer.com/10.1007/BF00370482.  .. autofunction:: box_quad # {{{ quadrature on rectangular box with no singularity Computes a (tensor product) double integral.      Integrate func on [a, b]X[c, d] using Gaussian quadrature with absolute     tolerance tol.      :param func: A double variable Python function or method to integrate.     :type func: function.     :param a: Lower-left corner of integration region.     :type a: float.     :param b: Lower-right corner of integration region.     :type b: float.     :param c: Upper-left corner of integration region.     :type c: float.     :param d: Upper-right corner of integration region.     :type d: float.     :param args: Extra arguments to pass to function.     :type args: tuple, optional.     :param tol: rtol Iteration stops when error between last two iterates is                 less than tol OR the relative change is less than rtol.     :type tol: float, optional.     :param rtol: Iteration stops when error between last two iterates is less                 than tol OR the relative change is less than rtol.     :type rtol: float, optional.     :param maxitero: Maximum order of outer Gaussian quadrature.     :type maxitero: int, optional.     :param maxiteri: Maximum order of inner Gaussian quadrature.     :type maxiteri: int, optional.     :param vec_func: True if func handles arrays as arguments (is a "vector"                 function). Default is True.     :type vec_func: bool, optional.     :param minitero: Minimum order of outer Gaussian quadrature.     :type minitero: int, optional.     :param miniteri: Minimum order of inner Gaussian quadrature.     :type miniteri: int, optional.      :returns:         - **val**: Gaussian quadrature approximation (within tolerance) to integral.         - **err**: Difference between last two estimates of the integral.     :rtype: tuple(float,float). # Using lambda for readability # NOQA # NOQA # Is there a simple way to retrieve err info from the inner quad calls? # Gauss quadrature with orders equal to maxiters # Using lambda for readability # NOQA # NOQA # Is there a simple way to retrieve err info from the inner quad calls? # }}} # {{{ affine mappings Computes the affine map and its inverse that maps the source_tria to        target_tria.         :param source_tria: The triangle to be mapped.        :type source_tria:             tuple(tuple(float,float),tuple(float,float),tuple(float,float)).        :param target_tria: The triangle to map to.        :type target_tria:             tuple(tuple(float,float),tuple(float,float),tuple(float,float)).         :returns:         - **mapping**: the forward map.         - **J**: the Jacobian.         - **invmap**: the inverse map.         - **invJ**: the Jacobian of inverse map.        :rtype:         tuple(lambda, float, lambda, float) # DOFs: A11, A12, A21, A22, b1, b2 # x, residuals, _, _ = np.linalg.lstsq(coef, rhs) # assert (np.allclose(residuals, 0)) # Using default value is the idiomatic way to "capture by value" # NOQA # NOQA # }}} # {{{ standard-triangle-to-rectangle mappings 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.        :rtype: tuple(lambda, lambda, lambda, lambda) # (x,y) --> (rho, theta): T --> R # x = rho * cos^2(theta), y = rho * sin^2(theta) # J = rho * sin(2*theta) Checks if a point is in the template triangle T.         :param pt: The point to be checked.        :type pt: tuple(float,float).         :returns: True if pt is in T.        :rtype: bool. Checks if a point is in the (template) rectangle R.         :param pt: The point to be checked.        :type pt: tuple(float,float).         :returns: True if pt is in R.        :rtype: bool. # }}} # {{{ quadrature on arbitrary triangle # v1 = p0 --> p1 # v2 = p0 --> p2 # v1 cross v2 == 0 <==> collinearity # v1 = p0 --> p1 # v2 = p0 --> p2 # v1 cross v2 > 0 <==> is positive 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 virtex of the tiangle. It is okay if func does not evaluate at the     singular point. This function handles that automatically.      :param func: A double variable Python function or method to integrate.     :type func: function.     :param tria: The triangular region to do quadrature.     :type tria:         tuple(tuple(float,float), tuple(float,float), tuple(float,float)).     :param args: Extra arguments to pass to function.     :type args: tuple, optional.     :param tol: rtol Iteration stops when error between last two iterates is         less than tol OR the relative change is less than rtol.     :type tol: float, optional.     :param rtol: Iteration stops when error between last two iterates is less         than tol OR the relative change is less than rtol.     :type rtol: float, optional.     :param maxiter: Maximum order of Gaussian quadrature.     :type maxiter: int, optional.     :param vec_func: True if func handles arrays as arguments         (is a "vector" function). Default is True.     :type vec_func: bool, optional.     :param miniter: Minimum order of Gaussian quadrature.     :type miniter: int, optional.      :returns:         - **val**: Gaussian quadrature approximation (within tolerance)             to integral.         - **err**: Difference between last two estimates of the integral.     :rtype: tuple(float,float). # Handle degenerate triangles # The function must be regular at the last two vertices # Solve for transforms # tria --> rect # rect --> tria # Transformed function is defined on [0,1]X[0,pi/2] # Transformed function, when multiplied by jacobian, should have no # singularity (numerically special treatment still needed) # integrand = func * jacobian # If something blows up, it is near the singular point # }}} # {{{ quadrature on a 2d box with a singular point inside Computes a (tensor product) double integral, with the integrand being         singular at some point inside the region.      Integrate func on [a, b]X[c, d] using transformed Gaussian quadrature with     absolute tolerance tol.      :param func: A double variable Python function or method to integrate.     :type func: function.     :param a: Lower-left corner of integration region.     :type a: float.     :param b: Lower-right corner of integration region.     :type b: float.     :param c: Upper-left corner of integration region.     :type c: float.     :param d: Upper-right corner of integration region.     :type d: float.     :param singular_point: The singular point of the integrand func.     :type singular_point: tuple(float, float).     :param args: Extra arguments to pass to function.     :type args: tuple, optional.     :param tol: rtol Iteration stops when error between last two iterates is                 less than tol OR the relative change is less than rtol.     :type tol: float, optional.     :param rtol: Iteration stops when error between last two iterates is less                 than tol OR the relative change is less than rtol.     :type rtol: float, optional.     :param maxiter: Maximum order of Gaussian quadrature.     :type maxiter: int, optional.     :param vec_func: True if func handles arrays as arguments (is a "vector"                 function). Default is True.     :type vec_func: bool, optional.     :param miniter: Minimum order of Gaussian quadrature.     :type miniter: int, optional.      :returns:         - **val**: Gaussian quadrature approximation (within tolerance) to integral.         - **err**: Difference between last two estimates of the integral.     :rtype: tuple(float,float). # When singular point is outside, project it onto the box bounday # This can import speed by not integrating around the actual singularity # when not necessary. (The splitting is still needed since it can be quite # close to singular). # quadrilateral Computes a double integral over a (non-twisted) quadrilateral, with the         integrand being singular at some point inside the region.      Integrate func on [a, b]X[c, d] using transformed Gaussian quadrature with     absolute tolerance tol.      :param func: A double variable Python function or method to integrate.     :type func: function.     :param quadrilateral: The integration region.     :type quadrilateral: tuple(tuple(float,float),         tuple(float,float), tuple(float,float), tuple(float,float)).     :param singular_point: The singular point of the integrand func.     :type singular_point: tuple(float, float).     :param args: Extra arguments to pass to function.     :type args: tuple, optional.     :param tol: rtol Iteration stops when error between last two iterates is                 less than tol OR the relative change is less than rtol.     :type tol: float, optional.     :param rtol: Iteration stops when error between last two iterates is less                 than tol OR the relative change is less than rtol.     :type rtol: float, optional.     :param maxiter: Maximum order of Gaussian quadrature.     :type maxiter: int, optional.     :param vec_func: True if func handles arrays as arguments (is a "vector"                 function). Default is True.     :type vec_func: bool, optional.     :param miniter: Minimum order of Gaussian quadrature.     :type miniter: int, optional.      :returns:         - **val**: Gaussian quadrature approximation (within tolerance) to integral.         - **err**: Difference between last two estimates of the integral.     :rtype: tuple(float,float). # split the quadrilateral into four triangles # }}} class DesingularizationMapping:     def __init__(self, nquad_points_1d):  def build_singular_box_quadrature(         kernel,         desing_mapping,          ):     """     :arg kernel: an instance of :class:`sumpy.kernel.Kernel`     :arg desing_mapping: an instance of :class:`sumpy.kernel.Kernel`     """ # vim: filetype=pyopencl.python:fdm=marker