A High-Order Isogeometric Collocation Method with Adaptive Refinement and Multi-patch Geometric Modeling

Abstract

Purpose

Isogeometric analysis (IGA) seamlessly integrates computer-aided design (CAD) with finite element analysis (FEA), streamlining the transition from geometric modeling to structural analysis.

Methods

This paper introduces a high-order isogeometric collocation method (IGA-C) specifically designed for analyzing complex multi-patch geometric structures. The method effectively resolves two primary challenges: the optimal selection of collocation points in high-order elements and ensuring stability in computations under complex geometric boundary conditions.

Results

Our contributions are threefold: first, we develop high-order basis function elements featuring local adaptive refinement tailored for IGA-C. Second, we investigate the optimal placement of collocation points across elements of varying orders, emphasizing their strategic distribution within non-uniform grids. Third, we present the Gaussian collocation method (IGA-GC) that employs advanced PHT-spline elements to facilitate the analysis of intricate multipatch structures.

Conclusions

Additionally, this work expands the algorithm to support elements of arbitrarily high order within the IGA-GC framework. Our numerical experiments validate that the proposed method not only achieves optimal convergence rates but also adeptly navigates the complexities associated with multi-patch geometric structures.

Introduction

Many structural problems in engineering ultimately require solving various partial differential equations (PDEs). Although these PDEs may appear similar, the specific challenges involved can vary significantly, including differences in geometry, material properties, loads, and boundary conditions. If an analytical solution were available, one could directly derive high-order derivatives from it. However, most engineering structures feature irregular geometries and complex boundary conditions, and often, no analytical solution exists. Consequently, these problems must be addressed using numerical methods. Finite element methods (FEM) [1], finite difference methods (FDM) [2], and other computational techniques play crucial roles in this context, offering solutions where traditional analytical methods fall short.

The traditional FEM process encompasses seven critical steps to effectively address engineering problems. Initially, the geometric configuration of the structure is defined using computer-aided design (CAD) tools. Subsequently, a mesh is generated based on the CAD geometry to accurately represent the structure. Analytical calculations are facilitated by defining Lagrangian polynomial functions on this newly generated mesh. Following this, essential matrices, including the stiffness and mass matrices, are assembled, and boundary conditions are applied. The resulting system of linear equations is then solved to determine the unknown quantities. An error analysis is conducted afterward to evaluate the accuracy of the solution. Based on this error analysis, the mesh is refined, and the solution process is iteratively repeated to enhance the accuracy of the results. This methodical sequence ensures comprehensive analysis and refinement, as depicted in Fig. 1a, which illustrates the main algorithms utilized at each stage.

[See PDF for image]

Fig. 1

Comparison of computational frameworks in traditional FEM and IGA. a Illustrates the traditional FEM pipeline, detailing the typical workflow from geometric modeling to mesh generation and solution computation. b Depicts the IGA mapping process, highlighting the seamless integration of geometric design and analysis by using a common set of basis functions for both tasks, thus streamlining the analysis pipeline

Research from Sandia National Laboratories in the United States indicates a significant inefficiency in industrial product design and structural analysis processes. According to their findings, as depicted in Fig. 1a, 60% of the total time is consumed by generating geometric models suitable for analytical calculations, with an additional 20% dedicated to mesh reconstruction. Collectively known as geometric configuration, these stages account for 80% of the workflow [3]. This distribution of effort is sub-optimal for most engineers, who would prefer to allocate more time to the actual structural analysis to achieve more functional designs. Furthermore, the current process involves building a CAD model that, once created, must be reconstructed with a mesh. This step not only consumes 60% of the total project time as a one-off task but also introduces geometric errors. Moreover, once the geometric configuration is set, making subsequent improvements to the structure requires reverting to the initial CAD design phase, thus highlighting the irreversible nature of this phase of the workflow.

In 2005, Professor T. J. R. Hughes, a distinguished academician of the United States Academy of Sciences, introduced the concept of isogeometric analysis (IGA) [4]. The core innovation of IGA is its use of spline functions, specifically non-uniform rational B-splines (NURBS) [5], to simultaneously perform geometric design and structural analysis, as illustrated in Fig. 1b. This approach establishes a seamless transition from design to analysis without necessitating geometric reconstruction, thereby avoiding the introduction of geometric errors and rendering the process reversible. Such improvements significantly reduce the costs associated with structural design. IGA allows for the construction of complex and accurate geometric structures with fewer elements and provides high-order continuity, which is particularly challenging to achieve with traditional FEM. Moreover, the stability of the method improves with the increased order of the elements [3]. Over the years, IGA has evolved and found extensive applications in various fields including the analysis of plate and shell structures [6, 7], structural optimization [8, 9], biomedical applications [10, 11], and acoustics [12], demonstrating substantial potential across these domains.

Currently, IGA methods can be categorized into two primary types: those based on the traditional isogeometric Galerkin method and those utilizing the isogeometric collocation method. The latter, particularly the traditional isogeometric collocation method (IGA-C) [13], offers significant advantages in terms of computational efficiency and cost, making it highly valuable for engineering applications. However, the research on IGA-C is still in its nascent stages, with several theoretical and algorithmic challenges that need addressing. Key issues include: (1) the development of locally refined functions and the determination of collocation points in non-uniform elements; (2) the selection of collocation points in high-order elements; and (3) addressing the continuity across multiple pieces of geometric structures and managing the boundary conditions of complex geometric structures. These challenges require further exploration to enhance the efficacy and applicability of IGA-C in practical engineering scenarios.

The selection of collocation points is critical in the collocation method, as it directly influences the accuracy and stability of computations. Traditionally, IGA-C employs Greville abscissae [14] as collocation points, which are calculated from the knot vector. This approach ensures that the number of basis functions matches the number of collocation points, resulting in a square coefficient matrix. However, one drawback is the relatively lower computational accuracy, particularly the error convergence rate of odd-order elements, which falls below theoretical expectations [13]. To address these shortcomings, Cosmin et al. proposed utilizing super-convergence points as collocation points, which are known to reduce collocation error [15]. Nevertheless, this approach leads to an over-determined coefficient matrix since the number of super-convergence points surpasses the number of basis functions, resulting in only an approximate solution. Additionally, it fails to accommodate non-uniform mesh refinement. Addressing non-uniform refinement necessitates research into spline elements that are amenable to collocation methods and capable of local refinement. Jia et al. introduced a Gaussian collocation method that employs Gaussian points—proven by Carl de Boor and Blair Schwartz to be super-convergence in one-dimensional B-spline elements with $C^1$ continuity [16]. These elements also exhibit $C^1$ continuity, and initial two-dimensional numerical experiments demonstrated optimal convergence rates for non-uniform elements [17]. However, this method currently only applies to third-order PHT spline elements and has not been extended to higher-order elements. Future work should focus on exploring collocation points for high-order PHT spline elements to advance the field further.

Initially, traditional IGA primarily employed NURBS, which only allowed for uniform mesh refinement. This limitation was later addressed by extending the algorithm to include other types of splines, such as T-splines [18], LRB-splines [19], THB-splines [20], and PHT-splines [21], all of which support adaptive refinement. PHT-splines, introduced by Deng et al. [21], are particularly notable for features like the convex hull property, affine invariance, and local support. The development of PHT-spline functions involves a three-step process: initially, mesh refinement is conducted through cross-insertion, followed by the redefinition of element basis functions via Bézier extraction [22, 23]. Lastly, to maintain the linear independence of the basis function space, some functions near T-nodes are truncated, a process that can lead to decay issues in traditional PHT-splines. This decay, characterized by non-convex domains and loss of the bell-shaped structure with excessive refinement, deteriorates the approximation quality of the basis functions as the element order increases [24]. To overcome these challenges, Deng et al. proposed an improved PHT-spline function that eliminates the need for truncation during refinement [24]. This improved version involves selecting function nodes (both interior and boundary), defining the support for each node as the smallest rectangular region encompassing it, and using local knot vectors to establish the basis functions for each node interval. These enhancements make improved PHT-splines more suitable for IGA-C, especially for achieving high-order continuity and handling complex geometric features. A key aspect of ongoing research is applying these improved PHT-splines to IGA-C with local refinement capabilities, aimed at addressing challenges presented by complex geometric structures. Previous studies primarily focused on cubic PHT spline elements, leaving high-order elements less explored. This paper aims to define and examine high-order PHT spline elements theoretically and explore their potential applications in IGA-C.

