# Copyright (C) 2016 Thomas H. Gibson
#
# This file is part of FIAT (https://www.fenicsproject.org)
#
# SPDX-License-Identifier: LGPL-3.0-or-later
import numpy as np
from FIAT.discontinuous_lagrange import DiscontinuousLagrange
from FIAT.dual_set import DualSet
from FIAT.finite_element import FiniteElement
from FIAT.functional import PointEvaluation
from FIAT.polynomial_set import mis
from FIAT.reference_element import (ufc_simplex, POINT,
LINE, QUADRILATERAL,
TRIANGLE, TETRAHEDRON,
TENSORPRODUCT)
from FIAT.tensor_product import TensorProductElement
# Numerical tolerance for facet-entity identifications
epsilon = 1e-10
[docs]
class TraceError(Exception):
"""Exception caused by tabulating a trace element on the interior of a cell,
or the gradient of a trace element."""
def __init__(self, msg):
super(TraceError, self).__init__(msg)
self.msg = msg
[docs]
class HDivTrace(FiniteElement):
"""Class implementing the trace of hdiv elements. This class
is a stand-alone element family that produces a DG-facet field.
This element is what's produced after performing the trace
operation on an existing H(Div) element.
This element is also known as the discontinuous trace field that
arises in several DG formulations.
"""
def __init__(self, ref_el, degree):
"""Constructor for the HDivTrace element.
:arg ref_el: A reference element, which may be a tensor product
cell.
:arg degree: The degree of approximation. If on a tensor product
cell, then provide a tuple of degrees if you want
varying degrees.
"""
sd = ref_el.get_spatial_dimension()
if sd in (0, 1):
raise ValueError("Cannot take the trace of a %d-dim cell." % sd)
# Store the degrees if on a tensor product cell
if ref_el.get_shape() == TENSORPRODUCT:
try:
degree = tuple(degree)
except TypeError:
degree = (degree,) * len(ref_el.cells)
assert len(ref_el.cells) == len(degree), (
"Number of specified degrees must be equal to the number of cells."
)
else:
if ref_el.get_shape() not in [TRIANGLE, TETRAHEDRON, QUADRILATERAL]:
raise NotImplementedError(
"Trace element on a %s not implemented" % type(ref_el)
)
# Cannot have varying degrees for these reference cells
if isinstance(degree, tuple):
raise ValueError("Must have a tensor product cell if providing multiple degrees")
# Initialize entity dofs and construct the DG elements
# for the facets
facet_sd = sd - 1
dg_elements = {}
entity_dofs = {}
topology = ref_el.get_topology()
for top_dim, entities in topology.items():
cell = ref_el.construct_subelement(top_dim)
entity_dofs[top_dim] = {}
# We have a facet entity!
if cell.get_spatial_dimension() == facet_sd:
dg_elements[top_dim] = construct_dg_element(cell, degree)
# Initialize
for entity in entities:
entity_dofs[top_dim][entity] = []
# Compute the dof numbering for all facet entities
# and extract nodes
offset = 0
pts = []
for facet_dim in sorted(dg_elements):
element = dg_elements[facet_dim]
nf = element.space_dimension()
num_facets = len(topology[facet_dim])
for i in range(num_facets):
entity_dofs[facet_dim][i] = list(range(offset, offset + nf))
offset += nf
# Run over nodes and collect the points for point evaluations
for dof in element.dual_basis():
facet_pt, = dof.get_point_dict()
transform = ref_el.get_entity_transform(facet_dim, i)
pts.append(tuple(transform(facet_pt)))
# Setting up dual basis - only point evaluations
nodes = [PointEvaluation(ref_el, pt) for pt in pts]
dual = DualSet(nodes, ref_el, entity_dofs)
# Degree of the element
deg = max([e.degree() for e in dg_elements.values()])
super(HDivTrace, self).__init__(ref_el, dual, order=deg,
formdegree=facet_sd,
mapping="affine")
# Set up facet elements
self.dg_elements = dg_elements
# Degree for quadrature rule
self.polydegree = deg
[docs]
def degree(self):
"""Return the degree of the (embedding) polynomial space."""
return self.polydegree
[docs]
def get_nodal_basis(self):
"""Return the nodal basis, encoded as a PolynomialSet object,
for the finite element."""
raise NotImplementedError("get_nodal_basis not implemented for the trace element.")
[docs]
def get_coeffs(self):
"""Return the expansion coefficients for the basis of the
finite element."""
raise NotImplementedError("get_coeffs not implemented for the trace element.")
[docs]
def tabulate(self, order, points, entity=None):
"""Return tabulated values of derivatives up to a 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``,
tabulated values are computed by geometrically
approximating which facet the points are on.
.. note ::
Performing illegal tabulations on this element will result in either
a tabulation table of `numpy.nan` arrays (`entity=None` case), or
insertions of the `TraceError` exception class. This is due to the
fact that performing cell-wise tabulations, or asking for any order
of derivative evaluations, are not mathematically well-defined.
"""
sd = self.ref_el.get_spatial_dimension()
facet_sd = sd - 1
# Initializing dictionary with zeros
phivals = {}
for i in range(order + 1):
alphas = mis(sd, i)
for alpha in alphas:
phivals[alpha] = np.zeros(shape=(self.space_dimension(), len(points)))
evalkey = (0,) * sd
# If entity is None, identify facet using numerical tolerance and
# return the tabulated values
if entity is None:
# NOTE: Numerical approximation of the facet id is currently only
# implemented for simplex reference cells.
if self.ref_el.get_shape() not in [TRIANGLE, TETRAHEDRON]:
raise NotImplementedError(
"Tabulating this element on a %s cell without providing "
"an entity is not currently supported." % type(self.ref_el)
)
# Attempt to identify which facet (if any) the given points are on
vertices = self.ref_el.vertices
coordinates = barycentric_coordinates(points, vertices)
unique_facet, success = extract_unique_facet(coordinates)
# If not successful, return NaNs
if not success:
for key in phivals:
phivals[key] = np.full(shape=(self.space_dimension(), len(points)), fill_value=np.nan)
return phivals
# Otherwise, extract non-zero values and insertion indices
else:
# Map points to the reference facet
new_points = map_to_reference_facet(points, vertices, unique_facet)
# Retrieve values by tabulating the DG element
element = self.dg_elements[facet_sd]
nf = element.space_dimension()
nonzerovals, = element.tabulate(order, new_points).values()
indices = slice(nf * unique_facet, nf * (unique_facet + 1))
else:
entity_dim, _ = entity
# If the user is directly specifying cell-wise tabulation, return
# TraceErrors in dict for appropriate handling in the form compiler
if entity_dim not in self.dg_elements:
for key in phivals:
msg = "The HDivTrace element can only be tabulated on facets."
phivals[key] = TraceError(msg)
return phivals
else:
# Retrieve function evaluations (order = 0 case)
offset = 0
for facet_dim in sorted(self.dg_elements):
element = self.dg_elements[facet_dim]
nf = element.space_dimension()
num_facets = len(self.ref_el.get_topology()[facet_dim])
# Loop over the number of facets until we find a facet
# with matching dimension and id
for i in range(num_facets):
# Found it! Grab insertion indices
if (facet_dim, i) == entity:
nonzerovals, = element.tabulate(0, points).values()
indices = slice(offset, offset + nf)
offset += nf
# If asking for gradient evaluations, insert TraceError in
# gradient slots
if order > 0:
msg = "Gradients on trace elements are not well-defined."
for key in phivals:
if key != evalkey:
phivals[key] = TraceError(msg)
# Insert non-zero values in appropriate place
phivals[evalkey][indices, :] = nonzerovals
return phivals
[docs]
def value_shape(self):
"""Return the value shape of the finite element functions."""
return ()
[docs]
def dmats(self):
"""Return dmats: expansion coefficients for basis function
derivatives."""
raise NotImplementedError("dmats not implemented for the trace element.")
[docs]
def get_num_members(self, arg):
"""Return number of members of the expansion set."""
raise NotImplementedError("get_num_members not implemented for the trace element.")
[docs]
@staticmethod
def is_nodal():
return True
[docs]
def construct_dg_element(ref_el, degree):
"""Constructs a discontinuous galerkin element of a given degree
on a particular reference cell.
"""
if ref_el.get_shape() in [LINE, TRIANGLE]:
dg_element = DiscontinuousLagrange(ref_el, degree)
# Quadrilateral facets could be on a FiredrakeQuadrilateral.
# In this case, we treat this as an interval x interval cell:
elif ref_el.get_shape() == QUADRILATERAL:
dg_a = DiscontinuousLagrange(ufc_simplex(1), degree)
dg_b = DiscontinuousLagrange(ufc_simplex(1), degree)
dg_element = TensorProductElement(dg_a, dg_b)
# This handles the more general case for facets:
elif ref_el.get_shape() == TENSORPRODUCT:
assert len(degree) == len(ref_el.cells), (
"Must provide the same number of degrees as the number "
"of cells that make up the tensor product cell."
)
sub_elements = [construct_dg_element(c, d)
for c, d in zip(ref_el.cells, degree)
if c.get_shape() != POINT]
if len(sub_elements) > 1:
dg_element = TensorProductElement(*sub_elements)
else:
dg_element, = sub_elements
else:
raise NotImplementedError(
"Reference cells of type %s not currently supported" % type(ref_el)
)
return dg_element
# The following functions are credited to Marie E. Rognes:
[docs]
def barycentric_coordinates(points, vertices):
"""Computes the barycentric coordinates for a set of points relative to a
simplex defined by a set of vertices.
:arg points: A set of points.
:arg vertices: A set of vertices that define the simplex.
"""
# Form mapping matrix
T = (np.asarray(vertices[:-1]) - vertices[-1]).T
invT = np.linalg.inv(T)
points = np.asarray(points)
bary = np.einsum("ij,kj->ki", invT, (points - vertices[-1]))
last = (1 - bary.sum(axis=1))
return np.concatenate([bary, last[..., np.newaxis]], axis=1)
[docs]
def map_from_reference_facet(point, vertices):
"""Evaluates the physical coordinate of a point using barycentric
coordinates.
:arg point: The reference points to be mapped to the facet.
:arg vertices: The vertices defining the physical element.
"""
# Compute the barycentric coordinates of the point relative to the reference facet
reference_simplex = ufc_simplex(len(vertices) - 1)
reference_vertices = reference_simplex.get_vertices()
coords = barycentric_coordinates([point, ], reference_vertices)[0]
# Evaluates the physical coordinate of the point using barycentric coordinates
point = sum(vertices[j] * coords[j] for j in range(len(coords)))
return tuple(point)
[docs]
def map_to_reference_facet(points, vertices, facet):
"""Given a set of points and vertices describing a facet of a simplex in n-dimensional
coordinates (where the points lie on the facet), map the points to the reference simplex
of dimension (n-1).
:arg points: A set of points in n-D.
:arg vertices: A set of vertices describing a facet of a simplex in n-D.
:arg facet: Integer representing the facet number.
"""
# Compute the barycentric coordinates of the points with respect to the
# full physical simplex
all_coords = barycentric_coordinates(points, vertices)
# Extract vertices of the reference facet
reference_facet_simplex = ufc_simplex(len(vertices) - 2)
reference_vertices = reference_facet_simplex.get_vertices()
reference_points = []
for (i, coords) in enumerate(all_coords):
# Extract the correct subset of barycentric coordinates since we know
# which facet we are on
new_coords = [coords[j] for j in range(len(coords)) if j != facet]
# Evaluate the reference coordinate of a point in barycentric coordinates
reference_pt = sum(np.asarray(reference_vertices[j]) * new_coords[j]
for j in range(len(new_coords)))
reference_points += [reference_pt]
return reference_points