## 1 Introduction  
The use of fixed underlying discretizations on which an immersed boundary evolves is a common approach in both finite element analysis and isogeometric analysis (IGA [1, 2]) to address many moving boundary problems, including fluid-structure interaction [3, 4], phase evolution [5, 6, 7, 8], and crack propagation [9, 10, 11]. In these problems, the behavioral fields may vary rapidly or even discontinuously across the boundaries, thus necessitating refinement locally near the inserted boundary. There are many well-established approaches to realizing local refinement of cells. A simple and natural subdivision scheme is the space tree illustrated in Figure 0(a), that is, a quad-tree in two-dimensions (2D) or oct-tree in three-dimensions (3D). Subdivision surfaces that converge to smooth limit surfaces are one solution to local refinement of complex free-form parametric surfaces, first proposed by [PERSON] and [PERSON] [12]. Subdivision surfaces were later used for finite element analysis [13, 14] to recursively refine the initial mesh in order to achieve local refinement as well as good continuity at vertices. In the isogeometric framework, the tensor product Non-Uniform Rational B-spline (NURBS) surfaces can be refined (\(h\)-refinement [15]) through knot insertion. Anisotropic refinement has been studied in [16]. However, these methods do not address local refinement. Extensions to tensor product representations, such as the Hierarchical B-splines [17], T-splines [18; 19], and Truncated Hierarchical B-splines (THB-splines) [20] allow for local refinement. These geometry modeling techniques have been used for isogeometric analysis as well [20; 21; 22]. In general, the quadrature cells are intersected by the boundaries regardless of the field refinement near the immersed boundaries, and therefore cell-wise numerical integration will be inaccurate. As illustrated in Figure 1b, the immersed boundaries necessitate a separate integration on each side of the boundary; however, an accurate integration over such cut cells is challenging. This inaccuracy is independent of any decision as to the optimal choice of the number of quadrature points discussed in reference [23]. In general, an adaptive quadrature strategy in which the quadrature cells are locally sub-divided in accordance with the immersed boundaries is necessary. [PERSON] and [PERSON] [24] described an adaptive quadrature scheme on trimmed regions of a parametric NURBS surface. Space tree quadrature cell refinement has been extensively used in isogeometric analysis [25; 26; 27; 28; 29; 30]. However, quad/oct-tree subdivision introduces redundant quadrature cells and significantly increases the computational time during matrix assembly. For example, in the quad-tree quadrature cells shown in Figure 1b, the total number of quadrature cells increases by 276% after a three-level subdivision. The computational performance is even worse if more cells are intersected by the immersed boundaries or if a higher level of subdivision is desired. To refine cut quadrature cells, adaptive quadrature algorithms that rely on subdivision have been developed to generate boundary-fitted quadrature sub-cells in the cut elements. [PERSON] and [PERSON] [31] planarized the immersed surfaces and then subdivided the hexahedral elements into boundary-fitted tetrahedral sub-cells. [PERSON] and [PERSON] [32] approximated the boundary segment within each cut element with a polynomial, followed by a subdivision of the cut element into a few serendipity elements. [PERSON] et al. [33] developed a directly integrable blending formulation for triangles with curved sides or quadrilaterals resulting from the cut elements. While the curved element methods yields very few quadrature sub-cells and preserves geometric accuracy, a robust implementation accounting for all special cases remains non-trivial. Furthermore, these methods are not easily extendable to three-dimensional problems.