Source code for FIAT.bernstein

# -*- coding: utf-8 -*-
#
# Copyright (C) 2018 Miklós Homolya
#
# This file is part of FIAT (https://www.fenicsproject.org)
#
# SPDX-License-Identifier:    LGPL-3.0-or-later

import math
import numpy

from FIAT.finite_element import FiniteElement
from FIAT.dual_set import DualSet
from FIAT.polynomial_set import mis


[docs] class BernsteinDualSet(DualSet): """The dual basis for Bernstein elements.""" def __init__(self, ref_el, degree): # Initialise data structures topology = ref_el.get_topology() entity_ids = {dim: {entity_i: [] for entity_i in entities} for dim, entities in topology.items()} # Calculate inverse topology inverse_topology = {vertices: (dim, entity_i) for dim, entities in topology.items() for entity_i, vertices in entities.items()} # Generate triangular barycentric indices dim = ref_el.get_spatial_dimension() kss = mis(dim + 1, degree) # Fill data structures nodes = [] for i, ks in enumerate(kss): vertices, = numpy.nonzero(ks) entity_dim, entity_i = inverse_topology[tuple(vertices)] entity_ids[entity_dim][entity_i].append(i) # Leave nodes unimplemented for now nodes.append(None) super(BernsteinDualSet, self).__init__(nodes, ref_el, entity_ids)
[docs] class Bernstein(FiniteElement): """A finite element with Bernstein polynomials as basis functions.""" def __init__(self, ref_el, degree): dual = BernsteinDualSet(ref_el, degree) k = 0 # 0-form super(Bernstein, self).__init__(ref_el, dual, degree, k)
[docs] def degree(self): """The degree of the polynomial space.""" return self.get_order()
[docs] def value_shape(self): """The value shape of the finite element functions.""" return ()
[docs] def tabulate(self, order, points, entity=None): """Return tabulated values of derivatives up to given order of basis functions at given points. :arg order: The maximum order of derivative. :arg points: An iterable of points. :arg entity: Optional (dimension, entity number) pair indicating which topological entity of the reference element to tabulate on. If ``None``, default cell-wise tabulation is performed. """ # Transform points to reference cell coordinates ref_el = self.get_reference_element() if entity is None: entity = (ref_el.get_spatial_dimension(), 0) entity_dim, entity_id = entity entity_transform = ref_el.get_entity_transform(entity_dim, entity_id) cell_points = list(map(entity_transform, points)) # Construct Cartesian to Barycentric coordinate mapping vs = numpy.asarray(ref_el.get_vertices()) B2R = numpy.vstack([vs.T, numpy.ones(len(vs))]) R2B = numpy.linalg.inv(B2R) B = numpy.hstack([cell_points, numpy.ones((len(cell_points), 1))]).dot(R2B.T) # Evaluate everything deg = self.degree() dim = ref_el.get_spatial_dimension() raw_result = {(alpha, i): vec for i, ks in enumerate(mis(dim + 1, deg)) for o in range(order + 1) for alpha, vec in bernstein_Dx(B, ks, o, R2B).items()} # Rearrange result space_dim = self.space_dimension() dtype = numpy.array(list(raw_result.values())).dtype result = {alpha: numpy.zeros((space_dim, len(cell_points)), dtype=dtype) for o in range(order + 1) for alpha in mis(dim, o)} for (alpha, i), vec in raw_result.items(): result[alpha][i, :] = vec return result
[docs] def bernstein_db(points, ks, alpha=None): """Evaluates Bernstein polynomials or its derivative at barycentric points. :arg points: array of points in barycentric coordinates :arg ks: exponents defining the Bernstein polynomial :arg alpha: derivative tuple :returns: array of Bernstein polynomial values at given points. """ points = numpy.asarray(points) ks = numpy.array(tuple(ks)) N, d_1 = points.shape assert d_1 == len(ks) if alpha is None: alpha = numpy.zeros(d_1) else: alpha = numpy.array(tuple(alpha)) assert d_1 == len(alpha) ls = ks - alpha if any(k < 0 for k in ls): return numpy.zeros(len(points)) elif all(k == 0 for k in ls): return numpy.ones(len(points)) else: # Calculate coefficient coeff = math.factorial(ks.sum()) for k in ls: coeff //= math.factorial(k) return coeff * numpy.prod(points**ls, axis=1)
[docs] def bernstein_Dx(points, ks, order, R2B): """Evaluates Bernstein polynomials or its derivatives according to reference coordinates. :arg points: array of points in BARYCENTRIC COORDINATES :arg ks: exponents defining the Bernstein polynomial :arg alpha: derivative order (returns all derivatives of this specified order) :arg R2B: linear mapping from reference to barycentric coordinates :returns: dictionary mapping from derivative tuples to arrays of Bernstein polynomial values at given points. """ points = numpy.asarray(points) ks = tuple(ks) N, d_1 = points.shape assert d_1 == len(ks) # Collect derivatives according to barycentric coordinates Db_map = {alpha: bernstein_db(points, ks, alpha) for alpha in mis(d_1, order)} # Arrange derivative tensor (barycentric coordinates) dtype = numpy.array(list(Db_map.values())).dtype Db_shape = (d_1,) * order Db_tensor = numpy.empty(Db_shape + (N,), dtype=dtype) for ds in numpy.ndindex(Db_shape): alpha = [0] * d_1 for d in ds: alpha[d] += 1 Db_tensor[ds + (slice(None),)] = Db_map[tuple(alpha)] # Coordinate transformation: barycentric -> reference result = {} for alpha in mis(d_1 - 1, order): values = Db_tensor for d, k in enumerate(alpha): for _ in range(k): values = R2B[:, d].dot(values) result[alpha] = values return result