Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

UFCHypercube / refactor reference_element.py #68

Open
wants to merge 1 commit into
base: master
Choose a base branch
from
Open

UFCHypercube / refactor reference_element.py #68

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
231 changes: 81 additions & 150 deletions FIAT/reference_element.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -587,7 +587,7 @@ def construct_subelement(self, dimension):

def contains_point(self, point, epsilon=0):
"""Checks if reference cell contains given point
(with numerical tolerance as given by the L1 distance (aka 'manhatten',
(with numerical tolerance as given by the L1 distance (aka 'manhattan',
'taxicab' or rectilinear distance) to the cell.

Parameters
Expand All @@ -606,7 +606,7 @@ def contains_point(self, point, epsilon=0):

def distance_to_point_l1(self, point):
# noqa: D301
"""Get the L1 distance (aka 'manhatten', 'taxicab' or rectilinear
"""Get the L1 distance (aka 'manhattan', 'taxicab' or rectilinear
distance) to a point with 0.0 if the point is inside the cell.

Parameters
Expand All @@ -617,7 +617,7 @@ def distance_to_point_l1(self, point):
Returns
-------
float
The L1 distance, also known as taxicab, manhatten or rectilinear
The L1 distance, also known as taxicab, manhattan or rectilinear
distance, of the cell to the point. If 0.0 the point is inside the
cell.

Expand Down Expand Up @@ -1058,7 +1058,7 @@ def compute_reference_normal(self, facet_dim, facet_i):

def contains_point(self, point, epsilon=0):
"""Checks if reference cell contains given point
(with numerical tolerance as given by the L1 distance (aka 'manhatten',
(with numerical tolerance as given by the L1 distance (aka 'manhattan',
'taxicab' or rectilinear distance) to the cell.

Parameters
Expand All @@ -1083,7 +1083,7 @@ def contains_point(self, point, epsilon=0):
True)

def distance_to_point_l1(self, point):
"""Get the L1 distance (aka 'manhatten', 'taxicab' or rectilinear
"""Get the L1 distance (aka 'manhattan', 'taxicab' or rectilinear
distance) to a point with 0.0 if the point is inside the cell.

For more information see the docstring for the UFCSimplex method."""
Expand All @@ -1103,6 +1103,8 @@ def cell_orientation_reflection_map(self):
return make_cell_orientation_reflection_map_tensorproduct(self.cells)

def compare(self, op, other):
"""Parent-based comparison between simplicial complexes.
This is done dimension by dimension."""
if hasattr(other, "product"):
other = other.product
if isinstance(other, type(self)):
Expand All @@ -1123,59 +1125,21 @@ def __le__(self, other):
return self.compare(operator.le, other)


class UFCQuadrilateral(Cell):
r"""This is the reference quadrilateral with vertices
(0.0, 0.0), (0.0, 1.0), (1.0, 0.0) and (1.0, 1.0).

Orientation of a physical cell is computed systematically
by comparing the canonical orderings of its facets and
the facets in the FIAT reference cell.

As an example, we compute the orientation of a
quadrilateral cell:

+---3---+ +--57---+
| | | |
0 1 43 55
| | | |
+---2---+ +--42---+
FIAT canonical Mapped example physical cell
class UFCHypercube(Cell):
"""Abstract class for a reference unit hypercube [0, 1]^d with vertices in
lexicographical order."""
dimension = None
shape = None

Suppose that the facets of the physical cell
are canonically ordered as:

C = [55, 42, 43, 57]

FIAT index to Physical index map must be such that
C[0] = 55 is mapped to a vertical facet; in this
example it is:

M = [43, 55, 42, 57]

C and M are decomposed into "vertical" and "horizontal"
parts, keeping the relative orders of numbers:

C -> C0 = [55, 43], C1 = [42, 57]
M -> M0 = [43, 55], M1 = [42, 57]

Then the orientation of the cell is computed as the
following:

C0.index(M0[0]) = 1; C0.remove(M0[0])
C0.index(M0[1]) = 0; C0.remove(M0[1])
C1.index(M1[0]) = 0; C1.remove(M1[0])
C1.index(M1[1]) = 0; C1.remove(M1[1])

o = 2 * 1 + 0 = 2
"""
def __init__(self):
product = TensorProductCell(UFCInterval(), UFCInterval())
cells = [UFCInterval()] * self.dimension
product = TensorProductCell(*cells)
pt = product.get_topology()

verts = product.get_vertices()
topology = flatten_entities(pt)

super(UFCQuadrilateral, self).__init__(QUADRILATERAL, verts, topology)
super(UFCHypercube, self).__init__(self.shape, verts, topology)

self.product = product
self.unflattening_map = compute_unflattening_map(pt)
Expand All @@ -1191,14 +1155,13 @@ def construct_subelement(self, dimension):

:arg dimension: subentity dimension (integer)
"""
if dimension == 2:
sd = self.get_spatial_dimension()
if dimension > sd:
raise ValueError("Invalid dimension: %d" % (dimension,))
elif dimension == sd:
return self
elif dimension == 1:
return UFCInterval()
elif dimension == 0:
return Point()
else:
raise ValueError("Invalid dimension: %d" % (dimension,))
return ufc_hypercube(dimension)

def get_entity_transform(self, dim, entity_i):
"""Returns a mapping of point coordinates from the
Expand All @@ -1216,13 +1179,14 @@ def volume(self):

def compute_reference_normal(self, facet_dim, facet_i):
"""Returns the unit normal in infinity norm to facet_i."""
assert facet_dim == 1
sd = self.get_spatial_dimension()
assert facet_dim == sd - 1
d, i = self.unflattening_map[(facet_dim, facet_i)]
return self.product.compute_reference_normal(d, i)

def contains_point(self, point, epsilon=0):
"""Checks if reference cell contains given point
(with numerical tolerance as given by the L1 distance (aka 'manhatten',
(with numerical tolerance as given by the L1 distance (aka 'manhattan',
'taxicab' or rectilinear distance) to the cell.

Parameters
Expand All @@ -1240,14 +1204,14 @@ def contains_point(self, point, epsilon=0):
return self.product.contains_point(point, epsilon=epsilon)

def distance_to_point_l1(self, point):
"""Get the L1 distance (aka 'manhatten', 'taxicab' or rectilinear
"""Get the L1 distance (aka 'manhattan', 'taxicab' or rectilinear
distance) to a point with 0.0 if the point is inside the cell.

For more information see the docstring for the UFCSimplex method."""
return self.product.distance_to_point_l1(point)

def symmetry_group_size(self, dim):
return [1, 2, 8][dim]
return [1, 2, 8, 48][dim]

def cell_orientation_reflection_map(self):
"""Return the map indicating whether each possible cell orientation causes reflection (``1``) or not (``0``)."""
Expand All @@ -1266,109 +1230,61 @@ def __le__(self, other):
return self.product <= other


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)

def get_dimension(self):
"""Returns the subelement dimension of the cell. Same as the
spatial dimension."""
return self.get_spatial_dimension()

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,))

def get_entity_transform(self, dim, entity_i):
"""Returns a mapping of point coordinates from the
`entity_i`-th subentity of dimension `dim` to the cell.

:arg dim: entity dimension (integer)
:arg entity_i: entity number (integer)
"""
d, e = self.unflattening_map[(dim, entity_i)]
return self.product.get_entity_transform(d, e)
class UFCQuadrilateral(UFCHypercube):
r"""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 volume(self):
"""Computes the volume in the appropriate dimensional measure."""
return self.product.volume()
Orientation of a physical cell is computed systematically
by comparing the canonical orderings of its facets and
the facets in the FIAT reference cell.

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)
As an example, we compute the orientation of a
quadrilateral cell:

def contains_point(self, point, epsilon=0):
"""Checks if reference cell contains given point
(with numerical tolerance as given by the L1 distance (aka 'manhatten',
'taxicab' or rectilinear distance) to the cell.
+---3---+ +--57---+
| | | |
0 1 43 55
| | | |
+---2---+ +--42---+
FIAT canonical Mapped example physical cell

Parameters
----------
point : numpy.ndarray, list or symbolic expression
The coordinates of the point.
epsilon : float
The tolerance for the check.
Suppose that the facets of the physical cell
are canonically ordered as:

Returns
-------
bool : True if the point is inside the cell, False otherwise.
C = [55, 42, 43, 57]

"""
return self.product.contains_point(point, epsilon=epsilon)
FIAT index to Physical index map must be such that
C[0] = 55 is mapped to a vertical facet; in this
example it is:

def distance_to_point_l1(self, point):
"""Get the L1 distance (aka 'manhatten', 'taxicab' or rectilinear
distance) to a point with 0.0 if the point is inside the cell.
M = [43, 55, 42, 57]

For more information see the docstring for the UFCSimplex method."""
return self.product.distance_to_point_l1(point)
C and M are decomposed into "vertical" and "horizontal"
parts, keeping the relative orders of numbers:

def symmetry_group_size(self, dim):
return [1, 2, 8, 48][dim]
C -> C0 = [55, 43], C1 = [42, 57]
M -> M0 = [43, 55], M1 = [42, 57]

def cell_orientation_reflection_map(self):
"""Return the map indicating whether each possible cell orientation causes reflection (``1``) or not (``0``)."""
return self.product.cell_orientation_reflection_map()
Then the orientation of the cell is computed as the
following:

def __gt__(self, other):
return self.product > other
C0.index(M0[0]) = 1; C0.remove(M0[0])
C0.index(M0[1]) = 0; C0.remove(M0[1])
C1.index(M1[0]) = 0; C1.remove(M1[0])
C1.index(M1[1]) = 0; C1.remove(M1[1])

def __lt__(self, other):
return self.product < other
o = 2 * 1 + 0 = 2
"""
dimension = 2
shape = QUADRILATERAL

def __ge__(self, other):
return self.product >= other

def __le__(self, other):
return self.product <= other
class UFCHexahedron(UFCHypercube):
"""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)."""
dimension = 3
shape = HEXAHEDRON


def make_affine_mapping(xs, ys):
Expand Down Expand Up @@ -1407,6 +1323,21 @@ def make_affine_mapping(xs, ys):
return A, b


def ufc_hypercube(spatial_dim):
"""Factory function that maps spatial dimension to an instance of
the UFC reference hypercube of that dimension."""
if spatial_dim == 0:
return Point()
elif spatial_dim == 1:
return UFCInterval()
elif spatial_dim == 2:
return UFCQuadrilateral()
elif spatial_dim == 3:
return UFCHexahedron()
else:
raise RuntimeError("Can't create UFC hypercube of dimension %s." % str(spatial_dim))


def default_simplex(spatial_dim):
"""Factory function that maps spatial dimension to an instance of
the default reference simplex of that dimension."""
Expand Down
Loading