Conforming patches are essential in FEM and IGA applications, as they share the same nodes and edges along their boundaries, ensuring continuity and smoothness. These patches can be created using various techniques such as subdivision [25], refinement [26], and optimization [3]. Conforming patches are particularly advantageous for multi-patch geometric modeling and structural analysis due to their ability to maintain element continuity seamlessly. Nonetheless, significant research also focuses on non-conforming patches, which do not share common nodes and edges and thus present challenges in continuity across patch boundaries. To address this, several methods are employed to ensure continuity, including penalty methods [27], Lagrange multiplier methods [28], and Nitsche’s methods [29]. Each method offers different benefits and is chosen based on the specific requirements of structural analysis. Some important relations and conclusions related to this work are available in [30, 31, 32, 33, 34, 35, 36–37]. Sayyad et al. [38] studied the effect of transverse out of plane normal strain on the bending results of a multi-layered arch subjected to a point load. Mahajan and Sharma [39] employed a novel higher order shear and normal deformation theory for studying static and dynamic responses of a multi-layered composite arch based on exponentially varying the out of plane shear strain. Ren and Lin [40] extended a review work on the isogeometric collocation methods with some applications to mechanical and mathematical simulations. Mahajan and Sharma [41] developed an analytical method for static analysis of multi-layered composite arches subjected to mechanical loadings. The cantilever beams can be used in the technical equipment and as an element of electrical and mechanical systems and structures. Jiang et al. [42] studied effect of woven fibers on the vibration responses of a composite cantilever beam. Pathak et al. [43] studied impact of a crack on the vibrational responses of cantilever beam using the surface methodology. Impact of nonlinear strain was studied on the large amplitude nonlinear vibration analysis of composite cantilever beam in moving mode [44]. Kurt and Orhan [45] developed energy harvesting applications of the cantilever beams.

This work is dedicated to advancing the field of IGA by focusing on the development of a high-order IGA-C method tailored for the intricate analysis and computational needs of complex multi-patch geometric structures. Our primary research objectives are articulated across three pivotal areas: enhancing the adaptive IGA-C methodology using advanced PHT spline functions to extend the capability of the algorithm for high-order continuity calculations, which improves accuracy and computational efficiency; investigating the optimization of collocation points within high-order elements, particularly non-uniform elements, to identify strategies that minimize computational errors and maximize stability; and applying an advanced IGA-GC method based on PHT spline elements to demonstrate the practical viability and robustness of our theoretical developments in solving real-world multi-patch complex geometry analysis problems. Additionally, this article presents a significant enhancement over our previous work [17], by broadening the adaptive capabilities of the IGA-C to incorporate high-order elements more effectively [46, 47–48].

The structure of this paper is organized into six sections. Initially, we provide an overview of the fundamental principles of IGA. Following this, we delve into various techniques used in constructing PHT splines. Subsequently, the isogeometric collocation method is discussed. Afterwards, we elucidate the core concepts of the Gaussian collocation method. Subsequent to this, we conduct several numerical experiments to validate our methods. The paper concludes with a summary of our findings and directions for future research.

Isogeometric Analysis

IGA represents an advanced computational approach that utilizes high-order continuous finite element calculations. This method employs a basis function to simultaneously perform geometric modeling and structural analysis, thereby establishing a seamless transition from geometric modeling to structural analysis without the need for intervening geometric reconstruction. As a result, it eliminates the introduction of geometric errors and the necessity for mesh reconstruction. Traditionally, IGA is based on NURBS, a fundamental geometric modeling function in CAD. NURBS facilitate the computation of arbitrary high-order continuity. Unlike traditional finite element methods, IGA defines both parameter space and physical space, with a one-to-one correspondence between elements in parameter space and those in physical space (illustrated in Fig. 1b). The transformation of parameters from parameter space to physical space elements is efficiently achieved using a specific transformation Eq. 1, streamlining the process and enhancing the accuracy of the analysis.

1

$$\mathit{F}(\xi )=\sum _{i=1}^n{N}_i^p(\xi )·{\mathit{C}}_i$$where $N_i^p$ represents the $i$th B-spline basis function with polynomial degree $p$. ${\mathit{C}}_{\mathit{i}}$ indicates the control points corresponding to the basis functions, which are defined in physical space. When employing B-spline functions for geometric modeling, having precise geometric parameter information is essential because B-splines typically do not satisfy the interpolation property. Therefore, to model sharp features accurately, one can strategically reduce the continuity at specific nodes. This approach allows for the precise delineation of detailed geometric features. Figure 2a displays a generated cubic B-spline curve, where control points are indicated by blue solid dots, and the resulting curve is depicted as a blue line. This curve utilizes an open knot vector, characterized by the repeated occurrence of the first and last nodes in the node vector, which impacts the curve's continuity. Specifically, each repetition of a node in the knot vector decreases the continuity at that point; for instance, node 5 is repeated twice, resulting in the curve being $C^0$ continuous at this node. The corresponding basis functions are illustrated in Fig. 2b. Figure 2c, d provide 2D examples, demonstrating how B-splines can be utilized to model complex surfaces. Figure 2c shows a control net, highlighting the influence of control points on the shape of the surface. Figure 2d depicts the corresponding 2D basis functions that define the surface shape, underscoring the flexibility and precision of B-spline modeling in capturing intricate geometric details.

[See PDF for image]

Fig. 2

Overview of B-spline representations and basis functions. This set of four figures demonstrates the fundamental elements of B-spline geometry. a A 1D B-spline curve showcasing its smooth and continuous nature. b The 1D B-spline basis functions that define the shape of the curve and properties. c A 2D control net illustrating the structural framework for surface modeling. d The corresponding 2D basis functions that influence the surface shape and continuity

Figure 3 highlights the distinctions between traditional FEM and IGA in terms of geometric modeling. In traditional FEM, as depicted in Fig. 3a, the geometric model from CAD is used for mesh reconstruction, shown in Fig. 3b. However, this figure also reveals a noticeable geometric discrepancy between the FEM mesh boundary and the boundary of the actual model. To minimize this geometric error, FEM often requires defining a large number of elements at feature boundaries, which not only complicates the calculations but also consumes significant storage space. In contrast, IGA employs a basis function for both geometric modeling and structural analysis, eliminating the need for intermediate mesh reconstruction. Figure 3d illustrates a problem domain constructed using four spline elements, with Fig. 3c displaying the corresponding control points. This streamlined approach in IGA enhances accuracy and efficiency by maintaining the integrity of the geometric model throughout the analysis process.

[See PDF for image]

Fig. 3

Workflow from CAD model to physical domain representation. a Displays the initial CAD model, representing the original geometric design. b Shows the FEM mesh generation, which discretizes the CAD model for analysis. c presents the control net, used in IGA to define the geometry of the model. d Illustrates the resulting physical domain, as defined by the control net, demonstrating the final shape used for analysis

Techniques in PHT Spline Construction

PHT spline functions, polynomial splines defined on hierarchical T grids, were initially utilized for geometric modeling [21]. These functions possess several advantageous properties over traditional B-spline functions, including non-negativity, normalization, and a local definition domain—essential attributes for basis functions in FEM. Additionally, PHT splines offer local subdivision capabilities, inheriting all the favorable characteristics of NURBS [30, 31]. As a result, these functions have been effectively applied in adaptive geometric analysis algorithms [26, 32] and have been used to address engineering challenges in areas such as vibration [33, 34] and fracture mechanics [35]. The construction of traditional PHT spline functions involves three main steps: firstly, selecting the unit for subdivision and dividing the quadrilateral unit into four sub-units through cross-insertion. Secondly, the basis function for each node is expressed using Bézier coordinates. Thirdly, a process known as truncation is applied, where Bézier coordinates associated with newly generated nodes are zeroed out to redefine the domain of the original basis function. Subsequently, the B-spline basis function corresponding to the new node is defined. However, traditional PHT splines exhibit two major limitations: each local subdivision leads to truncation around the T node, potentially causing basis function decay—a consequence of severe function deformation, as noted in previous studies [24]. Moreover, there has been no effective algorithm to extend traditional PHT splines to high-order elements comprehensively, with basis functions in such elements being particularly susceptible to decay. To address these issues, Deng et al. introduced an improved PHT spline function [24], aiming to enhance performance and applicability in complex applications.

