Advanced expression manipulation

Tensor-product representation

TensorExpr lowers an integrated form and expresses it through one- dimensional symbolic factors. This is useful to inspect the Kronecker structure consumed by tensor-product discretizations.

"""Expose the tensor-product structure of a scalar bilinear form."""

from sympde.api import BilinearForm, Constant, Domain, ScalarFunctionSpace
from sympde.api import TensorExpr, dot, elements_of, grad, integral


domain = Domain('Omega', dim=2)
space = ScalarFunctionSpace('V', domain)
trial, test = elements_of(space, names='u, v')
coefficient = Constant('mu', real=True)

form = BilinearForm(
    (trial, test),
    integral(
        domain,
        coefficient * trial * test + dot(grad(trial), grad(test)),
    ),
)

tensor_expression = TensorExpr(form, domain=domain)
print(tensor_expression)

The example is a complete script, which the documentation workflow runs before Sphinx builds the page.

Logical expressions

LogicalExpr pulls expressions on a mapped physical domain back to its logical domain. The analytical-mapping example demonstrates symbolic and numerical use of this transformation:

"""Work with built-in and user-defined analytical mappings."""

import numpy as np
from sympy import pi, sin, symbols

from sympde.api import CollelaMapping2D, Constant, Domain, LogicalExpr, Mapping


# Inspect a built-in mapping symbolically.
mapping = CollelaMapping2D('M', dim=2)
mapped_domain = mapping(Domain('Omega', dim=2))
x1, x2 = symbols('x1, x2', real=True)
eps, k1, k2 = [Constant(name) for name in ('eps', 'k1', 'k2')]

expected_x = 2.0 * (
    x1 + eps * sin(2.0 * pi * k1 * x1) * sin(2.0 * pi * k2 * x2)
) - 1.0
assert LogicalExpr(mapping[0], mapped_domain) == expected_x


# Fix its parameters and evaluate it numerically.
numerical_mapping = CollelaMapping2D(
    'M_numeric',
    dim=2,
    eps=0.1,
    k1=1.0,
    k2=1.0,
).get_callable_mapping()
np.testing.assert_allclose(numerical_mapping(0.25, 0.5), (-0.5, 0.0))


# Define a custom mapping by specifying its coordinate expressions.
class ComplexSquareMapping(Mapping):
    """
    Map (x1, x2) to (x, y) such that x + i*y = A * (x1 + i*x2)**2.

    Expanding the square gives x = A * (x1**2 - x2**2) and y = 2 * A * x1 * x2.
    """

    _expressions = {
        'x': 'A * (x1**2 - x2**2)',
        'y': '2 * A * x1 * x2',
    }


custom_mapping = ComplexSquareMapping('F', dim=2, A=1.0).get_callable_mapping()
np.testing.assert_allclose(custom_mapping(0.5, 0.25), (0.1875, 0.25))

Nonlinear forms

linearize differentiates a nonlinear residual with respect to a field and returns its symbolic Jacobian form:

"""Linearize the residual of a nonlinear Poisson problem."""

from sympy import exp

from sympde.api import BilinearForm, Domain, LinearForm, ScalarFunctionSpace
from sympde.api import dot, elements_of, grad, integral, linearize


domain = Domain('Omega', dim=2)
space = ScalarFunctionSpace('V', domain)
field, increment, test = elements_of(space, names='u, delta_u, v')


residual = LinearForm(
    test,
    integral(domain, dot(grad(test), grad(field)) - 4.0 * exp(-field) * test),
)

jacobian = linearize(residual, field, trials=increment)
assert isinstance(jacobian, BilinearForm)
assert jacobian(increment, test) == integral(
    domain,
    dot(grad(test), grad(increment))
    + 4.0 * exp(-field) * increment * test
)