Mixed poisson equation with Nitsche method to apply boundary conditions

1

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.

2

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
3

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)?

4

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.”