Isogeometric analysis based on geometric reconstruction models

In isogeometric analysis (IGA), the boundary representation of computer-aided design (CAD) and the tensor-product non-uniform rational B-spline structure make the analysis of three-dimensional (3D) problems with irregular geometries difficult. In this paper, an IGA method for complex models is presented by reconstructing analysis-suitable models. The CAD model is represented by boundary polygons or point cloud and is embedded into a regular background grid, and a model reconstruction method is proposed to obtain the level set function of the approximate model, which can be directly used in IGA. Three 3D examples are used to test the proposed method, and the results demonstrate that the proposed method can deal with complex engineering parts reconstructed by boundary polygons or point clouds.

The most important part of IGA is the parameterization of analysis domain with spline models. In practice, representing a complex model by a complete tensorproduct non-uniform rational B-spline (NURBS) patch, which impedes IGA to be utilized in engineering models with complex topological structures, is quite difficult. Although the use of multiple patches and multiple point constraints in IGA can deal with a model with complex geometry [18,19], it is very complicated and cannot guarantee the high continuity between patches. Moreover, CAD systems usually use boundary representation (Brep) as an industry standard to represent geometric objects [20]. Thus, the 3D splines for 3D problems do not exist in CAD systems. To solve the above problems, researchers attempt to find solutions from different aspects and obtain some achievements. Zhu et al. [21] proposed a surface parameterization method using the one-step inverse forming and used Coons surface parameterization to rebuild the CAD models for IGA. Jaxon and Qian [22] presented the IGA on triangulations, where the geometry and physical fields were triangulated by rational triangular Bernstein-Bézier splines, to achieve the automated meshing of complex objects. Martin et al. [23] used discrete volumetric harmonic functions to parameterize a 3D model, and then a single trivariate Bspline was used to fit the parameterized model and made it suitable for IGA. However, the above methods are all based on discrete models, which cannot avoid meshing techniques. Xu et al. [24] investigated the volume parameterization problem of the multi-block computational domain for IGA, where a constraint optimization problem of minimizing quadratic energy functions was solved to generate the multi-block volume parameterization without self-intersections. Zuo et al. [25] proposed an IGA method for constructive solid geometry (CSG) models, which decomposed a complex model into simple primitive models represented by NURBS and then constructed an analysis-suitable CSG model by Boolean operators, such as union, subtraction, and intersection operators. In addition, the use of other splines that easily represent complex geometries (e.g., T-splines [26,27] and U-splines [28]), is also an alternative for IGA. Although these approaches can alleviate the difficulty of IGA model generation, they have to introduce complicated data structures incompatible to general CAD systems or solve extra problems, thereby imposing extra difficulties and burdens on the preprocessing of IGA.
Another type of method to deal with complex geometric models in IGA is the use of trimming techniques. Kim et al. [29] proposed an IGA for trimmed CAD surfaces, where the trimmed elements were decomposed into triangular or quadrilateral elements with only one trimmed edge, and the integration of trimmed elements was implemented by coordinate transformation. Wang et al. [30] further subdivided the types of trimmed elements, proposed corresponding integration strategies, and implemented the IGA with multiple intersected surfaces. Nagy et al. [31,32] proposed a numerical algorithm to construct efficient quadrature rules for trimmed elements by designing the locations and weights of quadrature points and applied it to isogeometric finite element analysis (FEA) and isogeometric boundary element analysis. Ruess et al. [33] proposed a reduced integration method for trimmed elements where a quadtree-based decomposition was used to locally subdivide trimmed elements, and a recursive bisection approach was used to identify the integration points inside the trimmed elements. Marussig et al. [34] presented extended B-splines as a stable basis for IGA with trimmed elements, in which the B-splines may cause illconditioned system matrices substituted by extensions of untrimmed B-splines. Guo et al. [35] proposed an isogeometric shell analysis with trimmed surfaces based on the standard for the exchange of product model data (STEP) format, where mapping and subdivision techniques were used to implement the integration of trimmed elements. More information about trimming techniques in IGA can be found in a recent comprehensive review written by Marussig and Hughes [36], which includes the design, data exchange and computational simulation, and challenges of trimming in IGA.
Apart from the aforementioned approaches, the embedded domain method that immerses the analysis model into a simple non-boundary-fitted background mesh is also an appropriate candidate for applying IGA to complex models. The finite cell method (FCM) [37,38] embeds the original analysis domain in a geometrically larger domain with simple geometry as the approximate solution domain, and such approximate domain can be easily represented by standard tensor-product NURBS basis functions. In addition, for geometric models represented by voxel models (e.g., the models from medical imaging technologies), the IGA cannot be used because no CAD model is created, and such a problem can be solved by combining the FCM and IGA. For a comprehensive review of FCM and its application in IGA, details can be found in the article of Schillinger and Ruess [39]. Similar to the FCM, another embedded domain method called immersogeometric analysis was proposed by Kamensky et al. [40], which directly analyzes a spline-based structure by immersing it into a non-boundary-fitted analysis mesh. This method was first applied to the FSI simulation of bioprosthetic heart valves and further extended to analyze turbulent flow around complex geometries by a tetrahedral FCM [41]. In embedded domain methods, capturing geometry, which is the core for classifying the quadrature points that determine the accuracy of analysis, is significant. For example, the quadrature points in a trimmed element can be divided into two types: (1) points in the physical domain for numerical integration, and (2) points in the fictitious domain that should be discarded. To improve the performance of integration, a local refinement scheme is used to subdivide the trimmed elements into subcells, and the quadrature point classification will be repeatedly applied to such subcells. For complex CAD models, the point membership classification and a large number of integration points are time-consuming. Furthermore, since embedded domain methods are non-boundaryfitted, special approaches will be used to enforce Dirichlet boundary conditions and visualize analysis results. Constructing geometric models from draft to solid models in CAD systems is called forward modeling, where all the geometric parameters must be known before the model construction. In many industries, the geometric models are directly obtained from the physical models, which is called reverse modeling (i.e., reverse engineering). For example, the automotive designers usually work on physical prototypes, such as clay or wood models, to achieve a favorable shape in automobile industries [42], and the geometric models of medical models (e.g., organs and bones) are created by reverse modeling based on the digital images from computed tomography, magnetic resonance imaging, and other medical scanning processes [43]. In reverse engineering, the procedure of generating geometric models is called model reconstruction, and such models are called reconstruction models. Model reconstruction can be divided into two types: (1) surfaced-based methods and (2) feature-based methods [44]. Comparing to the featurebased methods, the surfaced-based methods are more mature and adopted by commercial software, such as CATIA, Geomagic, and Rapidform. The most important part of the surfaced-based methods is the isosurface generation. In 1987, marching cube (MC) algorithm was proposed by Lorenson and Cline, which has become a standard of isosurface generation because of its simplicity and high speed. However, the original MC algorithm may cause ambiguity problems (i.e., different isosurfaces may be obtained from the same voxel data) [45]. To eliminate the ambiguous problems, many new algorithms were proposed by different research groups (e.g., extended MC algorithms [46][47][48] and marching tetrahedra algorithms [49][50][51]). Nowadays, reconstruction models can accurately represent the physical models. Hence, the numerical analysis (e.g., FEA and IGA) based on such reconstruction models can also accurately predict the performance of the physical objects. Zander et al. [52] directly used a 3D finite cell toolbox FCMLab to work on a femur defined by voxel model. Verhoosel et al. [53] proposed an image-based adaptive IGA and applied it to the micromechanical modeling of trabecular bone, where the goaladaptive FCM was used to compute elastic properties. Recently, Kudela et al. [54] completed the direct structural analysis of domains defined by oriented point clouds based on the FCM. However, the adaptive FCM will increase the complexity of the analysis, as well as the computational cost.
In this paper, an IGA framework based on reconstruction models is proposed, which can be directly used to analyze complex models with arbitrary geometries, as well as the point cloud data. The outline of this paper is organized as follows: In Section 2, a brief overview of NURBS-based IGA is given; in Section 3, a model reconstruction method suitable for IGA is established; in Section 4, numerical quadrature for elements in detail is presented; in Section 5, numerical examples to verify the efficiency and accuracy of the proposed method are modeled and analyzed; finally, conclusions are summarized in Section 6. consists of a sequence of non-decreasing real numbers as where p is the order of the B-spline, ( ) is variable for the knot coordinate in parameter space, and n is the number of basis functions. The whole interval is a patch, and the knot interval is a span. Referred to the Cox-de Boor formula [56], the B-spline basis functions are recursively defined as where is the ith B-spline basis function of p-degree for the coordinate in parameter space.
One important property of the B-spline functions is that they constitute a partition of unity.
and another important property is that the continuity between knot spans achieves where k is the multiplicity of the knots.
Two-dimensional NURBS basis functions from the knot vectors and are constructed as where is the weight value corresponding to the tensor product . In the NURBS-based IGA, the above NUBRS basis functions are used for representing the physical fields as x(ξ, η) where can be a displacement, force, or other physical values at the position whose coordinates in the parameter space is , and and are the basis function vector and the physical value vector corresponding to the control points acting on the position , respectively. K e In this work, the NURBS-based IGA is used to solve static linear elasticity problems. The element stiffness matrix may be written as [14] where B and D are the strain-displacement matrix and the stress-strain matrix, respectively, and , , and are the physical, NURBS parametric, and integration domains of the element, respectively. Jacobian and

