Solve Poisson with discontinuous Galerkin (DG=0)

Question

Hi,
I try to solve Poisson equation with Dirichlet boundary conditions and discontinuous Galerkin element (DG=0)
The BC are:x[0]=0 c=500 and x[0]=1 c=100.
There is no source term.

I have used the demo_dg-poisson.py to start
https://github.com/FEniCS/dolfin/blob/master/demo/undocumented/dg-poisson/python/demo_dg-poisson.py

Here is my code (FEniCS 2016.1)

When i plot the solution field at the end, the concentrations show a non-linear gradient i do not understand.
I guess the BC or the variational formulation are ill-posed.
If someone can help.

Thanks.

from dolfin import *

tolerance_boundary_location = DOLFIN_EPS*1000;
class DirichletBoundary(SubDomain):
  def inside(self, x, on_boundary):
    return ( abs(x[0])  < tolerance_boundary_location or abs(x[0]-1)  < tolerance_boundary_location)

# Create mesh and define function space
mesh = UnitCubeMesh(10, 10, 10)
V = FunctionSpace(mesh, 'DG', 0)

# Define test and trial functions
u = TrialFunction(V)
v = TestFunction(V)

# Define normal vector and mesh size
n = FacetNormal(mesh)
h = CellSize(mesh)
h_avg = (h('+') + h('-'))/2

# Define the source term f, Dirichlet term u0 and Neumann term g
u0 = Expression('-400*x[0] + 500', degree=2)
f = Constant(0)

# Mark facets of the mesh
boundaries = FacetFunction('size_t', mesh, 0)
DirichletBoundary().mark(boundaries, 1)

# Define outer surface measure aware of Dirichlet and Neumann boundaries
ds = Measure('ds', domain=mesh, subdomain_data=boundaries)

# Define parameters
alpha = 4.0
gamma = 8.0

plot(boundaries, title='Dirichlet BC are applied here', interactive=True)

# Define variational problem
a = dot(grad(v), grad(u))*dx \
   - dot(avg(grad(v)), jump(u, n))*dS \
   - dot(jump(v, n), avg(grad(u)))*dS \
   + alpha/h_avg*dot(jump(v, n), jump(u, n))*dS \
   - dot(grad(v), u*n)*ds(1) \
   - dot(v*n, grad(u))*ds(1) \
   + (gamma/h)*v*u*ds(1)
L = v*f*dx - u0*dot(grad(v), n)*ds(1) + (gamma/h)*u0*v*ds(1)

# Compute solution
u = Function(V)
solve(a == L, u)

# Plot solution
plot(u, interactive=True)

Comments

What do you exactly mean with the 'non-linear gradient' in the results for DG=0? I have run your results and except from the stepwise behavior of the solution (which you can't expect to be better when using DG=0 ), the solution looks fine.

---

You are right. I have re-tested the code after reading your message and the solution is linear. I must have modified something when i have edited this post that has corrected my initial error.
Well, glad it works!

---

Does it make sense to solve with DG0? The discretization of the elliptic operator will be zero:
dot(grad(v), grad(u))*dx = 0 The finite element space is not in H_1.

---

Hmm, motivation is (i) how to have a working code with DG and (ii) i am curious about the solution with a field constant element by element.