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

Errors
1

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.

2

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?

3

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.

4

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.

5

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

6

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.

7

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