Source code for fe_utils.reference_elements
import numpy as np
[docs]class ReferenceCell(object):
def __init__(self, vertices, topology, name):
"""An object storing the geometry and topology of the reference cell.
:param vertices: a list of coordinate vectors corresponding to
the coordinates of the vertices of the cell.
:param topology: a dictionary of dictionaries such that topology[d][i]
is the list of vertices incident to the `i`-th entity of dimension
`d`.
"""
#: The vertices making up each topological entity of the reference
#: cell.
self.topology = topology
#: The list of coordinate veritices of the reference cell.
self.vertices = np.array(vertices, dtype=np.double)
self.name = name
#: The geometric and topological dimension of the reference cell.
self.dim = self.vertices.shape[1]
if self.dim != len(topology) - 1:
raise ValueError(
"Dimension mismatch between vertices and topology."
)
#: The number of entities of each dimension.
self.entity_counts = np.array([len(d) for d in topology.values()])
[docs] def point_in_entity(self, x, e):
""" Return true if the point x lies on the entity e.
:param x: The coordinate vector of the point.
:param e: The (d, i) pair describing the entity of dimension `d` and
index `i`.
"""
vertices = self.topology[e[0]][e[1]]
# Offset from first vertex.
dx = np.subtract(x, self.vertices[vertices[0]])
# Project onto space spanned by remaining vertices.
for v in vertices[1:]:
dv = np.subtract(self.vertices[v], self.vertices[vertices[0]])
dx -= dv * np.dot(dx, dv) / np.dot(dv, dv)
return (np.round(dx, 12) == 0).all()
def __repr__(self):
return self.name
#: A :class:`ReferenceCell` storing the geometry and topology of the
#: interval [0, 1].
ReferenceInterval = ReferenceCell(vertices=[[0.], [1.]],
topology={0: {0: [0],
1: [1]},
1: {0: [0, 1]}},
name="ReferenceInterval")
#: A :class:`ReferenceCell` storing the geometry and topology of the triangle
#: with vertices [[0., 0.], [1., 0.], [0., 1.]].
ReferenceTriangle = ReferenceCell(vertices=[[0., 0.], [1., 0.], [0., 1.]],
topology={0: {0: [0],
1: [1],
2: [2]},
1: {0: [1, 2],
1: [0, 2],
2: [0, 1]},
2: {0: [0, 1, 2]}},
name="ReferenceTriangle")