{'fraction_non_alphanumeric': 0.03044370532870993, 'fraction_numerical': 0.012355813459769301, 'mean_word_length': 4.775852203057304, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Direct Immersogeometric Fluid Flow and Heat Transfer Analysis of Objects Represented by Point Clouds

Immersogeometric analysis (IMGA) is a geometrically flexible method that enables one to perform multiphysics analysis directly using complex computer-aided design (CAD) models. In this paper, we develop a novel IMGA approach for simulating incompressible and compressible flows around complex geometries represented by point clouds. The point cloud object's geometry is represented using a set of unstructured points in the Euclidean space with (possible) orientation information in the form of surface normals. Due to the absence of topological information in the point cloud model, there are no guarantees for the geometric representation to be watertight or 2-manifold or to have consistent normals. To perform IMGA directly using point cloud geometries, we first develop a method for estimating the inside-outside information and the surface normals directly from the point cloud. We also propose a method to compute the Jacobian determinant for the surface integration (over the point cloud) necessary for the weak enforcement of Dirichlet boundary conditions. We validate these geometric estimation methods by comparing the geometric quantities computed from the point cloud with those obtained from analytical geometry and tessellated CAD models. In this work, we also develop thermal IMGA to simulate heat transfer in the presence of flow over complex geometries. The proposed framework is tested for a wide range of Reynolds and Mach numbers on benchmark problems of geometries represented by point clouds, showing the robustness and accuracy of the method. Finally, we demonstrate the applicability of our approach by performing IMGA on large industrial-scale construction machinery represented using a point cloud of more than 12 million points.


Introduction
In immersogeometric analysis (IMGA), a solid object is immersed into a non-boundary-fitted discretization of the background fluid domain that is used to solve flow physics using finite-element-based computational fluid dynamics (CFD) [1,2]. IMGA alleviates the labor-intensive and time-consuming geometry cleanup process needed to create a boundary-conforming fluid mesh to perform traditional CFD. For example, in boundary-fitted CFD mesh generation, small and thin geometric features often need to be manually defeatured, and the resulting gaps in the boundary representation (B-rep) of the solid model need to be filled to create a watertight surface representation of the solid object. This process can take up to several months for extensive industrial-scale simulations of flow over complex geometries such as tractors and trucks [3][4][5][6].
IMGA has been shown to be a practical method [2,[7][8][9] for the simulation of incompressible flow (both laminar and turbulent) around geometrically complex objects in the context of a tetrahedral finite cell approach [10]. Since the fluid domain is meshed independently, defeaturing or removing small geometric features from the immersed object is no longer necessary. The method was extended by Xu et al. [11,12] to handle moving particles, by Zhu et al. [13] to consider free-surface flows, and by Saurabh et al. [14,15] to perform industrial scale large eddy simulations using adaptive octree meshes. Xu et al. [16] recently proposed a compressible-flow version of the IMGA formulation for the simulation of aircraft aerodynamics, and Hoang et al. [17] developed a skeleton-stabilized IMGA technique for the simulation of fluid flow through a porous medium. The immersogeometric approach has also been shown as an efficient method to solve computational fluid-structure interaction (FSI) problems for heart valve applications [18][19][20][21][22][23][24][25][26]. It is also flexible enough to be automated and placed in an optimization loop that searches for an optimal design [27]. A subset of the immersogeometric FSI functionality was recently implemented as the open-source library CouDALFISh [28,29].
We have previously developed the immersogeometric method to directly use the B-rep of the computeraided design (CAD) model to perform fluid flow simulations. We were able to perform IMGA of flow over B-rep models represented using triangles [19], trimmed non-uniform rational B-splines (NURBS) [7], and finally analytic surfaces [8]. These approaches were able to handle complex geometry without needing any geometric defeaturing. However, they still worked based on the boundary representation of the well-defined CAD model (2-manifold, watertight, and with consistent normals). This paper extends the work to perform IMGA directly with point cloud representation of the solid geometry, which relaxes these restrictions (manifoldness and watertightness of the boundary representation).
In a point cloud representation, the boundary of the object is represented using a set of unstructured points in the Euclidean space with (possible) orientation information in the form of surface normals. Many geometric data acquisition methods, including optical laser-based scanners, LiDAR scanners, and even passive methods such as multi-view stereo, produce a point cloud representation of the geometry [30]. Further, there has been a high interest in performing flow simulations over as-manufactured or in-use objects rather than ideal as-designed objects for performing analysis (and to provide feedback) using digital twins [31][32][33][34][35]. In practice, we need the geometric information corresponding to the in-use physical object to analyze the flow over them. While it is possible to perform IMGA by first reconstructing the boundary surfaces and then using them for the analysis, reconstructing the surfaces to generate watertight and manifold solid models is by itself very challenging. In addition, such an approach would relax the restrictions on the B-rep CAD models obtained directly after design. The geometry no longer needs to be cleaned up to ensure it is watertight; points can be directly sampled from boundary surfaces to perform IMGA.
There are two main geometry processing steps for performing IMGA directly on point cloud representations. First, we need to perform an inside-outside test to conduct a point membership classification (PMC) on the background fluid mesh based only on the point cloud representation of the solid model. Second, we need to perform surface integration over the point cloud to impose the weak Dirichlet boundary conditions [36]. Surface integration requires estimating the surface normals and area for each point in the point cloud to be used as the Jacobian determinant for each surface element during integration. In addition, the fluid mesh needs to be adequately refined near the surface locations, i.e., the region around the points of the point cloud. In this paper, we exploit the prior art in geometric processing methods to achieve these steps.
Inside-outside evaluation for the background fluid mesh is necessary to identify points inside the solid geometry. Traditionally, a point membership classification is performed over the solid CAD model. However, these approaches are restricted to two-manifold solids with a watertight surface representation. Since point clouds are not manifold and do not have any inherent order, we propose using the generalized winding-number-based inside-outside testing approach to perform the point membership classification. This method relaxes the manifold and watertight requirements of the surface representation.
As mentioned earlier, the immersed objects used in previous immersogeometric methods consisted of tessellations, trimmed NURBS, or analytic surfaces. In the case of triangles, the Gaussian quadrature points were directly generated on the planar triangular surfaces [2]. In the case of trimmed NURBS, the NURBS parameterization was used to generate the quadrature points [7]. Similarly, the surface parametric equations were used to generate the Gaussian quadrature points for trimmed analytic surfaces [8]. In this work, it is convenient for the integration points to be co-located with the points of the point cloud. However, this approach requires calculating the Jacobian of the integration for each point in the point cloud. In this paper, we have developed a Voronoi projection-based area to compute the Jacobian determinant (or the effective area) associated with each point of the point cloud. In addition, we also need the surface normal corresponding to each point in the point cloud. We use a local hyperplane fitting algorithm to compute the surface normals.
In many industrial applications, fluid flow analyses are employed to verify and validate the efficacy of thermal control systems. The thermal analysis predicts the surface and ambient temperatures of critical components in an industrial product assembly. Previous work demonstrated the capability to simulate heat transfer using a variational multiscale method with weakly enforced Dirichlet conditions on conforming boundaries [37]. In this paper, we develop a thermal IMGA formulation to apply fixed temperature boundary conditions weakly on immersed point cloud surfaces, showing convection and conduction of thermal quantities in the flow. We first validate the proposed formulation using benchmark problems and later demonstrate the utility of thermal IMGA on a large vehicle assembly.
To summarize, the specific contributions in this paper include: 1) Methods to compute the surface Jacobian determinant and surface normals from the point cloud representation of a solid model. 2) Methods to compute the inside-outside information for classifying the quadrature points of the background fluid mesh. 3) A thermal IMGA formulation for modeling the heat transfer in flow over an immersed object. 4) Validation studies on the proposed methods to understand their efficacy in performing IMGA with a point cloud. 5) Demonstration of coupled fluid and thermal flow analysis on a large industrial-scale object represented using a point cloud.
This paper is organized as follows. We provide the mathematical formulation of IMGA for incompressible and compressible flows in Section 2. In Section 3, we provide details on processing the point cloud for computing the normals, Jacobian determinant, and inside-outside information. We also perform validation studies for the point cloud processing methods. Next, we provide the details of the flow and thermal validation studies in Section 4 and the application of the method to industrial-scale parts represented as point clouds in Section 5. Finally, in Section 6, we conclude our work and provide a few directions for future work.

