Isogeometric convolution hierarchical deep-learning neural network: isogeometric analysis with versatile adaptivity

Abstract

A method of Convolution-Hierarchical Deep-learning Neural Network isogeometric Analysis (C-IGA) of a geometric shape of a subject matter comprises performing a first mapping process between a physical domain of a subject matter and a parametric domain of the subject matter; and performing a second mapping process between the parametric domain of the subject matter and a parent domain of the subject matter; wherein the first mapping process comprises constructing an C-IGA interpolation in the parametric domain based on CAD data of the subject matter.

Description

Attorney Docket No.0116936.296WO2 ISOGEOMETRIC CONVOLUTION HIERARCHICAL DEEP-LEARNING NEURAL NETWORK: ISOGEOMETRIC ANALYSIS WITH VERSATILE ADAPTIVITY CROSS-REFERENCE TO RELATED PATENT APPLICATIONS

[0001] This application claims priority to and the benefit of U.S. Provisional Patent Application Serial No.63 / 534,609, filed August 25, 2023, which is incorporated herein in its entirety by reference. FIELD OF THE INVENTION

[0002] The present disclosure relates generally to the field of Hierarchical Deep Learning Neural Networks which pertains to a method, algorithm, and computer program employed in engineering analysis, and more particularly, to isogeometric analysis with versatile adaptivity. BACKGROUND OF THE INVENTION

[0003] The background description provided herein is for the purpose of generally presenting the context of the present invention. The subject matter discussed in the background of the invention section should not be assumed to be prior art merely as a result of its mention in the background of the invention section. Similarly, a problem mentioned in the background of the invention section or associated with the subject matter of the background of the invention section should not be assumed to have been previously recognized in the prior art. The subject matter in the background of the invention section merely represents different approaches, which in and of themselves may also be inventions.

[0004] The field of computational science and engineering is witnessing a paradigm shift in both algorithms and software. Computational scientists are developing methods that rely on training neural networks from experimental or computational data. Andrej Karpathy defined this new generation of software as software 2.0. Later, Erik Meijer elaborated on the idea by explaining two drastically different approaches of software development: software 1.0 where the developers need to specify the problem, design the algorithms, and break down a complex task into multiple smaller components, and software 2.0 where the developers use data, design an optimization problem, and train a neural network to solve any problem.Attorney Docket No.0116936.296WO2

[0005] For real-time control and optimization of the manufacturing processes, one needs both a fast and differentiable tool, as it involves computing derivatives of the Lagrangian many times. Currently, many works rely on reduced-order method to avoid computing a more accurate but computationally intensive solution of the governing equations. However, minimization of Lagrangian form of governing equation using Convolution Hierarchical Deep-learning Neural Network Artificial Intelligence (C-HiDeNN-AI), and applying the C-HiDeNN-AI in the field of isogeometric analysis are still absent.

[0006] Therefore, a heretofore unaddressed need exists in the art to address the aforementioned deficiencies and inadequacies. SUMMARY OF THE INVENTION

[0007] In light of the foregoing, this invention discloses a method of Convolution- Hierarchical Deep-learning Neural Network Isogeometric Analysis (C-IGA) of a geometric shape of a subject matter. The method comprises performing a first mapping process between a physical domain of a subject matter and a parametric domain of the subject matter; and performing a second mapping process between the parametric domain of the subject matter and a parent domain of the subject matter; wherein the first mapping process comprises constructing a C-IGA interpolation in the parametric domain based on CAD data of the subject matter.

[0008] In one embodiment, the CAD data of the subject matter comprises at least one of knot vectors, polynomial orders for each parametric direction, control points, and weights for spline basis functions.

[0009] In one embodiment, the spline basis functions comprise Non-Uniform Rational B- Splines (NURBS) basis functions.

[0010] In one embodiment, the step of the first mapping process comprises mapping nodal points though the geometric mapping.

[0011] In one embodiment, the nodal points are discrete.

[0012] In one embodiment, the C-IGA interpolation in the parametric domain is constructed using: xC−IGA(ξ) = ∑k∈AesNk(ξ)xGeo(ξk) .freedom are attached to the nodal points.

[0014] In one embodiment, the C-IGA interpolation is constructed based on at least oneAttorney Docket No.0116936.296WO2 calculation parameter.

[0015] In one embodiment, the calculation parameter comprises a mesh coordinate r, a mesh size h, a patch size s, a dilation parameter a, and reproducing order p.

[0016] In one embodiment, the reproducing order p and spline basis functions corresponding to the reproducing order p comprises the spline basis functions used in the CAD data.

[0017] In one embodiment, the CAD data comprises multiple CAD patches.

[0018] In one embodiment, the multiple CAD patches are linked via at least one linking method.

[0019] In one embodiment, the at least one linking method comprises a nodal compatibility method. The nodal compatibility method comprises building a conforming mesh in the physical domain; and conducting an inverse IGA mapping to construct a parametric mesh.

[0020] In one embodiment, the at least one linking method comprises a G0compatibility method. The G0compatibility method comprises building a conforming mesh in the physical domain; conducting an inverse IGA mapping to construct a parametric mesh; and reconstructing a regular mesh in vicinity of an interface of one of the multiple CAD patches in the parametric domain; wherein the C-IGA interpolant is constructed in the regular mesh.

[0021] In one embodiment, the at least one linking method comprises a nodal compatibility method, wherein the nodal compatibility method comprises: building a conforming mesh in the physical domain; and conducting an inverse IGA mapping to construct a parametric mesh.

[0022] In one embodiment, the at least one linking method comprises a G0compatibility method, wherein the G0compatibility method comprises building a conforming mesh in the physical domain; conducting an inverse IGA mapping to construct a parametric mesh; and reconstructing a regular mesh in the parametric domain; wherein the C-IGA interpolants is constructed in the regular mesh.

[0023] In another aspect of the invention, a method for graphics processing units (GPU) programming of Convolution-Hierarchical Deep-learning Neural Network Isogeometric Analysis (C-IGA) of a geometric shape of a subject matter comprises (1) providing a mesh to obtain a Hierarchical Deep-learning Neural Network – Finite Element Method (HiDeNN-FEM) interpolant based on the mesh; (2) computing an elemental patch information and a nodal connectivity based on the mesh and graph theory; (3) computing an assembled moment matrix G in a parametric domain using the elemental patch information and a calculation parameter; (4)Attorney Docket No.0116936.296WO2 obtaining a convolutional patch function in the parametric domain based on the steps (2) and (3); (5) obtaining a shape function in a physical domain based on the convolutional patch function; (6) integrating the shape function with the HiDeNN-FEM interpolant obtained in the step (1) to generate an isoparametric C-HiDeNN interpolant; (7) building a global stiffness matrix and force vector system of equations using the C-HiDeNN interpolant; and (8) solving the global stiffness matrix and force vector system of equations using an equation solver.

[0024] In one embodiment, the convolution patch function is obtained according to nodal connectivity and a calculation parameter.

[0025] In one embodiment, the calculation parameter comprises a mesh coordinate r, a mesh size h, a patch size s, a dilation parameter a, and reproducing order p.

[0026] In one embodiment, the steps (1) and (3)-(8) are parallelized.

[0027] In one embodiment, the steps (3)-(7) are serially conducted and stacked in a GPU memory having a memory limit; wherein the parallelized steps (3)-(7) comprises at least one variable having a variable size; and wherein when the variable size exceeding the memory limit of GPU, the steps (3)-(7) are serially conducted and stacked in a CPU memory.

[0028] In another aspect of the invention, a non-transitory tangible computer readable medium storing a program which, when executed by one or more processors, causes a computer to execute a process for GPU programming of Convolution-Hierarchical Deep-learning Neural Network Isogeometric Analysis (C-IGA) of a geometric shape of a subject matter comprises (1) providing a mesh to obtain a Hierarchical Deep-learning Neural Network – Finite Element Method (HiDeNN-FEM) interpolant based on the mesh; (2) computing an elemental patch information and a nodal connectivity based on the mesh and graph theory; (3) computing an assembled moment matrix G in a parametric domain using elemental patch information and a calculation parameter; (4) obtaining a convolutional patch function in the parametric domain based on the steps (2) and (3); (5) obtaining a shape function in a physical domain based on the convolutional patch function; (6) integrating the shape function with the HiDeNN-FEM interpolant obtained in the step (1) to generate a C-IGA interpolant; (7) building a global stiffness matrix and force vector system of equations using the C-IGA interpolant; and (8) solving the global stiffness matrix and force vector system of equations using an equation solver. BRIEF DESCRIPTION OF THE DRAWINGSAttorney Docket No.0116936.296WO2

[0029] The accompanying drawings illustrate one or more embodiments of the invention and together with the written description, serve to explain the principles of the invention. Wherever possible, the same reference numbers are used throughout the drawings to refer to the same or like elements of an embodiment.

[0030] Fig.1 shows a diagram of components of Convolution Hierarchical Deep-learning Neural Network Artificial Intelligence system for engineering software 2.0.

[0031] Fig.2 shows a diagram of 2D-C-HiDeNN isoparametric mapping between the parent domain and the physical domain.

[0032] Fig.3 shows a diagram of neural network structure of C-HiDeNN for a 1D scientific problem; panel (a) shows nodal coordinates in the parent domain and physical domain whose mapping is defined by the C-HiDeNN interpolation in an isoparametric manner; panel (b) shows elemental neural network of C-HiDeNN, with the first two hidden layers have linear or ReLU activation functions while the third hidden layer has Gaussian or cubic spline activation functions that enable nonlinear interpolation; panel (c) shows the C-HiDeNN structure for the entire problem, with the variational energy is used as a loss function of the problem.

[0033] Fig.4 shows an overview of Isogeometric Analysis (IGA) for a 3D toroidal surface.

[0034] Fig.5 shows a diagram of C-IGA mappings between the physical domain, parametric domain, and the parent domain. The mapping from the parent domain to the parametric domain is the standard FEM mapping used for numerical integration like IGA. The parametric domain shows the knot space defined in IGA framework.

[0035] Fig.6 shows a diagram of neural network structure of C-IGA for a 1D problem; panel (a) shows nodal coordinates in the parent domain, parametric domain, and physical domain; panels (b) and (c) show the elemental neural network of C-IGA (b) and the global neural network of C-IGA (c).

[0036] Fig.7 shows a diagram of comparison between IGA and C-IGA mappings; panel (a) shows IGA mapping is built on a regular mesh because the B-splines and NURBS basis functions are constructed by tensor product; panel (b) shows C-IGA mapping can be built on an irregular mesh, enabling the automatic r-adaptivity. Like conventional FEM or C-HiDeNN, the DoFs of C-IGA are attached to nodes.

[0037] Fig.8 shows charts reflecting the body force b(x) and analytical solution u(x) for 1D numerical model.Attorney Docket No.0116936.296WO2

[0038] Fig.9 shows a diagram of 1D numerical model and illustrations for C-IGA interpolants.

[0039] Fig.10 shows charts reflecting the L2-norm error comparison for C-IGA with varying s and fixed p = 2; panel (a) shows cubic spline radial function and panel (b) Gaussian radial function in Eq. (6).640 elements are used. hξis the element size in the parametric domain.

[0040] Fig.11 shows charts reflecting the L2-norm error convergence rates against the mesh refinement for different methods and different sets of parameters for C-IGA with cubic spline radial function: panel (a) p = 2 (L2-norm error); panel (b) p = 3 (L2-norm error); panel (c) p = 4 (L2-norm error).

[0041] Fig.12 shows charts reflecting optimized mesh for C-IGA with r-adaptivity: initial mesh size (a) h = L / 40; (b) h = L / 80. The parameters are set to be s = 2, p = 2, a / hξ = 1.8.

[0042] Fig.13 shows charts reflecting numerical results for C-IGA with r-adaptivity: (a) p = 2; (b) p = 3; (c) p = 4. The parameters are set to be s = 2, a / hξ = 1.8 for p = 2, s = 3, a / hξ = 1.86 for p = 3, and s = 4, a / hξ = 9.9 for p = 4.

[0043] Fig.14 shows charts reflecting numerical results for C-IGA with a-adaptivity: (a) p = 2; (b) p = 3; (c) p = 4.

[0044] Fig.15 shows illustrations numerical results for p-adaptivity: panel (a) illustrations of 1D model for p-adaptivity; panel (b) convergence study for p-adaptivity. The fourth order C-IGA interpolants are used in the blue region (x ∈ [1.5, 3.5] ∪ [6.5, 8.5]), and the second order interpolants for other regions.

[0045] Fig.16 shows 2D model with a body force centered at (5, 15) and geometric mapping; panel (a) Let IGA mapping with NURBS basis functions be the geometric mapping xGeo(ξ); panel (b) Isoparametric mapping of C-IGA transforms the physical mesh to the parametric mesh.

[0046] Fig.17 shows mesh comparisons among linear FEM, IGA and C-IGA (p = 2); IGA and C-IGA meshes can reproduce the geometry exactly, whereas linear FEM mesh approximates the geometry with straight lines.

[0047] Fig.18 shows exact solution (a) versus C-IGA solution with 160 by 160 elements (b).

