shenfun.forms package

Submodules

shenfun.forms.arguments module

class shenfun.forms.arguments.Array(space, val=0, buffer=None, reltol=1e-12, abstol=1e-15)[source]

Bases: ShenfunBaseArray

Numpy array for TensorProductSpace

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

The Function represents in 1D

\[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:

  • space (TensorProductSpace or SpectralBase)

  • 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

Examples

>>> from shenfun import comm, FunctionSpace, TensorProductSpace, Array
>>> K0 = FunctionSpace(8, 'F', dtype='D')
>>> K1 = FunctionSpace(8, 'F', dtype='d')
>>> FFT = TensorProductSpace(comm, [K0, K1])
>>> u = Array(FFT)
forward(output_array=None, kind=None)[source]

Return Function used to evaluate Array

Parameters:

  • output_array (array, optional) – Expansion coefficients

  • kind (str, optional) –

    • ‘fast’ - use fast transform if implemented

    • ‘vandermonde’ - Use Vandermonde matrix

property forward_output

Return whether self is the result of a forward transform

get_cartesian_vector()[source]

Return self as a Cartesian vector.

Note

This method is only relevant for curvilinear coordinates, where the vector uses different basis functions from the Cartesian.

offset()[source]

Return offset of this basis

The offset is the number of scalar TensorProductSpace es ahead of this Arrays space in a CompositeSpace.

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

Bases: object

Base class for arguments to shenfun’s Expr

Parameters:
property argument

Return argument of basis

property base

Return base

property dimensions

Return dimensions of function space

expr_rank()[source]

Return rank of 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 TensorProductSpace es ahead of this space in a CompositeSpace.

property tensor_rank

Return tensor 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:

  • basis (BasisFunction) – TestFunction, TrialFunction or Function

  • 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 ndim, which is the number of dimensions, and for second order tensors it equals ndim*ndim.

    • 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 CompositeSpace. Only used when the basis of the form is composite. 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

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

get_contravariant_component(i, j=None)[source]

Return contravariant vector component

Parameters:

  • i (int or None) – First vector component. In 1D set i to None

  • j (int, optional) – Second component of tensor

Note

The contravariant vector components are the components to the covariant basisvectors. If the basis vectors are normalized, i.e., if the setting is config[‘basisvectors’] == ‘normal’, then the normal vector components must be scaled in order to get to the contravariant components.

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

Return indices of Expr

property is_basis_function
num_components()[source]

Return number of components in Expr

num_terms()[source]

Return number of terms in Expr

scales()[source]

Return scales of Expr

set_basis(u)[source]

Change the basis function this Expr applies to

simplify()[source]

Join terms, that are otherwise equal, using scale

subs(a, b)[source]

Replace a with b in scales

Parameters:

  • a (Sympy Symbol)

  • b (Number)

property tensor_rank

Return rank of Expr’s BasisFunction

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, reltol=1e-12, abstol=1e-15)[source]

Bases: ShenfunBaseArray, BasisFunction

Spectral Galerkin function for given TensorProductSpace or SpectralBase

The Function in one dimension is

\[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 or FunctionSpace())

  • 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.

  • reltol (number, optional) – Relative tolerance for adaptively finding dimension of space from the buffer.

  • abstol (number, optional) – Absolute tolerance for adaptively finding dimension of space from the buffer

Note

For more information, see numpy.ndarray

Example

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

You get the array the space is planned for. This is probably what you expect:

>>> u = Function(K0)
>>> print(u.shape)
(8,)

But make the TensorProductSpace change K0 and K1 inplace:

>>> T = TensorProductSpace(comm, [K0, K1], modify_spaces_inplace=True)
>>> u = Function(K0)
>>> print(u.shape)
(8, 5)

If you want to preserve K0 as a 1D function space, then use modify_spaces_inplace=False, which is the default behaviour.

Note

The array returned will have the same shape as the arrays space is planned for. So if you want a Function on a 1D FunctionSpace, then make sure that FunctionSpace is not planned for a TensorProductSpace.

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, kind=None, mesh=None, padding_factor=None)[source]

Return Function evaluated on some quadrature mesh