The improved PHT spline function retains all the advantageous properties of the traditional PHT spline function while successfully avoiding the decay phenomenon associated with the latter. The construction of the improved PHT spline unit involves three distinct steps: first, selecting the basis function nodes, which are categorized into nodes within the problem domain and nodes on the problem boundary. Second, the minimum rectangular scope for each basis function node is defined. Finally, the B-spline basis function for each node is established by a tensor product based on the local node vector. This method ensures that the scope of each basis function node remains a rectangular area, resulting in all associated basis functions maintaining a bell-shaped curve, thereby preventing the decay phenomenon caused by severe deformation. Moreover, the use of tensor products to define basis functions facilitates easy extension to arbitrarily high-order units. These significant improvements make the enhanced PHT spline function a promising candidate for integration with the IGA algorithm, offering substantial benefits for research and application in this area.

To construct PHT spline basis functions, we need to first determine the type and domain of the basis function nodes. The basis function nodes can be divided into interior nodes and boundary nodes, and the boundary nodes can be further divided into left and right boundary nodes and upper and lower boundary nodes. We can find the rectangular domain range of each intermediate basis function node according to the coordinates and marked elements of the basis function nodes, and define the local node vectors. For example, for a basis function node ${(u}_i,v_j)$, we can define the following local node vectors:

2

$$LocalKnotU=\left[u_i^{\mathit{left}},,,u_i^{\mathit{left}},,,u_i,,,u_i,,,u_i^{\mathit{right}},,,u_i^{\mathit{right}}\right],$$

3

$$LocalKnotV=\left[v_j^{\mathit{down}},,,v_j^{\mathit{down}},,,v_j,,,v_j,,,v_j^{\mathit{up}},,,v_j^{\mathit{up}}\right].$$

If $LocalKnotU(1)=0$, $LocalKnotU(2)=0$, $LocalKnotU(3)=0$, $LocalKnotU(4)=0$, then the node is a left boundary node. Similarly, criteria for right, upper, and lower boundary nodes are applied. Following this classification, it is essential to utilize the Bézier extraction operator to represent the spline basis functions. This operator enables the expression of spline basis functions as linear combinations of Bernstein polynomial functions. It is necessary to calculate the coefficients, also known as Bézier extraction operators, that precede the Bernstein polynomial functions in these linear combinations. Given that PHT spline elements exhibit only $C^1$ continuity between them, each node is associated with only four basis functions. Consequently, it is crucial to identify which of these basis functions are boundary basis functions and to further categorize the type of boundary each represents. This step is critical as it simplifies the process of imposing boundary conditions later in the analysis, ensuring that the mathematical and computational framework is robust and effective for subsequent applications.

To efficiently manage spline functions within the framework of Bézier extraction, an exhaustive traversal of all basis function nodes is conducted. This traversal involves identifying all basis functions associated with each node and subsequently pinpointing all elements linked to that node. The Bézier extraction operator is then employed to represent these basis functions within their respective elements. It is important to note that the basis functions represented might only constitute portions of the original spline function. This representation process is extended to adjacent elements as well, where segments of basis functions are amalgamated to reconstruct the complete original spline function, which is a central concept of the Bézier extraction technique. As part of the refining process, which is hierarchical, a comprehensive traversal of all elements is necessary to ensure the retrieval of the lowest level elements without any overlap. Each element is encapsulated within a structured variable that holds various property variables, facilitating the assignment of relevant data. This structure allows for each function Bernstein basis function coefficients and function number within the element to be assigned, alongside the attributes of each function (identifying whether it is a boundary function or an internal function). By systematically assigning all pertinent information to each element, this method lays a robust foundation for constructing both the local and global matrices of the linear solving system. Additionally, given the importance of the order of basis functions when imposing boundary conditions, it proves beneficial to systematically number the functions along the boundary in sequence. This ordering facilitates the precise imposition of boundary conditions, enhancing the overall integrity and effectiveness of the analysis.

IGA Collocation Method

Isogeometric analysis methods are primarily classified into two categories: the traditional IGA-G and the traditional IGA-C. Unlike IGA-G, which addresses the weak form of PDEs, IGA-C tackles the strong form directly, as shown in the general expression of the boundary value problem:

4

$$\mathcal{L}\left(u\left(x\right)\right)=f\left(x\right)\text{in}\mathrm{\Omega }$$

5

$$\mathrm{\Gamma }\left(u\left(x\right)\right)=g\left(x\right)\text{on}\partial \mathrm{\Omega }$$where $\mathrm{\Omega }$ and $\partial \mathrm{\Omega }$ represent the problem domain and its boundaries, respectively, $u(x)$ is the unknown solutions, $\mathcal{L}$ and $\mathrm{\Gamma }$ are differential operators, and $f(x)$ and $g(x)$ are known functions. While many engineering problems are reducible to PDEs, obtaining analytical solutions is often infeasible, necessitating numerical approximation methods such as:

6

$$u_n=\sum N_i{a}_i,$$where $N_i$ denotes the basis function and $a_i$ the coefficient to be determined. The IGA-G method involves approximating the solutions into Eqs. 4 and 5, multiplying by a test function, integrating, applying the Gaussian divergence theorem, and integrating by parts to derive the weak form of the PDE, subsequently solving the linear system to approximate $u$. In contrast, IGA-C simplifies this process by directly substituting the approximate solution into the strong forms of the PDEs and selecting collocation points within $\mathrm{\Omega }$ and on $\partial \mathrm{\Omega }$ to formulate and solve the system of linear equations, where the number of rows in the coefficient matrix equals the number of collocation points, matching the number of basis functions, thus ensuring a square matrix. The solution of IGA-C satisfies the following relationship.

7

$$\mathcal{L}\left(u_h\left(x_i^{\mathit{int}}\right)\right)=f\left(x_i^{\mathit{int}}\right)\text{in}\mathrm{\Omega }i=1\dots N^{\mathit{int}},$$

8

$$\mathrm{\Gamma }\left(u_h\left(x_j^{\mathit{bnd}}\right)\right)=g\left(x_j^{\mathit{bnd}}\right)\text{on}j=1\dots N^{\mathit{bnd}}.$$

where $\{x_i^{\mathit{int}}\},i=1,\cdots ,N^{\mathit{int}}$ and $\{x_j^{\mathit{bnd}}\}=1,\cdots ,N^{\mathit{bnd}}$ are the set of the collocation points inside the problem domain and on the problem boundary. Traditional IGA-C typically employs "Greville abscissae" [14] as collocation points, calculated from the given knot vector.

9

$$\overline{{\xi }_i}=\frac{{\xi }_{i+1}+{\xi }_{i+2}+\cdots +{\xi }_{i+p}}{p}$$where $\{{\xi }_1,{\xi }_2,\cdots ,{\xi }_{n+p+1}\}$ is given knot vector. The primary advantage of using Greville abscissae as collocation points in IGA-C is that their quantity matches the number of basis functions. This correspondence ensures the resulting system of collocation linear equations forms a square matrix, which theoretically enables precise solutions. However, empirical findings indicate that the computational accuracy achieved with the Greville abscissae method is suboptimal. Specifically, the error convergence rate for odd-order elements often falls below anticipated theoretical values [13]. Consequently, the determination of optimal collocation points for IGA-C remains a significant and unresolved area of research.

Gaussian Collocation Method

Compared to traditional IGA-C, the Gaussian collocation method (IGA-GC) exhibits several distinct advantages. Firstly, it utilizes improved PHT spline functions as the basis for adaptive, non-uniform mesh refinement, guided by criteria based on the recovered solution and its error [26]. Secondly, Gaussian points are employed as collocation points. This selection is strategic to match the number of basis functions with the number of collocation points, ensuring the coefficient matrix of the linear equation system remains square. In areas surrounding T nodes, excess collocation points are eliminated, and the remaining points are repositioned to optimize placement. Thirdly, this method supports extension to calculations involving high-order elements, achieving optimal convergence rates. Fourthly, IGA-GC is particularly effective in addressing complex geometric problems, such as structures with holes and multiple re-entrant corner points. These complexities often require the definition of multiple spline layers stitched together, where continuity not only in the geometric representation but also in the basis functions across layers and the alignment of control points are critical.

