Stiffness Matrix Formulation

Preliminary Steps for Analysis

The main difference between isogeometric analysis and finite element analysis is the type of the shape functions used for the approximation of the solution field. The substitution of Lagrange polynomials with NURBS (widely utilized in the CAD community) leads to several special characteristics of isogeometric analysis.

Shape Functions

Finite element method shape functions are usually polynomials (e.g. Lagrange polynomials). However, they hold certain disadvantages. FEM Shape functions are interpolatory at all nodes, internal and external, and have C-1 continuity at the edge. This results in the incompetence of FEM to define stresses and strains at any boundary of any element, leading in the introduction of other corrective methods such as extrapolation, in order to achieve that.
NURBS as shape functions hold an overlapping that is useful, as nearby elements are strongly connected, thus the simulation provides an improved approach to the natural problem. Continuous shape function derivatives lead to a continuous stress and strain field, minimizing the need for application of corrective methods.

Control Points

Classical FEM downsizes the natural problem of infinite unknowns to a finite number of unknowns, which are the degrees of freedom of the nodes. The position of the nodes depends on the element type. As a general rule, the nodes can be usually found in the corners and middle of the sides of the elements. They are part of the element and therefore part of the physical model. Displacements in other areas of the model can be approximated by a linear combination of displacements on the degrees of freedom. Distribution in the model is evaluated via the corresponding shape functions. In isoparametric elements, shape functions and their respective nodes are also used to approximate the geometry, thus enabling relatively complex shapes to be approximated with Lagrange polynomials.

In isogeometric analysis, NURBS are chosen as shape functions. The isoparametric concept is reversed, as geometrical mapping now defines the solution approximation. The geometrical representation is achieved through a combination of control points and their corresponding shape functions. degrees of freedom at the control points are now the unknowns.


Two ingredients of the isogeometric universe correspond to the essence of FEM “element”. The isogeometric element could be either the patch or the knot span of the patch. In order to perform exact numerical integration [4], a certain number of Gauss points are chosen for the domain of every piece of the polynomial basis functions. The domain of this piece is the knot span, which resembles the finite element of FEM. This is the reason why in this thesis knot spans and not patches are considered isogeometric elements. In IGA, shape functions are not restricted to the interior of each element (knot span). Instead, they are non-zero across  p+1 knot spans and overlap with more shape functions.
This overlapping results in a denser stiffness matrix than the classical finite element matrix with the same degrees of freedom. Apart from that, the fact that B-SPLine functions are defined in the whole domain allows for integration throughout the Patch without building local element matrices separately. Of course, this would be time-consuming and it is not advised for advanced software technologies, but serves well for research purposes, where a flexible quadrature code is needed in order to test and discover new methods and ideas.

Gauss Points- Parametric Coordinates

Figure 1. Gauss Points
(a) on Reference Knot Span and (b) transferred to Parameter Space.
​Full tensor product properties apply here as well, leading in similar equations for the other two parametric directions:

Gauss Point Number

For 2D and 3D problems, the maximum order of the deformation matrix is determined by partial derivation of piecewise polynomial shape functions. Therefore the maximum degree is p for derivation in the remaining directions.
Conclusively per knot span:
  • For 1D problems,p  Gauss points per knot span are required.
  • For 2D and 3D problems,p+1  Gauss points per knot span are required.


As mentioned before, patches are used when a change in geometry type or material occurs. Patches can also be used in any case C-1 or C0 continuity is required. If separate knot vectors are used for every patch, stiffness matrices will be produced for every patch. The separate matrices are combined into one via connectivity arrays. If the connection is watertight, the procedure is the same as the one used in classical FEM to combine local element matrices to a single, global stiffness matrix.

Elasticity Matrix

Elasticity matrix types vary depending on the stress and strain field for each case. Data is obtained straight from classical FEM applications. Elasticity matrices for 1D elasticity, plane strain, plane stress and 3D elasticity are presented as follows, where E is the Young’s modulus and ν is the Poisson’s ratio.
1D elasticity:
2D elasticity, Plane Stress:
2D elasticity, Plane Strain:
3D elasticity:
The corresponding stress and strain vectors are:
1D elasticity:
2D elasticity, Plane Stress and Plane Strain:
3D elasticity:

Stiffness Matrix Assembly

The general process for the global stiffness matrix assembly, as obtained from finite element method, is shown in the following flow chart:
Figure 3. Stiffness Matrix Assembly in Finite Element Method.
There are three loops:
  • Patch loop
  • Element loop
  • Gauss Point loop
It is worth mentioning that the element loop in IGA can be avoided. In this case, the stiffness contribution of each control point pair is added directly to the stiffness matrix of the patch. The reason for this is that parameter space is local to patch rather than elements.

Therefore, an engineer accustomed to the methods of isogeometric analysis knows that the element loop can be completely avoided, leading to the following flow chart:
Figure 4. Stiffness Matrix Assembly in Isogeometric Analysis.  

