Welcome to PySOFE‘s documentation!¶

Warning

This documentation is under construction and several aspects might be missing or are not covered in detail. Especially the sections on the theoretical and implementational background are still very “work in progress”.

PySOFE User’s Guide¶

This section provides an overview on what Py$\small\mathbb{SOFE}$ does and how to use it.

To jump right in go to the Usage Example or have a look at the other Examples. For a more thorough introduction follow the Tutorial. If you want an explanation of a specific class, method, etc. please have a look at the Modules Reference.

About PySOFE¶

Py$\small\mathbb{SOFE}$ is a Python software package for solving partial differential equations using the finite element method in 1D, 2D or 3D.

Furthermore, the aim of this project is to provide background information on both the theoretical and implementational aspects of the finite element method. It is intended to serve as a starting point for people who are interested in the finite element method and especially its implementation.

Py$\small\mathbb{SOFE}$ as a project emerged from the lecture FEM: Practical Aspects given by Dr. Lars Ludwig at Dresden University of Technology during the winter term 2014/15 and is based in parts on the Matlab framework $\small\mathbb{SOFE}$ he developped for his dissertation ([Lud13]).

Although the underlying concepts are similar, Py$\small\mathbb{SOFE}$ is a separate project.

Contacts¶

If you face any problems using Py$\small\mathbb{SOFE}$, feel free to contact me:

andreas.kunze <at> mailbox.tu-dresden.de

License¶

Py$\small\mathbb{SOFE}$ is published under the BSD license.

Copyright (c) 20014-2016 Andreas Kunze

All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:

  a. Redistributions of source code must retain the above copyright notice,
     this list of conditions and the following disclaimer.
  b. Redistributions in binary form must reproduce the above copyright
     notice, this list of conditions and the following disclaimer in the
     documentation and/or other materials provided with the distribution.
  c. Neither the name of PySOFE nor the names of its contributors
     may be used to endorse or promote products derived from this software
     without specific prior written permission.


THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE FOR
ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH
DAMAGE.

Getting Started¶

This section provides an overview on how to setup and use Py$\small\mathbb{SOFE}$.

Contents

Installation
- PyPI
- Source
Usage

Installation ¶

Please note that Py$\small\mathbb{SOFE}$ has been developped and tested under Linux and Python 2.7.

A prerequisite for the installation of Py$\small\mathbb{SOFE}$ is that you already have installed Python on your computer as well as the modules setuptools and pip. If this is not the case, please follow this guide to do so.

If you have troubles installing Py$\small\mathbb{SOFE}$, please don’t hesitate to contact me.

PyPI ¶

Py$\small\mathbb{SOFE}$ is hosted at the Python Packaging Index so the easiest way of installing it would be by running:

$ pip install pysofe

in a terminal, provided that you have installed the Python module pip

Source ¶

If you would like to install Py$\small\mathbb{SOFE}$ from the source files please download the latest release from here. After the download is complete open the archive and change directory into the extracted folder. Then run the following command:

$ python setup.py install

Usage ¶

This sections is intended to provide a minimalistic worked example to show the basic usage of Py$\small\mathbb{SOFE}$.

Consider the linear Poisson equation on the unit square in 2D with homogeneous Dirichlet boundary conditions

\[\begin{align*} a\Delta u(x) &= f(x) &&\text{ in }\Omega=(0,1)^2 \\ u(x) &= 0 &&\text{ on }\partial\Omega \end{align*}\]

with the constant coefficient $a$ and the right hand site $f \in L^2(\Omega)$. To keep things simple let $a = 1$ and define $f(x) = x_0^2 - x_1^2$ where $x = (x_0,x_1)\in\Omega\subset\mathbb{R}^2$.

First, we create the mesh that discretizes the spatial domain $\Omega$ of our problem. To do so we import the predefined class UnitSquareMesh and instantiate a mesh with $100$ nodes on each axis

>>> from pysofe import UnitSquareMesh
>>> mesh = UnitSquareMesh(100, 100)

Then, we create the reference element that provides the basis functions. We will use linear basis functions on triangles, implemented in the class P1 which takes the spatial dimension of our problem as an argument

>>> from pysofe import P1
>>> element = P1(dimension=2)

Next, we create the function space in which we look for a solution. We do this by creating an instance of the FESpace class that brings together the mesh and reference element

>>> from pysofe import FESpace
>>> fe_space = FESpace(mesh, element)

For the formulation of our problem we also need to specify boundary conditions the solution should comply with. This part takes a bit more effort than the previous steps.

The first thing to do for this is to define a function specifying the boundary where the conditions should hold. Py$\small\mathbb{SOFE}$ expects a function that decides for a given array of points which of them lie on the boundary. In our case the boundary points are determined by their $x_0$ and $x_1$ coordinate which are equal to $0$ or $1$ on the boundary of the unit square

>>> from numpy import logical_or as or_
>>> def dirichlet_domain(x):
...     x0_is_0_or_1 = or_(x[0] == 0., x[0] == 1.)
...     x1_is_0_or_1 = or_(x[1] == 0., x[1] == 1.)
...     return or_(x0_is_0_or_1, x1_is_0_or_1)

Then, since we want to impose a Dirichlet boundary condition, we need to define the function that the solution of our problem should be equal to on the boundary. In the simple case of a homogeneous boundary condition it suffices to define this function as a constant

>>> g = 0.

Now, we can create the boundary condition implemented in the DirichletBC class and pass the arguments defined above

>>> from pysofe import DirichletBC
>>> dirichlet_bc = DirichletBC(fe_space, dirichlet_domain, g)

What remains is to define the actual boundary value problem. Therefore, we will use the predefined Poisson class. But first we need to define the parameters for this class which are the factor $a$ and the right hand site function $f$. Since we set the constant factor $a$ equal to $1$ we can simply do so in the code as well

>>> a = 1.

and we define the right hand site as

>>> def f(x):
...     return x[0]*x[0] - x[1]*x[1]

So, now we have all things we need to create the object that represent our boundary value problem

>>> from pysofe import Poisson
>>> pde = Poisson(fe_space, a, f, dirichlet_bc)

Finally, to solve it we call

>>> u = pde.solve()

which returns a callable function object that represents our approximate solution and can be visualized with

>>> import pysofe
>>> pysofe.show(u)

which should produce the following graphics.

Tutorial¶

This section intends to explain how to use the different parts of the Py$\small\mathbb{SOFE}$ library more thoroughly than the Usage part of the Getting Started section.

to be continued...

Examples¶

This section provides some specific examples of problems involving partial differential equations that will be solved using Py$\small\mathbb{SOFE}$.

to be continued...

Change Log¶

Unreleased¶

Added¶

Neumann and periodic boundary conditions
Mesh generator using signed distance functions

Release 0.1.0¶

Added¶

Lagrange elements P1, P2
Basic mesh data structures and uniform refinement method
Gauss-Legendre quadrature for simplicial domains
FE space class
Mass matrix, L2 product, Laplacian and L2 projection operators
Dirichlet boundary conditions
Poisson equation class
Basic documentation and tests

PySOFE Background Notes¶

