Quadrature-free immersed isogeometric analysis

This paper presents a novel method for solving partial differential equations on three-dimensional CAD geometries by means of immersed isogeometric discretizations that do not require quadrature schemes. It relies on a newly developed technique for the evaluation of polynomial integrals over spline boundary representations that is exclusively based on analytical computations. First, through a consistent polynomial approximation step, the finite element operators of the Galerkin method are transformed into integrals involving only polynomial integrands. Then, by successive applications of the divergence theorem, those integrals over B-Reps are transformed into the first surface and then line integrals with polynomials integrands. Eventually, these line integrals are evaluated analytically with machine precision accuracy. The performance of the proposed method is demonstrated by means of numerical experiments in the context of 2D and 3D elliptic problems, retrieving optimal error convergence order in all cases. Finally, the methodology is illustrated for 3D CAD models with an industrial level of complexity.


Introduction
The integration of Computer-Aided Design (CAD) and Computer-Aided Engineering has gained interest during the last two decades with the introduction of new numerical approaches as, for instance, the isogeometric paradigm [1,2] or meshfree strategies [3]. Particularly, spline-based geometric models have been found to present excellent performance for numerical simulations [4][5][6][7]. This opens the door to the formation of all-in-one design frameworks where a single geometric model is simultaneously used for parameterizing the shape of the object of interest and performing advanced numerical analyses [8][9][10][11]. The combination into one single model of both high-fidelity geometrical properties and efficient analysis performances is, however, far from trivial in general. Indeed, generating analysis-suitable geometric models for complex industrial designs requires advanced numerical tools. To achieve this goal, two different strategies can be undertaken: The first one consists in generating a fully boundary-conforming and matching geometric model such that standard analysis procedures can be directly applied. Generating those spline meshes is however a quite challenging task for geometries with complex topologies, especially when only tensor-product splines are considered [12][13][14][15]. For those cases, unstructured spline meshes [16][17][18][19][20] constitute an appealing alternative. On the contrary, the second approach aims to directly use standard CAD models which may contain non-conforming and trimmed surfaces and present geometric defects, such as water leaks or surface overlaps, and require the use of high-end analysis procedures [21][22][23][24][25][26]. Interested readers may refer to [27], and the many references therein, for an extensive review in the context of isogeometric methods. The present work falls into this second category.
A major ingredient that is commonly required to perform numerical analyses over CAD models is an efficient integration procedure which enables to evaluate integrals over complex domains such as curved polyhedrons. This is, for instance, the case when employing non-conformal analysis methods, where the geometric representation is decoupled from the discretization of the solution [28][29][30][31][32][33][34].
In this context of immersed and enriched FEM, there exist several integration approaches. In 3D, among the most common ones is worth highlighting octree subdivision [35][36][37][38] which consists in adaptively subdividing the domain of integration into sub-cells (voxels in 3D, or simple pixels in 2D). The obtained piecewise constant approximation of the underlying geometry can be improved by performing a local boundary reparameterization at the finest level of this recursion procedure via a (low-order) tessellation method [39,40]. Despite the beneficial simplicity and robustness of this decomposition-based method, it may suffer from high computational cost due to the large number of integration sub-cells, especially in three-dimensional and high-order methods.
For problems where the geometric representation of the boundary is of major importance, alternative approaches are considered. They consist in generating boundary-conforming sub-meshes which are generally non-analysis-suitable (due to the presence of hanging nodes, missing connectivity, singularities, etc.) but which are handy for integration purposes. The high-fidelity representation of the geometry boundaries, even for complex geometries, yields a highaccuracy in the evaluation of integrals. Nonetheless, even if the difficulty of generating such a high-order mesh is lower than building fully analysis-suitable boundary-conforming parameterizations, it still remains a challenging and time consuming task for complex 3D geometries. On the other hand, for two-dimensional geometries, the problem can be usually solved in a more accurate way: We refer the interested reader to the extensive survey [27].
An appealing alternative to these two approaches is the use of moment fitting techniques [41][42][43][44] in which coarse, but accurate, quadrature rules are generated for complex integration domains by tuning the positions and/or weights of the quadrature points. Nevertheless, these methods come at a price: The creation of tailored quadrature rules requires the computation of polynomial integrals over complex domains at a pre-processing stage, which calls for the use of alternative integration techniques.
Finally, there exists a fourth group of strategies for computing integrals over curved polyhedrons that lies in deriving dedicated integration rules for specific classes of integrands, as for instance polynomial functions. Indeed, it is known that integrating polynomials and other homogeneous functions over (curved) polyhedrons can be done more efficiently by invoking the divergence theorem [45][46][47][48][49]. These results can be exploited in several ways: One can perform a polynomial approximation of the integrands of interest such that the integration can be done straightforwardly [50][51][52]; those specific rules can be applied at the pre-processing stage of moment-fitting methods [42,43,53]; or can be used for creating quadrature schemes on the edges or faces of the polyhedrons for integrating the involved operators [54,55].
Within this category, worth mentioning are the recent works [49,55], where the divergence theorem is used for transforming volumetric integrals into either surface or line integrals. In [55], the authors reduced 3D integrals of general functions to 1D integrals, that are finally evaluated using fine quadrature rules. This extends the previous work [56] for the case of 2D geometries. Similarly, in [49] the complexity of 3D integrals is reduced to just vertices evaluations in the case of planar polyhedra. For the case of B-reps composed of Bézier triangles or nontrimmed B-splines patches, the authors in [49] applied the divergence theorem just once, transforming 3D integrals into 2D ones, which are approximated through standard quadrature rules.
Aligned with these ideas, in this work we present a fully quadrature-free method for integrating polynomials over general B-rep models enclosed by trimmed spline surfaces. The procedure is based on two successive applications of the divergence theorem, reducing volumetric integrals to the first surface and then line integrals, that are computed analytically up to machine precision. Hence, this can be seen as a generalization of those previous works, eliminating the need of quadrature rules. Such an approach is particularly well suited to B-Rep models as it only uses a description of the boundaries. On the other hand, handling B-Rep models with octree subdivision methods may be cumbersome as they have to evaluate if a point in the Euclidean space lies inside or outside the body for every single quadrature point, what is not always trivial.
Furthermore, we show how this integration procedure, combined with a consistent polynomial approximation step, leads to a new analysis tool for immersed isogeometric methods that skips the need of complex quadrature rules. This new integration procedure is highly-accurate (up to surface-surface intersection errors), and thus enables to handle analysis over high-order discretizations. In comparison, it is known that low-order approaches, as for instance, octree methods, require many quadrature points to keep the consistency error below the discretization error such that optimal convergence rates are retained in simulations. Consequently, it leads to high computational costs in general, which drastically reduces the benefits of employing highorder discretizations.
The developed approach is presented as follows: We firstly introduce in Sect. 2 the basics regarding immersed isogeometric analysis to further detail the scope of application of this work, and describe a consistent approximation step required for transforming the involved integrands into polynomials. Then, in Sect. 3, we discuss the geometric modeling via splines, trimming, and boundary-representation, as commonly undertaken in CAD. In Sect. 4, the proposed quadrature-free integration over B-Reps is presented. Finally, in Sect. 5, we solve elliptic PDEs and perform several numerical experiments to confirm the accuracy of the approach. Lastly, concluding remarks are summarized in Sect. 6.

