where *N*_{z} is the local shape function, $Nnodese$ is the number of nodes of the *e*th element, and $ui,mq,z$ is the displacement degree-of-freedom (DOF) at node *z* for phase $q=[1,2]$ in the *i*th direction. The enrichment level is *m*, with *M* being the maximum number of enrichment levels. The Heaviside function, *H*, turns on/off two sets of shape functions associated with the phases 1 and 2. The Kronecker delta, $\delta rmz$, selects the enrichment level, *r*, such that the solution at a point, *x*, is interpolated only by a single-shape function defined at node *z*. For each phase, multiple enrichment levels, i.e., sets of shape functions, are necessary if the DOFs interpolate the solution in multiple, physically disconnected regions; see Refs. [14], [24], and [25]. This generalization prevents spurious coupling and load transfer between disconnected regions of the same phase [14].