# 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")