Mixed poisson equation with Nitsche method to apply boundary conditions

In the mixed formulation of the poisson equation, the fluxes boundary conditions are strong boundary conditions and are applied using DirichletBC as in https://fenicsproject.org/docs/dolfin/latest/python/demos/mixed-poisson/demo_mixed-poisson.py.html.

On the other hand, the Nitsche method allows you to apply strong boundary conditions in a weak fashion through the variational form (for instance https://github.com/MiroK/fenics-nitsche/blob/master/poisson/poisson_circle_dirichlet.py)

I am looking for a way to combine both so I can apply the flux boundary conditions in the mixed formulation using the Nitsche method. Any help is appreciated.

I haven’t worked much with the mixed Poisson problem, so there’s probably something smarter in the porous media literature, but my quick attempt of

#bc = DirichletBC(W.sub(0), G, boundary)

# Mark boundary sub-domains:
class Gamma_N(SubDomain):
    def inside(self,x,on_boundary):
        return boundary(x)
boundary_marker = MeshFunction("size_t",mesh,mesh.topology().dim()-1,0)
Gamma_N().mark(boundary_marker,1)
ds = Measure("ds",domain=mesh,subdomain_data=boundary_marker)

# Formulation:
n = FacetNormal(mesh)
h = CellDiameter(mesh)
C_pen = Constant(1e1)
res = (dot(sigma, tau) + div(tau)*u + div(sigma)*v + f*v)*dx \
      - dot(tau,n)*u0*ds(0) - dot(tau,n)*u*ds(1) - dot(sigma-G,n)*v*ds(1) \
      + C_pen*h*inner(dot(sigma-G,n),dot(tau,n))*ds(1)
a = lhs(res)
L = rhs(res)

seems to work reasonably well, with \Vert u^h-u\Vert_{L^2} converging at first order, as expected for DG0 spaces.

My “derivation” of this is to do the following:

  • Keep the boundary term from div(tau)*u*dx on \Gamma_N, since \tau\cdot n is no longer zero.
  • Integrate div(sigma)*v*dx by parts, replace dot(sigma,n) with dot(G,n) in the boundary term, then integrate by parts again, but without replacing dot(sigma,n).
  • Add a penalty term, where the factor of h compensates for the difference in units between ds and dx (comparing units with the dot(sigma,tau)*dx term).
1 Like

Thank you. Not sure if you have a typo when you say: div(tau)*u*dx on \Gamma_N? Did you mean dot(tau, n)*u*ds(1)?

Not a typo, but now that I look at it again, I see that the wording is unclear. What I meant was: “Keep the boundary term -dot(tau,n)*u*ds(1) on \Gamma_N, which comes from the integration by parts used to obtain the term div(tau)*u*dx.”