# shenfun.forms package¶

## shenfun.forms.arguments module¶

class shenfun.forms.arguments.Array(space, val=0, buffer=None)[source]

Bases: shenfun.forms.arguments.ShenfunBaseArray

Numpy array for TensorProductSpace

The Array is the result of a Function evaluated on its quadrature mesh.

The Function is the product of all 1D basis expansions, that for each dimension is defined like

$u(x) = \sum_{k \in \mathcal{K}} \hat{u}_k \psi_k(x),$

where $$\psi_k(x)$$ are the trial functions and $$\{\hat{u}_k\}_{k\in\mathcal{K}}$$ are the expansion coefficients. Here an index set $$\mathcal{K}=0, 1, \ldots, N$$ is used to simplify notation.

For an M+1-dimensional TensorProductSpace with coordinates $$x_0, x_1, \ldots, x_M$$ we get

$u(x_{0}, x_{1}, \ldots, x_{M}) = \sum_{k_0 \in \mathcal{K}_0}\sum_{k_1 \in \mathcal{K}_1} \ldots \sum_{k_M \in \mathcal{K}_M} \hat{u}_{k_0, k_1, \ldots k_M} \psi_{k_0}(x_0) \psi_{k_1}(x_1) \ldots \psi_{k_M}(x_M),$

where $$\mathcal{K}_j$$ is the index set for the wavenumber mesh along axis $$j$$.

Note that for a Cartesian mesh in 3D it would be natural to use coordinates $$(x, y, z) = (x_0, x_1, x_2)$$ and the expansion would be the simpler and somewhat more intuitive

$u(x, y, z) = \sum_{l \in \mathcal{K}_0}\sum_{m \in \mathcal{K}_1} \sum_{n \in \mathcal{K}_2} \hat{u}_{l, m, n} \psi_{l}(x) \psi_{m}(y) \psi_{n}(z).$

The Array’s values (the Numpy array) represent the left hand side, evaluated on the Cartesian quadrature mesh. With this we mean the $$u(x_i, y_j, z_k)$$ array, where $$\{x_i\}_{i=0}^{N_0}$$, $$\{y_j\}_{j=0}^{N_1}$$ and $$\{z_k\}_{k=0}^{N_2}$$ represent the mesh along the three directions. The quadrature mesh is then

$(x_i, y_j, z_k) \quad \forall \, (i, j, k) \in [0, 1, \ldots, N_0] \times [0, 1, \ldots, N_1] \times [0, 1, \ldots, N_2]$

The entire spectral Galerkin function can be obtained using the Function class.

Parameters
• val (int or float) – Value used to initialize array

• buffer (Numpy array, Function or sympy Expr) – If array it must be of correct shape. A sympy expression is evaluated on the quadrature mesh and the result is used as buffer.