Model reconstruction method
In CAD systems, the B-rep is used to represent a geometric object, which means that only boundary information is included in the CAD models (e.g., a 3D model represented by its boundary with only surface information). The B-rep representation is quite different from the computational mesh of finite-element-based computer-aided engineering (CAE), and Hughes et al. [1] illustrated that up to 80% of overall analysis time is devoted to mesh generation in the automotive, aerospace, and ship building industries. Additionally, the 3D splines for 3D IGA cannot be directly obtained from CAD systems due to the B-rep representation.
In this section, we propose a model reconstruction method that can construct an analysis-suitable model of arbitrary geometry for IGA. To reconstruct an analysissuitable model, the boundary of the original model is first discretized, similar to the boundary discretization of the boundary element method. Subsequently, the model is embedded into a regular domain, which can be represented by complete NURBS basis functions, and the embedded domain is divided into subcells (i.e., the elements for IGA). Fortunately, the triangular mesh is applied in CAD systems for visualization as the boundary discretization mesh to avoid the use of additional mesh generation approaches [57]. Based on the boundary discretization and model division, the elements are classified into different types, and the empty elements are compressed. The following subsections will give the details of the proposed method.

Elements of embedded domain
In embedded domain methods, such as FCM and immersogeometric analysis, a simple unfitted structural mesh is used to approximate the solution field to avoid generating boundary conforming meshes. As shown in Fig. 1, the elements of the embedded domain are divided into three types: 1) fictitious elements, 2) solid elements, and 3) trimmed elements. The identification of element types is one of the key issues for embedded domain methods.
The identification of element types is actually a series of inside/outside tests for element positions. One of the simple and efficient approaches to implement such tests is the ray casting algorithm [58,59], which casts rays from a point and count intersections with a closed geometry and has been successfully applied to irregular geometries [40,60]. For 2D models, the inside/outside test can be efficiently implemented by subdividing the closed curve into segments and by applying the ray casting algorithm. However, for 3D models, the boundary surfaces will be subdivided into polygons, and the ray casting algorithm will be more tedious and time-consuming.
Taking a bracket as an example (see Fig. 2), the embedded domain is a regular cuboid. For an element with n nodes ( in Fig. 2), we can define a rule to identify the element type as follows: (i) fictitious element: All the nodes are outside the bracket domain; (ii) solid element: All the nodes are inside the bracket domain; (iii) trimmed element: Some nodes are inside the bracket domain and the other nodes are outside the domain. If the bracket surfaces are subdivided into m polygons (i.e., triangles in Fig. 2), then the ray casting algorithm has to implement times of point-polygon intersection tests to identify the element type for one element. When the number of elements and the number of polygons are large, it is very time-consuming to implement such element classification.

