Generate the density boundary conditions based on the phases' ones For all the phases, their boundary condition have to match at any point. It is either homogeneous Neumann or Dirichlet. In the first case: it's also homogeneous Neumann for density. In the latter: it's an interpolation of density values of each phases.