Refinement and Basis Function Organization in High-Order PHT-Splines

The definition of basis functions for cubic PHT spline elements is straightforward due to their $C^1$ continuity and the presence of 16 basis functions per element. These functions are categorized into four groups, each corresponding to a basis function node, which as shown in Fig. 4b. However, defining basis functions for higher-order PHT spline elements is more complex due to the varying number and types of functions within each element. In addition to the node basis functions, there are edge basis functions associated with the element edges and center basis functions linked to the center of the element. To organize the basis functions in high-order elements, we partition them into nine regions as depicted in Fig. 4a. This partition includes four node regions, four edge regions, and one center region. Subsequent Fig. 4b–f illustrate the number of basis functions associated with each region for polynomial degrees ranging from 3 to 7. It is noteworthy that the number of basis functions at the corners remains constant at four across all polynomial degrees, maintaining $C^1$ continuity along the elements. Conversely, the number of basis functions associated with the edges and the center region increases as the polynomial degree rises, reflecting the growing complexity of structure of splines in high-order elements.

[See PDF for image]

Fig. 4

Distribution of basis functions across regions in high-order elements. a Depicts the partitioning of basis functions into nine regions within a high-order element. b–f Illustrate the number of basis functions associated with each region, with the polynomial degrees ranging from 3 to 7

The connection of elements and the update of basis vertices are critical considerations after the cross-insertion process, which modifies the connectivity of the elements and alters the support of the basis vertices. As depicted in Fig. 5, cross insertion subdivides a parent element into four child elements, introducing five new vertices. However, not all of these vertices automatically qualify as basis vertices. For example, T-junctions are excluded as basis vertices since they do not introduce new basis functions. In Fig. 5, the central vertex, marked by the red dots and with a valence of four, qualifies as a basis vertex due to its full connectivity. Conversely, the status of the four vertices on the edges of the refined element requires additional evaluation. A new vertex that lies on the edge of the problem boundary is considered a basis vertex, as shown by the green dot in Fig. 5. If the vertex is not on the problem boundary, its classification as a basis vertex depends on the connectivity with neighboring elements. For instance, the orange triangle in Fig. 5 represents a T-junction and is not considered a basis vertex due to its insufficient connectivity and lack of contribution to new basis functions. This selective approach ensures that only vertices contributing to the structural and functional integrity of the mesh are designated as basis vertices. Additionally, as elements are further refined, T-junctions may be eliminated, altering their status in the mesh. For example, a T-junction between elements 2 and 4 could be removed with further refining of an element to the east. In such cases, the node transitions from a T-junction to a new basis vertex, necessitating the definition of all required information for the vertex and updating all originally connected elements accordingly. This ongoing adjustment is vital for maintaining the accuracy and integrity of the mesh structure.

[See PDF for image]

Fig. 5

Analysis of vertex classification in element subdivision by cross insertion. a Illustrates the subdivision of a parent element into four child elements through cross insertion, introducing five new vertices. b Focuses on the evaluation of edge vertices, depicting a vertex on the boundary with a green dot, indicating its status as a basis vertex, and a T-junction vertex marked by an orange triangle, which is not considered a basis vertex due to its limited connectivity and absence of new basis functions

PHT elements are refined through cross insertion, which necessitates updating the connections between basis vertices and elements. The Bézier extraction method is an invaluable tool for computing spline functions, particularly effective for both cubic and higher-order elements. This method involves expressing the spline function as a linear combination of Bernstein basis functions and determining the corresponding coefficient matrix for each element. The process of cross insertion and the calculation of Bézier coefficients are systematically outlined in four detailed steps to ensure clarity and efficiency in implementation.

Step 1: Initially, it is essential to identify and label the basis vertices. The support for each basis vertex, typically denotes as $U_{\mathit{left}}$,$U_{\mathit{right}}$,$V_{\mathit{down}}$,$V_{\mathit{up}}$, is determined based on its location relative to other elements, as illustrated in Fig. 6a. This involves starting from the center of the basis vertex and extending outwards in the left, right, down, and up directions to define the minimum boundaries of the rectangular region that will constitute the support of the basis vertex. If a basis vertex is positioned on the problem boundary, its support is adjusted to match the coordinates of the vertex, ensuring that the boundary representation is precise, as depicted in Fig. 6b.

[See PDF for image]

Fig. 6

Determining the support for basis vertices. a Demonstrates the method for defining the support of a basis vertex by extending boundaries outward in the left, right, down, and up directions from the center of the vertex to determine the rectangular support region. b Illustrates the adjustment of support for a basis vertex located on the problem boundary

Step 2: Each element is defined by four vertices. It is crucial to establish a relationship between these elements and the basis vertices. This involves identifying all elements that are supported by the basis vertices.

Step 3: To define edge basis vertices, it is first necessary to determine whether an edge is vertical or horizontal. Once classified, the coordinates of the edge basis vertices can be calculated. The support for these edge basis vertices is then defined, ensuring it does not extend beyond the boundaries of the element itself; hence, the corner points of the element are used to define the limits of the edge basis support. Finally, as elements undergo refinement, the support for the original basis must be updated, which typically results in the support dimensions being reduced by half.

Step 4: For cubic elements, Bézier extraction involves identifying node basis functions. Typically, the knot vector for these elements features four repeated values at each endpoint and two repeated values at each interior node, which facilitates the definition of these functions. For elements of the fourth order or higher, the Bézier extraction process not only includes node basis functions but also extends to boundary and center basis functions within each element. For a corner basis vertex, the local knot vector configuration is:

10

$$LocalKnotU=\left[u_i^{\mathit{left}},,,u_i^{\mathit{left}},,,u_i^{\mathit{left}},,,u_i,,,u_i,,,u_i,,,u_i^{\mathit{right}},,,u_i^{\mathit{right}},,,u_i^{\mathit{right}}\right],$$

11

$$LocalKnotV=\left[v_j^{\mathit{down}},,,v_j^{\mathit{down}},,,v_j^{\mathit{down}},,,v_j,,,v_j,,,v_j,,,v_j^{\mathit{up}},,,v_j^{\mathit{up}},,,v_j^{\mathit{up}}\right].$$

In fourth-order elements, the knot vector for the center basis function contains three repeated values at each end, as shown:

12

$$LocalKnotU=\left[u_i^{\mathit{left}},,,u_i^{\mathit{left}},,,u_i^{\mathit{left}},,,u_i^{\mathit{right}},,,u_i^{\mathit{right}},,,u_i^{\mathit{right}}\right],$$

13

$$LocalKnotV=\left[v_j^{\mathit{down}},,,v_j^{\mathit{down}},,,v_j^{\mathit{down}},,,v_j^{\mathit{up}},,,v_j^{\mathit{up}},,,v_j^{\mathit{up}}\right].$$

Edge basis vertices are differentiated into vertical and horizontal types. For horizontal type, the local knot vector is:

14

$$LocalKnotU=\left[u_i^{\mathit{left}},,,u_i^{\mathit{left}},,,u_i^{\mathit{left}},,,u_i^{\mathit{right}},,,u_i^{\mathit{right}},,,u_i^{\mathit{right}}\right],$$

15

$$LocalKnotV=\left[v_j^{\mathit{down}},,,v_j^{\mathit{down}},,,v_j^{\mathit{down}},,,v_j,,,v_j,,,v_j,,,v_j^{\mathit{up}},,,v_j^{\mathit{up}},,,v_j^{\mathit{up}}\right].$$

For the vertical type, the local knot vector configuration is:

16

$$LocalKnotU=\left[u_i^{\mathit{left}},,,u_i^{\mathit{left}},,,u_i^{\mathit{left}},,,u_i,,,u_i,,,u_i,,,u_i^{\mathit{right}},,,u_i^{\mathit{right}},,,u_i^{\mathit{right}}\right],$$

17

$$LocalKnotV=\left[v_j^{\mathit{down}},,,v_j^{\mathit{down}},,,v_j^{\mathit{down}},,,v_j^{\mathit{up}},,,v_j^{\mathit{up}},,,v_j^{\mathit{up}}\right].$$

Applying de Casteljaus algorithm within these configurations allows for the precise calculation of Bézier coordinates in each element, ensuring accurate and efficient spline function generation.

Gaussian Collocation Points: Elimination and Relocation