A high-efficient element classification method
As mentioned in Section 3.1, it is inefficient to use ray casting algorithm for element classification of large-scale problems. To solve this problem, we will propose a highefficient element classification method in this section. The element classification method consists of three main steps: 1) embedded domain construction, 2) identification of trimmed elements, and 3) element classification.

CAD surface model
For a CAD model, the main steps of the element classification method are shown in Fig. 3, and a description is given as follows: (1) A polygon model is first created by dividing the boundary of the CAD model into polygons. In this work, the triangles for the CAD visualization are used as polygons, which can be directly applied to the visualization of results. Based on such boundary polygon model, a bounding box can be generated as the embedded domain, which is divided into regular subdomains (i.e., the elements of the embedded domain). The element size should be larger than the polygon size.
(2) According to the relative position to boundary polygons, the elements of the embedded domain are initially classified into two types: (i) trimmed elements (i.e., the elements cut by the boundary polygons), and (ii) untrimmed elements. As shown in Fig. 4, if any vertex or the center of a polygon is in an element, then the element is defined as a trimmed element; otherwise, it is an untrimmed element. (3) The minimum distance (MD) from a vertex of a trimmed element to the boundary polygons can be calculated using the algorithm of the MD between a point and a triangle [61], which will be given in Section 3.3. The trimmed element and its adjacent elements are considered an element group. Considering that the element size is larger than the polygon size, we only need to calculate the distances from the vertex to the polygons in the element group and select the minimum one as the MD. The MD can be converted into a sign minimum distance (SMD) by comparing the directions between the vector from the vertex to the polygon plane and the outward vector of the polygon . If , then the vertex is in the solid domain, and the sign distance is positive; otherwise, , the vertex is in the fictitious domain, and the sign distance is negative.
After the sign distances are obtained, the elements can be divided into three types by iteratively expanding the signs of the MDs to all vertices of untrimmed elements (Fig. 5). The sign expanding rule is: (i) If the arbitrary vertex of a untrimmed element is in the solid domain (i.e., SMD > 0), then all vertices of this element should be in the solid domain, and this element is a solid element; and (ii) if the arbitrary vertex of a untrimmed element is in the fictitious domain (i.e., SMD < 0), then all vertices of this element should be in the fictitious domain, and this element is a fictitious element.