Immersed isogeometric analysis
With the aim of introducing immersed methods, the used notation, and the main ideas behind this work, let us first introduce a classical Poisson's problem as our driving example. Even if the problem is presented in a 3D context, the same ideas are directly applicable to 2D problems.
Let Ω ⊂ ℝ 3 be the computational domain whose boundar y is par titioned as Γ N ∪ Γ D = Ω and Γ N ∩ Γ D = � . We also define a functional space , such that the Poisson's problem reads: find u ∈ H 1 D (Ω) solution of: where K ∈ L 2 (Ω) 3×3 is the symmetric diffusivity operator; f ∈ L 2 (Ω) and g ∈ H −1∕2 (Γ N ) are the source and Neumann terms, respectively; and n ∈ ℝ 3 ( Ω) is the outward pointing unit normal on the boundary. For the sake of clarity, and without constituting any limitation, in the problem (1) and hereinafter we assume homogeneous Dirichlet boundary conditions on Γ D . The associated weak problem can be written as: find

Immersed methods
The philosophy behind immersed methods is depicted in Fig. 1. It consists in embedding the computational domain Ω into a grid T h (Ω 0 ) of a larger domain Ω 0 , such that Ω ⊂ Ω 0 ⊂ ℝ 3 . The solution of the weak problem (3) is then discretized over a subset of the grid T h (Ω 0 ) , which allows a decoupling of the solution discretization from the actual geometry. This simple and rather straightforward procedure is the one and only mesh generation task to undertake within immersed-like approaches, making this class of methods very appealing. Indeed, this can largely ease the design-to-analysis workflow since the computational domain can be directly prescribed as a geometric model with any representation commonly used in CAD, as for instance the Boundary-Representation (B-Rep), detailed in Sect. 3. In return, the price to pay during the analysis lies in the introduction of so-called cut or trimmed elements, as illustrated in Fig. 1. This requires the integration of quantities over cut elements (as discussed in the introduction, see Sect. 1). This work focuses on this particular challenge one would face when dealing with enriched or unfitted finite element methods over B-Rep models.  As the computational domain is Ω and not Ω 0 , the partition T h (Ω 0 ) is restricted to a subset T h (Ω) as: Indeed, the grid T h (Ω 0 ) naturally splits the domain Ω 0 into three complementary partitions of elements: Fig. 1, the elements belonging to these three subsets are denoted as cut, non-cut, and inactive elements, respectively.
In this work, we limit our discussion to the case of 3D immersed isogeometric methods, nevertheless, the presentation is kept rather general and can be easily adapted to generic immersed methods [28] or particular cases as, for instance, CutFEM [31] or Finite Cell Methods [57], among others.
To solve numerically the weak problem (3) we construct a discrete spline space h (Ω 0 ) over the grid T h (Ω 0 ) as: where N p i denotes generic spline basis functions of degree p > 0 and arbitrary continuity (up to p − 1 ), and I 0 is the set of indices of those basis functions. In this work we use tensor-product B-splines, but the extensions to other cases as, e.g., hierarchical splines [58] or T-splines [59], is straightforward. For the sake of simplicity, henceforward we drop the superscript p from N p i and assume that the spline degree p is constant along the three parametric directions.
The support of some basis functions of the space h (Ω 0 ) may not intersect the domain Ω and, consequently, they do not contribute to the solution of the problem (3). Therefore, we trim the space h (Ω 0 ) as: that, as already studied in [9], holds optimal approximation properties. It is a well-known fact that the active support of some basis functions in h (Ω) ( supp{N i } ∩ Ω ) may be small, which could yield ill-conditioned operators. This is an active research topic [27,[60][61][62] that exceeds the scope of this work.
Henceforward, we assume the Dirichlet boundary Γ D to be such that Γ D ⊂ Ω 0 ∩ Ω , what grants the strong enforcement of Dirichlet boundary conditions. The opposite case ( Γ D ⊄ Ω 0 ) entails the imposition of Dirichlet conditions in a weak sense. We refer the interested reader to [63][64][65] for a dedicated discussion and to [62] for a study, in the case of spline spaces, of the inherent stability issues.
Thus, by means of the assumption Γ D ⊂ Ω 0 ∩ Ω , we can define the space: that allows us to discretize the continuous weak problem (3) as: find u h ∈ D h (Ω) solution of: where the discrete versions of the bilinear form a and the linear form b are decomposed as: The computation of the integrals over non-cut elements Q ∈ T int h (Ω) is straightforward and can be performed using classical quadrature schemes. However, the evaluation of integrals over cut elements Q ∈ T Γ h (Ω) is a challenging problem and one of the Achilles' heels of isogeometric immersed methods in 3D (see the related discussion in Sect. 1). The main contribution of this article regards the computation of those integrals through a quadrature-free approach for the case of cut elements defined as B-Rep models. This procedure is presented in Sect. 4. Nonetheless, this method is only applicable to the case in which the integrands are polynomial functions. Thus, before introducing it, in the next section the integrals in (10) are transformed such as they only rely on polynomial integrands.

Polynomial approximation of finite element operators
When considering spline discretizations over the grid T h (Ω) , the terms ∇u h , ∇v h , and v h in the operators (10) take polynomial forms ∀Q ∈ T h (Ω) . On the contrary, the datum quantities involved (i.e., K , f, and g) may not be polynomials in general. Hence, to work with integrals that only present polynomial integrands, we seek to exploit a key result introduced in [66]: It is possible to perform a polynomial approximation of the integrands in (10) without deteriorating the solution. More specifically, instead of solving the problem (9), we consider the following approximate problem: where the discrete forms in (10) are replaced by: that involves the following polynomial approximations: In the approximations above, the projection spaces must be chosen carefully, such that the introduced consistency errors do not pollute the numerical solution. Thus, by recalling [66,Theorem 13], we know that the projection of K , f, and g into spline spaces of degree q ≥ p − 1 yields a solution ū h that approximates optimally the true solution u, presenting convergence order p for the error measured in the H 1 semi-norm when the mesh size h → 0 . In [66], the authors also observed, through numerical experiments, that a projection degree q > p − 1 yields optimal convergence order also respect to the L 2 norm of the error (rate p + 1).

Remark 1
The non-polynomial nature of the quantities K , f, and g may derive from an additional mapping that further deforms the domain Ω 0 (see, e.g., [67]). A numerical example addressing this case is presented in Sect. 5.2.1 (the multi-perforated quarter of annulus). On the contrary, these quantities might be low-order polynomials (even zero-order polynomials) by construction and it is therefore not necessary to project them into polynomial spaces.
In [66], the projections (13) are performed patch-wise. Nevertheless, the same error estimates hold in the case they are carried out in an element-wise way, that is the case of this work. This results in polynomial approximations that are element-wise discontinuous. Thus, for each element Q ∈ T h (Ω) we introduce a local L 2 -projector: where ℚ q 1 ,q 2 ,…,q m denotes the space of tensor-product polynomials with degrees (q 1 , q 2 , … , q m ) along the m parametric directions.
By employing a tensor-product Bernstein basis, the projected quantities K , f , and ḡ restricted to element Q can be expressed as: ∈ ℝ , and ḡ (Q) k ∈ ℝ are the projection coefficients, and B k are tensor-product Bernstein polynomials defined over Q and with degrees = (q, q, q) such that We refer the interested reader to the Sect. 1 of Appendix A for a discussion about tensor-product Bernstein polynomials.

Operators assembly through lookup tables
In what follows, we detail the assembly of the elemental stiffness matrix and the right-hand-side vector associated to the operators (12). Thus, plugging the projections (15) into (12), a single entry of the elemental matrix and vector can be computed as: where N i , N j ∈ (Ω) are test and trial basis functions, respectively. In the expressions above it is easy to realize that all the integrands restricted to a single element Q are polynomials: Notice also that the functions N i , N j , and B k are naturally defined over the full support of each element Q, and not only over its active part Q ∩ Ω.
Finally, by exploiting their polynomial nature, the element integrals in (17) can be computed as: where B and B are tensor-product Bernstein polynomials w i t h d e g r e e s = (2p + q, 2p + q, 2p + q) a n d i,k, ∈ ℝ are element dependent constant coefficients that can be calculated by means of the Bézier extraction operators [68][69][70] associated to the spline space h (Ω).
Then, the assembly of the operators (17) reduces to the computation of the coefficients (Q) i,j,k, , (Q) i,k, , and (Q) i,k, , as well as the integrals 1 : Thus, the integrals 3D Q, and 2D Q, can be precomputed for every element Q and stored in lookup tables, that will be accessed along the assembly process to create the elemental operators, in a similar way as proposed in [66].
Nevertheless, as discussed in Sect. 1, the computation of the integrals (20) is a challenging task. In the case of non-cut elements, their evaluation is straightforward: It can be precomputed analytically for a single unit cube and subsequently adapted to every non-cut element's domain through simple transformations (translations and scalings). But in the case of cut elements the evaluation of the integrals 3D Q, and (18a) 2D Q, is far from simple. For that purpose, in Sect. 4 we propose a quadrature-free approach for the common case in which the active part of elements ( Q ∩ Ω ) can be defined through a B-Rep, discussed in Sect. 3.

Geometric modeling via boundary representation
In this section, we introduce the notation and some basic concepts about splines and geometric modeling. Hence, we provide a mathematical way of describing the active part of the cut elements Q ∩ Ω , discussed in the previous section, by means of B-Rep representations. This constitutes the basis for the integration method presented in Sect. 4.

Spline representation
Splines are considered a de facto standard in Computer-Aided Design and have been extensively studied in the literature, see for instance [71][72][73]. Among the different representation techniques available, in this work we focus on the use of polynomial mappings, and more specifically, B-spline and Bézier curves and surfaces. A B-spline or Bézier curve c can be expressed in the form: where N p i are univariate basis functions, either B-splines or Bernstein polynomials, of degree p, and P i ∈ ℝ d are their associated control points, where d is the number of spatial dimensions. In Appendix A we provide further details about Bernstein polynomials (Appendix A.1 and A.2) and Bézier geometries (Appendix A.3), that are extensively used in this work. For an in-depth discussion about B-splines, we refer the interested reader to the existing literature [71][72][73].
Using tensor-product combinations of those basis functions, B-spline and Bézier surfaces S can be constructed as: where N p 1 i and N p 2 j are univariate B-spline or Bernstein basis functions of degrees p 1 and p 2 , respectively, and P i,j ∈ ℝ d are the associated control points. For the sake of simplicity, we assumed the parametric domains of the mappings (21) and (22), Dom(c) and Dom(S) , to be [0, 1] and [0, 1] 2 , respectively.

Trimmed surfaces and boundary representations
Simple spline mappings (21) and (22) cannot represent complex real-world geometries. Instead, the multitude of these geometric objects are usually combined for such a purpose. More specifically, Boolean operations (namely, unions, differences, and/or intersections) of several geometrical entities are commonly adopted in Computer-Aided Design [71]. By means of these operations, volumetric geometries are often represented in an implicit way: the volume enclosed by a set of, possibly trimmed, boundaries surfaces. This paradigm, known as Boundary Representation (B-Rep) [74,75] and extensively used in industrial modeling tools, is considered throughout this work. As illustrated in Fig. 2, we consider a domain V ⊂ ℝ 3 , nonsimply connected in general, whose boundary V is defined by a set of connected faces F i , i = 1, … , n F , such as: The domain V may correspond to the active part of the cut elements Q ∩ Ω discussed in Sect. 2.1.
We consider the faces F i to be defined as trimmed B-spline or Bézier surfaces that are piecewise smooth. Every trimmed face F i is composed of two elements: an underlying spline surface mapping S i of the form (22), and a group of connected curvilinear segments ̂i ,j ⊂ Dom(S i ), j = 1, … , n c,i , that delimit the active region of Dom(S i ) (see Figs. 3 and 4). We denote this active region as F i ⊂ Dom(S i ).  Each segment ̂i ,j is the image of a spline curve mapping ĉ i,j ∶ [0, 1] →̂i ,j of the form (21). Thus, the boundary of the active region F i is: therefore, we can define F i as: We again refer to Fig. 3 where all the introduced quantities are depicted for an illustrative example.

Remark 2
To work exclusively with pure polynomial representations, instead of (rational) piecewise polynomials, in this work we only consider non-rational Bézier curves and surfaces. Using only Béziers does not constitute any limitation: By refining at its internal knots, any face F i , defined by means of B-spline curves and surfaces, can be easily split into a set of trimmed Bézier faces, whose underlying curves and surfaces are Béziers (see Fig. 5). On the other hand, the exclusive use of non-rational polynomials may be a limiting factor as it turns impossible the creation of exact conic curves and surfaces. This limitation can be circumvented in the case of the resolution of elliptic PDEs using immersed IGA. As discussed in [67], in those cases it is possible to approximate the geometry of the cut elements Q ∩ Ω ∀Q ∈ T Γ h (Ω) by means of Bézier curves and surfaces of degree p, the same as the solution's discretization, and still preserve optimal approximation properties.

Quadrature-free integration of polynomials over B-Reps
In this section, we deal with the integration of polynomials over a domain V whose bounding faces F i are represented as trimmed Bézier surfaces, as described in the previous section. More specifically, we seek to compute the integral: where a ∶ V → ℝ is a polynomial function. This addresses the computation of the integrals 3D Q, over cut elements Q ∩ Ω as described in (20).
The approach presented in this section consists in the successive application of the divergence theorem, as similarly done, for instance, in [45,51,54,76]. Let us first recall here the classical divergence theorem, also known as Gauss-Ostrogradsky's theorem.

Theorem 1 Let V be a subset of ℝ 3 which is compact and has a piecewise smooth boundary
where ∇⋅ is the divergence operator and n ∶ V → ℝ 3 is the outward pointing unit normal on the boundary V.
By applying the divergence theorem, the three-dimensional integral (26) is transformed into, first, surface, and then line integrals that can be evaluated analytically with machine precision accuracy. This is possible in the present context due to the polynomial nature of the successive integrands which ease the formation of the antiderivatives involved in the integration process.

From volume integral to surface integrals
To apply the divergence theorem, let us first rewrite the initial integral (26) in the same form as the one in (27): The vector field A ∶ V → ℝ 3 can be expressed as: with e i as the Cartesian unit vectors and Q i ∶ V → ℝ as the antiderivatives of a, computed by: Here 1 , 2 , 3 , 1 , 2 , and 3 are real constants, such that Since a is a polynomial function, the computation of the antiderivatives in (30) is straightforward [see Appendix A, Eq. (59)]. Furthermore, due to this polynomial nature, the continuity requirements of the divergence theorem are granted for the vector field A.
Applying the divergence theorem to (28) we obtain: where we recall that n ∶ V → ℝ 3 is the outward pointing unit normal on the boundary V . Recalling the definition of the boundary V in (23), the integral (31) can be split as: where n i are the outward pointing unit normals of the surfaces S i , i = 1, … , n F . Exploiting the parametric representation of the surfaces S i , these unit normal vector fields can be expressed as: where the normal vectors N i are computed as: In (34) we assumed that the surface parameterization is oriented such that the cross-product N i points out of V. Plugging (33) into the expression of the surface integrals I 2D in (32), they become: for i = 1, … , n F . And pulling back these integrals to the parametric domain of S i , we obtain: where the integrands r i are defined as: Interestingly, the normalization and the inversion involved in the definition of the unit normal vectors (33) vanish after the pull-back, as observed in [46], for instance. Furthermore, as the surface S i is assumed to be polynomial, then the composition A•S i is also a polynomial bivariate, but with a higher degree. Additionally, the non-normalized normal vector field N i is also a polynomial since it is computed as the product of polynomial terms (the partial derivatives of S i are polynomials). Finally, the scalar product of two polynomial vector fields, A•S i and N i , is a polynomial scalar field. Consequently, r i is a polynomial. In the case of Bernstein polynomials we refer the interested reader to Appendix A: see Eq. (75) for the details of the composition A•S i between a trivariate and a surface; and Eq. (71) for the multiplications of multivariate polynomials involved in the cross and scalar products of Eqs. (34) and (37), respectively.

Remark 3 The integrals I 2D
i in (36) are equivalent to the boundary integrals 2D Q, depicted in (20) and required for the assembly of boundary conditions in immersed methods (see Sect. 2).

Remark 4
In the case of non-trimmed Bézier surfaces, like the one depicted in Fig. 5, the integrals (36) can be easily evaluated analytically using Eq. (69).

Remark 5
In some situations the normal fields n i of the surfaces S i may be aligned with one of three the Cartesian axes. This occurs quite often in the case of immersed methods for solving PDEs, presented in Sect. 2, in which the integration domains V correspond to the cut elements Q ∩ Ω ∀Q ∈ T h (Ω) of the grid embedded in a B-Rep geometry. In that particular situation many faces F i will be planar trimmed surfaces parallel to the Cartesian axes. For those cases, a wise choice of the coefficients 1 , 2 , and 3 in the antiderivatives (30) will make the scalar product A ⋅ n i vanish, minimizing the number of two-dimensional integrals to be computed. For instance, in the case of a face F i that is perpendicular to the z Cartesian axis, choosing 3 = 0 will make the term A ⋅ n i vanish. Nevertheless, for a given domain V the coefficients 1 , 2 , and 3 must be set once and for all, and cannot be independently chosen for every face Thus, an optimal strategy may be to set 1 , 2 , and 3 independently for every V such that the largest number of surface integrals vanish for that specific domain.

Evaluating the surface boundary integrals
Applying again the divergence theorem (27), we can transform the two-dimensional integrals I 2D i in (36) into line integrals as: where m i ∶F i → ℝ 2 is the outward pointing unit normal on the boundary F i . The vector field R i ∶ Dom(S i ) → ℝ 2 is defined such that r i =∇ ⋅R i , as for instance: and 1 , 2 , 1 , and 2 are real constants, such that 1 + 2 = 1.
Splitting the boundary F i according to (24) we obtain: where m i,j ∶ Img(ĉ i,j ) → ℝ 2 are the outward pointing unit normals of the curves ĉ i,j , i = 1, … , n c,i . Exploiting the parametric representation of the curves ĉ i,j , these unit normal vector fields can be expressed as, where the normal vectors M i,j are computed as: In the previous expression, we assume that the curves ĉ i,j are oriented such as the external boundaries of F i present a counter-clockwise orientation, while the internal ones are clockwise oriented (see Fig. 4).
Plugging (41) into the expression of the line integrals I 1D involved in (40), they become: Finally, pulling back these integrals to the parametric domain of the underlying curves ĉ i,j , we obtain: where, as for the two-dimensional case, the normalization and the inversion involved in the definition of the unit normal vectors (41) vanish after the pull-back. We gather all the integrand terms together as: where As the curve ĉ i,j is a Bézier, the composition R i •ĉ i,j is a higher degree univariate polynomial. Additionally, the nonnormalized normal vector field M i,j is also a polynomial since it is computed from Bézier derivatives. Finally, the scalar product of two polynomial vector fields, R i •ĉ i,j and M i,j , is a polynomial scalar field. Consequently, t i,j is a polynomial. Therefore, the integrals (45) can be easily evaluated in an analytic way, with machine precision accuracy, without the need for quadrature schemes.
Further details for the case of Bernstein polynomials are provided in Appendix A: the composition R i •ĉ i,j between a bivariate and a curve is detailed in Eq. (75); the derivative involved in (42) is easily determined by computing the derivatives of the Bernstein basis functions as described in (58); the scalar product in (46) can be evaluated by computing the product of the individual components (Eq. (66)) and then summing the resulting expressions (Eq. (65)); finally, the 1D integrals (45) can be analytically determined using the expression (63).

Remark 6
The Remark 5 is extensible to the line integrals detailed above. In some situations (see for instance Fig. 4), some boundaries ̂i ,j may be aligned with the Cartesian axes. In those cases, the constants 1 and 2 arising in the antiderivatives (39) can be chosen such as the product R i ⋅m i,j vanishes in some of those boundaries. These constants can be chosen independently for every face integral I 2D i such as the number of 1D integrals to be evaluated is minimized.

Polynomial degree
The reader may have noticed that due to the involved compositions, A•S i and R i •ĉ i,j , as well as the products of Bézier curves and surfaces, the resulting polynomial term t i,j can potentially present a very high degree. In this section, we detail the computation of this degree, as well as the order of other terms involved in the intermediate steps.
For the sake of simplicity, hereinafter we assume that the polynomial a to integrate, as well as the Bézier mappings S i and ĉ i,j , have constant degrees along all their parametric directions and for all their components: with r ≥ 0 , s > 0 , and c > 0 , and where the polynomial spaces ℚ follow the notation introduced in Sect. 2.2. According to the definitions (34) and (42)  Analogously to the case of A , the degree of R i (Eq. 39), and its composition R i •ĉ i,j , are simply computed as: Finally, the polynomial term t i,j presents a degree: The degree of t i,j can be potentially very high what may induce numerical instabilities. Nevertheless, in the examples of Sect. 5.2.2 very high order polynomials were involved (in the order of hundreds) but no instabilities were noticed. This is due to the fact that we use Bézier curves and surfaces that are expressed in terms of Bernstein polynomials, known to be numerically more stable than other choices, as, for instance, monomial or Lagrange bases. Along with this work, we compute derivatives, integrals, additions, and multiplications of Bernstein polynomials, that are stable operations, but we never evaluate polynomials. See Appendix A for further details.

Numerical experiments
In this section, we show the performance of the presented quadrature-free approach by means of numerical experiments. In a first set of examples, in Sect. 5.1, we apply the method to the computation of simple integrals in 2D and 3D domains and compare them with standard methods based on the use of boundary-conforming quadrature schemes. Afterwards, in Sect. 5.2 we apply it to the solution of elliptic PDEs using the immersed isogeometric framework presented in Sect. 2. Figures 6 and 7 present two numerical studies used to validate the presented integration strategy. The two-dimensional case, described in Fig. 6, consists of a quadratic Bézier surface which is trimmed by three holes and a vertical curved slice. The three-dimensional case, described in Fig. 7, involves a trimmed domain defined by the intersection of a cube and a free-form cubic trivariate. We compute the mass M and the center of gravity C M of these two geometries, defined by:

Computation of integrals over B-reps
where the density is considered to be constant = 1.
Reference values of (53) are obtained through boundaryconformal quadrature schemes created by reparameterizing the interior of V with a technique similar to the one presented in [14]. This approach subdivides the domain of integration and leads to integration sub-cells. Standard quadrature rules can then be used to integrate numerically. For the sake of comparison, an overkill number of quadrature points ( 64 × 64 × 64 ) were used within each integration cell for both examples.
The obtained results are presented in Table 1. For the 2D-geometry (Fig. 6), the computed relative differences, compared with the reparameterization approach, are below 10 −15 , i.e., close to machine precision. Nevertheless, for the 3D-geometry (Fig. 7), relative differences of the order of 10 −7 were noticed.

Remark 7
We associate the larger differences in the 3D case to the intrinsic tolerances involved in some geometric operations. In this work we employ algorithms provided by Open CASCADE Technology [77] which is an open source C++ library designed for geometric modeling applications. For instance, in the specific case of surface-surface intersections between B-spline or Bézier surfaces, Open CASCADE limits the lowest tolerance to 10 −7 , which truncates the achievable accuracy and agrees with the results reported in Table 1. Similar tolerances apply to other non-linear operations. These limitations are not exclusive of Open CASCADE, as

Fig. 6
The two-dimensional trimmed geometry for the validation of the quadrature-free integration procedure Fig. 7 The three-dimensional trimmed geometry for the validation of the quadrature-free integration procedure similar issues can be found in other commercial and noncommercial geometric kernels available: Tolerances of the order of 10 −7 are more than enough for most of the applications these tools are designed for. On the other hand, we use Irit [78], an open-source geometric modeler, for other 2D operations, as it is the case of the computation of intersections between planar spline curves. The involved tolerances in Irit can be tuned according to our needs, which allows us to reach a higher accuracy for the 2D problem. In addition, it is important to remark that these limitations pollute the geometrical approximation not just for the presented quadrature-free method, but as well for other approaches, as for instance, for surface and volumetric untrimming, as previously discussed in [67]. Nevertheless, we believe that the obtained results confirm the viability of the quadraturefree integration strategy for 3D geometries.

Remark 8
For computing the quantities (53) in the case of the 2D-geometry (Fig. 6), the integration procedure can be directly started from Eq. (36), by replacing

Immersed isogeometric analysis
In this section, we demonstrate the effectiveness of the quadrature-free approach for solving PDEs in the context of the immersed isogeometric framework presented in Sect. 2.
In particular, we perform a series convergence analyses for Poisson's problem in different 2D (Sect. 5.2.1) and 3D (Sect. 5.2.2) immersed domains. Optimal error convergence rates are retrieved in all the cases. Finally, in Sect. 5.2.3, the flexibility and robustness of the proposed approach is demonstrated in the case of geometries that present a level complexity analogous to the ones found in real industrial applications. For all the studied cases, we consider the approximated Poisson's problem (11), previously discussed in Sect. 2. We adopt manufactured solutions: except for the complex geometries in Sect. 5.2.3. Accordingly, the source and Neumann terms, f and g, are defined as: The Dirichlet boundary Γ D will be defined for each particular case, and, consequently, Neumann boundary conditions will be applied on Γ N = Ω ⧵ Γ D .
The choice of such regular functions as target solutions (Eq. 54) is motivated by the aim of focusing our study on the consistency error, mainly controlled by numerical integration and geometric representation errors, while keeping the discretization error small. The approximation properties of trimmed spline spaces for the solution of elliptic PDEs have been previously studied in [67].

Poisson's problem for 2D trimmed-geometries
Let us first tackle the Poisson's problem for several twodimensional problems: -a square with a circular hole (Fig. 8), -a square with a free-form hole (Fig. 9), -a multi-perforated quarter annulus (Fig. 10).
Several solution degrees are considered: i.e., from p = 1 for the trimmed squares, and p = 2 for the annulus, to p = 4 . Importantly, the presence of conic sections require to perform some geometric approximations such that the integrals in the finite element operators involve only nonrational polynomials. As already discussed in Remark 2, to do so we rely on the results proven in [67] which reveal that approximating the elements' geometry using degree p leads to optimal numerical results. Therefore, Béziers of degree p are used to approximate the rational geometrical quantities at the element level.
In addition, it is important to remark the presence of a non-identity mapping in the problem depicted in Fig. 10. This leads to the introduction of an extra non-polynomial term in the bilinear form (see Remark 1) that is approximated through a local polynomial projection, as discussed in Sect. 2.2.
The H 1 and L 2 relative norms of the solution errors are evaluated along with the analyses. Optimal convergence rates, p and p + 1 , respectively, are retrieved for the three cases, see Figs. 8, 9, and 10. In the case of the plate with a hole case (Fig. 8), the L ∞ norm was also studied observing an optimal convergence behavior 2 .
In addition, for that particular test case, the results of the proposed immersed approach were compared against the ones obtained using a boundary-fitted method. As it can be seen in Fig. 8, for a fixed element size h, both results are comparable in terms of accuracy for all the computed norms.
In Figs. 8, 9, and 10, the H 1 and L 2 norms were computed using tensor-product Gauss-Legendre quadrature rules with p + 6 points per direction for the active non-cut elements, including the elements of the boundary-fitted method. For the cut-elements, the norms were evaluated by means of 2 The L ∞ norm of a quantity f ∈ L ∞ (Ω) is known to be lower bounded by the L 2 norm as ||f || 2 L 2 (Ω) ≤ C||f || 2 L ∞ (Ω) , where C is a constant equal to the volume of the domain Ω. the reparameterization approach already employed during the validation of the integrals computed in Sect. 5.1, using p + 6 points per direction for every integration sub-cell. On the other hand, in Fig. 8, the L ∞ norm was computed using 64 equally distributed points along each direction for every non-cut element and integration sub-cell.
The numerical solutions obtained with the quadraturefree approach enable to validate the present methodology for two-dimensional cases.
Nevertheless, it is important to remark that for the finest discretizations in the case p = 4 , the error reaches a plateau (around 10 −10 for the relative L 2 error norms, and around 10 −8 for the relative H 1 and L ∞ norms). For those cases, the discretization error becomes lower than the error induced by geometrical operations as, for instance, the slicing of the domain Ω into elements. See the related discussion in Remark 7. Similar plateaux were observed in [67,79].

Poisson's problem for simple 3D trimmed-geometries
To go one step further, we perform several analyses on three-dimensional trimmed domains. We consider four trimmed domains with several levels of complexity. Each of them consists in a cube with length L = 1000 with different trimmed regions:  R=2L/3 L=1000 -a simple planar cut (Fig. 11), -a free-form cut (Fig. 12) which is defined by a bi-quadratic surface with the following control points: -one-quarter of a cylinder (Fig. 13), -one-eighth of a sphere (Fig. 14).
As for the 2D-cases, we study the convergence rate in both H 1 and L 2 relative norms for several spline degrees. The norms are again evaluated via a reparameterization procedure. The obtained results confirm the theoretical expectations: Optimal convergence rates are confirmed. The curves created by Open CASCADE [77] during the surface-surface intersections are represented as B-splines of high degree and possibly rational. Such high order curves may lead to very high degrees during the polynomial compositions, as detailed in Sect. 4.3. As for the 2D-cases, and according to Remark 2, it is always possible to approximate at the element level those geometrical entities with Béziers of degree equal to the solution degree. This turns out to be mandatory in the case of rational curves and surfaces. For all trimmed cubes included in this section, the curves arising from surface-surface intersections were approximated at element level using Bézier curves of degree p. In the same way, for the cases in Figs. 13 and 14, the underlying rational surfaces were also approximated at element level with Bézier surfaces of degree p along both parametric directions.
Thus, for the case of the planar cut described in Fig. 11, the accuracy of the surface-surface intersections is very good due to the simplicity of the underlying geometric slices (just straight lines). Consequently, in this particular example, the convergence rates are optimal, even for p = 4 over the finest mesh (see again Fig. 11). The cube with the cylindrical removal also presents optimal convergence rates as shown in Fig. 13. We observe that for this geometry, the surface-surface intersections are also precisely computed using Open CASCADE [77]: The intersection curves are straight lines in the parametric domain of the cylindrical surface. However, in the case of the free-form trimmed cube (Fig. 12) and spherical trimmed cube (Fig. 14), the optimal convergence rates start to deteriorate for the finest discretization and p = 4 . This is due to the fact that the intersection curves are no longer straight lines in the parametric domain of the trimming surfaces, thus, they are strongly influenced by the used tolerance values, as already discussed in Remark 7. Let us mention that similar results were previously observed in [67].
Let us also study the involved polynomial degrees for the four examples included in this section according to the estimation detailed in Sect. 4.3. Applying the quadrature-free Finally, regarding the computational cost, we observe that the slicing process as well as the posterior approximation step are not particularly expensive operations. Thus, for instance, for the example of Fig. 14 with 16 × 16 × 16 elements, the intersection of the geometry with the Cartesian background grid took around 10.8 seconds running in a single core of an Intel i7-8559U 2.7 GHz processor. For that specific case, the approximation stage, took 3.3, 3.9, 4.6, and 5.8 seconds, for degrees p = {1, 2, 3, 4} , respectively. It is important to note that these operations (slicing and approximation) are easily parallelizable, and therefore, the total time can be significantly reduced by using all the cores available on modern processor architectures.

Poisson's problem on complex 3D trimmed-geometries
To show the viability of the quadrature-free approach to handle complex 3D geometries, we consider the two CAD models shown in Figs. 15 and 16. These B-Rep geometries have been extracted from the Open CASCADE database [77]. Generating a boundary-conforming volumetric parameterization of these geometries is far from a simple task. Instead, the B-Rep models are immersed into Cartesian grids (see Sect. 2). The solutions are discretized with C 1 -continuous quadratic B-spline basis functions. Again, we solve Poisson's (56) t i,j ∈ ℚ w , with w = 12p 3 + 6p 2 − 1.  15 and 16). In order to build the finite element operators, the presented quadrature-free approach is applied. The obtained solutions are depicted in Figs. 15 and 16.
The stiffness matrices associated with these examples were ill-conditioned due to the presence of small cut elements. In order to solve this issue, the associated linear systems were preconditioned using a Jacobi preconditioner as described in [61].
Regarding the computational cost of the geometric operations, for the example in Fig. 15 the slicing process and subsequent approximation stage took 55.0 and 84.8 seconds, respectively; while 46.2 and 68.5 seconds were measured, respectively, for the test in Fig. 16. As in the previous sections, these times were obtained using a single core of an Intel i7-8559U 2.7 GHz processor.
We believe that these two complex geometries highlight the viability of the developed approach to deal with designs of industrial complexity level.

Conclusions
We have presented a novel approach for the solution of partial differential equations on B-Rep geometries by means of immersed isogeometric discretizations that do not require quadrature schemes. For such purpose, we developed a new quadrature-free technique for the evaluation of integrals with polynomial integrands over B-Reps enclosed by trimmed non-rational spline surfaces.
This technique is based on two successive applications of the divergence theorem, transforming 3D integrals into line integrals that are eventually computed analytically. The involved steps require the creation and manipulation of (potentially) very high-degree polynomials. Nevertheless, we do not perform explicit evaluation of such functions, but just operations as additions or multiplications (using Bernstein bases), that are known to be more stable. The accuracy of this integration method has been verified numerically by evaluating integrals of low order polynomials over 2D and 3D domains and comparing the obtained results against reference solutions computed through boundary-conformal quadrature schemes.
To apply such an integration method to the resolution of PDEs over CAD models using immersed Galerkin discretizations, we transform the integrands of the finite element operators into polynomials. Thus, relying on [66] we create local polynomial approximations of those integrands for every element. In addition, according to [67], we also approximate at element level the rational B-splines, that may define the geometry, as non-rational Bézier curves and surfaces. This opens the door to the application of the method to B-Reps enclosed by rational splines.
The combination of the results in [66,67] theoretically guarantees the optimal approximation properties of the proposed method for elliptic problems. This approach is directly extendable to other non-elliptic problems, however, suitable approximation properties are not backed up by theoretical evidences.
A series of numerical experiments support our claims. Thus, the method's performance is illustrated by a series of elliptic problems on immersed 2D and 3D geometries, some of which present rational geometries. Optimal convergence rates were confirmed in all the cases. Finally, and to prove the potential of the method, its real applicability is demonstrated with a couple of 3D B-Rep models with an industrial level of geometrical complexity.
In this work, we particularize our study to the case of isogeometric discretizations. Nevertheless, the ideas behind are straightforwardly extendable to other immersed methods as, for instance, the finite cell method or CutFEM/IGA [31][32][33], or to other discretization techniques like XFEM or high-order virtual element [80,81] methods. In addition, the quadrature-free integration could be also handy for the evaluation of the right-hand-side integrals involved in the moment fitting techniques [53].

Appendix: A Bernstein polynomials
In this Appendix we discuss the construction of polynomials using Bernstein bases. We first introduce, in A.1, the Bernstein basis, its main properties, and the construction of univariate polynomials. Afterwards, in A.2, we discuss its generalization to the case of tensor-product polynomials. And finally, in A.3 we present the case of multi-dimensional vector polynomials. Most of the constructions detailed in this Appendix are rather classical and can be found, for instance, in [82].

A.1 Bernstein basis and univariate polynomials
Let us first introduce the Bernstein polynomials basis for a degree p ≥ 0: It is well-known that this basis constitutes an appealing alternative to monomials in terms of numerical stability when it comes to floating-point operations [82].
In addition, the Bernstein basis presents some handy properties that simplify the manipulation of polynomials. For instance, their derivatives can be easily computed as a function of lower degree polynomials. Thus, for p > 0 : In the same way, their primitives can be computed using polynomials of higher degree: that yields: In addition, any Bernstein polynomial of degree p − 1 , with p > 0 , can be expressed as a combination of polynomials of degree p as: Using the Bernstein basis, a univariate real polynomial f(t) of degree p can be expressed as: where f i ∈ ℝ . Applying (58), (59), (60), and (61) to each Bernstein basis function of the polynomial f(t), it is straightforward to compute the derivative of f(t), its antiderivative, integrate it over the domain [0, 1], and express it using a basis of degree p + 1 , respectively. In particular, due to its particular interest in this work, the integral of f(t) over the domain [0, 1] is detailed: This result can be directly applied to the computation of the integral (45), in Sect. 4. We remark that in this operation no polynomial evaluations are involved, simply the linear combination of the coefficients f i , what makes this computation stable even for high degree polynomials. Let us know introduce now a second polynomial g(t) of degree q ≥ 0:    In the case q = p , the addition (subtraction) of f(t) and g(t) it is easily computed by adding (subtracting) their coefficients: On the other hand, if q < p , g(t) must be firstly written in the Bernstein basis of degree p, applying (61) p − q times, and then the expression (65) can be directly used.
The multiplication of polynomials is another operation that is extensively used in Sect. 4. The product f(t)g(t) yields a new polynomial of degree p + q that can be computed as: Based on that, the composition of two polynomials f •g(t) is expressed as: where the terms g(t) i (1 − g(t)) p−i , i = 1, … , p , can be evaluated by means of the polynomials product expression (66).

A.2 Multivariate polynomials
The univariate construction (62) can be extended to the case of m-dimensional tensor-product polynomials as: where (p 1 , p 2 , … , p m ) are the non-negative degrees along the m parametric directions, and i = (i 1 , i 2 , … , i m ) the multiindex accounting for all the univariate indices. Operations defined for univariate polynomials, as derivatives (58), primitives (59), or degree raising (61), can now be applied for every parametric direction independently. For instance, the computation of antiderivatives along different directions is required in Eqs. (30) and (39). On the other hand, the integral of h(t 1 , … , t m ) over a domain [0, 1] m can be easily computed as: (66) where the products are computed performing multiplications between multi-dimensional scalar polynomials, detailed in Eq. (71).