Assembling custom forms

pyiga makes it possible to define custom linear or bilinear forms and have highly efficient code for assembling them automatically generated, similar to the facilities in FEniCS or NGSolve.

Assembling variational forms from their textual description

The most convenient interface for assembling custom forms is by calling pyiga.assemble.assemble() with a string describing the problem as the first argument. Here is a simple but complete example:

from pyiga import bspline, geometry, assemble

kvs = (bspline.make_knots(3, 0.0, 1.0, 30),
       bspline.make_knots(3, 0.0, 1.0, 30))
geo = geometry.quarter_annulus()
A = assemble.assemble('inner(grad(u), grad(v)) * dx', kvs, geo=geo)

Here we define a cubic 2D tensor product spline space with 30 intervals per coordinate direction and then assemble a standard Laplace variational form over a quarter annulus geometry. The result is the sparse stiffness matrix A.

The string representing the integrand occurring in this example should be mostly self-explanatory: it represents the inner product of the gradients of the trial and test functions u and v, multiplied by the volume measure dx.

In a similar way we can assemble linear functionals, such as we would need for computing a right-hand side for a discretized problem:

def f(x, y): return np.cos(x) * np.exp(y)
rhs = assemble.assemble('f * v * dx', kvs, geo=geo, f=f)

Note that we specified the additional input function f as a keyword argument to assemble(). The assemble() function automatically determined that we are assembling a linear functional from the fact that we are referring to only one basis function, v, in the description of our form. The result is not a 1D vector as one might expect, but a 2D array whose shape is

rhs.shape == (kvs[0].numdofs, kvs[1].numdofs)

However, it is easy to convert it to a 1D vector by calling rhs.ravel().

Vector-valued problems are realized by informing assemble() of the number of components of the basis functions via the bfuns argument. For instance, the stiffness matrix for the vector Laplace problem could be realized as follows:

A = assemble.assemble('inner(grad(u), grad(v)) * dx', kvs,
        bfuns=[('u',2), ('v',2)], geo=geo)

Here grad() computes the Jacobian matrices of the vector-valued basis functions u and v, and inner() computes their Frobenius product.

In addition to the Supported functions described below, the following keywords can be used in the string description:

dx volume measure
ds surface measure (switches to Surface integrals)
u, v default trial and test function (only when not using the bfuns argument of assemble())
x current position in physical coordinates (i.e., the value of the geometry map)
n the unit normal vector (only for surface integrals)
gw the Gauss quadrature weight at the current node (usually not needed since it is included in dx and ds)
jac the Jacobian matrix of the geometry map with shape geo_dim x dim

The assemble() function also works in hierarchical spline spaces; simply pass an instance of pyiga.hierarchical.HSpace as the second argument instead of the tuple of knot vectors kvs. The result will be a sparse matrix or a 1D vector, and in each case the degrees of freedom are ordered according to the canonical order defined in the pyiga.hierarchical module.

For most problems, this simple string-based interface is sufficient to generate all required matrices and vectors. However, in some edge cases it may be necessary to define a custom vform by hand. The remainder of this document describes the internals of how to create such objects. It also gives a more formal description of how the expressions occurring in the strings above should be interpreted. In essence, these strings are nothing but Python expressions evaluated in a context where the members of the pyiga.vform module are made available.

Programmatically defining VForms

The pyiga.vform module contains the tools for describing variational forms, and pyiga.vform.VForm is the main class used to represent abstract variational forms. Its constructor has one required argument which describes the space dimension. For instance, to initialize a VForm for a three-dimensional problem:

from pyiga import vform

vf = vform.VForm(3)

In order to create expressions for our form, we need objects which represent the functions that our form operates on. By default VForm assumes a bilinear form, and therefore we can obtain objects for the trial and the test function using VForm.basisfuns() like this:

u, v = vf.basisfuns()

The objects that we work with when specifying vforms are abstract expressions (namely, instances of Expr) and all have certain properties such as a shape, Expr.shape, which is a tuple of dimensions just like a numpy array has. By default, VForm assumes a scalar-valued problem, and therefore both the trial function u and the test function v are scalar:

>>> u.shape
()
>>> v.shape
()

We can now start building expressions using these functions. Let’s first import some commonly needed functions from the vform module.