[0048] Fig.19 shows charts reflecting convergence study results for 2D Poisson’s problem with cubic spline radial function and a / h = 50. (a-c) Energy norm error for p = 2 vs. (a) mesh sizeAttorney Docket No.0116936.296WO2 hξ, (b) DoFs, and (c) computation time. (d-f) Energy norm error for p = 3 vs. (d) mesh size hξ, (e) DoFs, and (f) computation time. Convergence rates are denoted in (a, d). hξ is the element size in the parametric domain. In (b,c), horizontal lines are drawn at 10−4 error to compare the efficiency of C-IGA (s = 5, p = 2) with IGA and FEM (p = 2). FEM (p =FEM (p = 3) use serendipity elements: 8-node and 12-node elements, respectively.

[0049] Fig.20 shows charts reflecting convergence study results for 2D Poisson’s problem with Gaussian radial function and a / hξ = 3.2; (a-c) Energy norm error for p = 2 vs. (a) mesh size hξ, (b) DoFs, and (c) computation time. Degrees of freedom (DoFs) are marked. In (b,c), horizontal lines are drawn at 10−4error to compare the efficiency of C-IGA (s = 5, p = 2) with IGA and FEM (p = 2).

[0050] Fig.21 shows charts reflecting C-IGA mapping from (b) the irregular parametric mesh to (a) the physical mesh for 5 by 5 elements; (c) Solution field of 128 by 128 irregular mesh in the parametric domain; (d) Convergence plot for IGA (p = 2), C-IGA (s = 2, p = 2) with regular and irregular mesh. The cubic spline radial function is used with a / hξ = 50.

[0051] Fig.22 shows charts reflecting C-IGA mapping from the parametric mesh with triangular elements (b) to the physical mesh (a). (c) Error vs. DoFs plot for IGA (p = 2), C-IGA (s = 3 and s = 5) with quadrilateral and triangular elements. The cubic spline radial function is used with a / hξ= 50.

[0052] Fig.23 shows charts reflecting 2D linear elastic example: infinite plate with a circular hole; panel (a) shows problem definition with parameters. Convergence study results: energy norm error vs. (b) hξ and (c) DoFs. Gaussian radial function is used for the convergence study with a / hξ= 8. C-IGA mapping from (e) the parametric mesh to (d) the physical mesh. The parametric mesh size hξ is denoted. Two IGA elements from the lowest order are denoted in (e): element 1 - [0, 0.5] X [0, 1]; element 2 - [0.5, 1] X [0, 1]. Panel (f) shows σxx plot (unitless) of 64 by 64 regular elements in the parametric domain. The Gaussian radial function is used with a / hξ= 8, s = 2, p = 2.

[0053] Fig.24 shows a chart reflecting convergence study of 2D linear elasticity with plane strain and Poisson’s ratio of 0.4999. Convergence rates are denoted. A severe locking is observed from linear FEM (p = 1). The Gaussian radial function uses a / hξ= 8 with varying p and s. hξ is the element size in the parametric domain.

[0054] Fig.25 shows a diagram of 3D linear elastic example: solid circular cylinderAttorney Docket No.0116936.296WO2 subjected to internal pressure loading. Physical mesh (a) and Parametric mesh (b) of 32 x 16 x 8 elements. Ri= 1, Ro= 2, and L = 5. The top and bottom surfaces of the cylinder are given symmetric boundary conditions. (c) Convergence study results for L2 norm error. hξ is the element size in the parametric domain. The cubic spline radial function is used with a / hξ = 50.

[0055] Fig.26 shows a diagram of the geometrical mapping of isogeometric C-HiDeNN-TD.

[0010] proposed a 2-step mapping to transform an irregular physical domain to a regular reference domain by FEM mapping and an affine. The regular domain is intrinsically separable for TD. C- IGA uses geometric mapping instead, e.g., IGA mapping, which is mature and straightforward. (a) 2-step mapping in

[0010] . (b) Geometric mapping of isogeometric C-HiDeNN-TD.

[0056] Fig.27 shows a diagram of a numerical example where a 3D hyper-elastic plate with a hole is under tension in the x direction; panel (a) shows tension of 3D hyperelastic plate with a hole; and panel (b) shows Maximum Mises stress as a function of the number of DoFs.

[0057] Fig.28 shows a diagram of the resulting optimized structure for nested TO design ([n] = 7); (a) Initial design region and the optimized nested scales topology design showing a bamboo-like structure, (b) seven-scale self-similar topology design with the same volume fraction in each scale, and (c) seven-scale non self-similar design topology design with a varying volume fraction at different scales.

[0058] Fig.29 shows a diagram of performance of the nested lattice structure design.

[0059] Fig.30 shows a diagram of 3D-printed bamboo-like structure for (a) self-similar nested lattice design, (b) non self-similar nested lattice design. A scaled lattice structure is printed for [n]=2,..,7 to demonstrate the optimized structure.

[0060] Fig.31 shows a diagram of neural network structure of HiDeNN for a 1D scientific problem. Panel (a) shows nodal coordinates in the physical domain; panel (b) shows elemental neural network that mimics the linear shape function; panel (c) shows the HiDeNN structure for the entire problem.

[0061] Fig.32 shows charts reflecting the H1-norm error convergence rates against the mesh refinement for different methods and different sets of parameters for C-IGA with cubic spline radial function: (a) p = 2 (H1-norm error); (b) p = 3 (H1-norm error); (c) p = 4 (H1-norm error).

[0062] Fig.33 shows charts reflecting condition number of stiffness matrix for different meshes; panel (a) shows condition number vs. h; panel (b) shows accuracy vs. condition number.

[0063] Fig.34 shows charts reflecting anchors in IGA; panel (a) shows quadratic basisAttorney Docket No.0116936.296WO2 functions (p = 2) and anchors for open knot vector Ξ = {0, 0, 0, 1, 2, 3, 4, 5, 5, 5}; panel (b) shows cubic basis functions (p = 3) and anchors for open knot vector Ξ = {0, 0, 0, 0, 1, 2, 3, 4, 5, 5, 5, 5}; anchors t1, t2 locates at the same parametric position; panel (c) shows modified anchors by redistributing repeating anchors. t2 is moved to the middle point of t1 and t3.

[0064] Fig.35 shows charts reflecting numerical results of C-IGA (Eq. (G.2)) and IGA; in which C-IGA results completely coincide with IGA results.

[0065] Fig.36 shows a flow chart of steps of C-IGA algorithm implemented on a GPU.

[0066] Fig.37 shows two-dimensional (top) and three-dimensional cases (bottom) for a two- CAD-patch system. Physical space (left) and parametric space (right). Yellow box represents patch ^^ ( ^^( ^^)), and red box represents patch ^^ ( ^^( ^^)).

[0067] Fig.38 provides illustration for compatibility conditions: panel (a) Physical space and ^^(parametric space; the nodal convolution patches ^^^^1)^^(2)^^, ^^^^^^for shared node ^^ are colored by ^^(^1)^^(2)^^(1)^^(2)green. ^^^^^, ^^^^^^are divided into the index sets ^^^^^^, ^^^^^^of internal nodes and the index set ^^^^^^ ^^of shared nodes along the interface; panel (b) Illustration for condition (C1): FEM shape ^^(^^1), ^^(^^2)for shared node ^^ match at the interface; panel (c) Illustration for condition () ( )(C2):patch functions ^^^^ 1 ^^ 2^^^^, ^^^^,^^, ^^, ^^^^, ^^^^,^^, ^^for node ^^ ∈ ^^^^^^match at the interface; ion (C3): patch^^^^( ^^)^^( ^^)panel (d) Illustration for condit^^, ^^^^,^^, ^^^^1 , ^^^^, ^^, ^^, ^^2for node^^(1) (2)^^^^^^1 ∈ ^^^^, ^^2∈ ^^^^^^vanish at the interface.C-IGA shape functions which are constructed by the product rule around the interface; in the region represented by a blue dashed box, regular meshes are established in the parametric domain ^^ where the FEM shape functions, and the convolution patch functions are constructed by the product rule defined in Eqs. ( 81 )-( 82 ). The nodal index ^^ is represented by ( ^^, ^^) with ^^, ^^ representing the index in ^^-direction and ^^-direction, respectively.

[0069] Fig.40 provides illustrations for 1D interpolation problem to simulate the interfaces of two patches. Two different mappings ^^(1)and ^^(2)generates two different interpolations ^^(1)and ^^(2), corresponding to interpolations in two neighboring patches. ^^(1)and ^^(2)share the same physical domain. The physical nodes are denoted by [ ^^ , ^^2, ^^3, … ],Attorney Docket No.0116936.296WO2 parametric nodes [ ^^(11), ^^(21), ^^(31), … ] for the first mapping ^^(1), and [ ^^(11), ^^(22), ^^(32), … ] for the second mapping ^^

[0070] Fig. reflecting convergence study of between ^^(1)(^^)and ^^(2)(^^). (a) Case 1: ^^(1)(^^)and ^^(2)(^^)are interpolations of a smooth function, ^^( ^^) = ^^ ^^ ^^ ^^. (b) Case 2: nodal values are assigned to be [ ^^1, ^^2, ^^3, ^^4, … ] = [1,2,1,2, … ]. In Case 1, the ^^2norm relative deviation reaches the ( ^^ + 1)-th order theoretical convergence rate. In Case 2, the convergence rate is only a first order.

[0071] Fig.42 shows charts reflecting convergence study of multi-patch C-IGA for interpolating a continuous function: panel (a) Geometrical description of a 2-patch problem; a plate with a circular hole; panel (b) Convergence rate with distinct polynomial order(^^)for each patch.

[0072] Fig.43 shows charts reflecting convergence study of multi-patch C-IGA with nodal compatibility for solving a PDE (linear elasticity equation): panel (a) Energy norm error of the entire domain ^^; panel (b) ^^2norm deviation at the interface.

[0073] Fig.44 shows the geometry of the problem domain in panel (a) and the exact solution in panel (b). The interpolated solution field in these two patches are denoted by ^^(1)(^^, ^^)and ^^(2)(^^, ^^), respectively. The Gaussian hump locates near the interfacered dashed curve in panel (a).

[0074] Fig.45 shows charts in panel (a) In the coarsest mesh, two IGA elements are defined in patch 1 and one IGA element is needed in patch 2; panel (b) Patch 2 is divided into two parts for G0compatibility.

[0075] Fig.46 shows charts in panel (a) Approximate solutions ^^(1)(patch 1) and ^^(2)(patch 2) along the interface in a polar coordinate system; panel (b)rate with distinct polynomial order ( ^^) for each patch. ^^ = ^^ is set for every case.

[0076] Fig.47 illustrates C-IGA shape functions with ^^0compatibility at the interface: panel (a) FEM shape functions ^^(1)and ^^(2)with ^^ = 2, ^^ = 2; panel (b) convolution ch functions ^^^^(1)(pat^^=8,^^, ^^, ^^,( ^^=9, ^^=15) ^^, ^^, ^^,( ^^=9, ^^=15)with ^^ = 2, ^^ = 2; panel (c) C-IGA shape functions ^^(1)and ^^(2)with ^^ = 2, ^^ = 2; panel (d) FEM shape functionsAttorney Docket No.0116936.296WO2 ^^(1)(^^=8, ^^=15)and ^^(2)(^^=8, ^^=15)with ^^ = 3, ^^ = 3; panel (e) convolution patch functions

[0077] The invention will now be described more fully hereinafter with reference to the accompanying drawings, in which exemplary embodiments of the invention are shown. This invention may, however, be embodied in many different forms and should not be construed as limited to the embodiments set forth herein. Rather, these embodiments are provided so that this invention will be thorough and complete, and will fully convey the scope of the invention to those skilled in the art. Like reference numerals refer to like elements throughout.

[0078] The terms used in this specification generally have their ordinary meanings in the art, within the context of the invention, and in the specific context where each term is used. Certain terms that are used to describe the invention are discussed below, or elsewhere in the specification, to provide additional guidance to the practitioner regarding the description of the invention. For convenience, certain terms may be highlighted, for example using italics and / or quotation marks. The use of highlighting has no influence on the scope and meaning of a term; the scope and meaning of a term is the same, in the same context, whether or not it is highlighted. It will be appreciated that same thing can be said in more than one way. Consequently, alternative language and synonyms may be used for any one or more of the terms discussed herein, nor is any special significance to be placed upon whether or not a term is elaborated or discussed herein. Synonyms for certain terms are provided. A recital of one or more synonyms does not exclude the use of other synonyms. The use of examples anywhere in this specification including examples of any terms discussed herein is illustrative only, and in no way limits the scope and meaning of the invention or of any exemplified term. Likewise, the invention is not limited to various embodiments given in this specification.

[0079] One of ordinary skill in the art will appreciate that starting materials, biological materials, reagents, synthetic methods, purification methods, analytical methods, assay methods, and biological methods other than those specifically exemplified can be employed in the practiceAttorney Docket No.0116936.296WO2 of the invention without resort to undue experimentation. All art-known functional equivalents, of any such materials and methods are intended to be included in this invention. The terms and expressions which have been employed are used as terms of description and not of limitation, and there is no intention that in the use of such terms and expressions of excluding any equivalents of the features shown and described or portions thereof, but it is recognized that various modifications are possible within the scope of the invention claimed. Thus, it should be understood that although the invention has been specifically disclosed by preferred embodiments and optional features, modification and variation of the concepts herein disclosed may be resorted to by those skilled in the art, and that such modifications and variations are considered to be within the scope of this invention as defined by the appended claims.

[0080] Whenever a range is given in the specification, for example, a temperature range, a time range, or a composition or concentration range, all intermediate ranges and subranges, as well as all individual values included in the ranges given are intended to be included in the invention. It will be understood that any subranges or individual values in a range or subrange that are included in the description herein can be excluded from the claims herein.

[0081] It will be understood that, as used in the description herein and throughout the claims that follow, the meaning of “a”, “an”, and “the” includes plural reference unless the context clearly dictates otherwise. Thus, for example, reference to “a cell” includes a plurality of such cells and equivalents thereof known to those skilled in the art. As well, the terms “a” (or “an”), “one or more” and “at least one” can be used interchangeably herein. It is also to be noted that the terms “comprising”, “including”, and “having” can be used interchangeably.

[0082] It will be understood that when an element is referred to as being “on”, “attached” to, “connected” to, “coupled” with, “contacting”, etc., another element, it can be directly on, attached to, connected to, coupled with or contacting the other element or intervening elements may also be present. In contrast, when an element is referred to as being, for example, “directly on”, “directly attached” to, “directly connected” to, “directly coupled” with or “directly contacting” another element, there are no intervening elements present. It will also be appreciated by those of skill in the art that references to a structure or feature that is disposed “adjacent” another feature may have portions that overlap or underlie the adjacent feature.

[0083] It will be understood that, although the terms first, second, third etc. may be used herein to describe various elements, components, regions, layers and / or sections, these elements,Attorney Docket No.0116936.296WO2 components, regions, layers and / or sections should not be limited by these terms. These terms are only used to distinguish one element, component, region, layer or section from another element, component, region, layer or section. Thus, a first element, component, region, layer or section discussed below could be termed a second element, component, region, layer or section without departing from the teachings of the invention.

[0084] Furthermore, relative terms, such as “lower” or “bottom” and “upper” or “top,” may be used herein to describe one element’s relationship to another element as illustrated in the figures. It will be understood that relative terms are intended to encompass different orientations of the device in addition to the orientation depicted in the figures. For example, if the device in one of the figures is turned over, elements described as being on the “lower” side of other elements would then be oriented on “upper” sides of the other elements. The exemplary term “lower”, can therefore, encompasses both an orientation of “lower” and “upper,” depending of the particular orientation of the figure. Similarly, if the device in one of the figures is turned over, elements described as “below” or “beneath” other elements would then be oriented “above” the other elements. The exemplary terms “below” or “beneath” can, therefore, encompass both an orientation of above and below.

[0085] It will be further understood that the terms “comprises” and / or “comprising”, or “includes” and / or “including”, or “has” and / or “having”, or “carry” and / or “carrying”, or “contain” and / or “containing”, or “involve” and / or “involving”, “characterized by”, and the like are to be open-ended, i.e., to mean including but not limited to. When used in this disclosure, they specify the presence of stated features, regions, integers, steps, operations, elements, and / or components, but do not preclude the presence or addition of one or more other features, regions, integers, steps, operations, elements, components, and / or groups thereof.

[0086] Unless otherwise defined, all terms (including technical and scientific terms) used herein have the same meaning as commonly understood by one of ordinary skill in the art to which this invention belongs. It will be further understood that terms, such as those defined in commonly used dictionaries, should be interpreted as having a meaning that is consistent with their meaning in the context of the relevant art and the invention, and will not be interpreted in an idealized or overly formal sense unless expressly so defined herein.

[0087] As used in the disclosure, “around”, “about”, “approximately” or “substantially” shall generally mean within 20 percent, preferably within 10 percent, and more preferably within 5Attorney Docket No.0116936.296WO2 percent of a given value or range. Numerical quantities given herein are approximate, meaning that the term “around”, “about”, “approximately” or “substantially” can be inferred if not expressly stated.

[0088] As used in the disclosure, the phrase “at least one of A, B, and C” should be construed to mean a logical (A or B or C), using a non-exclusive logical OR. As used herein, the term “and / or” includes any and all combinations of one or more of the associated listed items.

[0089] As used in the disclosure, the term “wearable” refers to articles, adornments or items designed to be worn by a user, incorporated into another item worn by a user, act as an orthosis for the user, or interfacing with the contours of a user's body.

[0090] As used in the disclosure, “biocompatible” material is a material that is compatible with living tissue or a living system by not being toxic or injurious and not causing immunological rejection.

[0091] As used in the disclosure, the term “therapy” refers to any protocol, method, and / or agent that can be used in the management, treatment, and / or amelioration of a given disease, or a symptom related thereto. In certain embodiments, the terms “therapies” and “therapy” refer to a biological therapy, supportive therapy, and / or other therapies known to one of skill in the art, such as medical personnel, useful in the management or treatment of a given disease, or symptom related thereto.

[0092] As used in the disclosure, “treat”, “treatment”, and “treating” refer to the reduction or amelioration of the progression, severity, and / or duration of a given disease resulting from the administration of one or more therapies (including, but not limited to, the administration of microspheres disclosed herein).

[0093] The present system described herein features a computational framework combining the C-HiDeNN theory with isogeometric analysis (IGA), called Convolution IGA (C-IGA). C- IGA of the present invention has five key features that advance IGA: 1) arbitrarily high-order smoothness and convergence rates without increasing degrees of freedom; 2) a Kronecker delta property that enables direct imposition of Dirichlet boundary conditions; 3) automatic and flexible global / local mesh-adaptivity with built-in length scale control and adjustable radial basis functions; 4) ability to handle irregular meshes and triangular / tetrahedral elements; and 5) GPU implementation that speeds up the program as fast as finite element method (FEM).

[0094] Mathematically, the present invention discloses that both IGA and C-IGA mappingsAttorney Docket No.0116936.296WO2 are equivalent, and by taking a special design and modified anchors as nodes, C-IGA degenerates to IGA. The present invention demonstrates the accuracy, convergence rates, mesh-adaptivity, and performance of C-IGA with several 1D, 2D, and 3D numerical examples. The present invention also discloses applications of C-IGA from topology optimization to product manufacturing with multi-GPU programming are discussed.

[0095] Fig.1 shows the main components of the C-HiDeNN-AI.

[0096] As shown in Fig.1, the first component of C-HiDeNN-AI is multifidelity data. This data is required for an AI to learn the appropriate physics of the system. For example, in the case of designing a feedback control system for metal additive manufacturing, one needs to monitor the temperature data coming from the specimen or the melt pool. However, the data can come from multiple sources and can be used to train the entire network or some part of it while the rest of the network solves partial differential equations.

[0097] The second building block is discrete calculus. The discrete calculus formulations, especially kernel and convolution-based equations, form the basis of the HiDeNN solvers. Kernel-based operator learning methods and deep learning interpretation of discrete calculus algorithms have demonstrated that neural networks can essentially degenerate into integral equations with each layer representing one particular integral transform. This approach is already established in the literature of partial differential equation solvers. Based on this method, anyfunctional relationship between an input ^^( ^^) and an u( ^^) can be expressed asu( ^^) = ∫ ^^( ^^, ^^′) ^^( ^^) ^^ ^^′ (1)

[0098] where the vector ^^vector and P is the total number of parameters. ^^(⋅,⋅) is the kernel function for interpolation. The HiDeNN solvers convert this kernel into custom neural networks. The parameter vector can include spatial variables such as p1 =x, p2 =y, p3 =z; temporal variables such as p4 =t; and system parameters such as scan speed for an additive manufacturing process p5 = v. For large-scale engineering problems, p can be extended to even higher dimensions. To counter the curse of dimensionality and improve the speed of calculations, the kernel can be split into multiple directionalconvolution kernels by the following formulation (more details will be in the next sections):u( ^^) = ∫ ^^ ′ ′ ′1( ^^1, ^^1 ) ∫ ^^2( ^^2, ^^2 ) … ∫ c ^^ ^^( ^^ ^^ , ^^ ^^ ) ^^( ^^) ^^ ^^′ ′ ′1 ^^ ^^2 .. d ^^ ^^ (2)Attorney Docket No.0116936.296WO2

[0099] where elements ofpare independent of each other. Based on this idea, a family of HiDeNN solvers has been developed over recent years (see, the third building block in Fig.1). One feature of these HiDeNN solvers is that they are designed to leverage modern hardware such as a Graphics Processing Unit (GPU) and a Tensor Processing Unit (TPU). The TPUs are specifically designed to train sufficiently deep neural networks. The C-HiDeNN-AI has already demonstrated its superior efficiency over serial programming when implemented with a single GPU. However, multi-GPU programming has become a standard for extremely large-scale scientific computing and generative AI training. For example, ChatGPT is known to be trained with over 10K NVIDIA H100 GPUs. Therefore, the development of C-HiDeNN-AI software that can fully leverage the unprecedentedly fast-growing GPU technology is now more essential than ever. For the convenience of readers, a summary of this evolution is given below and shown Table 1.

[0100] HiDeNN: HiDeNN consists of a custom neural network where a group of neurons from each layer mimics the interpolation of finite element shape functions that can be either trained or solved. If a HiDeNN is trained, a single dataset of independent variables is required and nodal adaptivity (r-adaptivity in FEM) is achieved through training by minimizing the loss function. When solved, no initial dataset is required. The HiDeNN minimizes the energy of the physical system with nodal adaptivity. Applications of HiDeNN-FEM to nonlinear problems have been discussed in.

[0101] HiDeNN-AI: Saha et al. proposed an artificial intelligence system that can enrich the HiDeNN with computational or experimental data and perform a parametric study. The article proposed three general classes of problems: type 1 for which a lot of data is available with little physical insight; type 2 where some physics is missing that can be recovered from data, and type 3 where a fast reduced order model is required to solve large scale problems.

[0102] HiDeNN-Tensor Decomposition (TD): To improve the speed of the HiDeNN, a tensor decom-position method is introduced where the kernels of equation (2) are approximated by tensor decomposition of the space-time-parameter matrix and superposition of the modes.

[0103] Convolution HiDeNN (C-HiDeNN): The HiDeNN method is extended to convolution HiDeNN (C-HiDeNN) to increase the smoothness and accuracy of the solutionAttorney Docket No.0116936.296WO2 without increasing the degrees of freedom (DoFs) to solve. This smoothness is achieved by designing an interpolation that accounts for additional layers of neighborhood elements.

[0104] Convolution HiDeNN (C-HiDeNN)-Tensor Decomposition (TD): To improve the speed of C-HiDeNN and solve problems spanning multiple length scales, C-HiDeNN-TD is developed. With this method, an unprecedented speed up is achieved for multi-scale and multi- material topology optimization. Combining topology optimization with multi-scale structures holds promise for multi-functional design with superior performance.

[0105] The present invention extends to Isogeometric Convolution HiDeNN (C-IGA). Compared to the previous methods, in which integrates spatial, temporal, and parameter spaces, the present invention focuses on geometric accuracy, which is an important parameter for manufacturing. Modern 3D printers can materialize virtually any design if supplied with a Computer-Aided Design (CAD) geometry. The quality of the final build depends largely on how smooth and geomterically accurate the supplied CAD geometry is. Isogeometric analysis (IGA) can capture the complex geometry of a problem, which inspired the current development of C- IGA of the present invention that integrates isogeometric shape functions with C-HiDeNN. With the flexibility of the shape functions used in C-IGA, the method of present invention achieves more versatile adaptivity in the solution, which is currently not available in FEM or IGA.

[0106] Table 1. A Summary of the HiDeNN Solvers based on Eq. (2). Parameters Physical Quantity Method(s) that Consider the Parameter in the Kernel P1Spatial X-direction HiDeNN, HiDeNN-AI, HiDeNN-TD, C-HiDeNN, C-HiDeNN-TD P2Spatial Y-direction HiDeNN, HiDeNN-AI, HiDeNN-TD, C-HiDeNN, C-HiDeNN-TD P3Spatial Z-directionHiDeNN, HiDeNN-AI, HiDeNN-TD,C-HiDeNN, C-HiDeNN-TD P4 Time HiDeNN-AI, HiDeNN-TD, C-HiDeNN, C-HiDeNN-TD P5 Manufacturing Parameter(s)HiDeNN-AI, HiDeTNDN-TD, C-HiDeNN-Attorney Docket No.0116936.296WO2 P6 Physical Parameter(s)HiDeNN-AI, HiDeTNDN-TD, C-HiDeNN-P7GeomCeotrnidcitBioonusndaryC-IGA

[0107] The HiDeNN-TD solvers use Eq. (2) to counter the curse of dimensionality. In software 1.0, for three space dimensions, one time dimension, and (P − 4) parameters where P ≥ 4, the total number of degrees of freedom (DoFs) increase exponentially with problem dimensions: (n3× n × n(P−4)) ≈ O(nP), assuming n DoFs for each dimension. In HiDeNN-TD, with Q modes, the total DoFs linearly: Q×(n × 3 + n + n × (P − 4)) ≈ O(Q × n ×P).

[0108] In all of these solvers, the physical equations, such as the minimum of potential energy (when it exists), are converted to an optimization problem via loss function and solved for the solutions. This loss function can be designed to match a desired control response from the system in real-time. Thus, the C-HiDeNN-AI aims to combine all the elements shown in Fig.1 to produce a software environment that can reduce the design to the final product cycle in the manufacturing plant with orders of magnitude speed up.

[0109] The final building block in Fig.1 is the real-time control. To ensure real-time control, one needs a fast, reliable, and accurate solver that can be integrated directly in manufacturing hardware. The proposed C-HiDeNN-AI provides the flexibility of training a neural network from data while keeping the structure of the network similar as solving. By generating the data using the solver, and training the network with the data, one can achieve the real-time control of manufacturing system for better product.

