  This work introduces a novel, fully robust and highly-scalable, $h$-adaptive
aggregated unfitted finite element method for large-scale interface elliptic
problems. The new method is based on a recent distributed-memory implementation
of the aggregated finite element method atop a highly-scalable Cartesian
forest-of-trees mesh engine. It follows the classical approach of weakly
coupling nonmatching discretisations at the interface to model internal
discontinuities at the interface. We propose a natural extension of a
single-domain parallel cell aggregation scheme to problems with a finite number
of interfaces; it straightforwardly leads to aggregated finite element spaces
that have the structure of a Cartesian product. We demonstrate, through
standard numerical analysis and exhaustive numerical experimentation on several
complex Poisson and linear elasticity benchmarks, that the new technique enjoys
the following properties: well-posedness, robustness with respect to cut
location and material contrast, optimal ($h$-adaptive) approximation
properties, high scalability and easy implementation in large-scale finite
element codes. As a result, the method offers great potential as a useful
finite element solver for large-scale interface problems modelled by partial
differential equations.