The objective of the discussed method is to minimize the residual between analytical and approximate solutions, where the latter is derived by substituting collocation points into the strong form of PDEs and solving the problems directly, bypassing variational methods. The placement of collocation points is crucial for the effectiveness of collocation methods, yet the optimization of these points remains an unresolved challenge. The dimensionality of the coefficient matrices is crucial; the number of rows is determined by the total number of collocation points, while the number of columns corresponds to the total number of basis functions. To ensure determinacy of the coefficient matrices, it is essential that the number of collocation points matches the number of basis functions.

PHT-splines are employed as basis functions, and the corresponding collocation points are typically Gaussian points, known for their super-convergence properties. However, challenges arise with non-uniform refinement: maintaining the same number of collocation points in each element, including those around T-junctions, leads to an excess of collocation points relative to basis functions, resulting in over-determined coefficient matrices that require iterative solutions.

This paper explores strategies for reducing and rearranging collocation points around T-junctions in adaptively refined elements. Due to the absence of new basis functions at T-junctions following cross insertion, collocation points adjacent to T-junctions can be eliminated to decrease the total number of collocation points. We demonstrate this process using a fourth-order element example: typically, nine Gaussian collocation points per element ensure equality between the total number of basis functions and collocation points under uniform refinement. However, adaptive refinement introduces complexities such as T-junction nodes, depicted in Fig. 7a. In elements affected by T-junctions (colored green and yellow in the figure), the nearest collocation points to the T-junctions are removed first, as shown in Fig. 7b. If the number of collocation points still exceeds that of the basis functions, further reductions are necessary, as depicted in Fig. 7c where four additional collocation points are deleted. Finally, to balance the equation, two extra collocation points are added along the common middle edge of elements adjacent to T-junction nodes, as illustrated in Fig. 7d. This process of elimination and relocation ensures that the number of collocation points precisely matches the number of total basis functions in fourth-order elements.

[See PDF for image]

Fig. 7

Managing collocation points in adaptive refinement near T-junctions. a Identifies T-junction nodes, highlighted by the green and yellow colored elements. b Shows the initial removal of the nearest collocation points to the T-junctions to avoid over-determination. c Illustrates further reduction by deleting four additional collocation points when the number still exceeds that of the basis functions. d Depicts the addition of two extra collocation points along the common middle edge of elements adjacent to T-junction nodes to maintain balance in the equation

Currently, there are no definitive guidelines for the elimination and relocation of collocation points around T-junction areas. Unlike IGA-G, where quadrature integration points evaluate the basis functions, IGA-C requires that the basis functions be evaluated directly at the collocation points. Therefore, in uniform refinements, collocation points should be uniformly distributed within the element, whereas in adaptively refined areas, the collocation points tend to be concentrated. Additionally, since T-junctions do not introduce new basis functions, it is feasible to eliminate several collocation points near the T-junction until the number of collocation points matches the number of basis functions, as illustrated in Fig. 8. In some instances, such as when the polynomial degree is 4 or 6, as shown in Fig. 8b, d, a strategic relocation of certain points becomes necessary. A more precise mathematical analysis of collocation points within higher-order collocation elements requires further investigation to optimize their distribution and effectiveness.

[See PDF for image]

Fig. 8

Gaussian collocation points distributions in high-order elements. This figure consists of five sub-figures, each containing a pair of plots: the left plot in each pair demonstrates uniform refinement with collocation points, while the right plot showcases adaptive refinement with collocation points. a–e Progressively display refinements from cubic to seventh order, illustrating the differences in collocation point distribution as the polynomial order increases

Multi-patches Geometric Modeling

In geometric modeling, constructing multi-patch structures is essential for representing complex geometries. Each patch is typically defined independently, with information localized within its boundaries. However, when considering the entire structure, it is crucial to seamlessly connect these patches. This is achieved by ensuring that the nodes and control points along the boundaries of adjacent patches are identical. This method simplifies maintaining continuity at patch-to-patch boundaries, even when local refinement of elements occurs along these boundaries. We employ a strategy where the control points and element nodes are precisely aligned along the patch interfaces, and both maintain the same elements and control points along the adjoining lines. In our approach, nodes, control points, and basis functions are assigned a global index. This global indexing, such as ‘NodesGlobal’, uniquely identifies nodes within elements and is especially useful for delineating adjacent elements in multi-patch structures. Although each patch element maintains its own set of locally defined node variables, these must be uniquely characterized by global identifiers to facilitate the integration of all patches. Adjacent patches share common edges, and it is critical that the nodes on these shared edges are not only coincident but also consistently positioned to ensure smooth and continuous transitions between patches. During each refinement iteration, it is necessary to verify that the division of elements along the patch boundaries is consistent. If discrepancies are observed, cross-insertion refinement techniques are applied to the sparser side until a uniform element structure without T-nodes is achieved along the common edges. This process guarantees the continuity and smoothness required for high-quality geometric modeling of complex structures.

Recovery Error Estimation

In adaptive refinement processes, efficiently determining the necessity for refinement is paramount. Commonly, this involves comparing the approximate solution to the analytical solution; however, deriving an analytical solution for complex scenarios is often not feasible. This work introduces an alternative approach: using a recovery solution. By interpolating the approximate solution at points of super-convergence, a recovered solution can be derived, which generally demonstrates reduced errors compared to the original approximation. The discrepancy between these solutions, referred to as the 'recovery error', can be quantitatively assessed for each element. By evaluating this recovery error against a predetermined threshold, the methodology enables precise subdivision of elements where improvement is most needed. This approach builds on the principles established by Zienkiewicz and Zhu in their foundational work on patch recovery for finite element methods, which has proven effective in error estimation and mesh refinement [26, 27]. Recent studies have also applied the recovery estimation method to IGA solvers [7, 22], demonstrating promising results. Nevertheless, both from a theoretical and a numerical perspective, there remains significant scope for further research in this area to enhance the efficacy and applicability of recovery-based refinement strategies.

Recovery error estimation is predicated on the assumption that the solution to the problem is smooth and continuous. Despite this, practical observations reveal that this technique retains its efficacy even in scenarios where these conditions are unmet, such as with non-smooth or discontinuous solutions. This robustness renders recovery error estimation an invaluable tool for pinpointing regions within the model that necessitate finer meshing. Once the coefficients of the approximate solution are ascertained, the corresponding function values at these points can be substituted, capitalizing on super-convergence theory to minimize errors. This creates a recovered solution that acts almost as an analytical benchmark, setting an upper limit on the accuracy for the approximate solution and facilitating precise error quantification and the requisite mesh refinement, even in the absence of the exact solution.

Numerical Examples

In this section, we apply IGA-GC method to solve five distinct numerical examples, ranging from a single annulus model to a complex multi-patch honeycomb model comprising 90 patches. Among these, certain models, such as the cantilever beam, feature singularities including re-entrant corners. Additionally, other models like the L-shaped brackets and the perforated beam are characterized by multiple circular holes.

Annulus Model

The annulus model is commonly used for method validation, and due to its symmetric geometry, only a quarter of the structure is considered for analysis. Detailed geometric specifications are provided in Fig. 9, where the inner area is subjected to pressure, and both the y-coordinates of the bottom edge and the x-coordinates of the left end edge are fixed. Figure 10 compares the outcomes of uniform refinement with adaptive refinement using cubic Gaussian collocation elements, demonstrating that adaptive refinement more effectively identifies stress concentration areas with more elements. Figure 11 shows the use of higher-order elements, ranging from fourth to seventh order, and highlights that using high-order elements increases computational costs. Figure 12a–c depict three stages of mesh refinement, with the remained figures in Fig. 12 presenting the computed stress results in three directions. Since an analytical solution is available for this problem, it is possible to compute convergence rates for error estimation. Figure 13 presents error estimations for elements from the cubic to the seventh order. A comparison of the results from IGA-GC and IGA-G indicates that both methods achieve optimal convergence rates.

[See PDF for image]

Fig. 9

Annulus model

[See PDF for image]

Fig. 10

Progression of refinement levels in uniform and adaptive methods. a–d Illustrate the uniform refinement process from level 2 to level 5, showing a consistent increase in mesh density at each successive level. e–h Depict adaptive refinement from level 3 to level 6, highlighting how refinement is targeted to specific regions

[See PDF for image]

Fig. 11