[0110] All mathematical notations and symbols are defined in Table 2 below. Table 2. Mathematical notations and symbols Symbol Description Ω, Ω , ΩPhysical, parametric, parent domain ^^, ^^ ^^Global nodal coordinates and nodal coordinates of node i in physical domain ^^, ^^ ^^Global nodal coordinates and nodal coordinates of node $i$ in parametric domainAttorney Docket No.0116936.296WO2 d Dimension of problem s, a, p Patch size, dilation parameter, polynomial order ^^^^Set of nodes in element e ^^^^^^ ^^, ^^^^^^^^ Set of nodes in a patch domain of node i for C-HiDeNN and C-IGA ^^ ^^ ^^Set of nodes in a patch domain of element e, ^^^^^^= ⋃^^∈ ^^^^^^^^^^ ^^^^^^^^^^Gauss point of element e ^^ ^^ ^^ , ^^ ^^ ^^Gauss points and weights in parent domain ^^, ^^^^Global nodal variables and nodal variable of node i ^^ ^^ ^^Column vector of nodal variables corresponding to ^^^^^^ ^^or ^^ ^^^^^^^^ ^^( ^^ ), ^^ ^^( ^^)Compact supported interpolation function at node i ^^^^^^, ^^( ^^ ), ^^^^^^,^^^^ ( ^^) Convolution patch function of node j, constructed for ^^^^^^^^or ^^^^^^^^with a dilation parameter of a ^^^^^^, ^^( ^^ ), ^^^^^^^,^^^( ^^) Vector of convolution patch functions constructed for ^^^^^^ ^^or ^^^^^^ ^^with a dilation parameter of a ^^ ^^ ( ^^ ), ^^ ^^(^^) C-HiDeNN and C-IGA shape functions ^^( ^^ ), ^^( ^^)Vector of basis functions to be reproduced ^^ ^^( ^^ ), Ψ ^^, ^^( ^^ )Vector of radial basis functions and its one component for node j ^^( ^^ ), ^^( ^^)Jacobian matrix for physical domain in C-HiDeNN and C-IGA ^^ ( ^^ )Jacobian matrix for parametric domain in C-IGA Π, Π^^Variational energy of global domain and element e ^^Elasticity matrix ^^^^^^Strain-displacement matrix of elemental patch ^^^^^^^^^^^^C-HiDeNN or C-IGA shape functions for elemental patch ^^^^^^^^^^( ^^), ^^^^Basis functions and coefficients for general geometric mapping ^^ ^^ Number of global nodes ^^ ^^Number of nodes in the nodal patch ^^ ^^^^^^ℎ, ℎ^^Element size in physical, parametric domain ^^^^ ^^ ^^, ^^ℎExact solution, numerical solution with a element size hAttorney Docket No.0116936.296WO2 C-HiDeNN framework C-HiDeNN approximation

[0111] C-HiDeNN borrows the same concept as HiDeNN that constructs a general approximation using neural networks with predefined coefficients and activation functions, but it has one additional convolution layer to enable nonlinear activation and versatile mesh adaptivity. Considering an element-wise approximation, the C-HiDeNN approximation can be written as ^^^^− ^^ ^^ ^^ ^^ ^^ ^^(ξ ) =∑^^∈ ^^^^^^ ^^(ξ )∑^^∈ ^^ξξ^^ ^^ (ξ ) ^^ ^^ =^^ ^^, ^^∑^^∈ ^^^^^^^^ ^^ (ξ ) ^^ ^^ . (3)

[0112] the problem and thephysical domain can be define as isoparametric mapping, as illustrated in Fig.2. Ni( ^^ ) is a function with compact support associated with node ^^ i of an element, and ^^ξ^^, ^^(ξ ) is a local convolution patch function with a dilation parameter a that interpolatesξI(As, called nodal patch) of the node ξı. For example, in Fig.2, the element ^^ has four{ξ1, ξ2, ξ3, ξ4}. The function ^^1) is the interpolant of node ξ1among the four nodes. Thepatch domain of nodeξ1,isAξs1, where ^^ is the size of the nodalindicates the number of layers of nodes surrounding the center node. In general, ^^ can be different in different directions (see Fig.2 where s =[3,2]in two different directions respectively), similarly for the dilation parameter a = [aξ, aη]. The present invention restricts the discussion to a uniform case with the same scalar ^^ and ^^ for all directions. The elemental patch domain of element ^^ is defined by Aes= ⋃i∈AeAξsı(blue dashed box in Fig.2). The final approximation can be written as a product of the convolution shape functions and the nodal coefficients in the elemental patch.

[0113] The convolution patch function ^^ξ^^,^^^^plays a similar role as the kernel function in convolution neural networks, which connects the neighboring nodes and is expected to improve the solution smoothness. The function ^^ ^^(ξ ) provides a compact support for the final interpolantNk (ξ ). Various choices of functions^^ ^^ξ ^^^^, ^^ have been discussed. In general, ^^ ^^can be anyfunction with compact support, such as finite element and isogeometric basis functions.Attorney Docket No.0116936.296WO2 Meshfree type approximations can be used for the local convolution patch function ^^ξ^^,^^^^. The present invention designs the C-HiDeNN approximation to satisfy the Kronecker and reproducing properties.

[0114] Among all possible combinations of the two interpolants ( ^^^^and ^^^^^^,^^^^ ) of C- HiDeNN, the present invention considera a special case, C-FEM, where FEM shape functions are used for ^^^^^^,^^^^ . The radial basis function is chosen for ^^^^^^,^^^^ , defined by ^^ξ^^, ^^= Ψa(ξ )αj+ P(4) where Ψa(ξ ) is a vectorwith a dilation of ^^. Here, ^^ normalizes the radial distance of interest ξ and the center node ξı. For example,Ψa,j(ξ ) can be written asΨa,j(ξ ) = ψ(|ξ − ξ^| / a) (5) with a cubic spline function or a Gaussianψ(z) = e−z2,

[0115] Where ^^ = | ^^ − ^^^^| / ^^. P(ξ ) is a vector of functions that are expected to be reproduced, such aspolynomials. Αjand κjare coefficients.

[0116] The neural network structure of C-HiDeNN approximation is shown in Fig.3. The solution field uC−HiDeNNis a function of nodal solutions ^^, nodal coordinates ^^ , and variable parameters a, s, p.

[0117] The nodal solutions of the partial differential equation can be found by solving the following optimization problem:Attorney Docket No.0116936.296WO2

[0118] where be viewed as anoperator that computes the variational energy. Here, the present invention fixes all parametersX, a, s, p except the nodal solution U. This is mathematically equivalent to solving the variationalequation like standard FEM using the matrix form. Train (9)

[0119] wherethe set of positive real numbers and natural numbers, respectively. The optimization of X is referred to r-adaptivity in FEM literature, but the other parameters do not have counterparts in standard FEM. There exists a difference in the optimization of parameters a, s, p: they are optimized patch by patch, while nodal coordinates X are optimized node by node. That is to say, all convolution patch functions ^^ξ^^,^^^^in the nodal patch AξsIalways share the same parameters a, s, p in the optimization, and these parameters can be different in different nodal patches. In addition, unlike most adaptive FEM exploits the a posteriori error estimator, HiDeNN theory only utilizes automatic differentiation and objective function (i.e., variational energy) for various adaptivity. Discussions on parameters: p, s and a

[0120] There are three controlling parameters (p,s,a) in C-HiDeNN interpolation, which can be optimized to enhance solution efficiency for a specific problem.

[0121] Positive integer p refers to the order of P(ξ ) in the convolution patch functions in Eq. (4), which directly relates to the order of reproducing properties. The reproducing property will directly affect the convergence rate.

[0122] Positive integer s is the size of the nodal patch AξsIthat is used as an influence domain for computing the convolution patch functions Waξ,ıs . patch AξsIin d-dimensional space comprises(2s + 1)dnodes centered at node ξIif the mesh isthat one needs enough nodes e.g., at least (p+1) nodes along one dimension) to guarantee the p-th orderAttorney Docket No.0116936.296WO2 reproducing property, which leads to:s ≥ p / 2. Increasing s means enhancing the nodal connectivity, but this also leads to a larger bandwidth of the global stiffness matrix and computational cost.

[0123] Positive real number a is the radial dilation of the radial basis function Ψa(ξ ) in the convolution patch functions in Eq. (4). It controls the characteristic length scale of the radial function, and further influences the shape of the convolution patch function. One can adjust this parameter to change the shape and characteristic length scale of shape functions without changing the mesh and the reproducing property. The selection of a is independent of s and p.

[0124] The conditions for the three parameters are summarized in Table 3. Table 3. conditions for the three parameters a, s, p Parameter p s a Condition^^ ∈ ^^+^^ ≥ ^^ / 2, ^^ ∈ ^^+^^ ∈ ^^+Example 1: Single-Patch C-IGA 1. Isogeometric C-HiDeNN (C-IGA) with two mappings 1.1 Isoparametric mapping with exact geometry

[0125] Like traditional FEM, the FE mesh used in C-HiDeNN approximation does not guarantee the geometric accuracy in the numerical analysis. The present invention borrows the idea of IGA and propose the so-called isogeometric C-HiDeNN, shortly C-IGA, to link directly from CAD (Computer Aided Design) models to numerical analysis. Compared to C-HiDeNN, there is an additional intermediate domain called the parametric domain between the physical and parent domains in C-IGA, which is similar to that of IGA. The three domains of IGA and C- IGA are illustrated in Fig.4 (IGA) and Fig.5 (C-IGA). The domains of element e in the physical, parametric, and parent domains are denoted by Ωe, Ωe, Ωerespectively, and corresponding coordinate systems are denoted by x, ξ, and ξ , respectively.

[0126] The key ideas of C-IGA are: 1) construct the isoparametric C-HiDeNN interpolation in the parametric domain, and 2) reproduce the geometric mapping xGeo(ξ) from CAD / IGA by mapping discrete nodal points through the geometric .Attorney Docket No.0116936.296WO2

[0127] Eqs. (10) and (11) shown below are the isoparametric C-HiDeNN interpolation in the parametric domain, which are the same as those of C-HiDeNN, except that the parent coordinateξ is changed to the parametric coordinates ξ.xC−IGA(ξ) = ∑i∈Ae ^^ ^^(ξ) ∑ ξ ^^ ej∈Aξis ^^ ^^, ^^ (ξ)xj = ∑k∈As Nk (ξ)xk (10)and^^^^− ^^ ^^ ^^(ξ) =^^∈ ^^^^ ^^^^(ξ)^^∈ ^^ξ ^^^^^^^^^^, ^^(ξ) ^^^^=^^∈ ^^^^^^^^^^(ξ) ^^^^(11)

[0128]

[0129] present equals the geometric mapping xGeo(ξ) from CAD / IGA. To achieve this, there are two points to consider:

[0130] The physical coordinates xk should come from the geometric mapping, xk =xGeo(ξk), with given parametric coordinates ξk in the parametric domain.

[0131] Substituting xk = xGeo(ξk) into Eq. (10), one can rewrite the coordinate mapping as:xC−IGA(ξ) = ∑k∈Aes Nk (ξ)xGeo(ξk) (12).

[0132] The C-IGA shape functions Nk(ξ)should be able to reproduce the basis functions ofthe geometric mappingxGeo(ξ). In general, the IGA geometric mapping is written in the form:xGeo(ξ) = ∑l Rl(ξ)Bl (13)

[0133] where Rl(ξ) are basis functions, and Blcoefficients. For example, Rl(ξ) are B-spline or NURBS basis functions, and Blare control points in the IGA mapping xIGA(ξ) as illustrated in Fig.4.

[0134] To reproducexGeo(ξ), the present invention starts from the reproducing condition forRl(ξ), that is,∑k∈AesNk(ξ)Rl(ξk)= Rl(ξ)(14)

[0135] This property of C-IGA interpolation is realized by sepcifically designing the polynomial basis P(ξ) in Eq. (4).Attorney Docket No.0116936.296WO2

[0136] Since IGA can use either B-spline basis or NURBS basis function for xIGA(ξ), there are two choices for P(ξ): 1) general polynomial basis functions, or 2) polynomials over the weight function (weighted polynomials).

[0137] Based on the above two points, one can prove that the C-IGA mapping is precisely the xGeo(ξ) and can reproduce the exact geometry of CAD models.

[0138] Similar to IGA, for the convenience of numerical integration, the present invention needs a second mapping from the parent domain to the parametric domain in one element, defined by ξ(ξ ) =∑i∈AeNi(ξ )ξi(15)

[0139] The DNN form of C-in Fig.6. Like C-HiDeNN, the elemental neural network of C-IGA as shown in Fig.6 panel(b) also has multiple neurons in the third hidden layer that connect neighboring nodes in the nodal patch domain AξsIwith convolution patch functions ^^^ξ^,^^^^as weights. The difference between C-HiDeNN and C-IGA in the neural network form isIGA has two mappings (the geometric mapping and the FEM mapping) and the main interpolation takes place in the parametric domain.

[0140] With the above two mappings, the variational energy for an elastic system in one element e is given:

[0141] force. Useis the unknown nodal displacement vector with a length of nis the strain-displacement matrix whose size is [3, n(Aes) × d] for a 2D element. D is the modulus matrix.}Attorney Docket No.0116936.296WO2

[0142] Unlike C-HiDeNN, C-IGA has two Jacobian matrices corresponding to the two mappings. J is for the first mapping defined by

[0143] Because of the identity ofGeox(ξ)= x(ξ)and the C-IGA mapping in Eq. (12), J(ξ) can be obtained from either of them. The second Jacobian matrix J corresponds to the second mapping defined in Eq. (15):

[0144] 1.2 Proof of exact geometry

[0145] The present invention assumes that the geometric mapping consists of basis functionsRl(ξ) and coefficients Bl as shown in Eq. (13).

[0146] By designingP(ξ), e.g., taking P(ξ) = Rl(ξ), the convolution patch functions in the nodal patch can reproduce each basis Rl(ξ) exactly, i.e., ∑ξ ^^j∈Aξis^^^^, ^^(ξ)Rl(ξj) = Rl(ξ) (19)

[0147] Combining with the partition of unity for ^^^^(ξ), C-IGA interpolations have the same reproducing property Eq. (14). The proof is ∑k Nk(ξ)Rl(ξk)=∑i ^^ ^^(ξ) ∑ξi ^^ξ^^,^^^^(ξ)Rl(ξj) = Rl(ξ) ∑i ^^ ^^(ξ)= Rl(ξ)(20)

[0148] C-IGA mapping (Eq. (12)) and the geometric mapping xGeo(ξ):Attorney Docket No.0116936.296WO2

[0149] Remark that the present invention mainly talks about B-spline and NURBS basis functions in this paper. The proof is still valid for T-spline if the convolution patch functions in the nodal patch can reproduce T-spline basis functions. Complex CAD geometries with extraordinary nodes and arbitrary genus number can also be handled by C-HiDeNN / C-IGA interpolation (Eq.(3)). The C-HiDeNN / C-IGA formulation uses both finite element shape functions and meshfree type basis functions, which are not influenced by extraordinary nodes. More details on the implementation will be investigated in our future work. 1.3 Degeneration of C-IGA to IGA

[0150] Although C-IGA can reproduce the IGA mapping, the present C-IGA can not degenerate to IGA directly, as C-IGA has less DoFs that IGA. However, if one takes the modified anchors {ti} of IGA as the supporting nodes instead of non-repeating knots {ξi} to construct C-IGA interpolation, and use IGA basis functions as the convolution patch functions,the new C-IGA will degenerate to IGA (they are equal in fact).uC−IGA(ξ) = ∑i Ni(ξ; t∗i) ∑j∈Atis Nj,p(ξ)dj (22)

[0151] where t∗∗i function NI, such as ti =[ti−1, tI, ti+1] in 1D linear case. The nodal patch AtIs is the set of nodes (anchors) in a patchdomain of node (anchor) ti. Nj,p(ξ)is the p-th order basis function for node (anchor) tj, and djis the corresponding nodal coefficient. Note that djis not a nodal value of the solution like uj, for there is no Kronecker delta property in the new formulation. 1.4 ^^-adaptivityAttorney Docket No.0116936.296WO2

[0152] In IGA, the multivariate B-spline and NURBS basis functions are constructed by tensor product, which requires a regular mesh in the parametric domain. The T-spline enables local mesh refinement but it is still limited by regular meshes; IGA cannot handle arbitrary irregular meshes. On the other hand, both C-IGA and C-HiDeNN allow an irregular mesh because they can generate any type of mesh in the parametric (C-IGA) or parent (C-HiDeNN) domain. The ability to handle an irregular mesh is one of the key advantages of C-IGA over IGA, as illustrated in Fig.7. r-adaptivity is achieved by minimizing the variational energy with respect to the nodal parametric coordinates, that is,

[0153] Note that χ = [ξ1, ξ2, … , ξnp] denotes a vector of all nodal parametric coordinates with np number of globalNodal solution U is obtained after given χ, a, s, p, which can be interpreted that U is the function of χ in the optimization process. The present invention limits the movement of boundary nodes to keep the exact geometry unchanged. The difference from the r-adaptivity of C-HiDeNN is to optimize nodal parametric coordinates instead of nodal physical coordinates. The physical coordinates xkwill also follow the migration of parametric coordinates ξkaccording to the mappingxk= xGeo(ξk).

[0154] There are several points to be considered here.

[0155] In one embodiment, one can leverage automatic gradient functionalities for gradient- based optimization in many popular scientific computation packages such as PyTorch and JAX.

[0156] In one embodiment, the convexity of optimization is not guaranteed, meaning that one might not be able to find the optimal nodal coordinates even after prolonged iteration. Furthermore, it may result in an unphysical mesh having a negative Jacobian determinant. A careful selection of learning rate and early stopping criteria is required.

[0157] In one embodiment, the global mesh update may suffer from an immense cost of optimization. A more efficient strategy might be to update locally. In-depth studies on this local refinement are needed.Attorney Docket No.0116936.296WO2 1.5 ^^-adaptivity

[0158] In addition to r-adaptivity, both C-HiDeNN and C-IGA can adjust the dilation parameter a to get better accuracy. The so-called a-adaptivity refers to the dilation parameter a varying from one nodal patch to another nodal patch. As discussed previously, the parameter a controls the characteristic length scale of the radial function, and further influences the shape of the convolution patch function. For a complex solution, the optimal characteristic length should vary with spatial positions. The present invention seeks the optimal spatial distribution of the dilation parameter a by minimizing the variational energy:

[0159] Theby a = [aξ1, aξ2, … , aξnp] with aξIrepresenting the dilation parameter of convolution patch functions Aξsi. Unlike the r-adaptivity where the global nodal coordinates arebe performed patch-wise, rendering a wide choice of optimization strategies.

[0160] Unlike r-adaptivity which might result in unphysical meshes, a-adaptivity is free from the same issue because the adjustment of ^^ does not influence the mesh. The only constraint for ^^ is to be a positive real number, as given in Table 3. Additionally, although the present invention adopts an isoparametric mapping, the mapping is precisely the geometrical mapping xGeo(ξ) according to the equivalence Eq. (21). That is, the mapping is independent of the dilation parameter ^^. This leads to a simpler derivation of analytical derivatives with respect to ^^ in ^^-adaptivity compared to r-adaptivity. 1.6 ^^- ^^-adaptivity

[0161] According to Table 3, s ≥ p / 2. However, the preliminary numerical tests of the present invention suggest s ≥ p for good accuracy, efficiency, and numerical stability. Small s compared to p might cause stability issues when one imposes the Kronecker delta property on the convolution patch functions. A larger ^^ leads to a higher order convergence rate while a larger ^^ leads to a more accurate solution and a higher computational cost. Spatially varying ^^Attorney Docket No.0116936.296WO2 and ^^ can be considered similarly to ^^ -adaptivity, called s-p-adaptivity. That is, like ^^- adaptivity, the s-p-adaptivity can be achieved nodal patch-wise. For efficiency, in one embodiment, one can adopt a hybrid strategy, i.e., use large ^^ or ^^ for important regions and small ones for the other regions. This idea is verified in the following 1D numerical examples. Table 4: Comparison between FEM, IGA, C-HiDeNN, and C-IGA. FEM IGA C-HiDeNN C-IGA DoFs Attached to Attached to Attached to nodes Attached to nodes nodes control points Interpolation Interpolates Does not Interpolates nodal Interpolates nodal nodal points interpolate points points control points Reproducing Exact Order Exact Exact Order Exact geometry property geometry Adaptivity ^^−, ℎ−, ^^ − ℎ−, ^^ − ^^−, ℎ−, ^^−, ^^−, ^^ − ^^−, ℎ−, ^^−, ^^−, ^^ − adaptivity refinement, adaptivity. adaptivity. need to update Versatile Versatile control points 1.7 Discussion

[0162] C-IGA is an extension of C-HiDeNN that realizes exact geometry like IGA. The present invention compares different aspects of FEM, IGA, C-HiDeNN, and C-IGA in Table 4.

