Recommendations on FEM software for implementing Nitsche's method on interfaces between matching meshes?

Suppose:

• I have two domains, $\Omega_{1} = [0, 1/2] \times [0, 1]$ and $\Omega_{2} = [1/2, 1] \times [0, 1]$.
• The domains share an interface $\Gamma = \{1/2\} \times [0, 1] = \partial\Omega_{1} \cap \partial\Omega_{2}$.
• I mesh both domains in such a way that nodes and faces on the interface match up between the two meshes.
• For convenience, $\Omega_{0} = \Omega_{1} \cup \Omega_{2} = [0, 1] \times [0, 1]$.
• The boundary of $\Omega_{0}$ is divided into two disjoint sets, $\partial{\Omega}_{0, D}$ (for Dirichlet conditions) and $\partial{\Omega}_{0, N}$ (for Neumann conditions)

I'd like to find $u: \Omega_{0} \rightarrow \mathbb{R}$ such that

\begin{align} -\nabla \cdot (\kappa_{i}(x)\nabla{u_{i}}) = f_{i}, i = 1, 2, \\ [\kappa(x)\nabla{u} \cdot n] = 0, \textrm{on $\Gamma$},\\ [u] = 0, \textrm{on $\Gamma$},\\ u = u_{D}, \textrm{on $\partial{\Omega}_{0, D}$}, \\ \kappa(x)\nabla{u} \cdot n = g, \textrm{on $\partial{\Omega_{0, N}}$}, \end{align}

where $x \in \Omega_{0}$, $u_{i} = u|_{\Omega_{i}}$ is the restriction of $u$ to $\Omega_{i}$, $[\cdot]$ is the jump operator, and $u_{D}$ and $g$ are constants.

I would like to solve this Poisson problem with interfacial conditions using Nitsche's method, as shown in 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 (paywalled link, unfortunately).

Within each subdomain, I plan on using a continuous Galerkin method with Lagrange elements of order 1 to start. If it is not possible to glue together two matching meshes to implement this method, an equivalent formulation uses XFEM and enriches the continuous spaces on $\Omega_{1}$ and $\Omega_{2}$ by using discontinuous Galerkin on the interface $\Gamma$.

Is there a library (or collection thereof) that is particularly amenable to setting up a problem like the one above? Eventually, I plan on changing the jump conditions to become more complicated (e.g., nonzero jumps; one of them will satisfy a nonlinear relationship in $[u]$ because it is a surface reaction).

I have tried using FEniCS (see my similar question on the FEniCS Q&A site), and I have gotten nowhere. DOLFIN-OLM seems to be a proof-of-concept implementation, but will not have the capability to change the jump conditions easily (or at least, not without me forking the code, and it does not look to be maintained). I don't care about the shape of the elements (triangles vs. quads does not matter). I do care about having access to the assembled stiffness matrix (or its action) to be able to perform sensitivity analyses. Parallelism and good solver backends (e.g., PETSc, Trilinos) are highly desirable, as is the ability to rapidly develop a code that solves the problem.

If there's an inelegant way to kludge this sort of thing into an existing framework (e.g., use discontinuous Galerkin everywhere and enforce the jumps at $\Gamma$), that would also be great, because I care more about getting something working than its efficiency. (I can work on the efficiency later.) If someone answering the question chooses that route, you should explain the weak form I should be using instead of the one in the Hansbo-Hansbo paper above.

• What problem are you encountering with your attempts so far? – Bill Barth Apr 10 '14 at 3:42
• The basic problem right now is trying to figure out the right primitives in DOLFIN. Meshes can't be glued together in DOLFIN, and any constructive solid geometry tools will give me a single mesh with no duplicate nodes at interfaces. I also can't define a finite element space over the interface I'd like, so I can't use the XFEM approach. Since I don't have much of an FEM background -- a collaborator has recommended this approach -- I have been trying to find a way to enforce the jump conditions on the internal interface in another way, more amenable to posing the problem in FEniCS. – Geoff Oxberry Apr 10 '14 at 4:06
• Have you looked into FreeFem++? You can either do DG directly‌​; or assemble the 2-by-2 block matrices corresponding to the two subdomains, to have DG/Nitsche face terms only across the subdomain interface. FreeFem++ has good interfaces to parallel solvers and can do the processor-local assembly in parallel. The only thing that might not be good is that every processor must have the full matrix which means a call to MPIAllReduce or MPIBroadcast. – Hui Zhang Apr 11 '14 at 8:45
• @BillBarth: I figured it out. There were some FEniCS quirks that made things tricky. Also, apparently, a cG basis can be used to define a weak form that includes dG-like terms (jumps and averages); I used these terms to penalize a flux jump over an interface. Time permitting, I'll write it up. – Geoff Oxberry Apr 11 '14 at 11:52

My suggestion is to just write your own code using PETSc instead of using an existing FE library. Parallel assembly/solve is the most complicated part of an FE code and PETSc takes care of it. The rest of it is simple anyway, maybe a few hundred lines of straightforward C/Fortran. Besides PETSc now has DMPLEX for managing meshes (though it is still in development).

If you don't want AMR and are not working on academic type problems then you will be much better off with your own code built directly on top of PETSc.

Like you, I explored the possibility of using an existing FE lib for my own code a few years back but each FE library came with its own set of quirks/limitations/issues and doing complicated things sometimes is either impossible or much more difficult (and time consuming) than doing it yourself.