Add experimental.lab

[kinnala]

Work on experimental and optional symbolic "language" for defining weak forms. This requires autograd and is continuation of my past work on automatic differentiation.

Some examples:

"""Navier-Stokes lid driven cavity."""
from skfem.experimental.lab import *
import numpy as np

m = MeshTri().refined(4)

elem = ElementVector(ElementTriP2()) * ElementTriP1()
basis = Basis(m, elem)

u, p, v, q = symbols(elem)

form = (ddot(grad(u), grad(v))
        + dot(mul(grad(u), u), v)
        - p * div(v)
        - q * div(u)
        - 1e-3 * p * q)

x = basis.zeros()
for itr in range(20):
    x[basis.get_dofs('top').all('u^1^1')] = 100.
    dt = solve(*condense(*form.assemble(basis, x=x),
                         D=basis.get_dofs().all(['u^1^1', 'u^2^1'])))
    x += dt
    print(np.linalg.norm(dt))

(u, ubasis), (p, pbasis) = basis.split(x)

pbasis.plot(p).show()
ubasis.plot(u).show()

"""Nonlinear elasticity."""
from skfem.experimental.lab import *
import numpy as np

m = MeshQuad.init_tensor(np.linspace(0, 5, 20),
                         np.linspace(0, 0.5, 5))

elem = ElementVector(ElementQuad1())
basis = Basis(m, elem)

u, _ = symbols(elem)

epsu = .5 * (grad(u) + grad(u).T + grad(u).T @ grad(u))
sigu = epsu + eye(tr(epsu), 2)
form = ddot(sigu, epsu) - 1e-5 * u[1]

x = basis.zeros()
for itr in range(20):
    x[basis.get_dofs('right').all('u^1')] = -np.minimum(itr, 10) / 20
    x[basis.get_dofs('right').all('u^2')] = 0
    dt = solve(*condense(*form.assemble(basis, x=x),
                         D=basis.get_dofs({'left', 'right'}).all()))
    x += dt
    print(np.linalg.norm(dt))

mdefo = m.translated(x[basis.nodal_dofs])
ax = m.draw(color='r')
mdefo.draw(ax=ax).show()