Immersogeometric analysis
The immersogeometric flow analysis methodology consists of three main components: 1) The thermal fluid system is modeled using stabilized finite element methods for incompressible [37][38][39] and compressible [40][41][42] flows. 2) The Dirichlet boundary conditions imposed on the immersed objects are enforced weakly in the sense of Nitsche's method [2,16,43]. 3) To accurately capture the geometry of the flow domain, the concept of the Finite Cell Method (FCM) is employed in which the quadrature rules are adaptively refined [2,44,45]. These numerical ingredients are presented in this section. Let Ω (subsets of R d , d ∈ {2, 3}) denote the spatial domain, and Γ be its boundary. In the context of immersogeometric flow analysis, the computational domain Ω consists of two exclusive parts, the physical domain Ω phys , i.e., the fluid domain, and the fictitious domain Ω fict , i.e., the domain enclosed by solid objects. Ω phys and Ω fict are separated by the immersed boundary Γ, as shown in Figure 1. Consider a collection of disjoint elements {Ω e }, ∪ e Ω e ⊂ R d , with closures covering the computational domain, Ω ⊂ ∪ e Ω e , and let Γ be discretized into a collection of boundary elements {Γ b }. In what follows, Ω e phys = Ω e Ω phys and Ω e fict = Ω e Ω fict . Note that the finite-element discretization of the domain Ω is created without conforming to the geometry Γ, which greatly simplifies the mesh generation process. In what follows, a superscript h indicates that the variable is evaluated in the discrete space. The variational multiscale (VMS) discretization of the incompressible thermal fluid flow problem can be stated as: find the pressure p h , fluid velocity Γ Ω = Ω phys + Ω fict Ω fict Ω phys where and In the above equations, ρ, c, and κ are the density, specific heat capacity, and thermal conductivity of the fluid, respectively. σ σ σ(u h , p h ) = −p h I + 2µε ε ε(u h ) and ε ε ε(u h ) = 1 2 (∇u h + (∇u h ) T ) are the Cauchy stress and strain-rate tensors of incompressible fluid, respectively, with µ being the dynamic viscosity. I is the d × d identity matrix if not otherwise specified. f buoy is the buoyancy force per unit mass modeled using Boussinesq approximation as −β (T − T ref ) g, with β being the fluid thermal expansion coefficient, T ref being the reference temperature, and g being the gravitational acceleration. h and h T contain the prescribed traction and heat flux conditions, respectively, applied on Γ N u and Γ N T , which are the portions of Γ where the corresponding Neumann boundary conditions are applied. The fine-scale velocity, pressure, and temperature are defined by u = −τ M r M /ρ, p = −ρτ C r C , and T = −τ E r E /(ρc), respectively, where In the above equations, the stabilization parameters are given by where ∆t is the time step size, the constants C t = 4 and C I = 3 are chosen from an appropriate element-wise inverse estimation [46,47], ν = µ/ρ is the kinematic viscosity, α = κ/(ρc) is the thermal diffusivity, and G contains the information about the element size derived from the element geometric mapping from the parametric parent element to physical coordinates x(¸). The components of G are defined as The term κ IC DC formulates a discontinuity capturing (DC) operator, which provides additional numerical stability to locations where temperature gradients are large. The DC stabilization parameter [48] is given by where C IC DC is a positive constant scaling the strength of the DC operator and is set to 0.5 in this paper. The standard way of imposing Dirichlet boundary conditions is to enforce them strongly by ensuring that these conditions are satisfied by all trial solution functions, which is not feasible in immersed methods. Instead, the strong enforcement is replaced by weakly enforced Dirichlet boundary conditions originally introduced by Bazilevs et al. [36,49,50] for incompressible isothermal flows and later extended in Refs. [37,51,52] for thermal fluid flows. Let Γ D = Γ D u ∪ Γ D T be the portions of Γ on which the velocity and temperature Dirichlet boundary conditions are applied. The semi-discrete problem of the thermal fluid system can now be stated as follows: In the above, u D is the prescribed velocity on Γ D u , T D is the prescribed temperature on Γ D T , n is the unit outward normal vector, and Γ D,− is the inflow part of Γ D , on which u h · n < 0. The value ofγ can be selected as 1 or −1 which determines whether Eq. (12) is a symmetric or non-symmetric type of Nitsche's method, respectively [16,36,53]. Finally, τ IC µ and τ IC κ are stabilization parameters that need to be estimated element-wise as a compromise between the conditioning of the stiffness matrix, variational consistency, and the stability of the formulation. The choice ofγ influences the performance of the weak boundary condition operator and the selection of τ IC µ and τ IC κ , which will be discussed in detail in Section 2.1.3. An important advantage of using weakly enforced Dirichlet boundary conditions is the release of the pointwise no-slip condition at the boundary of the fluid domain. This, in turn, allows the flow to slip slightly on the solid surface and imitates the presence of the thin boundary layer that typically needs to be resolved with spatial refinement. It was shown in Bazilevs et al. [50] and Hsu et al. [54] that weak boundary conditions allow for an accurate overall flow solution even if the mesh size in the wall-normal direction is relatively large. In the immersogeometric method, the fluid mesh is arbitrarily cut by the object boundary, leaving a boundary layer discretization of inferior quality compared to the boundary-fitted counterpart. However, it was shown in Xu et al. [2] that accurate flow solutions were obtained using the immersogeometric method with a mesh resolution and refinement pattern comparable to the boundary-fitted mesh used to obtain the reference values. We believe this is partially due to the use of weak-boundary-condition formulation.