[0163] First, unlike IGA where the control points have DoFs (as known as control variables ), all FEM, C-HiDeNN, and C-IGA have DoFs attached to nodal points. This is the reason why C-IGA can handle an irregular mesh, as discussed in the Subsection 1.4 and Fig.7. The three methods (FEM, C-HiDeNN, and C-IGA) are also able to exactly interpolate nodal points because the convolution patch functions have a Kronecker delta property, whereas IGA cannot interpolate control points except for repeating knots. In terms of the reproducing property, FEM and C-HiDeNN can reproduce polynomial orders up to a given p because they use polynomial basis functions. On the other hand, both IGA and C-IGA use NURBS basis which are rational functions (neglecting the case when the NURBS weights are all equal to 1, i.e., B-splines with polynomials). With NURBS basis functions, they can reproduce exact geometry, but not the polynomials like FEM or C-HiDeNN. Additionally, C-HiDeNN uses a parent domain[−1,1]d(with some extended neighborhood depending on the patch size s) like FEM which focuses onAttorney Docket No.0116936.296WO2 one element, as shown in Fig.2. In contrast, C-IGA uses the parametric domain which is defined globally, as shown in Fig.5.

[0164] In terms of adaptivity, C-HiDeNN and C-IGA have versatile ^^−, ℎ−, ^^−, ^^−, ^^ −adaptivities in contrast with FEM or IGA. Detailed comparison of adaptivities between IGA and C-IGA can be found in Table 5. Some key points from the table are listed below.

[0165] In one embodiment, C-IGA follows the same framework of C-HiDeNN or FEM. There is no control net in the interpolations, although the geometric mapping may come from the IGA mapping with a coarse control net. B-spline, NURBS, and T-spline in IGA are all developed based on a regular mesh, since their basis functions are obtained by tensor product. Compared with IGA, C-IGA allows an irregular mesh and nodal adaptivity in all directions except for the boundary nodes. In short, the r-adaptivity of C-IGA can improve accuracy without additional DoFs.

[0166] In one embodiment, the r -adptivity of C-HiDeNN is realized in the physical domain. In contrast, C-IGA introduces an additional intermediate domain, parametric domain, to realize the exact geometry. The r -adaptivity is realized in such a parametric domain, and the physical coordinates move accordingly.

[0167] In one embodiment, C-IGA and IGA both have the parameter p to denote the order. The IGA can elevate the order p by p-refinement or k-refinement. The former adopts the order elevation element by element after knot insertion, which is analogous to p-refinement in FEM. The latter adopts the order elevation followed by the knot insertion. The $k$-refinement significantly lowers the DoFs increment compared with the p-refinement, but still requires some additional DoFs and basis functions. In contrast, the order elevation of C-IGA occurs in the convolution patch function, specifically with P(ξ) in Eq. (4). Hence, one can elevate the order p without changing the mesh and DoFs.

[0168] In one embodiment, C-IGA has two additional parameters compared with IGA: the patch size s and the dilation parameter a. Patch size s denotes the size of the patch domain involved in building shape functions. It controls the nodal connectivity and affects the bandwidth of the global stiffness matrix. As discussed previously, the dilation parameter a controls theAttorney Docket No.0116936.296WO2 characteristic length scale of the radial function. C-IGA allows the a-adaptivity to seek an optimal distribution of patch-wise a for improved solution accuracy without changing the mesh and DoFs. Table 5. Comparison of adaptivity between IGA and C-IGA. Parameters ICA C-IGA Mesh (r) - Automatic r-adaptivity due to HiDeNN Mesh (h) - Control net - No control net - Allows nonuniform but - Allows irregular mesh regular mesh * Order (p) - p-refinement knot -Elevates p without increasing insertion followed by DoFs. order elevation - k-refinement order elevation followed by knot insertion - Both cases need extra DoFs and recomputation of control points and knot vectors. Patch size (s) - Controls solution accuracy through modifying connectivity and bandwidth of the stiffness matrix Dilation (a)- ControlsCharacteristic length scales - Does not change the mesh * B-spline, NURBS, and T-spline are all developed based on the regular mesh, because their basis functions are obtained by tensor product. 2. Continuity and convergence rate 2.1 High-order continuityAttorney Docket No.0116936.296WO2

[0169] The present invention analyzes the continuity of shape functions Nk(ξ)of C- HiDeNN / C-IGA interpolations.

[0170] The C-HiDeNN / C-IGA convolution shape function Nk (ξ) is defined byN(ξ) ∑ ^^| ^^∈As ^^ ( ^^) ^^ ^^k ξi ^^ ^^ ^^, ^^ (25)

[0171] According to theof Nk(ξ) isofNk (ξ) depends on the continuities of ^^ ^^(ξ) and ^^ξ ^^^^, ^^ . Assuming ^^ ^^(ξ) is Cn continuous and^^ξ^^,^^^^is Cmcontinuous, DNNk(ξ) is continuous up to ^^ = min(n, m).

[0173] There exists a special case for 1D interpolation: Nk (ξ)has CN continuity with ^^ =min(n + 1, m), if the Dn+1 ^^ ^^(ξ) is discontinuous across nodes, e.g., FEM shape functions. Thatis due to the Kronecker delta property of ^^ξ^^,^^^^and partition of unity property of ^^^^(ξ). The first term in Eq. (26) at the left and right limits of ξjisterm and the final term.

[0175] If adopting the radial basis functions (Eq. (4)) as convolution patch function ^^^^^^,^^^^ , the continuity of ^^^^^^,^^^^ depends on two parts: the radial function ^^^^(^^)and reproduced basisAttorney Docket No.0116936.296WO2 functions ^^(^^). Generally, ^^(^^)is ^^∞, such as polynomials. ^^^^(^^)can also be ^^∞when taking the Gaussian radial functions Eq. (7). As a result, the overall continuity is depending on the selection of ^^^^( ^^) in general. Careful design of ^^^^( ^^) and ^^^^^^,^^^^ can lead to arbitrary high order ^^^^continuity of C-HiDeNN / C-IGA interpolants, which is difficult for traditional FEM. Note that in one embodiment, the present invention only requires partition of unity and compact support for ^^^^(^^). Thus B-spline basis functions and meshfree basis functions satisfying partition of unity can be good selections for ^^^^(ξ). 2.2 Theoretical convergence rate

[0176] If taking ^^(ξ)as p-th order polynomials in the reproducing conditions, C-HiDeNN / C- IGA interpolants satisfy the reproducing property as well, e.g., the following reproducing conditions for 1D p-th order polynomials ∑k∈AesNk(ξ)(ξk)q= ξq, ^^ = 0,1,2, ⋯ , ^^. (29)

[0177] Theinterpolants onthe upper bound of the error of C-HiDeNN for elliptic equations|uC-HiDeNN(x) − uExt(x)|L2 ≤ C p+1 Ext0 h |u |Hp+1 (30)and|uC-HiDeNN(x) − uExt(x)|H1 ≤ C1 hp|uExt|Hp+1 (31)

[0178] where C0, C1are constants independent of h, h is the mesh size, p is the reproduced polynomial order, and uExtdenotes the exact solution. It indicates that C-HiDeNN / C-IGA achieves the (p+1)-th order convergence rate in terms of L2-norm and the p-th order convergence rate in terms of H1-norm when the interpolation functions satisfy the p-th order reproducing conditions in Eq. (29). The present invention also verifies this conclusion for C-IGA in numerical examples, and further numerically study convergence rates with adaptivities of various parameters. 3. Numerical Examples 3.11D example and implementationsAttorney Docket No.0116936.296WO2

[0179] The present invention first uses a 1D example to illustrate C-IGA and study the convergence rate as well as all kinds of parameter adaptivity. For illustrative purposes, the present invention considers a rod fixed at both ends under a body force ^^(x) as shown in Fig.9. The strong form reads with Dirichlet boundaryu(0)=0, u(L)=0 (33)

[0180] Here, ^^(x)is the displacement field, ^^ the stiffness of the rod defined as E=175, A the section area set to be A=1, and L the length of the rod set to be L=10. All parameters are unitless. The body force ^^(x)is given by

[0181] 1D C-IGA interpolation is defined in the parametric domain:uh(ξ) = ∑npi=1 ^^ ^^(ξ) ∑ ξξ ^^ npi ^^ ^^, ^^ (ξ)uj = ∑k=1 Nk (ξ)uk (36)

[0183] inventionadopts the cubic spline function in Eq. (6) as the radial function, p-th order polynomials[1, ξ, ξ2, … , ξp] ^^ ^^ ^^(ξ) in the convolution patch function ^^ξ ^^^^, ^^ (Eq. (4)).

[0184] For IGA, the mapping is defined by a knot vector Ξ = [0,0,1,1] with two ends as control points B1= 0, B2= ^^. This gives the simplest geometric mapping, an affine, i.e.,Attorney Docket No.0116936.296WO2 xIGA(ξ) = ∑2l=1Rl(ξ)Bl= ^^ (37)withR1(ξ) = 1 − ^^, R2(ξ) = ^^ (38)

[0185] In one embodiment, the present invention adopts the same geometric mappingxGeo(ξ) = ^^ as the first mapping between the parametric and the physical domains. Note that itis also an isoparametric mappingxC−IGA(ξ) = ∑k Nk (ξ)xk = ∑k Nk (ξ) ⋅ ^^ξk (39)

[0186] by taking xk= xGeo(ξk)= ^^ξkand ^^(ξ)=[1, ξ, ξ2, … , ξp]in the reproducing condition. Since the C-HiDeNN interpolation reproduces p-th polynomials including thelinear basis, the present invention has the following equivalencexC−IGA(ξ) = ∑k Nk (ξ)xk = ∑k Nk (ξ) ⋅ ^^ξk = ^^ ⋅ ∑k Nk (ξ)ξk = ^^ ^^. (40)

[0187] The second mapping is also an affine in the element e, i.e.,

[0188]

[0189] with 1D interpolants uhdefined in Eq. (36).

[0190] According to above simple mappings Eqs. (40) and (41), the present invention has3.1.1 Convergence studyAttorney Docket No.0116936.296WO2

[0191] The present invention starts with a uniform mesh to study convergence rates. Theelement sizes in the physical domain and the parametric domain are denoted by h=L / (np-1) andhξ=1 / (np-1), respectively, where np-1 is the number of elements. The accuracy of numericalsolution uhis measured by the relative L2-norm error defined by hE^^ ^^ ^^ ^^ ^^ =|u −u xt|L2|uExt|L2(44)

[0192] Before theinvestigates the accuracy of C- IGA results with varying dilation parameter a and the patch size s=2,3,4,5 for cubic spline and Gaussian radial functions as shown in Fig.10. For the cubic spline radial function in Fig.10 panel (a), the error oscillates for low a, but converges to a specific value when a is large enough. The oscillation observed in small a is attributed to the piecewise form of the cubic spline function with a compact support size 2. The lowest error occurs at some specific values of a, but it requires optimization of a which is called a-adpativity. On the other hand, the error converges for large a with moderate accuracy. Note that the value a at which the error converges increases as the patch size s enlarges. When s=2, the error converges for a / hξ> 7 while for s=5, the error converges for a / hξ> 19. Therefore, when one uses the cubic spline radial function, it is recommended to use a large a if one wants to skip the a-adpativity.

[0193] Fig.10 panel(b) shows the case when the Gaussian radial basis is used. Unlike the Cubic spline, the error curves diverge for large a and there is a global optimum for a that yields considerably accurate results. In general, for those optimal a, the Gaussian function is more accurate than the cubic spline. Therefore, when the Gaussian radial function is used, one needs to optimize a.

[0194] Based on the above discussion, the present invention chooses the optimal dilation parameter a for different patch sizes to study their ^^2-norm error convergence rates and compare the C-IGA solutions with linear FEM and IGA results. Here, the cubic spline radial function is used. The results are shown in Fig.11. As expected, the ^^2-norm error of linear FEM solution has a theoretical convergence rate of 2, and the p-th order IGA solution reaches the (p+1)-th order of convergence. For p=2,3, the convergence rates of C-IGA are asymptotically p+1 but the solution accuracy is better than linear FEM and IGA. For p=4, C-IGA (s=4, a / hξ= 9.9) almostAttorney Docket No.0116936.296WO2 coincides with IGA results, and the other solutions also reach the fifth-order convergence rate. These observations verify the theoretical analysis of the present invention. When optimal dilation parameters are used, the patch size s has less influence on the convergence and accuracy. However, large s broadens the bandwidth of the stiffness matrix, which leads to more computational costs both in time and memory. Hence, small s, e.g., s=p, might be a more effective option. Furthermore, the present invention discloses that both FEM and IGA need additional DoFs when the order p is elevated, although such DoFs increases for IGA are relatively small. In contrast, the C-IGA can achieve high order convergence rates by simply adjusting the polynomial order in the convolution patch function without increasing the number of nodes and DoFs.

[0195] In one embodiment, the present invention further performs a numerical study on the condition number of the stiffness matrix. The numerical results demonstrate that the condition number of C-IGA for different solution orders (p=2,3,4) is similar to that of linear FEM, slightly larger than IGA for orders (p=2,3,4), and significantly lower than the quadratic FEM. For fine mesh, C-IGA gives the best accuracy at a given condition number compared to IGA and FEM. 3.1.2 r-adaptivity

[0196] For the r-adaptivity, C-IGA can adjust nodal coordinates in the parametric domain then the physical mesh is obtained by the geometric mapping xGeo(ξ). To avoid the singularity due to extremely small element size, the present invention limits the element size to be bigger than hξ / 5 where ℎ^^is the original element size in the parametric domain.

[0197] two kinds of optimization algorithms: gradient-based methods and global optimization methods. Here, since the problem size is moderate and as one aims for the global optimum, the present invention adopts one of the global optimization methods, particle swarm optimization (PSO), to optimize nodal positions except for the boundary nodes. In one embodiment, the present invention showed that PSO achieves more accurate results than the Adam optimizer since PSO is designed to search for the global minima from many different random starting points.

[0198] Fig.12 illustrates the final optimized coarse mesh of C-IGA (s=2, p=2, a / hξ= 1.8) for h=L / 40 and h=L / 80. As seen in the figure, the nodesof two Gaussian humps, since the solution curve at these places needs more nodes compared to otherAttorney Docket No.0116936.296WO2 flat parts. The L2-norm error is plotted in Fig.13. The r-adaptivity improves the accuracy considerably for coarse mesh. C-IGA errors with r-adaptivity are merely 0.5-5% of the original C-IGA errors for meshes h=L / 40 and L / 80. While for fine mesh, the improvement from r- adaptivity becomes small because there are already enough nodes to capture the solution exactly. 3.1.3 a-adaptivity

[0199] The C-IGA has an additional parameter a, called the dilation parameter, compared to FEM and IGA. The present invention has observed a significant effect of a on the accuracy in previous discussions. However, the previous analyses are based on a constant a, i.e., all nodal patches share the same a. Here, a-adaptivity allows different aξIfor each nodal patch AξsI, and aims to seek for the optimal distribution of aξIin the global domain. Again, the present invention adopts the PSO algorithm because of the oscillatory nature of the cubic spline radial function with small a as shown in Fig.10 panel (a), which brings a big challenge for gradient-based optimizers; Gradient-based methods are more suitable for Gaussian radial functions. One of the advantages of a-adaptivity over the r-adaptivity is that one does not change the mesh during the optimization, and avoid the mesh distortion issue of the r-adaptivity. Here, for numerical convenience, the present invention limits aξi / hξ∈(1,20]. The numerical results from a- adaptivity are shown in Fig.14. The a-adaptivity helps this method achieve better results than those for the best constant a, and also better than IGA results. 3.1.4 p-adaptivity

[0200] Larger p results in a faster convergence rate and better accuracy, but this also leads to additional computational costs. To balance accuracy and efficiency, the present invention expects to use high-order interpolants in high-gradient regions such as the two Gaussian humps in Fig.8, and low-order interpolants in the remaining regions. This is the spirit of p-adaptivity. The flexibility of the C-HiDeNN framework allows such operation by simply elevating the order p of the convolution patch functions corresponding to the nodes in the high-gradient regions without adding nodes or changing mesh. Note that the patch size s should be adjusted correspondingly to guarantee p-th order reproducing property according to the requirements in Table 3. For illustrative purposes, the computational domain illustrated in Fig.15 panel (a) is divided into two parts: the high-gradient region where the two Gaussian humps are located, and the remaining regions for the flat part of the solution. The first region is covered by the blue boxAttorney Docket No.0116936.296WO2 ( ^^ ∈[1.5,3.5]∪[6.5,8.5]) in the figure. The present invention uses the fourth order C-IGA interpolants (s=4, p=4, a / hξ= 9.9) in the first region, and the second order interpolants (s=2, p=2, a / hξ= 1.8) for other regions. As shown in Fig.15 panel (b), one can observe that the numerical results with p-adaptivity reach the same order of accuracy as the complete fourth- order interpolants. This verifies the efficiency of p-adaptivity. 3.2 2D Poisson's problem on a quarter-ring domain

[0201] Starting from 2D examples, the present invention uses GPU programming to speed up both the construction of convolution patch functions and the iterative solver. The inventors have already demonstrated orders of magnitude speed up with GPU programming. Note that JAX is the high-level Python library for GPU programming developed by Google. As C-IGA formulation is similar to C-HiDeNN, the parallel program of C-IGA is developed according to that of C-HiDeNN. The present invention uses an NVIDIA A6000 GPU with 48 GB of GPU memory. 3.2.1 Convergence study

[0202] In one embodiment, the present invention considers a quarter ring as illustrated in Fig.16. The ring is the region between two concentric circles whose inner and outer radii are 10 and 20 (unitless), respectively. The present invention solves a Poisson’s problem given by: Δ ^^(x, y) + ^^(x, y) = 0, (x, y) ∈ Ω (45)

[0203] Let the body force b(x,y) be centered at (5,15):

[0204] uExt(x, y) = e−π(x−5)2e−π(y−15)2(47)

[0205] The

[0206] which is almostAttorney Docket No.0116936.296WO2

[0207] In one embodiment, the present invention uses an IGA mapping as the geometricmapping xGeo(ξ), given byxGeo(ξ) = ∑3i=1 ∑3j=1 Ri,j(ξ, η)Bi,j (49)where RI,j(ξ, η)areThe IGA geometric mapping from the parametric domain to the physical domain is shown in Fig.16 panel (a) with control points. Fig.16 panel (b) illustrates how the nodes and elements defined in the parametric domain are mapped into the physical domain using the IGA mapping. It is important to note that once the nodes in the parametric domain are mapped into the physical domain with IGA mapping(xk= xGeo(ξk)), the field mapping constructed by C-IGA ( ^^(ξ) = ∑kNk(ξ)ξk) becomes exactly as the IGA mapping according to the proof in Subsection 1.2. Thisrenders the isoparametric feature of C-IGA ( ^^(ξ) = ∑kNk(ξ)uk). The present invention defines the relative energy-norm error to measure the accuracy:mesh based on the nodes and elements in the physical domain. Simply put, the nodes in the physical domain are linearly connected like standard FEM meshes as illustrated in Fig.17 panel (a). On the other hand, in Fig. 17 panel (b), one can see the equivalence between IGA and C-IGA mappings. The difference takes place when the order p is increased: IGA has additional repeated knots at the boundary with more basis functions and DoFs while C-IGA has no repeated nodes and keeps DoFs unchanged.

[0209] Fig.18 illustrates the C-IGA solution (p=2) compared with the analytical one uExt. The results of convergence studies are plotted in Fig.19 for cubic spline radial functions and Fig. 20 for Gaussian radian functions. In those figures, the present invention compares the energy- norm error of FEM, IGA, and C-IGA with different s and fixed p. Fig.19 panels (a,b,c) and Fig. 20 panels (a,b,c) are for p=2 and Fig.19 panels (d,e,f) are for p=3. Note that the largest C-IGA problem has around 4 million DoFs.Attorney Docket No.0116936.296WO2

[0210] In Fig.19 panel (a), IGA (p=2), FEM (p=2), and C-IGA (s=2, p=2) have almost the same error for a given mesh size. However, as s increases, C-IGA shows faster convergence around hξ= 10−2, thus resulting in better accuracy for finer meshes. Note that the convergence rates of all cases are the same as 2, except for the linear FEM (p=1). Although the FEM (p=2) has a similar error compared to IGA (p=2) and C-IGA (s=2, p=2) at the given mesh size, it requires more DoFs than both IGA and C-IGA because of the higher order element. As shown in Fig.19 panel(b), to achieve 10−4error, C-IGA (s=5, p=2) only requires around 6 × 104DoFs while IGA (p=2) and FEM (p=2) needs around 2 × 106and 6 × 106DoFs, respectively, that is, around 30 times and 100 times more DoFs than C-IGA. The accuracy of C-IGA improves when one uses the Gaussian radial function with a properly chosen dilation parameter. As shown in Fig.20 panel(b), IGA (p=2) and FEM (p=2) require around 60 and 200 times more DoFs than C-IGA (s=5, p=2).

