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.