Adaptive Poisson problem

[Abelsylowlie]

Hi, @kinnala, you gave me an example of adaptive Poisson problem in #844. However, when #829 fixed, I found something wrong when run the following codes.

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()))) 
    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()

This refinements are as follows

[kinnala]

I'll look into it. I don't know how #966 can affect this code because Basis.split is not used.

[kinnala]

There is something wrong with basis.find_dofs() which has been also deprecated for about one year. Change it to basis.get_dofs() and it should work.

[kinnala]

Ah yes, I added some default boundaries to all domains and this seems to have caused basis.find_dofs() to work differently.

[Abelsylowlie]

Thanks!