Input Data

Information necessary for the whole process of analysis is given as input at this point. The information essential for the formulation of the stiffness matrix is divided into two categories:

  • Structural analysis and material
  • Computational geometry

Stiffness Matrix 1D

Let us consider a simple 1D application as an example. These applications are utilized only in truss systems with axial tension, but their significance is mostly academic. Simplifying the problem, without loss of generality, can lead to a better understanding of the principles involved. In our research career, we often address 1D analogies for the solution of complex problems and it has always been extremely helpful.
In such a case, only axial deformation for each point of the truss exists. This deformation is u(x)=u(C(ξ))=u(ξ). The respective strain matrix consists of one value:
In order to calculate the derivative of u, we must first establish a transition from physical space to parameter space:
In finite element analysis, the reverse transformation is utilized. This is why the inverse matrix [J]-1 is needed.
Special care has to be taken in order for the Jacobian to be correct. The positive direction of the axes in parameter and physical space must coincide, or the determinant of the Jacobian will be negative and the matrix [J] irreversible. Numerical integration on points of singularity, such as two points on parameter space mapped into the same point on physical space, has to be avoided as well.
After that, the deformation matrix is evaluated.
Deformation matrix produces strain values anywhere in the model, by utilizing nodal displacements.
The Stiffness matrix for the patch is evaluated as shown:
Direct integration is almost never applicable. Numerical integration is used instead, looping through all the Gauss points of a patch and their respective weights:

Stiffness Matrix 2D

2D elasticity problems have many applications in modern analysis. The logic is exactly the same as in 1D problems. The main difference is, obviously, the utilization of one more dimension. Parameter space is defined on (ξ, η) and physical space on (χ, y). Displacements per x, y at any point in the entire domain are defined as u(x,y)=u(S(ξ,η))=u(ξ,η) and v(x,y)=v(ξ,η) respectively.
The strain vector is defined as:
The transformation of a function φ between parameter and physical space yields:
Thus, the 2D Jacobian matrix can be defined as:
and the inverse mapping:
The Jacobian matrix can be evaluated as shown:
The inverse Jacobian matrix is used in stiffness matrix calculation:
The determinant of the Jacobian matrix is also required and is equal to:
In order to evaluate the stiffness matrix, integration is required.
Numerical integration procedures for ξ, η lead to integration for tensor product Gauss points.

Stiffness Matrix 3D

3D elasticity is merely the extension of 2D elasticity in all directions, with a complete stress field. Every other problem can be created by downgrading 3D problems into 2D and 1D problems.
The displacement field for each point in physical space is now defined for x, y, z by u(x,y,z)=u(S(ξ,η,ζ)=u(ξ,η,ζ), v(ξ,η,ζ), w(ξ,η,ζ) respectively. The strain field can now be defined as:
, which leads to the definition of the Jacobian matrix for 3D:
and the inverse Jacobian matrix as well:
Jacobian matrix can be calculated from the derivatives of the shape functions.
The inverse of the Jacobian matrix is:
As a result, the deformation matrix for 3D elasticity is calculated as:
The corresponding stiffness matrix is produced by integration
Numerical integration is used in 3D as well

Stiffness Matrix Examples

  Figure 7. (a) Stiffness Matrix of the model where 1470 non-zero elements are produced
and (b) its corresponding basis functions.  
    Figure 8. (a) Stiffness Matrix of the model where 1762 non-zero elements are produced
and (b) its corresponding basis functions.
In Figure 4.8, the model is analyzed for a degree of p=2 and C1 Continuity with the same 66 degrees of freedom. The interconnectivity between the elements has increased and therefore provides a better approximation of the model than the one in Figure 4.7. In this case, the configuration of the matrix requires more time, but the overlapping between the B-SPLines leads to better results.

A Stiffness Matrix coefficient is non-zero when shared support applies for the corresponding B-SPLine and extensively between Control Points.
Figure 9. (a) Stiffness Matrix of the model where 2340 non-zero elements are produced and
(b) its corresponding basis functions.
In Figure 4.9, the model is analyzed for p=3 and C2 Continuity and for the same 66 degrees of freedom.
Figure 4.11 represents the influence of the continuity in the analysis of the model in Figure 4.10. As the continuity increases, the interconnectivity of the elements affects more degrees of freedom. This leads to a better approximation of the physical problem, but also the creation of the Stiffness Matrix becomes more time-consuming. For the same polynomial order, more non-zero elements are created at the Stiffness Matrix.

This is a result of the shared support of the B-SPLine basis functions, that exists for 2p B-SPLines, which leads to overlapping between the B-SPLines and interconnectivity of the degrees of freedom in the Stiffness Matrix. Therefore, as Continuity increases, the number of non-zero elements increases as well.

  Figure 10. 3D Structure for Continuity investigation.