Compressible flow formulation
The compressible-flow governing equations are discretized using a streamline upwind Petrov-Galerkin (SUPG) formulation [55][56][57][58][59][60][61][62][63][64][65][66][67][68][69] augmented by a DC operator [70][71][72][73][74][75][76][77][78][79][80][81]. In what follows, Roman indices take on values {1, ..., d}, and summation convention on repeated indices is applied. In addition, we use (·) ,t to denote a partial time derivative, and we use (·) ,i to denote the spatial gradient. Let Y = [p u T] T denote the solution vector of pressure, velocity and temperature, and W = [w p w w T ] T denote the test function vector of their respective test functions. The problem can be stated as follows: find Y h ∈ S h such that for all W h ∈ V h , where and In the above, A's and K ij are the Euler Jacobian matrices and the diffusivity matrix, respectively, whose specific definitions can be found in A. Note that the superscript h for A's and K ij are dropped for the clarity of notation, even though they are evaluated in the discrete space. H contains the prescribed fluid traction and heat flux boundary conditions, and Γ H is the subset of Γ where H is specified. Res is the residual of the compressible-flow equations defined as The stabilization matrixτ τ τ SUPG is defined aŝ Note that in Eq (17), the identity matrix I is (d + 2) × (d + 2), and the notation( · ) on top of the matrices indicates that they are evaluated based on the governing equations using conservation variables. The term associated withκ C DC in Eq. (14) is a DC operator, where the DC stabilization parameter is given bŷ In the above, C C DC is a O(1) positive constant (C C DC = 0.5 in this work), andÃ −1 0 is the inverse of the zeroth Euler Jacobian of the transformation between the conservation and entropy variables (see Rajanna et al. [82,Appendix A]). Equation (18) is an extension of the δ 91 definition designed by Tezduyar and colleagues [58,83], where only the convective part of the full residual operator Res(Y h ) was employed. While the SUPG terms provide the necessary stability across a wide range of Reynolds numbers, the DC operator provides the necessary additional dissipation in the shock regions. Finally, the compressible-flow version of the weak-boundary-condition formulation is added to the weak form: where λ is the second coefficient of viscosity (λ = −2µ/3 based on Stokes' hypothesis), and τ C µ , τ C λ , and τ C κ are the stabilization parameters of Nitsche's method. In this work, the compressible gas is assumed to be calorically perfect and the specific heats at constant volume and constant pressure can be defined as c v = R/(γ − 1) and c p = γR/(γ − 1), respectively, where R is the ideal gas constant and γ is the heat capacity ratio. The pressure, density, and temperature are related through the ideal gas equation of state, p = ρRT. In addition, thermal conductivity can be calculated by κ = c p µ/Pr, where Pr is the Prandtl number.

Stabilization parameter selection for weak boundary conditions
Eqs. (12) and (19) can be either the symmetric or non-symmetric form of Nitsche's method, depending on the value ofγ being selected as 1 or −1, respectively. The symmetric Nitsche method provides excellent accuracy and robustness when the stabilization parameters are properly estimated, which sometimes requires the solution of a local eigenvalue problem [84][85][86][87][88]. However, the arbitrarily intersected elements in immersed methods make obtaining accurate solutions to the eigenvalue problem challenging. An intuitive way to circumvent this difficulty is to use a uniform value that can be determined by a trial-and-error approach to simultaneously satisfy the requirements for numerical stability and conditioning of the system matrix [1,2]. Some delicately designed algorithms have been proposed for the local selection of stabilization parameters in intersected elements [86,89] and for preconditioning of the linear system [90][91][92]. However, even with the added algorithmic complexity, careful estimations of τ's do not necessarily result in improved accuracy in L 2 errors when compared with the uniform-valued stabilization parameters [53]. Therefore, in our IMGA CFD for incompressible flows, we set τ IC µ = 10 3 and τ IC κ = cτ IC µ to take advantage of this most straightforward but effective strategy.
Our previous numerical experiments suggest that, compared with incompressible flows, IMGA CFD for compressible flows is more sensitive to the oscillations around immersed boundaries caused by the uniform, large-value stabilization parameters. With this observation, the non-symmetric Nitsche type weak boundary condition operator [16,53,[93][94][95][96] becomes an attractive alternative. The non-symmetric Nitsche formulation can be parameter-free [97][98][99][100][101], or with a stabilization added to reduce the oscillations near interfaces and to improve the L 2 accuracy [102][103][104][105]. Since the stabilization parameter in this case is not required to be larger than a specific lower-bound value to ensure stability, the estimation of the parameter becomes much simpler. In addition, the value of stabilization does not need to be very large, so it is less likely to overshadow the consistency term. Inspired by Wu et al. [27], we scale the stabilization parameters as τ C µ = τ C λ = 4ρh e n /∆t and τ C κ = c v τ C µ for compressible-flow IMGA with non-symmetric Nitsche-type weak boundary conditions. Note that here h e n = (n · Gn) −1/2 is calculated from the full element, which greatly simplifies the evaluation in intersected elements.

Adaptive quadrature near point cloud surface
Accurate numerical integration over the physical domain Ω phys with the presence of intersected elements is crucial for obtaining accurate simulation results in immersogeometric flow analysis. This is achieved in this paper via a subdivision-based adaptive quadrature rule in the context of the tetrahedral finite cell method [2]. Without changing the original background elements that support the basis functions, the method adaptively refines the quadrature distributions to create an aggregation of quadrature points in the vicinity of the immersed boundary. We refer to the elements used for generating adaptive quadrature rules as cells to avoid confusion with elements used for constructing basis functions. An intersected cell is recursively split into sub-cells until the sub-cell is completely inside Ω phys or Ω fict , or when the subdivision reaches a prescribed refinement level. Standard quadrature rules are applied to determine the set of quadrature points and weights within the sub-cells. For clarity of visualization, Figure 2 shows an illustration of the method in 2D based on adaptive sub-cells for triangles. The adaptive quadrature scheme needs two phases of point membership classification: the first classification is performed on the cell vertices to determine whether this cell is intersected, and the second classification is performed on the quadrature points to determine whether they should be active (inside Ω phys ) or inactive (inside Ω fict ). In the context of point cloud surface representation, a point membership classification method based on winding number will be presented in Section 3.

Point cloud processing
Before we extend the IMGA for point cloud representation, we begin with a few notations for representing the point cloud and its associated differential quantities. A solid geometry is represented as a smooth manifold S isometrically embedded in the Euclidean space R d , and the boundary surface of the solid geometry is represented as Γ (interacting with the domain Ω). A point cloud P ⊂ Γ could be defined as a set of points sampled from the boundary surface, i.e. P = {p ∈ Γ}. We assume that Γ is two-manifold (to ensure that S is solid) and piecewise smooth (for estimating the normals and curvatures). Further, P is sampled randomly from Γ, with the number of points denoted by |P |. For performing IMGA using point clouds, we present the following methods for (i) estimating the normal or the orientation of each point in the point cloud, (ii) enforcing Dirichlet boundary conditions, and (iii) performing point membership classification for any given domain point with respect to the boundary surface represented by point clouds for adaptive quadrature computations.

Normal computation
The raw point cloud P does not contain any information of the orientation (or the normal vector for each point in the point cloud). We estimate the normals from the raw (unoriented) point cloud by constructing a local surface proxy for each point and then estimating the normals for the surface proxy as illustrated in Figure 3. For constructing this surface, we first identify the neighborhood of each point. For a given query point p i ∈ P, where i ∈ [1, |P |], we define the neighborhood N i to be the k-nearest neighbors of the point in the set P. For fitting the local surface for each point p i , we first transform all the neighboring points, q ∈ N i , to a local coordinate system with the origin at p i . Considering the points lie in R 3 , we now define a height function for f (x, y): with Here, B (·,·) are different constants for a n-degree surface fit for the height function, z = f (x, y). The neighboring points are fit with this height function by minimizing objective function of To simplify this problem, we write it in matrix representation, where the objective function is to minimize ||MB − Z|| 2 . Here, M is the Vandermonde matrix 1, x l , y l , . We then find the solution to the linear system MB = Z, which is equivalent to fitting the best-fit surface. The guarantee for the existence of a unique solution and the accuracy of the solution to B as O(h n−j+1 ) are proven in Cazals and Pouget [106]. Once the solution for B is obtained, the normal can be computed as follows: Note that when fitting a surface, we use a different and simpler approach by constructing a covariance matrix C, which can be defined as: Neighboring points, Delaunay triangulation of Voronoi area for First, we identify the k-nearest neighbors for each point. Then we fit a local surface using the neighbors, and the computed normal of the local surface is considered as the normal for that corresponding point. On the right, we show the estimation of the area (Jacobian) corresponding to each point. Using the k-nearest neighbors computed earlier, we first use Delaunay triangulation to obtain the local triangulation surrounding the query point. Then using the barycentric coordinates, we compute the Voronoi area corresponding to each point.
Here,p represents the centroid of the neighboring points cluster N i . We can estimate the normal vector for the plane at point p i by computing the eigenvalues and eigenvectors of the covariance matrix C. The normal is defined as the eigenvector corresponding to the smallest eigenvalue [107]. Since the eigenvector corresponding to the smallest eigenvalue is also a principal component for the point cloud, this method is also known as principal component analysis (PCA)-based normal estimation.
Note that the method above does not guarantee that the calculated normals consistently point inward or outward with respect to the point cloud geometry. For a given point, we construct a minimum spanning tree consisting of the nearest neighbors and flip the calculated normal vector if its orientation is inconsistent with the neighbors [107].
While the PCA-based normal estimation performs well (using plane fitting), for better accuracy, we can fit a quadratic surface to obtain better curvature and smoothness [106,108] (also commonly known as 2-jet fitting, or simply, jet fitting). In this paper, we compare both methods (PCA and jet fitting) to understand the advantages and disadvantages of both methods.

Point cloud surface integration rules
To weakly enforce the Dirichlet boundary conditions on Γ D , we first discretize the boundary surface into elements. We then compute the weak-boundary-condition (or Nitsche) terms in Eq. (12) using the points of the point cloud and their corresponding Jacobian for each discrete element of the surface. For a point cloud represented using P, we do not need any additional discretization, and hence we consider each point p i ∈ P as a surface element 1 and its associated (single) quadrature point for performing the integration of the Nitsche terms over Γ D . However, computing the Jacobian determinant for each quadrature point for performing the integration is not straightforward. While the surface area or the Jacobian determinant is not defined in a strict sense for a point-cloud-based representation, we define an effective surface area a i for each point in the point cloud (i.e., point area or Voronoi area). This area can be associated with the geodesic Voronoi area of the point p i on the boundary surface Γ. For some point cloud acquisition processes, this is already available from the sensors [109], and for others, this can be approximated as shown in Belkin et al. [110] and Barill et al. [111]. Similar to the approach used for computing the normals of the point cloud, we construct a local neighborhood for each point p i ∈ P to obtain the geodesic Voronoi area a i as shown in Figure 3b.
The complete process can be detailed in three steps: (i) For each point p i in the set of point cloud P, we compute the neighborhood N i to be the k-nearest neighbors. (ii) Using each set of neighborhood points N i , we perform Delaunay triangulation to result in a local surface T i . (iii) The dual graph of a Voronoi diagram V i is the Delaunay triangulation T i . Therefore, computing the Voronoi area of p i can be easily performed by taking the contribution area of each triangle in T i to a i . Once we obtain the area a i corresponding to each quadrature point, we can integrate and obtain the weak-boundary-condition terms in Eq. (12).

Point membership classification
For IMGA, the attribution of a given point q in the domain Ω to the physical domain Ω phys or fictitious domain Ω fict is fundamental. In other words, we would have to determine if the point q is inside the solid S or outside it. Traditional B-rep CAD representations perform inside-outside testing by using ray-intersectionbased methods. However, these approaches assume that the surface Γ representing the boundary of solid S is 2-manifold and watertight. When the CAD geometry is represented using a non-watertight surface representation or badly oriented normals, the ray-intersection-based methods fail to determine the membership correctly. Particularly, ray-intersection-based methods do not work for representations such as point clouds (which are non-manifold).
There has been recent work on defining inside-outside tests for arbitrary non-manifold geometries such as triangle soups and point clouds using generalized winding numbers [111,112]. The winding number w(q, Γ) ∈ R can be computed as: Here, A refers to the total solid angle subtended by Γ on q, dA refers to the differential solid angle subtended by each element of Γ on q. We then project the surface onto a unit sphere to obtain the surface integral and further simplify it for a discrete point cloud represented by P as follows: p i , n i , and a i refer to the point coordinates, the normal, and the geodesic Voronoi area, respectively, for the ith point in the point cloud P. Although the winding number is defined as a real number, for a 2-manifold, watertight representation of the geometry, the possible winding number value at any given point q ∈ R 3 are positive integers. Further, if the surface does not have self-intersections, the winding number is in the binary set of {0, 1}. In other words, for any point q ∈ R 3 : However, in the case of a non-manifold boundary representation such as a point cloud representation P, the winding number w(q, Γ) is still zero for all q ∈ Ω phys with a continuous value between 0 and s + 1, where s is the number of self-intersections of the surface Γ. We illustrate how the winding number field for a non-manifold and manifold boundary representation differs in Figure 4. Therefore, for attribution of a given point q to be inside or outside of Ω phys , we threshold the winding number at 0.5, and any point with a smaller winding number is considered outside the physical domain, Ω phys .

Validation of geometric quantities
Before performing IMGA using point clouds, we first validate the estimated geometric quantities to understand their accuracy in representing the geometry. We consider three different geometries represented by point clouds for our study: (i) an icosahedron-based sampling of a sphere; (ii) a randomly sampled sphere; and (iii) a Fandisk model (a benchmark CAD geometry) containing sharp and smooth features [113].

Validation on icosphere and randomly sampled sphere
While our proposed method for estimating the geometric quantities required for IMGA works independent of the point sampling, for the first study, we consider the case of an icosphere. An icosphere is a set of points obtained by tessellating an icosahedron inscribed inside a sphere of unit radius [114]. Icosphere provides an ideal set of points sampled on a sphere such that there is a uniform distance between the points lying on the sphere. This geometric setup provides an ideal example for understanding the behavior of different parameters used to estimate geometric quantities. A more realistic geometry is a randomly sampled point cloud. We randomly generate points on a sphere by generating a set of random points and then projecting them to the sphere by dividing each point by its distance to the sphere center. Using this method, we can obtain an arbitrarily large number of points on the spherical surface. 1XPEHURISRLQWV Figure 5: The x, y, and z components of the normals estimated from 163,842 points sampled from an icosphere and the histogram of the deviation in the estimated normals (right). 'HYLDWLRQ H 1XPEHURISRLQWV For some methods used for acquiring the point cloud, the normal information is not available. Hence, we estimate the normals using the 2-jet fitting (mentioned above) for the point cloud. It takes 4.06 s to estimate normals for 163,842 points on the icosphere. Figure 5 shows the visualization of the estimated normals for the icosphere. Here, we visualize the three components of the normal vector individually. We observe that each normal component is aligned with the corresponding primary axes. A similar visualization for a randomly sampled sphere with 200,000 points is shown in Figure 6.
We also compute the deviation between the estimated and theoretical spherical normals. The theoretical normal of a given point on a unit sphere can be calculated by normalizing the vector from the sphere center to the point. The deviation is computed by subtracting the dot product between the estimated normal vector n esti and the theoretical normal vector n theo from unity (1 − n theo · n esti ). The maximum deviation in estimating the normals for 163,842 points is 1.15 × 10 −11 . A similar behavior can be observed for 200,000 randomly sampled in Figure 6 with a slightly larger deviation of 2.00 × 10 −10 .
Since the method for normal estimation relies on the local geometry, naturally, this deviation reduces with the number of points used to represent the geometry (hence capturing more intricate features of the geometry). In Figure 7, we compare the mean deviation in normal estimation and deviation in total area estimation. For normal estimation, we use 1 − n theo · n esti for estimating the deviation at each point and then use the mean deviation for comparison. However, we compare the sum of each point area for area estimation. The analytical value of the surface area of a sphere is used as the baseline for comparing the deviation. We see a linear trend for deviation of the geometric quantities with the increase in the number of points on the log scale. Also, the difference in mean deviation of normals obtained from the icosphere and randomly sampled sphere reduces with an increase in the number of points. The total area deviation for randomly sampled points is almost similar to the icosphere. Another parameter affecting the estimation of geometric quantities is the number of neighbors used to estimate them. Figure 8 shows the trend in the mean deviation of the estimated normals and the deviation in the estimated total surface area of the point cloud for different numbers of neighbors. The deviation is the minimum with an optimal number of neighbors (empirically close to 18-20 neighboring points, subject to change based on the number of points in the point cloud) and remains almost constant with a slight increase in the mean deviation of the estimated normals. This increase in deviation with an increase in the number of neighbors is because of adding more points that do not meaningfully contribute to the normal estimation. While a similar behavior is observed in normal estimation for a randomly sampled number of points, we also observe that random sampling increases the optimal number of neighbors to account for non-uniform distances between the points, specifically in point area estimation.
Another parameter to study is the degree of the surface used for fitting from the nearest neighbors. Recall that the linear plane fitting is equivalent to PCA-based fitting and is faster than the jet-fitting approach. While in terms of speed, the PCA-based plane fitting takes 0.79 s for normal estimation (compared to 4.06 s for 2-jet fitting), there is a difference of an order of magnitude in the estimation of the normals and consequently even in the area estimation. Figure 9 shows the histogram of deviation in normals estimation for icosphere and randomly sampled points using both PCA and 2-jet fitting. While there is a significant difference in the computational time and the deviation, the maximum deviation in the estimation is still less than 10 −4 , and hence is still a viable method for estimating the quantities for IMGA. Both methods, with appropriately chosen neighbors and point cloud density, have very little deviation for the sphere.

Validation of normals estimation on a point cloud sampled from the Fandisk model
As mentioned above, the normal estimation depends on the nearest neighbors for a given point. This means there could be significant deviation when the geometry has sharp edges or very fine features. To test this, we use a benchmark CAD model (Fandisk model, shown in Figure 10a). The Fandisk model consists of smooth and sharp features, thus enabling us to test the robustness of the normal and area estimation methods. While the experiments elucidated above can also be performed for the Fandisk model, we show just a few key results that provide insights for brevity.
We use the normals computed from the triangle faces of the benchmark CAD model to set a baseline for the normal vectors. This baseline preserves sharp features since the normals are completely different for adjacent triangles. However, such sharp features cannot be replicated using point clouds due to a lack of connectivity information. In Figure 10b, we show the histogram of deviations in normal estimation. While most points have very little deviations, about 20% of the total points have some significant deviation. We visualize the deviations of the normals at each point on the surface in Figure 10c. While most of the points at the center of each face have minimal deviation, all the points near the edges show significant normal deviation, which is expected. Note that the deviation of 0.30 (the maximum deviation for a point) refers to a deviation of 45 • between the two normal vectors, the average direction between perpendicular planes.

Validation of winding number computation
Since the winding number is a scalar field, we analyze the behavior by plotting the winding number at one point along the x axis. The ideal winding number for a sphere looks similar to the line-cut representing 163, 842 points as shown in Figure 11a. The winding number at x < −0.5 and x > 0.5 for a sphere centered around origin with a diameter 1.0 is 0, and 1 for −0.5 ≤ x ≤ 0.5. When fewer points are sampled, the dipole moments (see Barill et al. [111] for a detailed discussion on dipole moments) do not cancel in the local region around the boundary, which leads to local self-intersections. However, we can resolve a watertight geometry without local self-intersections by using 0.5 as a threshold to classify a given point to be within the geometry.

Robustness to noise
Often point clouds acquired from actual sensors have significant noise associated with them. While understanding the efficacy of noise removal approaches is not in the scope of this paper, we still would like to ensure that the proposed method is robust to noise. To understand the effect of noise, we add a specified Gaussian noise at each point. Figure 12 shows the increase in the deviation with an increase in the % noise added to the original point cloud. To see how this affects the geometric reconstruction, we obtain line cuts of the winding number as earlier in Figure 11b. As would be expected, adding more noise makes the reconstruction non-manifold. Figure 13 shows the renderings of the reconstructed surface mesh of the sphere. With the increase in noise, many dimples are created on the surface that could change the aerodynamics of the shape; nevertheless, the shapes generated are all 2-manifold and watertight.

Incompressible flow around a point cloud sphere
We perform mesh convergence and validation studies using icosphere point clouds representing spheres. We create an icosphere by subdividing the faces of an icosahedron and projecting the resulting nodes onto the enclosing sphere. Though the icosphere does not achieve an ideal pattern of equilateral triangles, the generated geometry is equivalent to a well-defined 2-manifold surface mesh. From this tessellated icosphere, points are spawned at the center of each triangle element and collected into a point cloud. We begin with a description of the problem setup and background meshes used for the convergence studies and validations.

Problem setup
We follow the flow around a sphere study undertaken by Xu et al. [2] with the same domain dimensions and very similar background meshes. Note that the problem is non-dimensional, and the icosphere has a diameter of one. The computational domain, boundary conditions, and immersed icosphere are shown in Figure 14. A uniform inflow velocity of one in the x-direction is set at the inlet, and the no-penetration condition is set at the lateral walls, both enforced strongly. The outflow boundary is traction free. The noslip condition is enforced weakly on the sphere point cloud using the methods described in the previous sections. The density of the fluid is set to one, and the Reynolds number (Re = µ −1 ) is defined as the inverse of the viscosity. Non-dimensional mesh settings used for this problem are summarized in Table 1. IM0, IM1, and IM2 are immersogeometric tetrahedral meshes of increasing mesh density. Note that IM2 has a comparable mesh resolution to the boundary-fitted mesh BM2 used in Xu et al. [2], which we use as a reference. We consider Re = 100, 300 and 3700 in this study, and the time step sizes in the simulations are set to 1.0 × 10 −2 , 1.0 × 10 −3 and 1.0 × 10 −3 , respectively.  As point clouds lack continuous surfaces to recursively refine surface quadrature points, a certain spatial density of point cloud relative to the background mesh's cut element size is required to prevent flow leakage. Hsu et al. [7] recommends a maximum ratio of 2 between the sizes of a tessellated surface element and the volume element cut by it. Before conducting validation studies, we derive an analogous ratio for point clouds with an icosphere geometry. This is found by converging flow quantities with increasingly dense point clouds on a fixed background mesh, specifically IM0.

Point cloud convergence with fixed mesh
We simulate flow over six different densities of point clouds on IM0 for Re = 100, 300, and 3700. The coarsest geometry has an average point spacing equal to about five times the cut element size, while the finest geometry has a point spacing of approximately 1/6 of IM0's cut element size. This tabulation aims to find the point cloud density at which our measured flow heuristics gain independence from point spacing. In Table 2, drag coefficient (C D ) and recirculation length (L/d) are listed for the steady-state Re = 100 solutions, whereas we report time-averaged drag coefficient (C D ) and Strouhal number (St) for vortex shedding cases of the two other Reynolds numbers. The drag coefficient is computed as C D = 2F D /(ρU 2 A), where F D is the drag force, ρ is the fluid density, U is the inflow speed, and A is the frontal area of the sphere. The drag force is evaluated using the variationally consistent conservative definition of traction [2,115]. The recirculation bubble length is computed as L/d, where d is the diameter of the sphere, and L is the length from the rear end of the sphere to the point where the velocity in x-direction changes sign. The Strouhal number is computed as St = f d/U, where f is the frequency of vortex shedding.
As Re increases, we observe an increase in the minimum point cloud density required to achieve stable flow solutions. As indicated by the blank cells in Table 2, the coarsest point clouds exhibit divergent solutions at higher Reynolds numbers, even with small time step sizes. As Re increases, it becomes more necessary that each cut element encloses at least one icosphere point so that each cut element adequately "feels" the forcing contribution from the immersed point cloud to prevent flow leakage.
The results in Table 2 show that convergence is consistently achieved when using a point cloud of 40,962 points. This density translates to an average point spacing valued at about 2/3 of IM0's cut element size. This agrees with the spatial density recommendations of Xu et al. [2], as this point cloud is the coarsest permutation that maintains a point cloud spacing smaller than IM0's cut element size. In subsequent simulations, we prescribe a point cloud density that ensures at least one point per cut element.

Mesh convergence and flow validation
Here, we perform convergence studies on mesh density and adaptive quadrature level and validate the flow quantities of interest at Re = 100, 300, and 3700 against those reported in Xu et al. [2]. The adaptive quadrature refinement cases are denoted as AQ followed by the level number. One icosphere cloud, with its average point spacing less than the cut element size of the finest mesh IM2, is used in all simulations in this section. The same flow qualities as in Section 4.2.2 are evaluated and the results are shown in Tables 3-5. For all cases, the convergence under adaptive quadrature refinement is clearly shown.
For Re = 100 and 300, data obtained above and including IM1 and AQ2 refinement levels are in excellent agreement with the reference values (BM2) reported in Xu et al. [2]. The results clearly show mesh convergence and demonstrate the importance of capturing the point cloud geometry in cut cells when integrating the background elements. For Re = 3700, the results show that IM2 is essential to obtaining accurate solutions. For this configuration and flow condition, there occurs a laminar flow separation near the equator of the sphere and a transition to turbulence in the wake of the object [116]. Compared to the drag coefficient, the Strouhal number appears particularly sensitive to mesh density, unable to reach a closer value until IM2. In Figure 15, the time-averaged pressure coefficient (C p ) along the sphere's upper crown line is plotted for the IM2 case. θ = 0 • corresponds to the stagnation point and θ = 180 • is the trailing point of the sphere. AQ0 shows oscillatory behavior in C p , signifying that further refined quadrature points are required in cut elements. AQ1 remedies the oscillation but largely follows the same pressure coefficient curve as AQ0. Finally, AQ2 produces a result that is in excellent agreement with the boundary-fitted reference [2]. It should be noted that the refinement levels, solution accuracy, and convergence behavior presented in this section are in full agreement with those in Xu et al. [2].

Buoyancy-driven flow
After the validation of isothermal fluid flow, we validate the thermal IMGA by simulating buoyancy-driven flows. Following the problem of natural convection in an enclosure with a heated sphere investigated by Yoon et al. [117], a static spherical boundary is situated within a sealed cubical enclosure filled with air, as shown in Figure 16. Buoyant flow is driven by a temperature differential between the sphere and the enclosure walls. This case is chosen as it requires simulation of heat transfer involving immersed Dirichlet boundaries and accurate resolution and application of buoyant forces on fluids. The sphere's boundary, with a non-dimensional radius of R = 0.2, is immersed within the cubical simulation mesh. This cube has non-dimensional edge lengths of 1, with the coordinate system at the domain's center and coordinate axes parallel to the domain's edges. The gravity acts in the −z direction. All boundaries are treated as no-slip walls. A high temperature T h = 1 is applied on the sphere, while a low temperature T c = 0 is applied on the cube domain surfaces. The normalized temperature can be defined as Θ = (T − T c )/(T h − T c ).
We simulate two cases at Rayleigh number Ra = 1000 with δ = 0 and δ = 0.25. Rayleigh number is defined as Ra = gβL 3 (T h − T c )/(να), where g = 1 is the gravitational acceleration, β = 1 is the thermal expansion coefficient, L = 1 is the length of the enclosure, ν is the kinematic viscosity, and α is the thermal  diffusivity. δ refers to the z coordinate of the sphere where x = 0 and y = 0. For both cases, the Prandtl number, Pr = ν/α, of air is used (Pr = 0.7). With the aforementioned problem setup, one can adjust the value of ν or α to achieve the desired Ra number. Three meshes are tested with each of the two δ cases, with non-dimensional element sizes listed in Table 6. The mesh sizes on the cube domain surfaces and near the immersed sphere are the "base elements" and "cut elements" sizes, respectively, and a smooth transition is achieved in the space between the boundaries. A time step size of 5.0 × 10 −3 is used in all simulations.
For each of the six simulation case permutations, the Nusselt number (Nu) is calculated at the top surface of the enclosure, referring to the wall perpendicular to the z-axis located at z = 0.5. The Nusselt number can be calculated by Nu = ∇Θ · n, where n is the unit normal vector to the wall. In Figure 17, we plot this value along a line between (0.0, 0.0, 0.5) and (0.5, 0.0, 0.5) and compare our results with the reference values reported by Yoon et al. [117].
In the case of δ = 0, for the most part, even the coarse mesh is sufficient in reproducing the temperature gradient curve. All meshes accurately simulate boundary values relative to the reference solution. Unsurprisingly, the coarse mesh produces some notable solution oscillations. Grid independence is apparent after the medium mesh, where the medium and fine meshes improve the solution accuracy and track the reference data to a strong degree. Similarly, all three meshes with δ = 0.25 produce comparable solutions near the edge of the top surface. The coarse mesh solution noticeably deviates from the reference solution at the center of the top surface. Again, the coarse mesh exhibits some oscillations throughout the solution, whereas the medium and fine meshes are successful in reproducing the reference curve.

Compressible flow over a torpedo-shaped body
In this section, we simulate laminar flow around a torpedo-shaped body represented by the point cloud data in both the subsonic and supersonic regimes. The dimension of the geometry and the computational domain can be found in Wang et al. [8]. The point cloud representation of the torpedo-shaped body is shown in Figure 18a. boundary is set to have the same total traction as the inlet. On the point cloud object, the velocity is set to zero, and the temperature is set as the stagnation temperature determined by T D = (1 + 0.5(γ − 1)M 2 )T; both conditions are enforced weakly. The heat capacity ratio γ is 1.4, the ideal gas constant R = 288.293, and the Prandtl number Pr is 0.72. Note that all quantities are dimensionless.
We perform a mesh refinement study to assess the performance of the compressible-flow point-cloud IMGA formulation. Simulations are carried out on three meshes named IM0, IM1, and IM2, from coarser mesh to finer mesh, respectively (see Wang et al. [8] for the notation and statistics of these meshes). To illustrate the mesh design, we plot a planar cut through the center of the coarsest mesh IM0 in Figure 18b. The simulations are performed using a time step size of 0.005 until a steady state is reached. Figure 19a and Figure 20a show the Mach number contour plots computed on IM2 for the subsonic and supersonic cases, respectively. The pressure coefficient distributions along the upper crown line of the torpedo-shaped body as a function of the streamwise coordinate for different meshes for the subsonic and supersonic cases are plotted in Figure 19b and Figure 20b, respectively. The pressure coefficient results are also compared with the boundary-fitted computations using a comparable mesh resolution for both cases. The results demonstrate that IM1 and IM2 meshes produce converged solutions and show excellent agreement with the boundaryfitted computations. They are also in excellent agreement with those reported in Xu et al. [16]. Note that a two-level recursive adaptive quadrature rule is employed to faithfully capture the immersed geometry and produce an essentially converged solution.

Flow over an industrial vehicle
Extending upon the practical application of IMGA to industrial geometry in B-reps and tessellated formats [2,7,8,16], we demonstrate the compatibility of this method with remarkably detailed geometry in a point cloud format. The workflow for computational analysis of an industrial product design typically begins with a B-rep CAD model representing an industrial object to be manufactured. This B-rep model primarily exists for purposes other than CFD, meaning that the geometry is often unsuitable for boundaryfitted analysis. Significant effort is invested to either "clean up" the existing B-rep or create a surrogate version with limited complexity and airtight topology. Only then may boundary-fitted mesh generation begin, accompanied by its own set of difficulties. Hsu et al. [7] and Wang et al. [8] subverted these computational analysis roadblocks using IMGA on an agricultural tractor and a tractor-trailer truck in trimmed NURBS and analytic surface formats. Compared to the tractor and truck, the geometry of this demonstration is unique in that it is overwhelmingly composed of finite-thickness shells, one type of topology targeted by CAD cleanup operations.
The geometry in this demonstration represents a medium-sized construction vehicle: a John Deere 544K Wheel Loader. Given that this is a low-speed utility vehicle, there is little concern for aerodynamic forces experienced by the vehicle, which is contrary to common vehicular CFD analysis, like that of a tractor-trailer truck. For this type of vehicle in an actual product analysis application, fluid simulation seeks to forecast how components, especially the engine and electronics, remain within thermal constraints throughout their continuous operation. Physics-coupled simulations are particularly useful to predict temperatures that the key components of the vehicle experience. Regarding the 544K demonstration that we subsequently exhibit, fluid flow similarly transpires within this low-speed incompressible flow regime, though we purposely avoid boundary conditions expressing resemblance to conditions within the operating range of the real-life vehicle. In this section, we demonstrate the utility of point-cloud IMGA applied to an industrial analysis workflow.
The generality of the point cloud format affords flexibility regarding the types of CAD formats that can be analyzed. The as-manufactured analysis advantage of point-cloud IMGA is fully realized with the virtualization of physical objects, but true as-designed analysis still presents worthwhile advantages over conventional analysis of simplified geometries. Lacking a practical 3D scan, we begin with the comprehensive CAD assembly used by Deere & Company for design and manufacturing purposes. This B-rep is ill-suited for boundary-fitted analysis as it is non-manifold with gaps, intersections, and collocations. Instead of converting it to watertight NURBS, analytic, or tessellated surfaces, the assembly is sampled into points with normal vectors and area scalars as in Figure 21 in a matter of minutes. A significant amount of manual labor is saved by circumventing the geometry simplification process typically required for boundary-fitted analysis of such complex geometry. All impermeable features within the vehicle are retained; in the ab-  sence of support for sub-scale porosity, this preliminary demonstration omits porous screens and radiators. It should be noted that all parameters utilized to simulate the 544K are unrepresentative of reality, and presented simulation results make no claims regarding the performance of the actual product as it exists physically. This exercise demonstrates point-cloud IMGA's compatibility with multiphysics simulation using complex point clouds.

Geometric pre-processing accuracy
The vehicle's point cloud is discretized by sampling its B-rep surface model. We begin with a 0.025 m spacing between points, which returns a total of 4,041,875 point cloud elements. This operation filters small geometric features, automating an otherwise intensive manual task. Similar to the icosphere, sampling of the surface geometry provides normal vectors and area scalars as analytical values to compare the quality of pre-processing approximations. Point spacing is decreased to produce three more point clouds with greater element populations of 8,135,975 points, 11,855,474 points, and 27,937,410 points. Normals estimation is performed using 4 neighboring points, and a box plot of each point cloud's collection of normal vector deviations is displayed in Figure 22.
Here, the performance of jet-fitting resembles the normal estimation on the Fandisk. A deviation of two indicates flipped normals, suggesting that the point cloud with four million elements contains an excess of incorrectly aligned normal vectors. Cazals and Pouget [106] warn that the ambiguity of normal vectors on flat surfaces is a weakness of normal vector calculation in its current form. Using a larger number of neighboring points for the normal estimation increases the error, which we theorize is due to the ubiquity of thin plates in this geometry. As the number of neighbors increases for a point on a thin plate, points on one side of the plate are grouped with points on the other side, malevolently affecting calculations for points on both sides. As the density of the point cloud increases, a point on one side of a plate is more likely to be in a neighborhood exclusively with points also on the same side of the plate. As will be discussed later, we find that the two coarsest point clouds are too sparse to produce reliable normals for stable flow simulations of John Deere 544K using point-cloud IMGA. The vehicle point cloud with 12 million elements, which has an average distance of 0.00625 m between points, is therefore used for the subsequent studies.

Immersed mesh generation
Only in rare cases is an open-source meshing tool such as Gmsh [118] suitable for boundary-fitted volume meshing of production-ready product assemblies. However, owing to the topological simplicity of generating a parallelepiped volume mesh, notable examples of open-source mesh generation software typically provide more than enough features to set up an IMGA simulation case such as this example. Here, the functionalities accessed through Gmsh allow us to create the simulation domain in a streamlined manner. 2D and 3D element size parameters help control mesh density in the flow regions of interest. For targeted refinement of tetrahedra near the John Deere 544K point cloud, Gmsh's "Attractor" field trivializes the import of point coordinates to define smooth refinement regions surrounding every point cloud element. Continuing the theme of open-source workflow, sampling discrete points from CAD surfaces is similarly straightforward with Open CASCADE [119] and freely available Python packages libigl [120] and trimesh [121].
Visible in Figure 23, the meshing strategy includes two rectangular prism refinement regions with the "outer" region completely enclosing the "inner" region, and elements adjacent to the floor have element size equal to that of the "inner" region's elements. A large majority of elements within the interior compartment must be refined, which increases the overall element count.
As mentioned earlier, the average distance between the points sampled on the B-rep model for the coarsest case is 0.025 m. Point spacing is decreased to produce three more point clouds with denser point distributions. For IMGA flow simulations, we observe that the estimated normals on the two coarsest point clouds are not sufficiently accurate to produce stable flow solutions, even though they satisfy the one-point-perelement criterion defined earlier. In the evaluation of weak-boundary-condition formulation, oscillatory normals can lead to significant instability. We find that the point cloud with 12 million points, which has an average spacing between points to be around 0.00625 m, has enough accuracy in the normal estimation and can produce a stable flow solution for our demonstration purposes. It should be noted again that requiring higher point cloud density here is not due to the point-cloud IMGA formulation but because of the need to obtain better normals estimation.
The parallelization strategy proposed by Hsu et al. [122] is employed for the IMGA simulations presented in this paper. In this strategy, the problem mesh is partitioned into subdomains by balancing the number of elements in each partition; each subdomain is then assigned to a processing core. However, this approach can sometimes create a highly unbalanced distribution of quadrature points in each partition since quadrature points aggregate in the cut elements due to the use of adaptive quadrature. In this work, we use the strategy proposed by Xu et al. [16] and weigh each element by the number of quadrature points it contains. This is then used as the metric in the graph/mesh partitioning package METIS [123] for determining mesh partitions by balancing the summation of user-defined weights. This approach produces a more balanced computational load on each processing core.

Simulation results
We solve for the airflow around the construction vehicle with the aforementioned boundary conditions and mesh. In this section, we observe an instantaneous snapshot of the flow solution. We focus on the velocity field, confirming that external airflow conforms to the immersed boundaries constructing the outer shell of the vehicle. We also look further inside the vehicle, especially in the engine compartment, to reveal the treatment of complex thin shell boundaries by inside-outside testing based on winding number. The latter portion of this section examines the vehicle interior, focusing on heat conduction and convection in the flow temperature field.  As is expected at Re = 3 × 10 6 , external flow is fully turbulent across the vehicle. Figure 24 shows the visualization of the instantaneous vortical structures colored by the velocity magnitude of turbulent flow around the John Deere 544K. The construction vehicle is visualized using a point cloud that was used to directly perform immersogeometric fluid flow and heat transfer analysis. To distinguish internal airflow from external airflow in subsequent visualizations, we substitute the surface CAD model in place of the point cloud utilized for simulation. Figure 25 shows that the front bucket attachment forces the brunt of flow movement, simulating a recirculation bubble above and downstream of the bucket. The driver cabin appears to direct some flow into the lower engine compartment, whereas flow over the top of the cabin immediately separates and heightens the large wake region behind the 544K. Thicker boundary layers are  apparent at the rear of the vehicle due to the adverse pressure gradients present there. Finally, the fan extracts air from the engine compartment and adds axial flow to the turbulent wake.
Looking at another sample of the same velocity field snapshot in Figure 26, the rotating wheel boundary conditions are evident. The surface velocity of each wheel's contact patch is equal to the freestream velocity, applied as a weak boundary condition to elements containing pertinent elements of the point cloud. The tires create minor recirculation bubbles downstream, but vortex shedding appears to be deterred by the presence of mudguards. We see another dimension of flow movement by the front bucket attachment to the degree that the vehicle's entire width is engulfed by the bucket's wake region. From this angle of the velocity field slice, the fan shows a greater level of flow extraction and interacts with the weak vortex street trailing the vehicle. Moving onto the interior compartment in Figure 27, we feature the same velocity field sample as in Figure 25 but with 2D slicing of the 544K and visualization of individual mesh elements. The bulk of airflow enters the interior section through the gap between the cabin and the hydraulic assembly. We observe moving air in the inner region, desirable for convecting heat away from critical components into the colder ambient air. We conclude this solution analysis by looking at the temperature field in Figure 28, which is purely demonstrative of point-cloud IMGA's heat transfer capabilities and unrepresentative of the real-life vehicle's actual performance. A temperature boundary condition of value T = 100 • C is weakly imposed on points belonging to the front axle hub, electronics housing, transmission case, rear axle hub, and engine. The long diagonal plate above the front axle hub impedes heat convection in that direction. As expected, the group of components in the primary engine compartment emits a greater heat load into the vehicle's wake region. The fan extracts high-temperature air that mixes with cooler air around the extremities of the compartment, pushing this air into the vehicle's wake. Again, the plates above the engine demonstrate heat conduction without air penetration across the boundaries.

Conclusions
We have presented a new method for immersogeometric fluid flow and heat transfer analysis that directly uses point cloud representation of objects. Analogous to previous work using analytic surface equations to generate surface Gaussian quadrature points, our method here repurposes point cloud elements into quadrature points. We employed computationally cheaper methods of normal vector and area generation using spatially local neighbors of points, proving that the results are appreciably accurate compared to full surface reconstruction. Prominent difficulties associated with geometry cleanup are avoided entirely, and the point cloud format trivializes immersed mesh generation using even basic meshing codes. Simulation results obtained using immersogeometric point clouds are in excellent agreement with reference solutions. Our method is painlessly applied to perform flow analysis of an incredibly complex industrial geometry, a large construction vehicle, incorporating moving point clouds and thermal fluid boundary conditions to demonstrate the utility of point-cloud IMGA in commercial and industrial applications.

Appendix A Definitions of Euler Jacobian and diffusivity matrices
In this appendix, we use the strong form of governing equations of compressible flows in 3D space (i.e. d = 3) to illustrate the definitions of Euler Jacobian and diffusivity matrices. The solution variables of compressible flows can be written using conservation variables U or pressure-primitive variables Y, defined as where ρ is the density, u i is the i th velocity component, i = 1, ..., d with d = 2 or 3 being the space dimension, e is the specific internal energy, p is the pressure, and T is the temperature. The pressure, density, and temperature are related through the ideal gas equation of state, p = ρRT, where R is the ideal gas constant. Furthermore, we assume a calorically perfect gas in this work and define the specific internal energy as e = c v T, where c v = R/(γ − 1) is the specific heat at constant volume, and γ is the heat capacity ratio. The governing equations of compressible flows can then be written as where F adv i and F diff i are the vectors of convective and diffusive fluxes, respectively, defined as F sp is the contribution of stress-power in the energy equation, defined as and S is the source term. In Eqs. (29)- (30), δ ij is the Kronecker delta, and τ ij and q i are the viscous stress and heat flux, respectively, given by where µ is the dynamic viscosity, λ is the second coefficient of viscosity (λ = −2µ/3 based on Stokes' hypothesis), and κ is the thermal conductivity. We further split the convective flux into F adv  [40].