Source code for FIAT.functional

# Copyright (C) 2008 Robert C. Kirby (Texas Tech University)
#
# Modified 2020 by the same from Baylor University
#
# This file is part of FIAT (https://www.fenicsproject.org)
#
# SPDX-License-Identifier:    LGPL-3.0-or-later

# functionals require:
# - a degree of accuracy (-1 indicates that it works for all functions
#   such as point evaluation)
# - a reference element domain
# - type information

from collections import OrderedDict
from itertools import chain
import numpy
import sympy

from FIAT import polynomial_set
from FIAT.quadrature import GaussLegendreQuadratureLineRule, QuadratureRule
from FIAT.reference_element import UFCInterval as interval


[docs]def index_iterator(shp): """Constructs a generator iterating over all indices in shp in generalized column-major order So if shp = (2,2), then we construct the sequence (0,0),(0,1),(1,0),(1,1)""" if len(shp) == 0: return elif len(shp) == 1: for i in range(shp[0]): yield [i] else: shp_foo = shp[1:] for i in range(shp[0]): for foo in index_iterator(shp_foo): yield [i] + foo
[docs]class Functional(object): r"""Abstract class representing a linear functional. All FIAT functionals are discrete in the sense that they are written as a weighted sum of (derivatives of components of) their argument evaluated at particular points. :arg ref_el: a :class:`Cell` :arg target_shape: a tuple indicating the value shape of functions on the functional operates (e.g. if the function eats 2-vectors then target_shape is (2,) and if it eats scalars then target_shape is () :arg pt_dict: A dict mapping points to lists of information about how the functional is evaluated. Each entry in the list takes the form of a tuple (wt, comp) so that (at least if the deriv_dict argument is empty), the functional takes the form :math:`\ell(f) = \sum_{q=1}^{N_q} \sum_{k=1}^{K_q} w^q_k f_{c_k}(x_q)` where :math:`f_{c_k}` indicates a particular vector or tensor component :arg deriv_dict: A dict that is similar to `pt_dict`, although the entries of each list are tuples (wt, alpha, comp) with alpha a tuple of nonnegative integers corresponding to the order of partial differentiation in each spatial direction. :arg functional_type: a string labeling the kind of functional this is. """ def __init__(self, ref_el, target_shape, pt_dict, deriv_dict, functional_type): self.ref_el = ref_el self.target_shape = target_shape self.pt_dict = pt_dict self.deriv_dict = deriv_dict self.functional_type = functional_type if len(deriv_dict) > 0: per_point = list(chain(*deriv_dict.values())) alphas = [tuple(foo[1]) for foo in per_point] self.max_deriv_order = max([sum(foo) for foo in alphas]) else: self.max_deriv_order = 0
[docs] def evaluate(self, f): """Obsolete and broken functional evaluation. To evaluate the functional, call it on the target function: functional(function) """ raise AttributeError("To evaluate the functional just call it on a function.")
def __call__(self, fn): raise NotImplementedError("Evaluation is not yet implemented for %s" % type(self))
[docs] def get_point_dict(self): """Returns the functional information, which is a dictionary mapping each point in the support of the functional to a list of pairs containing the weight and component.""" return self.pt_dict
[docs] def get_reference_element(self): """Returns the reference element.""" return self.ref_el
[docs] def get_type_tag(self): """Returns the type of function (e.g. point evaluation or normal component, which is probably handy for clients of FIAT""" return self.functional_type
[docs] def to_riesz(self, poly_set): r"""Constructs an array representation of the functional so that the functional may be applied to a function expressed in in terms of the expansion set underlying `poly_set` by means of contracting coefficients. That is, `poly_set` will have members all expressed in the form :math:`p = \sum_{i} \alpha^i \phi_i` where :math:`\{\phi_i\}_{i}` is some orthonormal expansion set and :math:`\alpha^i` are coefficients. Note: the orthonormal expansion set is always scalar-valued but if the members of `poly_set` are vector or tensor valued the :math:`\alpha^i` will be scalars or vectors. This function constructs a tensor :math:`R` such that the contraction of :math:`R` with the array of coefficients :math:`\alpha` produces the effect of :math:`\ell(f)` In the case of scalar-value functions, :math:`R` is just a vector of the same length as the expansion set, and :math:`R_i = \ell(\phi_i)`. For vector-valued spaces, :math:`R_{ij}` will be :math:`\ell(e^i \phi_j)` where :math:`e^i` is the canonical unit vector nonzero only in one entry :math:`i`. """ es = poly_set.get_expansion_set() ed = poly_set.get_embedded_degree() nexp = es.get_num_members(ed) pt_dict = self.get_point_dict() pts = list(pt_dict.keys()) npts = len(pts) bfs = es.tabulate(ed, pts) result = numpy.zeros(poly_set.coeffs.shape[1:], "d") # loop over points for j in range(npts): pt_cur = pts[j] wc_list = pt_dict[pt_cur] # loop over expansion functions for i in range(nexp): for (w, c) in wc_list: result[c][i] += w * bfs[i, j] if self.deriv_dict: dpt_dict = self.deriv_dict # this makes things quicker since it uses dmats after # instantiation es_foo = polynomial_set.ONPolynomialSet(self.ref_el, ed) dpts = list(dpt_dict.keys()) dbfs = es_foo.tabulate(dpts, self.max_deriv_order) ndpts = len(dpts) for j in range(ndpts): dpt_cur = dpts[j] wac_list = dpt_dict[dpt_cur] for i in range(nexp): for (w, alpha, c) in wac_list: result[c][i] += w * dbfs[tuple(alpha)][i, j] return result
[docs] def tostr(self): return self.functional_type
[docs]class PointEvaluation(Functional): """Class representing point evaluation of scalar functions at a particular point x.""" def __init__(self, ref_el, x): pt_dict = {x: [(1.0, tuple())]} Functional.__init__(self, ref_el, tuple(), pt_dict, {}, "PointEval") def __call__(self, fn): """Evaluate the functional on the function fn.""" return fn(tuple(self.pt_dict.keys())[0])
[docs] def tostr(self): x = list(map(str, list(self.pt_dict.keys())[0])) return "u(%s)" % (','.join(x),)
[docs]class ComponentPointEvaluation(Functional): """Class representing point evaluation of a particular component of a vector function at a particular point x.""" def __init__(self, ref_el, comp, shp, x): if len(shp) != 1: raise Exception("Illegal shape") if comp < 0 or comp >= shp[0]: raise Exception("Illegal component") self.comp = comp pt_dict = {x: [(1.0, (comp,))]} Functional.__init__(self, ref_el, shp, pt_dict, {}, "ComponentPointEval")
[docs] def tostr(self): x = list(map(str, list(self.pt_dict.keys())[0])) return "(u[%d](%s)" % (self.comp, ','.join(x))
[docs]class PointDerivative(Functional): """Class representing point partial differentiation of scalar functions at a particular point x.""" def __init__(self, ref_el, x, alpha): dpt_dict = {x: [(1.0, tuple(alpha), tuple())]} self.alpha = tuple(alpha) self.order = sum(self.alpha) Functional.__init__(self, ref_el, tuple(), {}, dpt_dict, "PointDeriv") def __call__(self, fn): """Evaluate the functional on the function fn. Note that this depends on sympy being able to differentiate fn.""" x = list(self.deriv_dict.keys())[0] X = sympy.DeferredVector('x') dX = numpy.asarray([X[i] for i in range(len(x))]) dvars = tuple(d for d, a in zip(dX, self.alpha) for count in range(a)) return sympy.diff(fn(X), *dvars).evalf(subs=dict(zip(dX, x)))
[docs]class PointNormalDerivative(Functional): """Represents d/dn at a point on a facet.""" def __init__(self, ref_el, facet_no, pt): n = ref_el.compute_normal(facet_no) self.n = n sd = ref_el.get_spatial_dimension() alphas = [] for i in range(sd): alpha = [0] * sd alpha[i] = 1 alphas.append(alpha) dpt_dict = {pt: [(n[i], tuple(alphas[i]), tuple()) for i in range(sd)]} Functional.__init__(self, ref_el, tuple(), {}, dpt_dict, "PointNormalDeriv")
[docs]class PointNormalSecondDerivative(Functional): """Represents d^/dn^2 at a point on a facet.""" def __init__(self, ref_el, facet_no, pt): n = ref_el.compute_normal(facet_no) self.n = n sd = ref_el.get_spatial_dimension() tau = numpy.zeros((sd*(sd+1)//2,)) alphas = [] cur = 0 for i in range(sd): for j in range(i, sd): alpha = [0] * sd alpha[i] += 1 alpha[j] += 1 alphas.append(tuple(alpha)) tau[cur] = n[i]*n[j] cur += 1 self.tau = tau self.alphas = alphas dpt_dict = {pt: [(n[i], alphas[i], tuple()) for i in range(sd)]} Functional.__init__(self, ref_el, tuple(), {}, dpt_dict, "PointNormalDeriv")
[docs]class IntegralMoment(Functional): """Functional representing integral of the input against some tabulated function f. :arg ref_el: a :class:`Cell`. :arg Q: a :class:`QuadratureRule`. :arg f_at_qpts: an array tabulating the function f at the quadrature points. :arg comp: Optional argument indicating that only a particular component of the input function should be integrated against f :arg shp: Optional argument giving the value shape of input functions. """ def __init__(self, ref_el, Q, f_at_qpts, comp=tuple(), shp=tuple()): self.Q = Q qpts, qwts = Q.get_points(), Q.get_weights() pt_dict = OrderedDict() self.comp = comp for i in range(len(qpts)): pt_cur = tuple(qpts[i]) pt_dict[pt_cur] = [(qwts[i] * f_at_qpts[i], comp)] Functional.__init__(self, ref_el, shp, pt_dict, {}, "IntegralMoment") def __call__(self, fn): """Evaluate the functional on the function fn.""" pts = list(self.pt_dict.keys()) wts = numpy.array([foo[0][0] for foo in list(self.pt_dict.values())]) result = numpy.dot([fn(p) for p in pts], wts) if self.comp: result = result[self.comp] return result
[docs]class IntegralMomentOfNormalDerivative(Functional): """Functional giving normal derivative integrated against some function on a facet.""" def __init__(self, ref_el, facet_no, Q, f_at_qpts): n = ref_el.compute_normal(facet_no) self.n = n self.f_at_qpts = f_at_qpts self.Q = Q sd = ref_el.get_spatial_dimension() # map points onto facet fmap = ref_el.get_entity_transform(sd-1, facet_no) qpts, qwts = Q.get_points(), Q.get_weights() dpts = [fmap(pt) for pt in qpts] self.dpts = dpts dpt_dict = OrderedDict() alphas = [tuple([1 if j == i else 0 for j in range(sd)]) for i in range(sd)] for j, pt in enumerate(dpts): dpt_dict[tuple(pt)] = [(qwts[j]*n[i]*f_at_qpts[j], alphas[i], tuple()) for i in range(sd)] Functional.__init__(self, ref_el, tuple(), {}, dpt_dict, "IntegralMomentOfNormalDerivative")
[docs]class IntegralLegendreDirectionalMoment(Functional): """Moment of v.s against a Legendre polynomial over an edge""" def __init__(self, cell, s, entity, mom_deg, comp_deg, nm=""): sd = cell.get_spatial_dimension() assert sd == 2 shp = (sd,) quadpoints = comp_deg + 1 Q = GaussLegendreQuadratureLineRule(interval(), quadpoints) legendre = numpy.polynomial.legendre.legval(2*Q.get_points()-1, [0]*mom_deg + [1]) f_at_qpts = numpy.array([s*legendre[i] for i in range(quadpoints)]) fmap = cell.get_entity_transform(sd-1, entity) mappedqpts = [fmap(pt) for pt in Q.get_points()] mappedQ = QuadratureRule(cell, mappedqpts, Q.get_weights()) qwts = mappedQ.wts qpts = mappedQ.pts pt_dict = OrderedDict() for k in range(len(qpts)): pt_cur = tuple(qpts[k]) pt_dict[pt_cur] = [(qwts[k] * f_at_qpts[k, i], (i,)) for i in range(2)] super().__init__(cell, shp, pt_dict, {}, nm)
[docs]class IntegralLegendreNormalMoment(IntegralLegendreDirectionalMoment): """Moment of v.n against a Legendre polynomial over an edge""" def __init__(self, cell, entity, mom_deg, comp_deg): n = cell.compute_scaled_normal(entity) super().__init__(cell, n, entity, mom_deg, comp_deg, "IntegralLegendreNormalMoment")
[docs]class IntegralLegendreTangentialMoment(IntegralLegendreDirectionalMoment): """Moment of v.t against a Legendre polynomial over an edge""" def __init__(self, cell, entity, mom_deg, comp_deg): t = cell.compute_edge_tangent(entity) super().__init__(cell, t, entity, mom_deg, comp_deg, "IntegralLegendreTangentialMoment")
[docs]class IntegralLegendreBidirectionalMoment(Functional): """Moment of dot(s1, dot(tau, s2)) against Legendre on entity, multiplied by the size of the reference facet""" def __init__(self, cell, s1, s2, entity, mom_deg, comp_deg, nm=""): # mom_deg is degree of moment, comp_deg is the total degree of # polynomial you might need to integrate (or something like that) sd = cell.get_spatial_dimension() shp = (sd, sd) s1s2T = numpy.outer(s1, s2) quadpoints = comp_deg + 1 Q = GaussLegendreQuadratureLineRule(interval(), quadpoints) # The volume squared gets the Jacobian mapping from line interval # and the edge length into the functional. legendre = numpy.polynomial.legendre.legval(2*Q.get_points()-1, [0]*mom_deg + [1]) * numpy.abs(cell.volume_of_subcomplex(1, entity))**2 f_at_qpts = numpy.array([s1s2T*legendre[i] for i in range(quadpoints)]) # Map the quadrature points fmap = cell.get_entity_transform(sd-1, entity) mappedqpts = [fmap(pt) for pt in Q.get_points()] mappedQ = QuadratureRule(cell, mappedqpts, Q.get_weights()) pt_dict = OrderedDict() qpts = mappedQ.pts qwts = mappedQ.wts for k in range(len(qpts)): pt_cur = tuple(qpts[k]) pt_dict[pt_cur] = [(qwts[k] * f_at_qpts[k, i, j], (i, j)) for (i, j) in index_iterator(shp)] super().__init__(cell, shp, pt_dict, {}, nm)
[docs]class IntegralLegendreNormalNormalMoment(IntegralLegendreBidirectionalMoment): """Moment of dot(n, dot(tau, n)) against Legendre on entity.""" def __init__(self, cell, entity, mom_deg, comp_deg): n = cell.compute_normal(entity) super().__init__(cell, n, n, entity, mom_deg, comp_deg, "IntegralNormalNormalLegendreMoment")
[docs]class IntegralLegendreNormalTangentialMoment(IntegralLegendreBidirectionalMoment): """Moment of dot(n, dot(tau, t)) against Legendre on entity.""" def __init__(self, cell, entity, mom_deg, comp_deg): n = cell.compute_normal(entity) t = cell.compute_normalized_edge_tangent(entity) super().__init__(cell, n, t, entity, mom_deg, comp_deg, "IntegralNormalTangentialLegendreMoment")
[docs]class IntegralMomentOfDivergence(Functional): """Functional representing integral of the divergence of the input against some tabulated function f.""" def __init__(self, ref_el, Q, f_at_qpts): self.f_at_qpts = f_at_qpts self.Q = Q sd = ref_el.get_spatial_dimension() qpts, qwts = Q.get_points(), Q.get_weights() dpts = qpts self.dpts = dpts dpt_dict = OrderedDict() alphas = [tuple([1 if j == i else 0 for j in range(sd)]) for i in range(sd)] for j, pt in enumerate(dpts): dpt_dict[tuple(pt)] = [(qwts[j]*f_at_qpts[j], alphas[i], (i,)) for i in range(sd)] super().__init__(ref_el, tuple(), {}, dpt_dict, "IntegralMomentOfDivergence")
[docs]class IntegralMomentOfTensorDivergence(Functional): """Like IntegralMomentOfDivergence, but on symmetric tensors.""" def __init__(self, ref_el, Q, f_at_qpts): self.f_at_qpts = f_at_qpts self.Q = Q qpts, qwts = Q.get_points(), Q.get_weights() nqp = len(qpts) dpts = qpts self.dpts = dpts assert len(f_at_qpts.shape) == 2 assert f_at_qpts.shape[0] == 2 assert f_at_qpts.shape[1] == nqp sd = ref_el.get_spatial_dimension() dpt_dict = OrderedDict() alphas = [tuple([1 if j == i else 0 for j in range(sd)]) for i in range(sd)] for q, pt in enumerate(dpts): dpt_dict[tuple(pt)] = [(qwts[q]*f_at_qpts[i, q], alphas[j], (i, j)) for i in range(2) for j in range(2)] super().__init__(ref_el, tuple(), {}, dpt_dict, "IntegralMomentOfDivergence")
[docs]class FrobeniusIntegralMoment(Functional): def __init__(self, ref_el, Q, f_at_qpts): # f_at_qpts is (some shape) x num_qpts shp = tuple(f_at_qpts.shape[:-1]) if len(Q.get_points()) != f_at_qpts.shape[-1]: raise Exception("Mismatch in number of quadrature points and values") qpts, qwts = Q.get_points(), Q.get_weights() pt_dict = {} for i, (pt_cur, wt_cur) in enumerate(zip(map(tuple, qpts), qwts)): pt_dict[pt_cur] = [] for alfa in index_iterator(shp): qpidx = tuple(alfa + [i]) pt_dict[pt_cur].append((wt_cur * f_at_qpts[qpidx], tuple(alfa))) super().__init__(ref_el, shp, pt_dict, {}, "FrobeniusIntegralMoment")
[docs]class PointNormalEvaluation(Functional): """Implements the evaluation of the normal component of a vector at a point on a facet of codimension 1.""" def __init__(self, ref_el, facet_no, pt): n = ref_el.compute_normal(facet_no) self.n = n sd = ref_el.get_spatial_dimension() pt_dict = {pt: [(n[i], (i,)) for i in range(sd)]} shp = (sd,) super().__init__(ref_el, shp, pt_dict, {}, "PointNormalEval")
[docs]class PointEdgeTangentEvaluation(Functional): """Implements the evaluation of the tangential component of a vector at a point on a facet of dimension 1.""" def __init__(self, ref_el, edge_no, pt): t = ref_el.compute_edge_tangent(edge_no) self.t = t sd = ref_el.get_spatial_dimension() pt_dict = {pt: [(t[i], (i,)) for i in range(sd)]} shp = (sd,) super().__init__(ref_el, shp, pt_dict, {}, "PointEdgeTangent")
[docs] def tostr(self): x = list(map(str, list(self.pt_dict.keys())[0])) return "(u.t)(%s)" % (','.join(x),)
[docs]class IntegralMomentOfEdgeTangentEvaluation(Functional): r""" \int_e v\cdot t p ds p \in Polynomials :arg ref_el: reference element for which e is a dim-1 entity :arg Q: quadrature rule on the face :arg P_at_qpts: polynomials evaluated at quad points :arg edge: which edge. """ def __init__(self, ref_el, Q, P_at_qpts, edge): t = ref_el.compute_edge_tangent(edge) sd = ref_el.get_spatial_dimension() transform = ref_el.get_entity_transform(1, edge) pts = tuple(map(lambda p: tuple(transform(p)), Q.get_points())) weights = Q.get_weights() pt_dict = OrderedDict() for pt, wgt, phi in zip(pts, weights, P_at_qpts): pt_dict[pt] = [(wgt*phi*t[i], (i, )) for i in range(sd)] super().__init__(ref_el, (sd, ), pt_dict, {}, "IntegralMomentOfEdgeTangentEvaluation")
[docs]class PointFaceTangentEvaluation(Functional): """Implements the evaluation of a tangential component of a vector at a point on a facet of codimension 1.""" def __init__(self, ref_el, face_no, tno, pt): t = ref_el.compute_face_tangents(face_no)[tno] self.t = t self.tno = tno sd = ref_el.get_spatial_dimension() pt_dict = {pt: [(t[i], (i,)) for i in range(sd)]} shp = (sd,) Functional.__init__(self, ref_el, shp, pt_dict, {}, "PointFaceTangent")
[docs] def tostr(self): x = list(map(str, list(self.pt_dict.keys())[0])) return "(u.t%d)(%s)" % (self.tno, ','.join(x),)
[docs]class IntegralMomentOfFaceTangentEvaluation(Functional): r""" \int_F v \times n \cdot p ds p \in Polynomials :arg ref_el: reference element for which F is a codim-1 entity :arg Q: quadrature rule on the face :arg P_at_qpts: polynomials evaluated at quad points :arg facet: which facet. """ def __init__(self, ref_el, Q, P_at_qpts, facet): P_at_qpts = [[P_at_qpts[0][i], P_at_qpts[1][i], P_at_qpts[2][i]] for i in range(P_at_qpts.shape[1])] n = ref_el.compute_scaled_normal(facet) sd = ref_el.get_spatial_dimension() transform = ref_el.get_entity_transform(sd-1, facet) pts = tuple(map(lambda p: tuple(transform(p)), Q.get_points())) weights = Q.get_weights() pt_dict = OrderedDict() for pt, wgt, phi in zip(pts, weights, P_at_qpts): phixn = [phi[1]*n[2] - phi[2]*n[1], phi[2]*n[0] - phi[0]*n[2], phi[0]*n[1] - phi[1]*n[0]] pt_dict[pt] = [(wgt*(-n[2]*phixn[1]+n[1]*phixn[2]), (0, )), (wgt*(n[2]*phixn[0]-n[0]*phixn[2]), (1, )), (wgt*(-n[1]*phixn[0]+n[0]*phixn[1]), (2, ))] super().__init__(ref_el, (sd, ), pt_dict, {}, "IntegralMomentOfFaceTangentEvaluation")
[docs]class MonkIntegralMoment(Functional): r""" face nodes are \int_F v\cdot p dA where p \in P_{q-2}(f)^3 with p \cdot n = 0 (cmp. Peter Monk - Finite Element Methods for Maxwell's equations p. 129) Note that we don't scale by the area of the facet :arg ref_el: reference element for which F is a codim-1 entity :arg Q: quadrature rule on the face :arg P_at_qpts: polynomials evaluated at quad points :arg facet: which facet. """ def __init__(self, ref_el, Q, P_at_qpts, facet): sd = ref_el.get_spatial_dimension() weights = Q.get_weights() pt_dict = OrderedDict() transform = ref_el.get_entity_transform(sd-1, facet) pts = tuple(map(lambda p: tuple(transform(p)), Q.get_points())) for pt, wgt, phi in zip(pts, weights, P_at_qpts): pt_dict[pt] = [(wgt*phi[i], (i, )) for i in range(sd)] super().__init__(ref_el, (sd, ), pt_dict, {}, "MonkIntegralMoment")
[docs]class PointScaledNormalEvaluation(Functional): """Implements the evaluation of the normal component of a vector at a point on a facet of codimension 1, where the normal is scaled by the volume of that facet.""" def __init__(self, ref_el, facet_no, pt): self.n = ref_el.compute_scaled_normal(facet_no) sd = ref_el.get_spatial_dimension() shp = (sd,) pt_dict = {pt: [(self.n[i], (i,)) for i in range(sd)]} super().__init__(ref_el, shp, pt_dict, {}, "PointScaledNormalEval")
[docs] def tostr(self): x = list(map(str, list(self.pt_dict.keys())[0])) return "(u.n)(%s)" % (','.join(x),)
[docs]class IntegralMomentOfScaledNormalEvaluation(Functional): r""" \int_F v\cdot n p ds p \in Polynomials :arg ref_el: reference element for which F is a codim-1 entity :arg Q: quadrature rule on the face :arg P_at_qpts: polynomials evaluated at quad points :arg facet: which facet. """ def __init__(self, ref_el, Q, P_at_qpts, facet): n = ref_el.compute_scaled_normal(facet) sd = ref_el.get_spatial_dimension() transform = ref_el.get_entity_transform(sd - 1, facet) pts = tuple(map(lambda p: tuple(transform(p)), Q.get_points())) weights = Q.get_weights() pt_dict = OrderedDict() for pt, wgt, phi in zip(pts, weights, P_at_qpts): pt_dict[pt] = [(wgt*phi*n[i], (i, )) for i in range(sd)] super().__init__(ref_el, (sd, ), pt_dict, {}, "IntegralMomentOfScaledNormalEvaluation")
[docs]class PointwiseInnerProductEvaluation(Functional): """ This is a functional on symmetric 2-tensor fields. Let u be such a field, p be a point, and v,w be vectors. This implements the evaluation v^T u(p) w. Clearly v^iu_{ij}w^j = u_{ij}v^iw^j. Thus the value can be computed from the Frobenius inner product of u with wv^T. This gives the correct weights. """ def __init__(self, ref_el, v, w, p): sd = ref_el.get_spatial_dimension() wvT = numpy.outer(w, v) pt_dict = {p: [(wvT[i][j], (i, j)) for i, j in index_iterator((sd, sd))]} shp = (sd, sd) super().__init__(ref_el, shp, pt_dict, {}, "PointwiseInnerProductEval")
[docs]class TensorBidirectionalMomentInnerProductEvaluation(Functional): r""" This is a functional on symmetric 2-tensor fields. Let u be such a field, f a function tabulated at points, and v,w be vectors. This implements the evaluation \int v^T u(x) w f(x). Clearly v^iu_{ij}w^j = u_{ij}v^iw^j. Thus the value can be computed from the Frobenius inner product of u with wv^T. This gives the correct weights. """ def __init__(self, ref_el, v, w, Q, f_at_qpts, comp_deg): sd = ref_el.get_spatial_dimension() wvT = numpy.outer(w, v) qpts, qwts = Q.get_points(), Q.get_weights() pt_dict = {} for k, pt in enumerate(map(tuple(qpts))): pt_dict[pt] = [] for i, j in index_iterator((sd, sd)): pt_dict[pt].append((qwts[k] * wvT[i][j] * f_at_qpts[i, j, k]), (i, j)) shp = (sd, sd) super().__init__(ref_el, shp, pt_dict, {}, "TensorBidirectionalMomentInnerProductEvaluation")
[docs]class IntegralMomentOfNormalEvaluation(Functional): r""" \int_F v\cdot n p ds p \in Polynomials :arg ref_el: reference element for which F is a codim-1 entity :arg Q: quadrature rule on the face :arg P_at_qpts: polynomials evaluated at quad points :arg facet: which facet. """ def __init__(self, ref_el, Q, P_at_qpts, facet): # scaling on the normal is ok because edge length then weights # the reference element quadrature appropriately n = ref_el.compute_scaled_normal(facet) sd = ref_el.get_spatial_dimension() transform = ref_el.get_entity_transform(sd - 1, facet) pts = tuple(map(lambda p: tuple(transform(p)), Q.get_points())) weights = Q.get_weights() pt_dict = OrderedDict() for pt, wgt, phi in zip(pts, weights, P_at_qpts): pt_dict[pt] = [(wgt*phi*n[i], (i, )) for i in range(sd)] super().__init__(ref_el, (sd, ), pt_dict, {}, "IntegralMomentOfScaledNormalEvaluation")
[docs]class IntegralMomentOfTangentialEvaluation(Functional): r""" \int_F v\cdot n p ds p \in Polynomials :arg ref_el: reference element for which F is a codim-1 entity :arg Q: quadrature rule on the face :arg P_at_qpts: polynomials evaluated at quad points :arg facet: which facet. """ def __init__(self, ref_el, Q, P_at_qpts, facet): # scaling on the tangent is ok because edge length then weights # the reference element quadrature appropriately sd = ref_el.get_spatial_dimension() assert sd == 2 t = ref_el.compute_edge_tangent(facet) transform = ref_el.get_entity_transform(sd - 1, facet) pts = tuple(map(lambda p: tuple(transform(p)), Q.get_points())) weights = Q.get_weights() pt_dict = OrderedDict() for pt, wgt, phi in zip(pts, weights, P_at_qpts): pt_dict[pt] = [(wgt*phi*t[i], (i, )) for i in range(sd)] super().__init__(ref_el, (sd, ), pt_dict, {}, "IntegralMomentOfScaledTangentialEvaluation")
[docs]class IntegralMomentOfNormalNormalEvaluation(Functional): r""" \int_F (n^T tau n) p ds p \in Polynomials :arg ref_el: reference element for which F is a codim-1 entity :arg Q: quadrature rule on the face :arg P_at_qpts: polynomials evaluated at quad points :arg facet: which facet. """ def __init__(self, ref_el, Q, P_at_qpts, facet): # scaling on the normal is ok because edge length then weights # the reference element quadrature appropriately n = ref_el.compute_scaled_normal(facet) sd = ref_el.get_spatial_dimension() transform = ref_el.get_entity_transform(sd - 1, facet) pts = tuple(map(lambda p: tuple(transform(p)), Q.get_points())) weights = Q.get_weights() pt_dict = OrderedDict() for pt, wgt, phi in zip(pts, weights, P_at_qpts): pt_dict[pt] = [(wgt*phi*n[i], (i, )) for i in range(sd)] super().__init__(ref_el, (sd, ), pt_dict, {}, "IntegralMomentOfScaledNormalEvaluation")