Distribution of collocation points in adaptive refinement across polynomial orders. This figure includes four sub-figures, each representing the placement of collocation points in adaptive refinement for polynomial orders from the fourth to the seventh

[See PDF for image]

Fig. 12

Sequential refinement and stress distribution in multiple directions. The first row, a–c, illustrates a sequence of refinements from level 1 to level 3. The second row, d–f, displays the stress contour plots in the xx-direction. The third row, g–i, depicts the stress contours in the xy-direction. The fourth row, j–l, presents the stress contours in the yy-direction

[See PDF for image]

Fig. 13

Comparative analysis of convergence rates using IGA-GC and IGA-G across polynomial orders. The left column, a–e, displays the convergence rates achieved by IGA-GC method for polynomial orders from the cubic to the seventh. The right column, f–j, illustrates the convergence rates achieved by the isogeometric Galerkin (IGA-G) method across the same range of polynomial orders

Cantilever Beam

A cantilever beam is a common structural element with one end fixed, often referred to as the fixed support, while the other end is free, as depicted in Fig. 14. This structure extends outward and typically bears concentrated or distributed loads. Cantilever beams are primarily used in applications requiring overhanging structures without additional supports, such as house balconies and bridges. This design allows the beam to extend into space without intermediate supports, providing an economical and practical solution in architectural design.

[See PDF for image]

Fig. 14

Cantilever beam model

Re-entrant corners are the internal corners of a building that have an angle of less than 180 degrees. These corners are stress concentration points, making them vulnerable to cracking and structural damage. To mitigate these issues, re-entrant corners often require enhanced reinforcement or specific design modifications. In numerical simulations, re-entrant corners present as singularities due to the abrupt changes in direction or curvature at these corners, which can affect the continuity and smoothness of the solution or its derivatives. Given that IGA relies on structured meshes with tensor-product bases, it struggles to effectively address these singularities because mesh refinement in one area tends to affect the entire mesh. Consequently, adaptive refinement becomes a practical solution for managing such singularities. This method facilitates localized refinement around the re-entrant corner while maintaining a coarser mesh in other areas, as illustrated in Fig. 15a–c. This approach ensures more accurate solutions without incurring unnecessary computational expenses, as demonstrated in the computed stress plots in Fig. 15. The benefits of adaptive refinement are clearly evident through error estimation, as depicted in Fig. 16. Even in the absence of an analytical solution, it is possible to compute the recovery errors between the recovered solution and the approximate solution. Adaptive refinement enables the retrieval of optimal convergence rates, whereas uniform refinement suffers due to the singularity region around the re-entrant corners.

[See PDF for image]

Fig. 15

Sequential refinement and stress distribution in multiple directions. The first row, a–c, illustrates a sequence of refinements from level 1 to level 3. The second row, d–f, displays the stress contour plots in the xx-direction. The third row, g–i, depicts the stress contours in the xy-direction. The fourth row, j–l, presents the stress contours in the yy-direction

[See PDF for image]

Fig. 16

Comparative analysis of convergence rates using IGA-GC and IGA-G across polynomial orders. The left column, a and b displays the convergence rates achieved by the isogeometric Gaussian collocation (IGA-GC) method for polynomial orders with the cubic and the fourth. The right column, c and d, illustrates the convergence rates achieved by the isogeometric Galerkin (IGA-G) method across the same range of polynomial orders

L-shape Bracket

The third example features an L-shaped bracket, as illustrated in Fig. 17. One end of the structure is fixed, while the other end is subject to a tensile force. This bracket is characterized not only by its L-shape but also by the presence of four circular holes on the plate. Based on prior experience, it is anticipated that stress will primarily concentrate around the circular holes and near the fixed end, a hypothesis that has been validated by numerical experiments. Figure 18 depicts the mesh, which has been divided using third-order PHT spline elements. A total of 18 spline patches were utilized in this structure. To preserve the geometric shape of the circular holes and ensure continuity with the surrounding spline patches, four patches were allocated to each circular hole. Additionally, two spline patches were used at the central right-angle turn. Figure 19 presents a comparison of the stress calculations between the second and fourth mesh levels. The stress contour plots indicate minimal changes in stress results beyond improvements in smoothness. However, as shown in Fig. 20, error calculations and the corresponding error cloud diagrams demonstrate that refining the mesh in areas of stress concentration significantly reduces errors.

[See PDF for image]

Fig. 17

L-shape bracket model

[See PDF for image]

Fig. 18

Progression of adaptive refinement in cubic order elements. This figure includes six sub-figures, each illustrating the adaptive refinement process from level 1 to level 6 using cubic order elements

[See PDF for image]

Fig. 19

Computed stresses in multiple directions at different refinement levels. The first row, a–c, displays the computed stresses in three directions (xx, xy, yy) for the second refinement level. The second row, d–f, illustrates the same stress directions for the fourth refinement level

[See PDF for image]

Fig. 20

Analysis of computed recovery errors across refinement levels in three directions. The first row, a–c, presents the computed recovery errors in the xx, xy, and yy directions for the sixth refinement level, illustrating the accuracy of the model at this stage. The second row, d–f, shows the errors in the same three directions for the ninth refinement level, allowing for an evaluation of error reduction as refinement progresses. The third row, g–i, depicts the recovery errors for the tenth refinement level, providing a further analysis of error trends at higher levels of refinement

Perforated Beam

The perforated beam is a widely used structure in various engineering applications. Its primary advantage is weight and material reduction while maintaining sufficient safety for the structure itself. The main issue is illustrated in Figs. 21, and 22 shows a sequence of adaptive refinements. As observed, stress concentrations primarily occur around the hole regions and the fixed corners, aligning with mechanical principles. For this model, a total of 40 patches are applied to define it. Unlike the holes in the L-shaped bracket from the previous example, two holes in this scenario require six patches to account for the features of neighboring patches. Figure 23 displays the computed recovery errors across three levels, demonstrating the efficiency of adaptive refinement. Figure 24 presents the various orders of Gaussian collocation elements. Figure 25 compares the convergence rates based on recovery errors. It is evident that while uniform refinement levels off as refinement progresses, adaptive refinement may exhibit some fluctuations, yet the overall trend remains optimal across different orders.

[See PDF for image]

Fig. 21

Perforated beam model

[See PDF for image]

Fig. 22

Sequence of adaptive refinement meshes from level 1 to level 8. This figure contains eight sub-figures, each depicting the mesh at a different stage of adaptive refinement, ranging from level 1 to level 8

[See PDF for image]

Fig. 23

Analysis of computed recovery errors across refinement levels in three directions. The first row, a–c, presents the computed recovery errors in the xx, xy, and yy directions for the third refinement level, illustrating the accuracy of the model at this stage. The second row, d–f, shows the errors in the same three directions for the fourth refinement level, allowing for an evaluation of error reduction as refinement progresses. The third row, g–i, depicts the recovery errors for the fifth refinement level, providing a further analysis of error trends at higher levels of refinement

[See PDF for image]

Fig. 24

Comparison of Gaussian collocation elements across different polynomial orders. The figure includes 5 sub-figures, each comparing characteristics of Gaussian collocation elements from cubic to seventh order

[See PDF for image]

Fig. 25

Comparative analysis of convergence rates using IGA-GC and IGA-G across polynomial orders. The left column, a–e, displays the convergence rates achieved by IGA-GC method for polynomial orders from third to seventh. The right column, f–j, illustrates the convergence rates achieved by IGA-G method across the same range of polynomial orders

Honeycomb

The honeycomb structure, notable for its porous architecture comprising multiple cells. In the scenario described, 90 PHT spline patches are utilized to construct the structure, with the bottom layer secured and the top layer subject to stretching. To analyze this configuration, both IGA-GC and IGA-G are employed. Figure 26a–f illustrates various mesh divisions using cubic elements. Figure 27a–f presents the convergence rates.

[See PDF for image]

Fig. 26

Sequence of adaptive refinement meshes from level 1 to level 6. This figure contains six sub-figures, each depicting the mesh at a different stage of adaptive refinement, ranging from level 1 to level 6

[See PDF for image]

Fig. 27

Comparative analysis of convergence rates using IGA-GC and IGA-G across polynomial orders. The left column, a–c, displays the convergence rates achieved by IGA-GC method for polynomial orders from third to fifth. The right column, d–f, illustrates the convergence rates achieved by IGA-G method across the same range of polynomial orders

