Skip to content

Commit

Permalink
remove Duffy, tidy up
Browse files Browse the repository at this point in the history
  • Loading branch information
@pbrubeck
pbrubeck committed Oct 31, 2023
1 parent bac180c commit d213473
Showing 1 changed file with 95 additions and 234 deletions.
329 changes: 95 additions & 234 deletions FIAT/expansions.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -32,201 +32,10 @@ def jrc(a, b, n):
return an, bn, cn


def _tabulate_dpts(tabulator, D, n, order, pts):
X = sympy.DeferredVector('x')

def form_derivative(F):
'''Forms the derivative recursively, i.e.,
F -> [F_x, F_y, F_z],
[F_x, F_y, F_z] -> [[F_xx, F_xy, F_xz],
[F_yx, F_yy, F_yz],
[F_zx, F_zy, F_zz]]
and so forth.
'''
out = []
try:
out = [sympy.diff(F, X[j]) for j in range(D)]
except (AttributeError, ValueError):
# Intercept errors like
# AttributeError: 'list' object has no attribute
# 'free_symbols'
for f in F:
out.append(form_derivative(f))
return out

def numpy_lambdify(X, F):
'''Unfortunately, SymPy's own lambdify() doesn't work well with
NumPy in that simple functions like
lambda x: 1.0,
when evaluated with NumPy arrays, return just "1.0" instead of
an array of 1s with the same shape as x. This function does that.
'''
try:
lambda_x = [numpy_lambdify(X, f) for f in F]
except TypeError: # 'function' object is not iterable
# SymPy's lambdify also works on functions that return arrays.
# However, use it componentwise here so we can add 0*x to each
# component individually. This is necessary to maintain shapes
# if evaluated with NumPy arrays.
lmbd_tmp = sympy.lambdify(X, F)
lambda_x = lambda x: lmbd_tmp(x) + 0 * x[0]
return lambda_x

def evaluate_lambda(lmbd, x):
'''Properly evaluate lambda expressions recursively for iterables.
'''
try:
values = [evaluate_lambda(l, x) for l in lmbd]
except TypeError: # 'function' object is not iterable
values = lmbd(x)
return values

# Tabulate symbolically
symbolic_tab = tabulator(n, X)
# Make sure that the entries of symbolic_tab are lists so we can
# append derivatives
symbolic_tab = [[phi] for phi in symbolic_tab]
#
data = (order + 1) * [None]
for r in range(order + 1):
shape = [len(symbolic_tab), len(pts)] + r * [D]
data[r] = numpy.empty(shape)
for i, phi in enumerate(symbolic_tab):
# Evaluate the function numerically using lambda expressions
deriv_lambda = numpy_lambdify(X, phi[r])
data[r][i] = \
numpy.array(evaluate_lambda(deriv_lambda, pts.T)).T
# Symbolically compute the next derivative.
# This actually happens once too many here; never mind for
# now.
phi.append(form_derivative(phi[-1]))
return data


def xi_triangle(eta):
"""Maps from [-1,1]^2 to the (-1,1) reference triangle."""
eta1, eta2 = eta
xi1 = 0.5 * (1.0 + eta1) * (1.0 - eta2) - 1.0
xi2 = eta2
return (xi1, xi2)


def xi_tetrahedron(eta):
"""Maps from [-1,1]^3 to the -1/1 reference tetrahedron."""
eta1, eta2, eta3 = eta
xi1 = 0.25 * (1. + eta1) * (1. - eta2) * (1. - eta3) - 1.
xi2 = 0.5 * (1. + eta2) * (1. - eta3) - 1.
xi3 = eta3
return xi1, xi2, xi3


def eta_square(xi):
"""Maps from the (-1,1) reference triangle to [-1,1]^2."""
xi1, xi2 = xi
with numpy.errstate(divide='ignore', invalid='ignore'):
eta1 = 2. * (1. + xi1) / (1. - xi2) - 1.
eta2 = xi2
if eta1.dtype != object:
eta1[numpy.logical_not(numpy.isfinite(eta1))] = 1.
return eta1, eta2


