Source code for FIAT.expansions

# 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
"""Principal orthogonal expansion functions as defined by Karniadakis
and Sherwin.  These are parametrized over a reference element so as
to allow users to get coordinates that they want."""

import numpy
import math
import sympy
from FIAT import reference_element
from FIAT import jacobi


[docs] def jrc(a, b, n): an = (2*n+1+a+b)*(2*n+2+a+b) / (2*(n+1)*(n+1+a+b)) bn = (a*a-b*b) * (2*n+1+a+b) / (2*(n+1)*(2*n+a+b)*(n+1+a+b)) cn = (n+a)*(n+b)*(2*n+2+a+b) / ((n+1)*(n+1+a+b)*(2*n+a+b)) 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
[docs] 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)
[docs] 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
[docs] class PointExpansionSet(object): """Evaluates the point basis on a point reference element.""" def __init__(self, ref_el): if ref_el.get_spatial_dimension() != 0: raise ValueError("Must have a point") self.ref_el = ref_el self.base_ref_el = reference_element.Point()
[docs] def get_num_members(self, n): return 1
[docs] def tabulate(self, n, pts): """Returns a numpy array A[i,j] = phi_i(pts[j]) = 1.0.""" assert n == 0 return numpy.ones((1, len(pts)))
[docs] def tabulate_derivatives(self, n, pts): """Returns a numpy array of size A where A[i,j] = phi_i(pts[j]) but where each element is an empty tuple (). This maintains compatibility with the interfaces of the interval, triangle and tetrahedron expansions.""" deriv_vals = numpy.empty_like(self.tabulate(n, pts), dtype=tuple) deriv_vals.fill(()) return deriv_vals
[docs] class LineExpansionSet(object): """Evaluates the Legendre basis on a line reference element.""" def __init__(self, ref_el): if ref_el.get_spatial_dimension() != 1: raise Exception("Must have a line") self.ref_el = ref_el self.base_ref_el = reference_element.DefaultLine() v1 = ref_el.get_vertices() v2 = self.base_ref_el.get_vertices() self.A, self.b = reference_element.make_affine_mapping(v1, v2) self.mapping = lambda x: numpy.dot(self.A, x) + self.b self.scale = numpy.sqrt(numpy.linalg.det(self.A))
[docs] def get_num_members(self, n): return n + 1
[docs] def tabulate(self, n, pts): """Returns a numpy array A[i,j] = phi_i(pts[j])""" if len(pts) > 0: ref_pts = numpy.array([self.mapping(pt) for pt in pts]) psitilde_as = jacobi.eval_jacobi_batch(0, 0, n, ref_pts) results = numpy.zeros((n + 1, len(pts)), type(pts[0][0])) for k in range(n + 1): results[k, :] = psitilde_as[k, :] * math.sqrt(k + 0.5) return results else: return []
[docs] def tabulate_derivatives(self, n, pts): """Returns a tuple of length one (A,) such that A[i,j] = D phi_i(pts[j]). The tuple is returned for compatibility with the interfaces of the triangle and tetrahedron expansions.""" ref_pts = numpy.array([self.mapping(pt) for pt in pts]) psitilde_as_derivs = jacobi.eval_jacobi_deriv_batch(0, 0, n, ref_pts) # Jacobi polynomials defined on [-1, 1], first derivatives need scaling psitilde_as_derivs *= 2.0 / self.ref_el.volume() results = numpy.zeros((n + 1, len(pts)), "d") for k in range(0, n + 1): results[k, :] = psitilde_as_derivs[k, :] * numpy.sqrt(k + 0.5) vals = self.tabulate(n, pts) deriv_vals = (results,) # Create the ordinary data structure. dv = [] for i in range(vals.shape[0]): dv.append([]) for j in range(vals.shape[1]): dv[-1].append((vals[i][j], [deriv_vals[0][i][j]])) return dv
[docs] class TriangleExpansionSet(object): """Evaluates the orthonormal Dubiner basis on a triangular reference element.""" def __init__(self, ref_el): if ref_el.get_spatial_dimension() != 2: raise Exception("Must have a triangle") self.ref_el = ref_el self.base_ref_el = reference_element.DefaultTriangle() v1 = ref_el.get_vertices() v2 = self.base_ref_el.get_vertices() self.A, self.b = reference_element.make_affine_mapping(v1, v2) self.mapping = lambda x: numpy.dot(self.A, x) + self.b # self.scale = numpy.sqrt(numpy.linalg.det(self.A))
[docs] def get_num_members(self, n): return (n + 1) * (n + 2) // 2
[docs] def tabulate(self, n, pts): if len(pts) == 0: return numpy.array([]) else: return numpy.array(self._tabulate(n, numpy.array(pts).T))
def _tabulate(self, n, pts): '''A version of tabulate() that also works for a single point. ''' m1, m2 = self.A.shape ref_pts = [sum(self.A[i][j] * pts[j] for j in range(m2)) + self.b[i] for i in range(m1)] def idx(p, q): return (p + q) * (p + q + 1) // 2 + q results = ((n + 1) * (n + 2) // 2) * [None] results[0] = 1.0 \ + pts[0] - pts[0] \ + pts[1] - pts[1] if n == 0: return results x = ref_pts[0] y = ref_pts[1] f1 = (1.0 + 2 * x + y) / 2.0 f2 = (1.0 - y) / 2.0 f3 = f2**2 results[idx(1, 0)] = f1 for p in range(1, n): a = (2.0 * p + 1) / (1.0 + p) # b = p / (p+1.0) results[idx(p+1, 0)] = a * f1 * results[idx(p, 0)] \ - p/(1.0+p) * f3 * results[idx(p-1, 0)] for p in range(n): results[idx(p, 1)] = 0.5 * (1+2.0*p+(3.0+2.0*p)*y) \ * results[idx(p, 0)] for p in range(n - 1): for q in range(1, n - p): (a1, a2, a3) = jrc(2 * p + 1, 0, q) results[idx(p, q+1)] = \ (a1 * y + a2) * results[idx(p, q)] \ - a3 * results[idx(p, q-1)] for p in range(n + 1): for q in range(n - p + 1): results[idx(p, q)] *= math.sqrt((p + 0.5) * (p + q + 1.0)) return results # return self.scale * results
[docs] def tabulate_derivatives(self, n, pts): order = 1 data = _tabulate_dpts(self._tabulate, 2, n, order, numpy.array(pts)) # Put data in the required data structure, i.e., # k-tuples which contain the value, and the k-1 derivatives # (gradient, Hessian, ...) m = data[0].shape[0] n = data[0].shape[1] data2 = [[tuple([data[r][i][j] for r in range(order+1)]) for j in range(n)] for i in range(m)] return data2
[docs] def tabulate_jet(self, n, pts, order=1): return _tabulate_dpts(self._tabulate, 2, n, order, numpy.array(pts))
[docs] class TetrahedronExpansionSet(object): """Collapsed orthonormal polynomial expanion on a tetrahedron.""" def __init__(self, ref_el): if ref_el.get_spatial_dimension() != 3: raise Exception("Must be a tetrahedron") self.ref_el = ref_el self.base_ref_el = reference_element.DefaultTetrahedron() v1 = ref_el.get_vertices() v2 = self.base_ref_el.get_vertices() self.A, self.b = reference_element.make_affine_mapping(v1, v2) self.mapping = lambda x: numpy.dot(self.A, x) + self.b self.scale = numpy.sqrt(numpy.linalg.det(self.A))
[docs] def get_num_members(self, n): return (n + 1) * (n + 2) * (n + 3) // 6
[docs] def tabulate(self, n, pts): if len(pts) == 0: return numpy.array([]) else: return numpy.array(self._tabulate(n, numpy.array(pts).T))
def _tabulate(self, n, pts): '''A version of tabulate() that also works for a single point. ''' m1, m2 = self.A.shape ref_pts = [sum(self.A[i][j] * pts[j] for j in range(m2)) + self.b[i] for i in range(m1)] def idx(p, q, r): return (p + q + r)*(p + q + r + 1)*(p + q + r + 2)//6 + (q + r)*(q + r + 1)//2 + r results = ((n + 1) * (n + 2) * (n + 3) // 6) * [None] results[0] = 1.0 \ + pts[0] - pts[0] \ + pts[1] - pts[1] \ + pts[2] - pts[2] if n == 0: return results x = ref_pts[0] y = ref_pts[1] z = ref_pts[2] factor1 = 0.5 * (2.0 + 2.0 * x + y + z) factor2 = (0.5 * (y + z))**2 factor3 = 0.5 * (1 + 2.0 * y + z) factor4 = 0.5 * (1 - z) factor5 = factor4**2 results[idx(1, 0, 0)] = factor1 for p in range(1, n): a1 = (2.0 * p + 1.0) / (p + 1.0) a2 = p / (p + 1.0) results[idx(p+1, 0, 0)] = a1 * factor1 * results[idx(p, 0, 0)] \ - a2 * factor2 * results[idx(p-1, 0, 0)] # q = 1 for p in range(0, n): results[idx(p, 1, 0)] = results[idx(p, 0, 0)] \ * (p * (1.0 + y) + (2.0 + 3.0 * y + z) / 2) for p in range(0, n - 1): for q in range(1, n - p): (aq, bq, cq) = jrc(2 * p + 1, 0, q) qmcoeff = aq * factor3 + bq * factor4 qm1coeff = cq * factor5 results[idx(p, q+1, 0)] = qmcoeff * results[idx(p, q, 0)] \ - qm1coeff * results[idx(p, q-1, 0)] # now handle r=1 for p in range(n): for q in range(n - p): results[idx(p, q, 1)] = results[idx(p, q, 0)] \ * (1.0 + p + q + (2.0 + q + p) * z) # general r by recurrence for p in range(n - 1): for q in range(0, n - p - 1): for r in range(1, n - p - q): ar, br, cr = jrc(2 * p + 2 * q + 2, 0, r) results[idx(p, q, r+1)] = \ (ar * z + br) * results[idx(p, q, r)] \ - cr * results[idx(p, q, r-1)] for p in range(n + 1): for q in range(n - p + 1): for r in range(n - p - q + 1): results[idx(p, q, r)] *= \ math.sqrt((p+0.5)*(p+q+1.0)*(p+q+r+1.5)) return results
[docs] def tabulate_derivatives(self, n, pts): order = 1 D = 3 data = _tabulate_dpts(self._tabulate, D, n, order, numpy.array(pts)) # Put data in the required data structure, i.e., # k-tuples which contain the value, and the k-1 derivatives # (gradient, Hessian, ...) m = data[0].shape[0] n = data[0].shape[1] data2 = [[tuple([data[r][i][j] for r in range(order + 1)]) for j in range(n)] for i in range(m)] return data2
[docs] def tabulate_jet(self, n, pts, order=1): return _tabulate_dpts(self._tabulate, 3, n, order, numpy.array(pts))
[docs] def get_expansion_set(ref_el): """Returns an ExpansionSet instance appopriate for the given reference element.""" if ref_el.get_shape() == reference_element.POINT: return PointExpansionSet(ref_el) elif ref_el.get_shape() == reference_element.LINE: return LineExpansionSet(ref_el) elif ref_el.get_shape() == reference_element.TRIANGLE: return TriangleExpansionSet(ref_el) elif ref_el.get_shape() == reference_element.TETRAHEDRON: return TetrahedronExpansionSet(ref_el) else: raise Exception("Unknown reference element type.")
[docs] def polynomial_dimension(ref_el, degree): """Returns the dimension of the space of polynomials of degree no greater than degree on the reference element.""" if ref_el.get_shape() == reference_element.POINT: if degree > 0: raise ValueError("Only degree zero polynomials supported on point elements.") return 1 elif ref_el.get_shape() == reference_element.LINE: return max(0, degree + 1) elif ref_el.get_shape() == reference_element.TRIANGLE: return max((degree + 1) * (degree + 2) // 2, 0) elif ref_el.get_shape() == reference_element.TETRAHEDRON: return max(0, (degree + 1) * (degree + 2) * (degree + 3) // 6) else: raise ValueError("Unknown reference element type.")