Okay, there are two categories when implementing jump conditions.

**Conforming meshes**: here you basically need to add additional degrees of freedom for the nodes that are located at the interface. ElmerSolver does not support that this is done for the interface only. You could use Discontinuous Galerkin, for example, and introduce discontinuity between elements for all interfaces between elements. This seems like an overkill though. Now for this reason there is the '-discont' hack in ElmerGrid that actually adds second set of nodes to the mesh file. This means that all fields can be discontinuous and actually need a BC that can deal with the discontinuity. As Elmer wants to know the new connections to build the matrix structure there are some special "connection elements" that implement the connections between the discontinuous boundaries.

Now there are two reasons why I'm not a fan of the '-discont' option. First of all, it hasn't been written to deal with generic n-body problems where a single node could give rise to, for example, eight nodes to implement the discontinuity between all bodies. Also, currently the discontinuity exists only for some equation. For others one should implement a continuity condition to treat multiphysics problems.

**Nonconforming meshes**: Now this category gives of course more flexibility on meshing the different bodies. In the current mortar finite element implementation the connection between elements are searched usin octree and there can be even many discontinuous boundaries. As said, the downside here is that the matrix structure is ruined and the resulting matrix equation may be somewhat more difficult to solve. From my experience it is quite typical that the number of iteration is doubled. We are working with this field currently and continuity over nonconforming interfaces is now a library function that can be used for many equations. Unfortunately here we are still lacking a generic discontinuity model. That would actually involve only rather simple matrix entry but hans't yet been done. There may be some applications coming up where this could be a good solution and I see this generic path more interesting.

Of course the conforming meshes can be dealt with the routines for the nonconforming meshes too. Also the linear systems for the nonconforming meshes could be solved in various ways. Now the contraint matrix is added into the system matrix which somewhat limits the available linear solvers. One could also use Dirichlet-Neumann etc. kind of DD-methods where the individual bodies would be solved by any method and the additionally cost would be the DD-iteration needed in gluing of the different domains.

-Peter