def eta_cube(xi):
"""Maps from the (-1,1) reference tetrahedron to [-1,1]^3."""
xi1, xi2, xi3 = xi
with numpy.errstate(divide='ignore', invalid='ignore'):
eta1 = 2. * (1. + xi1) / (-xi2 - xi3) - 1.
eta2 = 2. * (1. + xi2) / (1. - xi3) - 1.
eta3 = xi3
if eta1.dtype != object:
eta1[numpy.logical_not(numpy.isfinite(eta1))] = 1.
if eta2.dtype != object:
eta2[numpy.logical_not(numpy.isfinite(eta2))] = 1.
return eta1, eta2, eta3


def dubiner_1d(order, dim, x):
"""Returns a tabulation of the orthonormal Dubiner 1D polymoials, defined as
c_i P_i^(0,0)(x) if dim == 0,
c_ij (1-x)^i P_j^(2i + dim, 0)(x) otherwise."""
if dim == 0:
scale = numpy.sqrt(0.5 + numpy.arange(order + 1))
results = jacobi.eval_jacobi_batch(0, 0, order, x[:, None])
return numpy.multiply(scale[:, None], results, out=results)
sd = (order + 1) * (order + 2) // 2
phi = numpy.zeros((sd, x.size), dtype=x.dtype)
x1 = (1. - x) * 0.5
for i in range(order + 1):
n = order - i
alpha = 2 * i + dim
results = jacobi.eval_jacobi_batch(alpha, 0, n, x[:, None])
if i > 0:
results *= x1 ** i
scale = numpy.sqrt(0.5*(alpha + 1) + numpy.arange(n + 1))
numpy.multiply(scale[:, None], results, out=results)
indices = [morton_index2(i, j) for j in range(n + 1)]
phi[indices, :] = results
return phi


def dubiner_deriv_1d(order, dim, x):
"""Returns a tabulation of the first derivatives of the orthonormal Dubiner
1D polynomials."""
def recurrence(dim, n, factors, phi, dfactors=None, dphi=None):
if dim == 0:
scale = numpy.sqrt(0.5 + numpy.arange(order + 1))
results = jacobi.eval_jacobi_deriv_batch(0, 0, order, x[:, None])
return numpy.multiply(scale[:, None], results, out=results)
sd = (order + 1) * (order + 2) // 2
dphi = numpy.zeros((sd, x.size), dtype=x.dtype)
x1 = (1. - x) * 0.5
for i in range(order + 1):
n = order - i
alpha = 2 * i + dim
derivs = jacobi.eval_jacobi_deriv_batch(alpha, 0, n, x[:, None])
if i > 0:
results = jacobi.eval_jacobi_batch(alpha, 0, n, x[:, None])
derivs *= x1
derivs += results * (-0.5 * i)
if i > 1:
derivs *= x1 ** (i - 1)
scale = numpy.sqrt(0.5*(alpha + 1) + numpy.arange(n + 1))
numpy.multiply(scale[:, None], derivs, out=derivs)
indices = [morton_index2(i, j) for j in range(n + 1)]
dphi[indices, :] = derivs
return dphi


def duffy_chain_rule(A, eta, tabulations):
"""Applies the chain rule associated with an affine transformation onto the
default simplex on (-1, 1), followed by the Duffy transformation onto [-1, 1]^d.
A: the Jacobian of the affine transformation
eta: the points in [-1, 1]^d
tabulations: On entry, the tabulations of the reference gradient with
respect to eta. On exit, the tabulations of the gradient in physical space.
"""
dim = len(eta)
dphi_dxi = [tabulations[alpha] for alpha in sorted(tabulations, reverse=True) if sum(alpha) == 1]
if len(dphi_dxi) < dim:
return
eta1 = [(1. - x) * 0.5 for x in eta]
for i in range(dim):
for j in range(i):
Jij = -0.5 * (1. + eta[j])
for k in range(j + 1, dim):
if k != i:
Jij *= eta1[k]
dphi_dxi[i] -= dphi_dxi[j] * Jij
for j in range(i + 1, dim):
dphi_dxi[i] /= eta1[j]
j = 0
for alpha in sorted(tabulations, reverse=True):
if sum(alpha) == 1:
tabulations[alpha] = sum(dphi_dxi[i] * A[i][j] for i in range(dim))
j += 1


