
By employing a different formulation to enforce boundary conditions, the stability problems caused by Nitsche's method are circumvented. Such alternative methodologies to enforce boundary conditions are discussed in Section REF  (e.g., the nonsymmetric Nitsche method [1]}, [2]}, [3]} or the shifted boundary method [4]}, [5]}, [6]}, [7]}). For some applications, however, these approaches lead to other complications, such as loss of consistency, symmetry, or adjoint consistency.

A fictitious domain stiffness is commonly employed with high-order discretizations in the finite cell method [8]}, [9]} (see [10]}, [11]} for reviews). This approach is mathematically analyzed in [12]}. With a fictitious domain stiffness, the volume integrals of the bilinear operator are extended into the fictitious domain multiplied by a small parameter (in practical applications this only holds for the fictitious parts of cut elements, i.e., \(\Omega _h\setminus \Omega \) ). This provides a stiffness on the fictitious part of cut elements, providing some control over the approximate solution on these. This approach has been applied to many real-world applications, notable examples of which are implant-vertebra models [13]} and additively manufactured structures [14]}. The fictitious domain stiffness is not consistent with the original problem, however. For this reason, a drawback of this approach is that it can be challenging (or for certain applications impossible) to set the parameter low enough to keep the consistency error at an acceptable level, but high enough to obtain sufficient control over the fictitious parts of elements. Besides the stability, the fictitious domain stiffness also improves the conditioning, but generally not to a level that permits the application of iterative solvers.

A minimal stabilization procedure is developed in [15]}. In this procedure, the elements are first classified as good (resp. bad) elements, i.e., elements with a sufficiently (resp. insufficiently) large overlap with the physical domain. Next, for each bad element, the normal derivative of the discrete function on the Dirichlet boundary that appears in the boundary integral of the bilinear form, is replaced by the normal derivative of an extension of the function on a nearby good element. This gives rise to a similar trace theorem as in boundary-fitted discretizations and, accordingly, the weak formulation is coercive with a Nitsche parameter that is inversely proportional to the (untrimmed) element size of the background grid such that optimal error estimates can be derived. This technique has mainly been applied to trimmed patches in isogeometric analysis, and can also be employed for coupling conditions on unfitted interfaces in overlapping multi-patch discretizations [16]}. Additionally, in [17]} it is shown that this methodology can be used to obtain inf-sup stable mixed formulations. It is to be noted that this approach does not provide control of the gradients in the (fictitious part of) cut elements themselves, and therefore only provides coercivity with respect to \(H^1(\Omega )\) , not extended coercivity with respect to \(H^1(\Omega _h)\) . As such, this approach yields stability and provides a well-posed imposition of boundary conditions, but it does not solve conditioning problems.

A recovery-based stabilization technique is commonly applied in the Cartesian grid finite element method (cgFEM, [18]}, [19]}), and is similar to the ghost penalty with patch control [20]}. In this stabilization procedure, the solution on a small cut element is weakly constrained by adding an \(L^2\) -projection term over the entire background element to the bilinear form. Instead of constraining the solution on the small cut elements to a smooth extension of the solution on neighboring elements, Richardson iteration is applied to constrain the solution to a smooth extension of the approximate solution in the previous iteration. As a benefit of this technique, [20]} mentions that it does not require operators that are not generally available in (immersed) finite element codes. Furthermore, [20]} establishes the convergence of the Richardson iteration and a bound on the condition number similar to that of boundary-fitted finite element methods, and also illustrates these aspects numerically.

A least squares stabilization term can be applied when the employed function space is (at least) \(C^1\) -continuous. In [23]} and [24]} it is demonstrated that a stable and coercive formulation can also be achieved by applying a least squares finite element term in the vicinity of the Dirichlet boundary. Because in \(H^2(\Omega )\)  the normal gradient on the boundary can be controlled by volumetric terms, coercivity is even achieved in the full space instead of only in the discrete space. A drawback of this approach is that it only applies to \(C^1\) -continuous bases. The least squares stabilization does not repair the conditioning, such that in [23]} this approach is combined with the removal of basis functions with very small supports in the physical domain and in [24]} it is combined with a fictitious domain stiffness.

