diff --git a/modepy/matrices.py b/modepy/matrices.py index 5c0fbc494745baf78fdc7c5f94c77d8fa36e320a..14a77134457212e6910b826adda2ab34af9c6a56 100644 --- a/modepy/matrices.py +++ b/modepy/matrices.py @@ -235,10 +235,8 @@ class _FaceMap(object): self.face_dim = vol_dim - 1 def __call__(self, points): - return ( - self.origin[:, np.newaxis] - + - np.einsum("ad,dn->an", self.span, points*0.5 + 0.5)) + return (self.origin[:, np.newaxis] + + np.einsum("ad,dn->an", self.span, points*0.5 + 0.5)) def modal_face_mass_matrix(trial_basis, order, face_vertices, test_basis=None): diff --git a/modepy/nodes.py b/modepy/nodes.py index 911640c9572d6c0b86c03727175531cf444c7bb1..3ee11feb7f39ff7bac99445e5006c1a0d36778c4 100644 --- a/modepy/nodes.py +++ b/modepy/nodes.py @@ -146,7 +146,7 @@ def warp_and_blend_nodes_2d(n, node_tuples=None): # {{{ 3D nodes _alpha_opt_3d = [ - 0, 0, 0, 0.1002, 1.1332, 1.5608, 1.3413, 1.2577, 1.1603, + 0, 0, 0, 0.1002, 1.1332, 1.5608, 1.3413, 1.2577, 1.1603, 1.10153, 0.6080, 0.4523, 0.8856, 0.8717, 0.9655] diff --git a/modepy/quadrature/__init__.py b/modepy/quadrature/__init__.py index b53cc34d69f685750d5dc26a51644cfae1944a41..8f3c733eacb99d5466771490057bd7d842e67ddd 100644 --- a/modepy/quadrature/__init__.py +++ b/modepy/quadrature/__init__.py @@ -45,7 +45,12 @@ class Quadrature(object): .. attribute :: weights A vector of length *nnodes*. + + .. attribute :: exact_to + + Total polynomial degree up to which the quadrature is exact. """ + def __init__(self, nodes, weights): self.nodes = nodes self.weights = weights @@ -61,7 +66,7 @@ class Quadrature(object): class Transformed1DQuadrature(Quadrature): - """A quadrature rule on an arbitrary interval :math:`(a,b)`. """ + """A quadrature rule on an arbitrary interval :math:`(a, b)`.""" def __init__(self, quad, left, right): """Transform a given 1D quadrature rule *quad* onto the @@ -70,9 +75,10 @@ class Transformed1DQuadrature(Quadrature): self.left = left self.right = right - length = right-left - assert length > 0 + length = right - left half_length = length / 2 + assert length > 0 + Quadrature.__init__(self, - left + (quad.nodes+1)/2*length, - quad.weights*half_length) + left + (quad.nodes+1) / 2 * length, + quad.weights * half_length) diff --git a/modepy/quadrature/grundmann_moeller.py b/modepy/quadrature/grundmann_moeller.py index 09c2a73c3606faad34562d863ba885c4f09b9776..1e2d6ca4b1a7a4052e73b0ddcc22e5f35a845544 100644 --- a/modepy/quadrature/grundmann_moeller.py +++ b/modepy/quadrature/grundmann_moeller.py @@ -49,7 +49,7 @@ def _extended_euclidean(q, r): quot, t = divmod(q, r) T = Q[0] - quot*R[0], Q[1] - quot*R[1] # noqa: N806 q, r = r, t - Q, R = R, T + Q, R = R, T # noqa: N806 return q, Q[0], Q[1] @@ -68,7 +68,7 @@ class GrundmannMoellerSimplexQuadrature(Quadrature): r"""Cubature on an *n*-simplex. This cubature rule has both negative and positive weights. - It is exact for polynomials up to order :math:`2s+1`, where + It is exact for polynomials up to order :math:`2s + 1`, where :math:`s` is given as *order*. The integration domain is the unit simplex. (see :ref:`tri-coords` @@ -82,7 +82,7 @@ class GrundmannMoellerSimplexQuadrature(Quadrature): * A. Grundmann and H.M. Moeller, Invariant integration formulas for the n-simplex by combinatorial methods, - SIAM J. Numer. Anal. 15 (1978), 282--290. + SIAM J. Numer. Anal. 15 (1978), 282--290. http://dx.doi.org/10.1137/0715019 """ @@ -96,7 +96,7 @@ class GrundmannMoellerSimplexQuadrature(Quadrature): """ s = order n = dimension - d = 2*s+1 + d = 2*s + 1 self.exact_to = d @@ -115,17 +115,17 @@ class GrundmannMoellerSimplexQuadrature(Quadrature): points_to_weights = {} - for i in range(s+1): + for i in range(s + 1): weight = (-1)**i * 2**(-2*s) \ - * (d + n-2*i)**d \ + * (d + n - 2*i)**d \ / factorial(i) \ - / factorial(d+n-i) + / factorial(d + n - i) - for t in generate_decreasing_nonnegative_tuples_summing_to(s-i, n+1): + for t in generate_decreasing_nonnegative_tuples_summing_to(s - i, n + 1): for beta in generate_unique_permutations(t): - denominator = d+n-2*i + denominator = d + n - 2*i point = tuple( - _simplify_fraction((2*beta_i+1, denominator)) + _simplify_fraction((2*beta_i + 1, denominator)) for beta_i in beta) points_to_weights[point] = \ @@ -133,17 +133,17 @@ class GrundmannMoellerSimplexQuadrature(Quadrature): from operator import add - vertices = [-1 * np.ones((n,))] \ + vertices = ([-1 * np.ones((n,))] + [np.array(x) - for x in wandering_element(n, landscape=-1, wanderer=1)] + for x in wandering_element(n, landscape=-1, wanderer=1)]) nodes = [] weights = [] dim_factor = 2**n for p, w in six.iteritems(points_to_weights): - real_p = reduce(add, (a/b*v for (a, b), v in zip(p, vertices))) + real_p = reduce(add, (a/b * v for (a, b), v in zip(p, vertices))) nodes.append(real_p) - weights.append(dim_factor*w) + weights.append(dim_factor * w) Quadrature.__init__(self, np.array(nodes).T, np.array(weights)) diff --git a/modepy/quadrature/jacobi_gauss.py b/modepy/quadrature/jacobi_gauss.py index ab268e14a826818d40b9b56e0c1e5aef1a78885a..91d3f3afcecc67a5bf2668aec9bde5091df0a4e8 100644 --- a/modepy/quadrature/jacobi_gauss.py +++ b/modepy/quadrature/jacobi_gauss.py @@ -31,25 +31,33 @@ from modepy.quadrature import Quadrature class JacobiGaussQuadrature(Quadrature): r"""An Gauss quadrature of order *N* associated with the - Jacobi weight :math:`(1-x)^\alpha(1+x)^\beta`. + Jacobi weight :math:`(1 - x)^\alpha (1 + x)^\beta`. - Assumes :math:`\alpha,\beta > -1` with - :math:`(\alpha,\beta)\not\in\{(-1/2,-1/2)\}`. + Assumes :math:`\alpha, \beta > -1` with + :math:`(\alpha, \beta) \not \in \{(-1/2, -1/2)\}`. - Integrates on the interval (-1,1). - The quadrature rule is exact up to degree :math:`2N+1`. + Integrates on the interval :math:`(-1, 1)`. + The quadrature rule is exact up to degree :math:`2N + 1`. Inherits from :class:`modepy.Quadrature`. See there for the interface to obtain nodes and weights. """ + def __init__(self, alpha, beta, N): # noqa x, w = self.compute_weights_and_nodes(N, alpha, beta) + self.exact_to = 2*N + 1 + Quadrature.__init__(self, x, w) @staticmethod def compute_weights_and_nodes(N, alpha, beta): # noqa - """Return (nodes, weights) for an n-th order Gauss quadrature - with the Jacobi polynomials of type (alpha, beta). + """ + :arg N: order of the Gauss quadrature (the order of exactly + integrated polynomials is :math:`2 N + 1`). + :arg alpha: power of :math:`1 - x` in the Jacobi polynomial weight. + :arg beta: power of :math:`1 + x` in the Jacobi polynomial weight. + + :return: a tuple ``(nodes, weights)`` of quadrature notes and weights. """ # FIXME: This could stand to be upgraded to the Rokhlin algorithm. @@ -68,42 +76,28 @@ class JacobiGaussQuadrature(Quadrature): # see Appendix A of Hesthaven/Warburton for these formulas def a(n): - return ( - 2/(2*n+apb) - * - sqrt( - (n*(n+apb)*(n+alpha)*(n+beta)) - / - ((2*n+apb-1)*(2*n+apb+1)) - ) - ) + return (2 / (2*n + apb) + * sqrt((n * (n+apb) * (n+alpha) * (n+beta)) + / ((2*n + apb - 1) * (2*n + apb + 1)))) def b(n): if n == 0: - return ( - -(alpha-beta) - / - (apb+2) - ) + return -(alpha-beta) / (apb+2) else: - return ( - -(alpha**2-beta**2) - / - ((2*n+apb)*(2*n+apb+2)) - ) + return -(alpha**2 - beta**2) / ((2*n + apb) * (2*n + apb + 2)) - T = np.zeros((N+1, N+1)) # noqa + T = np.zeros((N + 1, N + 1)) # noqa: N806 - for n in range(N+1): + for n in range(N + 1): T[n, n] = b(n) if n > 0: - T[n, n-1] = current_a # noqa + T[n, n-1] = current_a # noqa: F821 if n < N: - next_a = a(n+1) + next_a = a(n + 1) T[n, n+1] = next_a - current_a = next_a # noqa + current_a = next_a # noqa: F841 - assert la.norm(T-T.T) < 1e-12 + assert la.norm(T - T.T) < 1.0e-12 eigval, eigvec = la.eigh(T) assert la.norm(np.dot(T, eigvec) - np.dot(eigvec, np.diag(eigval))) < 1e-12 @@ -113,7 +107,7 @@ class JacobiGaussQuadrature(Quadrature): p0 = partial(jacobi, alpha, beta, 0) # that's a constant, sure nodes = eigval weights = np.array( - [eigvec[0, i]**2 / p0(nodes[i])**2 for i in range(N+1)]) + [eigvec[0, i]**2 / p0(nodes[i])**2 for i in range(N + 1)]) return nodes, weights @@ -121,8 +115,8 @@ class JacobiGaussQuadrature(Quadrature): class LegendreGaussQuadrature(JacobiGaussQuadrature): """An Gauss quadrature associated with weight 1. - Integrates on the interval (-1,1). - The quadrature rule is exact up to degree :math:`2N+1`. + Integrates on the interval :math:`(-1, 1)`. + The quadrature rule is exact up to degree :math:`2N + 1`. Inherits from :class:`modepy.Quadrature`. See there for the interface to obtain nodes and weights. @@ -137,25 +131,23 @@ def jacobi_gauss_lobatto_nodes(alpha, beta, N): # noqa nodes corresponding to the :class:`JacobiGaussQuadrature` with the same parameters. - Exact to degree :math:`2N-3`. + Exact to degree :math:`2N - 3`. """ - x = np.zeros((N+1,)) + x = np.zeros((N + 1,)) x[0] = -1 x[-1] = 1 if N == 1: return x - x[1:-1] = np.array( - JacobiGaussQuadrature(alpha+1, beta+1, N-2).nodes - ).real + x[1:-1] = np.array(JacobiGaussQuadrature(alpha + 1, beta + 1, N - 2).nodes).real return x def legendre_gauss_lobatto_nodes(N): # noqa """Compute the Legendre-Gauss-Lobatto quadrature nodes. - Exact to degree :math:`2N-1`. + Exact to degree :math:`2N - 1`. """ return jacobi_gauss_lobatto_nodes(0, 0, N) diff --git a/modepy/quadrature/vr_quad_data_tet.py b/modepy/quadrature/vr_quad_data_tet.py index df8ba46ce5e11c677fdba68b51f6303b28c8fd52..a41d11f5ab928a1b24b1e92615bb82bb6a30ef6e 100644 --- a/modepy/quadrature/vr_quad_data_tet.py +++ b/modepy/quadrature/vr_quad_data_tet.py @@ -23,6 +23,7 @@ def process_rule(table): return result + tetrahedron_data = process_rule({ # noqa 0: {points: [[0.00000000000000000000000000000000e+00], [0.00000000000000000000000000000000e+00], diff --git a/modepy/quadrature/vr_quad_data_tri.py b/modepy/quadrature/vr_quad_data_tri.py index 5c58da5013363ad0af4514c19f5b0cdf4a776e18..2e890f8ac9b5edb81ab13315c7f9fe0d093242c0 100644 --- a/modepy/quadrature/vr_quad_data_tri.py +++ b/modepy/quadrature/vr_quad_data_tri.py @@ -23,6 +23,7 @@ def process_rule(table): return result + triangle_data = process_rule({ # noqa 0: {points: [[0.00000000000000000000000000000000e+00], [0.00000000000000000000000000000000e+00], diff --git a/modepy/tools.py b/modepy/tools.py index 4989c00cbcb69a8d2d22d5561ee7a122550440d0..1656192a92a98d02e6f98a656a0ed8decaf2250a 100644 --- a/modepy/tools.py +++ b/modepy/tools.py @@ -100,10 +100,9 @@ class Monomial: from pytools import factorial from operator import mul - return (self.factor*2**len(self.exponents) * - reduce(mul, (factorial(alpha) for alpha in self.exponents)) - / - factorial(len(self.exponents)+sum(self.exponents))) + return (self.factor * 2**len(self.exponents) + * reduce(mul, (factorial(alpha) for alpha in self.exponents)) + / factorial(len(self.exponents)+sum(self.exponents))) def diff(self, coordinate): diff_exp = list(self.exponents) @@ -146,14 +145,14 @@ class AffineMap: EQUILATERAL_TO_UNIT_MAP = { 1: AffineMap([[1]], [0]), 2: AffineMap([ - [1, -1 / sqrt(3)], - [0, 2 / sqrt(3)]], - [-1/3, -1/3]), + [1, -1/sqrt(3)], + [0, 2/sqrt(3)]], + [-1/3, -1/3]), 3: AffineMap([ [1, -1/sqrt(3), -1/sqrt(6)], [0, 2/sqrt(3), -1/sqrt(6)], [0, 0, sqrt(6)/2]], - [-1/2, -1/2, -1/2]) + [-1/2, -1/2, -1/2]) } @@ -206,14 +205,14 @@ EQUILATERAL_VERTICES = { ]), 2: np.array([ [-1, -1/sqrt(3)], - [1, -1/sqrt(3)], - [0, 2/sqrt(3)], + [1, -1/sqrt(3)], + [0, 2/sqrt(3)], ]), 3: np.array([ [-1, -1/sqrt(3), -1/sqrt(6)], - [1, -1/sqrt(3), -1/sqrt(6)], - [0, 2/sqrt(3), -1/sqrt(6)], - [0, 0, 3/sqrt(6)], + [1, -1/sqrt(3), -1/sqrt(6)], + [0, 2/sqrt(3), -1/sqrt(6)], + [0, 0, 3/sqrt(6)], ]) } @@ -232,7 +231,7 @@ def pick_random_simplex_unit_coordinate(rng, dims): r = np.zeros(dims, np.float64) for j in range(dims): rn = rng.uniform(0, remaining) - r[j] = base+rn + r[j] = base + rn remaining -= rn return r diff --git a/setup.cfg b/setup.cfg index 7d2560e5b39a70a6c64d1da77c6ce1543ac0f5c2..9328bab7fdc59d71722ef389044edb65f327c10a 100644 --- a/setup.cfg +++ b/setup.cfg @@ -1,3 +1,3 @@ [flake8] -ignore = E126,E127,E128,E123,E226,E241,E242,E402,W503 +ignore = E123,E126,E127,E128,E226,E241,E242,E402,W503 max-line-length=85 diff --git a/test/test_quadrature.py b/test/test_quadrature.py index 4f395686834555d6aad67b32eba65047f43c627c..07ba925479c56ca57faf8d3bd045d43a22ccd105 100644 --- a/test/test_quadrature.py +++ b/test/test_quadrature.py @@ -33,31 +33,33 @@ def test_transformed_quadrature(): """Test 1D quadrature on arbitrary intervals""" def gaussian_density(x, mu, sigma): - return 1/(sigma*np.sqrt(2*np.pi))*np.exp(-(x-mu)**2/(2*sigma**2)) + return 1 / (sigma * np.sqrt(2*np.pi)) * np.exp(-(x-mu)**2 / (2 * sigma**2)) from modepy.quadrature import Transformed1DQuadrature from modepy.quadrature.jacobi_gauss import LegendreGaussQuadrature mu = 17 sigma = 12 - tq = Transformed1DQuadrature(LegendreGaussQuadrature(20), mu-6*sigma, mu+6*sigma) + tq = Transformed1DQuadrature(LegendreGaussQuadrature(20), + mu - 6*sigma, mu + 6*sigma) result = tq(lambda x: gaussian_density(x, mu, sigma)) - assert abs(result - 1) < 1e-9 + assert abs(result - 1) < 1.0e-9 def test_gauss_quadrature(): from modepy.quadrature.jacobi_gauss import LegendreGaussQuadrature - for s in range(9+1): - cub = LegendreGaussQuadrature(s) - for deg in range(2*s+1 + 1): + for s in range(9 + 1): + quad = LegendreGaussQuadrature(s) + for deg in range(quad.exact_to + 1): def f(x): return x**deg - i_f = cub(f) - i_f_true = 1/(deg+1)*(1-(-1)**(deg+1)) + + i_f = quad(f) + i_f_true = 1 / (deg+1) * (1 - (-1)**(deg + 1)) err = abs(i_f - i_f_true) - assert err < 2e-15 + assert err < 2.0e-15, (s, deg, err, i_f, i_f_true) @pytest.mark.parametrize(("quad_class", "highest_order"), [ @@ -77,7 +79,7 @@ def test_simplex_quadrature(quad_class, highest_order, dim): try: quad = quad_class(order, dim) except mp.QuadratureRuleUnavailable: - print(("UNAVAIL", quad_class, order)) + print(("UNAVAILABLE", quad_class, order)) break if isinstance(quad_class, mp.VioreanuRokhlinSimplexQuadrature):