Adaptive Poisson equation

[Abelsylowlie]

I try to use adaptive method to solve poisson equation with quadratic (or cubic) element.
Comparing with ex22, Delta u_h need to interpolate, and I try to
caculate the Hessian matrix value of the bases function, but it returns None.

I guess those elements are defined by using lbasis (phi and dphi),but the second derivatives of phi are not given.

How to get the term Delta u_h of in interior residual.

[kinnala]

Thanks, it is a good question. There are several options. However, I suggest you employ some finite element knowledge.

Partial derivative of ElementTriP<k> is ElementDG(ElementTriP<k-1>). So you could first find the derivative and then calculate the residual. This is an example:

from skfem import *
mesh = MeshTri().refined(3)
basis = Basis(mesh, ElementTriP3())
grad_basis = basis.with_element(ElementVector(ElementDG(ElementTriP2())))
u = basis.project(lambda x: x[0] ** 3)
grad_u = grad_basis.project(basis.interpolate(u).grad)

@Functional
def residual(w):
    (uxx, uxy), (uyx, uyy) = w['grad_u'].grad
    return (uxx + uyy) ** 2

out = residual.elemental(grad_basis, grad_u=grad_u)

from skfem.visuals.matplotlib import *
plot(mesh, out, colorbar=True)
show()

I hope it works out.

[kinnala]

You probably don't want to redefine u inside eval_estimator like you've done.

[kinnala]

I can try it out later, I'm now on mobile so cannot test properly.

[kinnala]

I had to fix one line in load_func and comment out another line:

from skfem import *
from skfem.models.poisson import laplace
from skfem.helpers import grad
import numpy as np

m = MeshTri.init_lshaped().refined()
e = ElementTriP2()

def load_func(x, y):
    return 1. + 0. * x  # fixed

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

def eval_estimator(m, u):    
    # interior residual
    basis = Basis(m, e)

    grad_basis = basis.with_element(ElementVector(ElementDG(ElementTriP1()))) 
    #u = basis.project(lambda x: x[0] ** 2 + x[1] ** 2) 
    w = {'grad_u': grad_basis.project(basis.interpolate(u).grad) }
 
    @Functional 
    def interior_residual(w): 
        h = w.h
        x, y = w.x
        (uxx, uxy), (uyx, uyy) = w['grad_u'].grad 
        return h ** 2 * (uxx + uyy + load_func(x, y)) ** 2 

    eta_K = interior_residual.elemental(grad_basis, **w) 

    # facet jump
    fbasis = [InteriorFacetBasis(m, e, side=i) for i in [0, 1]]
    w = {'u' + str(i + 1): fbasis[i].interpolate(u) for i in [0, 1]}
    
    @Functional
    def edge_jump(w):
        h = w.h
        n = w.n
        dw1 = grad(w['u1'])
        dw2 = grad(w['u2'])
        return h * ((dw1[0] - dw2[0]) * n[0] +
                    (dw1[1] - dw2[1]) * n[1]) ** 2

    eta_E = edge_jump.elemental(fbasis[0], **w)
    
    tmp = np.zeros(m.facets.shape[1])
    np.add.at(tmp, fbasis[0].find, eta_E)
    eta_E = np.sum(.5 * tmp[m.t2f], axis=0)
    
    return eta_K + eta_E

if __name__ == "__main__":
    from skfem.visuals.matplotlib import draw, plot, show
    draw(m)

for itr in reversed(range(6)):
        
    basis = Basis(m, e)
    
    K = asm(laplace, basis)
    f = asm(load, basis)

    u = solve(*condense(K, f, D=basis.find_dofs()))
        
    if itr > 0:
        m = m.refined(adaptive_theta(eval_estimator(m, u)))
        draw(m)
        show()

[kinnala]

[Abelsylowlie]

Happy new year to you！ And thanks so much for this library and your help!