Using DOLFIN (in Python, preferably), I'd like to implement Nitsche's method based on A. Hansbo, P. Hansbo, An unfitted finite element method, based on Nitsche's method, for elliptic interface problems, Computer Methods in Applied Mechanics and Engineering, Volume 191, Issues 47-48, , p. 5537-5552, 2002.

I have two interfaces corresponding to three meshes. For the purposes of description, let's suppose they are:

- $[0, 1/3] \times [0, 1]$,
- $[1/3, 2/3] \times [0, 1]$,
- $[2/3, 1] \times [0, 1]$

Mesh nodes should match up on interfaces. Put another way, I should be able to join the three meshes to get one mesh without any hanging nodes, cut cells, or deformities. I want to mesh these three regions so that on the interfaces between regions, there is an extra set of nodes (sort of like dG) that I can use to enforce a jump in the flux (or the field itself). Is there a way to set up a mesh in DOLFIN in this fashion?

Within each region, I plan on using cG with Lagrange elements of order 1. An equivalent formulation should be enriching the continuous finite element space with functions that admit jumps in the field and flux across interfaces. If I can't rig the meshing to implement Nitsche's method, is it possible to use DOLFIN to enrich the finite element space instead?

If it matters, I'm solving a nonlinear Poisson equation in 2D.

From what I've gathered, there are kludgey ways to do this in DOLFIN that are slow (see Christian Waluga's question on mortar methods on the old FEniCS Launchpad site; my case also involves no overlap, but the meshes match at interfaces), there are related methods for weakly enforcing Dirichlet boundary conditions (see Andre Massing's answer to this question on the old FEniCS Launchpad site), and there is DOLFIN-OLM (I am unsure of the status of this project, but if it's usable, and has enough documentation for me to get a Poisson equation up and running in a few days, by all means, recommend it).

Thanks! I was able to get things to work using standard DOLFIN & FFC. Basically, I was trying to implement interior jump conditions across interfaces, and I could set up my mesh so that the interfaces aligned with mesh facets. It took some experimentation to get it right, but it seems to work. I'm under a deadline right now, but would be happy to post an example later to show what I mean.