MWE_nonlinear_DG

import ngsolve as ng
from netgen.geom2d import unit_square
from ngsolve.webgui import Draw

def SimpleNewtonSolve(gfu,a,tol=1e-13,maxits=25):
    res = gfu.vec.CreateVector()
    du = gfu.vec.CreateVector()
    fes = gfu.space
    for it in range(maxits):
        print ("Iteration {:3}  ".format(it),end="")
        a.Apply(gfu.vec, res)
        a.AssembleLinearization(gfu.vec)
        du.data = a.mat.Inverse(fes.FreeDofs()) * res
        gfu.vec.data -= du

        #stopping criteria
        stopcritval = ng.sqrt(abs(ng.InnerProduct(du,res)))
        print ("<A u",it,", A u",it,">_{-1}^0.5 = ", stopcritval)
        if stopcritval < tol:
            break

T_ex = ng.x*(ng.x-1)*ng.y*(ng.y-1)
f = -2*(ng.y*(ng.y-1) + ng.x*(ng.x-1))

mesh = ng.Mesh(unit_square.GenerateMesh(maxh=0.05))

order=1
fes = ng.L2(mesh, order=order, dgjumps=True)
u,v = fes.TnT()
gfu = ng.GridFunction(fes, name="uDG")

jump_u = u-u.Other()
jump_v = v-v.Other()
n = ng.specialcf.normal(2)
mean_dudn = 0.5*n * (ng.grad(u)+ng.grad(u.Other()))
mean_dvdn = 0.5*n * (ng.grad(v)+ng.grad(v.Other()))

alpha = 4
h = ng.specialcf.mesh_size
a = ng.BilinearForm(fes)
a += ng.SymbolicBFI(ng.grad(u)*ng.grad(v))
a += ng.SymbolicBFI(alpha*order**2/h*jump_u*jump_v, skeleton=True)
a += ng.SymbolicBFI(alpha*order**2/h*u*v, ng.BND, skeleton=True)
a += ng.SymbolicBFI(-mean_dudn*jump_v -mean_dvdn*jump_u, skeleton=True)
a += ng.SymbolicBFI(-n*ng.grad(u)*v-n*ng.grad(v)*u, ng.BND, skeleton=True)

#==========================
# Version 1. This works.
#a.Assemble()

#L = ng.LinearForm(fes)
#L += ng.SymbolicLFI(f*v)
#L.Assemble()

#gfu.vec.data = a.mat.Inverse() * L.vec
# =========================

#==========================
# Version 2. This does NOT work, and I don't know why.
a += ng.SymbolicBFI(-f*v)
SimpleNewtonSolve(gfu, a)
#==========================

Draw(T_ex, mesh)
Draw(gfu, mesh)
Draw(gfu - T_ex, mesh)