Point cloud model
The main steps of the element classification method for point cloud models are given as follows: (1) Similar to the embedded domain construction of the CAD surface model, a bounding box that contains all the points is generated as the embedded domain, as shown in Fig. 6. The elements are obtained by dividing the embedded domain into regular subdomains. Normals of the point cloud can be directly obtained from the 3D scanner.
(2) The identification of trimmed elements for point cloud models is straightforward. If an element contains at least one point, then this element is a trimmed element, otherwise, it is an untrimmed element.
(3) The MD from a vertex of a trimmed element to the point cloud is obtained by calculating vertex-to-point distance for the points in the trimmed element and its neighbor elements, and then the minimum one is selected as the MD from the vertex to point cloud. Since the outward normals of the points can be directly obtained from general 3D scanners, the SMD is obtained by the vector from the vertex to the nearest point and the outward vector of the point, and the elements can be also divided into three types similar to CAD surface model cases, as that in Section 3.2.1.

Minimum distance between a point and a triangle in 3D
T(s, t) = A + sl In Section 3.2, the MD between a point and a triangle plays an important role in element classification. In this study, we used the algorithm proposed in Ref. [61], which is briefly summarized herein. As shown in Fig. 7, a triangle △ACE (vertices are A, C, E) can be parametrically represented by while , where and are the vector from A to C and the vector from A to E. The MD between a point P and this triangle is to find the values so that is the triangle point closest to P. The distance from P to any point on the triangle ( ) is calculated by domain, as shown in Fig. 8. Assuming , , and , the criteria determining which region contains the closest point are given in Table 1. When the region of the closest point is obtained, the MD is computed in terms of the region, and the detail algorithm can be found in Ref. [61].

Model reconstruction based on sign minimum distances
Taking the SMDs of vertices as level set function (LSF) values, we reconstruct an approximate model that is implicitly represented by the zero level set isosurface, and When the level set model is obtained, we can construct the explicit geometric model by the marching tetrahedra algorithm [49]. In the marching tetrahedra algorithm, a cube is split into 6 tetrahedrons, as shown in Fig. 9, and a tetrahedron isosurface is constructed for each tetrahedron to reconstruct the model. For a tetrahedron, if symbol 1/0 indicates that a vertex is inside/outside the boundary polygon, then a binary number of 4 digits can be used to represent different cases of a tetrahedron cut by a planar facet, as shown in Fig. 10. Since the SMDs of the vertices are known, the symbol 1/0 can be directly obtained. In marching tetrahedra, the planar facet approximate to the zero-distance isosurface is calculated for each tetrahedron. The facet vertices are the points whose MD to the boundary polygon is zero, and these zero-distance points can be obtained by linearly interpolating the SMDs of vertices along each edge as Eq. (9): where d 1 and d 2 are the SMDs of the two vertices of an edge, x 1 and x 2 are their corresponding coordinates, and x 0 is the coordinate of the zero-distance point. Figure 11 gives an example of tetrahedron trimming based on SMDs.
The trimmed tetrahedron in Fig. 11 can be regarded as a reconstruction model constructed by marching tetrahedra algorithm based on the MD. Notably, the explicit geometric model is unnecessary in the IGA method proposed in this paper, but the explicit model is necessary for conventional FEA or manufacturing.