[0211] General trends of p=3 shown in Fig.19 panel (d,e) are similar to those of p=2. C- IGA accuracy improves as s goes up. Although it seems FEM (p=3) has the lowest error among all others given the mesh size, it requires much larger DoFs (almost 5 times more than 4-node elements) because of the 12-node serendipity element. In Fig.19 panel (e), C-IGA with higher s still gives the lowest error given the same DoFs.

[0212] Computation time is plotted against the error in Fig.19 panel (c,f) for the cubic spline radial function and in Fig.20 panel (c) for the Gaussian radial function. When p=2, C-IGA with a large s is faster than that with a small s when the error is10−4, indicating that the increased computation cost for higher s is compensated by enhanced accuracy. A similar trend is observed in Fig.19 panel(f) when p=3. In general, FEM is slightly faster than C-IGA when the error level is high (problems with small DoFs), but the gap is reduced or even reversed as the problem size grows. At the error level of 10−4in Fig.19 panel (c), C-IGA (s=5, p=2) is slightly slower than FEM (p=2). On the other hand, in Fig.20 panel (c), C-IGA (s=5, p=2) is slightly faster than FEM (p=2). The present invention discloses that most computational costs of C-IGA come from the construction of convolution patch functions, which can be further parallelized with multiple GPUs. Implementation of C-IGA programming on multiple GPUs is therefore the next step.

[0213] Although the Gaussian radial function is more accurate than the cubic spline radial function, there is an unstable moment matrix issue with the Gaussian radial function as the meshAttorney Docket No.0116936.296WO2 refines and s increases. This is the reason why C-IGA plots in Fig.20 have fewer points than those in Fig.19. 3.2.2 Irregular and triangular parametric mesh

[0214] As discussed in Subsection 1.7 and Fig.7, C-IGA can handle irregular meshes in the parametric domain. Fig.21 shows that the present invention can create an irregular mesh in the parametric domain (b) and directly map it to the physical domain (a) using both IGA and C-IGA mappings (they are equivalent). However, IGA cannot create shape functions based on the irregular parametric mesh shown in Fig.21 panel (b). Fig.21 panel (c) is the solution field computed from a 128 by 128 irregular mesh in the parametric domain, which still accurately captures the exact solution field displayed in Fig.18 panel (a). The convergence plot in Fig.21 panel (d) also confirms the irregular mesh can reliably capture the solution field. This naturally leads to the capability of r-adaptivity of C-IGA, similar to standard FEM.

[0215] C-IGA can also handle triangular elements in IGA framework. IGA formulated with triangular elements has never been reported in the literature. Fig.22 panel (a,b) illustrates how the triangular elements in the parametric domain are mapped into the physical domain. Again, the C-IGA mapping is identical to the IGA mapping but the IGA formulation cannot solve the problem with triangular elements. In Fig.22 panel (c), the energy norm error is plotted with DoFs for p=2 cases: IGA (Quadrilateral elements), C-IGA (s=3, Quadrilateral and triangular elements), and C-IGA (s=5, Quadrilateral and triangular elements). It is clearly shown that triangular elements can achieve the same level of accuracy as quadrilateral elements for the same DoFs. 3.32D linear elasticity: infinite plate with a circular hole

[0216] An infinite plate with a circular hole problem is solved. The problem is defined in Fig.23 panel (a) and the exact solution is given below:Attorney Docket No.0116936.296WO2

