The error in using the DG method to solve an elliptic problem

Hi, everyone.

I am using the DG method to solve a quasilinear elliptic equation, but I have encountered difficulties in implementation. The following is my code, and the error in the execution results is likely caused by the definition of the variable Nonlinear2.

How should I modify my implementation?

from fenics import *
from mshr import *
import numpy as np
import sympy as sp
import math

x, y = sp.symbols('x[0] x[1] ')

p = 4.0
pdeg = 1
size = 10

def vecnorm(vec):
  return pow(inner(vec,vec),0.5)

mesh = RectangleMesh(Point(1, 1), Point(2, 2), size, size)
r = Expression("pow((pow(x[0],2)+pow(x[1],2)),0.5)", degree=3)
uex = Expression("pow(r,(p-2)/(p-1))",degree=3,r=r,p=p)
f = Constant(0.0)
f = Expression(sp.printing.ccode(f), degree=3)


W = VectorFunctionSpace(mesh , "DG", pdeg )
V = FunctionSpace(mesh , "DG", pdeg )
u = TrialFunction(V)
v = TestFunction(V)

n = FacetNormal(mesh)
h = CellDiameter(mesh)
h_avg = (h('+') + h('-'))/2


# Boundary
uD = uex
bc = DirichletBC(V, uD, "on_boundary")
uD = project(uD, V)

# Parameters
delta = Constant(1e-14)
alpha = Constant(10*pdeg)
sigma_b = alpha/h
sigma_o = alpha/h_avg

# Initialization 
u0 = project(uex, V)

eps = 1.0
tol = 1.0e-14
iter = 0
maxiter = 10
while eps > tol and iter < maxiter:

  u = TrialFunction(V)
  v = TestFunction(V)
  gradu0 = Function(W)
  gradu0 = project(grad(u0),W)

  Nonlinear1 = pow((delta + vecnorm(gradu0)),p-2)
  Nonlinear2 = pow((delta + vecnorm(1/h_avg*jump(u0,n))),p-2)
  Nonlinear3 = pow((delta + vecnorm(1/h*(u0-uD))),p-2)

  b0 = -inner(avg(Nonlinear1*grad(u)),jump(v,n))*dS - inner(avg(Nonlinear2*grad(v)),jump(u,n))*dS \
        + sigma_o*inner(jump(u,n),jump(v,n))*dS 
  b1 = -inner(Nonlinear1*grad(u),v*n)*ds - inner(Nonlinear3*grad(v),(u-uD)*n)*ds \
        + sigma_b*inner(v*n,(u-uD)*n)*ds 

  N = inner(Nonlinear1*grad(u),grad(v))*dx

  L = f*v*dx


  Lhs = N + b0 + b1
  Rhs = L

  form = Lhs-Rhs
  Lhs = lhs(form)
  Rhs = rhs(form)


  u = Function(V)
  solve(Lhs == Rhs, u, bc)


  print ("iter=%d:" %iter)
  u_er = project(u-u0, V)
  u_L2iter = sqrt(assemble(inner(u_er,u_er)*dx))

  eps_u = u_L2iter
  eps = eps_u

  print ("iter=%d: norm=%g " %(iter, eps))

  u0 = u
  iter = iter + 1

UFLException: Cannot restrict an expression twice.

As stated by the error message, this is because you use multiple instances of jump and avg.

You can take the avg on top of this in:

What are you actually trying to implement here? Could you show the mathematics?

Hi, dokken

Thank you for your response.

This is the paper I referred to.
https://scispace.com/pdf/discontinuous-galerkin-finite-element-approximation-of-225ik17iy1.pdf
Equation (2.2) is the DG scheme I referred to.

In this scheme, there is not a \langle \left[\cdot \right]\rangle which is what you have when you have jump(u0, n) inside an avg. Please carefully write out your formulation.

It looks like this is the term Chen is referring to:
image

I didn’t spot this, sorry about that. But what does the average of a jump mean?
avg(jump(a))= avg(a("+") - a("-"))= 0.5 * ( (a("+")-a("-"))("+") + (a("+")-a("-"))("-") ) where some pf these restrictions do not really make sense.

Yes, I agree it would make a lot more sense to extract the single valued \mu expression from the average operator. Then Dolfinx can handle it properly:

\int_\Gamma \mu(h^{-1}|[\![w]\!]|)\langle\nabla v\cdot n\rangle [\![w]\!] ds

So, in terms of implementation for @Chen-03:

inner(Nonlinear2*avg(grad(v)),jump(u,n))*dS