Specify boundary in dual-mixed formulation

Question

Hi friends,

I can't seem to find out how to specify boundaries correctly for the dual mixed Poisson problem.
I want to give left and right a fixed potential (1 and 0) and on the lower and upper boundary zero outward flux.
I used this example as guide here

The Dirichlet Boundary Conditions are here natural (due to the formulation). I tried to put g = 1-x into the formulation. Is it correct that when I specify upper and lower edge as DirichletBC (although physically Neumann) , then left and right remain to be evaluated as natural boundary domains? Hence I thought my function g would only be evaluated at left and right bound to 1 and 0 respectively. Hence the right side $g\cdot vm\cdot ds$. Is that correct?

Then I tried to make the outward flux at upper and lower edge 0, which I tried to express like in the example with the overloaded eval_cell function.

Yet neither seems to work, I cannot get the $\sigma$ on the upper and lower bound to be tangential to the bound (always normal, why?), nor does the left boundary seem to be 1.
What am I doing wrong?

As solution we expect a linearly decreasing potential from 1 to 0.

from dolfin import *
nn = 150
mesh = UnitSquareMesh(nn,nn)

HDIV = FunctionSpace(mesh,"RT",1)
DG = FunctionSpace(mesh,"DG",0)
W = HDIV * DG

(sigma,um) = TrialFunctions(W)
(tau,vm) = TestFunctions(W)

J = Expression("0.0")

def neumann(x):
    return near(x[1], 0.0) or near(x[1],1.0)

class BoundarySource(Expression):
    def __init__(self, mesh):
        self.mesh = mesh
    def eval_cell(self, values, x, ufc_cell):
        cell = Cell(self.mesh, ufc_cell.index)
        n = cell.normal(ufc_cell.local_facet)
        g = 0
        values[0] = g*n[0]
        values[1] = g*n[1]
    def value_shape(self):
        return (2,)

G = BoundarySource(mesh)
bc_neumann = DirichletBC(W.sub(0), G, neumann)

g= Expression("1.0-x[0]")

am = (dot(sigma,tau) + div(tau)*um + div(sigma)*vm)*dx
fm = -g*vm*dx   
w = Function(W)
solve(am == fm, w,bc_neumann)
(sigma,um) = w.split()

plot(um)
interactive()

from dolfin import *

nn = 150
mesh = UnitSquareMesh(nn,nn)

HDIV = FunctionSpace(mesh,"RT",1)
DG = FunctionSpace(mesh,"DG",0)
W = HDIV * DG

(sigma,um) = TrialFunctions(W)
(tau,vm) = TestFunctions(W)

J = Expression("0.0")

def neumann(x):
    return near(x[1], 0.0) or near(x[1],1.0)

class BoundarySource(Expression):
    def __init__(self, mesh):
        self.mesh = mesh
    def eval_cell(self, values, x, ufc_cell):
        cell = Cell(self.mesh, ufc_cell.index)
        n = cell.normal(ufc_cell.local_facet)
        g = 0
        values[0] = g*n[0]
        values[1] = g*n[1]
    def value_shape(self):
        return (2,)

G = BoundarySource(mesh)
bc_neumann = DirichletBC(W.sub(0), G, neumann)

g= Expression("1.0-x[0]")

am = (dot(sigma,tau) + div(tau)*um + div(sigma)*vm)*dx
fm = -g*vm*ds   
w = Function(W)
solve(am == fm, w,bc_neumann)
(sigma,um) = w.split()

plot(um)
interactive()

Answer 1

The right hand side is incorrect. You are correct in that the Dirichlet boundary condition becomes a natural boundary boundary condition, but it should be integrated against a testfunction in the $H(\operatorname{div})$ space. Also, the integral should only cover the Dirichlet part of the boundary.

If the Dirichlet boundary condition reads

$$ u = g\, \mbox{ on }\Gamma_D, $$

then the right hand side should have the form

$$ L = \int_{\Gamma_D} (q\cdot n)g \, ds $$

where $q\in H(\operatorname{div})$. See also example below.

from dolfin import *

N = 32
mesh = UnitSquareMesh(N, N)

G_N = CompiledSubDomain("(x[1]<DOLFIN_EPS)||(x[1]>1-DOLFIN_EPS)")
G_D = CompiledSubDomain("(x[0]<DOLFIN_EPS)||(x[0]>1-DOLFIN_EPS)")

facet_domains = FacetFunction("size_t", mesh)
G_N.mark(facet_domains, 1)
G_D.mark(facet_domains, 2)

ds = Measure("ds", domain = mesh, subdomain_data = facet_domains)

RT = FiniteElement("RT", triangle, 1)
DG = FiniteElement("DG", triangle, 0)
P1 = FiniteElement("CG", triangle, 1)
W  = FunctionSpace(mesh, RT * DG)

n = FacetNormal(mesh)
U = Expression("1-x[0]", P1)
f = Constant(0)
g = Expression("1-x[0]", P1)
h = Constant((0., 0.))

(p, u) = TrialFunctions(W)
(q, v) = TestFunctions(W)

a = (inner(p,q) + u * div(q) + v * div(p)) * dx
L = f * v * dx + g * inner(q, n) * ds(2)

bc = DirichletBC(W.sub(0), h, facet_domains, 1)

w = Function(W); p, u = w.split()
solve(a == L, w, bc)

print errornorm(U, u)