"""Quadrature schemes on cells
This module generates quadrature schemes on reference cells that integrate
exactly a polynomial of a given degree using a specified scheme.
Scheme options are:
scheme="default"
scheme="canonical" (collapsed Gauss scheme)
Background on the schemes:
Keast rules for tetrahedra:
Keast, P. Moderate-degree tetrahedral quadrature formulas, Computer
Methods in Applied Mechanics and Engineering 55(3):339-348, 1986.
http://dx.doi.org/10.1016/0045-7825(86)90059-9
"""
# Copyright (C) 2011 Garth N. Wells
# Copyright (C) 2016 Miklos Homolya
#
# This file is part of FIAT (https://www.fenicsproject.org)
#
# SPDX-License-Identifier: LGPL-3.0-or-later
#
# First added: 2011-04-19
# Last changed: 2011-04-19
# NumPy
from numpy import array, arange, float64
# FIAT
from FIAT.reference_element import QUADRILATERAL, HEXAHEDRON, TENSORPRODUCT, UFCTriangle, UFCTetrahedron
from FIAT.quadrature import QuadratureRule, make_quadrature, make_tensor_product_quadrature
[docs]
def create_quadrature(ref_el, degree, scheme="default"):
"""
Generate quadrature rule for given reference element
that will integrate an polynomial of order 'degree' exactly.
For low-degree (<=6) polynomials on triangles and tetrahedra, this
uses hard-coded rules, otherwise it falls back to a collapsed
Gauss scheme on simplices. On tensor-product cells, it is a
tensor-product quadrature rule of the subcells.
:arg cell: The FIAT cell to create the quadrature for.
:arg degree: The degree of polynomial that the rule should
integrate exactly.
"""
if ref_el.get_shape() == TENSORPRODUCT:
try:
degree = tuple(degree)
except TypeError:
degree = (degree,) * len(ref_el.cells)
assert len(ref_el.cells) == len(degree)
quad_rules = [create_quadrature(c, d, scheme)
for c, d in zip(ref_el.cells, degree)]
return make_tensor_product_quadrature(*quad_rules)
if ref_el.get_shape() in [QUADRILATERAL, HEXAHEDRON]:
return create_quadrature(ref_el.product, degree, scheme)
if degree < 0:
raise ValueError("Need positive degree, not %d" % degree)
if scheme == "default":
# TODO: Point transformation to support good schemes on
# non-UFC reference elements.
if isinstance(ref_el, UFCTriangle):
return _triangle_scheme(degree)
elif isinstance(ref_el, UFCTetrahedron):
return _tetrahedron_scheme(degree)
else:
return _fiat_scheme(ref_el, degree)
elif scheme == "canonical":
return _fiat_scheme(ref_el, degree)
elif scheme == "KMV": # Kong-Mulder-Veldhuizen scheme
return _kmv_lump_scheme(ref_el, degree)
else:
raise ValueError("Unknown quadrature scheme: %s." % scheme)
def _fiat_scheme(ref_el, degree):
"""Get quadrature scheme from FIAT interface"""
# Number of points per axis for exact integration
num_points_per_axis = (degree + 1 + 1) // 2
# Create and return FIAT quadrature rule
return make_quadrature(ref_el, num_points_per_axis)
def _kmv_lump_scheme(ref_el, degree):
"""Specialized quadrature schemes for P < 6 for KMV simplical elements."""
sd = ref_el.get_spatial_dimension()
# set the unit element
if sd == 2:
T = UFCTriangle()
elif sd == 3:
T = UFCTetrahedron()
else:
raise ValueError("Dimension not supported")
if degree == 1:
x = ref_el.vertices
w = arange(sd + 1, dtype=float64)
if sd == 2:
w[:] = 1.0 / 6.0
elif sd == 3:
w[:] = 1.0 / 24.0
else:
raise ValueError("Dimension not supported")
elif degree == 2:
if sd == 2:
x = list(ref_el.vertices)
for e in range(3):
x.extend(ref_el.make_points(1, e, 2)) # edge midpoints
x.extend(ref_el.make_points(2, 0, 3)) # barycenter
w = arange(7, dtype=float64)
w[0:3] = 1.0 / 40.0
w[3:6] = 1.0 / 15.0
w[6] = 9.0 / 40.0
elif sd == 3:
x = list(ref_el.vertices)
x.extend(
[
(0.0, 0.50, 0.50),
(0.50, 0.0, 0.50),
(0.50, 0.50, 0.0),
(0.0, 0.0, 0.50),
(0.0, 0.50, 0.0),
(0.50, 0.0, 0.0),
]
)
# in facets
x.extend(
[
(0.33333333333333337, 0.3333333333333333, 0.3333333333333333),
(0.0, 0.3333333333333333, 0.3333333333333333),
(0.3333333333333333, 0.0, 0.3333333333333333),
(0.3333333333333333, 0.3333333333333333, 0.0),
]
)
# in the cell
x.extend([(1 / 4, 1 / 4, 1 / 4)])
w = arange(15, dtype=float64)
w[0:4] = 17.0 / 5040.0
w[4:10] = 2.0 / 315.0
w[10:14] = 9.0 / 560.0
w[14] = 16.0 / 315.0
else:
raise ValueError("Dimension not supported")
elif degree == 3:
if sd == 2:
alpha = 0.2934695559090401
beta = 0.2073451756635909
x = list(ref_el.vertices)
x.extend(
[
(1 - alpha, alpha),
(alpha, 1 - alpha),
(0.0, 1 - alpha),
(0.0, alpha),
(alpha, 0.0),
(1 - alpha, 0.0),
] # edge points
)
x.extend(
[(beta, beta), (1 - 2 * beta, beta), (beta, 1 - 2 * beta)]
) # points in center of cell
w = arange(12, dtype=float64)
w[0:3] = 0.007436456512410291
w[3:9] = 0.02442084061702551
w[9:12] = 0.1103885289202054
elif sd == 3:
x = list(ref_el.vertices)
x.extend(
[
(0, 0.685789657581967, 0.314210342418033),
(0, 0.314210342418033, 0.685789657581967),
(0.314210342418033, 0, 0.685789657581967),
(0.685789657581967, 0, 0.314210342418033),
(0.685789657581967, 0.314210342418033, 0.0),
(0.314210342418033, 0.685789657581967, 0.0),
(0, 0, 0.685789657581967),
(0, 0, 0.314210342418033),
(0, 0.314210342418033, 0.0),
(0, 0.685789657581967, 0.0),
(0.314210342418033, 0, 0.0),
(0.685789657581967, 0, 0.0),
]
) # 12 points on edges of facets (0-->1-->2)
x.extend(
[
(0.21548220313557542, 0.5690355937288492, 0.21548220313557542),
(0.21548220313557542, 0.21548220313557542, 0.5690355937288492),
(0.5690355937288492, 0.21548220313557542, 0.21548220313557542),
(0.0, 0.5690355937288492, 0.21548220313557542),
(0.0, 0.21548220313557542, 0.5690355937288492),
(0.0, 0.21548220313557542, 0.21548220313557542),
(0.5690355937288492, 0.0, 0.21548220313557542),
(0.21548220313557542, 0.0, 0.5690355937288492),
(0.21548220313557542, 0.0, 0.21548220313557542),
(0.5690355937288492, 0.21548220313557542, 0.0),
(0.21548220313557542, 0.5690355937288492, 0.0),
(0.21548220313557542, 0.21548220313557542, 0.0),
]
) # 12 points (3 points on each facet, 1st two parallel to edge 0)
alpha = 1 / 6
x.extend(
[
(alpha, alpha, 0.5),
(0.5, alpha, alpha),
(alpha, 0.5, alpha),
(alpha, alpha, alpha),
]
) # 4 points inside the cell
w = arange(32, dtype=float64)
w[0:4] = 0.00068688236002531922325120561367839
w[4:16] = 0.0015107814913526136472998739890272
w[16:28] = 0.0050062894680040258624242888174649
w[28:32] = 0.021428571428571428571428571428571
else:
raise ValueError("Dimension not supported")
elif degree == 4:
if sd == 2:
alpha = 0.2113248654051871 # 0.2113248654051871
beta1 = 0.4247639617258106 # 0.4247639617258106
beta2 = 0.130791593829745 # 0.130791593829745
x = list(ref_el.vertices)
for e in range(3):
x.extend(ref_el.make_points(1, e, 2)) # edge midpoints
x.extend(
[
(1 - alpha, alpha),
(alpha, 1 - alpha),
(0.0, 1 - alpha),
(0.0, alpha),
(alpha, 0.0),
(1 - alpha, 0.0),
] # edge points
)
x.extend(
[(beta1, beta1), (1 - 2 * beta1, beta1), (beta1, 1 - 2 * beta1)]
) # points in center of cell
x.extend(
[(beta2, beta2), (1 - 2 * beta2, beta2), (beta2, 1 - 2 * beta2)]
) # points in center of cell
w = arange(18, dtype=float64)
w[0:3] = 0.003174603174603175 # chk
w[3:6] = 0.0126984126984127 # chk 0.0126984126984127
w[6:12] = 0.01071428571428571 # chk 0.01071428571428571
w[12:15] = 0.07878121446939182 # chk 0.07878121446939182
w[15:18] = 0.05058386489568756 # chk 0.05058386489568756
else:
raise ValueError("Dimension not supported")
elif degree == 5:
if sd == 2:
alpha1 = 0.3632980741536860e-00
alpha2 = 0.1322645816327140e-00
beta1 = 0.4578368380791611e-00
beta2 = 0.2568591072619591e-00
beta3 = 0.5752768441141011e-01
gamma1 = 0.7819258362551702e-01
delta1 = 0.2210012187598900e-00
x = list(ref_el.vertices)
x.extend(
[
(1 - alpha1, alpha1),
(alpha1, 1 - alpha1),
(0.0, 1 - alpha1),
(0.0, alpha1),
(alpha1, 0.0),
(1 - alpha1, 0.0),
] # edge points
)
x.extend(
[
(1 - alpha2, alpha2),
(alpha2, 1 - alpha2),
(0.0, 1 - alpha2),
(0.0, alpha2),
(alpha2, 0.0),
(1 - alpha2, 0.0),
] # edge points
)
x.extend(
[(beta1, beta1), (1 - 2 * beta1, beta1), (beta1, 1 - 2 * beta1)]
) # points in center of cell
x.extend(
[(beta2, beta2), (1 - 2 * beta2, beta2), (beta2, 1 - 2 * beta2)]
) # points in center of cell
x.extend(
[(beta3, beta3), (1 - 2 * beta3, beta3), (beta3, 1 - 2 * beta3)]
) # points in center of cell
x.extend(
[
(gamma1, delta1),
(1 - gamma1 - delta1, delta1),
(gamma1, 1 - gamma1 - delta1),
(delta1, gamma1),
(1 - gamma1 - delta1, gamma1),
(delta1, 1 - gamma1 - delta1),
] # edge points
)
w = arange(30, dtype=float64)
w[0:3] = 0.7094239706792450e-03
w[3:9] = 0.6190565003676629e-02
w[9:15] = 0.3480578640489211e-02
w[15:18] = 0.3453043037728279e-01
w[18:21] = 0.4590123763076286e-01
w[21:24] = 0.1162613545961757e-01
w[24:30] = 0.2727857596999626e-01
else:
raise ValueError("Dimension not supported")
# Return scheme
return QuadratureRule(T, x, w)
def _triangle_scheme(degree):
"""Return a quadrature scheme on a triangle of specified order. Falls
back on canonical rule for higher orders."""
if degree == 0 or degree == 1:
# Scheme from Zienkiewicz and Taylor, 1 point, degree of precision 1
x = array([[1.0/3.0, 1.0/3.0]])
w = array([0.5])
elif degree == 2:
# Scheme from Strang and Fix, 3 points, degree of precision 2
x = array([[1.0/6.0, 1.0/6.0],
[1.0/6.0, 2.0/3.0],
[2.0/3.0, 1.0/6.0]])
w = arange(3, dtype=float64)
w[:] = 1.0/6.0
elif degree == 3:
# Scheme from Strang and Fix, 6 points, degree of precision 3
x = array([[0.659027622374092, 0.231933368553031],
[0.659027622374092, 0.109039009072877],
[0.231933368553031, 0.659027622374092],
[0.231933368553031, 0.109039009072877],
[0.109039009072877, 0.659027622374092],
[0.109039009072877, 0.231933368553031]])
w = arange(6, dtype=float64)
w[:] = 1.0/12.0
elif degree == 4:
# Scheme from Strang and Fix, 6 points, degree of precision 4
x = array([[0.816847572980459, 0.091576213509771],
[0.091576213509771, 0.816847572980459],
[0.091576213509771, 0.091576213509771],
[0.108103018168070, 0.445948490915965],
[0.445948490915965, 0.108103018168070],
[0.445948490915965, 0.445948490915965]])
w = arange(6, dtype=float64)
w[0:3] = 0.109951743655322
w[3:6] = 0.223381589678011
w = w/2.0
elif degree == 5:
# Scheme from Strang and Fix, 7 points, degree of precision 5
x = array([[0.33333333333333333, 0.33333333333333333],
[0.79742698535308720, 0.10128650732345633],
[0.10128650732345633, 0.79742698535308720],
[0.10128650732345633, 0.10128650732345633],
[0.05971587178976981, 0.47014206410511505],
[0.47014206410511505, 0.05971587178976981],
[0.47014206410511505, 0.47014206410511505]])
w = arange(7, dtype=float64)
w[0] = 0.22500000000000000
w[1:4] = 0.12593918054482717
w[4:7] = 0.13239415278850616
w = w/2.0
elif degree == 6:
# Scheme from Strang and Fix, 12 points, degree of precision 6
x = array([[0.873821971016996, 0.063089014491502],
[0.063089014491502, 0.873821971016996],
[0.063089014491502, 0.063089014491502],
[0.501426509658179, 0.249286745170910],
[0.249286745170910, 0.501426509658179],
[0.249286745170910, 0.249286745170910],
[0.636502499121399, 0.310352451033785],
[0.636502499121399, 0.053145049844816],
[0.310352451033785, 0.636502499121399],
[0.310352451033785, 0.053145049844816],
[0.053145049844816, 0.636502499121399],
[0.053145049844816, 0.310352451033785]])
w = arange(12, dtype=float64)
w[0:3] = 0.050844906370207
w[3:6] = 0.116786275726379
w[6:12] = 0.082851075618374
w = w/2.0
else:
# Get canonical scheme
return _fiat_scheme(UFCTriangle(), degree)
# Return scheme
return QuadratureRule(UFCTriangle(), x, w)
def _tetrahedron_scheme(degree):
"""Return a quadrature scheme on a tetrahedron of specified
degree. Falls back on canonical rule for higher orders"""
if degree == 0 or degree == 1:
# Scheme from Zienkiewicz and Taylor, 1 point, degree of precision 1
x = array([[1.0/4.0, 1.0/4.0, 1.0/4.0]])
w = array([1.0/6.0])
elif degree == 2:
# Scheme from Zienkiewicz and Taylor, 4 points, degree of precision 2
a, b = 0.585410196624969, 0.138196601125011
x = array([[a, b, b],
[b, a, b],
[b, b, a],
[b, b, b]])
w = arange(4, dtype=float64)
w[:] = 1.0/24.0
elif degree == 3:
# Scheme from Zienkiewicz and Taylor, 5 points, degree of precision 3
# Note: this scheme has a negative weight
x = array([[0.2500000000000000, 0.2500000000000000, 0.2500000000000000],
[0.5000000000000000, 0.1666666666666666, 0.1666666666666666],
[0.1666666666666666, 0.5000000000000000, 0.1666666666666666],
[0.1666666666666666, 0.1666666666666666, 0.5000000000000000],
[0.1666666666666666, 0.1666666666666666, 0.1666666666666666]])
w = arange(5, dtype=float64)
w[0] = -0.8
w[1:5] = 0.45
w = w/6.0
elif degree == 4:
# Keast rule, 14 points, degree of precision 4
# Values taken from http://people.sc.fsu.edu/~jburkardt/datasets/quadrature_rules_tet/quadrature_rules_tet.html
# (KEAST5)
x = array([[0.0000000000000000, 0.5000000000000000, 0.5000000000000000],
[0.5000000000000000, 0.0000000000000000, 0.5000000000000000],
[0.5000000000000000, 0.5000000000000000, 0.0000000000000000],
[0.5000000000000000, 0.0000000000000000, 0.0000000000000000],
[0.0000000000000000, 0.5000000000000000, 0.0000000000000000],
[0.0000000000000000, 0.0000000000000000, 0.5000000000000000],
[0.6984197043243866, 0.1005267652252045, 0.1005267652252045],
[0.1005267652252045, 0.1005267652252045, 0.1005267652252045],
[0.1005267652252045, 0.1005267652252045, 0.6984197043243866],
[0.1005267652252045, 0.6984197043243866, 0.1005267652252045],
[0.0568813795204234, 0.3143728734931922, 0.3143728734931922],
[0.3143728734931922, 0.3143728734931922, 0.3143728734931922],
[0.3143728734931922, 0.3143728734931922, 0.0568813795204234],
[0.3143728734931922, 0.0568813795204234, 0.3143728734931922]])
w = arange(14, dtype=float64)
w[0:6] = 0.0190476190476190
w[6:10] = 0.0885898247429807
w[10:14] = 0.1328387466855907
w = w/6.0
elif degree == 5:
# Keast rule, 15 points, degree of precision 5
# Values taken from http://people.sc.fsu.edu/~jburkardt/datasets/quadrature_rules_tet/quadrature_rules_tet.html
# (KEAST6)
x = array([[0.2500000000000000, 0.2500000000000000, 0.2500000000000000],
[0.0000000000000000, 0.3333333333333333, 0.3333333333333333],
[0.3333333333333333, 0.3333333333333333, 0.3333333333333333],
[0.3333333333333333, 0.3333333333333333, 0.0000000000000000],
[0.3333333333333333, 0.0000000000000000, 0.3333333333333333],
[0.7272727272727273, 0.0909090909090909, 0.0909090909090909],
[0.0909090909090909, 0.0909090909090909, 0.0909090909090909],
[0.0909090909090909, 0.0909090909090909, 0.7272727272727273],
[0.0909090909090909, 0.7272727272727273, 0.0909090909090909],
[0.4334498464263357, 0.0665501535736643, 0.0665501535736643],
[0.0665501535736643, 0.4334498464263357, 0.0665501535736643],
[0.0665501535736643, 0.0665501535736643, 0.4334498464263357],
[0.0665501535736643, 0.4334498464263357, 0.4334498464263357],
[0.4334498464263357, 0.0665501535736643, 0.4334498464263357],
[0.4334498464263357, 0.4334498464263357, 0.0665501535736643]])
w = arange(15, dtype=float64)
w[0] = 0.1817020685825351
w[1:5] = 0.0361607142857143
w[5:9] = 0.0698714945161738
w[9:15] = 0.0656948493683187
w = w/6.0
elif degree == 6:
# Keast rule, 24 points, degree of precision 6
# Values taken from http://people.sc.fsu.edu/~jburkardt/datasets/quadrature_rules_tet/quadrature_rules_tet.html
# (KEAST7)
x = array([[0.3561913862225449, 0.2146028712591517, 0.2146028712591517],
[0.2146028712591517, 0.2146028712591517, 0.2146028712591517],
[0.2146028712591517, 0.2146028712591517, 0.3561913862225449],
[0.2146028712591517, 0.3561913862225449, 0.2146028712591517],
[0.8779781243961660, 0.0406739585346113, 0.0406739585346113],
[0.0406739585346113, 0.0406739585346113, 0.0406739585346113],
[0.0406739585346113, 0.0406739585346113, 0.8779781243961660],
[0.0406739585346113, 0.8779781243961660, 0.0406739585346113],
[0.0329863295731731, 0.3223378901422757, 0.3223378901422757],
[0.3223378901422757, 0.3223378901422757, 0.3223378901422757],
[0.3223378901422757, 0.3223378901422757, 0.0329863295731731],
[0.3223378901422757, 0.0329863295731731, 0.3223378901422757],
[0.2696723314583159, 0.0636610018750175, 0.0636610018750175],
[0.0636610018750175, 0.2696723314583159, 0.0636610018750175],
[0.0636610018750175, 0.0636610018750175, 0.2696723314583159],
[0.6030056647916491, 0.0636610018750175, 0.0636610018750175],
[0.0636610018750175, 0.6030056647916491, 0.0636610018750175],
[0.0636610018750175, 0.0636610018750175, 0.6030056647916491],
[0.0636610018750175, 0.2696723314583159, 0.6030056647916491],
[0.2696723314583159, 0.6030056647916491, 0.0636610018750175],
[0.6030056647916491, 0.0636610018750175, 0.2696723314583159],
[0.0636610018750175, 0.6030056647916491, 0.2696723314583159],
[0.2696723314583159, 0.0636610018750175, 0.6030056647916491],
[0.6030056647916491, 0.2696723314583159, 0.0636610018750175]])
w = arange(24, dtype=float64)
w[0:4] = 0.0399227502581679
w[4:8] = 0.0100772110553207
w[8:12] = 0.0553571815436544
w[12:24] = 0.0482142857142857
w = w/6.0
else:
# Get canonical scheme
return _fiat_scheme(UFCTetrahedron(), degree)
# Return scheme
return QuadratureRule(UFCTetrahedron(), x, w)