# Copyright (C) 2008 Robert C. Kirby (Texas Tech University)
#
# This file is part of FIAT (https://www.fenicsproject.org)
#
# SPDX-License-Identifier: LGPL-3.0-or-later
#
# Modified by Marie E. Rognes (meg@simula.no), 2012
# Modified by David A. Ham (david.ham@imperial.ac.uk), 2015
import itertools
import math
import numpy
from FIAT import reference_element, expansions, jacobi, orthopoly
[docs]
class QuadratureRule(object):
"""General class that models integration over a reference element
as the weighted sum of a function evaluated at a set of points."""
def __init__(self, ref_el, pts, wts):
if len(wts) != len(pts):
raise ValueError("Have %d weights, but %d points" % (len(wts), len(pts)))
self.ref_el = ref_el
self.pts = pts
self.wts = wts
[docs]
def get_points(self):
return numpy.array(self.pts)
[docs]
def get_weights(self):
return numpy.array(self.wts)
[docs]
def integrate(self, f):
return sum([w * f(x) for (x, w) in zip(self.pts, self.wts)])
[docs]
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):
# this gives roots on the default (-1,1) reference element
# (xs_ref, ws_ref) = compute_gauss_jacobi_rule(a, b, m)
(xs_ref, ws_ref) = compute_gauss_jacobi_rule(0., 0., m)
Ref1 = reference_element.DefaultLine()
A, b = reference_element.make_affine_mapping(Ref1.get_vertices(),
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)
[docs]
class GaussLobattoLegendreQuadratureLineRule(QuadratureRule):
"""Implement the Gauss-Lobatto-Legendre quadrature rules on the interval using
Greg von Winckel's implementation. This facilitates implementing
spectral elements.
The quadrature rule uses m points for a degree of precision of 2m-3.
"""
def __init__(self, ref_el, m):
if m < 2:
raise ValueError(
"Gauss-Labotto-Legendre quadrature invalid for fewer than 2 points")
Ref1 = reference_element.DefaultLine()
verts = Ref1.get_vertices()
if m > 2:
# Calculate the recursion coefficients.
alpha, beta = orthopoly.rec_jacobi(m, 0, 0)
xs_ref, ws_ref = orthopoly.lobatto(alpha, beta, verts[0][0], verts[1][0])
else:
# Special case for lowest order.
xs_ref = [v[0] for v in verts[:]]
ws_ref = (0.5 * (xs_ref[1] - xs_ref[0]), ) * 2
A, b = reference_element.make_affine_mapping(Ref1.get_vertices(),
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)
[docs]
class GaussLegendreQuadratureLineRule(QuadratureRule):
"""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)
[docs]
class RadauQuadratureLineRule(QuadratureRule):
"""Produce the Gauss--Radau quadrature rules on the interval using
an adaptation of Winkel's Matlab code.
The quadrature rule uses m points for a degree of precision of 2m-1.
"""
def __init__(self, ref_el, m, right=True):
assert m >= 1
N = m - 1
# Use Chebyshev-Gauss-Radau nodes as initial guess for LGR nodes
x = -numpy.cos(2 * numpy.pi * numpy.linspace(0, N, m) / (2 * N + 1))
P = numpy.zeros((N + 1, N + 2))
xold = 2
free = numpy.arange(1, N + 1, dtype='int')
while numpy.max(numpy.abs(x - xold)) > 5e-16:
xold = x.copy()
P[0, :] = (-1) ** numpy.arange(0, N + 2)
P[free, 0] = 1
P[free, 1] = x[free]
for k in range(2, N + 2):
P[free, k] = ((2 * k - 1) * x[free] * P[free, k - 1] - (k - 1) * P[free, k - 2]) / k
x[free] = xold[free] - ((1 - xold[free]) / (N + 1)) * (P[free, N] + P[free, N + 1]) / (P[free, N] - P[free, N + 1])
# The Legendre-Gauss-Radau Vandermonde
P = P[:, :-1]
# Compute the weights
w = numpy.zeros(N + 1)
w[0] = 2 / (N + 1) ** 2
w[free] = (1 - x[free])/((N + 1) * P[free, -1])**2
if right:
x = numpy.flip(-x)
w = numpy.flip(w)
xs_ref = x
ws_ref = w
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)
[docs]
class CollapsedQuadratureTriangleRule(QuadratureRule):
"""Implements the collapsed quadrature rules defined in
Karniadakis & Sherwin by mapping products of Gauss-Jacobi rules
from the square to the triangle."""
def __init__(self, ref_el, m):
ptx, wx = compute_gauss_jacobi_rule(0., 0., m)
pty, wy = compute_gauss_jacobi_rule(1., 0., m)
# map ptx , pty
pts_ref = [expansions.xi_triangle((x, y))
for x in ptx for y in pty]
Ref1 = reference_element.DefaultTriangle()
A, b = reference_element.make_affine_mapping(Ref1.get_vertices(),
ref_el.get_vertices())
mapping = lambda x: numpy.dot(A, x) + b
scale = numpy.linalg.det(A)
pts = tuple([tuple(mapping(x)) for x in pts_ref])
wts = [0.5 * scale * w1 * w2 for w1 in wx for w2 in wy]
QuadratureRule.__init__(self, ref_el, tuple(pts), tuple(wts))
[docs]
class CollapsedQuadratureTetrahedronRule(QuadratureRule):
"""Implements the collapsed quadrature rules defined in
Karniadakis & Sherwin by mapping products of Gauss-Jacobi rules
from the cube to the tetrahedron."""
def __init__(self, ref_el, m):
ptx, wx = compute_gauss_jacobi_rule(0., 0., m)
pty, wy = compute_gauss_jacobi_rule(1., 0., m)
ptz, wz = compute_gauss_jacobi_rule(2., 0., m)
# map ptx , pty
pts_ref = [expansions.xi_tetrahedron((x, y, z))
for x in ptx for y in pty for z in ptz]
Ref1 = reference_element.DefaultTetrahedron()
A, b = reference_element.make_affine_mapping(Ref1.get_vertices(),
ref_el.get_vertices())
mapping = lambda x: numpy.dot(A, x) + b
scale = numpy.linalg.det(A)
pts = tuple([tuple(mapping(x)) for x in pts_ref])
wts = [scale * 0.125 * w1 * w2 * w3
for w1 in wx for w2 in wy for w3 in wz]
QuadratureRule.__init__(self, ref_el, tuple(pts), tuple(wts))
[docs]
class UFCTetrahedronFaceQuadratureRule(QuadratureRule):
"""Highly specialized quadrature rule for the face of a
tetrahedron, mapped from a reference triangle, used for higher
order Nedelecs"""
def __init__(self, face_number, degree):
# Create quadrature rule on reference triangle
reference_triangle = reference_element.UFCTriangle()
reference_rule = make_quadrature(reference_triangle, degree)
ref_points = reference_rule.get_points()
ref_weights = reference_rule.get_weights()
# Get geometry information about the face of interest
reference_tet = reference_element.UFCTetrahedron()
face = reference_tet.get_topology()[2][face_number]
vertices = reference_tet.get_vertices_of_subcomplex(face)
# Use tet to map points and weights on the appropriate face
vertices = [numpy.array(list(vertex)) for vertex in vertices]
x0 = vertices[0]
J = numpy.vstack([vertices[1] - x0, vertices[2] - x0]).T
# This is just a very numpyfied way of writing J*p + x0:
points = numpy.einsum("ij,kj->ki", J, ref_points) + x0
# Map weights: multiply reference weights by sqrt(|J^T J|)
detJTJ = numpy.linalg.det(numpy.dot(J.T, J))
weights = numpy.sqrt(detJTJ) * ref_weights
# Initialize super class with new points and weights
QuadratureRule.__init__(self, reference_tet, points, weights)
self._reference_rule = reference_rule
self._J = J
[docs]
def reference_rule(self):
return self._reference_rule
[docs]
def jacobian(self):
return self._J
[docs]
def make_quadrature(ref_el, m):
"""Returns the collapsed quadrature rule using m points per
direction on the given reference element. In the tensor product
case, m is a tuple."""
if isinstance(m, tuple):
min_m = min(m)
else:
min_m = m
msg = "Expecting at least one (not %d) quadrature point per direction" % min_m
assert (min_m > 0), msg
if ref_el.get_shape() == reference_element.POINT:
return QuadratureRule(ref_el, [()], [1])
elif ref_el.get_shape() == reference_element.LINE:
return GaussJacobiQuadratureLineRule(ref_el, m)
elif ref_el.get_shape() == reference_element.TRIANGLE:
return CollapsedQuadratureTriangleRule(ref_el, m)
elif ref_el.get_shape() == reference_element.TETRAHEDRON:
return CollapsedQuadratureTetrahedronRule(ref_el, m)
elif ref_el.get_shape() == reference_element.QUADRILATERAL:
line_rule = GaussJacobiQuadratureLineRule(ref_el.construct_subelement(1), m)
return make_tensor_product_quadrature(line_rule, line_rule)
elif ref_el.get_shape() == reference_element.HEXAHEDRON:
line_rule = GaussJacobiQuadratureLineRule(ref_el.construct_subelement(1), m)
return make_tensor_product_quadrature(line_rule, line_rule, line_rule)
else:
raise ValueError("Unable to make quadrature for cell: %s" % ref_el)
[docs]
def make_tensor_product_quadrature(*quad_rules):
"""Returns the quadrature rule for a TensorProduct cell, by combining
the quadrature rules of the components."""
ref_el = reference_element.TensorProductCell(*[q.ref_el
for q in quad_rules])
# Coordinates are "concatenated", weights are multiplied
pts = [list(itertools.chain(*pt_tuple))
for pt_tuple in itertools.product(*[q.pts for q in quad_rules])]
wts = [numpy.prod(wt_tuple)
for wt_tuple in itertools.product(*[q.wts for q in quad_rules])]
return QuadratureRule(ref_el, pts, wts)
# rule to get Gauss-Jacobi points
[docs]
def compute_gauss_jacobi_points(a, b, m):
"""Computes the m roots of P_{m}^{a,b} on [-1,1] by Newton's method.
The initial guesses are the Chebyshev points. Algorithm
implemented in Python from the pseudocode given by Karniadakis and
Sherwin"""
x = []
eps = 1.e-8
max_iter = 100
for k in range(0, m):
r = -math.cos((2.0 * k + 1.0) * math.pi / (2.0 * m))
if k > 0:
r = 0.5 * (r + x[k - 1])
j = 0
delta = 2 * eps
while j < max_iter:
s = 0
for i in range(0, k):
s = s + 1.0 / (r - x[i])
f = jacobi.eval_jacobi(a, b, m, r)
fp = jacobi.eval_jacobi_deriv(a, b, m, r)
delta = f / (fp - f * s)
r = r - delta
if math.fabs(delta) < eps:
break
else:
j = j + 1
x.append(r)
return x
[docs]
def compute_gauss_jacobi_rule(a, b, m):
xs = compute_gauss_jacobi_points(a, b, m)
a1 = math.pow(2, a + b + 1)
a2 = math.gamma(a + m + 1)
a3 = math.gamma(b + m + 1)
a4 = math.gamma(a + b + m + 1)
a5 = math.factorial(m)
a6 = a1 * a2 * a3 / a4 / a5
ws = [a6 / (1.0 - x**2.0) / jacobi.eval_jacobi_deriv(a, b, m, x)**2.0
for x in xs]
return xs, ws