Source code for FIAT.barycentric_interpolation
# Copyright (C) 2021 Pablo D. Brubeck
#
# This file is part of FIAT (https://www.fenicsproject.org)
#
# SPDX-License-Identifier: LGPL-3.0-or-later
#
# Written by Pablo D. Brubeck (brubeck@protonmail.com), 2021
import numpy
[docs]def barycentric_interpolation(xsrc, xdst, order=0):
"""Return tabulations of a 1D Lagrange nodal basis via the second barycentric interpolation formula
See Berrut and Trefethen (2004) https://doi.org/10.1137/S0036144502417715 Eq. (4.2) & (9.4)
:arg xsrc: a :class:`numpy.array` with the nodes defining the Lagrange polynomial basis
:arg xdst: a :class:`numpy.array` with the interpolation points
:arg order: the integer order of differentiation
:returns: dict of tabulations up to the given order (in the same format as :meth:`~.CiarletElement.tabulate`)
"""
# w = barycentric weights
# D = spectral differentiation matrix (D.T : u(xsrc) -> u'(xsrc))
# I = barycentric interpolation matrix (I.T : u(xsrc) -> u(xdst))
D = numpy.add.outer(-xsrc, xsrc)
numpy.fill_diagonal(D, 1.0E0)
w = 1.0E0 / numpy.prod(D, axis=0)
D = numpy.divide.outer(w, w) / D
numpy.fill_diagonal(D, numpy.diag(D) - numpy.sum(D, axis=0))
I = numpy.add.outer(-xsrc, xdst)
idx = numpy.argwhere(numpy.isclose(I, 0.0E0, 1E-14))
I[idx[:, 0], idx[:, 1]] = 1.0E0
I = 1.0E0 / I
I *= w[:, None]
I[:, idx[:, 1]] = 0.0E0
I[idx[:, 0], idx[:, 1]] = 1.0E0
I = (1.0E0 / numpy.sum(I, axis=0)) * I
derivs = {(0,): I}
for k in range(0, order):
derivs[(k+1,)] = numpy.matmul(D, derivs[(k,)])
return derivs