  We consider the reliable implementation of high-order unfitted finite element
methods on Cartesian meshes with hanging nodes for elliptic interface problems.
We construct a reliable algorithm to merge small interface elements with their
surrounding elements to automatically generate the finite element mesh whose
elements are large with respect to both domains. We propose new basis functions
for the interface elements to control the growth of the condition number of the
stiffness matrix in terms of the finite element approximation order, the number
of elements of the mesh, and the interface deviation which quantifies the mesh
resolution of the geometry of the interface. Numerical examples are presented
to illustrate the competitive performance of the method.
