Source code for fe_utils.function_spaces

import numpy as np
from . import ReferenceTriangle, ReferenceInterval
from .finite_elements import LagrangeElement, lagrange_points
from matplotlib import pyplot as plt
from matplotlib.tri import Triangulation

[docs]class FunctionSpace(object): def __init__(self, mesh, element): """A finite element space. :param mesh: The :class:`~.mesh.Mesh` on which this space is built. :param element: The :class:`~.finite_elements.FiniteElement` of this space. Most of the implementation of this class is left as an :ref:`exercise <ex-function-space>`. """ #: The :class:`~.mesh.Mesh` on which this space is built. self.mesh = mesh #: The :class:`~.finite_elements.FiniteElement` of this space. self.element = element raise NotImplementedError # Implement global numbering in order to produce the global # cell node list for this space. #: The global cell node list. This is a two-dimensional array in #: which each row lists the global nodes incident to the corresponding #: cell. The implementation of this member is left as an #: :ref:`exercise <ex-function-space>` self.cell_nodes = None #: The total number of nodes in the function space. self.node_count =, mesh.entity_counts) def __repr__(self): return "%s(%s, %s)" % (self.__class__.__name__, self.mesh, self.element)
[docs]class Function(object): def __init__(self, function_space, name=None): """A function in a finite element space. The main role of this object is to store the basis function coefficients associated with the nodes of the underlying function space. :param function_space: The :class:`FunctionSpace` in which this :class:`Function` lives. :param name: An optional label for this :class:`Function` which will be used in output and is useful for debugging. """ #: The :class:`FunctionSpace` in which this :class:`Function` lives. self.function_space = function_space #: The (optional) name of this :class:`Function` = name #: The basis function coefficient values for this :class:`Function` self.values = np.zeros(function_space.node_count)
[docs] def interpolate(self, fn): """Interpolate a given Python function onto this finite element :class:`Function`. :param fn: A function ``fn(X)`` which takes a coordinate vector and returns a scalar value. """ fs = self.function_space # Create a map from the vertices to the element nodes on the # reference cell. cg1 = LagrangeElement(fs.element.cell, 1) coord_map = cg1.tabulate(fs.element.nodes) cg1fs = FunctionSpace(fs.mesh, cg1) for c in range(fs.mesh.entity_counts[-1]): # Interpolate the coordinates to the cell nodes. vertex_coords = fs.mesh.vertex_coords[cg1fs.cell_nodes[c, :], :] node_coords =, vertex_coords) self.values[fs.cell_nodes[c, :]] = [fn(x) for x in node_coords]
[docs] def plot(self, subdivisions=None): """Plot the value of this :class:`Function`. This is quite a low performance plotting routine so it will perform poorly on larger meshes, but it has the advantage of supporting higher order function spaces than many widely available libraries. :param subdivisions: The number of points in each direction to use in representing each element. The default is :math:`2d+1` where :math:`d` is the degree of the :class:`FunctionSpace`. Higher values produce prettier plots which render more slowly! """ fs = self.function_space d = subdivisions or ( 2 * ( + 1) if > 1 else 2 ) if fs.element.cell is ReferenceInterval: fig = plt.figure() fig.add_subplot(111) # Interpolation rule for element values. local_coords = lagrange_points(fs.element.cell, d) elif fs.element.cell is ReferenceTriangle: fig = plt.figure() ax = fig.gca(projection='3d') local_coords, triangles = self._lagrange_triangles(d) else: raise ValueError("Unknown reference cell: %s" % fs.element.cell) function_map = fs.element.tabulate(local_coords) # Interpolation rule for coordinates. cg1 = LagrangeElement(fs.element.cell, 1) coord_map = cg1.tabulate(local_coords) cg1fs = FunctionSpace(fs.mesh, cg1) for c in range(fs.mesh.entity_counts[-1]): vertex_coords = fs.mesh.vertex_coords[cg1fs.cell_nodes[c, :], :] x =, vertex_coords) local_function_coefs = self.values[fs.cell_nodes[c, :]] v =, local_function_coefs) if fs.element.cell is ReferenceInterval: plt.plot(x[:, 0], v, 'k') else: ax.plot_trisurf(Triangulation(x[:, 0], x[:, 1], triangles), v, linewidth=0)
@staticmethod def _lagrange_triangles(degree): # Triangles linking the Lagrange points. return (np.array([[i / degree, j / degree] for j in range(degree + 1) for i in range(degree + 1 - j)]), np.array( # Up triangles [np.add(np.sum(range(degree + 2 - j, degree + 2)), (i, i + 1, i + degree + 1 - j)) for j in range(degree) for i in range(degree - j)] # Down triangles. + [ np.add( np.sum(range(degree + 2 - j, degree + 2)), (i+1, i + degree + 1 - j + 1, i + degree + 1 - j)) for j in range(degree - 1) for i in range(degree - 1 - j) ]))
[docs] def integrate(self): """Integrate this :class:`Function` over the domain. :result: The integral (a scalar).""" raise NotImplementedError