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.