Parameters:

  • output_array (array, optional) – Return array, Function evaluated on mesh

  • kind (dict or str, optional) – Can be used to overload config[‘transforms’][‘kind’], using a dictionary with family as key and values either one of the strings - ‘fast’ - ‘recursive’ - ‘vandermonde’ For example kind={‘chebyshev’: ‘vandermonde’} Note that for one-dimensional problems one can use just the string value and no dict

  • mesh (str, list or functionspace, optional) –

    • ‘quadrature’ - use quadrature mesh of self

    • ‘uniform’ - use uniform mesh

    • A function space with the same mesh distribution as self

  • padding_factor (None or number or tuple, optional) – If padding_factor is a number different from 1, then perform a padded backward transform, using a padded work array ‘self._backward_work_array’. Use a tuple for different padding in different directions.

Note

If output_array is not provided, then a new Array is created and returned. If padding_factor is used, then the returned Array is cached and reused for subsequent calls using the same padding_factor.

property base

Base object if memory is from some other object.

Examples

The base of an array that owns its memory is None:

>>> x = np.array([1,2,3,4])
>>> x.base is None
True

Slicing creates a view, whose memory is shared with x:

>>> y = x[2:]
>>> y.base is x
True
copy_from_flattened(f, j=(), dims=None, sl=None)[source]

Copy dofs to self from a flattened array

Parameters:

  • f (array 1D) – The array to copy from

  • j (int or tuple, optional) – Index into diagonal axes (Fourier)

  • dims (array 1D, optional) – The flat dim of each variable in the composite basis

  • sl (sequence of slices, optional) – list of slices into non-diagonal axes

Note

All Functions have allocated Numpy arrays assuming all dofs are being used. However, for non-orthogonal bases with boundary conditions we use only the first N-nb dofs, where nb is the number of boundary conditions. This function does not copy the boundary dofs.

copy_to_flattened(f=None, j=(), dims=None, sl=None)[source]

Copy dofs of self to a flattened array

Parameters:

  • f (array 1D, optional) – The array to copy to

  • j (int or tuple, optional) – Index into diagonal axes (Fourier). Not used if sl is provided

  • dims (array 1D, optional) – The dim of each variable in the composite basis

  • sl (sequence of slices, optional) – list of slices into non-diagonal axes

Note

All Functions have allocated Numpy arrays assuming all dofs are being used. However, for non-orthogonal bases with boundary conditions we use only the first N-nb dofs, where nb is the number of boundary conditions. This function does not copy the boundary dofs.

eval(x, output_array=None)[source]

Evaluate Function at points x

Parameters:

  • x (float or array of floats)

  • 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

get_dealiased_space(padding_factor)[source]

Return a function space to be used for padded transforms

Parameters:

padding_factor (number or tuple) – The padding factor of the backward transform A tuple can be used for padding differently in each direction.

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, basis=None, 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,

    • Chebyshevu or U,

    • Legendre or L,

    • Fourier or F,

    • Laguerre or La,

    • Hermite or H

    • Jacobi or J

    • Ultraspherical or Q

  • bc (tuple or dict, optional) – See BoundaryConditions

  • 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=Chebyshevu:

      • GU - Chebyshevu-Gauss

      • 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

    • For family=Ultraspherical:

      • QG - Jacobi-Gauss

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

  • scaled (bool) – Whether to use scaled basis

  • basis (str) – Name of basis to use, if there are more than one possible basis for a given boundary condition. For example, there are three Dirichlet bases for the Chebyshev family: Heinrichs, ShenDirichlet and Phi1

  • 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 coordinates. Used for curves parametrized with one parameter in \(\mathbb{R}^2\) or \(\mathbb{R}^3\). The parameter is the first item. Second item is a tuple of length 2 or 3 (for \(\mathbb{R}^2\) or \(\mathbb{R}^3\)) for the Cartesian position vector. The sympy assumptions are there to help Sympy compute basis vectors and metrics, see Coordinates for options. Example circle of radius 1:

    theta = sp.Symbol('x', real=True)
rv = (sp.cos(theta), sp.sin(theta))
FunctionSpace(16, 'L', coordinates=(theta, rv), domain=(0, 2*np.pi))

Examples

>>> from shenfun import FunctionSpace
>>> F0 = FunctionSpace(16, 'F')
>>> C1 = FunctionSpace(32, 'C', quad='GC')
class shenfun.forms.arguments.TestFunction(space, index=0, offset=0, base=None)[source]

Bases: BasisFunction

Test function - BasisFunction with argument = 0

Parameters:
property argument

Return argument of basis

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

Bases: BasisFunction

Trial function - BasisFunction with argument = 1

Parameters:
property argument

Return argument of basis

shenfun.forms.inner module

class shenfun.forms.inner.Inner(v, uh)[source]

