ElementDG() for ElementQuadP()

[Abelsylowlie]

Hi, all. I try to use example 07 on quadrilateral mesh. It is OK for ElementQuad1() and ElementQuad2(). But for ElementQuadP(), it gets this

I would like to know what causes this error and how to fix it. Thanks in advance.

[kinnala]

Thanks for the report. I think there is a bug in the DOF indexing. I should be able to fix it in few days.

[Abelsylowlie]

Thank you! I want to know the difference between ElementQuadP() and spectral element with nodal basis, since both of them use Legendre polynomial.

[kinnala]

Can you report which order you used for ElementQuadP?

[kinnala]

For me this works just fine:

"""Discontinuous Galerkin method."""

from skfem import *
from skfem.helpers import grad, dot, jump
from skfem.models.poisson import laplace, unit_load

m = MeshQuad().refined(2)
e = ElementDG(ElementQuadP(7))
alpha = 1e-3

ib = Basis(m, e)
bb = FacetBasis(m, e)
fb = [InteriorFacetBasis(m, e, side=i) for i in [0, 1]]

@BilinearForm
def dgform(u, v, p):
    ju, jv = jump(p, u, v)
    h = p.h
    n = p.n
    return ju * jv / (alpha * h) - dot(grad(u), n) * jv - dot(grad(v), n) * ju

@BilinearForm
def nitscheform(u, v, p):
    h = p.h
    n = p.n
    return u * v / (alpha * h) - dot(grad(u), n) * v - dot(grad(v), n) * u

A = asm(laplace, ib)
B = asm(dgform, fb, fb)
C = asm(nitscheform, bb)
b = asm(unit_load, ib)

x = solve(A + B + C, b)

M, X = ib.refinterp(x, 4)

def visualize():
    from skfem.visuals.matplotlib import plot, draw
    ax = draw(M, boundaries_only=True)
    return plot(M, X, shading="gouraud", ax=ax, colorbar=True)

if __name__ == "__main__":
    visualize().show()

[Abelsylowlie]

It did not work using the above snippet for order >= 3 on my laptop with Windows 10 and scikit-fem-7.0.1.

[Abelsylowlie]

It works on the current unreleased version. Thanks.