This section is intended to provide information on both the theoretical background of Py$\small\mathbb{SOFE}$ and the implementational aspects.

Implementational Aspects¶

This section is intended to show how the theoretical concepts are actually implemented in the software.

to be continued...

The Mesh¶

In the finite element environment the purpose of the mesh is to discretize and approximate the spatial domain $\Omega$ of the considered partial differential equation.

In Py$\small\mathbb{SOFE}$ the mesh is defined by its geometry and a topology implemented in the respective classes MeshGeometry and MeshTopology which is a concept similar to the one presented in [Log09]. Furthermore, it stores a family of reference maps, implemented in the ReferenceMap class, that connect the physical mesh entities to a reference domain.

The Mesh Geometry¶

The MeshGeometry class provides geometrical information about the mesh which currently amounts in storing the spatial coordinate of the mesh nodes.

The Mesh Topology¶

The MeshTopology class provides topological information about the mesh, i.e. it gives access to its entities such as edges, faces or cells and neighborly relation, e.g. the neighboring cells to each edge.

All the information is stored using incidence matrices. So, for every pair $(d,d')$ of topological dimensions $0 \leq d,d' \leq D$ where $D \in \{ 1,2,3 \}$ is the spatial dimension of the mesh, there exists a matrix $M_{d,d'}$ that for the $d$-dimensional entities stores their incident entities of topological dimension $d'$. The entries $m_{ij}$ are either $1$ or $0$ depending on whether the $j$-th $d'$-dimensional mesh entity is incident to the $i$-th mesh entity of topological dimension $d$.

To explain how to get to these matrices we will consider the following simple sample mesh in 2D.

$\coordinate [label={below left:{\small $1$}}] (1) at (0,0); \coordinate [label={below right:{\small $2$}}] (2) at (2,0); \coordinate [label={above left:{\small $3$}}] (3) at (0,2); \coordinate [label={above right:{\small $4$}}] (4) at (2,2); \draw (3) -- (1) -- (2) -- (3) -- (4) -- (2); \draw (0.6,0.6) node[circle, draw, inner sep=1pt, minimum size=0pt] {\small 1}; \draw (1.4,1.4) node[circle, draw, inner sep=1pt, minimum size=0pt] {\small 2}; \draw (1.0,-0.2) node {\tiny\underline{1}}; \draw (-0.1,1.0) node {\tiny\underline{2}}; \draw (1.2,1.0) node {\tiny\underline{3}}; \draw (2.1,1.0) node {\tiny\underline{4}}; \draw (1.0,2.2) node {\tiny\underline{5}};$

which consists of 4 vertices, 5 edges and 2 triangular cells. For this mesh the incidence relation $D \to 0$, i.e. the vertices incident to each cell $(D=2)$, would be stored in a matrix

\[M_{2,0} = \begin{pmatrix} 1 & 1 & 1 & 0 \\ 0 & 1 & 1 & 1 \end{pmatrix}\]

This matrix shows that e.g. the first cell, which corresponds to the first row in the matrix, is incident to the vertices with indices $1,2,3$, as those are the non-zero column indices which represent the mesh vertices.

Note that this mesh is very simple and coarse. If the meshes get finer the number of non-zero entries will become very small compared to the total number of entries, i.e. the incidence matrices will be sparse. To make use of this fact we will utilize the sparse matrix capabilities of the scipy package which provides sparse matrix classes that only store the non-zero entries, hence saving lots of memory.

Similarly to [Log09] the construction of these matrices relies on the three algorithms build, transpose and intersection, which will be explained in the following sections.

Contents

Initialize
Build
Transpose
Intersection
- The case $d = d'$
- The case $d > d'$
Computing Incidence Relations

Initialize ¶

A prerequisite for the algorithms to work is, that we already have the incidence matrix for the relation $D \to 0$. This matrix is constructed using the following auxiliary intialization method.

First, we need the connectivity array that defines the mesh cells (entities of topological dimension $D$) row-wise via their node indices. For the sample mesh this array would be

\[\begin{bmatrix} 1 & 2 & 3 \\ 2 & 3 & 4 \end{bmatrix},\]

i.e. the first cell is define by the node indices $1,2,3$ and the second cell is defined by the nodes $2,3,4$. This array has to be supplied by the user and passed as an argument when creating an instance of the MeshTopology class.

In the end we will represent the incidence relation using a sparse matrix in coordinate format (COO) which stores its non-zero entries by their ‘cordinates’ $(i,j)$, i.e. their row and column indices. To do so, three 1-dimensional arrays are used. The first one is the data array which contains the values of the non-zero entries. In our case this will be an array of $1$‘s. The second array will contain the row indices and the third one the respective column indices.

To create the row index array we need to know the number of vertices that define the cells of the mesh. For a mesh of triangles this number is equal to $nV = 3$ (in 3D for tetrahedra it would be $nV = 4$). Remember that we want to construct the incidence relation between the mesh cells and the mesh vertices, so if each cell is defined by $nV$ vertices there will be that many non-zero entries in every row of the matrix. Hence, for each row index $i$ (which corresponds to one cell in the mesh) there will be $nV$ coordinate tuples $(i,j)$ that have this index as a first component. Therefore, the number of vertices tells us how often each cell index appears in the row array.

For our sample mesh the row index array would be:

\[rows = \begin{bmatrix} 1 & 1 & 1 & 2 & 2 & 2 \end{bmatrix},\]

containing each row index $nV = 3$ times.

The columns of the incidence matrix we create correspond to the vertices of the mesh. So, each index in the column array is the column index of the non-zero entry corresponding to the respective index in the row array. As those non-zero entries represent the incident vertices for each cell the column indices are given by the connectivity array that has been passed as an argument.

Again, for our sample mesh this column index array would be:

\[cols = \begin{bmatrix} 1 & 2 & 3 & 2 & 3 & 4 \end{bmatrix}.\]

As stated before, the remaining data array contains the values of the non-zero entries in the matrix. For us those values are all equal to $1$, so the data array would be

\[data = \begin{bmatrix} 1 & 1 & 1 & 1 & 1 & 1 \end{bmatrix}.\]

and we can now create the incidence matrix for the relation $D \to 0$ by passing the three array to the scipy class for a sparse matrix in COO format.

Build ¶

The build algorithm computes the incidence relation $d \to 0$ for a topological dimension $0 < d < D$, i.e. it determines the incident vertices for every $d$-dimensional mesh entity.

Starting from the incidence relation $D \to 0$ that was constructed using the initialization method above, we first have to determine the local vertex sets for every $d$-dimensional subentity of the mesh cells, i.e. we need for each subentity of topological dimension $d$ of every mesh cell to vertex indices that define it.

For example, assume we want to build the incidence relation $1 \to 0$, i.e. the incident vertices of every edge (topological dimension $d = 1$). As every triangular cell is incident to three edges and each edge is defined by two vertex indices, the local vertex set for the edges of every cell will consist of three 2-tuples of indices.

For our sample mesh the local vertex set for the edges of the first cell would be

\[\begin{bmatrix} (1,2) & (1,3) & (2,3) \end{bmatrix}\]

because the first edge of the cell is defined by the vertex indices $1$ and $2$, the second edge is defined by the indices $1,3$ and the third edge by the indices $2,3$. Analogously, the local vertex for the edges of the second cell would be

\[\begin{bmatrix} (2,3) & (2,4) & (3,4) \end{bmatrix}.\]

Obviously, there are duplicates in the local vertex sets because edges will be incident to two neighboring triangular cells if they are not part of the boundary of the domain. So, the next step is to determine the unique vertex tuples in all the local vertex sets. For us this would result in the following array, where the tuples are now written as the rows.

\[\begin{bmatrix} 1 & 2 \\ 1 & 3 \\ 2 & 3 \\ 2 & 4 \\ 3 & 4 \end{bmatrix}.\]

This array is now a connectivity array that defines the edges of our sample mesh just like the one we had in the initialization method for the cells. That means we can use the same procedure to construct the sparse matrix for the incidence relation $d \to 0$.

As there are five edges in the mesh and each edge is incident to two vertices the row index array for the edges will contain each index from $1$ to $5$ twice.

\[rows = \begin{bmatrix} 1 & 1 & 2 & 2 & 3 & 3 & 4 & 4 & 5 & 5 \end{bmatrix}\]

The column index array just consists of the indices given by the connectivity array flattened row-wise

\[cols = \begin{bmatrix} 1 & 2 & 1 & 3 & 2 & 3 & 2 & 4 & 3 & 4 \end{bmatrix}\]

and the data array is an array of $1$‘s with the same length as the row and column index arrays.

\[data = \begin{bmatrix} 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 \end{bmatrix}\]

Transpose ¶

The transpose algorithm computes the incidence relation $d \to d'$ for $d < d'$ from $d' \to d$.

This is the simplest and most straight forward of the algorithms. Assuming we already have the incidence matrix $M_{d',d}$ for the relation $d' \to d$ we get the incidence matrix $M_{d,d'}$ simply by transposing

\[M_{d,d'} = M_{d',d}^{T}\]

For example, we already know the incident vertices to every mesh edge

\[M_{1,0} = \begin{pmatrix} 1 & 1 & 0 & 0 \\ 1 & 0 & 1 & 0 \\ 0 & 1 & 1 & 0 \\ 0 & 1 & 0 & 1 \\ 0 & 0 & 1 & 1 \end{pmatrix}.\]

If we now want to know the incident edges for every mesh vertex we just have to transpose this matrix to get

\[M_{1,0}^{T} = M_{0,1} = \begin{pmatrix} 1 & 1 & 0 & 0 & 0 \\ 1 & 0 & 1 & 1 & 0 \\ 0 & 1 & 1 & 0 & 1 \\ 0 & 0 & 0 & 1 & 1 \\ \end{pmatrix}.\]

where we can see for each of the four mesh vertices the edges that are incident to them.

Intersection ¶

The intersection algorithm computes the incidence relation $d \to d'$ for $d \geq d'$ from the two relations $d \to d''$ and $d'' \to d'$.

The auxiliary topological dimension $d''$ is set to $0$, so we assume that we already have the incidence matrices $M_{d,0}, M_{0,d'}$ for the incidence relations $d \to 0$ and $0 \to d'$.

Then, the intersection is computed as the dot product of these two matrices. The resulting intersection matrix $I_{d,d'}$ describes for every $d$-dimensional mesh entity (represented by the rows) how many vertices ($d'' = 0$-dimensional entities) it shares with each of the mesh entities of topological dimension $d$ (represented by the columns).

Next, we have to distinguish between two cases to define when the incidence relation $d \to d'$ holds for two entities.

The case $d = d'$¶

In case $d = d'$, e.g. if we want the neighboring (incident) cells for every mesh cell ($2 \to 2$), we define two entities of topological dimension $d$ incident to each other if they share a $d - 1 = \tilde{d}$-dimensional mesh entity. So, for two triangles ($d = 2$) to be incident they must share an edge ($\tilde{d} = 1$).

Simplicial entities of dimension $d$, i.e. triangles for $d = 2$ or tetrahedra for $d = 3$, which are the only type of entities that currently are supported by Py$\small\mathbb{SOFE}$, are defined by $d + 1$ vertices. So, as the computed intersection matrix $I_{d,d'}$ tells us the number of vertices shared between each of the $d$-dimensional entities and every entity of topological dimension $d'$ the incidence matrix $M_{d,d'}$ for $d = d'$ is defined by

\[(M_{d,d})_{ij} = m_{ij} = \begin{cases} 1 & (I_{d,d})_{ij} = d \\ 0 & \text{else} \end{cases}\]

Suppose, we want to know the neighboring edges for every edge of our sample mesh, i.e. we are interested in the incidence relation $1 \to 1$. The necessary incidence matrices for computing the intersection are

\[M_{1,0} = \begin{pmatrix} 1 & 1 & 0 & 0 \\ 1 & 0 & 1 & 0 \\ 0 & 1 & 1 & 0 \\ 0 & 1 & 0 & 1 \\ 0 & 0 & 1 & 1 \end{pmatrix} \quad M_{0,1} = \begin{pmatrix} 1 & 1 & 0 & 0 & 0 \\ 1 & 0 & 1 & 1 & 0 \\ 0 & 1 & 1 & 0 & 1 \\ 0 & 0 & 0 & 1 & 1 \\ \end{pmatrix}.\]

Then the intersection matrix is

\[M_{1,0} \cdot M_{0,1} = I_{1,1} = \begin{pmatrix} 2 & 1 & 1 & 1 & 0 \\ 1 & 2 & 1 & 0 & 1 \\ 1 & 1 & 2 & 1 & 1 \\ 1 & 0 & 1 & 2 & 1 \\ 0 & 1 & 1 & 1 & 2 \end{pmatrix}\]

and therefore the resulting incidence matrix would be

\[M_{1,1} = \begin{pmatrix} 0 & 1 & 1 & 1 & 0 \\ 1 & 0 & 1 & 0 & 1 \\ 1 & 1 & 0 & 1 & 1 \\ 1 & 0 & 1 & 0 & 1 \\ 0 & 1 & 1 & 1 & 0 \end{pmatrix}.\]

A special case is the relation $0 \to 0$ for which we define that vertices are incident to themselves only. So the incidence matrix $M_{0,0}$ would be the identity matrix.

The case $d > d'$¶

In case $d > d'$ we define a entity of topological dimension $d'$ to be incident to a $d$-dimensional one if all its vertices are incident to this entity of dimension $d$.

Since the $d'$-dimensional mesh entities are defined by $d' + 1$ vertices the incidence matrix $M_{d,d'}$ is defined similarly to the previous section via the intersection matrix

\[(M_{d,d'})_{ij} = m_{ij} = \begin{cases} 1 & (I_{d,d'})_{ij} = d' + 1 \\ 0 & \text{else} \end{cases}\]

Suppose, we want to know the incident edges for every cell of our sample mesh, i.e. we are interested in the incidence relation $2 \to 1$. The necessary incidence matrices for computing the intersection are

\[M_{2,0} = \begin{pmatrix} 1 & 1 & 1 & 0 \\ 0 & 1 & 1 & 1 \end{pmatrix} \quad M_{0,1} = \begin{pmatrix} 1 & 1 & 0 & 0 & 0 \\ 1 & 0 & 1 & 1 & 0 \\ 0 & 1 & 1 & 0 & 1 \\ 0 & 0 & 0 & 1 & 1 \\ \end{pmatrix}.\]

Then the intersection matrix is

\[M_{2,0} \cdot M_{0,1} = I_{2,1} = \begin{pmatrix} 2 & 2 & 2 & 1 & 1 \\ 1 & 1 & 2 & 2 & 2 \end{pmatrix}\]

and therefore the resulting incidence matrix would be

\[M_{2,1} = \begin{pmatrix} 1 & 1 & 1 & 0 & 0 \\ 0 & 0 & 1 & 1 & 1 \end{pmatrix}.\]

Computing Incidence Relations ¶

Now that we have the algorithms build, transpose and intersection we can construct the incidence matrix for any relation $d \to d'$ by a combination of those algorithms as follows

Compute relation (d -> d')
--------------------------

if relation (d -> 0) does not exist
   build relation (d -> 0)

if relation (d' -> 0) does not exist
   build relation (d' -> 0)

if d < d'
   compute relation (d' -> d)
   transpose relation (d' -> d)
else
   intersect relation (d -> dd)

Mesh Generation¶

One of the first steps in the practical application of the finite element method is creating a mesh to discretize the domain of the considered differential equation. Every subsequent step of the method strongly depends on the quality of the generated mesh. It is therefore indispensable to use a mesh of good quality to get an accurate solution.

A mesh consists of node positions and connectivity information. While it is possible for simple geometries like a square to explicitly define or generate the nodes and cells this approach becomes infeasible once the geometries get more complicated.

In Py$\small\mathbb{SOFE}$ we implement a method proposed by Per-Olof Persson and Gilbert Strang in their paper [PS04]. They make use of the analogy between a simplicial mesh and a truss structure and by assuming a force-displacement function for the bars of the structure, i.e. the edges of the mesh, iteratively generate meshes by solving for a force equilibrium.

The Algorithm¶

to be continued...

The Reference Maps¶

Within the finite element method there are certain operations that have to be done on every mesh element. Examples include the evaluation of basis functions and their derivatives or the evaluation of a function in quadrature points. Doing this for each mesh entity $K$ separately becomes infeasable when the number of elements increases.

Hence, those operations are rather pulled back on a single reference element $\hat{K}$ using a socalled reference map $\Phi_{K}:\hat{K}\to K$ for each mesh element $K \in \mathcal{T}_h$ where $\mathcal{T}_h$ is our triangulation, i.e. the union of all the elements that discretize the spatial domain of the considered partial differential equation.

\draw (4,-1) -- (5,-1) -- (6,-1) -- (6,0) -- (6,1) -- (5,1) -- (4,1) -- (4,0) -- (4,-1); \draw (4,0) -- (5,0) -- (6,0); \draw (5,-1) -- (5,0) -- (5,1); \draw (5,-1) -- (4,0); \draw (6,-1) -- (5,0) -- (4,1); \draw (6,0) -- (5,1); \draw (0,-0.5) -- (1,-0.5) -- (0,0.5) -- (0,-0.5); \node at (0.2,-0.2) {\small $\hat{x}$}; \node at (0.5,-0.8) {\small $\hat{K}$}; \draw [->, line width=1pt] (0.3,-0.1) to [out=50, in=130] (4.6,-0.2); \node at (2.5,1.0) {\small $\Phi_{K_{7}}$}; \node at (4.7,-0.4) {\small $x$}; \node at (5.7,-1.5) {\small $K_{7}$}; \draw (4.8,-0.6) -- (5.5, -1.3);

Now, every point $x \in K$ in some mesh element $K \in \mathcal{T}_h$ can be expressed as

\[x = \Phi_{K}(\hat{x})\]

for an appropriate point $\hat{x}\in\hat{K}$ on the reference element. An example where this will be used is numerical integration. Instead of storing or computing the quadrature points for the domain of every single mesh element we just have to store the quadrature points for the reference domain and can map them to each of the physical mesh elements when needed.

Another advantage of this approach is, that we don’t have to explicitly define all the basis functions corresponding to the degrees of freedom of the finite dimensional subspace $V_h$ (INSERT LINK TO THEORY). Consider herefore the above mesh with $9$ nodes and suppose we use linear Lagrange basis functions. In this case the degrees of freedom correspond to the nodes of the mesh and the basis functions $\mathcal{B} = \{\varphi_j, j=1,\ldots,9\}$ are defined to be equal to $1$ at exactly one node $x_k, k=1,\ldots,9$ and vanish at all other nodes, i.e.

\[\varphi_j(x_k) = \begin{cases} 1 & j = k \\ 0 & j \not= k \end{cases}\]

Since we pull back the evaluation of basis functions to the reference element we only have to define $3$ reference basis functions $\hat{\varphi}_i, i=1,2,3$ corresponding to the nodes of the reference element. Then, by means of the reference maps we can evaluate any basis function $\varphi \in \mathcal{B}$ in any point $x \in K$ by evaluating the appropriate reference basis function $\hat{\varphi}$ in the corresponding local point $\hat{x} \in \hat{K}$.

\[\varphi(x)\big|_{x\in K} = \varphi(\Phi_K(\hat{x})) = \hat{\varphi}(\hat{x})\]

Evaluation of the Reference Maps¶

The reference maps are used in many different ways. They have to provide zero order information, i.e. ordinary evaluation, to map local points on the reference domain to their global counterparts on the physical mesh elements. First order information is needed e.g. to compute Jacobians arising when we transform integrals to the reference domain or for the evaluation of derivatives of the global basis functions.

Furthermore it doesn’t suffice to be able to map points from the reference element domain to the mesh elements domains. We also need to be able to map local points to the mesh subentities, e.g. edges in 2D.

To show how the reference maps are described and evaluated we consider the 2-dimensional case and straight-sided triangular elements. Then the reference map $\Phi_K$ for $K \in \mathcal{T}_h$ can be defined as

\[\Phi_{K}(\mathbf{\hat{x}}) = \mathbf{c}_1 + (\mathbf{c}_2 - \mathbf{c}_1\ \mathbf{c}_3 - \mathbf{c}_1)\cdot\mathbf{\hat{x}}\]

where $\mathbf{c}_1, \mathbf{c}_2, \mathbf{c}_3 \in\mathbb{R}^2$ are the coordinates of the nodes of the mesh element $K$. An alternative way of describing this in terms of the local basis functions $\{ \hat{\varphi}_1, \hat{\varphi}_2, \hat{\varphi}_3 \}$ is

\[\Phi_{K}(\mathbf{\hat{x}}) = \mathbf{c}_1 \hat{\varphi}_1(\mathbf{\hat{x}}) + \mathbf{c}_2 \hat{\varphi}_2(\mathbf{\hat{x}}) + \mathbf{c}_3 \hat{\varphi}_3(\mathbf{\hat{x}}) = \sum_{i=1}^{3} \mathbf{c}_i \hat{\varphi}_i(\mathbf{\hat{x}})\]

Recall herefore that $\hat{\varphi}_i(\mathbf{\hat{c}}_j) = \delta_{ij}$, where $\{\mathbf{\hat{c}}_j, j=1,2,3\}$ are the coordinates of the reference element’s nodes. This means e.g. for $\mathbf{\mathbf{\hat{x}}} = \mathbf{\hat{c}}_1$ that

\[\hat{\varphi}_1(\mathbf{\hat{c}}_1) = 1 \qquad \hat{\varphi}_2(\mathbf{\hat{c}}_1) = 0 \qquad \hat{\varphi}_3(\mathbf{\hat{c}}_1) = 0\]

so the local node $\mathbf{\hat{c}}_1$ is correctly mapped to its global counterpart

\[\begin{align*} \Phi_{K}(\mathbf{\hat{c}}_1) &= \mathbf{c}_1 \hat{\varphi}_1(\mathbf{\hat{c}}_1) + \mathbf{c}_2 \hat{\varphi}_2(\mathbf{\hat{c}}_1) + \mathbf{c}_3 \hat{\varphi}_3(\mathbf{\hat{c}}_1) \\ &= \mathbf{c}_1 \cdot 1 + \mathbf{c}_2 \cdot 0 + \mathbf{c}_3 \cdot 0 \\ &= \mathbf{c}_1. \end{align*}\]

So the reference map $\Phi_K$ for each mesh element $K$ is evaluated by

\[\Phi_{K}(\mathbf{\hat{x}}) = \sum_{i=1}^{nB} \mathbf{c}_i \hat{\varphi}_1(\mathbf{\hat{x}})\]

where $\mathbf{c}_i, i=1,\ldots,nB$ are the coordinates of the mesh elements nodes and $nB$ is the number of local basis functions defined for the reference element ($2$ for edges, $3$ for cells in case of triangles).

In the code all reference maps are evaluated simultaneously for a given point or even for a given array of local points on the reference domain as follows.

Let $\mathbf{\hat{X}} \in \mathbb{R}^{nP \times nD}$ be the array of local points, so $nP$ is the number of points and $nD$ is their dimension (in our case $nD = 2$ for local points on the reference triangle).

First we evaluate the local basis functions in those points which will give us an array $B \in\mathbb{R}^{nB \times nP}$ with values between $0$ and $1$

\[B = \left[\hat{\varphi}_1(\mathbf{\hat{X}})\ \ldots\ \hat{\varphi}_{nB}(\mathbf{\hat{X}})\right]^{T} \quad\text{with } \hat{\varphi}_i(\mathbf{\hat{X}}) \in \mathbb{R}^{nP}, i=1,\ldots,nB.\]

Next, we need the coordinates of the nodes of every mesh element which will result in an array $C \in \mathbb{R}^{nE \times nB \times nD}$ where $nE$ is the number of mesh elements, so

\[C = [C_1\ \ldots\ C_{nE}] \quad\text{with } C_k \in\mathbb{R}^{nB \times nD}, k=1,\ldots,nE.\]

Then the reference map for each element is evaluated as the matrix product

\[\Phi_k = B^{T}C_k \in \mathbb{R}^{nP \times nD} \qquad k = 1,\ldots,nE.\]

which will give us the array of the family of reference maps evaluated in each given point $\Phi \in \mathbb{R}^{nE \times nP \times nD}$.

The Reference Element¶

In the sense of Ciarlet [Cia78] a finite element is defined as a triple $(\hat{K}, P(\hat{K}), \Sigma(\hat{K}))$ where $\hat{K} \subset \mathbb{R}^d$ is a closed bounded subset with non-empty interior, $P(\hat{K})$ is a finite dimensional vector space on $\hat{K}$ and $\Sigma(\hat{K})$ is a basis of the dual space $P'(\hat{K})$, i.e. the space of linear functionals on $P(\hat{K})$.

In most cases $\hat{K}$ is a polygonial domain and $P(\hat{K})$ is a polynomial function space. The elements $\hat{\varphi}_j, j=1,\ldots,n$ of a basis of $P(\hat{K})$, with $n\in\mathbb{N}$ being its dimension, are called shape functions and the linear functionals $N_i \in\Sigma(\hat{K}), i=1,\ldots,n$ define the socalled degrees of freedom.

The shape functions $\hat{\varphi}_j$ are determined by the linear functionals $N_i$ via

\[N_{i}(\hat{\varphi}_{j}) = \delta_{ij} = \begin{cases} 1 & i = j \\ 0 & \text{otherwise} \end{cases}\]

so the choice of these functionals provides different families of finite elements (see [LB13]).

The purpose of the reference element is to provide methods for the evaluation of these shape or basis functions $\hat{\varphi}$ defined on the reference domain and their derivatives.

Lagrange Elements¶

One of the most widely used family of finite elements are the Lagrange elements, also often called Courant elements, which were first defined in [Cou43] with use of Lagrange interpolation polynomials. Their defining functionals $N_{i}$ are given by

\[N_i(v) = v(\xi_i),\quad i=1,\ldots,n\]

where $\xi_i \in\mathbb{R}^d$ are specific node points. As each basis function is therefore associated with a particular node they form a socalled nodal basis of $P(\hat{K})$.

It follows from the equations (1) and (2) that

\[\hat{\varphi}_{i}(\xi_{j}) = \delta_{ij}\]

for the nodes $\xi_j$ on the elements’ domain. The number of nodes required for the definition of the shape or basis functions is determined by their polynomial order.

$\mathbb{P}_1$ Elements¶

The Lagrange elements $\mathbb{P}_1$ define linear shape functions on simplicial domains. The nodes $\xi_i$ each basis function is associated with are the vertices of the element.

$\mathbb{P}_2$ Elements¶

The Lagrange elements $\mathbb{P}_2$ define quadratic shape functions and in addition to the elements’ vertices use the midpoints of the elements’ edges as node for their definition.

References¶

This section contains references to the work and ideas of other authors that had an impact on Py$\small\mathbb{SOFE}$.

[Cia78]

Philippe G. Ciarlet. The Finite Element Method for Elliptic Problems. Studies in Mathematics and its Applications. Elsevier Science, 1978.

[Cou43]

Richard Courant. Variational methods for the solution of problems of equilibrium and vibrations. Bulletin of the American Mathematical Society, 49:1–23, 1943.

[LB13]

Mats G. Larson and Frederik Bengzon. The Finite Element Method: Theory, Implementation and Applications. Texts in Computational Science and Engineering. Springer Berlin Heidelberg, 2013.

[Log09]

Anders Logg. Efficient representation of computational meshes. International Journal of Computational Science and Engineering, 4(4):283–295, 2009.

[Lud13]

Lars Ludwig. Analytical investigations and numerical experiments for singularly perturbed convection-diffusion problems with layers and singularities using a newly developed FE-software. PhD thesis, Dresden University of Technology, 2013. URL: http://nbn-resolving.de/urn:nbn:de:bsz:14-qucosa-137301.

[PS04]

Per-Olof Persson and Gilbert Strang. A simple mesh generator in matlab. SIAM review, 46(2):329–345, 2004.

PySOFE Modules Reference¶

This section provides a detailed view of the data structures available in Py$\small\mathbb{SOFE}$.

`elements` Module¶

This module provides the data structures that represent finite elements. The purpose of these data structures is to evaluate the basis functions and their derivatives at points on the reference domain.

`elements.base` Module¶

Provides the base class for all finite element classes.

class pysofe.elements.base.Element(dimension, order, n_basis, n_verts)[source]¶

Provides an abstract base class for all reference elements.

Derived classes have to implement the methods

_eval_d0_basis(points)()
_eval_d1_basis(points)()
_eval_d2_basis(points)()

that evaluate the basis functions or their derivatives.

Parameters:	dimension (int) – The spatial dimension of the element order (int) – The polynomial order of the basis functions n_basis (iterable) – The number of basis functions related to the entities of every topological dimension n_verts (iterable) – The number of vertices for the entities of each topological dimension

dimension¶: The spatial dimension of the element.

order¶: The polynomial order of the element’s basis functions.

n_basis¶: The number of basis functions associated with the entities of each spatial (sub-) dimenson.

n_verts¶: The number of vertices of the element’s (sub-) entities of each spatial dimension.

dof_tuple¶: The number of degrees of freedom associated to the entities of every spatial dimension.

eval_basis(points, deriv=0)[source]¶

Evaluates the element’s basis functions or their derivatives at given local points on the reference domain.

Parameters:	points (array_like) – The local points at which to evaluate the basis functions deriv (int) – The derivation order

`elements.simple` Module¶

Provides classes for simple finite elements.

`elements.simple.lagrange` Module¶

Provides classes for some simple Lagrange type finite elements.

class pysofe.elements.simple.lagrange.P1(dimension)[source]¶

Linear Lagrange basis functions on simplicial domains.

Parameters:	dimension (int) – The spatial dimension of the element

class pysofe.elements.simple.lagrange.P2(dimension)[source]¶

Quadratic Lagrange basis functions on simplicial domains.

Parameters:	dimension (int) – The spatial dimension of the element

`meshes` Module¶

Provides the data structures for the finite element meshes.

The main purpose of the mesh in the finite element environment is to discretize the spatial domain of the considered partial differential equation. Furthermore, it gives access to its topological entities such as vertices, edges, faces and cells and provides methods for refinement.

The mesh data structure in Py$\small\mathbb{SOFE}$ is represented by the Mesh class and encapsulates the classes MeshGeometry and MeshTopology. It also has an instance of the ReferenceMap class as an attribute to connect the physical mesh entities with the reference domain.

`meshes.mesh` Module¶

Provides the data structure for approximating the spatial domain of the partial differential equation.

class pysofe.meshes.mesh.Mesh(nodes, connectivity)[source]¶

Provides a class for general meshes.

Basically clues together the MeshGeometry and MeshTopology classes and provides interfaces for mesh refinement and function evaluation.

Parameters:	nodes (array_like) – The coordinates of the mesh nodes connectivity (array_like) – The connectivity array defining the mesh cells via their vertex indices

dimension¶: The spatial dimension fo the mesh.

nodes¶: The coordinates of the mesh nodes.

cells¶: The incident vertex indices of the mesh cells.

facets¶: The incident vertex indices of the mesh facets.

faces¶: The incident vertex indices of the mesh faces.

edges¶: The incident vertex indices of the mesh edges.

boundary(fnc=None)[source]¶

Determines the mesh facets that are part of the boundary specified by the function fnc and returns a corresponding boolean array.

Parameters:	fnc (callable) – Function specifying some part of the boundary for which to return the corresponding facets

refine(method='uniform', **kwargs)[source]¶

Refines the mesh using the given method.

Parameters:	method (str) – A string specifying the refinement method to use

eval_function(fnc, points, local=True)[source]¶

Evaluates a given function in the global mesh points corresponding to the given local points on the reference domain.

Parameters:	fnc (callable) – The function to evaluate points (array_like) – The local points on the reference domain
Returns:	nE x nP [x ...]
Return type:	numpy.ndarray

search(points, return_preimage=True)[source]¶

Determines for every given global query point its containing mesh cell.

If return_preimage is True (default) it also returns the corresponding local points on the reference domain.

Parameters:	points (array_like) – The global query points return_preimage (bool) – Whether to return the corresponding local points on the reference domain

class pysofe.meshes.mesh.UnitSquareMesh(nx=10, ny=10)[source]¶

Discretizes the domain of the unit square in 2D.

Parameters:	nx (int) – Number of nodes in the x-direction ny (int) – Number of nodes in the y-direction

`meshes.geometry` Module¶

Provides the data structure that supplies geometrical information on the meshes.

class pysofe.meshes.geometry.MeshGeometry(nodes)[source]¶

Stores geometrical information of a mesh.

Parameters:	nodes (array_like) – The coordinates of the mesh grid points

nodes¶: The coordinates of the mesh grid points

`meshes.topology` Module¶

Provides the data structure that supplies topological information on meshes.

class pysofe.meshes.topology.MeshTopology(cells, dimension)[source]¶

Stores topological information on a mesh such as its entities and connectivity information.

Parameters:	cells (array_like) – The connectivity array of the mesh cells dimension (int) – The spatial dimension of the mesh cells

n_entities¶: The numbers of elements for every topological dimension

get_connectivity(d, dd, return_indices=False)[source]¶

Returns the incidence relation d -> dd, i.e. for each each d-dimensional mesh entity its incident entities of topological dimension dd.

Parameters:	dd (d,) – The topological dimensions of the entities involved return_indices (bool) – Whether to return the indices of the incident entities instead of the sparse incidence matrix

get_entities(d)[source]¶

Returns the array of all d-dimensional mesh entities represented by the indices of their incident vertices.

Parameters:	d (int) – The topological dimension of the mesh entities to return

get_boundary(d)[source]¶

Returns boolean array specifying the boundary entities of topological dimension d.

Parameters:	d (int) – The topological dimension of the desired boundary entities

`meshes.reference_map` Module¶

Provides the data structure for the family of reference maps.

class pysofe.meshes.reference_map.ReferenceMap(mesh)[source]¶

Establishes a connection between the reference domain and the physical mesh elements.

Parameters:	mesh (pysofe_light.meshes.mesh.Mesh) – The mesh instance

eval(points, deriv=0, mask=None)[source]¶

Evaluates each member of the family of reference maps or their derivatives at given local points.

Provides the following forms of information:

Zero order information, i.e. the ordinary evaluation, is needed to compute global points given their local counterparts
First order information is used to compute Jacobians, e.g. in integral transformations

Parameters:	points (array_like) – The local points at which to evaluate deriv (int) – The derivation order mask (array_like) – Integer or boolean mask specifying certain elements of which to evaluate the reference maps
Returns:	nE x nP x nD [x nD]
Return type:	numpy.ndarray

eval_inverse(points, hosts)[source]¶

Evaluates the inverse maps of the reference maps corresponding to the given host elements at given global points.

Parameters:	points (array_like) – The global points for which to evaluate the inverse mappings hosts (array_like) – The host element for each of the global points

jacobian_inverse(points, mask=None)[source]¶

Returns the inverse of the reference maps’ jacobians evaluated at given points.

Parameters:	points (array_like) – The local points at which to evaluate the jacobians mask (array_like) – Integer or boolean mask specifying certain elements of which to evaluate the inverse jacobians

jacobian_determinant(points, mask=None)[source]¶

Returns the determinants of the reference maps’ jacobians evaluated at given points.

Parameters:	points (array_like) – The local points at which to compute the determinants mask (array_like) – Integer or boolean mask specifying certain elements of which to evaluate the jacobian determinants

`quadrature` Module¶

Provides the data structure for numerical integration.

The purpose of these data structures is to give access to the quadrature points and weights for several spatial domains.

`quadrature.base` Module¶

Provides the base class for all quadrature rules.

class pysofe.quadrature.base.QuadRule(order, dimension)[source]¶

Provides an abstract base class for all quadrature rules.

Parameters:	order (int) – The polynomial order up to which the quadrature should be exact

`quadrature.gaussian` Module¶

Provides the data structures for Gauss-Legendre quadrature on simplicial domains.

class pysofe.quadrature.gaussian.GaussQuadSimp(order, dimension)[source]¶

Provides Gauss-Legendre quadrature points and weights for standard simplicial domains in several spatial dimensions.

Parameters:	order (int) – The polynomial order up to which the numerical integration should be exact dimension (int) – The spatial dimension up to which to provide quadrature data

`spaces` Module¶

Provides the data structures for finite element spaces as well as for functionals and operators on them.

`spaces.space` Module¶

Provides the data structures that represents a finite element space.

class pysofe.spaces.space.FESpace(mesh, element)[source]¶

Base class for all finite element spaces.

Connects the mesh and reference element via a degrees of freedom manager.

Parameters:	mesh (pysofe.meshes.mesh.Mesh) – The mesh used for approximating the pde domain element (pysofe.elements.base.Element) – The reference element

get_quadrature_data(d)[source]¶

Returns the quadrature points and weighths associated with the d-dimensional entities as well as the jacobian determinants of for integral transformation to the reference domain.

Parameters:	d (int) – The topological dimension of the entities for which to return the quadrature data

eval_global_derivatives(points, deriv=1)[source]¶

Evaluates the global basis functions’ derivatives at given local points.

Parameters:	points (array_like) – The local points on the reference element deriv (int) – The derivation order
Returns:	(nE x nB x nP x nD [x nD]) array containing for all elements (nE) the evaluation of all basis functions first derivatives (nB) in each point (nP)
Return type:	numpy.ndarray

eval_dofs(dofs, points, deriv=0, local=True)[source]¶

Evaluates the linear combination of the basis functions or their derivatives in the given points w.r.t. the given values for the degrees of freedom.

Parameters:	dofs (numpy.ndarray) – The degrees of freedom values points (array_like) – The points in which to evaluate the function deriv (int) – The derivation order local (bool) – Whether the points are local on the reference domain

`spaces.manager` Module¶

Provides the data structure for the degrees of freedom manager.

class pysofe.spaces.manager.DOFManager(mesh, element)[source]¶

Connects the finite element mesh and reference element through the assignment of degrees of freedom to individual elements.

Parameters:	mesh (pysofe.meshes.Mesh) – The mesh for which to provide connectivity information element (pysofe.elements.base.Element) – The reference element

n_dof¶: The total number of degrees of freedom

get_dof_map(d, mask=None)[source]¶

Returns the degrees of freedom mapping that connects the global mesh entities of topological dimension d to the local reference element.

Parameters:	d (int) – The topological dimension of the entities for which to return the degrees of freedom mapping mask (array_like) – An 1d array marking certain entities of which to get the dof map

extract_dofs(d, mask=None, return_indices=False)[source]¶

Returns an index or boolean (default) array specifying the degrees of freedom associated with the mesh entities of topological dimension d.

Parameters:	d (int) – The topological dimension of the mesh entities mask (array_like) – An 1d array marking certain entities of which to get the dofs return_indices (bool) – Whether to return an index array or not

`spaces.operators` Module¶

Provides data structures that represent weak forms differential operators.

class pysofe.spaces.operators.Operator(fe_space)[source]¶

Base class for all operators.

Derived classes have to implement the method _compute_entries()

Parameters:	fe_space (pysofe_light.spaces.space.FESpace) – A function space the operator works on

assemble(codim=0, mask=None)[source]¶

Assembles the discrete weak operator.

Parameters:	codim (int) – The codimension of the entities for which to assemble mask (array_like) – Boolean 1d array marking specific entities for assembling
Returns:
Return type:	scipy.sparse.lil_matrix

class pysofe.spaces.operators.MassMatrix(fe_space, c=1)[source]¶

Represents the operator

\[\int_{\Omega} c uv\]

where $u,v \in V$ and $c \in L^{2}$.

Parameters:	fe_space (pysofe_light.spaces.space.FESpace) – The function space the operator works on c (scalar, callable) – The function factor

class pysofe.spaces.operators.L2Product(fe_space, f=1)[source]¶

Represents the operator

\[\int_{\Omega} fv\]

where $v \in V$ and $f \in L^{2}$.

Parameters:	fe_space (pysofe_light.spaces.space.FESpace) – The function space the operator works on f (scalar, callable) – The function factor

class pysofe.spaces.operators.Laplacian(fe_space, a=1)[source]¶

Represents the operator

\[\int_{\Omega} a \nabla u \cdot \nabla v\]

where $u,v \in V$ and $a \in L^{2}$.

Parameters:	fe_space (pysofe_light.spaces.space.FESpace) – The function space the operator works on a (scalar, callable) – The function factor

class pysofe.spaces.operators.RDiffusion(fe_space, q=1.0)[source]¶

Provides the discrete operator for

\[\int_{\Omega} q \cdot \nabla v\]

Parameters:	fe_space (pysofe.spaces.FESpace) – The finite element space the operator lives on

class pysofe.spaces.operators.L2Projection[source]¶

Provides an object for the $L^{2}$-projection of a given function.

static project(fnc, fe_space, codim=0, mask=None)[source]¶

Projects the given function to the given finite element space.

Parameters:	fnc (scalar, callable) – The function to project fe_space – The function space onto which to project codim (int) – The codimension of the considered entities mask (array_like) – Boolean array marking entities onto which to project

`spaces.functions` Module¶

Provides convinience classes for functions in the fem framework.

class pysofe.spaces.functions.FEFunction(fe_space, dof_values)[source]¶

A finite element function defined via degrees of freedoms.

Parameters:	fe_space (pysofe.spaces.space.FESpace) – The function space dof_values (array_like) – Values for the degrees of freedom of the function space

order¶: The polynomial order of the function.

class pysofe.spaces.functions.PeriodicFEFunction(fe_space, dof_values, period)[source]¶

A periodic finite element function.

Parameters:	fe_space (pysofe.spaces.space.FESpace) – The function space dof_values (array_like) – Values for the degrees of freedom of the function space period (scalar, iterable) – Uniform period or period for each axis direction

class pysofe.spaces.functions.MeshFunction(fnc, mesh)[source]¶

Wrapper for the evaluation of a given function on a specific mesh.

Parameters:	fnc (callable) – The function mesh (pysofe.meshes.mesh.Mesh) – The considered mesh

`pde` Module¶

Provides the data structures that represent specific partial differential equations.

`pde.base` Module¶

Provides the base class for all pde objects.

class pysofe.pde.base.PDE(fe_space, bc)[source]¶

Base class for all partial differential equations.

Parameters:	fe_space (pysofe.spaces.space.FESpace) – The function space of the pde bc (pysofe.pde.conditions.BoundaryCondition, iterable) – The boundary condition(s)

assemble()[source]¶: Assembles all discrete operators that are part of the pde.

apply_conditions()[source]¶: Applies the boundary conditions to the stiffness matrix and load vector.

solve()[source]¶: Solves the partial differential equation.

`pde.poisson` Module¶

Provides the data structure that represents the Poisson equation.

class pysofe.pde.poisson.Poisson(fe_space, a=1.0, f=0.0, bc=None)[source]¶

Represents the linear Poisson equation

\[- a \Delta u = f\]

Parameters:	fe_space (pysofe.spaces.space.FESpace) – The considered function space a (callable) – The coefficient function f (callable) – The right hand site function bc (pysofe.pde.conditions.BoundaryCondition, iterable) – The boundary condition(s)

`pde.conditions` Module¶

Provides the data structures that represent the boundary conditions for the pdes.

class pysofe.pde.conditions.BoundaryCondition(fe_space, domain)[source]¶

Base class for all boundary conditions.

Parameters:	fe_space (pysofe.spaces.space.FESpace) – The considered function space domain (callable) – A function specifying the boundary part where the condition should hold

apply(A=None, b=None)[source]¶: Applies the boundary condition to the stiffness matrix A and/or load vector b.

class pysofe.pde.conditions.DirichletBC(fe_space, domain, g=0)[source]¶

This class represents Dirichlet boundary conditions of the form

\[u(x,t) = u_{D}|_{\Gamma_{D}}, \Gamma_{D}\subseteq\partial\Omega\]

Parameters:	fe_space (pysofe.spaces.space.FESpace) – The considered function space domain (callable) – A function specifying the boundary part where the condition should hold g (callable) – A function specifying the values at the boundary

class pysofe.pde.conditions.NeumannBC(fe_space, domain, h=0)[source]¶

Represents Neumann boundary conditions of the form

\[\frac{\partial u}{\partial n} = h|_{\Gamma_{N}\subseteq\partial\Omega}\]

where $n$ denotes the outer unit normal on $\Gamma_{N}$.

Parameters:	fe_space (pysofe.spaces.space.FESpace) – The considered function space domain (callable) – A function specifying the boundary part where the condition should hold h (callable) – A function specifying the values at the boundary

class pysofe.pde.conditions.PeriodicBC(fe_space, master_domain, slave_domain)[source]¶

Represents periodic boundary conditions for n-orthotopes.

Parameters:	fe_space (pysofe.spaces.FESpace) – The considered function space master_domain (callable) – Functions specifying the master facets slave_domain (callable) – Functions specifying the slave facets

`utils` Module¶

Provides some auxilliary functions that are used in various points in the software.

pysofe.utils.unique_rows(A, return_index=False, return_inverse=False)[source]¶

Returns B, I, J where B is the array of unique rows from the input array A and I, J are arrays satisfying A = B[J,:] and B = A[I,:].

Parameters:	A (numpy.ndarray) – The 2d array of which to determine the unique rows return_index (bool) – Whether to return I return_inverse (bool) – Whether to return J

pysofe.utils.lagrange_nodes(dimension, order)[source]¶

Returns the nodes that determine the Lagrange shape functions of given order on a simplicial domain.

Parameters:	dimension (int) – The spatial dimension of the points order (int) – The polynomial order of the shape functions

pysofe.utils.match_nodes(nodes0, nodes1, dim=0)[source]¶: Determines index sets I, J such that nodes0[I, d] == nodes1[J, d] for each d != dim.

pysofe.utils.int2bool(arr, size=None)[source]¶

Turns an integer index array into a boolean mask.

Parameters:	arr (array_like) – The 1d array with the indices of the later True values size (int) – The size of the boolean mask, if None the maximum integer index will be used

pysofe.utils.bool2int(mask)[source]¶

Returns the indices of the nonzero entries in the given 1D boolean mask.

Parameters:	mask (array_like) – The boolean mask

pysofe.utils.mesh_area(mesh)[source]¶

Computes the area of a 2D triangular mesh using Heron’s formula.

Given the lengths $a,b,c$ of the edges of a triangle its area is given by

\[\sqrt{s(s-a)(s-b)(s.c)}\]

where $s = \frac{a + b + c}{2}$.

Welcome to PySOFE‘s documentation!¶

PySOFE User’s Guide¶

About PySOFE¶

Contacts¶

License¶

Getting Started¶

Tutorial¶

Examples¶

Change Log¶

Unreleased¶

Added¶

Release 0.1.0¶

Added¶

PySOFE Background Notes¶

Implementational Aspects¶

The Mesh¶

The Mesh Geometry¶

The Mesh Topology¶

Mesh Generation¶

The Algorithm¶

The Reference Maps¶

Evaluation of the Reference Maps¶

The Reference Element¶

Lagrange Elements¶

\(\mathbb{P}_1\) Elements¶

\(\mathbb{P}_2\) Elements¶

References¶

PySOFE Modules Reference¶

elements Module¶

elements.base Module¶

elements.simple Module¶

elements.simple.lagrange Module¶

meshes Module¶

meshes.mesh Module¶

meshes.geometry Module¶

meshes.topology Module¶

meshes.reference_map Module¶

quadrature Module¶

quadrature.base Module¶

quadrature.gaussian Module¶

spaces Module¶

spaces.space Module¶

spaces.manager Module¶

spaces.operators Module¶

spaces.functions Module¶

pde Module¶

pde.base Module¶

pde.poisson Module¶

pde.conditions Module¶

utils Module¶

Indices and tables¶

`elements` Module¶

`elements.base` Module¶

`elements.simple` Module¶

`elements.simple.lagrange` Module¶

`meshes` Module¶

`meshes.mesh` Module¶

`meshes.geometry` Module¶

`meshes.topology` Module¶

`meshes.reference_map` Module¶

`quadrature` Module¶

`quadrature.base` Module¶

`quadrature.gaussian` Module¶

`spaces` Module¶

`spaces.space` Module¶

`spaces.manager` Module¶

`spaces.operators` Module¶

`spaces.functions` Module¶

`pde` Module¶

`pde.base` Module¶

`pde.poisson` Module¶

`pde.conditions` Module¶

`utils` Module¶