Examples

These examples are complete Python scripts that can be run on their own. The documentation workflow runs all of them before building the HTML pages, ensuring that their imports and demonstrated APIs remain current.

Variational forms

This example constructs a scalar Poisson problem and a Nitsche-type form on a two-patch domain.

 1"""Construct scalar variational forms on single- and multi-patch domains."""
 2
 3from sympy import cos
 4
 5from sympde.api import BilinearForm, Constant, Domain, LinearForm, ScalarFunctionSpace
 6from sympde.api import Square, Dn, dot, elements_of, find, grad, integral
 7from sympde.api import jump, latex, minus, plus
 8
 9
10# A Poisson problem on a single square.
11domain = Square('Omega')
12space = ScalarFunctionSpace('V', domain)
13trial, test = elements_of(space, names='u, v')
14
15bilinear_form = BilinearForm(
16    (trial, test),
17    integral(domain, dot(grad(trial), grad(test))),
18)
19
20x, y = domain.coordinates
21linear_form = LinearForm(test, integral(domain, cos(x + y) * test))
22equation = find(
23    trial,
24    forall=test,
25    lhs=bilinear_form(trial, test),
26    rhs=linear_form(test),
27)
28
29print(latex(equation))
30
31
32# A Nitsche-type interface contribution on two joined patches.
33patch_a = Square('A')
34patch_b = Square('B')
35domain = Domain.join(
36    [patch_a, patch_b],
37    [((0, 0, 1), (1, 0, -1), 1)],
38    'Omega',
39)
40
41space = ScalarFunctionSpace('V', domain)
42trial, test = elements_of(space, names='u, v')
43interface = domain.interfaces
44kappa = Constant('kappa')
45
46interface_integrand = (
47    -jump(trial) * jump(Dn(test))
48    + kappa * jump(trial) * jump(test)
49    + plus(Dn(trial)) * minus(test)
50    + minus(Dn(trial)) * plus(test)
51)
52
53multi_patch_form = BilinearForm(
54    (trial, test),
55    integral(domain, dot(grad(trial), grad(test)))
56    + integral(interface, interface_integrand),
57)
58
59print(latex(multi_patch_form))

Differential calculus

This example demonstrates the currently supported symbolic vector-calculus identities.

 1"""Verify representative symbolic vector-calculus identities."""
 2
 3from sympy import expand
 4
 5from sympde.api import Constant, Domain, ScalarFunctionSpace, VectorFunctionSpace
 6from sympde.api import curl, div, elements_of, grad, laplace
 7
 8
 9domain = Domain('Omega', dim=3)
10scalar_space = ScalarFunctionSpace('V', domain)
11vector_space = VectorFunctionSpace('W', domain)
12
13alpha, beta = [Constant(name) for name in ('alpha', 'beta')]
14f, g, h = elements_of(scalar_space, names='f, g, h')
15field_f, field_g, field_h = elements_of(vector_space, names='F, G, H')
16
17# Gradient linearity and product rules.
18assert grad(f + g) == grad(f) + grad(g)
19assert grad(alpha * f + beta * g) == alpha * grad(f) + beta * grad(g)
20assert grad(f * g) == f * grad(g) + g * grad(f)
21assert grad(f / g) == -f * grad(g) / g**2 + grad(f) / g
22assert expand(grad(f * g * h)) == (
23    f * g * grad(h) + f * h * grad(g) + g * h * grad(f)
24)
25
26# Linearity of curl, Laplace, and divergence operators.
27assert curl(field_f + field_g) == curl(field_f) + curl(field_g)
28assert curl(alpha * field_f + beta * field_g) == (
29    alpha * curl(field_f) + beta * curl(field_g)
30)
31assert laplace(f + g) == laplace(f) + laplace(g)
32assert laplace(alpha * f + beta * g) == alpha * laplace(f) + beta * laplace(g)
33assert div(field_f + field_g) == div(field_f) + div(field_g)
34assert div(alpha * field_f + beta * field_g) == (
35    alpha * div(field_f) + beta * div(field_g)
36)
37
38# Exact-sequence identities.
39assert curl(grad(h)) == 0
40assert div(curl(field_h)) == 0

Analytical mappings

This example evaluates a built-in analytical mapping and defines a custom mapping from coordinate expressions.

 1"""Work with built-in and user-defined analytical mappings."""
 2
 3import numpy as np
 4from sympy import pi, sin, symbols
 5
 6from sympde.api import CollelaMapping2D, Constant, Domain, LogicalExpr, Mapping
 7
 8
 9# Inspect a built-in mapping symbolically.