from pyiga.vform import grad, div, inner, dx

We will often need the gradient of a function, obtained via grad():

>>> gu = grad(u)
>>> gu.shape
(3,)

Note that grad(u) is itself an expression. As expected, the gradient of a scalar function is a three-component vector. We could take the divergence (div()) of the gradient and get back a scalar expression which represents the Laplacian \(\Delta u = \operatorname{div} \nabla u\) of u:

>>> Lu = div(gu)
>>> Lu.shape
()

However, in order to express the standard variational form for the Laplace problem, we only require the inner product \(\nabla u \cdot \nabla v\) of the gradients of our input functions, which can be computed using inner():

>>> x = inner(grad(u), grad(v))
>>> x.shape
()

Again, this is a scalar since inner() represents a contraction over all axes of its input tensors; for vectors, it is the scalar product, and for matrices, it is the Frobenius product.

Note

In general, the syntax for constructing forms sticks as closely as possible to that of the UFL language used in FEniCS, and therefore the UFL documentation is also a helpful resource.

Finally we want to represent the integral of this expression over the computational domain. We do this by multiplying with the symbol dx:

integral = inner(grad(u), grad(v)) * dx

Internally, dx is actually a scalar expression which represents the absolute value of the determinant of the geometry Jacobian, i.e., the scalar term that stems from transforming the integrand from the physical domain to the parameter domain.

We are now ready to add this expression to our VForm via VForm.add(), and since the expression is rather simple, we can skip all the intermediate steps and variables and simply do

vf.add(inner(grad(u), grad(v)) * dx)

Note that the expression passed to VForm.add() here is exactly the string we passed to assemble() in the first example in the previous section.

A simple example

It’s usually convenient to define vforms in their own functions so that we don’t pollute the global namespace with the definitions from the vform module. The Laplace variational form

\[a(u,v) = \int_\Omega \nabla u \cdot \nabla v \, dx\]

would be defined like this:

def laplace_vf(dim):
    from pyiga.vform import VForm, grad, inner, dx
    vf = VForm(dim)
    u, v = vf.basisfuns()
    vf.add(inner(grad(u), grad(v)) * dx)
    return vf

Calling this function results in a VForm object:

>>> laplace_vf(2)
<pyiga.vform.VForm at 0x7f0fdcf5c0f0>

Note

Currently, the predefined Laplace variational form in pyiga is defined in a different way which leads to a slightly higher performance of the generated code.

Vector-valued problems

By default, basis functions are assumed to be scalar-valued. To generate forms with vector-valued functions, simply pass the components argument with the desired sizes to VForm.basisfuns():

>>> vf = vform.VForm(2)
>>> u, v = vf.basisfuns(components=(2,2))

>>> u.shape, v.shape
((2,), (2,))

We can still compute gradients (Jacobians) using grad() as before:

>>> grad(u).shape
(2, 2)

As a simple example, the div-div bilinear form \(a(u,v) = \int_\Omega \operatorname{div} u \operatorname{div} v \,dx\) would be implemented using

vf.add(div(u) * div(v) * dx)

It is also possible to mix vector and scalar functions, e.g. for Stokes-like problems:

vf = vform.VForm(2)
u, p = vf.basisfuns(components=(2,1))

vf.add(div(u) * p * dx)

In this example, u is a vector-valued function and p is scalar-valued.

Working with coefficient functions

Often you will need to provide additional functions as input to your assembler, for instance to represent a diffusion coefficient which varies over the computational domain. A scalar input field is declared using the VForm.input() method as follows:

>>> vf = vform.VForm(2)
>>> coeff = vf.input('coeff')

>>> coeff.shape
()

The new variable coeff now represents a scalar expression that we can work with just as with the basis functions, e.g.

>>> grad(coeff).shape
(2,)

As a simple example, to use this as a scalar diffusion coefficient, we would do

u, v = vf.basisfuns()
vf.add(coeff * inner(grad(u), grad(v)) * dx)

Input fields can be declared vector- or matrix-valued simply by prescribing their shape. In this example, we declare a 2x2 matrix-valued coefficient function:

>>> vf = vform.VForm(2)
>>> coeff = vf.input('coeff', shape=(2,2))

>>> coeff.shape
(2, 2)