[0217] and R is theradius of as assumes a with the elastic modulus of 103and the Poisson`s ratio of 0.3. All values are unitless. One-quarter of the plate is simulated under the axis-symmetric boundary condition and the exact traction force is applied at x=-4 and y=4 surface.

[0218] The minimum order IGA mapping is defined by the knot vectors: Ξ = {0,0,0,0.5,1,1,1}, ^^ = {0,0,0,1,1,1} (54)

[0219] with order p=2. Two IGA elements are defined with these knot vectors: element 1 -[0,0.5] × [0,1]; element 2 -[0.5,1] × [0,1], which are denoted in Fig. 23 panel (e). These IGAelements are different from the elements on which the convolution patch functions are built. They are fixed elements once the minimum order IGA mapping is defined. It is important to note that when one construct convolution patch functions in the parametric domain, both the nodal patch AξsIand the elemental patch Aesmust be confined to one IGA element where the element e belongs to. This is because each IGA element has different basis functions. If a patch domain invades multiple IGA elements, the equivalence between the IGA and C-IGA mappings (stated in Eq. (21)) no longer holds because the reproducing condition is violated.

[0220] Since the present invention has the analytical solution, convergence studies are available. Fig.23 panel (b,c) are the results where the energy norm error is measured and the Gaussian radial function is used with a / hξ= 8. As shown in Fig.23 panel (b), the C-IGA error is slightly higher than IGA or FEM (p=2) for a given mesh size although all of them have the same theoretical convergence rate of 2. When it comes to the DoFs in Fig.23 panel (c), for a given level of error, IGA needs the smallest DoFs while C-IGA and FEM (p=2) need more than IGA. Note that we have not conducted any adaptivity for C-IGA. A careful adaptivity of C-IGA such as r,s,a,p-adaptivity might significantly reduce the error. Finally, Fig.23 panel (f) shows the σxxplot computed from a 64 by 64 regular element domain. As expected from the exact solution, there is a stress peak at r=1 and θ = ^^ / 2, which is 3.0. 3.3.1 Removal of lockingAttorney Docket No.0116936.296WO2

[0221] In one embodiment, the present invention changed the problem to plane strain elasticity and Poisson's ratio to 0.4999 to study the locking behavior. As shown in Fig.24, a severe locking is observed from linear FEM (p=1) while the other cases are locking-free. The IGA is expected to be locking-free because the minimum order of the NURBS basis of the current problem is 2. C-IGA is also locking-free because we can achieve an arbitrary reproducing order according to the C-HiDeNN interpolation theory. As long as p is greater than 1, C-IGA can remove locking similar to C-HiDeNN. Especially when p=2, C-IGA (s=3) converges faster than C-IGA (s=2). 3.43D linear elasticity: solid circular cylinder subjected to internal pressure loading

[0222] In one embodiment, a 3D linear elasticity problem is solved. The circular cylinder shown in Fig.25 panel (a) is constructed from the control points and weights provided according to one embodiment of the ppresent invention. Elastic modulus and Poisson's ratio are set as E=1000 and ^^ = 0.3, respectively. The internal pressure of P=100 is applied. All quantities are unitless. A 32 × 16 × 8 parametric mesh is generated in Fig.25 panel (b) while the corresponding physical mesh is shown in Fig.25 panel (a).

[0223] The exact solution to this problem under plane strain assumption is:

[0224] are the inner and outer radius of the cylinder, respectively. The ^^2norm error convergence study is plotted in Fig.25 panel (c). The convergence rate of C-IGA (s=2, p=2) is close to its theoretical rate of 3.0 for coarse meshes but is reduced to 2.0 for fine meshes. Increasing s or p to 3 did not improve the accuracy and convergence rate. One attributes these suboptimal convergence rates to several reasons: 1) a suboptimal convergence of error on the boundary; 2) a relatively smooth exact solution; and 3) the rational form of the ^^(ξ) basis function. The basis functions ^^(ξ) in C-IGA are set as polynomial functions divided by the NURBS weighting function to reproduce exact geometry. This might cause the suboptimal convergence for finer meshes. Further studies on these three aspects are needed. 3.5 Limitations of the current computational schemeAttorney Docket No.0116936.296WO2

[0225] Although the present invention has demonstrated the various capabilities and versatile adaptivity of C-IGA, there are some limitations of the current computational scheme that need further investigation.

[0226] When one imposes the Kronecker delta property to the convolution patch function, the moment matrix becomes unstable as the order p increases, especially when the Gaussian radial function is used. This issue occurs at large s or large a. Hence, it is recommended to take a smaller dilation parameter a for a larger patch size s, or a larger dilation parameter a for a smaller patch size s according to the order p. Considering that large s requires too much computation costs, the latter set might be a more efficient choice. Further mathematical analysis is needed to determine the optimal combination of s and a for a given p. %It is recommended to increase the patch size s according to the order p, but large s requires too much computation costs.

[0227] In the case of the Gaussian radial basis, when the dilation parameter a is suboptimal, the iterative solver does not converge. This might hinder the automatic a-adaptivity.

[0228] Most CAD models have only 3D surface information although the 3D IGA is defined on the 3D volumetric mesh. There are many approaches to resolve this issue within the IGA formulation, and some of these may be suitable with C-IGA.

[0229] As discussed in subsection 3.4, the suboptimal convergence especially for large p needs to be investigated. Example 2: Multi-Patch C-IGA Formulation without Compatibility Conditions

[0230] Multi-patch IGA is formulated by expanding the NURBS mapping to multiple non- overlapping and fully matching patches (B-splines follow the same procedure, but they are omitted here for brevity). Assume the two-dimensional physical domain Ω consists of ^^ NURBS patches whose mappings are denoted by: ^^( ^^): Ω → Ω( ^^), ^^ = 1, … , ^^(56)such thatΩ= ⋃^^^^=1 Ω( ^^).(57)

[0231] The superscript ( ^^) denotes functions or variables that differ by patches. Here, Ω(α)is a closed set. The boundary of a patch Ω(α)is denoted by Γ(α). The Γ represents the boundary of the whole computational domain Ω which is imposedboundary conditions.Attorney Docket No.0116936.296WO2

[0232] Each patch is connected to one or several neighboring patches. Patches Ω(α)and Ω(β)(α < β) are called neighboring patches if Ω(α)⋂ Ω(β)has positive ( ^^ − 1)-dimensional measure. Their interface is denoted by Γ(α,β)= Γ(α)⋂ Γ(β), ^^ < ^^. The set ℱ consists of indices of all neighboring patches, i.e., ℱ(α)and Ω(β)are ^^ ^^ ^^ ^^ℎ ^^ ^^ ^^ ^^ ^^ ^^ ^^ ^^ ^^ ^^ℎ ^^ ^^, ^^ < ^^}. Eachmapping has its own ^^1, ^^2. The discretized solution space in the physical)domain at patch(^^)is defined as ^^ℎ( ^^)= span{^^(^^^^)( ^^) ≔ ^^(^^^^)∘ ^^( ^^)−1( ^^), ^^ ∈ ℛ( ^^)}, (58)

[0233] where the patch-wise inverse mapping ^^( ^^) (^^)is uniquely defined from the Bijective Assumption.

[0234] Bijective Assumption: For each patch Ω( ^^), the absolute value of the Jacobian determinant assumed to be always greater than a certain positive constant. That is, there exists ^^ > 0 such that:|det (∇( ^^(α)))| > ^^, ^^ = 1,2, … , ^^.

[0235] In0^^of the solution field across patch interfaces (i.e., ^^0continuity – globally ^^0) is guaranteed if the patches are fully matching. To achieve fully matching geometrical mappings, the conformity assumption is required:

[0236] Conformity Assumption for IGA: Let Γ( ^^, ^^)= Γ(α)∩ Γ(β), ( ^^, ^^) ∈ ℱ be the interface between the neighboring patches Ω( ^^)and Ω( ^^).invention defines a set of indices of patch ^^ whose support domain intersects the interface Γ( ^^, ^^)as ℬ( ^^, ^^)={^^ ∈ ℛ( ^^): supp ( ^^( ^^)^^) ∩ Γ( ^^, ^^)≠ ∅}. The two patches fully match if the two following conditions are met.(1) The interface Γ( ^^, ^^)is the image of a full edge of the respective parameter domains. (2) For eachℬ( ^^, ^^), there exists a unique index ^^ ∈ ℬ( ^^, ^^), such that:^^( ^^)^^|Γ(^^, ^^)= ^^( ^^)^^(^^, ^^).

[0237] Given the fully matachieving ^^ continuity in the solution field is straightforward; the two matching basis functions ^^( ^^)and ^^( ^^)^^ ^^should have the same control variables ^^( ^^) ( ^^) ( ^^) ( ^^)^^= ^^^^and control ^^^^= ^^^^. Since the control variablesAttorney Docket No.0116936.296WO2 are the degrees of freedom (DoFs) of IGA, the imposition of ^^( ^^)^^= ^^( ^^)^^is simply conducted by generating a connectivity array like FEM.

[0238] Within a single CAD patch, C-IGA achieves ^^0at the element boundary and higher order of continuity inside each element. On the other hand, for a multi-patch system, C-IGA interpolation may suffer from compatibility issues at patch interfaces, similar to IGA without the Conformity Assumption. In this section, the present invention discusses the effect of this non-conforming patch interface on the convergence rate of a general interpolation.

[0239] First, the present invention defines a general interpolation operator ℐ. The interpolation of a continuous function in a Sobolev space ^^^^+1(Ω), ^^(^^)∈ ^^^^+1(Ω)( ^^ ≥ 0) reads: ^^ℎ( ^^) = ℐ ^^( ^^), (59)

[0240] where the interpolated function ^^ℎ(^^)is constructed with multiple non-overlapping patches, each of which is parameterized differently. Thus, ^^ℎ( ^^) is patch-wise continuous (i.e., ^^ℎ,( ^^)∈ ^^1(Ω( ^^))) and can be discontinuous at patch interfaces due to incompatiblepresent invention defines the following broken Sobolev space: ^^^^^^(Ω)= { ^^ ∈ ^^2(Ω): ^^( ^^)∈ ^^^^(Ω( ^^)), ^^ = 1,2, … , ^^}. (60)

[0241] where ^^ℎ ^^+1^^ ∈ ^^^^(Ω) and ^^^^(Ω) is a subset of this setup, the present invention investigates the convergence of the interpolation of a continuous function and the solution of a partial differential equation (PDE) in. 1.1 Interpolation Estimate on a Multi-Patch System without Compatibility Conditions

[0242] The present invention defines a norm in ^^^^^^(Ω) as: ^^^^^^= ^^(1)^^ 1 +⋯+ ^^( ^^) 2^^ ^^ .

[0243] If thepatch Ω( ^^), the interpolation of a continuous function ^^ throughout the whole domain Ω can also reach the ^^-th order convergence rate. See Lemma 1 below.

[0244] Lemma 1:

[0245] For ^^( ^^) ∈ ^^^^+1(Ω), assume that the ^^-th order interpolation estimate for ^^ℎ= ℐ ^^ holds in each single patch ^^( ^^), ^^ = 1,2, … , ^^, i.e.,Attorney Docket No.0116936.296WO2 ‖^^( ^^)−^^ℎ,( ^^)‖^^1(Ω(^^() ^^)) ≤ ^^^^ ^^ ^^ ^^ ^^ ^^ℎ^^‖^^‖^^ ^^+1(Ω(^^)), ^^ = 1,2, … , ^^, (62)

[0246] size ℎ. Then thefollowing interpolation estimate in the whole domain ^^ holds: ‖ ^^ − ^^ℎ‖^^ ^1^(Ω)≤ ^^^^ ^^ ^^ ^^ ^^ ^^ℎ ^^‖ ^^‖^^^^+1(Ω) (63)

[0249] In contrast to the interpolation estimate of a continuous function discussed in the previous section, the convergence rate of a solution to a PDE is suboptimal due to the compatibility issue at patch interfaces. Here, the present invention takes an elliptic PDE to illustrate the effect of compatibility.

[0250] Consider the following elliptic PDE: −Δ ^^ + ^^ = ^^ (64)

[0251] with ^^ being a source term. The boundary ^^ of the entire domain Ω consists of Dirichlet boundary ^^Dassigned with ^^|^^^^ = ^^^^, and Neumann boundary ^^Nassigned with ^^ ^^ ^^ ^^|^^= ^^. Note that ^^D≠ ∅ and ^^D, ^^Ndo not overlap. ^^The solution ^^ ∈ ^^^^+1(Ω) of Eq. (64) satisfies the following strong form PDE for any test function ^^ ∈ ^^^1^(Ω),

[0253] By^^ ^^, − ∑ ( ^^(^^, ^^),^^ ^^^(^^, ^^))= ^^, + ^^,^^^^,

[0254] The^^ ^Attorney Docket No.0116936.296WO2 1

[0255] Note that ( ^^( ^^, ^^)^^( ^^))2is precisely the ^^^1^ norm of ^^ over the whole domain Ω, that is,‖^^‖^^ ^1^(Ω). Here, the introduces the following^^ ^^^1^(Ω):〈^^〉 ( ^^, ^^) = ( ^^( ^^)| + ^^( ^^)| ) – mean v ( ^^, ^^)^^ 2 ^^( ^^, ^^) ^^( ^^, ^^) alue of the traces of ^^ at the interface ^^ .–jump of ^^ at the interface ^^( ^^, ^^).may find thenorm‖⋅‖in ^^^(Ω)like the one in the single patch, i.e., 1 ^^1‖ ^^‖ ≤ ( ^^( ^^, ^^))2≤ ^^2‖ ^^‖,(68)

[0257] where ^^ and ^^2arebe measured by this equivalent norm‖⋅‖.

[0258] The weak form in a multi-patch domain Ω reads: ^^( ^^ℎ, ^^ℎ) + ^^ ∑ ([^^ℎ]^^( ^^, ^^),[^^ℎ]^^( ^^, ^^))ℎ^^( ^^, ^^)= ( ^^ , ^^) + ( ^^ℎ, ^^)^^ ^^, ^^

[0259] conditions on ^^D, i.e., ^^ℎ= { ^^ℎ| ^^ℎ= ℐ ^^, ^^ ∈ ^^^^+1(Ω) and ^^|^^^^ = ^^^^}, (70)

[0260] and ^^ℎisD^^ , i.e., ^^ℎ={^^ℎ| ^^ℎ= ℐ ^^, ^^ ∈ ^^^^+1(Ω)and ^^|^^^^= 0}. (71)

[0261] The term ^^the uniqueness of the solution, where ^^ is apenalty coefficient.

[0262] The weak form of a single patch system will lead to ^^-th order convergence rate if the ^^-th order interpolation estimate is achieved. However, when solving a multi-patch system, the compatibility issue will limit the convergence according to Theorem 1. The ^^2norm at the interface Γ( ^^, ^^)is defined by ‖ ^^‖^^2(Γ( ^^, ^^))=∫Γ( ^^, ^^)^^2^^ ^^( ^^).

[0264] Let ^^ℎ∈ ^^ℎbe the trial function of the weak form:Attorney Docket No.0116936.296WO2 ^^( ^^ℎ, ^^ℎ) + ^^ ∑ ([^^ℎ]^^(^^, ^^),[^^ℎ]^^(^^, ^^))^^( ^^, ^^)= ( ^^ℎ, ^^) + ( ^^ℎ, ^^)^^ ^^, ^^( ^^, ^^)∈ℱ

[0265]

[0266] an estimate of the discontinuity at the interface Γ( ^^, ^^)satisfies: ‖[^^ℎ]^^(^^, ^^)‖ ^^2(Γ( ^^, ^^))

[0267] where ℎ^^( ^^, ^^)is the element size around the interface Γ( ^^, ^^), ^^ ∈ ℕ and ^^^^ ^^ ^^> 0 are constants, then the following error estimate holds:‖^^ℎ− ^^‖^^ ^1^(Ω))

[0268]

[0269] Corollary: If the interpolations always achieve compatibility at the interfaces, i.e.,the trial function ^^ℎ ∈ ^^ℎ and the test function ^^ℎ ∈ ^^ℎ are both continuous at the interface:[ ^^ℎ] ^^( ^^, ^^) = 0, [ ^^ℎ] ^^( ^^, ^^) = 0, then the coefficient ^^ ^^ ^^ ^^ in Eq. (73) becomes zero. Accordingly,hand side of inequality (74) vanishes, and the following error estimate holds: ‖ ^^ℎ− ^^‖ ^^1^^(Ω)≤ ^^2ℎ ^^‖ ^^ ^^ ^^ ^^‖^^^^+1(Ω),(75)

[0270] that leads to the

[0271] Theorem 1 indicates that the convergence rate is limited by the gap of interpolations at the interface. In general, ^^ is smaller than ^^. For C-IGA interpolations, the numerical tests provided in Section 3 show that ^^ = 1. Hence, the final convergence rate without compatibility is merely first order even though the reproducing order of interpolation for each patch is greater than one. Thus, one must satisfy the compatibility conditions for multi-patch systems to ensure ^^0continuity globally (i.e., ^^0continuity), and achieve the theoretical convergence rate. 2. Multi-Patch C-IGA Formulation with Compatibility ConditionsAttorney Docket No.0116936.296WO2

[0272] The single-patch C-IGA interpolation at patch Ω( ^^)is defined as ^^ℎ,( ^^)= ∑ ^^^(^^^)( ^^( ^^)) ∑ ^^^^( ^^)( ^^^, ^^^^,^^, ^^( ^^(= ∑ ^^^^^)( ^^( ^^)) ^^ ^^^^( ^^) ^ ^^ (76)^^^∈ ^^ ^^

[0273] ^^(^^( ^^))= ∑ ^^^^( ^^)(^^( ^^))∑ ^^ ^^ ( ^^)^^, ^^^^,^^, ^^ (^^( ^^))^^ ^^ = ∑ ^^^^( ^^( ^^)) ^^ ^^.^^( ^^) ^^∈ ^^ ^^^^ (77)

[0274] C-IGA (56) are equivalent. Moreover, the interpolation in Eq. (76) has two properties: ^^-th order reproducing property, and Kronecker delta property.

[0275] The former guarantees the consistency and the ^^-th order interpolation estimate in each patch. The latter imposes interpolating features to C-IGA shape functions, which is a crucial requirement for the compatibility conditions.

[0276] Theorem 1 indicates that the convergence rate of an elliptic PDE solution in a multi-patch system is suboptimal due to the incompatibility at the interface. To resolve this, it is required to have a ^^0compatibility defined as follows:

[0277] ^^^^compatibility : The solution field achieves global ^^0continuity over the entire multi- patch physical domain Ω, having ^^0continuity across patch interfaces and up to ^^^^( ^^ = min( ^^, ^^)) continuity within each patch, where FEM shape functions are ^^^^continuous and convolution patch functions are ^^^^continuous.

[0278] This is realized by two steps:

[0279] nodal compatibility condition to guarantee the uniqueness of solution. It replaces the penalty term in Eq. (69).

[0280] ^^^^compatibility conditions for C-IGA shape functions to achieve point-wise continuous interpolations ( ^^0continuity) along the interface. 2.1 Nodal Compatibility Condition

[0281] The nodal compatibility gives constraints on the mesh construction. Due to the Kronecker delta property of C-IGA interpolations, if the two neighboring patches share the same physical nodes, the nodal compatibility is achieved. That is, ^^( ^^)(^^^^)= ^^( ^^)(^^^^), ∀ ^^^^∈Attorney Docket No.0116936.296WO2 Γ(α,β), Γ(α,β)∈ ℱ, where the subscript ^^ being the nodal index along the interface. Thus, the of the solution is guaranteed, and the penalty term can be eliminated in Eq. (69).nodal compatibility is straightforward. First, construct a conforming physical mesh in the whole domain Ω (no hanging nodes at the interface) using a standard FEM meshing tool. Then find corresponding parametric mesh for each patch using the inverse mapping. The present invention devised an inverse mapping algorithm based on the neural network, which will be described below.

[0283] The inverse mapping of a multi-patch system is defined as: ^^( ^^)−1: ^^( ^^)→ ^^( ^^), ^^ ∈ Ω( ^^), ^^ = 1,2, … , ^^.(78)

[0284] According tomapping is non-zero in the given computational domain, i.e.,^^ ^^^^ ^^( ^^)≠ 0, the inverse mapping ^^( ^^)−1exists. This property is mentioned in the previous Bijective Assumption.

[0285] In the mesh generation step of C-IGA, one needs to find out the parametric nodal coordinates ^^(^^^^)that correspond to the given physical coordinates ^^(^^^^)(here, the nodal index ^^ refers to the nodes in the entire patch Ω( ^^), not only in the patch interface). One possible approach is to use a feedforward neural network (FFNN) to train the inverse IGA mapping. This is a regression problem to be solved by minimizing the loss function below: ^^^^^^ ^^ ^^ ^^1^^( ^^)−12 ^^( ^^) − ^^( ^^)

[0286] where themapping ^^(^^^^)= ^^( ^^)( ^^(^^^^)), and ^^(^^^^)−1^^is the FFNN inverse mapping. Once the neural network is trainedinvention uses the trained inverse mapping to generate parametric mesh ^^( ^^)according t( ^^)^^o the given physical mesh ^^^^. Detailed procedures of training a neural networkthe inverse IGA mapping are given in Algorithm 1 below.

[0287] Algorithm 1: Training procedure of IGA inverse mapping.

[0288] Given: physical mesh, ^^(^^^^), and IGA forward mapping, ^^( ^^)= ^^( ^^)( ^^( ^^)).

[0289] Goal: find parametric mesh ^^(^^^^)that satisfies the IGA forward mapping.Attorney Docket No.0116936.296WO2

[0290] Generate training data. Generate ( ^^(^^^^), ^^(^^^^)) pairs using the given IGA forward mapping. The training data can be mesh.

[0291] Train a neural network for the inverse mapping, ^^( ^^)= ^^(^^^^)−1^^( ^^( ^^)).

[0292] A standard FFNN with sigmoid activation error (MSE)loss function is a fair starting point.

[0293] For the given physical mesh, compute parametric mesh from the trained inverse mapping: ^^(^^^,^^) (^^^= ^^^^^^)−1 (^^( ^^^^^^)).

[0294]

[0295] If ‖ ^^( ^^)( ^^(^^^,^^) (^^^) − ^^^^^^)‖ < ^^ ( ^^ ^^ ^^ ^^ ^^ ^^ ^^ ^^ ^^), use ^^(^^,^^^) (^^^as ^^^^^^).

[0296] with ^^( ^^)^^, ^^ ^^being the initial condition: ^^(^^^^)= ^^ ^^ ^^m( ) ( ) ( ^^)^^( ^i^n)‖ ^^^^( ^^^^) − ^^^^‖. If a unique inverse mapping exists, the loss of themachine error).

[0297] Although C-IGA formulation with nodal compatibility is straightforward to implement, it does not guarantee ^^0continuity because the C-IGA shape functions of patch(^^)are constructed only with nodes inside that patch. Thus, additional constraints should be imposed to achieve ^^0continuity. 2.2 ^^^^Compatibility Condition

[0298] The ^^0compatibility requires further constraints on the C-IGA shape functions in addition to the nodal compatibility. The objective is to achieve point-wise continuity across the interfaces of every pair of neighboring patches Ω(α)and Ω(β), that is, ^^ℎ,( ^^)(^^)= ^^ℎ,( ^^)(^^), ∀ ^^ ∈ Γ(α,β).(80)

[0299] In Theorem 2 below,on the FEM shape functions, and the convolution patch functions whose support domains intersect with the interface are sufficient to establish the ^^0compatibility.

[0300] For ease of presentation, the present invention makes some notations. Assume that two neighboring patches Ω( ^^), Ω( ^^)satisfy the nodal compatibility, that is, share the same physical nodes along the interface Γ(α,β). Consider C-IGA interpolations in Ω( ^^)and Ω( ^^)use the same parameters ^^, ^^, ^^. For anode ^^^^∈ Γ(α,β), the present invention denotes the corresponding parametric nodes in Ω( ^^)and Ω( ^^)as ^^(^^^^), ^^(^^^^), respectively. Let the nodal convolution patches for node ^^Attorney Docket No.0116936.296WO2 ^^(be ^^^^^^)^^, ^^^^( ^^) ( ^^) ( ^^)^^^^, respectively. The index set of shared nodes in the convolution patches ^^^^^^^^, ^^^^^^^^is ^^( ^^) ^^( ^^) ^^( ^^) ^^( ^^) ^ ( ^^)denoted by ^^^^ ^^= ^^^^∩ ^^^^. The index set of internal no^^ ^^ ^ ^^ ^^^^ ^^ ^^des in ^^^^is ^^^^= ^^^^\ ^^^^, and that in ^^( ^^) ^^(^^^^i^^^^s ^^^^

[0301] Theorem 2: Let two neighboring patches Ω( ^^), Ω( ^^)satisfy the nodal compatibility condition, and the C-IGA interpolations on them use the same ^^, ^^, ^^. For each shared node ^^ at theinterface ^^( ^^, ^^), if the following compatibility conditions for C-IGA interpolations hold:

[0302] (C1) ^^^(^^^)( ^^)|^^( ^^, ^^) = ^^(^^^^)( ^^)|^^( ^^, ^^) , ∀ ^^ ∈ ^^( ^^, ^^)(FEM shape functions match at the interface)

[0303] (C2) For ^^ ∈ ^^^^^^ ^^, ^^^^^^, ^^^^,^^, ^^( ^^)|^^(α, ^^) = ^^^^^^, ^^^^,^^, ^^( ^^)|^^( ^^, ^^) , ∀ ^^ ∈ ^^( ^^, ^^)^^

[0305] (C3) For ^^1∈ ^^^^^^, ^^2∈ ^^^^^^^^, ^^^^^^, ^^^^,^^, ^^1( ^^)|^^(^^, ^^)= 0, ^^^^^^, ^^^^,^^, ^^2( ^^)|^^(^^, ^^)= 0, ∀ ^^ ∈ ^^( ^^, ^^)(convolution patchinterpolation reaches ^^0compatibility. That is, ^^ℎis continuous across the interface, ^^ℎ,( ^^)( ^^) = ^^ℎ,( ^^)( ^^), ∀ ^^ ∈ ^^( ^^, ^^). 2.3 Implementation of ^^^^Compatibility Condition via Product Rule

[0306] To satisfy the three conditions in Theorem 2, the present invention constructs FEM shape functions and convolution patch functions via product rule over a regular parametric mesh around the interface. Here, the present invention considers B-spline and NURBS cases. To simplify the notations, two 2D neighboring patches, patch 1 (Ω(1)) and patch 2 (Ω(2)), are considered for illustration, and ^^(1,2)is their interface. Without loss of generality, the present invention assumes that the interface of the two patches lies along the ^^-direction in the parametric domain Ω (where ^^( ^^)= ( ^^( ^^), ^^( ^^) denotes the vector of parametric coordinates in patch ^^).

[0307] For each shared node ^^ at the interface ^^(1,2), the FEM shape functions ^^^^( ^^), ^^ = 1,2, and (^^)convolution patch functions ^^^^ ^^^^are as

[0308] the index in ^^-direction and ^^-direction, respectively. Correspondingly, the nodal index ^^ is representedAttorney Docket No.0116936.296WO2 by ( ^^, ^^). The ^^^^and ^^^^are 1D FEM shape functions in ^^-direction and ^^-direction, respectively. The (^ ) ( ^^)^^^^ ^^^, ^^^^, ^^, ^^and ^^^^^^, ^^^^,^^, ^^are 1D convolution patch functions, which satisfy Kronecker delta and reproducing (^^), ^^( ^^)) is a factor coming from the weighting function ^^( ^^)( ^^( ^^), ^^( ^^)) in NURBS,1, B-^, ^^ ) = { ( ^ ) ( ) ( )spline^^( ^ ( ^^) ( ^^) ^^ ^ ( ^^ ^^ , ^^ ^^ )

[0309] Namely, the interface. For B- as aconstant 1, and thus ^^( ^^( ^^), ^^( ^^)) = 1. At node ^^, ^^( ^^( ^^), ^^( ^^)) by ^^ ( ^^( ^^), ^^( ^^)) to ensure^^ ^^^^)Kronecker delta, ^^, ^^. Based on^^(1,2), the^^^^( ^^)( ^^) reduces to ^^^^( ^^( ^^)the convolution patch function ^^^^( ^^)^^, ^^^^,^^, ^^reduces to ^^^^( ^^)^^, ^^^^, ^^, ^^( ^^( ^^)). Therefore, if ^^^^( ^^^ ( ^^)( ^^^, ^^^^, ^^, ^^( ^^( ^^)), ^^ = 1,2 are the same at the interface, the three compatibility conditions in 2 are satisfied automatically. See Theorem 3.

[0310] Theorem 3: C-IGA basis constructed via product rule. Let two neighboring patches Ω(1), Ω(2)satisfy the nodal compatibility condition with a patch interface ^^(1,2)that lies along the ^^-direction in the parametric domain Ω . Let C-IGA interpolations in two patches beover regular parametric meshes in Ω with the same parameters ^^, ^^, ^^.

[0311] Then for constructions of FEM shape functions ^^^(^^^), ^^ = 1,2, and convolution patch functions ^^^^( ^^)^^, ^^^^,^^, ^^given in Eq. (81) and Eq. (82), if ^^^^( ^^( ^^)) and ^^^^( ^^)^^^^( ^^( ^^)), ^^ = 1,2 match at the interface, compatibilityis attained at the interface ^^(1,2).

[0312] Next, the present invention provides detailed constructions to ensure ^^^^( ^^( ^^)) and ^^^^( ^^)^^^^^^( ^^)^^ = 1,2 match at the interface. Node ^^ = ( ^^, ^^) is at the interface.share the same ^^^^( ^^(1)), i.e., ^^^^( ^^(1)) in place of ^^^^( ^^(2)) (or ^^^^( ^^(2)) in place of ^^^^( ^^(1)) ).tches share the ^^^ ( ^^)^ Two pa^ ^^^^( ^^)) given by,Attorney Docket No.0116936.296WO2 ^^^^(1) ^ (2)^^^,^^^, ^^( ^^(1)) = ^^^^^^,^^^, ^^, ^^( ^^(2)) = ^^^^^^( ^^) ^^^^+ ^^^^( ^^(1), ^^(2)) ^^^^with new basis ^^, , ,

[0313] The ^^^^^^^,^^^, ^^represents a vector of all convolution patch functions in the nodal convolution ^^( ^^) patch ^^^^^^. The ^^^^( ^^) is a vector of radial basis functions defined in the physical domain. The ^^( ^^(1), ^^(2)) isof basis functions to be reproduced, which are ^^-th order polynomial basis over a ^^(1,2). Coefficients ^^^^and ^^^^are determined to ensure Kronecker delta and reproducingspline basis( ) (^^^^( ^^^^)) becomes a vector of polynomial basis.

[0314] In ^^-direction, ^^^^( ^^( ^^))( ^^( ^^)) are constructed in a similar way as 1D cases. ( Based on these constructions,convolution patch functions ^^^^ ^^)^^, ^^^^,^^, ^^preserve the Kronecker delta property and the ability to reproduce NURBS basis functions.

[0315] Notably, the regular meshes and special constructions in Eq. (81, 82) are only required for elements that are adjacent to the interface within ^^ layers. For elements away from the interface, C-IGA shape functions are constructed similar to the single-patch C-IGA where irregular meshes are allowed. 3. Numerical Examples 3.11D Interpolations with Nodal Compatibility

[0316] The present invention considers the following 1D interpolation problem to study the convergence of interface errors in a 2D multi-patch system. As illustrated in Fig.40, the presentinvention adopts two different mappings ^^(1), ^^(2) over the same physical domain of interface ^^(1,2): ^^ ∈[0,10] with the corresponding twointerpolations ^^(1)( ^^(1)), ^^(2)( ^^(2)). The nodalcompatibility suggests that ^^(1)( ^^(1)) and ^^(2)( ^^(2)) share the1^^2, ^^3, … ] shown as red points in Fig.40.variables are denoted by [ ^^1, ^^2, ^^3, … ]. Since different mappings for each patch lead to different constructions of shape functions, the discrepancy between the two C-IGA interpolations along the interface ‖ ^^(1)(^^)− ^^(2)( ^^)‖^^2 Γ(1,2))will not vanish.Attorney Docket No.0116936.296WO2

[0317] The two mappings are constructed with quadratic B-spline basis functions with control points ^^ = 0,3,10, i.e., ^^(1)= ^^(1)( ^^(1)) = 0 ⋅ (1 − ^^(1))2+ 3 ⋅ 2 ^^(1)(1 − ^^(1)) + 10 ⋅ ⋅1D function, and 2) when the nodal values are taken from a discontinuous 1D function:,2, ^^3, ^^4, … ] = [1,2,1,2, … ]. The former makes it as an interpolation of a smoothwhile the latter appears in the estimate of a solution of an elliptic PDE. The difference between the two interpolations is measured by the following relative ^^2-norm deviation: ‖ ^^(1)( ^^) − ^^(2)( ^^)‖^^2(Γ(1,2)) .

[0318] First, let ^^(1)( ^^), ^^(2)( ^^) be the interpolations of ^^( ^^) = sin ^^, i.e., ^^^^= sin ^^^^, ^^ =1,2,3, ⋯. The convergence of the relative ^^2-norm deviation as the mesh sizedecreases isplotted in Fig. 41 panel (a). As expected from Lemma 1, the ^^2 norm deviation reaches the ( ^^ +1)-th order convergence rate. This can be understood from the following inequality:‖ ^^(1)( ^^) − ^^(2)‖^^2(Γ(1,2))≤ ‖ ^^(1) − ^^‖^^2(Γ(1,2)(2))+ ‖ ^^ − ^^‖^^2(Γ(1,2)).(85)

[0319] C-rate for the two terms in the right-hand side. Thus, the present invention has the ( ^^ + 1)-th order ( ^^2norm) convergence rate for the interpolation of a smooth function in the multi-patch system.

[0320] Second, let the nodal values be [ ^^1, ^^2, ^^3, ^^4, … ] = [1,2,1,2, … ]. Obviously, the frequency of oscillations increases as the mesh refines. As illustrated in Fig.41 panel (b), the ^^2norm deviations exhibit an order 1 convergence rate rather than the ( ^^ + 1)-th order convergence. In the error estimate of a PDE solution, the suboptimal convergence of the discontinuity across the interface obstructs the theoretical ^^-th order (for ^^1norm) convergence rate as indicated in Theorem 1. 3.2.2D Interpolations with Nodal Compatibility

[0321] Next, the present invention expands the problem definition to 2D and provide numerical proof on Lemma 1. Consider a plate with circular hole of radius ^^ = 0.5, as shown in Fig.42 panel (a). The present invention interpolates a continuous function:Attorney Docket No.0116936.296WO2 ^^ ^^= 1 −^^2 3 4^^ ( ^^, ^^)3 ^^^^2(2cos 2 ^^ + cos 4 ^^) +2 ^^4cos 4 ^^,(86)

[0322] which is force is applied to the plate. Fig.42 panelthe interpolation in the whole physical domain, Ω. When ^^ = 2, ^^ = 2 are used for both patches, the convergence rateapproaches to 3.0, which is the theoretical rate for ^^ = 2. The same has been observed for ^^ =3, ^^ = 3, which numerically proves Lemma 1.

[0323] The present invention has shown that multi-patch C-IGA with nodal compatibility can achieve theoretical convergence when it is used for interpolating a continuous function. However, as pointed out in Theorem 1, the theoretical convergence is not guaranteed for solving a PDE. The present invention demonstrates this through an infinite plate with a circular hole problem, schematically illustrated in Fig.42 panel (a). The exact solution of this problem is given as: ^^ ^^(^, ^^)= 1 −^^2^^2(3 4^^ ^2cos 2 ^^ + cos 4 ^^) +3 ^^2^^4cos 4 ^^

[0324] where thepanel (a), the exact traction (Eq. ( 87 )) is applied to the top and left end of the quarter domain under symmetric boundary condition and the hole radius is ^^ = 0.5. The energy norm error is measured: 1 / 2^^^^ ^^ ^^ ^^ ^^ ^^1^^−1^^ℎ− ^^^^ ^^ ^^ ^^ ^^⋅ ^^ℎ− ^^^^ ^^ ^^ ^^ ^^

[0325] assumes a linear elasticity with the Young’s modulus and Poisson’s ratio being 103and 0.3, respectively. As shown in Fig.43 panel (a), the solution convergence is sub-optimal for higher orders ^^ ≥ 2. This is attributed to the discontinuity at the interface; even though the nodal solutions match at the interface, the ^^0continuity is not achieved. This is clearly shown in Fig.43 panel (b), a plot for the L2norm deviation between the two patches at the interface (Eq. ( 84 )). The convergence rate remains constant as 1.0 for different polynomial orders ^^. 3.3. Convergence Study of G0CompatibilityAttorney Docket No.0116936.296WO2

[0326] The ^^0compatibility condition is provided in Theorem 2, and detailed implementations using product rule are provided in Theorem 3. In this section, the present invention provides numerical proofs on those compatibility conditions.

[0327] Consider a Poisson’s equation within the same geometry Ω as in Fig 42 panel (a), given by ∆ ^^( ^^, ^^) + ^^( ^^, ^^) = 0, ( ^^, ^^) ∈ Ω. (89)

[0328] where ∆ is the Laplace operator. The present invention manufactures an exact solution as a Gaussian hump centered at ( ^^, ^^) = (−0.5,1): ^^^^ ^^ ^^(^^, ^^)= ^^− ^^( ^^+0.5)2− ^^( ^^−1)2.(90)

[0329] The source term ^^( ^^, ^^) is then derived from Eq. (86) using the exact solution. The boundary conditions are also determined from the exact solution, i.e., ^^(^^, ^^)|^^Ω= ^^^^ ^^ ^^(^^, ^^)|^^Ω. (91)

[0330] Fig.44 shows the geometry of the problem domain (a) and the exact solution (b). The interpolated solution field in these two patches are denoted by ^^(1)(^^, ^^)and ^^(2)(^^, ^^), respectively. The Gaussian hump locates near the interface denoted as the red dashed curve in Fig.44 panel (a).

[0331] It is worthwhile to note that in the coarsest IGA mesh (before any mesh refinement), as shown in Fig.45 panel (a), two IGA elements exist in patch 1, while only one IGA element is in patch 2. When C-IGA shape functions are computed in patch 1, the convolution patch domains (the support domains of C-IGA shape functions) should not invade this IGA element boundary (same for patch 2) because the two patches have different parameterizations. If the convolution patch domains invade the IGA element boundary, the equivalence of geometric mappings among IGA and C-IGA will be violated.

[0332] For nodal compatibility, it is straightforward to construct C-IGA shape functions in each IGA element. However, if one wants to achieve ^^0compatibility where the shape functions of the two patches coincide along the interface, it is required to divide patch 2 into two IGA elements similar to patch 1, as illustrated in Fig.45 panel (b). Then the shape functions around the interface are constructed using the product rule. Detailed constructions are provided in Subsection 2.3, which ensures the three ^^0compatibility conditions in Theorem 2.

[0333] In Fig.46 panel (a), the present invention plots ^^(1)( ^^, ^^) and ^^(2)( ^^, ^^) along the interface to evaluate the discontinuity; the two curveslisted in Table 6, the relative deviation between ^^(1)(^^, ^^)and ^^(2)(^^, ^^)at the interface (defined in Eq.Attorney Docket No.0116936.296WO2 (81)) are less than 5e-12, which mostly originates from the machine precision. In contrast, the nodal compatibility leads to non-negligible relative deviations at the interface, varying from 1e-5 to 1e-7 depending on mesh refinement.

[0334] C-IGA shape functions with the ^^0compatibility at the interface is provided in Fig. 47. The C-IGA shape functions encompass multiple C-IGA elements (not IGA elements) within a support domain of a normalized size 2(^^ + 1)× 2( ^^ + 1). They satisfy the Kronecker delta property, that is, for a given node, the corresponding shape function value becomes 1 at the node and 0 at the other nodes in the support domain of the shape function. As shown in Fig.47, C- IGA shape functions are allowed to have negative values within their support, which is different from IGA basis functions (i.e., splines).

[0335] The present invention defines the relative energy-norm error to measure the accuracy: ℎ^^^ ^^ ^^ ^^ ^^ =‖ ^^ − ^^ ^ ^^ ^^‖ ^^^^^^ ^^ ^^, ‖ ^^‖^^= ∫ |∇ ^^|2d ^^d ^^.(92)

[0336] The indifferent convergence rates, as presented in Fig.46 panel (b). The multi-patch C-IGA with ^^0compatibility can achieve theoretical convergence, i.e., ^^-th order energy-norm convergence rate. For ^^ = 3, ^^ = 3, C-IGA achieves a convergence rate of 3.61. However, the nodal compatibility results in a suboptimal convergence rate of 1.19 for ^^ = 2, ^^ = 2. As pointed out in Theorem 1, the solution accuracy is significantly influenced by the discontinuity across the interface. Table 6 Relative deviations ( ^^2norm defined by Eq. (81)) between ^^(1)(^^)and ^^(2)(^^)at the interface. Nodal compatibility ^^0compatibility Mesh ^^ = 2, ^^ = 2 ^^ = 3, ^^ = 3 ^^ = 2, ^^ = 2 ^^ = 3, ^^ = 340 × 40 1.76E-05 2.42E-05 1.18E-12 3.82E-1280 × 80 6.34E-06 9.91E-06 5.71E-13 4.03E-12160 × 160 2.30E-06 3.32E-06 1.27E-12 4.64E-12320 × 320 8.02E-07 1.16E-06 1.73E-12 1.98E-12640 × 640 - - 1.28E-12 4.63E-12Attorney Docket No.0116936.296WO2

[0337] Multi-CAD-patch (in short, multi-patch) C-IGA theory and its compatibility conditions are developed and validated with numerical examples. The multi-patch C-IGA formulation accompanies incompatibility of the solution field at patch interfaces, resulting in a suboptimal convergence of the solution accuracy when solving a partial differential equation (PDE). The two compatibility conditions are proposed and investigated: 1) nodal compatibility and 2) ^^0compatibility conditions. The former imposes the condition during the meshing stage, that is, the physical nodes of adjacent patches are set to be matched at the interface. This can be achieved by leveraging standard FEM meshing tools and inverse geometrical mapping of IGA. The nodal compatibility only guarantees ^^0continuity at those matching nodes. To achieve point-wise continuity across the interface, ^^0compatibility conditions are required. This is realized by constructing C-IGA shape functions with the product rule such that the shape functions at the interface are only affected by the nodes at the interface.

[0338] In terms of numerical implementation, the nodal compatibility is more practical because it only requires inverse mappings for each patch. The ^^0compatibility necessitates reparameterization of the near-interface region into a regular mesh. For most engineering applications where lots (even thousands) of CAD patches are involved for the geometrical representation, this reparameterization step will substantially impede the overall numerical analysis and complicate programming. In the future, the implementation of nodal compatibility to a realistic engineering problem with hundreds or thousands of CAD patches will be conducted.

[0339] Another critical and promising research direction is to leverage the multi-patch theory as a starting point of tensor decomposition (TD) – based reduced order modeling. In traditional TD or proper generalized decomposition (PGD) reduced order methods, a regular cartesian grid was indispensable, limiting their application to complex 3D geometries. The multi-patch IGA theories proposes a systematic reparameterization of a 3D complex CAD geometry into a regular parametric mesh. PGD has already been implemented to a single-patch domain. It is also recently reported that the C-HiDeNN-TD provides a powerful solver for space-time-parametric (S-T-P) PDEs where the spatial domain is represented with regular grid. If is successfully merged with the multi-patch theory that provides a mapping from a regular mesh to a complex 3D geometry, it will be able to handle S-T-P PDEs with complex geometries, which has not been accomplished by other numerical methods.Attorney Docket No.0116936.296WO2 Applications and Conclusions Multi-GPU programming

[0340] The recent surge in the development of generative AI, particularly due to the ChatGPT hype, has resulted in unprecedented attention to multi-GPU-based parallel computing (or accelerated computing). NVIDIA, which provided more than 10,000 A100 GPU computers for ChatGPT, has joined the $1 trillion valuation club in June 2023. Many other big tech companies are also competing to develop parallel supercomputers for AI.

[0341] It may seem like these amazing computers are still far off for ordinary researchers, but computers evolve very quickly, and models often become outdated within a few years of development. It is asserted that in 5 to 10 years, the most recent NVIDIA supercomputer will be available for use in at least university research labs.

[0342] Therefore, it is necessary to develop mechanistic machine learning tools suitable for rapidly evolving AI supercomputers.

[0343] One of the key features of the GPU solver is that it can solve a linear system of equations considerably fast if one can store all the variables in GPU memory. For example, solving a 1 million DoFs linear FEM takes around 10 seconds with GPU, while it takes more than 5 minutes with Abaqus running on CPU. However, as the size of the GPU memory is limited, one is limited by the problem size. The largest problem that one can solve with a single A6000 GPU (48 GB memory) is around 10 million DoFs. However, NVIDIA's DGX GH200 supercomputer, scheduled for release later this year, supports 144 TB of GPU memory, which means it can solve problems with up to 30 billion DoFs using highly smooth interpolation. It is not easy for existing FEM tools to solve 1 billion DoFs with linear elements. It is expected that C-HiDeNN and C-IGA to become the most widely used FEM analysis methods in the parallel computing era. Isogeometric C-HiDeNN-TD (C-IGA-TD)

[0344] In order to enhance the efficiency of solving PDEs, especially for extra-large scale problems with large numbers of DoFs, C-HiDeNN-TD has been developed. It is a reduced order model of C-HiDeNN using the so-called separation of variables % which enables a reduced order model of C-HiDeNN using the so-called separation of variables. The C-HiDeNN-TD approximation for a 3D problem is written asAttorney Docket No.0116936.296WO2 uC−HiDeNN−TD=∑Qq=1 (∑nξi=1 Nξ,ı(ξ)d(q) nη (q) nγ (q)ξ,i ) (∑j=1 Nη,^(η)dη,j ) (∑k=1 Nγ,k(γ)dγ,k ) (91)^^, ^^, ^^- directions respectively, d(q)ξ,i, d(q)η,j, d(q)γ,kare corresponding nodal coefficients for the q-th mode, Q is the number of modes, and nξ, nη, nγare the number of nodes in the corresponding directions.

[0346] By tensor decomposition, one only needs to solve 1D problems in each direction instead of solving the full 3D problem. Thus the number of DoFs of the final discretized system is reduced from nξ× nη× nγto (nξ+ nη+ nγ) × ^^ (i.e., exponential growth to linear growth), which greatly enhances the efficiency. In general, one can expect the number of modes Q is much fewer than the number of nodes in each direction, i.e., ^^ ≪ nξ, nη, nγ.

[0347] It can be seen that the above TD form has one limitation: when the computational domain is an arbitrary domain and is not intrinsically separable, the separated representation cannot be applied directly. Some techniques have been developed to resolve this issue. For example, immersing the non-separable domains onto a fully separable one and using geometrical mapping to deal with layered domains where interfaces are not planar. In one embodiment, one resolved this issue by using FEM mapping and an affine to map an arbitrary domain into a regular domain as illustrated in Fig.26 panel (a). In the present invention, C-IGA provides astraightforward way to map an arbitrary physical domain to a regular parametric domain (Ω =[0,1]dwith dimension d). The proposed method is able to adopt the IGA mapping as thegeometric mapping directly instead of the 2-step mapping, as shown in Fig.26 panel (b).

[0348] Compared with the original strategy, geometric mapping is more mature and convenient. The parametric domain facilitates the application of TD to C-HiDeNN, leading to isogeometric C-HiDeNN-TD. Nonlinear isogeometric C-HiDeNN

[0349] C-HiDeNN theory has been successfully applied to model nonlinear materials with large deformation using the total Lagrangian formulation. As discussed before, classical numerical solvers belong to the realm of engineering software 1.0, where the solvers have to be adapted for different material models. For example, in implicit analysis, the material tangent stiffness matrix has to be derived from different constitutive equations. Based on the idea ofAttorney Docket No.0116936.296WO2 \textit{Engineering Software 2.0}, the nonlinear C-HiDeNN solver is differentiable where it uses automatic differentiation to obtain the tangent stiffness matrix for a wide range of nonlinearmaterials, as shown in the equation below.KiIjJ = ^^ ^^(fiI, ujJ) (92)

[0350] where KiIjJis thenodal force, ujJis the nodal displacement. AD can be efficiently implemented using JAX.

[0351] AD can also make the C-HiDeNN nonlinear solver very powerful for inverse problems with multiple design variables since the gradient information can be obtained directly from the forward analysis.

[0352] To demonstrate the powerfulness of nonlinear C-HiDeNN differentiable solver, Fig. 27 shows a numerical example where a 3D hyperelastic plate with a hole is under tension in the x direction. Due to the symmetry of the problem, only a quarter plate is used for the currentanalysis where the left and bottom sides are fixed and a displacement boundary condition ux =0.1 is applied on the right. The Neo-Hookean model is used as the constitutive law for theanalysis where the strain energy density is shown as: ^^ =μ(I − 3)K( )221̅+2J − 1 (93) 2

[0353] where ^^ is theI1̅= J−3I1, ^^ = ^^ ^^ ^^(F) and I1is the first invariant of right Cauchy-Green tensor. Μ = 232kPa, K=500 kPa are adopted for the current analysis. The dimensions of the quarter plate is given in Fig.27 panel (a). Fig.27 panel (b) shows the obtained maximum von Mises stress as a function of the numbers of DoFs. When the maximum von Mises stress is close to 218 kPa, nonlinear FEM requires 2435830 DoFs. However, C-HiDeNN p=1,s=1 requires 15183 DoFs (6.2% total DoFs of nonlinear FEM), C-HiDeNN p=2, s=2 requires 5850 DoFs (0.24% total DoFs of nonlinear FEM). As can be observed, nonlinear C-HiDeNN converges much faster than classical nonlinear finite element analysis.

[0354] Based on automatic differentiation, one can extend the current C-IGA framework to nonlinear materials with large deformation.Attorney Docket No.0116936.296WO2

[0355] Fig.36 illustrates a general process of flow chart of steps of C-IGA algorithm implemented on a GPU. In particular, step (1) provides a mesh to obtain a Hierarchical Deep- learning Neural Network – Finite Element Method (HiDeNN-FEM) interpolant based on the mesh; step (2) computes an elemental patch information and a nodal connectivity based on the mesh and graph theory using ( ^^^^^^, ^^^^^^); step (3) computes an assembled moment matrix G in the parametric domain using the patch information and a calculation parameter; step (4)computes a convolutional patch function and derivative in the parametric domain ^^^^^^,^^^^, ^^, Ә ^^^^^^,^^^^, ^^ / Ә ^^; step (5) obtains a shape function derivate in a physical domain Ә / Ә ^^; step(6) integrates the HiDeNN-FEM shape function and derivative ^^^^, Ә ^^^^ / Ә builds aglobal stiffness matrix and force vector system of equations the C-HiDeNN interpolant; and step (8) solves the global stiffness matrix and force vector system of equations using a matrix solver.

[0356] The step (2) is based on graph theory, while steps (3)-(8) adopts GPU computation. Isogeometric C-HiDeNN-TD (C-IGA-TD) for topology design and manufacturing C-HiDeNN in nested topology design

[0357] C-HiDeNN theory has been applied for nested scales topology optimization in a bamboo-like structure with tensor decomposition (TD), to address a seven-scale concurrent design problem. Bamboo is a prevalent natural material with a nested structure, offering a lightweight design while retaining considerable strength. The presented results are based on mechanistic nested topology design theory, implemented through the C-HiDeNN-TD approach.

[0358] C-HiDeNN theory is an excellent choice for topology optimization problem as it provides a built-in convolution filter with arbitrary smoothness. The design goal of this example problem is to optimize a nested bamboo-like structure under applied axial, bending, and torsional force for minimum compliance. The equality constraints are the equilibrium conditions and the volume fractions of the materials. The present invention applies a Lagrange multiplier to enforce these constraints while minimizing the objective.

[0359] Three key parameters are considered for the nested TO design: nested scales ([n]), design in each scale ([d]), and material choice or design (m). The resulting optimized structure for [n]=7 is presented in Fig.28. With a constant volume fraction of 0.7 in each scale, a nested self-similar fractal design is obtained. The video of the 7-scale topology design can be found inAttorney Docket No.0116936.296WO2 the supplementary material. On the other hand, using different volume fractions in each scale yields a non-self-similar design. These results indicate that the method can generate fractal and beyond fractal type, opening new opportunities for multiscale topology optimization in product design. An animation of these structures is provided in the supplementary video. A summary of these designs is presented in Table 7 which demonstrates the capability of the C-HiDeNN-TD- TO (Topology Optimization) for an extremely large scale high-resolution (∼ 1030voxels) topology design problem.

[0360] The performance of the bamboo-like structure for the seven scales is presented in Fig. 29. Introducing new lattice scales increases the stiffness for axial, bending, and torsional loading up to the fifth lattice scale. However, after that, introducing additional scales scarcely contributes to the stiffness of the structure. The present invention printed the predicted structures using PLA filament as shown in Fig.30. The zoomed-in version of the structure highlights the method's capability to create a smooth design. Table 7 Design parameters, objectives, performance for the nested scales design of bamboo-like structure Design Parameters Objective Force Resolution Optimized performance (N.mm) Self-similar[n]=1,.Compliance Axial 2.27 × 1029Axial design[..,7^^] = 1, ∀[ ^^]Bending compliance: [m]=1 Torsion 140.9543 Bending compliance: 237.8368 Torsion compliance: 279.773 Non-self-[n]=1,...,7Compliance Axial 2.27 × 1029Axial similar[ ^^] = 1, ∀[ ^^]Bending compliance: design [m]=1 Torsion 124.5211 Bending compliance: 215.2245 Torsion compliance: 286.2291Attorney Docket No.0116936.296WO2 C-IGA-TD in nested topology design of advanced materials systems

[0361] C-IGA-TD can facilitate complex geometries and boundary conditions in engineering product design. By utilizing the isogeometric technique, the C-IGA-TD topology optimization (TO) method directly incorporates CAD files into the design process, thereby eliminating the tedious post-processing for manufacturing. Moreover, employing the same shape function for both the solution and optimization ensures an accurate representation of the CAD file within the optimization scheme. C-IGA-TD holds the potential to be a powerful tool for transitioning from traditional Engineering Software 1.0 to the more advanced Engineering Software 2.0 across various applications, including topology design and additive manufacturing of advanced materials systems. Conclusion

[0362] An isogeometric version of the C-HiDeNN framework, called C-IGA, has been developed. C-IGA introduces an intermediate domain, parametric domain, and corresponding geometric mapping like IGA. The given geometric mapping is able to transform the parametric domain into the exact geometry in the physical domain. C-IGA reproduces such geometric mapping in the isoparametric framework by 1) assigning discrete physical nodes according to the geometric mapping and 2) reproducing the basis of the mapping function. C-IGA inherits the advantages of C-HiDeNN in the aspects of high-fidelity solutions and versatile adaptivity. Compared to FEM and IGA, C-IGA can achieve arbitrarily high-order smoothness and convergence rates without increasing DoFs. In addition to the order elevation, C-IGA allows the adjustment of other parameters to improve the performance, including 1) the patch size s that controls the number of neighboring nodes contributing to the nodal shape function and 2) the dilation parameter a that controls characteristic length scales. Another benefit of C-IGA is the Kronecker delta property, which eases the enforcement of boundary conditions and avoids the loss of optimal convergence rates. Lastly, unlike IGA, C-IGA can handle irregular meshes or low-order elements such as 3-node triangles without sacrificing accuracy.

[0363] Both theoretical analysis and numerical tests show that C-IGA can achieve (p+1)-th order of convergence rates for L2norm error and p-th order of convergence rates for H1norm error, if the interpolants satisfy p-th order reproducing conditions. Importantly, the adaptive nature of C-IGA is significantly flexible. The r,a,p-adaptivities are investigated in the numericalAttorney Docket No.0116936.296WO2 examples. The r-adaptivity optimizes nodal positions for more accurate mesh by minimizing the variational energy. The a-adaptivity searches for the optimal distribution of patch-wise dilation parameter a of convolution patch functions associated with each node. Both r-adaptivity and a- adaptivity improve the solution accuracy, especially for the coarse mesh. The present invetion explores the idea of p-adaptivity that adopts high-order interpolations in the important regions and low-order ones in the background region, resulting in a more efficient solution scheme.

[0364] The future research direction lies in developing multi-GPU programming of C-IGA, and its reduced-order extension, Isogeometric C-HiDeNN-TD, aiming to achieve significantly improved performance for extremely large-scale optimization problems. Isogeometric C- HiDeNN-TD will be designed to effectively handle boundary conditions and integrate CAD files seamlessly into the design process. This will significantly enhance the smoothness and accuracy of optimized designs with dramatically reduced DoFs and orders of magnitude of speed up. The ongoing development and application of C-IGA methodology pave the way for the transition from Engineering Software 1.0 to Engineering Software 2.0, unlocking the immense potential for interconnectivity between simulation-based product development, industry-standard design, and real-time controlled additive manufacturing. HiDeNN approximation

[0365] Hierarchical Deep-learning Neural Network (HiDeNN) shown in Fig.31 is the foundation of the recently developed convolution HiDeNN (C-HiDeNN) theory. The key idea of HiDeNN is that one can rewrite finite element shape functions with partially connected neural networks and simple activation functions (e.g., ReLU and identity) whose weights and biases are functions of nodal coordinates and nodal field variables. These elemental neural networks (Fig. 31 panel (b)) hierarchically construct ``HiDeNN'' in Fig.31 panel (c), whose input is element- wise Gauss quadrature (or integration) points and loss function is the variational energy. If one fixes the nodal coordinates and optimize the nodal field variables, the solution is mathematically equivalent to FEM. If the nodal coordinates are also optimized, one can achieve r-adaptivity, which in general improves local accuracy because the nodes tend to move toward the high gradient region. 1. Radial Basis functionsAttorney Docket No.0116936.296WO2

[0366] The radial basis function is a specific option for the convolution patch function ^^ ξ^^^^, ^^that simultaneously satisfies Kronecker delta and reproducing conditions, defined in Eq. (4). Denote a vector Waξ,is = [ ^^ξ^^,^^1 , ^^ξ^^,^^2 , … , ^^^ξ^,^^^^ ^^] composed of all convolution patch functions in the nodal patch Aξsi. , ^^ = 1,2, … , ^^ ^^. ns is the number of nodes inAξsI,i.e., ^^ ^^ = ^^ (Aξsi). Here, ^^ refers to parent coordinates in parent domain for C-HiDeNN, or parametric coordinates in parametric domain for C-IGA. For a scalar patch size s in alldirections, ^^ ^^ = (2s + 1)d with dimension d. One can write Eq. (4) in the matrix formWξia,s (ξ) = Ψa(ξ) ^^ + ^^(ξ) ^^, (B.1)

[0367] where ^^ = [α1,matrices. The radial functions Ψa(ξ) = [Ψa,1(ξ), Ψa,2(ξ), … , Ψa,ns(ξ)] can be chosen as cubic spline functions Eq. (6) or , P (ξ), … , P (ξ)] ( ^^ ≤ ^^ ^^) are the functions to be1 mreproduced in the reproducing conditions, e.g., polynomials ( ^^(ξ)=[1, ξ, ξ2, … , ξp]in 1D for example, m=p+1).

[0368] To satisfy Kronecker delta and reproducing conditions, one sets coefficients to be ^^ = R−1(Ins×ns− QK), ^^ =(QTR−1Q)−1QTR−1(B.2)

[0369] ^^ and ^^ are both matrices composed of nodal values of Ψa,j(ξ)and Pk(ξ):Attorney Docket No.0116936.296WO2

[0370] Detailed

[0371] The Kronecker delta and reproducing conditions are (Kronecker delta) (B.5)are as (Kronecker delta) (B.7) (Reproducing conditions) (B.8)

[0373] , ^^ = R−1(Ins×ns− QK) (B.9)

[0374] Substituting it into the second equation Eq. (B.8), one has Ψa(ξ)R−1(Ins×ns− QK)^^ + ^^(ξ)^^ ^^ = ^^(ξ), (B.10)

[0375] and then Ψa(ξ)R−1^^(Im×m− KQ)= ^^(ξ)(Im×m− KQ). (B.11)

[0376] Finally, one has (Ψa(ξ)R−1Q − P(ξ))(Im×m− KQ)= 0. (B.12)

[0377] For, ^^ ^^ = Im×m. (B.13)Attorney Docket No.0116936.296WO2

[0378] Since ^^ is an ^^ ^^ × ^^ matrix, and might not be a square matrix, the coefficient matrix ^^ is the pseudo-inverse of ^^ which is not unique for ^^ < ^^ ^^. The radial basis function adopts one specific solution ^^ = (QTR−1Q)−1QTR−1. (B.14)

[0379] Remark that the least-square form ^^ = (QTQ)−1QT(B.15) is also the solution to Eq. (B.13). In fact, ^^ can be ^^ =(QTMQ)−1QT^^ (B.16)

[0380] with arbitrary ns-dimensional invertible matrix ^^. Remark that one can choose different ^^ to adjust the condition number of QT^^ ^^. 2. Elemental interpolation of C-HiDeNN in matrix form

[0381] The C-HiDeNN formulation becomes similar to FEM if one re-writes the interpolation in the matrix form: uh,e T(ξ ) = ∑Ti∈Ae (Nsi(ξ )) Usi= (Nse) Use(C.1)

[0382] s^^ s. in Eq. (C.1) is only over the elementthe first summation in Eq. (3). The i-th node in the ξeelement e is denoted by ξeı. UsIis a matrix of nodal displacements in the nodal patch Asıas eillustrated in Fig. 3(b). It has a dimension of [n (Aξıs ) , d] if the dimension of solution is d.Wξeı (ξ ) is a vector of convolution patch functionswith nodes ξe^ in the nodalξe e esı. Likewise, the length of the column vectors Waξ,ıs and Nsiequals ^^ (Aξs ), that is, it depends on the patch size s.

[0383] In FEM, nodal displacements are dealt as a column vector under Voigt notation. C- HiDeNN can do the same with the two column vectors, shape function vector ^^^^^^and nodal displacement vector ^^^^^^, with a length of ^^( ^^^^^^) × ^^. Note that the components of ^^^^^^are from ^^^^(ξ ) in Eq. (3).Attorney Docket No.0116936.296WO2

[0384] The variational energy (for the elastic system) in one element e is

[0385] body forc^^e. ^^^^is the strain-displacement matrix whose size is[3, ^^(^^^^^^)× ^^]^^ is the modulusmatrix while ^^ is the Jacobian matrix defined by ^^(ξ ) = ^^ ^^ / ^^ξ (C.3)

[0386] with the isoparametric mapping ^^(ξ ).

[0387] The integral in the parent domain Ω^^= [−1,1]^^can be realized by Gauss quadrature in the same way as FEM:3. Design of ^^( ^^) to reproduce B-spline and NURBS basis functions

[0389] In C-IGA, the present invention needs to design ^^(ξ) in the convolution patch functions Eq. (D.1) in order to reproduce the basis functions of geometric mapping. Especially for IGA mapping or CAD model as geometric mapping, B-spline and NURBS basis functions need to be reproduced accordingly. A straightforward way is to choose B-spline or NURBS basis functions ^^^^(ξ)as ^^(ξ)directly. An alternative way is to use polynomials or polynomials over a weighting function (denominator of NURBS basis functions) instead, which is convenient to be rescaled in the parent domain.

[0390] The first way is to use B-spline or NURBS basis functions ^^^^(ξ) as ^^(ξ) directly:Attorney Docket No.0116936.296WO2 ^^(ξ)=[R1(ξ), R2(ξ), …](D.1)

[0391] The second way is to use polynomials or polynomials over a weighting function instead. The present invention notices that B-spline basis functions are polynomials in their elements. Thus in one element, one can use polynomials instead, e.g., 1D p-th order polynomials for p -th order B-spline basis functions, ^^(ξ)=[1, ξ, ξ2, … , ξ^^](D.2)

[0392] or in 2D,^^(ξ)=[1, ξ, η, ξ2, η2, … , ξpηp](D.3)

[0393] It is important to note that the number of basis functions for C-IGA is determined by the tensor product while that of C-HiDeNN is from the complete polynomials of Pascal's triangle. For example, when p=2 and the problem dimension is 2D, the basis functions for C-HiDeNN are ^^(ξ) = [1, ξ, η, ξ2, ξη, η2] while those of C-IGA are ^^(ξ) =[1, ξ, η, ξ2,C-IGA requires more polynomials in general.

[0394] As for NURBS basis functions, they are represented by B-spline functions over a weighting function ^^(ξ). Thus one can use polynomials over the same weighting function ^^(ξ)instead, e.g., 2D rational form for p-th order NURBS basis functions is ^^(ξ, η)=[1, ξ, η, ξ2, η2, … , ξ^^η^^] / ^^(ξ, η)(D.4) 4. Derivatives with respect to dilation parameter a

[0395] This Subsection studies the derivatives of variational energy with respect to dilation parameter a for a -adaptivity or seeking for optimal global a by optimization. Although in an isoparametric framework, the first mapping is precisely given geometric mapping, which is independent of a. The second mapping Eq. (15) is also independent of a. Thus one only needs to consider the derivatives of interpolants with respect to a in the variational energy Eq. (16), i.e.,Attorney Docket No.0116936.296WO2

[0396] Note that the coefficient matrices ^^ and ^^ in the convolution patch functions ^^^^ξ,^^^^(ξ; ^^)(see Eq. (B.1) and Eq. (B.2)) are also function of a: ^^( ^^) = ^^−1( ^^)( ^^^^ ^^× ^^ ^^− ^^ ^^( ^^)), ^^ = ( ^^^^^^−1( ^^) ^^)−1^^^^^^−1( ^^). (E.2)in the matrix form is(E.5)rates and condition number of stiffness matrix

[0400] The ^^1-norm error convergence rates of C-IGA for 1D model, and comparisons with linear FEM, IGA are plotted in Fig.32. As expected, the convergence rates of C-IGA are p or even larger than p if satisfying p -th order reproducing conditions. This is consistent with theAttorney Docket No.0116936.296WO2 present invention theoretical analysis.

[0401] The condition number of stiffness matrix for 1D model are plotted in Fig.33. Among IGA (p =2,3,4), FEM (p =1,2) and C-IGA (p =2,3,4, s=p), IGA gives the smallest condition number, and quadratic FEM gives the largest one at the same mesh, as DoFs of quadratic FEM are almost twice those of other methods. The condition numbers of IGA, linear FEM and C-IGA have the same order of magnitude. The lines of different orders of C-IGA almost coincide, so one only plots one line to represent them. The present invention does the same thing for IGA. The condition number of C-IGA is close to that of linear FEM for the same mesh. In addition, it is noticed that all condition numbers are proportional to the square of element size no matter IGA (p =2,3,4), FEM (p =1,2) and C-IGA (p =2,3,4). The accuracy vs. condition number plot is shown in Fig.33 panel (b). When considering accuracy, C-IGA (p =2) gives the smallest condition number at the same level of error (e.g., 10−4energy-norm error) for fine mesh compared to FEM (p =1,2) and IGA (p =2). C-IGA with larger patch size s has a similar condition number, so only results for s=2 are presented. Remark that the condition number of C- IGA also varies with s and a altered like accuracy. The condition number shown in the Fig. might not be the smallest. 5. Degeneration of C-IGA to IGA

[0402] C-IGA can be degenerated to IGA by a special choice of nodes and convolution patch function. Note that knots {ξ^^} and basis functions { ^^^^, ^^} are not in a one-to-one correspondence in IGA, and the DoFs increase as the order elevates. Therefore, anchors { ^^^^} are used to locate basis functions in the parametric domain. For the j-th basis function ^^^^, ^^, the corresponding anchor, denoted ^^^^, is given by (G.1)

[0403] The are

[0404] If one uses these anchors as nodes to construct the C-IGA mesh in the parametric domain, C-IGA will have the same DoFs as IGA. The present invention notes that repeating anchors appear when knots are repeated and ^^ > 2, especially at the boundary of the computational domain, as illustrated in Fig.34 panel (b,c). In one embodiment, the presentAttorney Docket No.0116936.296WO2 invention modifies anchors by redistributing their positions locally in the region between the two adjacent non-repeating anchors. As shown in Fig.34 panel (c), the repeating anchors ^^1and ^^2are redistributed uniformly between ^^1and ^^3, i.e., ^^2=(^^1+ ^^3) / 2. After this modification, the present invention adopts anchors as nodes and construct C-IGA interpolation while using IGAbasis functions in the convolution function, given by ^^^^− ^^ ^^ ^^(ξ) = ∑^^^^^^(ξ; ^^^∗^) ∑^^∈ ^^^^^^^^^^^^, ^^(ξ) ^^^^. (G.2)

[0405] ^^∗ ^^ ∗^^^^, such as ^^ ^^ =[ ^^ ^^−1, ^^ ^^, ^^ ^^+1] in 1D linear case. ^^ ^^ ^^^^ called nodal patch is the set of nodes (anchors) in a patchnodal coefficient for ^^^^. The new interpolation does notsatisfy the Kronecker delta property.

[0406] In fact, it is precisely IGA interpolations, while written in the C-HiDeNN form, so the increase of patch size s does not influence the performance of interpolations. The only requirement is that patch size s should be larger than the support size of IGA basis functions.

[0407] Then the present invention uses the new C-IGA Eq. (G.2) and IGA to solve 1D problems in Subsection 3.1 of Example 1. The results plotted in Fig.35 clearly show that C-IGA results completely coincide with IGA results, which verifies the equivalence of C-IGA (Eq. G.2) and IGA. Thus one can conclude that C-IGA can degenerate to IGA by using modified anchors as nodes and IGA basis functions as convolution patch functions. 7. Geometric data of 2D model: a quarter-ring

[0408] The ring illustrated in Fig.16 has an inner radius 10, and an outer radius 20. The coarse mesh uses rational quadratic basis functions, i.e., NURBS, in two directions. The NURBS basis functions ^^^^, ^^(ξ, η)is defined by

[0409] ^^(ξ, η) = ∑3^^=1∑3^^=1^^^^ ^^^^^^, ^^=2(ξ) ^^^^, ^^=2(η) (H.2)Attorney Docket No.0116936.296WO2

[0410] where ^^^^, ^^=2(ξ), ^^^^, ^^=2(η)are 1D quadratic B-spline basis functions given in Table 8, and ^^^^ ^^are weights listed in Table 9. ^^(ξ, η) is the so-called weighting function. In this model, with above weights and B-spline basis, it is ^^(ξ, η) = 1 + (√2 − 2) ^^(1 − ξ) (H.3)

[0411] In onefunction to construct ^^(ξ) in C-IGA. In addition to NURBS basis functions, control points ^^^^, ^^in the geometric mapping Eq. (49) are listed in Table 10. Remark that the basis in ^^-direction, the radial direction, can be linear. Here, the present invention uses quadratic basis in both two directions for consistency. Table 8.1D B-spline basis functions in two directions i i=1 i=2 i=3 ^^^^, ^^= ^^(^^)( ^^ − ^^)^^^^ ^^(^^ − ^^)^^^^^^^^, ^^= ^^(^^)( ^^ − ^^)^^^^ ^^(^^ − ^^)^^^^Table 9 Weights. ^^^^ ^^i=1 i=2 i=3 j=1 1√2 / 21j=2 1√2 / 21j=3 1√2 / 21Table 10 Control points ^^^^, ^^i=1 i=2 i=3 j=1 (0, 10) (10, 10) (10, 0) j=2 (0, 15) (15, 15) (15, 0) j=3 (0, 20) (20, 20) (20, 0)

[0412] The foregoing description of illustrative embodiments of the invention has been presented for purposes of illustration and of description. It is not intended to be exhaustive or toAttorney Docket No.0116936.296WO2 limit the invention to the precise form disclosed, and modifications and variations are possible in light of the above teachings or may be acquired from practice of the invention. The embodiments were chosen and described in order to explain the principles of the invention and as practical applications of the invention to enable one skilled in the art to utilize the invention in various embodiments and with various modifications as suited to the particular use contemplated. It is intended that the scope of the invention be defined by the claims appended hereto and their equivalents.

Claims

Attorney Docket No.0116936.296WO2 CLAIMS What is claimed is:

1. A method of Convolution-Hierarchical Deep-learning Neural Network Isogeometric Analysis (C-IGA) of a geometric shape of a subject matter, comprising: performing a first mapping process between a physical domain of a subject matter and a parametric domain of the subject matter; and performing a second mapping process between the parametric domain of the subject matter and a parent domain of the subject matter; wherein the first mapping process comprises constructing a C-IGA interpolation in the parametric domain based on CAD data of the subject matter.

2. The method according to claim 1, wherein the CAD data of the subject matter comprises at least one of knot vectors, polynomial orders for each parametric direction, control points, and weights for spline basis functions.

3. The method according to claim 2, wherein the spline basis functions comprise Non-Uniform Rational B-Splines (NURBS) basis functions.

4. The method according to claim 2, wherein the step of the first mapping process comprises mapping nodal points though the geometric mapping.

5. The method according to claim 4, wherein the nodal points are discrete.

6. The method according to claim 4, wherein the C-IGA interpolation in the parametric domainis constructed using:xC−IGA(ξ) = ∑k∈Aes Nk (ξ)xGeo(ξk)7. The method according to claim 5, wherein degrees of freedom are attached to the nodal points.Attorney Docket No.0116936.296WO2 8. The method according to claim 1, wherein the C-IGA interpolation is constructed based on at least one calculation parameter.

9. The method according to claim 1, wherein the calculation parameter comprises a mesh coordinate r, a mesh size h, a patch size s, a dilation parameter a, and reproducing order p.

10. The method according to claim 9, wherein the reproducing order p and spline basis functions corresponding to the reproducing order p comprises the spline basis functions used in the CAD data.

11. The method of claim 1, wherein the CAD data comprises multiple CAD patches.

12. The method of claim 11, wherein the multiple CAD patches are linked via at least one linking method.

13. The method of claim 12, wherein the at least one linking method comprises a nodal compatibility method, wherein the nodal compatibility method comprises: building a conforming mesh in the physical domain; and conducting an inverse IGA mapping to construct a parametric mesh.

14. The method of claim 12, wherein the at least one linking method comprises a G0compatibility method, wherein the G0compatibility method comprises: building a conforming mesh in the physical domain; conducting an inverse IGA mapping to construct a parametric mesh; and reconstructing a regular mesh in vicinity of an interface of one of the multiple CAD patches in the parametric domain; wherein the C-IGA interpolant is constructed in the regular mesh.

15. A method for graphics processing units (GPU) programming of Convolution-Hierarchical Deep-learning Neural Network Isogeometric Analysis (C-IGA) of a geometric shape of a subject matter, comprising:Attorney Docket No.0116936.296WO2 (1) providing a parametric mesh to obtain a Hierarchical Deep-learning Neural Network – Finite Element Method (HiDeNN-FEM) interpolant based on the parametric mesh; (2) computing an elemental patch information and a nodal connectivity based on the parametric mesh and graph theory; (3) computing an assembled moment matrix G in a parametric domain using the elemental patch information and a calculation parameter; (4) obtaining a convolutional patch function in the parametric domain based on the steps (2) and (3); 5) integrating convolutional patch function with the HiDeNN-FEM interpolant obtained in the step (1) to generate a C-IGA interpolant; (6) building a global stiffness matrix and force vector to form a system of equations using the C-IGA interpolant; and (7) solving the system of equations from step (6) using a equation solver.

16. The method according to claim 15, wherein the convolution patch function is obtained according to nodal connectivity and a calculation parameter.

17. The method according to claim 16, wherein the calculation parameter comprises a mesh coordinate r, a mesh size h, a patch size s, a dilation parameter a, and reproducing order p.

18. The method according to claim 15, wherein the steps (1) and (3)-(8) are parallelized.

19. The method according to claim 18, wherein the steps (3)-(7) are serially conducted and stacked in a GPU memory having a memory limit; wherein the parallelized steps (3)-(7) comprises at least one variables having a variable size; and wherein when the variable size exceeding the memory limit of GPU, the steps (3)-(7) are serially conducted and stacked in a CPU memory.

20. A non-transitory tangible computer readable medium storing a program which, when executed by one or more processors, causes a computer to execute a process for GPU programming of Convolution-Hierarchical Deep-learning Neural Network Isogeometric Analysis (C-IGA) of aAttorney Docket No.0116936.296WO2 geometric shape of a subject matter, the process comprising: (1) providing a parametric mesh to obtain a Hierarchical Deep-learning Neural Network – Finite Element Method (HiDeNN-FEM) interpolant based on the parametric mesh; (2) computing an elemental patch information and a nodal connectivity based on the parametric mesh and graph theory; (3) computing an assembled moment matrix G in a parametric domain using the elemental patch information and a calculation parameter; (4) obtaining a convolutional patch function in the parametric domain based on the steps (2) and (3); 5) integrating convolutional patch function with the HiDeNN-FEM interpolant obtained in the step (1) to generate a C-IGA interpolant; (6) building a global stiffness matrix and force vector to form a system of equations using the C-IGA interpolant; and (7) solving the system of equations from step (6) using an equation solver.

21. The non-transitory tangible computer readable medium according to claim 20, wherein the convolution patch function is obtained according to nodal connectivity and a calculation parameter.

22. The non-transitory tangible computer readable medium according to claim 21, wherein the calculation parameter comprises a mesh coordinate r, a mesh size h, a patch size s, a dilation parameter a, and reproducing order p.

23. The non-transitory tangible computer readable medium according to claim 20, wherein the steps (1) and (3)-(7) are parallelized.

24. The non-transitory tangible computer readable medium according to claim 23, wherein the steps (3)-(6) are serially conducted and stacked in a GPU memory having a memory limit; wherein the parallelized steps (3)-(6) comprises at least one variables having a variable size; and wherein when the variable size exceeding the memory limit of GPU, the steps (3)-(6) are serially conducted and stacked in a CPU memory.

Images