forked from FEniCS/fiat
-
Notifications
You must be signed in to change notification settings - Fork 7
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
GL/GLL tabulation via barycentric interpolation (#65)
Implement numerical stable tabulation of (derivatives) of GL/GLL basis functions using the second barycentric interpolation formula of Berrut and Trefethen (2004).
- Loading branch information
Showing
3 changed files
with
92 additions
and
0 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,46 @@ | ||
# 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 ([email protected]), 2021 | ||
|
||
import numpy | ||
|
||
|
||
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 |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
|
@@ -5,9 +5,14 @@ | |
# SPDX-License-Identifier: LGPL-3.0-or-later | ||
# | ||
# Written by David A. Ham ([email protected]), 2015 | ||
# | ||
# Modified by Pablo D. Brubeck ([email protected]), 2021 | ||
|
||
import numpy | ||
|
||
from FIAT import finite_element, polynomial_set, dual_set, functional, quadrature | ||
from FIAT.reference_element import LINE | ||
from FIAT.barycentric_interpolation import barycentric_interpolation | ||
|
||
|
||
class GaussLegendreDualSet(dual_set.DualSet): | ||
|
@@ -31,3 +36,21 @@ def __init__(self, ref_el, degree): | |
dual = GaussLegendreDualSet(ref_el, degree) | ||
formdegree = ref_el.get_spatial_dimension() # n-form | ||
super(GaussLegendre, self).__init__(poly_set, dual, degree, formdegree) | ||
|
||
def tabulate(self, order, points, entity=None): | ||
# This overrides the default with a more numerically stable algorithm | ||
|
||
if entity is None: | ||
entity = (self.ref_el.get_dimension(), 0) | ||
|
||
entity_dim, entity_id = entity | ||
transform = self.ref_el.get_entity_transform(entity_dim, entity_id) | ||
|
||
xsrc = [] | ||
for node in self.dual.nodes: | ||
# Assert singleton point for each node. | ||
(pt,), = node.get_point_dict().keys() | ||
xsrc.append(pt) | ||
xsrc = numpy.asarray(xsrc) | ||
xdst = numpy.array(list(map(transform, points))).flatten() | ||
return barycentric_interpolation(xsrc, xdst, order=order) |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
|
@@ -5,9 +5,14 @@ | |
# SPDX-License-Identifier: LGPL-3.0-or-later | ||
# | ||
# Written by David A. Ham ([email protected]), 2015 | ||
# | ||
# Modified by Pablo D. Brubeck ([email protected]), 2021 | ||
|
||
import numpy | ||
|
||
from FIAT import finite_element, polynomial_set, dual_set, functional, quadrature | ||
from FIAT.reference_element import LINE | ||
from FIAT.barycentric_interpolation import barycentric_interpolation | ||
|
||
|
||
class GaussLobattoLegendreDualSet(dual_set.DualSet): | ||
|
@@ -31,3 +36,21 @@ def __init__(self, ref_el, degree): | |
dual = GaussLobattoLegendreDualSet(ref_el, degree) | ||
formdegree = 0 # 0-form | ||
super(GaussLobattoLegendre, self).__init__(poly_set, dual, degree, formdegree) | ||
|
||
def tabulate(self, order, points, entity=None): | ||
# This overrides the default with a more numerically stable algorithm | ||
|
||
if entity is None: | ||
entity = (self.ref_el.get_dimension(), 0) | ||
|
||
entity_dim, entity_id = entity | ||
transform = self.ref_el.get_entity_transform(entity_dim, entity_id) | ||
|
||
xsrc = [] | ||
for node in self.dual.nodes: | ||
# Assert singleton point for each node. | ||
(pt,), = node.get_point_dict().keys() | ||
xsrc.append(pt) | ||
xsrc = numpy.asarray(xsrc) | ||
xdst = numpy.array(list(map(transform, points))).flatten() | ||
return barycentric_interpolation(xsrc, xdst, order=order) |