10mapping = CollelaMapping2D('M', dim=2)
11mapped_domain = mapping(Domain('Omega', dim=2))
12x1, x2 = symbols('x1, x2', real=True)
13eps, k1, k2 = [Constant(name) for name in ('eps', 'k1', 'k2')]
14
15expected_x = 2.0 * (
16    x1 + eps * sin(2.0 * pi * k1 * x1) * sin(2.0 * pi * k2 * x2)
17) - 1.0
18assert LogicalExpr(mapping[0], mapped_domain) == expected_x
19
20
21# Fix its parameters and evaluate it numerically.
22numerical_mapping = CollelaMapping2D(
23    'M_numeric',
24    dim=2,
25    eps=0.1,
26    k1=1.0,
27    k2=1.0,
28).get_callable_mapping()
29np.testing.assert_allclose(numerical_mapping(0.25, 0.5), (-0.5, 0.0))
30
31
32# Define a custom mapping by specifying its coordinate expressions.
33class ComplexSquareMapping(Mapping):
34    """
35    Map (x1, x2) to (x, y) such that x + i*y = A * (x1 + i*x2)**2.
36
37    Expanding the square gives x = A * (x1**2 - x2**2) and y = 2 * A * x1 * x2.
38    """
39
40    _expressions = {
41        'x': 'A * (x1**2 - x2**2)',
42        'y': '2 * A * x1 * x2',
43    }
44
45
46custom_mapping = ComplexSquareMapping('F', dim=2, A=1.0).get_callable_mapping()
47np.testing.assert_allclose(custom_mapping(0.5, 0.25), (0.1875, 0.25))

Exterior calculus

This example demonstrates differential forms, exterior derivatives, adjoint operators, products, contractions, and Lie derivatives.

 1"""Apply exterior-calculus operators to symbolic differential forms."""
 2
 3from sympy import Symbol
 4
 5from sympde.api import Constant, DifferentialForm, Domain, VectorFunctionSpace
 6from sympde.api import d, delta, element_of, hodge, ip, jp, latex, ld, wedge
 7
 8
 9dimension = Symbol('n')
10alpha = Constant('alpha')
11
12u_0, v_0 = [
13    DifferentialForm(name, index=0, dim=dimension)
14    for name in ('u_0', 'v_0')
15]
16u_1, v_1 = [
17    DifferentialForm(name, index=1, dim=dimension)
18    for name in ('u_1', 'v_1')
19]
20u_2 = DifferentialForm('u_2', index=2, dim=dimension)
21u_n = DifferentialForm('u_n', index=dimension, dim=dimension)
22
23# The exterior derivative is linear, nilpotent, and zero on top forms.
24assert d(d(u_0)) == 0
25assert d(u_0 + v_0) == d(u_0) + d(v_0)
26assert d(alpha * u_0 + v_0) == alpha * d(u_0) + d(v_0)
27assert d(u_n) == 0
28
29# Its adjoint is linear and zero on zero-forms.
30assert delta(u_0) == 0
31assert delta(u_1 + v_1) == delta(u_1) + delta(v_1)
32assert delta(alpha * u_1 + v_1) == alpha * delta(u_1) + delta(v_1)
33
34domain = Domain('Omega', dim=3)
35vector_space = VectorFunctionSpace('W', domain)
36beta = element_of(vector_space, name='beta')
37
38expressions = (
39    d(u_1) + u_2,
40    wedge(d(u_1), u_2),
41    hodge(u_1),
42    ip(beta, u_1),
43    jp(beta, u_1),
44    ld(beta, u_1),
45)
46
47for expression in expressions:
48    print(latex(expression))

Expression trees

This example writes the expression trees of a differential operator and of an expanded derivative in the DOT language of Graphviz. The quickstart guide draws both trees.

 1"""Write the expression trees of symbolic expressions in the DOT language."""
 2
 3from sympy.printing.dot import dotprint
 4
 5from sympde.api import Domain, ScalarFunctionSpace, div, dx, dy, elements_of, grad
 6
 7
 8domain = Domain('Omega', dim=2)
 9space = ScalarFunctionSpace('V', domain)
10u, v = elements_of(space, names='u, v')
11
12
13# [operator-start]
14def L(w):
15    return -div(grad(w)) + w
16
17
18operator_graph = dotprint(L(u))
19# [operator-end]
20
21
22# [derivative-start]
23expr = dx(dy(u * v))
24assert expr == u * dx(dy(v)) + v * dx(dy(u)) + dx(u) * dy(v) + dx(v) * dy(u)
25
26derivative_graph = dotprint(expr)
27# [derivative-end]
28
29
30# Graphs drawn in the documentation, by file name
31graphs = {
32    'operator': operator_graph,
33    'derivative': derivative_graph,
34}

Linearization of nonlinear problems

This example constructs a nonlinear Poisson residual and derives its Jacobian form symbolically.

 1"""Linearize the residual of a nonlinear Poisson problem."""
 2
 3from sympy import exp
 4
 5from sympde.api import BilinearForm, Domain, LinearForm, ScalarFunctionSpace
 6from sympde.api import dot, elements_of, grad, integral, linearize
 7
 8
 9domain = Domain('Omega', dim=2)
10space = ScalarFunctionSpace('V', domain)
11field, increment, test = elements_of(space, names='u, delta_u, v')
12
13
14residual = LinearForm(
15    test,
16    integral(domain, dot(grad(test), grad(field)) - 4.0 * exp(-field) * test),
17)
18
19jacobian = linearize(residual, field, trials=increment)
20assert isinstance(jacobian, BilinearForm)
21assert jacobian(increment, test) == integral(
22    domain,
23    dot(grad(test), grad(increment))
24    + 4.0 * exp(-field) * increment * test
25)