To actually supply these functions to the assembler, they must be passed as keyword arguments to the constructor of the generated assembler class. It is possible to pass either standard Python functions (in which case differentiation is not supported) or instances of pyiga.bspline.BSplineFunc or pyiga.geometry.NurbsFunc. In fact, the predefined input geo for the geometry map is simply declared as a vector-valued input field. See the section Compiling and assembling for an example of how to pass these functions.

By default, input functions are considered to be defined in the coordinates of the parameter domain. If your input function is given in terms of physical coordinates, declare it as follows:

coeff = vf.input('coeff', physical=True)

Note

In the simple string-based interface described in the first section, functions passed as BSplineFunc or similar objects are assumed to be given in parametric coordinates, whereas standard Python functions are assumed to be given in physical coordinates. This simple heuristic usually produces the expected result.

For performance reasons, it is sometimes beneficial to be able to update a single input function without recreating the entire assembler class, for instance when assembling the same form many times with different coefficients in a Newton iteration for a nonlinear problem. In this case, we can declare the function as follows:

func = vf.input('func', updatable=True)

The generated assembler class then has an update() method which takes the function as a keyword argument and updates it accordingly, e.g.,

asm.update(func=F)

Defining constant values

If a needed coefficient function is constant, it is unnecessary to use the VForm.input() machinery. Instead, we can simply define constant values using the as_expr(), as_vector(), and as_matrix() functions as follows:

>>> coeff = vform.as_expr(5)
>>> coeff.shape
()

>>> vcoeff = vform.as_vector([2,3,4])
>>> vcoeff.shape
(3,)

>>> mcoeff = vform.as_matrix([[2,1],[1,2]])
>>> mcoeff.shape
(2, 2)

We can then work with these constants exactly as with any other expression.

For constant scalar values as well as tuples of constants or expressions, the coercion to expressions is performed implicitly, making as_expr() and as_vector() unnecessary in these cases. This means that we can directly write expressions such as

vf.add(inner(3 * grad(u), grad(v)) * dx)
vf.add(inner((2.0, 3.0), grad(u)) * dx)

The first example also shows that multiplication of a scalar with a vector works as expected, i.e., the vector is multiplied componentwise with the scalar.

The above approach has the disadvantage that the assembler needs to be recompiled every time a constant changes. It is possible to instead define constant parameters which are unknown at compile-time and specified only when assembling. Such constants are defined, analogously to the input fields in the previous section, using the VForm.parameter() method:

>>> vf = vform.VForm(2)
>>> b = vf.parameter('b', shape=(2,))

>>> b.shape
(2,)

Again, scalar, vector-valued and matrix-valued parameters are supported. Before assembling, these parameters must be set using the update_params() method of the assembler class, analogously to the update() method for input fields described in the previous section.

When using the string-based interface, such constants may be defined simply by passing them as keyword arguments into the assemble function:

f = assemble.assemble('inner(b, grad(v)) * dx',
                      kvs, geo=geo, b=(2.0, -1.0))

Defining linear (unary) forms

By default, VForm assumes the case of a bilinear form, i.e., having a trial function u and a test function v. For defining right-hand sides, we usually need linear forms which have only one argument. We can do this by passing arity=1 to the VForm constructor. The VForm.basisfuns() method returns only a single basis function in this case.

Below is a simple example for defining the linear form \(\langle F,v \rangle = \int_\Omega f v \,dx\) with a user-specified input function f:

vf = vform.VForm(3, arity=1)
v = vf.basisfuns()
f = vf.input('f')
vf.add(f * v * dx)

Working with parametric derivatives

By default, a VForm assumes that you will provide it with a geometry map under the input field name geo and automatically transforms the derivatives of the basis functions u and v, as well as gradients of any input fields, accordingly. If for some reason you need to work with untransformed gradients with respect to parametric coordinates, you can simply pass the keyword argument parametric=True to the derivative functions such as grad() like this:

vf = vform.VForm(2)
u, v = vf.basisfuns()
f = vf.input('f')
gu = grad(u, parametric=True)
gf = Dx(f, 1, parametric=True)

You can compute both parametric and physical derivatives of basis functions as well as input fields given in parametric coordinates. An input field that is given in terms of physical coordinates only supports physical derivatives.

