Surface integral for the mixed poisson problem

1

Hi all,

I am try to duplicate and validate the mixed poisson demo. However, when I get the solution and compare the gradient of solution given by grad(u) and sigma, I find the surface integral always gives me zero. Is there anyone can tell me the reason. I will appreciate it very much.

my code is

from dolfin import *
import matplotlib.pyplot as plt
mesh = UnitSquareMesh(32, 32)
BDM = FiniteElement("BDM", mesh.ufl_cell(), 1)
DG  = FiniteElement("DG", mesh.ufl_cell(), 0)
W = FunctionSpace(mesh, BDM * DG)
(sigma, u) = TrialFunctions(W)
(tau, v) = TestFunctions(W)
tol = 1e-14
class BoundaryX0(SubDomain):
      def inside(self, x, on_boundary):
          return on_boundary and near(x[0], 0, tol)

class BoundaryX1(SubDomain):
      def inside(self, x, on_boundary):
          return on_boundary and near(x[0], 1, tol)

class BoundaryY0(SubDomain):
      def inside(self, x, on_boundary):
          return on_boundary and near(x[1], 0, tol)

class BoundaryY1(SubDomain):
      def inside(self, x, on_boundary):
          return on_boundary and near(x[1], 1, tol)
boundary_markers = MeshFunction('size_t', mesh,1)
boundary_markers.set_all(9999)
bx0 = BoundaryX0()
bx1 = BoundaryX1()
by0 = BoundaryY0()
by1 = BoundaryY1()
bx0.mark(boundary_markers, 0)
bx1.mark(boundary_markers, 1)
by0.mark(boundary_markers, 2)
by1.mark(boundary_markers, 3)
# Redefine boundary integration measure
ds = Measure('ds', domain=mesh, subdomain_data=boundary_markers)

f = Expression("10*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)", degree=2)
a = (dot(sigma, tau) + div(tau)*u + div(sigma)*v)*dx
L = - f*v*dx
class BoundarySource(UserExpression):
    def __init__(self, mesh, **kwargs):
        self.mesh = mesh
        super().__init__(**kwargs)
    def eval_cell(self, values, x, ufc_cell):
        cell = Cell(self.mesh, ufc_cell.index)
        n = cell.normal(ufc_cell.local_facet)
        g = sin(5*x[0])
        values[0] = g*n[0]
        values[1] = g*n[1]
    def value_shape(self):
        return (2,)

G = BoundarySource(mesh, degree=2)
def boundary(x):
    return x[1] < DOLFIN_EPS or x[1] > 1.0 - DOLFIN_EPS
bc = DirichletBC(W.sub(0), G, boundary)
w = Function(W)
solve(a == L, w, bc)
(sigma, u) = w.split()
nn = FacetNormal(mesh)
print(assemble(dot(grad(u),nn)*ds))

and my output is
0.0

2

Your u is in a discontinuous Lagrange function of order 0. That means that it is piecewise constant on each element, and therefore grad(u) will always be 0, as the derivative of a constant is 0.

3

Thanks for your valuable response. You are right that if I use the 0 order element the derivative will be zero.

Have a nice day.