# Copyright (C) 2008 Robert C. Kirby (Texas Tech University)
#
# This file is part of FIAT (https://www.fenicsproject.org)
#
# SPDX-License-Identifier: LGPL-3.0-or-later
#
# Modified by David A. Ham (david.ham@imperial.ac.uk), 2014
# Modified by Lizao Li (lzlarryli@gmail.com), 2016
"""
Abstract class and particular implementations of finite element
reference simplex geometry/topology.
Provides an abstract base class and particular implementations for the
reference simplex geometry and topology.
The rest of FIAT is abstracted over this module so that different
reference element geometry (e.g. a vertex at (0,0) versus at (-1,-1))
and orderings of entities have a single point of entry.
Currently implemented are UFC and Default Line, Triangle and Tetrahedron.
"""
from itertools import chain, product, count
from functools import reduce
from collections import defaultdict
import operator
from math import factorial
import numpy
POINT = 0
LINE = 1
TRIANGLE = 2
TETRAHEDRON = 3
QUADRILATERAL = 11
HEXAHEDRON = 111
TENSORPRODUCT = 99
[docs]
def lattice_iter(start, finish, depth):
"""Generator iterating over the depth-dimensional lattice of
integers between start and (finish-1). This works on simplices in
1d, 2d, 3d, and beyond"""
if depth == 0:
return
elif depth == 1:
for ii in range(start, finish):
yield [ii]
else:
for ii in range(start, finish):
for jj in lattice_iter(start, finish - ii, depth - 1):
yield jj + [ii]
[docs]
def make_lattice(verts, n, interior=0):
"""Constructs a lattice of points on the simplex defined by verts.
For example, the 1:st order lattice will be just the vertices.
The optional argument interior specifies how many points from
the boundary to omit. For example, on a line with n = 2,
and interior = 0, this function will return the vertices and
midpoint, but with interior = 1, it will only return the
midpoint."""
vs = numpy.array(verts)
hs = (vs - vs[0])[1:, :] / n
m = hs.shape[0]
result = [tuple(vs[0] + numpy.array(indices).dot(hs))
for indices in lattice_iter(interior, n + 1 - interior, m)]
return result
[docs]
def linalg_subspace_intersection(A, B):
"""Computes the intersection of the subspaces spanned by the
columns of 2-dimensional arrays A,B using the algorithm found in
Golub and van Loan (3rd ed) p. 604. A should be in
R^{m,p} and B should be in R^{m,q}. Returns an orthonormal basis
for the intersection of the spaces, stored in the columns of
the result."""
# check that vectors are in same space
if A.shape[0] != B.shape[0]:
raise Exception("Dimension error")
# A,B are matrices of column vectors
# compute the intersection of span(A) and span(B)
# Compute the principal vectors/angles between the subspaces, G&vL
# p.604
(qa, _ra) = numpy.linalg.qr(A)
(qb, _rb) = numpy.linalg.qr(B)
C = numpy.dot(numpy.transpose(qa), qb)
(y, c, _zt) = numpy.linalg.svd(C)
U = numpy.dot(qa, y)
rank_c = len([s for s in c if numpy.abs(1.0 - s) < 1.e-10])
return U[:, :rank_c]
[docs]
class Cell(object):
"""Abstract class for a reference cell. Provides accessors for
geometry (vertex coordinates) as well as topology (orderings of
vertices that make up edges, facecs, etc."""
def __init__(self, shape, vertices, topology):
"""The constructor takes a shape code, the physical vertices expressed
as a list of tuples of numbers, and the topology of a cell.
The topology is stored as a dictionary of dictionaries t[i][j]
where i is the dimension and j is the index of the facet of
that dimension. The result is a list of the vertices
comprising the facet."""
self.shape = shape
self.vertices = vertices
self.topology = topology
# Given the topology, work out for each entity in the cell,
# which other entities it contains.
self.sub_entities = {}
for dim, entities in topology.items():
self.sub_entities[dim] = {}
for e, v in entities.items():
vertices = frozenset(v)
sub_entities = []
for dim_, entities_ in topology.items():
for e_, vertices_ in entities_.items():
if vertices.issuperset(vertices_):
sub_entities.append((dim_, e_))
# Sort for the sake of determinism and by UFC conventions
self.sub_entities[dim][e] = sorted(sub_entities)
# Build connectivity dictionary for easier queries
self.connectivity = {}
for dim0, sub_entities in self.sub_entities.items():
# Skip tensor product entities
# TODO: Can we do something better?
if isinstance(dim0, tuple):
continue
for entity, sub_sub_entities in sorted(sub_entities.items()):
for dim1 in range(dim0+1):
d01_entities = filter(lambda x: x[0] == dim1, sub_sub_entities)
d01_entities = tuple(x[1] for x in d01_entities)
self.connectivity.setdefault((dim0, dim1), []).append(d01_entities)
def _key(self):
"""Hashable object key data (excluding type)."""
# Default: only type matters
return None
def __eq__(self, other):
return type(self) == type(other) and self._key() == other._key()
def __ne__(self, other):
return not self.__eq__(other)
def __hash__(self):
return hash((type(self), self._key()))
[docs]
def get_shape(self):
"""Returns the code for the element's shape."""
return self.shape
[docs]
def get_vertices(self):
"""Returns an iterable of the element's vertices, each stored as a
tuple."""
return self.vertices
[docs]
def get_spatial_dimension(self):
"""Returns the spatial dimension in which the element lives."""
return len(self.vertices[0])
[docs]
def get_topology(self):
"""Returns a dictionary encoding the topology of the element.
The dictionary's keys are the spatial dimensions (0, 1, ...)
and each value is a dictionary mapping."""
return self.topology
[docs]
def get_connectivity(self):
"""Returns a dictionary encoding the connectivity of the element.
The dictionary's keys are the spatial dimensions pairs ((1, 0),
(2, 0), (2, 1), ...) and each value is a list with entities
of second dimension ordered by local dim0-dim1 numbering."""
return self.connectivity
[docs]
def get_vertices_of_subcomplex(self, t):
"""Returns the tuple of vertex coordinates associated with the labels
contained in the iterable t."""
return tuple([self.vertices[ti] for ti in t])
[docs]
def get_dimension(self):
"""Returns the subelement dimension of the cell. For tensor
product cells, this a tuple of dimensions for each cell in the
product. For all other cells, this is the same as the spatial
dimension."""
raise NotImplementedError("Should be implemented in a subclass.")
[docs]
def construct_subelement(self, dimension):
"""Constructs the reference element of a cell subentity
specified by subelement dimension.
:arg dimension: `tuple` for tensor product cells, `int` otherwise
"""
raise NotImplementedError("Should be implemented in a subclass.")
[docs]
class Simplex(Cell):
"""Abstract class for a reference simplex."""
[docs]
def compute_normal(self, facet_i):
"""Returns the unit normal vector to facet i of codimension 1."""
# Interval case
if self.get_shape() == LINE:
verts = numpy.asarray(self.vertices)
v_i, = self.get_topology()[0][facet_i]
n = verts[v_i] - verts[[1, 0][v_i]]
return n / numpy.linalg.norm(n)
# first, let's compute the span of the simplex
# This is trivial if we have a d-simplex in R^d.
# Not so otherwise.
vert_vecs = [numpy.array(v)
for v in self.vertices]
vert_vecs_foo = numpy.array([vert_vecs[i] - vert_vecs[0]
for i in range(1, len(vert_vecs))])
(u, s, vt) = numpy.linalg.svd(vert_vecs_foo)
rank = len([si for si in s if si > 1.e-10])
# this is the set of vectors that span the simplex
spanu = u[:, :rank]
t = self.get_topology()
sd = self.get_spatial_dimension()
vert_coords_of_facet = \
self.get_vertices_of_subcomplex(t[sd-1][facet_i])
# now I find everything normal to the facet.
vcf = [numpy.array(foo)
for foo in vert_coords_of_facet]
facet_span = numpy.array([vcf[i] - vcf[0]
for i in range(1, len(vcf))])
(uf, sf, vft) = numpy.linalg.svd(facet_span)
# now get the null space from vft
rankfacet = len([si for si in sf if si > 1.e-10])
facet_normal_space = numpy.transpose(vft[rankfacet:, :])
# now, I have to compute the intersection of
# facet_span with facet_normal_space
foo = linalg_subspace_intersection(facet_normal_space, spanu)
num_cols = foo.shape[1]
if num_cols != 1:
raise Exception("barf in normal computation")
# now need to get the correct sign
# get a vector in the direction
nfoo = foo[:, 0]
# what is the vertex not in the facet?
verts_set = set(t[sd][0])
verts_facet = set(t[sd - 1][facet_i])
verts_diff = verts_set.difference(verts_facet)
if len(verts_diff) != 1:
raise Exception("barf in normal computation: getting sign")
vert_off = verts_diff.pop()
vert_on = verts_facet.pop()
# get a vector from the off vertex to the facet
v_to_facet = numpy.array(self.vertices[vert_on]) \
- numpy.array(self.vertices[vert_off])
if numpy.dot(v_to_facet, nfoo) > 0.0:
return nfoo
else:
return -nfoo
[docs]
def compute_tangents(self, dim, i):
"""Computes tangents in any dimension based on differences
between vertices and the first vertex of the i:th facet
of dimension dim. Returns a (possibly empty) list.
These tangents are *NOT* normalized to have unit length."""
t = self.get_topology()
vs = list(map(numpy.array, self.get_vertices_of_subcomplex(t[dim][i])))
ts = [v - vs[0] for v in vs[1:]]
return ts
[docs]
def compute_normalized_tangents(self, dim, i):
"""Computes tangents in any dimension based on differences
between vertices and the first vertex of the i:th facet
of dimension dim. Returns a (possibly empty) list.
These tangents are normalized to have unit length."""
ts = self.compute_tangents(dim, i)
return [t / numpy.linalg.norm(t) for t in ts]
[docs]
def compute_edge_tangent(self, edge_i):
"""Computes the nonnormalized tangent to a 1-dimensional facet.
returns a single vector."""
t = self.get_topology()
(v0, v1) = self.get_vertices_of_subcomplex(t[1][edge_i])
return numpy.array(v1) - numpy.array(v0)
[docs]
def compute_normalized_edge_tangent(self, edge_i):
"""Computes the unit tangent vector to a 1-dimensional facet"""
v = self.compute_edge_tangent(edge_i)
return v / numpy.linalg.norm(v)
[docs]
def compute_face_tangents(self, face_i):
"""Computes the two tangents to a face. Only implemented
for a tetrahedron."""
if self.get_spatial_dimension() != 3:
raise Exception("can't get face tangents yet")
t = self.get_topology()
(v0, v1, v2) = list(map(numpy.array,
self.get_vertices_of_subcomplex(t[2][face_i])))
return (v1 - v0, v2 - v0)
[docs]
def compute_face_edge_tangents(self, dim, entity_id):
"""Computes all the edge tangents of any k-face with k>=1.
The result is a array of binom(dim+1,2) vectors.
This agrees with `compute_edge_tangent` when dim=1.
"""
vert_ids = self.get_topology()[dim][entity_id]
vert_coords = [numpy.array(x)
for x in self.get_vertices_of_subcomplex(vert_ids)]
edge_ts = []
for source in range(dim):
for dest in range(source + 1, dim + 1):
edge_ts.append(vert_coords[dest] - vert_coords[source])
return edge_ts
[docs]
def make_points(self, dim, entity_id, order):
"""Constructs a lattice of points on the entity_id:th
facet of dimension dim. Order indicates how many points to
include in each direction."""
if dim == 0:
return (self.get_vertices()[entity_id], )
elif 0 < dim < self.get_spatial_dimension():
entity_verts = \
self.get_vertices_of_subcomplex(
self.get_topology()[dim][entity_id])
return make_lattice(entity_verts, order, 1)
elif dim == self.get_spatial_dimension():
return make_lattice(self.get_vertices(), order, 1)
else:
raise ValueError("illegal dimension")
[docs]
def volume(self):
"""Computes the volume of the simplex in the appropriate
dimensional measure."""
return volume(self.get_vertices())
[docs]
def volume_of_subcomplex(self, dim, facet_no):
vids = self.topology[dim][facet_no]
return volume(self.get_vertices_of_subcomplex(vids))
[docs]
def compute_scaled_normal(self, facet_i):
"""Returns the unit normal to facet_i of scaled by the
volume of that facet."""
dim = self.get_spatial_dimension()
v = self.volume_of_subcomplex(dim - 1, facet_i)
return self.compute_normal(facet_i) * v
[docs]
def compute_reference_normal(self, facet_dim, facet_i):
"""Returns the unit normal in infinity norm to facet_i."""
assert facet_dim == self.get_spatial_dimension() - 1
n = Simplex.compute_normal(self, facet_i) # skip UFC overrides
return n / numpy.linalg.norm(n, numpy.inf)
[docs]
def get_dimension(self):
"""Returns the subelement dimension of the cell. Same as the
spatial dimension."""
return self.get_spatial_dimension()
# Backwards compatible name
ReferenceElement = Simplex
[docs]
class UFCSimplex(Simplex):
[docs]
def get_facet_element(self):
dimension = self.get_spatial_dimension()
return self.construct_subelement(dimension - 1)
[docs]
def construct_subelement(self, dimension):
"""Constructs the reference element of a cell subentity
specified by subelement dimension.
:arg dimension: subentity dimension (integer)
"""
return ufc_simplex(dimension)
[docs]
def contains_point(self, point, epsilon=0):
"""Checks if reference cell contains given point
(with numerical tolerance)."""
result = (sum(point) - epsilon <= 1)
for c in point:
result &= (c + epsilon >= 0)
return result
[docs]
class Point(Simplex):
"""This is the reference point."""
def __init__(self):
verts = ((),)
topology = {0: {0: (0,)}}
super(Point, self).__init__(POINT, verts, topology)
[docs]
def construct_subelement(self, dimension):
"""Constructs the reference element of a cell subentity
specified by subelement dimension.
:arg dimension: subentity dimension (integer). Must be zero.
"""
assert dimension == 0
return self
[docs]
class DefaultLine(Simplex):
"""This is the reference line with vertices (-1.0,) and (1.0,)."""
def __init__(self):
verts = ((-1.0,), (1.0,))
edges = {0: (0, 1)}
topology = {0: {0: (0,), 1: (1,)},
1: edges}
super(DefaultLine, self).__init__(LINE, verts, topology)
[docs]
def get_facet_element(self):
raise NotImplementedError()
[docs]
class UFCInterval(UFCSimplex):
"""This is the reference interval with vertices (0.0,) and (1.0,)."""
def __init__(self):
verts = ((0.0,), (1.0,))
edges = {0: (0, 1)}
topology = {0: {0: (0,), 1: (1,)},
1: edges}
super(UFCInterval, self).__init__(LINE, verts, topology)
[docs]
class DefaultTriangle(Simplex):
"""This is the reference triangle with vertices (-1.0,-1.0),
(1.0,-1.0), and (-1.0,1.0)."""
def __init__(self):
verts = ((-1.0, -1.0), (1.0, -1.0), (-1.0, 1.0))
edges = {0: (1, 2),
1: (2, 0),
2: (0, 1)}
faces = {0: (0, 1, 2)}
topology = {0: {0: (0,), 1: (1,), 2: (2,)},
1: edges, 2: faces}
super(DefaultTriangle, self).__init__(TRIANGLE, verts, topology)
[docs]
def get_facet_element(self):
return DefaultLine()
[docs]
class UFCTriangle(UFCSimplex):
"""This is the reference triangle with vertices (0.0,0.0),
(1.0,0.0), and (0.0,1.0)."""
def __init__(self):
verts = ((0.0, 0.0), (1.0, 0.0), (0.0, 1.0))
edges = {0: (1, 2), 1: (0, 2), 2: (0, 1)}
faces = {0: (0, 1, 2)}
topology = {0: {0: (0,), 1: (1,), 2: (2,)},
1: edges, 2: faces}
super(UFCTriangle, self).__init__(TRIANGLE, verts, topology)
[docs]
def compute_normal(self, i):
"UFC consistent normal"
t = self.compute_tangents(1, i)[0]
n = numpy.array((t[1], -t[0]))
return n / numpy.linalg.norm(n)
[docs]
class IntrepidTriangle(Simplex):
"""This is the Intrepid triangle with vertices (0,0),(1,0),(0,1)"""
def __init__(self):
verts = ((0.0, 0.0), (1.0, 0.0), (0.0, 1.0))
edges = {0: (0, 1),
1: (1, 2),
2: (2, 0)}
faces = {0: (0, 1, 2)}
topology = {0: {0: (0,), 1: (1,), 2: (2,)},
1: edges, 2: faces}
super(IntrepidTriangle, self).__init__(TRIANGLE, verts, topology)
[docs]
def get_facet_element(self):
# I think the UFC interval is equivalent to what the
# IntrepidInterval would be.
return UFCInterval()
[docs]
class DefaultTetrahedron(Simplex):
"""This is the reference tetrahedron with vertices (-1,-1,-1),
(1,-1,-1),(-1,1,-1), and (-1,-1,1)."""
def __init__(self):
verts = ((-1.0, -1.0, -1.0), (1.0, -1.0, -1.0),
(-1.0, 1.0, -1.0), (-1.0, -1.0, 1.0))
vs = {0: (0, ),
1: (1, ),
2: (2, ),
3: (3, )}
edges = {0: (1, 2),
1: (2, 0),
2: (0, 1),
3: (0, 3),
4: (1, 3),
5: (2, 3)}
faces = {0: (1, 3, 2),
1: (2, 3, 0),
2: (3, 1, 0),
3: (0, 1, 2)}
tets = {0: (0, 1, 2, 3)}
topology = {0: vs, 1: edges, 2: faces, 3: tets}
super(DefaultTetrahedron, self).__init__(TETRAHEDRON, verts, topology)
[docs]
def get_facet_element(self):
return DefaultTriangle()
[docs]
class IntrepidTetrahedron(Simplex):
"""This is the reference tetrahedron with vertices (0,0,0),
(1,0,0),(0,1,0), and (0,0,1) used in the Intrepid project."""
def __init__(self):
verts = ((0.0, 0.0, 0.0), (1.0, 0.0, 0.0), (0.0, 1.0, 0.0), (0.0, 0.0, 1.0))
vs = {0: (0, ),
1: (1, ),
2: (2, ),
3: (3, )}
edges = {0: (0, 1),
1: (1, 2),
2: (2, 0),
3: (0, 3),
4: (1, 3),
5: (2, 3)}
faces = {0: (0, 1, 3),
1: (1, 2, 3),
2: (0, 3, 2),
3: (0, 2, 1)}
tets = {0: (0, 1, 2, 3)}
topology = {0: vs, 1: edges, 2: faces, 3: tets}
super(IntrepidTetrahedron, self).__init__(TETRAHEDRON, verts, topology)
[docs]
def get_facet_element(self):
return IntrepidTriangle()
[docs]
class UFCTetrahedron(UFCSimplex):
"""This is the reference tetrahedron with vertices (0,0,0),
(1,0,0),(0,1,0), and (0,0,1)."""
def __init__(self):
verts = ((0.0, 0.0, 0.0), (1.0, 0.0, 0.0), (0.0, 1.0, 0.0), (0.0, 0.0, 1.0))
vs = {0: (0, ),
1: (1, ),
2: (2, ),
3: (3, )}
edges = {0: (2, 3),
1: (1, 3),
2: (1, 2),
3: (0, 3),
4: (0, 2),
5: (0, 1)}
faces = {0: (1, 2, 3),
1: (0, 2, 3),
2: (0, 1, 3),
3: (0, 1, 2)}
tets = {0: (0, 1, 2, 3)}
topology = {0: vs, 1: edges, 2: faces, 3: tets}
super(UFCTetrahedron, self).__init__(TETRAHEDRON, verts, topology)
[docs]
def compute_normal(self, i):
"UFC consistent normals."
t = self.compute_tangents(2, i)
n = numpy.cross(t[0], t[1])
return -2.0 * n / numpy.linalg.norm(n)
[docs]
class TensorProductCell(Cell):
"""A cell that is the product of FIAT cells."""
def __init__(self, *cells):
# Vertices
vertices = tuple(tuple(chain(*coords))
for coords in product(*[cell.get_vertices()
for cell in cells]))
# Topology
shape = tuple(len(c.get_vertices()) for c in cells)
topology = {}
for dim in product(*[cell.get_topology().keys()
for cell in cells]):
topology[dim] = {}
topds = [cell.get_topology()[d]
for cell, d in zip(cells, dim)]
for tuple_ei in product(*[sorted(topd)for topd in topds]):
tuple_vs = list(product(*[topd[ei]
for topd, ei in zip(topds, tuple_ei)]))
vs = tuple(numpy.ravel_multi_index(numpy.transpose(tuple_vs), shape))
topology[dim][tuple_ei] = vs
# flatten entity numbers
topology[dim] = dict(enumerate(topology[dim][key]
for key in sorted(topology[dim])))
super(TensorProductCell, self).__init__(TENSORPRODUCT, vertices, topology)
self.cells = tuple(cells)
def _key(self):
return self.cells
@staticmethod
def _split_slices(lengths):
n = len(lengths)
delimiter = [0] * (n + 1)
for i in range(n):
delimiter[i + 1] = delimiter[i] + lengths[i]
return [slice(delimiter[i], delimiter[i+1])
for i in range(n)]
[docs]
def get_dimension(self):
"""Returns the subelement dimension of the cell, a tuple of
dimensions for each cell in the product."""
return tuple(c.get_dimension() for c in self.cells)
[docs]
def construct_subelement(self, dimension):
"""Constructs the reference element of a cell subentity
specified by subelement dimension.
:arg dimension: dimension in each "direction" (tuple)
"""
return TensorProductCell(*[c.construct_subelement(d)
for c, d in zip(self.cells, dimension)])
[docs]
def volume(self):
"""Computes the volume in the appropriate dimensional measure."""
return numpy.prod([c.volume() for c in self.cells])
[docs]
def compute_reference_normal(self, facet_dim, facet_i):
"""Returns the unit normal in infinity norm to facet_i of
subelement dimension facet_dim."""
assert len(facet_dim) == len(self.get_dimension())
indicator = numpy.array(self.get_dimension()) - numpy.array(facet_dim)
(cell_i,), = numpy.nonzero(indicator)
n = []
for i, c in enumerate(self.cells):
if cell_i == i:
n.extend(c.compute_reference_normal(facet_dim[i], facet_i))
else:
n.extend([0] * c.get_spatial_dimension())
return numpy.asarray(n)
[docs]
def contains_point(self, point, epsilon=0):
"""Checks if reference cell contains given point
(with numerical tolerance)."""
lengths = [c.get_spatial_dimension() for c in self.cells]
assert len(point) == sum(lengths)
slices = TensorProductCell._split_slices(lengths)
return reduce(operator.and_,
(c.contains_point(point[s], epsilon=epsilon)
for c, s in zip(self.cells, slices)),
True)
[docs]
class UFCQuadrilateral(Cell):
"""This is the reference quadrilateral with vertices
(0.0, 0.0), (0.0, 1.0), (1.0, 0.0) and (1.0, 1.0)."""
def __init__(self):
product = TensorProductCell(UFCInterval(), UFCInterval())
pt = product.get_topology()
verts = product.get_vertices()
topology = flatten_entities(pt)
super(UFCQuadrilateral, self).__init__(QUADRILATERAL, verts, topology)
self.product = product
self.unflattening_map = compute_unflattening_map(pt)
[docs]
def get_dimension(self):
"""Returns the subelement dimension of the cell. Same as the
spatial dimension."""
return self.get_spatial_dimension()
[docs]
def construct_subelement(self, dimension):
"""Constructs the reference element of a cell subentity
specified by subelement dimension.
:arg dimension: subentity dimension (integer)
"""
if dimension == 2:
return self
elif dimension == 1:
return UFCInterval()
elif dimension == 0:
return Point()
else:
raise ValueError("Invalid dimension: %d" % (dimension,))
[docs]
def volume(self):
"""Computes the volume in the appropriate dimensional measure."""
return self.product.volume()
[docs]
def compute_reference_normal(self, facet_dim, facet_i):
"""Returns the unit normal in infinity norm to facet_i."""
assert facet_dim == 1
d, i = self.unflattening_map[(facet_dim, facet_i)]
return self.product.compute_reference_normal(d, i)
[docs]
def contains_point(self, point, epsilon=0):
"""Checks if reference cell contains given point
(with numerical tolerance)."""
return self.product.contains_point(point, epsilon=epsilon)
[docs]
class UFCHexahedron(Cell):
"""This is the reference hexahedron with vertices
(0.0, 0.0, 0.0), (0.0, 0.0, 1.0), (0.0, 1.0, 0.0), (0.0, 1.0, 1.0),
(1.0, 0.0, 0.0), (1.0, 0.0, 1.0), (1.0, 1.0, 0.0) and (1.0, 1.0, 1.0)."""
def __init__(self):
product = TensorProductCell(UFCInterval(), UFCInterval(), UFCInterval())
pt = product.get_topology()
verts = product.get_vertices()
topology = flatten_entities(pt)
super(UFCHexahedron, self).__init__(HEXAHEDRON, verts, topology)
self.product = product
self.unflattening_map = compute_unflattening_map(pt)
[docs]
def get_dimension(self):
"""Returns the subelement dimension of the cell. Same as the
spatial dimension."""
return self.get_spatial_dimension()
[docs]
def construct_subelement(self, dimension):
"""Constructs the reference element of a cell subentity
specified by subelement dimension.
:arg dimension: subentity dimension (integer)
"""
if dimension == 3:
return self
elif dimension == 2:
return UFCQuadrilateral()
elif dimension == 1:
return UFCInterval()
elif dimension == 0:
return Point()
else:
raise ValueError("Invalid dimension: %d" % (dimension,))
[docs]
def volume(self):
"""Computes the volume in the appropriate dimensional measure."""
return self.product.volume()
[docs]
def compute_reference_normal(self, facet_dim, facet_i):
"""Returns the unit normal in infinity norm to facet_i."""
assert facet_dim == 2
d, i = self.unflattening_map[(facet_dim, facet_i)]
return self.product.compute_reference_normal(d, i)
[docs]
def contains_point(self, point, epsilon=0):
"""Checks if reference cell contains given point
(with numerical tolerance)."""
return self.product.contains_point(point, epsilon=epsilon)
[docs]
def make_affine_mapping(xs, ys):
"""Constructs (A,b) such that x --> A * x + b is the affine
mapping from the simplex defined by xs to the simplex defined by ys."""
dim_x = len(xs[0])
dim_y = len(ys[0])
if len(xs) != len(ys):
raise Exception("")
# find A in R^{dim_y,dim_x}, b in R^{dim_y} such that
# A xs[i] + b = ys[i] for all i
mat = numpy.zeros((dim_x * dim_y + dim_y, dim_x * dim_y + dim_y), "d")
rhs = numpy.zeros((dim_x * dim_y + dim_y,), "d")
# loop over points
for i in range(len(xs)):
# loop over components of each A * point + b
for j in range(dim_y):
row_cur = i * dim_y + j
col_start = dim_x * j
col_finish = col_start + dim_x
mat[row_cur, col_start:col_finish] = numpy.array(xs[i])
rhs[row_cur] = ys[i][j]
# need to get terms related to b
mat[row_cur, dim_y * dim_x + j] = 1.0
sol = numpy.linalg.solve(mat, rhs)
A = numpy.reshape(sol[:dim_x * dim_y], (dim_y, dim_x))
b = sol[dim_x * dim_y:]
return A, b
[docs]
def default_simplex(spatial_dim):
"""Factory function that maps spatial dimension to an instance of
the default reference simplex of that dimension."""
if spatial_dim == 1:
return DefaultLine()
elif spatial_dim == 2:
return DefaultTriangle()
elif spatial_dim == 3:
return DefaultTetrahedron()
else:
raise RuntimeError("Can't create default simplex of dimension %s." % str(spatial_dim))
[docs]
def ufc_simplex(spatial_dim):
"""Factory function that maps spatial dimension to an instance of
the UFC reference simplex of that dimension."""
if spatial_dim == 0:
return Point()
elif spatial_dim == 1:
return UFCInterval()
elif spatial_dim == 2:
return UFCTriangle()
elif spatial_dim == 3:
return UFCTetrahedron()
else:
raise RuntimeError("Can't create UFC simplex of dimension %s." % str(spatial_dim))
[docs]
def ufc_cell(cell):
"""Handle incoming calls from FFC."""
# celltype could be a string or a cell.
if isinstance(cell, str):
celltype = cell
else:
celltype = cell.cellname()
if " * " in celltype:
# Tensor product cell
return TensorProductCell(*map(ufc_cell, celltype.split(" * ")))
elif celltype == "quadrilateral":
return UFCQuadrilateral()
elif celltype == "hexahedron":
return UFCHexahedron()
elif celltype == "vertex":
return ufc_simplex(0)
elif celltype == "interval":
return ufc_simplex(1)
elif celltype == "triangle":
return ufc_simplex(2)
elif celltype == "tetrahedron":
return ufc_simplex(3)
else:
raise RuntimeError("Don't know how to create UFC cell of type %s" % str(celltype))
[docs]
def volume(verts):
"""Constructs the volume of the simplex spanned by verts"""
# use fact that volume of UFC reference element is 1/n!
sd = len(verts) - 1
ufcel = ufc_simplex(sd)
ufcverts = ufcel.get_vertices()
A, b = make_affine_mapping(ufcverts, verts)
# can't just take determinant since, e.g. the face of
# a tet being mapped to a 2d triangle doesn't have a
# square matrix
(u, s, vt) = numpy.linalg.svd(A)
# this is the determinant of the "square part" of the matrix
# (ie the part that maps the restriction of the higher-dimensional
# stuff to UFC element
p = numpy.prod([si for si in s if (si) > 1.e-10])
return p / factorial(sd)
[docs]
def tuple_sum(tree):
"""
This function calculates the sum of elements in a tuple, it is needed to handle nested tuples in TensorProductCell.
Example: tuple_sum(((1, 0), 1)) returns 2
If input argument is not the tuple, returns input.
"""
if isinstance(tree, tuple):
return sum(map(tuple_sum, tree))
else:
return tree
[docs]
def is_hypercube(cell):
if isinstance(cell, (DefaultLine, UFCInterval, UFCQuadrilateral, UFCHexahedron)):
return True
elif isinstance(cell, TensorProductCell):
return reduce(lambda a, b: a and b, [is_hypercube(c) for c in cell.cells])
else:
return False
[docs]
def flatten_reference_cube(ref_el):
"""This function flattens a Tensor Product hypercube to the corresponding UFC hypercube"""
flattened_cube = {2: UFCQuadrilateral(), 3: UFCHexahedron()}
if numpy.sum(ref_el.get_dimension()) <= 1:
# Just return point/interval cell arguments
return ref_el
else:
# Handle cases where cell is a quad/cube constructed from a tensor product or
# an already flattened element
if is_hypercube(ref_el):
return flattened_cube[numpy.sum(ref_el.get_dimension())]
else:
raise TypeError('Invalid cell type')
[docs]
def flatten_entities(topology_dict):
"""This function flattens topology dict of TensorProductCell and entity_dofs dict of TensorProductElement"""
flattened_entities = defaultdict(list)
for dim in sorted(topology_dict.keys()):
flat_dim = tuple_sum(dim)
flattened_entities[flat_dim] += [v for k, v in sorted(topology_dict[dim].items())]
return {dim: dict(enumerate(entities))
for dim, entities in flattened_entities.items()}
[docs]
def compute_unflattening_map(topology_dict):
"""This function returns unflattening map for the given tensor product topology dict."""
counter = defaultdict(count)
unflattening_map = {}
for dim, entities in sorted(topology_dict.items()):
flat_dim = tuple_sum(dim)
for entity in entities:
flat_entity = next(counter[flat_dim])
unflattening_map[(flat_dim, flat_entity)] = (dim, entity)
return unflattening_map