Flexible weak form expressions

[adtzlr]

FElupe does only provide form assemblers for a given formal architecture of weak form expressions (double-dot products), for example a linear form with the gradient of v:

∫ f_ij : v_i,j dx

With some minor additions very flexible weak-form definitions similar to the interior-basis functionality of scikit-fem may be possible with FElupe. This is just a quick and dirty extension, notation and class methods may be improved in the future.

import felupe as fe
import numpy as np

class Basis:
    def __init__(self, field):
        self.field = field
        self.grad = np.einsum(
            "ij,akpe->aijkpe", np.eye(self.field.region.element.dim), self.field.region.dhdX
        )
        self.basis = np.einsum(
            "ij,ap,e->aijpe", 
            np.eye(self.field.region.element.dim), 
            self.field.region.h, 
            np.ones(self.field.region.mesh.ncells)
        )

class LinearForm:
    def __init__(self, v, grad_v=False):
        self.v = v
        self.grad_v = grad_v
        self.dx = v.field.region.dV
        self.form = fe.IntegralForm(None, v.field, self.dx, grad_v=grad_v)

    def integrate(self, weakform, *args, **kwargs):
        if self.grad_v:
            v = self.v.grad
        else:
            v = self.v.basis
        values = np.zeros((len(v), *v.shape[-3:]))
        for a, vbasis in enumerate(v):
            for i, vb in enumerate(vbasis):
                values[a, i] = weakform(vb, *args, **kwargs) * self.dx
        return values.sum(-2)

    def assemble(self, weakform, *args, **kwargs):
        values = self.integrate(weakform, *args, **kwargs)
        return self.form.assemble(values)

class BilinearForm:
    def __init__(self, v, u, grad_v=False, grad_u=False):
        self.v = v
        self.grad_v = grad_v
        self.u = u
        self.grad_u = grad_u
        self.dx = v.field.region.dV
        self.form = fe.IntegralForm(None, v.field, self.dx, u.field, grad_v, grad_u)

    def integrate(self, weakform, *args, **kwargs):

        if self.grad_v:
            v = self.v.grad
        else:
            v = self.v.basis

        if self.grad_u:
            u = self.u.grad
        else:
            u = self.u.basis

        values = np.zeros((len(v), v.shape[-3], len(u), *u.shape[-3:]))

        for a, vbasis in enumerate(v):
            for i, vb in enumerate(vbasis):

                for b, ubasis in enumerate(u):
                    for j, ub in enumerate(ubasis):

                        values[a, i, b, j] = weakform(vb, ub, *args, **kwargs) * self.dx

        return values.sum(-2)

    def assemble(self, weakform, *args, **kwargs):
        values = self.integrate(weakform, *args, **kwargs)
        return self.form.assemble(values)

Using these two classes in the following example is already quite intuitive. Note that no fourth-order tensor has to be created for the assembly of the stiffness matrix.

from felupe.math import ddot

def lform(v, F):
    return ddot(F, v)

def bform(v, u, F):
    return ddot(F, v) * ddot(F, u)

mesh = fe.Cube(n=3)
region = fe.RegionHexahedron(mesh)
field = fe.Field(region, dim=3)
basis = Basis(field)

L = LinearForm(v=basis, grad_v=True)
r = L.integrate(lform, F=field.extract())

a = BilinearForm(v=basis, u=basis, grad_v=True, grad_u=True)
K = a.integrate(bform, F=field.extract())

checks are fine for both linear and bilinear forms:

L.form.fun = field.extract()
check = np.allclose(r, L.form.integrate())
print(check)

a.form.fun = fe.math.dya(field.extract(), field.extract())
check = np.allclose(K, a.form.integrate())
print(check)

[adtzlr]

However, I'm not 100% sure if it is really useful to mimic a functionality of a package which does this close to perfection...

But I think this will improve the usage of FElupe anyway, so this and an additional (Bi)LinearFormMixed, and probably BasisMixed should be included in the future.

[adtzlr]

see progress in branch add-basis-flexible-weakform

[adtzlr]

How to deal with mixed formulations? A proposal for linear forms:

L = fe.LinearFormMixed(v=(basis_displacement, basis_pressure), grad_v=(True, False))

def linearform_mixed(v_, F, J):
    grad_v, q = v_
    return [fe.math.ddot(F, grad_v), fe.math.dot(J, q)]

F = displacement.extract()
L.integrate(linearform_mixed, F=F, J=fe.math.det(F))
L.assemble(linearform_mixed, F=F, J=fe.math.det(F))

And for bilinear forms use the upper triangle entries...?

[adtzlr]

All done, see d738e8b