def recurrence(dim, n, factors, phi, dfactors=None, dphi=None):
if dim == 1:
elif dim == 1:
idx = lambda p: p
elif dim == 2:
idx = morton_index2
Expand All @@ -244,7 +53,7 @@ def recurrence(dim, n, factors, phi, dfactors=None, dphi=None):

# p = 1
phi[idx(1)] = f1
if dphi is not None:
if not skip_derivs:
dphi[idx(1)] = df1

# general p by recurrence
Expand All @@ -257,8 +66,8 @@ def recurrence(dim, n, factors, phi, dfactors=None, dphi=None):
phi[inext] = a * f1 * phi[icur] - b * f2 * phi[iprev]
if skip_derivs:
continue
dphi[inext] = a * f1 * dphi[icur] - b * f2 * dphi[iprev] \
+ a * phi[icur] * df1 - b * phi[iprev] * df2
dphi[inext] = (a * f1 * dphi[icur] - b * f2 * dphi[iprev] +
a * phi[icur] * df1 - b * phi[iprev] * df2)
if dim < 2:
return

Expand All @@ -281,7 +90,7 @@ def recurrence(dim, n, factors, phi, dfactors=None, dphi=None):
iprev = idx(p, q - 1)
aq, bq, cq = jrc(2 * p + 1, 0, q)
g = aq * f3 + (bq - aq) * f4
h = cq * f5
h = cq * f5
phi[inext] = g * phi[icur] - h * phi[iprev]
if skip_derivs:
continue
Expand Down Expand Up @@ -322,6 +131,94 @@ def recurrence(dim, n, factors, phi, dfactors=None, dphi=None):
dphi[inext] = (ar * z + br) * dphi[icur] + ar * phi[icur] * dz - cr * dphi[iprev]


def _tabulate_dpts(tabulator, D, n, order, pts):
X = sympy.DeferredVector('x')

def form_derivative(F):
'''Forms the derivative recursively, i.e.,
F -> [F_x, F_y, F_z],
[F_x, F_y, F_z] -> [[F_xx, F_xy, F_xz],
[F_yx, F_yy, F_yz],
[F_zx, F_zy, F_zz]]
and so forth.
'''
out = []
try:
out = [sympy.diff(F, X[j]) for j in range(D)]
except (AttributeError, ValueError):
# Intercept errors like
# AttributeError: 'list' object has no attribute
# 'free_symbols'
for f in F:
out.append(form_derivative(f))
return out

def numpy_lambdify(X, F):
'''Unfortunately, SymPy's own lambdify() doesn't work well with
NumPy in that simple functions like
lambda x: 1.0,
when evaluated with NumPy arrays, return just "1.0" instead of
an array of 1s with the same shape as x. This function does that.
'''
try:
lambda_x = [numpy_lambdify(X, f) for f in F]
except TypeError: # 'function' object is not iterable
# SymPy's lambdify also works on functions that return arrays.
# However, use it componentwise here so we can add 0*x to each
# component individually. This is necessary to maintain shapes
# if evaluated with NumPy arrays.
lmbd_tmp = sympy.lambdify(X, F)
lambda_x = lambda x: lmbd_tmp(x) + 0 * x[0]
return lambda_x

def evaluate_lambda(lmbd, x):
'''Properly evaluate lambda expressions recursively for iterables.
'''
try:
values = [evaluate_lambda(l, x) for l in lmbd]
except TypeError: # 'function' object is not iterable
values = lmbd(x)
return values

