Source code for FIAT.pointwise_dual

# Copyright (C) 2020 Robert C. Kirby (Baylor University)
#
# This file is part of FIAT (https://www.fenicsproject.org)
#
# SPDX-License-Identifier:    LGPL-3.0-or-later
#
import numpy as np
from FIAT.functional import Functional
from FIAT.dual_set import DualSet
from collections import defaultdict
from itertools import zip_longest


[docs]def compute_pointwise_dual(el, pts): """Constructs a dual basis to the basis for el as a linear combination of a set of pointwise evaluations. This is useful when the prescribed finite element isn't Ciarlet (e.g. the basis functions are provided explicitly as formulae). Alternately, the element's given dual basis may involve differentiation, making run-time interpolation difficult in FIAT clients. The pointwise dual, consisting only of pointwise evaluations, will effectively replace these derivatives with (automatically determined) finite differences. This is exact on the polynomial space, but is an approximation if applied to functions outside the space. :param el: a :class:`FiniteElement`. :param pts: an iterable of points with the same length as el's dimension. These points must be unisolvent for the polynomial space :returns: a :class `DualSet` """ nbf = el.space_dimension() T = el.ref_el sd = T.get_spatial_dimension() assert np.asarray(pts).shape == (int(nbf / np.prod(el.value_shape())), sd) z = tuple([0] * sd) nds = [] V = el.tabulate(0, pts)[z] # Make a square system, invert, and then put it back in the right # shape so we have (nbf, ..., npts) with more dimensions # for vector or tensor-valued elements. alphas = np.linalg.inv(V.reshape((nbf, -1)).T).reshape(V.shape) # Each row of alphas gives the coefficients of a functional, # represented, as elsewhere in FIAT, as a summation of # components of the input at particular points. # This logic picks out the points and components for which the # weights are actually nonzero to construct the functional. pts = np.asarray(pts) for coeffs in alphas: pt_dict = defaultdict(list) nonzero = np.where(np.abs(coeffs) > 1.e-12) *comp, pt_index = nonzero for pt, coeff_comp in zip(pts[pt_index], zip_longest(coeffs[nonzero], zip(*comp), fillvalue=())): pt_dict[tuple(pt)].append(coeff_comp) nds.append(Functional(T, el.value_shape(), dict(pt_dict), {}, "node")) return DualSet(nds, T, el.entity_dofs())