import numpy as nm
try:
import sympy as sm
except ImportError:
sm = None
import pytest
from sfepy.base.base import assert_, ordered_iteritems
import sfepy.base.testing as tst
[docs]def symarray(prefix, shape):
"""
Copied from SymPy so that the tests pass for its different versions.
"""
arr = nm.empty(shape, dtype=object)
for index in nm.ndindex(shape):
arr[index] = sm.Symbol('%s_%s' % (prefix, '_'.join(map(str, index))))
return arr
[docs]def get_poly(order, dim, is_simplex=False):
"""
Construct a polynomial of given `order` in space dimension `dim`,
and integrate it symbolically over a rectangular or simplex domain
for coordinates in [0, 1].
"""
xs = symarray('x', dim)
opd = max(1, int((order + 1) / dim))
poly = 1.0
oo = 0
for ii, x in enumerate(xs):
if ((oo + opd) > order) or (ii == (len(xs) - 1)):
opd = max(order - oo, 0)
poly *= (x**opd + 1)
oo += opd
assert_(oo == order)
limits = [[xs[ii], 0, 1] for ii in range(dim)]
if is_simplex:
for ii in range(1, dim):
for ip in range(0, ii):
limits[ii][2] -= xs[ip]
integral = sm.integrate(poly, *reversed(limits))
return xs, poly, limits, integral
[docs]def test_weight_consistency():
"""
Test if integral of 1 (= sum of weights) gives the domain volume.
"""
from sfepy.discrete.quadratures import quadrature_tables
ok = True
for geometry, qps in ordered_iteritems(quadrature_tables):
tst.report('geometry:', geometry)
for order, qp in ordered_iteritems(qps):
diff = nm.abs(qp.weights.sum() - qp.volume)
_ok = diff < 1e-14
tst.report('order %d: %s (%.2e)' % (order, _ok, diff))
ok = ok and _ok
assert_(ok)
def _test_quadratures(quad, geometry, order):
dim = int(geometry[0])
n_v = int(geometry[2])
is_simplex = n_v == (dim + 1)
porder = order if is_simplex else dim * order
xs, poly, limits, integral = get_poly(porder, dim,
is_simplex=is_simplex)
tst.report(' polynomial:', poly)
tst.report(' limits:', limits)
tst.report(' integral:', integral)
def fun(coors):
vals = nm.empty(coors.shape[0], dtype=nm.float64)
subs = {}
for ir, cc in enumerate(coors):
for ic, x in enumerate(xs):
subs[x] = coors[ir,ic]
vals[ir] = float(poly.subs(subs))
return vals
val = quad.integrate(fun, order=order, geometry=geometry)
ok = nm.allclose(val, float(integral), rtol=0.0, atol=1e-11)
tst.report(' sym. == num.: %f == %f -> %s' %
(float(integral), val, ok))
return ok, integral, val
[docs]@pytest.mark.slow
def test_quadratures():
"""
Test if the quadratures have orders they claim to have, using
symbolic integration by sympy.
"""
from sfepy.discrete.quadratures import quadrature_tables, tp_geometries
from sfepy.discrete.integrals import Integral
if sm is None:
tst.report('cannot import sympy, skipping')
return
quad = Integral('test_integral')
ok = True
failed = []
for geometry, qps in ordered_iteritems(quadrature_tables):
tst.report('geometry:', geometry)
if geometry == '0_1':
continue
elif geometry in tp_geometries:
iter_qp = range(1, 11)
elif geometry == '2_3':
iter_qp = range(1, 25)
elif geometry == '3_4':
iter_qp = range(1, 12)
elif geometry == '3_6':
iter_qp = range(4, 9)
else:
iter_qp = sorted(qps.keys())
for order in iter_qp:
tst.report('order:', order)
_ok, integral, val = _test_quadratures(quad, geometry, order)
if not _ok:
failed.append((geometry, order, float(integral), val))
ok = ok and _ok
if not ok:
tst.report('failed:')
for aux in failed:
tst.report(aux, '%.1e' % abs(aux[2] - aux[3]))
assert_(ok)