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