# Tabulate symbolically
symbolic_tab = tabulator(n, X)
# Make sure that the entries of symbolic_tab are lists so we can
# append derivatives
symbolic_tab = [[phi] for phi in symbolic_tab]
#
data = (order + 1) * [None]
for r in range(order + 1):
shape = [len(symbolic_tab), len(pts)] + r * [D]
data[r] = numpy.empty(shape)
for i, phi in enumerate(symbolic_tab):
# Evaluate the function numerically using lambda expressions
deriv_lambda = numpy_lambdify(X, phi[r])
data[r][i] = \
numpy.array(evaluate_lambda(deriv_lambda, pts.T)).T
# Symbolically compute the next derivative.
# This actually happens once too many here; never mind for
# now.
phi.append(form_derivative(phi[-1]))
return data


def xi_triangle(eta):
"""Maps from [-1,1]^2 to the (-1,1) reference triangle."""
eta1, eta2 = eta
xi1 = 0.5 * (1.0 + eta1) * (1.0 - eta2) - 1.0
xi2 = eta2
return (xi1, xi2)


def xi_tetrahedron(eta):
"""Maps from [-1,1]^3 to the -1/1 reference tetrahedron."""
eta1, eta2, eta3 = eta
xi1 = 0.25 * (1. + eta1) * (1. - eta2) * (1. - eta3) - 1.
xi2 = 0.5 * (1. + eta2) * (1. - eta3) - 1.
xi3 = eta3
return xi1, xi2, xi3


class ExpansionSet(object):
point_set = RecursivePointSet(lambda n: GaussLegendreQuadratureLineRule(UFCInterval(), n + 1).get_points())

Expand Down Expand Up @@ -389,42 +286,6 @@ def _tabulate_derivatives(self, n, pts):
self._normalize(n, dphi)
return phi, dphi

def _tabulate_duffy(self, n, pts):
"""Returns a dict of tabulations of phi_i(pts[j]) and each component of
the gradient d/dx_k phi_i(pts[j]). Here we employ the Duffy transform,
and thus this tabulation mode is only recommended for use with interior
points.
"""
from FIAT.polynomial_set import mis
dim = self.ref_el.get_spatial_dimension()
sd = self.get_num_members(n)
xi = numpy.transpose(numpy.dot(pts, self.A.T) + self.b)
eta = (lambda x: x, lambda x: x, eta_square, eta_cube)[dim](xi)
basis = [dubiner_1d(n, k, eta[k]) for k in range(dim)]
derivs = [dubiner_deriv_1d(n, k, eta[k]) for k in range(dim)]
tabulations = {}
alphas = mis(dim, 0) + mis(dim, 1)
for alpha in alphas:
V = [v if a == 0 else dv for a, v, dv in zip(alpha, basis, derivs)]
phi = V[0]
if dim >= 2:
phi1 = phi
phi = numpy.copy(V[1])
for i in range(n + 1):
indices = [morton_index2(i, j) for j in range(n + 1 - i)]
phi[indices] *= phi1[i]
if dim >= 3:
phi2 = phi
phi = numpy.zeros((sd, V[0].shape[1]), dtype=V[0].dtype)
for i in range(n + 1):
for j in range(n + 1 - i):
Vij = phi2[morton_index2(i, j)]
for k in range(n + 1 - i - j):
phi[morton_index3(i, j, k)] = V[2][morton_index2(i + j, k)] * Vij
tabulations[alpha] = phi
duffy_chain_rule(self.A, eta, tabulations)
return tabulations

def make_dmats(self, degree):
"""Returns a numpy array with the expansion coefficients dmat[k, j, i]
of the gradient of each member of the expansion set:
Expand Down Expand Up @@ -650,7 +511,7 @@ def _normalize(self, n, phi):
for p in range(n + 1):
for q in range(n - p + 1):
for r in range(n - p - q + 1):
phi[idx(p, q, r)] *= math.sqrt((p + 0.5) * (p + q + 1.0) * (p + q + r +1.5))
phi[idx(p, q, r)] *= math.sqrt((p + 0.5) * (p + q + 1.0) * (p + q + r + 1.5))

def tabulate_derivatives(self, n, pts):
order = 1
Expand Down

0 comments on commit d213473

Please sign in to comment.