Element renumbering and reusing approaches
When an embedded domain is constructed, all the elements have to be numbered to complete the element   classification, as shown in Fig. 3. However, the fictitious elements that are fully outside the geometric model should be removed in IGA. Thus, the trimmed and solid elements (i.e., the reserved elements) need to be renumbered.
In this work, a simple approach is used to loop over the original elements and renumber the trimmed and solid elements (Fig. 12), and an element renumbering example is given in Fig. 13. In this approach, each trimmed/solid element has two different IDs whose mapping relation is recorded in a list so that the original data (e.g., coordinates of nodes and element connectivity) can be directly reused in IGA based on the mapping relation.
One advantage of the embedded domain method is that the domain can be divided into uniform elements. Therefore, the elements with the same property can be classified into a group, and only one stiffness matrix needs to be calculated. Considering Fig. 13 as an example, only one stiffness matrix, which can be applied to all the other elements, needs to be solved since all the solid elements have the same geometry and property (Fig.  14). In IGA, apart from the material property, the distribution pattern of control points is also an element     property, since the distribution pattern of the control points may be different for the elements that are located at different positions. Figure 15 shows a comparison of computational models (2D problem) between FEA and IGA, where the distribution pattern of nodes is the same for all elements, but the distribution pattern of control points may be different. Due to the uncertainty of trimming cases, the stiffness matrices of the trimmed elements are still calculated individually. Based on the aforementioned element reusing approach, the computational efficiency of element stiffness matrices will be improved, which is especially useful for large-scale problems.

Numerical quadrature for elements
In embedded domain methods, the untrimmed elements are integrated by standard Gauss quadrature method. The trimmed elements are usually integrated by the adaptive quadrature algorithm [39,40]. However, identifying the    Gauss points in the solid domain, especially for the complex trimming cases, is very difficult and timeconsuming. In this work, the Gauss points in the solid of the trimmed element can be easily identified due to the known LSF values of the element vertices. As such, the trimmed elements can be integrated efficiently.

Regular element integration
Since the embedded domain is regular, the regular elements can be cubes and the standard Gauss quadrature method can be used to complete the element integration. When the element is a regular cube, the Jacobian determinants and can be efficiently computed by and V Ω e VΩ e V Ω e where , , and are the volumes of the element in physical, NURBS parametric, and integration domains, respectively. For 2D problems, the element volume should be replaced by the element area.