Discussion

In this section, we conduct numerical experiments using five models, ranging from a single annulus model to a complex 90-patch honeycomb structure. Some models feature sharp elements such as re-entrant corners, which can introduce singularities, while others incorporate multiple holes that may lead to stress concentration around these areas. One effective strategy to address these challenges is adaptive refinement, which allocates more elements in regions with sharp features and stress concentrations but fewer elsewhere, thus optimizing computational resources. Additionally, adaptive refinement can help achieve optimal convergence rates in problems characterized by singularities.

Furthermore, we explore the performance of high-order PHT elements in numerical testing. Both uniform and adaptive refinements of these elements demonstrate optimal convergence rates. The use of Gaussian collocation points with high-order elements proves particularly effective, especially in adaptive scenarios, helping to balance the number of basis functions and collocation points. We also detail the rules for eliminating and relocating collocation points around T-junctions in high-order elements. These rules have proven efficient in testing, further enhancing the robustness and accuracy of our simulations.

Conclusion

In this paper, we present a novel high-order Gaussian collocation method designed to analyze complex multi-patch structures. Our primary objective is to refine PHT for IGA-C and adaptive refinement. Traditional PHT splines, which normalize functions via basis function truncation at T-junctions, often suffer from a decay phenomenon. Our advancements in PHT elements retain the benefits of traditional elements while mitigating this decay issue and facilitating the extension to high-order approximations. We delve into the strategic selection of collocation points tailored for these high-order elements, especially in scenarios necessitating adaptive refinement. By leveraging the intrinsic properties of PHT elements, we utilize Gaussian points as collocation points. To prevent the system from becoming over-determined, it is crucial to either remove or relocate certain points near T-junctions. Consequently, we have formulated comprehensive guidelines for the distribution of collocation points across arbitrary high-order elements. Our method also addresses the seamless integration of multi-patch structures, ensuring continuity between patches by aligning control points, knots, and basis functions, as well as by synchronizing the refinement levels along the boundaries of adjacent patches. This connectivity is facilitated by a global indexing system, which enhances the structural integrity of the model. Numerical examples substantiate the efficacy of our method across a spectrum of structures, ranging from single patches to intricate multi-patch geometries featuring singularities such as re-entrant corners and multiple circular holes. Our method consistently achieves optimal convergence rates. While high-order approximations inherently increase computational demands, our findings suggest a balanced approach, optimizing between computational efficiency and accuracy. Looking ahead, our future work will concentrate on developing a three-dimensional solver for high-order approximations and applying the method to more practical engineering challenges, aiming to further validate and enhance its robustness and applicability. The proposed method and strategy can be extended for solution of the novel nanocomposite structures made of new nanofillers such as graphene nanoplatelets and graphene origami reinforced structures. Furthermore, one can suggest application of this method for analysis of the more complicated structures such as structures with curved boundaries and including holes.

Acknowledgements

This work was supported by the National Natural Science Foundation of China (11902263) and the Degree and Graduate Education Research Fund of NPU (2023YMs005).

Data availability statement

Data sharing not applicable to this article as no datasets were generated or analyzed during the current study.

Declarations

Conflict of interest

The authors declare no conflict of interest.

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

References

1. Hughes, TJR. The finite element method: linear static and dynamic finite element analysis; 2000; New York, Dover Publications: 1191.74002

2. Smith, GD. Numerical solution of partial differential equations: finite difference methods; 1985; 3 Oxford, Oxford University Press: 0576.65089