Note that the symbol dx still includes the geometry Jacobian, and therefore you should not multiply your expression with it if you want to integrate over the parameter domain instead of the physical domain. In this case, you should multiply your expression with the attribute VForm.GaussWeight instead, which represents the weight for the Gauss quadrature. When using the textual description, use the keyword gw for this purpose. When computing integrals over the physical domain, the quadrature weight is automatically subsumed into dx and does not need to be specified explicitly.

Surface integrals

By default, a VForm will assume that the dimension of the image of the geometry map is the same as the dimension of the spline space over which we are integrating. For computing matrices and vectors over surfaces, we can specify the geo_dim argument of the VForm constructor to be one higher than the input dimension. Of course, the geo function we pass when assembling must match that output dimension. We also have to multiply with the surface measure ds instead of the volume measure dx when computing such integrals.

Here is a simple example which describes a linear functional over a surface:

def L2_surface_functional_vf(dim):
    V = VForm(dim, geo_dim=dim+1, arity=1)
    v = V.basisfuns()
    f = V.input('f')
    V.add(f * v * ds)
    return V

When called with dim=2, it represents a 2D surface integral in a 3D ambient space.

If you need to use the normal vector in your expression, you can access it via the VForm.normal attribute. It is oriented according to the standard right-hand rule and has unit length.

At the moment, transformations of derivatives on surfaces are not implemented, and therefore you can only use parametric derivatives.

The string-based interface described above will automatically switch to surface integration when it detects that ds was used instead of dx.

Supported functions

The following functions and expressions are implemented in pyiga.vform and have the same semantics as in the UFL language (see the UFL documentation):

dx Dx() grad() div() curl() as_vector() as_matrix() inner() dot() tr() det() inv() cross() outer() abs() sqrt() exp() log() sin() cos() tan()

In addition, all expressions have the members Expr.dx() for partial derivatives (analogous to the global function Dx()), Expr.T for transposing matrices, and Expr.dot() which is analogous to the global dot() function. Vector and matrix expressions can be indexed and sliced using the standard Python [] operator. Expressions also support the standard arithmetic operators +, -, *, /, **.

Compiling and assembling

Once a vform has been defined, it has to be compiled into efficient code and then invoked in order to assemble the desired matrix or vector. Currently, there is only one backend in pyiga which is based on Cython – the vform is translated into Cython code, compiled on the fly and loaded as a dynamic module. The compiled modules are cached so that compiling a given vform for a second time does not recompile the code.

Below is an example for assembling the Laplace variational form defined in the section A simple example:

from pyiga import assemble, bspline, geometry

# define the trial space
kv = bspline.make_knots(3, 0.0, 1.0, 20)
kvs = (kv, kv)   # 2D tensor product spline space

# define the geometry map
geo = geometry.quarter_annulus()   # NURBS quarter annulus

A = assemble.assemble_vf(stiffness_vf(2), kvs, geo=geo, symmetric=True)

Any further input functions the assembler requires can be passed as further keyword arguments in the call to assemble_vf(). The function will automatically detect whether the VForm has arity 1 or 2 and generate a vector or a matrix correspondingly.

Manual compilation of the variational form

Sometimes it may be necessary to directly work with the assembler class that results from compiling a VForm. The functions used for compilation are contained in the pyiga.compile module, and the resulting matrices can be computed using the pyiga.assemble.assemble_entries() functions.

Using these functions, the Laplace variational form defined above can be assembled as follows:

from pyiga import compile, assemble, bspline, geometry

# compile the vform into an assembler class
Asm = compile.compile_vform(laplace_vf(2))

# define the trial space
kv = bspline.make_knots(3, 0.0, 1.0, 20)
kvs = (kv, kv)   # 2D tensor product spline space

# define the geometry map
geo = geometry.quarter_annulus()   # NURBS quarter annulus

A = assemble.assemble_entries(Asm(kvs, geo=geo), symmetric=True)

The geometry map is passed using geo= to the constructor of the assembler class, and further input functions defined as described in the section Working with coefficient functions can be passed in the same way using their given name as the keyword.

The resulting object A is a sparse matrix in CSR format; different matrix formats can be chosen by passing the format= keyword argument to assemble_entries(). The argument symmetric=True takes advantage of the symmetry of the bilinear form in order to speed up the assembly.