Evaluating local functionals

[kinnala]

I did again some stuff related to adaptive finite element methods.

I ended up writing stuff like this for evaluating local error estimators:

    # interior estimator 
    basis = InteriorBasis(m, e)
    hK = basis.mesh_parameters()
    xK, yK = basis.global_coordinates()

    eta_K = np.sum(hK**2 / k(xK, yK) * load_func(xK, yK)**2 * basis.dx, axis=1)

and

    # facet estimator
    fbasis = [FacetBasis(m, e, side=i) for i in [0, 1]]
    u1, du1 = fbasis[0].interpolate(u, derivative=True)
    u2, du2 = fbasis[1].interpolate(u, derivative=True)
    hE = fbasis[0].mesh_parameters()
    xE, yE = fbasis[0].global_coordinates()
    eps = 1e-6
    n = fbasis[0].normals
    xE1 = xE - eps*n[0]
    yE1 = yE - eps*n[1]
    xE2 = xE + eps*n[0]
    yE2 = yE + eps*n[1]

    eta_E = np.sum(hE / k(xE, yE) *\
                   ((k(xE1, yE1)*du1[0] - k(xE2, yE2)*du2[0])*n[0] +\
                    (k(xE1, yE1)*du1[1] - k(xE2, yE2)*du2[1])*n[1])**2 * fbasis[0].dx, axis=1)

    # add eta_E to elements
    tmp = np.zeros(m.facets.shape[1])
    np.add.at(tmp, fbasis[0].find, eta_E)
    eta_E = np.sum(0.5*tmp[m.t2f], axis=0)

This looks pretty similar to what one does in skfem.assembly.asm. Maybe we could add skfem.assembly.evaluate (or similar) for evaluating functionals.

[kinnala]

Did some quasi-cool stuff yesterday. Finally added adaptive refinement of triangular meshes.

from skfem import *
from skfem.models.poisson import laplace
import numpy as np
%matplotlib inline

m = MeshTri.init_symmetric()
m.refine(2)
e = ElementTriP1()

def load_func(x, y):
    return (x > 0.5)*(y > 0.5)

@linear_form
def load(v, dv, w):
    x, y = w.x
    return load_func(x, y) * v

def eval_estimator(m, u):    
    # interior residual
    basis = InteriorBasis(m, e)
    hK = basis.mesh_parameters()
    x, y = basis.global_coordinates()

    eta_K = np.sum(hK**2 * load_func(x, y)**2 * basis.dx, axis=1)

    # facet jump
    fbasis = [FacetBasis(m, e, side=i) for i in [0, 1]]

    u1, du1 = fbasis[0].interpolate(u, derivative=True)
    u2, du2 = fbasis[1].interpolate(u, derivative=True)
    hE = fbasis[0].mesh_parameters()
    n = fbasis[0].normals

    eta_E = np.sum(hE * ((du1[0] - du2[0])*n[0] +\
                         (du1[1] - du2[1])*n[1])**2 * fbasis[0].dx, axis=1)

    tmp = np.zeros(m.facets.shape[1])
    np.add.at(tmp, fbasis[0].find, eta_E)
    eta_E = np.sum(0.5*tmp[m.t2f], axis=0)

    return eta_K + eta_E

m.draw()

for itr in range(13):

    if itr > 1:
        m.refine(adaptive_theta(eval_estimator(m, u)))

    basis = InteriorBasis(m, e)

    K = asm(laplace, basis)
    f = asm(load, basis)
    u = np.zeros_like(f)

    I = m.interior_nodes()
    u[I] = solve(*condense(K, f, I=I))

m.draw()
m.plot(u, smooth=True)

This solves Poisson equation adaptively.

Note that the function eval_estimator relates to this issue.

[gdmcbain]

Very nice! I will be using this a lot!

[kinnala]

I'm now playing around with the syntax related to this issue. This is what I have now:

@integral
def interior_residual(w):
    h, x, y = w.h, w.x
    return h**2 * load_func(x, y)**2

@integral(facets=True, interior=True)
def facet_residual(w):
    dw1, dw2, n = w.dw1, w.dw2, w.n
    return h * ((dw1[0] - dw2[0])*n[0] +\
                (dw1[1] - dw2[1])*n[1])**2

integrate(interior_residual, basis=basis)

[kinnala]

There is an example demonstrating this feature at

https://github.com/kinnala/scikit-fem/blob/master/docs/examples/ex22.py

[kinnala]

I guess this can be closed now.