3. Cottrell, JA; Hughes, TJR; Bazilevs, Y. Isogeometric analysis: toward integration of CAD and FEA; 2009; New York, Wiley: [DOI: https://dx.doi.org/10.1002/9780470749081] 1378.65009

4. Hughes, TJR; Cottrell, JA; Bazilevs, Y. Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement. Comput Methods Appl Mech Eng; 2005; 194, 39–41 pp. 4135-4195.2152382 [DOI: https://dx.doi.org/10.1016/j.cma.2004.10.008] 1151.74419

5. Piegl, L; Tiller, W. The NURBS book; 1997; 2 New York, Springer: [DOI: https://dx.doi.org/10.1007/978-3-642-59223-2] 0868.68106

6. Casquero, H; Wei, X; Toshniwal, D; Li, A; Hughes, TJR; Kiendl, J; Zhang, YJ. Seamless integration of design and Kirchhoff-Love shell analysis using analysis-suitable unstructured T-splines. Comput Methods Appl Mech Eng; 2020; 360, 112765.4045265 [DOI: https://dx.doi.org/10.1016/j.cma.2019.112765] 1441.74111

7. Wang, Z; Cao, J; Wei, X; Chen, Z; Casquero, H; Zhang, YJ. Kirchhoff-Love shell representation and analysis using trangle configuration B-splines. Comput Methods Appl Mech Eng; 2023; 416, 116316. [DOI: https://dx.doi.org/10.1016/j.cma.2023.116316] 1536.74133

8. Popez, J; Anitescu, C; Rabczuk, T. Isogeometric structural shape optimization using automatic sensitivity analysis. Appl Math Model; 2021; 89, pp. 1004-1024.4145485 [DOI: https://dx.doi.org/10.1016/j.apm.2020.07.027] 1481.74631

9. Popez, J; Anitescu, C; Rabczuk, T. CAD-compatible structural shape optimization with movable Bézier tetrahedral mesh. Comput Methods Appl Mech Eng; 2020; 367, 113066. [DOI: https://dx.doi.org/10.1016/j.cma.2020.113066] 1442.74231

10. Qian, K; Liao, AS; Gu, S; Webster-Wood, VA; Zhang, YJ. Biomimetic IGA neuron growth modeling with neurite morphometric features and CNN-based prediction. Comput Methods Appl Mech Eng; 2023; 417, 116213.4668228 [DOI: https://dx.doi.org/10.1016/j.cma.2023.116213] 1539.92029

11. Li, A; Zhang, YJ. Isogeometric analysis-based physics-informed graph neural network for studying traffic jam in neurons. Comput Methods Appl Mech Eng; 2023; 403, 115757.4512465 [DOI: https://dx.doi.org/10.1016/j.cma.2022.115757] 1536.76156

12. Videla, J; Anitescu, C; Khajah, T; Bordas, SPA; Atroshchenko, E. h- and p-adaptivity driven by recovery and residual-based error estimators for PHT-splines applied to time-harmonic acoustics. Comput Math Appl; 2019; 77, 9 pp. 2369-2395.3945119 [DOI: https://dx.doi.org/10.1016/j.camwa.2018.12.026] 1442.65409

13. Auricchio, F; Beirão da Veoga, L; Hughes, TJR; Reali, A; Sangalli, G. Isogeometric collocation methods. Math Models Methods Appl Sci; 2010; 20, 11 pp. 2075-2107.2740715 [DOI: https://dx.doi.org/10.1142/S0218202510004878] 1226.65091

14. Farin, G. Curves and surface for CAGD: a practical guide; 2002; 5 Burlington, Morgan Kaufmann: 0817.68131

15. Anitescu, C; Jia, Y; Zhang, YJ; Rabczuk, T. An isogeometric collocation method using superconvergent points. Comput Methods Appl Mech Eng; 2015; 284, pp. 1073-1097.3310318 [DOI: https://dx.doi.org/10.1016/j.cma.2014.11.038] 1425.65193

16. de Boor, C; Schwartz, B. Collocation at Gaussian points. SIAM J Numer Anal; 1973; 10, 4 pp. 582-606.373328 [DOI: https://dx.doi.org/10.1137/0710052] 0232.65065

17. Jia, Y; Anitescu, C; Zhang, YJ; Rabczuk, T. An adaptive isogeometric analysis collocation method with a recovery-based error estimator. Comput Methods Appl Mech Eng; 2019; 345, pp. 52-74.3880137 [DOI: https://dx.doi.org/10.1016/j.cma.2018.10.039] 1440.65248

18. Sederberg, TW; Zheng, J; Bakenov, A. T-splines and T-NURCCs. ACM Trans Graphics; 2003; 22, 3 pp. 477-484. [DOI: https://dx.doi.org/10.1145/882262.882295] 1069.65542

19. Dokken, T; Lyche, T; Pettersen, KF. Polynomial splines over locally refined Box-partition. Comput Aided Geomet Des; 2013; 30, 3 pp. 331-356.3019748 [DOI: https://dx.doi.org/10.1016/j.cagd.2012.12.005] 1264.41011

20. Giannelli, C; Jüttler, B; Speleers, H. THB-splines: The truncated basis for hierarchical splines. Comput Aided Geomet Des; 2012; 29, 7 pp. 485-498.2925951 [DOI: https://dx.doi.org/10.1016/j.cagd.2012.03.025] 1252.65030

21. Deng, J; Chen, F; Li, X; Tong, CHuW; Yang, Z; Feng, Y. Polynomial splines over hierarchical T-meshes. Graph Models; 2008; 70, 4 pp. 76-86. [DOI: https://dx.doi.org/10.1016/j.gmod.2008.03.001] 1417.65069

22. Borden, MJ; Scott, MA; Evan, JA; Hughes, TJR. Isogeometric finite element data structures based on Bézier extraction of NURBS. Int J Numer Methods Eng; 2011; 87, pp. 15-47. [DOI: https://dx.doi.org/10.1002/nme.2968] 1242.74097

23. Scott, MA; Borden, MJ; Verhoosel, CV; Sederberg, TW; Hughes, TJR. Isogeometric finite element data structures based on Bézier extraction of T-splines. Int J Numer Methods Eng; 2011; 88, pp. 126-156. [DOI: https://dx.doi.org/10.1002/nme.3167] 1242.65243

24. Kang, H; Xu, J; Chen, F; Deng, J. A new basis for PHT-splines. Graph Models; 2015; 82, pp. 149-159. [DOI: https://dx.doi.org/10.1016/j.gmod.2015.06.011] 1349.94163

25. Wei X, Li X, Zhang YJ, Hughes TJR (2021) Tuned hybrid nonuniform subdivision surfaces with optimal convergence rates. 122(9):2117–2144

26. Chan, C; Anitescu, C; Rabczuk, T. Volumetric parametrization from a level set boundary representation with PHT-splines. Comput Aided Des; 2017; 82, pp. 29-41.3575510 [DOI: https://dx.doi.org/10.1016/j.cad.2016.08.008] 1475.68422

27. Babuška, I. The finite element method with penalty. Math Comput; 1973; 27, pp. 221-228.351118 [DOI: https://dx.doi.org/10.1090/S0025-5718-1973-0351118-5] 0299.65057

28. Babuška, I. The finite element method with Lagrangian multipliers. Numer Math; 1973; 20, 3 pp. 179-192.359352 [DOI: https://dx.doi.org/10.1007/BF01436561] 0258.65108

29. Guo, Y; Ruess, M. Nitsche’s method for a coupling of isogeometric thin shells and blended shell structures. Comput Methods Appl Mech Eng; 2015; 284, pp. 881-905.3310311 [DOI: https://dx.doi.org/10.1016/j.cma.2014.11.014] 1423.74573

30. Li, X; Deng, J; Chen, F. Surface modeling with polynomial splines over hierarchical T-meshes. Vis Comput; 2007; 23, 12 pp. 1027-1033. [DOI: https://dx.doi.org/10.1007/s00371-007-0170-3] 1364.65141

31. Li, X; Deng, J; Chen, F. Polynomial splines over general T-meshes. Virtual Comput; 2010; 26, 4 pp. 277-286.0955.62581

32. Anitescu, C; Hossian, MN; Rabczuk, T. Recovery-based error estimation and adaptivity using high-order splines over hierarchical T-meshs. Comput Methods Appl Mech Eng; 2018; 328, pp. 638-662. [DOI: https://dx.doi.org/10.1016/j.cma.2017.08.032] 1439.65011

33. Videla J, Contreras F, Nguyen HX, Atroshchenko E (2020) Application of PHT-splines in bending and vibration analysis of cracked Kirchhoff-Love plates 361:112754

34. Yang, HS; Dong, CY; Qin, XC; Wu, YH. Viration and buckling analyses of FGM plates with multiple internal defects using XIGA-PHT and FCM under thermal and mechanical loads. Appl Math Model; 2020; 78, pp. 433-481.4029064 [DOI: https://dx.doi.org/10.1016/j.apm.2019.10.011] 1481.74330

35. Yang, HS; Dong, CY. Adaptive extended isogeometric analysis based on PHT-splines for thin cracked plates and shells with Kirchhoff-Love theory. Appl Math Model; 2019; 76, pp. 759-799.3980884 [DOI: https://dx.doi.org/10.1016/j.apm.2019.07.002] 1481.74682

36. Zienkiewicz, OC; Zhu, JZ. The superconvergence patch recovery and a posteriori error estimation in the finite element method, Part 1: the recovery technique. Int J Numer Methods Eng; 1992; 33, 7 pp. 1331-1364. [DOI: https://dx.doi.org/10.1002/nme.1620330702] 0769.73084

37. Zienkiewicz, OC; Zhu, JZ. The superconvergence patch recovery and a posteriori error estimation in the finite element method, Part 2: error estimates and adaptivity. Int J Numer Methods Eng; 1992; 33, 7 pp. 1365-1382. [DOI: https://dx.doi.org/10.1002/nme.1620330703] 0769.73085

38. Sayyad, AS; Mahajan, VM; Shinde, BM. Effects of transverse normal strain on the deformation of laminated and sandwich arches under the action of concentrated force. Mech Adv Mater Struct; 2023; [DOI: https://dx.doi.org/10.1080/15376494.2023.2229825] 07815808

39. Mahajan, VM; Sharma, A. Evaluation of static and dynamic responses considering thickness stretching effect for layered composite and sandwich arches using exponential shear and normal deformation theory. Forces Mech; 2023; 12, 100204. [DOI: https://dx.doi.org/10.1016/j.finmec.2023.100204] 1481.35053

40. Ren, J; Lin, H. A survey on isogeometric collocation methods with applications. Mathematics; 2023; 11, 469. [DOI: https://dx.doi.org/10.3390/math11020469] 1203.76062

41. Mahajan, VM; Sharma, A. Evaluation of static responses for layered composite arches. Curved Layer Struct; 2023; 10, 1 20220185. [DOI: https://dx.doi.org/10.1515/cls-2022-0185] 07866011

42. Jiang, X; Jiang, Y; Lu, S et al. Free vibration of three-dimensional angle-interlock woven composite multilayered cantilever beam. J Vib Eng Technol; 2023; 11, pp. 3387-3398. [DOI: https://dx.doi.org/10.1007/s42417-022-00755-x] 1443.34040

43. Pathak, D; Kushari, S; Maity, S et al. Vibration analysis of cracked cantilever beam using response surface methodology. J Vib Eng Technol; 2023; 11, pp. 2429-2452. [DOI: https://dx.doi.org/10.1007/s42417-022-00713-7] 1297.93145

44. Liu, Y. Nonlinear dynamic analysis of an axially moving composite laminated cantilever beam. J Vib Eng Technol; 2023; 11, pp. 3307-3319. [DOI: https://dx.doi.org/10.1007/s42417-022-00750-2] 1213.68499

45. Kurt, P; Orhan, S. Improving energy harvesting from cantilever-like structures based on beam geometry. J Vib Eng Technol; 2024; [DOI: https://dx.doi.org/10.1007/s42417-024-01326-y] 07906597

46. Xu, X; Fu, X; Zhao, H; Liu, M; Xu, A; Ma, Y. Three-dimensional reconstruction and geometric morphology analysis of lunar small craters within the patrol range of the Yutu-2 rover. Remote Sens; 2023; 15, 17 4251. [DOI: https://dx.doi.org/10.3390/rs15174251]

47. Wang, Z; Chen, M; Guo, Y; Li, Z; Yu, Qa. Bridging the domain gap in satellite pose estimation: a self-training approach based on geometrical constraints. IEEE Transact Aero Electro Syst; 2024; 60, 3 pp. 2500-2514. [DOI: https://dx.doi.org/10.1109/TAES.2023.3250385] 1484.92028

48. Wu, Y; Fan, Y; Zhou, S; Wang, X; Chen, Q; Li, X. Research on the cross-sectional geometric parameters and rigid skeleton length of reinforced concrete arch bridges: a case study of Yelanghu Bridge. Struct; 2024; 69, 107423. [DOI: https://dx.doi.org/10.1016/j.istruc.2024.107423] 1537.16022