Bases: object

Return an instance of a class that can perform the inner product of the linear form efficiently through a matrix-vector product

Parameters:

Note

This is an optimization only for linear forms, not bilinear. There is no need to use this class for regular scalar products, where uh is simply an Array.

shenfun.forms.inner.inner(expr0, expr1, output_array=None, assemble=None, kind=None, fixed_resolution=None, return_matrices=False)[source]

Return (weighted or unweighted) discrete inner product

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

where \(\mathcal{I}=0, 1, \ldots, N-1, N \in \mathbb{Z}^+\), \(f\) is a number or 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.

Note

If \(f\) is a number (typically one) and \(g\) an Array, then inner represents an unweighted, regular integral over the domain.

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, 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.

  • assemble (None or str, optional) – Determines how to perform the integration

    • ‘quadrature’ (default)

    • ‘exact’

    • ‘adaptive’

    Exact and adaptive should result in the same matrix. Exact computes the integral using Sympy integrate, whereas adaptive makes use of adaptive quadrature through scipy. For most bilinear forms quadrature will also return the same matrix as exact. With non-constant measures (like for curvilinear coordinates) there is more likely to be differences.

  • kind (None, str or dict, optional) – Alternative methods.

    For bilinear forms (matrices) kind is a string:

    • ‘implemented’ - Hardcoded implementations

    • ‘stencil’ - Use orthogonal bases and stencil-matrices

    • ‘vandermonde’ - Use Vandermonde matrix

    The default (for kind=None) is to first try to look for implemented kind, and if that fails try first ‘stencil’ and then finally fall back on vandermonde. Vandermonde creates a dense matrix of size NxN, so it should be avoided (e.g., by implementing the matrix) for large N.

    For linear forms (vectors) the kind keyword can be used to overload the default methods for transforms set in config[‘transforms’][‘kind’]. Hence, kind is a dictionary with family as key and values either one of the possible methods

    • ‘fast’ - Use FFT (only Fourier and Chebyshev first and second kind)

    • ‘recursive’ - Low-memory implementation (only for Jacobi polynomials)

    • ‘vandermonde’ - Use Vandermonde matrix

    E.g., kind={‘chebyshev’: ‘recursive’}. Note that for one-dimensional problems it is enough to use just the value of the dictionary.

  • fixed_resolution (Number or sequence of integers, optional) – A fixed number of quadrature points used to compute the inner product along each dimension of the domain. If ‘fixed_resolution’ is set, then assemble is set to ‘quadrature’ and kind is set to ‘vandermonde’. fixed_resolution is argument to TensorProductSpace.get_refined().

  • return_matrices (bool, optional) – For linear forms, whether to simply return the matrices that are used to compute the result with matrix-vector products.

Note

For most matrices all methods will lead to the same result. For bilinear forms with polynomial coefficients regular quadrature will become inaccurate and a ‘fixed_resolution’ higher than the regular number of quadrature points should be considered.

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

See also

project()

Example

Compute mass matrix of Shen’s Chebyshev Dirichlet basis:

>>> from shenfun import FunctionSpace, TensorProductSpace
>>> from shenfun import TestFunction, TrialFunction, Array
>>> 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]

# Compute unweighted integral >>> F = FunctionSpace(10, ‘F’, domain=(0, 2*np.pi)) >>> T = TensorProductSpace(comm, (SD, F)) >>> area = inner(1, Array(T, val=1)) >>> print(‘Area of domain =’, area) Area of domain = 12.56637061435917

shenfun.forms.operators module

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]

Return grad(test)

Parameters:

test (Expr or BasisFunction)

Note

Increases the rank of Expr by one

shenfun.forms.project module

class shenfun.forms.project.Project(uh, T, output_array=None)[source]

Bases: object

Return an instance of a class that can perform projections efficiently

Parameters:

Example

>>> from shenfun import TestFunction, Function, Dx, Project, FunctionSpace
>>> T = FunctionSpace(8, 'C')
>>> u = Function(T)
>>> u[1] = 1
>>> dudx = Project(Dx(u, 0, 1), T)
>>> dudx()
Function([1., 0., 0., 0., 0., 0., 0., 0.])

Note that u[1] = 1 sets the coefficient of the second Chebyshev polynomial (i.e., x) to 1. Hence the derivative of u is dx/dx=1, which is the result of dudx().

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) –

  • T (SpectralBase, TensorProductSpace or CompositeSpace)

  • 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

See also

inner()

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