• .. note:: For more information, see numpy.ndarray <https (//docs.scipy.org/doc/numpy/reference/generated/numpy.ndarray.html>_)

Examples

>>> from mpi4py import MPI
>>> from shenfun import Basis, TensorProductSpace, Function
>>> K0 = FunctionSpace(8, 'F', dtype='D')
>>> K1 = FunctionSpace(8, 'F', dtype='d')
>>> FFT = TensorProductSpace(MPI.COMM_WORLD, [K0, K1])
>>> u = Array(FFT)

forward(output_array=None)[source]

Return Function used to evaluate Array

property forward_output

Return whether self is the result of a forward transform

offset()[source]

Return offset of this basis

The offset is the number of scalar TensorProductSpacees ahead of this Arrays space in a :class:.MixedTensorProductSpace.

shenfun.forms.arguments.Basis(*args, **kwargs)[source]
class shenfun.forms.arguments.BasisFunction(space, index=0, basespace=None, offset=0, base=None)[source]

Bases: object

Base class for arguments to shenfun’s Exprs

Parameters
property argument

Return argument of basis

property base

Return base

property basespace

Return base space

property dimensions

Return dimensions of function space

expr_rank()[source]

Return rank of expression involving basis

function_space()[source]

Return function space of BasisFunction

index()[source]

Return index into base space

num_components()[source]

Return number of components in basis

offset()[source]

Return offset of this basis

The offset is the number of scalar TensorProductSpacees ahead of this space in a :class:.MixedTensorProductSpace.

property rank

Return rank of basis

class shenfun.forms.arguments.Expr(basis, terms=None, scales=None, indices=None)[source]

Bases: object

Class for spectral Galerkin forms

An Expression that is linear in TestFunction, TrialFunction or Function. Exprs are used as input to inner() or project().

Parameters
• terms (list of list of lists of length dimension) – Describes the differential operations performed on the basis function in the Expr

The triply nested terms lists are such that

• the outermost list represents a tensor component. There is one item for each tensor component. If the Expr is a scalar with rank = 0, then len(terms) = 1. For vectors it equals the number of dimensions and for second order tensors it equals ndim**2

• the second nested list represents the different terms in the form, that may be more than one. For example, the scalar valued div(grad(u)) has three terms in 3D:

$\partial^2u/\partial x^2 + \partial^2u/\partial y^2 + \partial^2u/\partial z^2$
• the last inner list represents the differential operations for each term and each tensor component, stored for each as a list of length = dim

The Expr div(grad(u)), where u is a scalar, is as such represented as a nested list of shapes (1, 3, 3), 1 meaning it’s a scalar, the first 3 because the Expr consists of the sum of three terms, and the last 3 because it is 3D. The entire representation is:

[[[2, 0, 0],
[0, 2, 0],
[0, 0, 2]]]


where the first [2, 0, 0] term has two derivatives in first direction and none in the others, the second [0, 2, 0] has two derivatives in second direction, and the last [0, 0, 2] has two derivatives in the last direction and none in the others.

• scales (list of lists) – Representing a scalar multiply of each inner product. Note that the scalar can also be a function of coordinates (using sympy). There is one scale for each term in each tensor component, and as such it is a list of lists.

• indices (list of lists) – Index into MixedTensorProductSpace. Only used when basis of form has rank > 0. There is one index for each term in each tensor component, and as such it is a list of lists.

Examples

>>> from shenfun import *
>>> from mpi4py import MPI
>>> comm = MPI.COMM_WORLD
>>> C0 = FunctionSpace(16, 'F', dtype='D')
>>> C1 = FunctionSpace(16, 'F', dtype='D')
>>> R0 = FunctionSpace(16, 'F', dtype='d')
>>> T = TensorProductSpace(comm, (C0, C1, R0))
>>> v = TestFunction(T)
>>> e = div(grad(v))
>>> e.terms()
[[[2, 0, 0], [0, 2, 0], [0, 0, 2]]]
>>> e2 = grad(v)
>>> e2.terms()
[[[1, 0, 0]], [[0, 1, 0]], [[0, 0, 1]]]


Note that e2 in the example has shape (3, 1, 3). The first 3 because it is a vector, the 1 because each vector item contains one term, and the final 3 since it is a 3-dimensional tensor product space.

property argument

Return argument of Expr’s basis

property base

Return base BasisFunction used in Expr

basis()[source]

Return basis of Expr

basis_rank()[source]

Return rank of Expr’s BasisFunction

property dimensions

Return ndim of Expr

eval(x, output_array=None)[source]

Return expression evaluated on x

Parameters

x (float or array of floats) – Array must be of shape (D, N), for N points in D dimensions

expr_rank()[source]

Return rank of Expr

function_space()[source]

Return function space of basis in Expr

index()[source]
indices()[source]

Return indices of Expr

num_components()[source]

Return number of components in Expr

num_terms()[source]

Return number of terms in Expr

property rank

Return rank of Expr’s BasisFunction

scales()[source]

Return scales of Expr

simplify()[source]

Join terms, that are otherwise equal, using scale

terms()[source]

Return terms of Expr

tolatex(symbol_names=None, funcname='u', replace=None)[source]
tosympy(basis=None, psi=x, y, z, r, s)[source]

Return self evaluated with a sympy basis

Parameters
• basis (sympy Expr or string) – if sympy Expr, then use this directly if str, then use that as name for a generic sympy Function

• psi (tuple of sympy Symbols)

Examples

>>> from shenfun import FunctionSpace, TensorProductSpace
>>> import sympy as sp
>>> theta, r = psi = sp.symbols('x,y', real=True, positive=True)
>>> rv = (r*sp.cos(theta), r*sp.sin(theta))
>>> F0 = FunctionSpace(8, 'F', dtype='d')
>>> T0 = FunctionSpace(8, 'C')
>>> T = TensorProductSpace(comm, (F0, T0), coordinates=(psi, rv))
>>> u = TrialFunction(T)
>>> ue = r*sp.cos(theta)
>>> du = div(grad(u))
>>> du.tosympy()
Derivative(u(x, y), (y, 2)) + Derivative(u(x, y), y)/y + Derivative(u(x, y), (x, 2))/y**2
>>> du.tosympy(basis=ue, psi=psi)
0

class shenfun.forms.arguments.Function(space, val=0, buffer=None)[source]

Bases: shenfun.forms.arguments.ShenfunBaseArray, shenfun.forms.arguments.BasisFunction

Spectral Galerkin function for given TensorProductSpace or Basis()

The Function is the product of all 1D basis expansions, that for each dimension is defined like

$u(x) = \sum_{k \in \mathcal{K}} \hat{u}_k \psi_k(x),$

where $$\psi_k(x)$$ are the trial functions and $$\{\hat{u}_k\}_{k\in\mathcal{K}}$$ are the expansion coefficients. Here an index set $$\mathcal{K}=0, 1, \ldots, N$$ is used to simplify notation.

For an M+1-dimensional TensorProductSpace with coordinates $$x_0, x_1, \ldots, x_M$$ we get

$u(x_{0}, x_{1}, \ldots, x_{M}) = \sum_{k_0 \in \mathcal{K}_0}\sum_{k_1 \in \mathcal{K}_1} \ldots \sum_{k_M \in \mathcal{K}_M} \hat{u}_{k_0, k_1, \ldots k_M} \psi_{k_0}(x_0) \psi_{k_1}(x_1) \ldots \psi_{k_M}(x_M),$

where $$\mathcal{K}_j$$ is the index set for the wavenumber mesh along axis $$j$$.

Note that for a Cartesian mesh in 3D it would be natural to use coordinates $$(x, y, z) = (x_0, x_1, x_2)$$ and the expansion would be the simpler and somewhat more intuitive

$u(x, y, z) = \sum_{l \in \mathcal{K}_0}\sum_{m \in \mathcal{K}_1} \sum_{n \in \mathcal{K}_2} \hat{u}_{l, m, n} \psi_{l}(x) \psi_{m}(y) \psi_{n}(z).$

The Function’s values (the Numpy array) represent the $$\hat{u}$$ array. The trial functions $$\psi$$ may differ in the different directions. They are chosen when creating the TensorProductSpace.

Parameters
• space (TensorProductSpace)

• val (int or float) – Value used to initialize array

• buffer (Numpy array, Function or sympy Expr) – If array it must be of correct shape. A sympy expression is evaluated on the quadrature mesh and forward transformed to create the buffer array.

Note

Examples

>>> from mpi4py import MPI
>>> from shenfun import FunctionSpace, TensorProductSpace, Function
>>> K0 = FunctionSpace(8, 'F', dtype='D')
>>> K1 = FunctionSpace(8, 'F', dtype='d')
>>> T = TensorProductSpace(MPI.COMM_WORLD, [K0, K1])
>>> u = Function(T)
>>> K2 = FunctionSpace(8, 'C', bc=(0, 0))
>>> T2 = TensorProductSpace(MPI.COMM_WORLD, [K0, K1, K2])
>>> v = Function(T2)

assign(u_hat)[source]

Assign self to u_hat of possibly different size

Parameters

u_hat (Function) – Function of possibly different shape than self. Must have the same function_space

backward(output_array=None, uniform=False)[source]

Return Function evaluated on quadrature mesh

eval(x, output_array=None)[source]

Evaluate Function at points

Parameters
• points (float or array of floats)

• coefficients (array) – Expansion coefficients

• output_array (array, optional) – Return array, function values at points

Examples

>>> import sympy as sp
>>> K0 = FunctionSpace(9, 'F', dtype='D')
>>> K1 = FunctionSpace(8, 'F', dtype='d')
>>> T = TensorProductSpace(MPI.COMM_WORLD, [K0, K1], axes=(0, 1))
>>> X = T.local_mesh()
>>> x, y = sp.symbols("x,y")
>>> ue = sp.sin(2*x) + sp.cos(3*y)
>>> ul = sp.lambdify((x, y), ue, 'numpy')
>>> ua = Array(T, buffer=ul(*X))
>>> points = np.random.random((2, 4))
>>> u = ua.forward()
>>> u0 = u.eval(points).real
>>> assert np.allclose(u0, ul(*points))

property forward_output

Return whether self is the result of a forward transform

mask_nyquist(mask=None)[source]

Set self to have zeros in Nyquist coefficients

refine(N, output_array=None)[source]

Return self with new number of quadrature points

Parameters

N (number or sequence of numbers) – The new number of quadrature points

Note

If N is smaller than for self, then a truncated array is returned. If N is greater than before, then the returned array is padded with zeros.

set_boundary_dofs()[source]
to_ortho(output_array=None)[source]

Project Function to orthogonal basis

shenfun.forms.arguments.FunctionSpace(N, family='Fourier', bc=None, dtype='d', quad=None, domain=None, scaled=None, padding_factor=1.0, dealias_direct=False, coordinates=None, **kw)[source]

Return function space

Parameters
• N (int) – Number of quadrature points

• family (str, optional) – Choose one of

• Chebyshev or C,

• Legendre or L,

• Fourier or F,

• Laguerre or La,

• Hermite or H

• Jacobi or J

• bc (str or tuple, optional) – Choose one of

• 2-tuple (a, b) - Dirichlet boundary condition with $$v(-1)=a$$ and $$v(1)=b$$.

• 4-tuple (a, b, 0, 0) - Biharmonic with the two non-zero Dirichlet conditions $$v(-1)=a$$ and $$v(1)=b$$.

• Dirichlet - Homogeneous Dirichlet.

• Neumann - Homogeneous Neumann.

• Biharmonic - Homogeneous Dirichlet and Neumann at both ends.

• UpperDirichlet - Homogeneous Dirichlet at x=1, nothing at x=-1.

• Polar - Homogeneous biharmonic, specifically used with polar coordinates.

• DirichletNeumann - Homogeneous Dirichlet at x=-1 and Neumann at x=1.

• NeumannDirichlet - Homogeneous Neumann at x=-1 and Dirichlet at x=1.

• dtype (str or np.dtype, optional) – The datatype of physical space (input to forward transforms)

• quad (str, optional) – Type of quadrature

• For family=Chebyshev:

• GL - Chebyshev-Gauss-Lobatto

• GC - Chebyshev-Gauss

• For family=Legendre:

• LG - Legendre-Gauss

• GL - Legendre-Gauss-Lobatto

• For family=Laguerre:

• LG - Laguerre-Gauss

• For family=Hermite:

• HG - Hermite-Gauss

• For family=Jacobi:

• JG - Jacobi-Gauss

• domain (two-tuple of floats, optional) – The computational domain

• scaled (bool) – Whether to use scaled basis (only Legendre)

• padding_factor (float, optional) – For padding backward transform (for dealiasing)

• dealias_direct (bool, optional) – Use 2/3-rule dealiasing

• coordinates (2- or 3-tuple (coordinate, position vector (, sympy assumptions)), optional) – Map for curvilinear coordinatesystem. The new coordinate variable in the new coordinate system is the first item. Second item is a tuple for the Cartesian position vector as function of the new variable in the first tuple. Example:

theta = sp.Symbols('x', real=True, positive=True)
rv = (sp.cos(theta), sp.sin(theta))


where theta and rv are the first and second items in the 2-tuple.

Examples

>>> from shenfun import FunctionSpace
>>> F0 = FunctionSpace(16, 'F')
>>> C1 = FunctionSpace(32, 'C', quad='GC')

class shenfun.forms.arguments.TestFunction(space, index=0, basespace=None, offset=0, base=None)[source]

Test function - BasisFunction with argument = 0

Parameters
• space (TensorProductSpace or MixedTensorProductSpace)

• index (int, optional) – Component of basis with rank > 0

• basespace (The base MixedTensorProductSpace if space is a) – subspace.

• offset (int) – The number of scalar spaces (i.e., :class:.TensorProductSpacees) ahead of this space

• base (The base TestFunction)

property argument

Return argument of basis

class shenfun.forms.arguments.TrialFunction(space, index=0, basespace=None, offset=0, base=None)[source]

Trial function - BasisFunction with argument = 1

Parameters
• space (TensorProductSpace or MixedTensorProductSpace)

• index (int, optional) – Component of basis with rank > 0

• basespace (The base MixedTensorProductSpace if space is a) – subspace.

• offset (int) – The number of scalar spaces (i.e., :class:.TensorProductSpacees) ahead of this space

• base (The base TrialFunction)

property argument

Return argument of basis

## shenfun.forms.inner module¶

This module contains the inner function that computes the weighted inner product.

shenfun.forms.inner.inner(expr0, expr1, output_array=None, level=0)[source]

Return (weighted) discrete inner product of linear or bilinear form

$(f, g)_w^N = \sum_{i\in\mathcal{I}}f(x_i) \overline{g}(x_i) w_i \approx \int_{\Omega} g\, \overline{f}\, w\, dx$

where $$\mathcal{I}=0, 1, \ldots, N, N \in \mathbb{Z}^+$$, $$f$$ is an expression linear in a TestFunction, and $$g$$ is an expression that is linear in TrialFunction or Function, or it is simply an Array (a solution interpolated on the quadrature mesh in physical space). $$w$$ is a weight associated with chosen basis, and $$w_i$$ are quadrature weights.

If the expressions are created in a multidimensional TensorProductSpace, then the sum above is over all dimensions. In 2D it becomes:

$(f, g)_w^N = \sum_{i\in\mathcal{I}}\sum_{j\in\mathcal{J}} f(x_i, y_j) \overline{g}(x_i, y_j) w_j w_i$

where $$\mathcal{J}=0, 1, \ldots, M, M \in \mathbb{Z}^+$$.

Parameters
• expr0, expr1 (Expr, BasisFunction, Array) – or number. Either one can be an expression involving a BasisFunction (TestFunction, TrialFunction or Function) an Array or a number. With expressions (Expr) on a BasisFunction we typically mean terms like div(u) or grad(u), where u is any one of the different types of BasisFunction. If one of expr0/expr1 involves a TestFunction and the other one involves a TrialFunction, then a tensor product matrix (or a list of tensor product matrices) is returned. If one of expr0/expr1 involves a TestFunction and the other one involves a Function, or a plain Array, then a linear form is assembled and a Function is returned. If one of expr0/expr1 is a number (typically 1), or a tuple of numbers for a vector, then the inner product represents a non-weighted integral over the domain. If a single number, then the other expression must be a scalar Array or a Function. If a tuple of numbers, then the Array/Function must be a vector.

• output_array (Function) – Optional return array for linear form.

• level (int) – The level of postprocessing for assembled matrices. Applies only to bilinear forms

• 0 Full postprocessing - diagonal matrices to scale arrays and add equal matrices

• 1 Diagonal matrices to scale arrays, but don’t add equal matrices

• 2 No postprocessing, return all assembled matrices

Returns

Function for linear forms.

SparseMatrix for bilinear 1D forms.

TPMatrix or list of TPMatrix for bilinear multidimensional forms.

Number, for non-weighted integral where either one of the arguments is a number.

Return type

Depending on dimensionality and the arguments to the forms

Example

Compute mass matrix of Shen’s Chebyshev Dirichlet basis:

>>> from shenfun import FunctionSpace
>>> from shenfun import TestFunction, TrialFunction
>>> SD = FunctionSpace(6, 'Chebyshev', bc=(0, 0))
>>> u = TrialFunction(SD)
>>> v = TestFunction(SD)
>>> B = inner(v, u)
>>> d = {-2: np.array([-np.pi/2]),
...       0: np.array([ 1.5*np.pi, np.pi, np.pi, np.pi]),
...       2: np.array([-np.pi/2])}
>>> [np.all(abs(B[k]-v) < 1e-7) for k, v in d.items()]
[True, True, True]


## shenfun.forms.operators module¶

This module contains the implementation of operators acting on arguments.

shenfun.forms.operators.Dx(test, x, k=1)[source]

Return k’th order partial derivative in direction x

Parameters
• test (Expr or BasisFunction)

• x (int) – axis to take derivative over

• k (int) – Number of derivatives

shenfun.forms.operators.curl(test)[source]

Return curl of test

Parameters

test (Expr or BasisFunction)

shenfun.forms.operators.div(test)[source]

Return div(test)

Parameters

test (Expr or BasisFunction) – Must be rank > 0 (cannot take divergence of scalar)

shenfun.forms.operators.grad(test)[source]

Parameters

test (Expr or BasisFunction)

Note

Increases the rank of Expr by one

## shenfun.forms.project module¶

shenfun.forms.project.project(uh, T, output_array=None, fill=True, use_to_ortho=True, use_assign=True)[source]

Project uh to tensor product space T

Find $$u \in T$$, such that

$(u - u_h, v)_w^N = 0 \quad \forall v \in T$
Parameters
• uh (Instance of either one of) –

• output_array (Function, optional) – Return array

• fill (bool, optional) – Whether to fill the output_array with zeros before projection

• use_to_ortho (bool, optional) – Whether to use fast to_ortho method for projection of Functions to orthogonal space.

• use_assign (bool, optional) – Whether to use fast assign method for projection of Function to a denser space of the same kind.

Returns

The projection of uh in T

Return type

Function

Example

>>> import numpy as np
>>> from mpi4py import MPI
>>> from shenfun import FunctionSpace, project, TensorProductSpace, Array, \
...     Function, Dx
>>> N = 16
>>> comm = MPI.COMM_WORLD
>>> T0 = FunctionSpace(N, 'C')
>>> K0 = FunctionSpace(N, 'F', dtype='d')
>>> T = TensorProductSpace(comm, (T0, K0))
>>> uj = Array(T)
>>> uj[:] = np.random.random(uj.shape)
>>> u = Function(T)
>>> u = project(uj, T, output_array=u) # Same as u = T.forward(uj, u)
>>> du = project(Dx(u, 0, 1), T)


## Module contents¶

shenfun.forms.extract_bc_matrices(mats)[source]

Extract boundary matrices from list of mats

Parameters

mats (list of list of :class:.TPMatrixes)

Returns

list of boundary matrices.

Return type

list

Note

The mats list is modified in place since boundary matrices are extracted.