# Copyright (C) 2008 Robert C. Kirby (Texas Tech University)
# Copyright (C) 2013 Andrew T. T. McRae
# Modified by Thomas H. Gibson, 2016
#
# This file is part of FIAT (https://www.fenicsproject.org)
#
# SPDX-License-Identifier: LGPL-3.0-or-later
import numpy
from FIAT.finite_element import FiniteElement
from FIAT.reference_element import TensorProductCell, UFCQuadrilateral, UFCHexahedron, flatten_entities, compute_unflattening_map
from FIAT.dual_set import DualSet
from FIAT.polynomial_set import mis
from FIAT import dual_set
from FIAT import functional
def _first_point(node):
return tuple(node.get_point_dict().keys())[0]
def _first_point_pair(node):
return tuple(node.get_point_dict().items())[0]
[docs]class TensorProductElement(FiniteElement):
"""Class implementing a finite element that is the tensor product
of two existing finite elements."""
def __init__(self, A, B):
# set up simple things
order = min(A.get_order(), B.get_order())
if A.get_formdegree() is None or B.get_formdegree() is None:
formdegree = None
else:
formdegree = A.get_formdegree() + B.get_formdegree()
# set up reference element
ref_el = TensorProductCell(A.get_reference_element(),
B.get_reference_element())
if A.mapping()[0] != "affine" and B.mapping()[0] == "affine":
mapping = A.mapping()[0]
elif B.mapping()[0] != "affine" and A.mapping()[0] == "affine":
mapping = B.mapping()[0]
elif A.mapping()[0] == "affine" and B.mapping()[0] == "affine":
mapping = "affine"
else:
raise ValueError("check tensor product mappings - at least one must be affine")
# set up entity_ids
Adofs = A.entity_dofs()
Bdofs = B.entity_dofs()
Bsdim = B.space_dimension()
entity_ids = {}
for curAdim in Adofs:
for curBdim in Bdofs:
entity_ids[(curAdim, curBdim)] = {}
dim_cur = 0
for entityA in Adofs[curAdim]:
for entityB in Bdofs[curBdim]:
entity_ids[(curAdim, curBdim)][dim_cur] = \
[x*Bsdim + y for x in Adofs[curAdim][entityA]
for y in Bdofs[curBdim][entityB]]
dim_cur += 1
# set up dual basis
Anodes = A.dual_basis()
Bnodes = B.dual_basis()
# build the dual set by inspecting the current dual
# sets item by item.
# Currently supported cases:
# PointEval x PointEval = PointEval [scalar x scalar = scalar]
# PointScaledNormalEval x PointEval = PointScaledNormalEval [vector x scalar = vector]
# ComponentPointEvaluation x PointEval [vector x scalar = vector]
nodes = []
for Anode in Anodes:
if isinstance(Anode, functional.PointEvaluation):
for Bnode in Bnodes:
if isinstance(Bnode, functional.PointEvaluation):
# case: PointEval x PointEval
# the PointEval functional just requires the
# coordinates. these are currently stored as
# the key of a one-item dictionary. we retrieve
# these by calling get_point_dict(), and
# use the concatenation to make a new PointEval
nodes.append(functional.PointEvaluation(ref_el, _first_point(Anode) + _first_point(Bnode)))
elif isinstance(Bnode, functional.IntegralMoment):
# dummy functional for product with integral moments
nodes.append(functional.Functional(None, None, None,
{}, "Undefined"))
elif isinstance(Bnode, functional.PointDerivative):
# dummy functional for product with point derivative
nodes.append(functional.Functional(None, None, None,
{}, "Undefined"))
else:
raise NotImplementedError("unsupported functional type")
elif isinstance(Anode, functional.PointScaledNormalEvaluation):
for Bnode in Bnodes:
if isinstance(Bnode, functional.PointEvaluation):
# case: PointScaledNormalEval x PointEval
# this could be wrong if the second shape
# has spatial dimension >1, since we are not
# explicitly scaling by facet size
if len(_first_point(Bnode)) > 1:
# TODO: support this case one day
raise NotImplementedError("PointScaledNormalEval x PointEval is not yet supported if the second shape has dimension > 1")
# We cannot make a new functional.PSNEval in
# the natural way, since it tries to compute
# the normal vector by itself.
# Instead, we create things manually, and
# call Functional() with these arguments
sd = ref_el.get_spatial_dimension()
# The pt_dict is a one-item dictionary containing
# the details of the functional.
# The key is the spatial coordinate, which
# is just a concatenation of the two parts.
# The value is a list of tuples, representing
# the normal vector (scaled by the volume of
# the facet) at that point.
# Each tuple looks like (foo, (i,)); the i'th
# component of the scaled normal is foo.
# The following line is only valid when the second
# shape has spatial dimension 1 (enforced above)
Apoint, Avalue = _first_point_pair(Anode)
pt_dict = {Apoint + _first_point(Bnode): Avalue + [(0.0, (len(Apoint),))]}
# The following line should be used in the
# general case
# pt_dict = {Anode.get_point_dict().keys()[0] + Bnode.get_point_dict().keys()[0]: Anode.get_point_dict().values()[0] + [(0.0, (ii,)) for ii in range(len(Anode.get_point_dict().keys()[0]), len(Anode.get_point_dict().keys()[0]) + len(Bnode.get_point_dict().keys()[0]))]}
# THE FOLLOWING IS PROBABLY CORRECT BUT UNTESTED
shp = (sd,)
nodes.append(functional.Functional(ref_el, shp, pt_dict, {}, "PointScaledNormalEval"))
else:
raise NotImplementedError("unsupported functional type")
elif isinstance(Anode, functional.PointEdgeTangentEvaluation):
for Bnode in Bnodes:
if isinstance(Bnode, functional.PointEvaluation):
# case: PointEdgeTangentEval x PointEval
# this is very similar to the case above, so comments omitted
if len(_first_point(Bnode)) > 1:
raise NotImplementedError("PointEdgeTangentEval x PointEval is not yet supported if the second shape has dimension > 1")
sd = ref_el.get_spatial_dimension()
Apoint, Avalue = _first_point_pair(Anode)
pt_dict = {Apoint + _first_point(Bnode): Avalue + [(0.0, (len(Apoint),))]}
# THE FOLLOWING IS PROBABLY CORRECT BUT UNTESTED
shp = (sd,)
nodes.append(functional.Functional(ref_el, shp, pt_dict, {}, "PointEdgeTangent"))
else:
raise NotImplementedError("unsupported functional type")
elif isinstance(Anode, functional.ComponentPointEvaluation):
for Bnode in Bnodes:
if isinstance(Bnode, functional.PointEvaluation):
# case: ComponentPointEval x PointEval
# the CptPointEval functional requires the component
# and the coordinates. very similar to PE x PE case.
sd = ref_el.get_spatial_dimension()
nodes.append(functional.ComponentPointEvaluation(ref_el, Anode.comp, (sd,), _first_point(Anode) + _first_point(Bnode)))
else:
raise NotImplementedError("unsupported functional type")
elif isinstance(Anode, functional.FrobeniusIntegralMoment):
for Bnode in Bnodes:
if isinstance(Bnode, functional.PointEvaluation):
# case: FroIntMom x PointEval
sd = ref_el.get_spatial_dimension()
pt_dict = {}
pt_old = Anode.get_point_dict()
for pt in pt_old:
pt_dict[pt+_first_point(Bnode)] = pt_old[pt] + [(0.0, sd-1)]
# THE FOLLOWING IS PROBABLY CORRECT BUT UNTESTED
shp = (sd,)
nodes.append(functional.Functional(ref_el, shp, pt_dict, {}, "FrobeniusIntegralMoment"))
else:
raise NotImplementedError("unsupported functional type")
elif isinstance(Anode, functional.IntegralMoment):
for Bnode in Bnodes:
if isinstance(Bnode, functional.PointEvaluation):
# case: IntMom x PointEval
sd = ref_el.get_spatial_dimension()
pt_dict = {}
pt_old = Anode.get_point_dict()
for pt in pt_old:
pt_dict[pt+_first_point(Bnode)] = pt_old[pt]
# THE FOLLOWING IS PROBABLY CORRECT BUT UNTESTED
shp = (sd,)
nodes.append(functional.Functional(ref_el, shp, pt_dict, {}, "IntegralMoment"))
else:
raise NotImplementedError("unsupported functional type")
elif isinstance(Anode, functional.Functional):
# this should catch everything else
for Bnode in Bnodes:
nodes.append(functional.Functional(None, None, None, {}, "Undefined"))
else:
raise NotImplementedError("unsupported functional type")
dual = dual_set.DualSet(nodes, ref_el, entity_ids)
super(TensorProductElement, self).__init__(ref_el, dual, order, formdegree, mapping)
# Set up constituent elements
self.A = A
self.B = B
# degree for quadrature rule
self.polydegree = max(A.degree(), B.degree())
[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")
[docs] def get_coeffs(self):
"""Return the expansion coefficients for the basis of the
finite element."""
raise NotImplementedError("get_coeffs not implemented")
[docs] def tabulate(self, order, points, entity=None):
"""Return tabulated values of derivatives up to given order of
basis functions at given points."""
if entity is None:
entity = (self.ref_el.get_dimension(), 0)
entity_dim, entity_id = entity
shape = tuple(len(c.get_topology()[d])
for c, d in zip(self.ref_el.cells, entity_dim))
idA, idB = numpy.unravel_index(entity_id, shape)
# Factor the entity argument to get entities of the component elements
entityA_dim, entityB_dim = entity_dim
entityA = (entityA_dim, idA)
entityB = (entityB_dim, idB)
pointsAdim, pointsBdim = [c.get_spatial_dimension()
for c in self.ref_el.construct_subelement(entity_dim).cells]
pointsA = [point[:pointsAdim] for point in points]
pointsB = [point[pointsAdim:pointsAdim + pointsBdim] for point in points]
Asdim = self.A.ref_el.get_spatial_dimension()
Bsdim = self.B.ref_el.get_spatial_dimension()
# Note that for entities other than cells, the following
# tabulations are already appropriately zero-padded so no
# additional zero padding is required.
Atab = self.A.tabulate(order, pointsA, entityA)
Btab = self.B.tabulate(order, pointsB, entityB)
npoints = len(points)
# allow 2 scalar-valued FE spaces, or 1 scalar-valued,
# 1 vector-valued. Combining 2 vector-valued spaces
# into a tensor-valued space via an outer-product
# seems to be a sensible general option, but I don't
# know how to handle the nestedness of the arrays
# if someone then tries to make a new "tensor finite
# element" where one component is already a
# tensor-valued space!
A_valuedim = len(self.A.value_shape()) # scalar: 0, vector: 1
B_valuedim = len(self.B.value_shape()) # scalar: 0, vector: 1
if A_valuedim + B_valuedim > 1:
raise NotImplementedError("tabulate does not support two vector-valued inputs")
result = {}
for i in range(order + 1):
alphas = mis(Asdim+Bsdim, i) # thanks, Rob!
for alpha in alphas:
if A_valuedim == 0 and B_valuedim == 0:
# for each point, get outer product of (A's basis
# functions f1, f2, ... evaluated at that point)
# with (B's basis functions g1, g2, ... evaluated
# at that point). This gives temp[point][f_i][g_j].
# Flatten this, so bfs are
# in the order f1g1, f1g2, ..., f2g1, f2g2, ...
# which is compatible with the entity_dofs order.
# We now have temp[point][full basis function]
# Transpose this to get temp[bf][point],
# and we are done.
temp = numpy.array([numpy.outer(
Atab[alpha[0:Asdim]][..., j],
Btab[alpha[Asdim:Asdim+Bsdim]][..., j])
.ravel() for j in range(npoints)])
result[alpha] = temp.transpose()
elif A_valuedim == 1 and B_valuedim == 0:
# similar to above, except A's basis functions
# are now vector-valued. numpy.outer flattens the
# array, so it's like taking the OP of
# f1_x, f1_y, f2_x, f2_y, ... with g1, g2, ...
# this gives us
# temp[point][f1x, f1y, f2x, f2y, ...][g_j].
# reshape once to get temp[point][f_i][x/y][g_j]
# transpose to get temp[point][x/y][f_i][g_j]
# reshape to flatten the last two indices, this
# gives us temp[point][x/y][full bf_i]
# finally, transpose the first and last indices
# to get temp[bf_i][x/y][point], and we are done.
temp = numpy.array([numpy.outer(
Atab[alpha[0:Asdim]][..., j],
Btab[alpha[Asdim:Asdim+Bsdim]][..., j])
for j in range(npoints)])
assert temp.shape[1] % 2 == 0
temp2 = temp.reshape((temp.shape[0],
temp.shape[1]//2,
2,
temp.shape[2]))\
.transpose(0, 2, 1, 3)\
.reshape((temp.shape[0], 2, -1))\
.transpose(2, 1, 0)
result[alpha] = temp2
elif A_valuedim == 0 and B_valuedim == 1:
# as above, with B's functions now vector-valued.
# we now do... [numpy.outer ... for ...] gives
# temp[point][f_i][g1x,g1y,g2x,g2y,...].
# reshape to temp[point][f_i][g_j][x/y]
# flatten middle: temp[point][full bf_i][x/y]
# transpose to temp[bf_i][x/y][point]
temp = numpy.array([numpy.outer(
Atab[alpha[0:Asdim]][..., j],
Btab[alpha[Asdim:Asdim+Bsdim]][..., j])
for j in range(len(Atab[alpha[0:Asdim]][0]))])
assert temp.shape[2] % 2 == 0
temp2 = temp.reshape((temp.shape[0], temp.shape[1],
temp.shape[2]//2, 2))\
.reshape((temp.shape[0], -1, 2))\
.transpose(1, 2, 0)
result[alpha] = temp2
return result
[docs] def value_shape(self):
"""Return the value shape of the finite element functions."""
if len(self.A.value_shape()) == 0 and len(self.B.value_shape()) == 0:
return ()
elif len(self.A.value_shape()) == 1 and len(self.B.value_shape()) == 0:
return (self.A.value_shape()[0],)
elif len(self.A.value_shape()) == 0 and len(self.B.value_shape()) == 1:
return (self.B.value_shape()[0],)
else:
raise NotImplementedError("value_shape not implemented")
[docs] def dmats(self):
"""Return dmats: expansion coefficients for basis function
derivatives."""
raise NotImplementedError("dmats not implemented")
[docs] def get_num_members(self, arg):
"""Return number of members of the expansion set."""
raise NotImplementedError("get_num_members not implemented")
[docs] def is_nodal(self):
# This element is nodal iff all factor elements are nodal.
return all([self.A.is_nodal(), self.B.is_nodal()])
[docs]class FlattenedDimensions(FiniteElement):
"""A wrapper class that flattens entity dimensions of a FIAT element defined
on a TensorProductCell to one with quadrilateral/hexahedron entities.
TensorProductCell has dimension defined as a tuple of factor element dimensions
(i, j) in 2D and (i, j, k) in 3D.
Flattened dimension is a sum of the tuple elements."""
def __init__(self, element):
nodes = element.dual.nodes
dim = element.ref_el.get_spatial_dimension()
if dim == 2:
ref_el = UFCQuadrilateral()
elif dim == 3:
ref_el = UFCHexahedron()
else:
raise ValueError("Illegal element dimension %s" % dim)
entity_ids = element.dual.entity_ids
flat_entity_ids = flatten_entities(entity_ids)
dual = DualSet(nodes, ref_el, flat_entity_ids)
super(FlattenedDimensions, self).__init__(ref_el, dual, element.get_order(), element.get_formdegree(), element._mapping)
self.element = element
# Construct unflattening map for passing correct values to tabulate()
self.unflattening_map = compute_unflattening_map(self.element.ref_el.get_topology())
[docs] def degree(self):
"""Return the degree of the (embedding) polynomial space."""
return self.element.degree()
[docs] def tabulate(self, order, points, entity=None):
"""Return tabulated values of derivatives up to given order of
basis functions at given points."""
if entity is None:
entity = (self.get_reference_element().get_spatial_dimension(), 0)
# Entity is provided in flattened form (d, i)
# Appropriate product entity is taken from the unflattening_map dict
entity_dim, entity_id = entity
product_entity = self.unflattening_map[(entity_dim, entity_id)]
return self.element.tabulate(order, points, product_entity)
[docs] def value_shape(self):
"""Return the value shape of the finite element functions."""
return self.element.value_shape()
[docs] def get_nodal_basis(self):
"""Return the nodal basis, encoded as a PolynomialSet object,
for the finite element."""
raise self.element.get_nodal_basis()
[docs] def get_coeffs(self):
"""Return the expansion coefficients for the basis of the
finite element."""
raise self.element.get_coeffs()
[docs] def dmats(self):
"""Return dmats: expansion coefficients for basis function
derivatives."""
raise self.element.dmats()
[docs] def get_num_members(self, arg):
"""Return number of members of the expansion set."""
raise self.element.get_num_members(arg)
[docs] def is_nodal(self):
# This element is nodal iff unflattened element is nodal.
return self.element.is_nodal()