Trimmed element integration
The adaptive quadrature algorithm [39,40] is one of the effective approaches for trimmed element integration, which divides the Gauss points of the embedded trimmed element into two parts: (1) Gauss points in the solid domain, and (2) Gauss points in the fictitious domain. Only the Gauss points in the solid domain will be used in the trimmed element integration. However, the Gauss point classification is very complicated, especially for 3D problems with complex trimming cases, since it is difficult to identify which Gauss point is inside/outside the solid domain if the trimming curves or trimming surfaces are irregular shapes.
When the level set reconstructed model is used, the above problem can be solved by calculating the LSF values at the Gauss points by using a simple interpolation as is the LSF value of the Gauss point whose , and and are the shape function (linear interpolation is used herein and the LSF value of the ith vertex of the regular embedded element, respectively. The above equation can be easily extended to 3D problems as The Gauss point is in the solid domain if the LSF value of the Gauss point is equal or larger than 0, otherwise, the Gauss point is out of the solid domain (i.e., in the fictitious domain). An illustration of the Gauss point classification is shown in Fig. 16. In this work, two levels of refinement (i.e., the number of levels is 3) and 3×3×3 Gauss quadrature points in each level are used.

Numerical examples
In this section, we use numerical examples, including three CAD models and a point cloud model, to investigate the performance of our IGA method in detail. All examples are run on a desktop: CPU is Intel Xeon E5-2620 2.1GHz, RAM is 64GB, OS is Windows 7 Professional, and the compiler is MATLAB 2014b.

Cylinder ν
The first example is a bottom-fixed cylinder subjected to the pressure caused by a heavy box as shown in Fig.  17(a), which can be converted into an analysis model as Fig. 17(b). The boundary polygon model of the cylinder is shown in Fig. 17(c), which can be embedded in a regular domain corresponding to a complete tensorproduct NURBS domain. In this model, Young's modulus, E, is 10 5 , the Poisson's ratio, , is 0.3, the bottom surface of the cylinder is fixed and the equivalent pressure on the square domain of the top surface is 100.
The cylinder is embedded into a regular domain with 7×7×14 quadratic 3D NURBS elements. The control points of the top surface are given in Fig. 18, and the pressure loading region is just a part of the surface patch. Imposing the accurate load at control points is not easy since the NURBS control points are not endpoint interpolatory. The load F to the solid control points that influences the loading region in Fig. 18 can be obtained by the following equation τ Ω where F is the force vector that should be added at the control points of the loading region, N is the NURBS shape function vector, which is related to the integration point, and is the traction of the loading region . More accurate loading approaches require more integration points, which can be obtained by the local refinement of the elements. The displacement and the stress along the pressure direction (i.e., Z direction) are given in Fig. 19. Two finite element reference results from ANSYS (Canonsburg, Pennsylvania, USA) using 1873 (Figs. 19(a) and 19(d)) and 72027 quadratic tetrahedral elements ( (Figs. 19(b) and 19(e)) are used to compare with the results of the proposed method ( (Figs. 19(c) and 19(f)). The displacement results are consistent for all the three models, but the large displacement region around the loading region is larger in the IGA result. The reason for this is that the forces are imposed at the control points   (solid control points in Fig. 18) which influence a larger region than the loading region of FEA. Comparing to the FEA with 1873 elements whose element size is similar to that of IGA, the stress distribution of IGA is closer to the FEA with 72027 elements which is regarded as the reference standard, especially for the stress concentration around the edge of the bottom surface. The above demonstrates the validity and accuracy of the proposed IGA method.

Automotive part
An automotive part subjected to tension loading is shown in Fig. 20, where the bottom is fixed; and the tension pressure is 10 4 , and Young's modulus and the Poisson's ratio are 10 9 and 0.3, respectively. The automotive part is embedded into a regular domain 16 × 32 × 8 = 4096 with quadratic 3D NURBS elements. Figs. 21(a) and 21(b) show the total deformation results of the FEA of ANSYS using 40143 quadratic tetrahedral elements and the proposed IGA, respectively. The two results are almost the same, that is, both results demonstrate the high accuracy of the proposed method for real engineering parts. Notably, the boundary triangular facets can be directly obtained from the CAD system, which usually uses triangular facets for visualization [57]. Thus, an additional meshing tool is not necessary.

Point cloud of a gear shaft
Different from the abovementioned examples whose reconstruction models are constructed based on boundary polygons, this gear shaft example will show how to directly implement the IGA based on point cloud data. The flowchart of the IGA in this paper is given in Fig. 22 which consists of the following main steps: (1) The 3D scanner JTscan-DS (Jeatech, Guangzhou, China) is used to scan the physical model of the gear shaft. The 360-degree auto scan mode is used to scan the model, and the oriented point cloud can be directly obtained by using the JTscan-DS software. Therefore, the reconstruction model can be constructed by applying the approach proposed in Section 3. (2) A total of 140959 points can be observed in the point cloud, and the point cloud is embedded into a regular domain with quadratic 3D NURBS elements. The analysis model is constructed by fixing the bottom and adding forces at the 24 control points (the force for each control point is 100) of a gear tooth as shown in Fig. 22(b), and Young's modulus and the Poisson's ratio are 105 and 0.3, respectively.
(3) After the IGA, the deformation of the control points is obtained, and the deformation of each point of the point cloud can be also obtained by the NURBS interpolation, which is shown in Fig. 22(c). The maximum deformation occurs at the force loading position while the minimum deformation occurs on the fixed bottom face, which is reasonable for the gear shaft model subjected to the boundary condition as Fig. 22(c). This example demonstrates that our method can be directly used to analyze point clouds.    In the future, a better method to impose Dirichlet boundary conditions will be explored to solve the lowfidelity loading problem, and a more accurate and efficient quadrature method (e.g., the optimal quadrature approach [31,32]) will be tested to replace the adaptive quadrature algorithm. Furthermore, the proposed method will be expanded to structure optimization [62,63] and other engineering problems. The images or other third-party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

Nomenclature
To view a copy of this license, visit http://creativecommons.org/licenses/ by/4.0/.