Skip to content

Commit

Permalink
Fix GaussJacobiLineQuadrature
Browse files Browse the repository at this point in the history
  • Loading branch information
pbrubeck committed Nov 26, 2023
1 parent f6287d4 commit c70914c
Showing 1 changed file with 4 additions and 20 deletions.
24 changes: 4 additions & 20 deletions FIAT/quadrature.py
Original file line number Diff line number Diff line change
Expand Up @@ -40,10 +40,10 @@ class GaussJacobiQuadratureLineRule(QuadratureRule):
"""Gauss-Jacobi quadature rule determined by Jacobi weights a and b
using m roots of m:th order Jacobi polynomial."""

def __init__(self, ref_el, m):
def __init__(self, ref_el, m, a=0, b=0):
# this gives roots on the default (-1,1) reference element
# (xs_ref, ws_ref) = gaussjacobi(m, a, b)
(xs_ref, ws_ref) = gaussjacobi(m, 0., 0.)
xs_ref, ws_ref = gaussjacobi(m, a, b)

Ref1 = reference_element.DefaultLine()
A, b = reference_element.make_affine_mapping(Ref1.get_vertices(),
Expand Down Expand Up @@ -96,31 +96,15 @@ def __init__(self, ref_el, m):
QuadratureRule.__init__(self, ref_el, xs, ws)


class GaussLegendreQuadratureLineRule(QuadratureRule):
class GaussLegendreQuadratureLineRule(GaussJacobiQuadratureLineRule):
"""Produce the Gauss--Legendre quadrature rules on the interval using
the implementation in numpy. This facilitates implementing
discontinuous spectral elements.
The quadrature rule uses m points for a degree of precision of 2m-1.
"""
def __init__(self, ref_el, m):
if m < 1:
raise ValueError(
"Gauss-Legendre quadrature invalid for fewer than 2 points")

xs_ref, ws_ref = numpy.polynomial.legendre.leggauss(m)

A, b = reference_element.make_affine_mapping(((-1.,), (1.)),
ref_el.get_vertices())

mapping = lambda x: numpy.dot(A, x) + b

scale = numpy.linalg.det(A)

xs = tuple([tuple(mapping(x_ref)[0]) for x_ref in xs_ref])
ws = tuple([scale * w for w in ws_ref])

QuadratureRule.__init__(self, ref_el, xs, ws)
super(GaussLegendreQuadratureLineRule, self).__init__(ref_el, m)


class RadauQuadratureLineRule(QuadratureRule):
Expand Down

0 comments on commit c70914c

Please sign in to comment.