Stiffness matrix
The stiffness matrix is a square matrix in structural mechanics that linearly relates the vector of applied nodal forces \mathbf{F} to the vector of corresponding nodal displacements \mathbf{u} in a discretized elastic structure, expressed by the equation \mathbf{K} \mathbf{u} = \mathbf{F}, where \mathbf{K} denotes the stiffness matrix itself.[1] This matrix encapsulates the inherent resistance of the structure to deformation, with its elements representing stiffness coefficients that quantify how forces at one node influence displacements at others.[1] It forms the cornerstone of computational methods for analyzing determinate and indeterminate structures, such as trusses, beams, and frames, by enabling the systematic solution of equilibrium equations.[2] In practice, the stiffness matrix is constructed using the direct stiffness method, a matrix-based approach where individual element stiffness matrices—derived from material properties, geometry, and boundary conditions—are transformed from local to global coordinates and assembled into a global stiffness matrix that enforces overall equilibrium.[2] This assembly process accounts for degrees of freedom at nodes (e.g., two per node for 2D trusses or three for 2D frames), allowing the solution for unknown displacements at free degrees of freedom before back-substituting to find reactions and internal forces.[2] The method's efficiency stems from its adaptability to computer implementation, making it indispensable for finite element analysis (FEA) in engineering disciplines like civil, mechanical, and aerospace.[2] Key mathematical properties of the stiffness matrix include symmetry (\mathbf{K} = \mathbf{K}^T), arising from Maxwell's reciprocity theorem, which states that the displacement at one point due to a unit force at another equals the displacement at the second point due to a unit force at the first.[1] For stable structures, \mathbf{K} is positive definite, with all eigenvalues positive, ensuring a unique solution to the displacement equations; zero or negative eigenvalues indicate instability or rigid body modes.[1] These properties facilitate numerical solving techniques, such as Gaussian elimination, while the matrix's sparsity (due to localized element interactions) optimizes computational performance in large-scale simulations.[1] The development of the stiffness matrix traces its roots to 19th-century structural theory, with early formulations appearing in Alfred Clebsch's 1862 work on three-dimensional pinned trusses, building on principles from Louis Navier (1826) and Barré de Saint-Venant (1838).[3] Iterative advancements, including Heinrich Manderla's prize-winning method for secondary stresses (1880) and Otto Mohr's modifications (1892), laid groundwork for matrix approaches in the early 20th century, alongside Hardy Cross's moment distribution method (1930).[3] The modern direct stiffness method emerged in the 1950s amid aerospace demands for analyzing complex aircraft frames, with pivotal contributions from Gabriel Kron's tensor-based formulations (1944) and John H. Argyris's matrix extensions (1954–1955).[3] Its formalization as part of the finite element method occurred in the seminal 1956 paper by M. J. Turner, R. W. Clough, H. C. Martin, and L. J. Topp, which introduced displacement-based stiffness assembly for continuum structures. This evolution transformed stiffness matrix applications from discrete frames to versatile FEA tools for simulating real-world phenomena like vibrations, buckling, and nonlinear behaviors.Fundamentals
Definition and basic properties
In linear elasticity and structural mechanics, the stiffness matrix, denoted as \mathbf{K}, is a square matrix that relates the applied force vector \mathbf{F} to the resulting displacement vector \mathbf{u} through the linear equation \mathbf{F} = \mathbf{K} \mathbf{u}.[1][4] This relation captures the system's resistance to deformation under external loads, assuming small displacements and linear material behavior where Hooke's law holds.[4] The stiffness matrix exhibits several fundamental properties. It is symmetric, satisfying \mathbf{K}^T = \mathbf{K}, which follows from Maxwell's reciprocity theorem stating that the displacement at one point due to a force at another equals the displacement at the second point due to the same force at the first.[1] For stable, physically realistic systems without buckling or instability, \mathbf{K} is positive definite, meaning all its eigenvalues are positive and the quadratic form \mathbf{u}^T \mathbf{K} \mathbf{u} > 0 for any nonzero \mathbf{u}, ensuring positive strain energy storage.[1][4] However, if the system is unconstrained (e.g., free-floating), \mathbf{K} is singular with zero eigenvalues corresponding to rigid body modes, allowing arbitrary translations or rotations without deformation.[1] The stiffness matrix arises from the principle of minimum potential energy in conservative systems. The total potential energy \Pi is the sum of the strain energy U (internal elastic energy) and the potential of external forces \Omega (negative work done by loads). Equilibrium displacements minimize \Pi, and setting the first variation \delta \Pi = 0 yields the force-displacement equations; the second variation corresponds to the Hessian of U, which is precisely \mathbf{K}.[5][6] A simple illustration is a single spring with constant k, fixed at one end and loaded at the other, representing a one-degree-of-freedom system where the stiffness is the scalar K = k, so F = k u. This extends to multi-degree-of-freedom systems, such as a spring connected between two nodes with displacements u_1 and u_2, yielding the element stiffness matrix \mathbf{k} = k \begin{bmatrix} 1 & -1 \\ -1 & 1 \end{bmatrix}, relating nodal forces to relative displacements via \mathbf{f} = \mathbf{k} \mathbf{u}.[6]Physical interpretation
The stiffness matrix generalizes the scalar stiffness concept from Hooke's law, which describes the linear relationship between force and displacement for a single-degree-of-freedom system as F = k x, to multi-dimensional systems where deformations occur in multiple directions simultaneously. In this matrix form, the relationship becomes \mathbf{F} = \mathbf{K} \mathbf{u}, where \mathbf{K} captures the coupled resistance of the structure to vector displacements \mathbf{u}, reflecting material properties, geometry, and connectivity that determine how loads induce deformations across degrees of freedom.[7] The entries of the stiffness matrix provide a direct physical insight into these interactions: the element K_{ij} denotes the force developed at the i-th degree of freedom due to a unit displacement imposed at the j-th degree of freedom, with all other displacements constrained to zero. Positive diagonal terms K_{ii} indicate the inherent stiffness against direct deformation in that direction, while off-diagonal terms K_{ij} (for i \neq j) reveal coupling effects, where the sign determines whether the induced force opposes or aids the displacement, arising from the structure's geometry and boundary conditions.[8] In the context of structural equilibrium, the stiffness matrix represents the collective internal resistance of the system, such that the product \mathbf{K} \mathbf{u} generates the nodal forces necessary to counteract applied external loads \mathbf{F} in the equation \mathbf{K} \mathbf{u} = \mathbf{F}. This formulation ensures that the deformed configuration balances internal elastic restoring forces against external actions, maintaining static equilibrium. For instance, in a simple 2D frame composed of beams, the diagonal entries embody direct stiffness contributions like axial rigidity along member lengths and flexural resistance to transverse loads, whereas off-diagonal entries account for geometric coupling, such as how a horizontal displacement at a joint induces vertical shear forces due to inclined or connected members.[8][9] The positive definiteness of the stiffness matrix further underscores its physical reliability, guaranteeing a unique and stable displacement response to any load without spontaneous deformation.[1]Formulation in mechanics
Discrete systems
In discrete systems, such as lumped-parameter models, the stiffness matrix assembles from individual element contributions to relate nodal forces to displacements across the structure.[6] These models approximate continuous structures by concentrating mass and stiffness at discrete nodes connected by elements like springs or bars, enabling efficient computation of static and dynamic responses.[10] The degrees of freedom in such systems correspond to the possible nodal displacements and associated forces. In one-dimensional models, like axial spring chains, each node has a single displacement degree of freedom along the axis, with forces acting in the same direction.[6] For two-dimensional truss structures, each node typically has two degrees of freedom—horizontal and vertical displacements—with corresponding force components.[10] In three dimensions, nodes in space trusses exhibit three translational degrees of freedom per node.[10] Lumped models assemble the global stiffness matrix by superimposing element-level matrices, ensuring compatibility at shared nodes. For a mass-spring chain, consider two springs in series connecting three nodes, with spring constants k_1 and k_2. The local stiffness matrix for each spring is \begin{bmatrix} k & -k \\ -k & k \end{bmatrix}, relating end forces to relative displacements.[6] Assembly yields the 3×3 global stiffness matrix: K = \begin{bmatrix} k_1 & 0 & -k_1 \\ 0 & k_2 & -k_2 \\ -k_1 & -k_2 & k_1 + k_2 \end{bmatrix}, where the off-diagonal terms reflect force equilibrium at nodes.[6] In truss structures, element stiffness matrices account for geometry and orientation. For a two-dimensional bar element with cross-sectional area A, modulus E, and length L at angle \theta to the global x-axis, the element stiffness matrix is K_e = \frac{AE}{L} \mathbf{T}^T \begin{bmatrix} 1 & -1 \\ -1 & 1 \end{bmatrix} \mathbf{T}, where \mathbf{T} is the transformation matrix incorporating \cos\theta and \sin\theta to rotate between local and global coordinates.[10] Explicitly, K_e = \frac{AE}{L} \begin{bmatrix} c^2 & cs & -c^2 & -cs \\ cs & s^2 & -cs & -s^2 \\ -c^2 & -cs & c^2 & cs \\ -cs & -s^2 & cs & s^2 \end{bmatrix}, with c = \cos\theta and s = \sin\theta.[10] For a two-bar truss, the global matrix forms by adding contributions from each element's K_e, expanded to the full degrees of freedom and condensed for connectivity.[10] Boundary conditions incorporate supports by modifying the stiffness matrix or the force vector. For fixed supports (zero displacement), partition the matrix by removing rows and columns corresponding to constrained degrees of freedom, reducing the system to free nodes.[6] For prescribed non-zero displacements, adjust the force vector by subtracting stiffness terms multiplied by the known displacements, maintaining equilibrium.[6] This ensures the solved displacements satisfy structural constraints without altering the matrix symmetry or positive definiteness.[6]Continuous systems via finite elements
In the finite element method (FEM), stiffness matrices for continuous systems are derived by discretizing the domain into finite elements and approximating the governing partial differential equations through variational principles. This approach transforms the continuous problem into a system of algebraic equations where the element stiffness matrix relates nodal displacements to forces, enabling the analysis of structures like beams, plates, and solids under elastic deformation.[11] The derivation begins with the weak form of the equilibrium equations, obtained by multiplying the strong form by a test function and integrating over the domain, often using integration by parts to reduce derivative orders and incorporate natural boundary conditions. For linear elasticity, the weak form expresses the principle of virtual work: the internal virtual work equals the external virtual work, leading to an integral equation over the element volume. The Galerkin method then approximates both the solution and test functions within this weak form using the same basis, ensuring symmetry in the resulting matrices for self-adjoint problems and promoting convergence.[11][11] Within each element, displacements are approximated as \mathbf{u} = \mathbf{N} \mathbf{d}, where \mathbf{N} is the matrix of shape functions and \mathbf{d} contains the nodal displacement values. Shape functions, typically polynomials like linear or quadratic, interpolate the field variables and ensure inter-element continuity (e.g., C^0 continuity for displacements). Strains are then derived as \boldsymbol{\epsilon} = \mathbf{B} \mathbf{d}, with the strain-displacement matrix \mathbf{B} formed from derivatives of the shape functions. The material constitutive relation links stresses to strains via \boldsymbol{\sigma} = \mathbf{D} \boldsymbol{\epsilon}, where \mathbf{D} is the material matrix (e.g., for isotropic elasticity, involving Young's modulus and Poisson's ratio). Substituting into the weak form yields the element stiffness matrix as \mathbf{K}_e = \int_{V_e} \mathbf{B}^T \mathbf{D} \mathbf{B} \, dV, which represents the resistance to deformation within the element.[11][11][11] For complex geometries, isoparametric elements map the physical element coordinates to a simple reference (parent) domain using the same shape functions for both geometry and displacements: \mathbf{x} = \sum N_i \mathbf{x}_i and \mathbf{u} = \sum N_i \mathbf{d}_i. This mapping involves the Jacobian matrix \mathbf{J} for coordinate transformation, adjusting the integration to the reference space as \mathbf{K}_e = \int_{\hat{V}_e} \mathbf{B}^T \mathbf{D} \mathbf{B} \det(\mathbf{J}) \, d\hat{V}, where d\hat{V} = d\xi \, d\eta \, d\zeta in natural coordinates (e.g., \xi \in [-1,1]). This formulation allows flexible handling of curved boundaries while maintaining computational efficiency through numerical quadrature.[11] A simple example is the 1D bar element under axial loading, with length L, cross-sectional area A, and Young's modulus E. Using linear shape functions N_1 = 1 - \xi and N_2 = \xi (where \xi = x/L), the displacement is u(x) = (1 - \xi) d_1 + \xi d_2, and the strain is constant: \epsilon = (d_2 - d_1)/L. The B matrix is \mathbf{B} = [-1/L, 1/L], and with \mathbf{D} = E, the stiffness matrix integrates to \mathbf{K}_e = \frac{EA}{L} \begin{bmatrix} 1 & -1 \\ -1 & 1 \end{bmatrix}. This 2x2 matrix captures the element's axial stiffness, with off-diagonals indicating load transfer between nodes. The physical interpretation aligns with the matrix embodying resistance to relative nodal displacements, analogous to spring constants in discrete systems.[11]Specific cases
Poisson equation
The Poisson equation, a fundamental elliptic partial differential equation, models scalar potential fields such as temperature in heat conduction or electrostatic potential, and is given by -\nabla \cdot (k \nabla u) = f in a domain \Omega, where k > 0 is a diffusion coefficient, u is the unknown scalar field, and f is a source term.[11] To apply the finite element method (FEM), the equation is first reformulated in its weak form by multiplying by a test function \phi \in H^1_0(\Omega) (the Sobolev space of functions with square-integrable first derivatives vanishing on the boundary for homogeneous Dirichlet conditions) and integrating by parts: \int_\Omega k \nabla \phi \cdot \nabla u \, d\Omega = \int_\Omega f \phi \, d\Omega + \int_{\partial \Omega_N} g \phi \, d\Gamma, where \partial \Omega_N denotes the Neumann boundary portion and g is the specified flux.[12] This bilinear form on the left-hand side ensures well-posedness under suitable coercivity conditions on k.[11] In the Galerkin FEM discretization, the solution is approximated as u_h = \sum_{j=1}^N u_j \phi_j, where \{\phi_j\} are basis functions spanning a finite-dimensional subspace V_h \subset H^1_0(\Omega), such as continuous piecewise polynomials on a mesh. Substituting into the weak form yields the discrete system K \mathbf{u} = \mathbf{F}, where \mathbf{u} collects the nodal coefficients u_j, the load vector has entries F_j = \int_\Omega f \phi_j \, d\Omega + \int_{\partial \Omega_N} g \phi_j \, d\Gamma (involving a mass matrix M_{ij} = \int_\Omega \phi_i \phi_j \, d\Omega when f is projected), and the stiffness matrix K is symmetric positive definite with entries K_{ij} = \int_\Omega k \nabla \phi_i \cdot \nabla \phi_j \, d\Omega.[12] This matrix encodes the diffusion operator and is the core of the FEM approximation for the Poisson problem, with its sparsity arising from local support of the basis functions.[11] For two-dimensional problems on triangular meshes, linear Lagrange basis functions are commonly used, where on each element T_e with vertices \mathbf{r}_1, \mathbf{r}_2, \mathbf{r}_3 and area A_e, the local basis \phi_i|_{T_e} = \frac{1}{2A_e} (a_i + b_i x + c_i y) has constant gradients \nabla \phi_i = \frac{1}{2A_e} (b_i, c_i). Assuming constant k for simplicity, the element stiffness matrix K_e is then explicit: K_e = \frac{k}{4 A_e} \begin{pmatrix} |\mathbf{r}_2 - \mathbf{r}_3|^2 & (\mathbf{r}_2 - \mathbf{r}_3) \cdot (\mathbf{r}_3 - \mathbf{r}_1) & (\mathbf{r}_2 - \mathbf{r}_3) \cdot (\mathbf{r}_1 - \mathbf{r}_2) \\ (\mathbf{r}_2 - \mathbf{r}_3) \cdot (\mathbf{r}_3 - \mathbf{r}_1) & |\mathbf{r}_3 - \mathbf{r}_1|^2 & (\mathbf{r}_3 - \mathbf{r}_1) \cdot (\mathbf{r}_1 - \mathbf{r}_2) \\ (\mathbf{r}_2 - \mathbf{r}_3) \cdot (\mathbf{r}_1 - \mathbf{r}_2) & (\mathbf{r}_3 - \mathbf{r}_1) \cdot (\mathbf{r}_1 - \mathbf{r}_2) & |\mathbf{r}_1 - \mathbf{r}_2|^2 \end{pmatrix}, which is assembled into the global K by summing contributions from shared nodes.[12] This form highlights the dependence on element geometry and ensures second-order accuracy in the H^1-norm for smooth solutions.[11] Boundary conditions influence the stiffness matrix differently based on type. Essential (Dirichlet) conditions, prescribing u = \bar{u} on \partial \Omega_D, are enforced by restricting V_h to functions satisfying them (e.g., setting corresponding rows and columns of K to zero except the diagonal to identity, and adjusting \mathbf{F}), preserving symmetry.[12] Natural (Neumann or Robin) conditions, specifying flux k \nabla u \cdot \mathbf{n} = g or a convective form, appear only in \mathbf{F} via boundary integrals and do not modify K directly, allowing the matrix to remain unchanged while incorporating inhomogeneous data.[11]Elasticity and structural problems
In solid mechanics, the stiffness matrix plays a central role in the finite element formulation of linear elasticity problems, which govern the deformation of deformable bodies under external loads. The fundamental equilibrium equation is \nabla \cdot \boldsymbol{\sigma} = \mathbf{f}, where \boldsymbol{\sigma} is the stress tensor, \mathbf{f} is the body force vector, and the stress-strain relation follows Hooke's law \boldsymbol{\sigma} = \mathbf{D} \boldsymbol{\varepsilon}, with the infinitesimal strain tensor defined as \boldsymbol{\varepsilon} = \frac{1}{2} (\nabla \mathbf{u} + (\nabla \mathbf{u})^T). Here, \mathbf{u} is the displacement vector, and \mathbf{D} is the constitutive (elasticity) matrix that encapsulates material properties such as Young's modulus E and Poisson's ratio \nu. This vector-valued framework distinguishes elasticity from scalar problems, as it accounts for coupled multi-directional deformations and stress components.[11] Within the finite element method, the element stiffness matrix \mathbf{K}_e for linear elasticity is derived from the principle of virtual work and expressed as the volume integral \mathbf{K}_e = \int_V \mathbf{B}^T \mathbf{D} \mathbf{B} \, dV, where \mathbf{B} is the strain-displacement matrix relating strains to nodal displacements. The matrix \mathbf{B} incorporates the gradients of the shape functions and explicitly includes terms for normal strains (\varepsilon_{xx}, \varepsilon_{yy}, \varepsilon_{zz}) and shear strains (\varepsilon_{xy}, \varepsilon_{xz}, \varepsilon_{yz}), ensuring symmetry in 2D and 3D formulations; for instance, in 2D, \mathbf{B} is a $3 \times 8 matrix for a 4-node element, with rows corresponding to \varepsilon_{xx} = \sum \frac{\partial N_i}{\partial x} u_i, \varepsilon_{yy} = \sum \frac{\partial N_i}{\partial y} v_i, and \varepsilon_{xy} = \frac{1}{2} \sum \left( \frac{\partial N_i}{\partial y} u_i + \frac{\partial N_i}{\partial x} v_i \right), where N_i are the shape functions and u_i, v_i are nodal displacements. The integral is typically evaluated numerically via Gauss quadrature to handle arbitrary geometries.[13] For two-dimensional approximations, plane stress and plane strain assumptions simplify the 3D problem by reducing the dimensionality while preserving key mechanics. In plane stress, applicable to thin plates where out-of-plane stresses vanish (\sigma_{zz} = \sigma_{xz} = \sigma_{yz} = 0), the reduced constitutive matrix is \mathbf{D} = \frac{E}{1 - \nu^2} \begin{bmatrix} 1 & \nu & 0 \\ \nu & 1 & 0 \\ 0 & 0 & \frac{1 - \nu}{2} \end{bmatrix}, relating in-plane stresses \{\sigma_{xx}, \sigma_{yy}, \sigma_{xy}\} to strains. Conversely, plane strain, suitable for long bodies constrained in the z-direction (\varepsilon_{zz} = \varepsilon_{xz} = \varepsilon_{yz} = 0), uses \mathbf{D} = \frac{E}{(1 + \nu)(1 - 2\nu)} \begin{bmatrix} 1 - \nu & \nu & 0 \\ \nu & 1 - \nu & 0 \\ 0 & 0 & \frac{1 - 2\nu}{2} \end{bmatrix}. These forms adjust the effective stiffness to reflect the assumptions, enabling efficient 2D simulations without full 3D meshing.[13] A representative example is the 4-node quadrilateral element for 2D linear elasticity, commonly used in plane stress or strain analyses due to its bilinear interpolation and ability to model irregular shapes via isoparametric mapping. The shape functions in natural coordinates (\xi, \eta) are N_i = \frac{1}{4} (1 + \xi_i \xi)(1 + \eta_i \eta), leading to a \mathbf{B} matrix that varies within the element. For a rectangular element of width $2a and height $2b (with coordinates aligned), the position-dependent \mathbf{B} matrix is \mathbf{B} = \begin{bmatrix} -\frac{b - y}{4ab} & 0 & \frac{b - y}{4ab} & 0 & \frac{b + y}{4ab} & 0 & -\frac{b + y}{4ab} & 0 \\ 0 & -\frac{a - x}{4ab} & 0 & -\frac{a + x}{4ab} & 0 & \frac{a + x}{4ab} & 0 & \frac{a - x}{4ab} \\ -\frac{a - x}{4ab} & -\frac{b - y}{4ab} & -\frac{a + x}{4ab} & \frac{b - y}{4ab} & \frac{a + x}{4ab} & \frac{b + y}{4ab} & \frac{a - x}{4ab} & -\frac{b + y}{4ab} \end{bmatrix}, where x, y are local coordinates within [-a, a] \times [-b, b]; the full \mathbf{K}_e = t \int_{-a}^{a} \int_{-b}^{b} \mathbf{B}^T \mathbf{D} \mathbf{B} \, dx \, dy can be computed analytically for uniform \mathbf{D} to yield an explicit $8 \times 8 matrix, or numerically via Gauss quadrature to account for the varying \mathbf{B}. This element captures shear locking if not integrated properly, typically requiring reduced (1-point) or selective integration.[14]Assembly and computation
Element-level construction
The construction of the local stiffness matrix for a finite element begins with evaluating the integral derived from the weak form of the governing equations, which relates the element's contributions to the overall system through shape functions and material properties. This process is performed in the element's local coordinate system to isolate computations before global assembly. For linear elastic materials, the element stiffness matrix \mathbf{k}^e is computed as \mathbf{k}^e = \int_{V^e} \mathbf{B}^T \mathbf{D} \mathbf{B} \, dV, where \mathbf{B} is the strain-displacement matrix, \mathbf{D} is the constitutive matrix, and the integration is over the element volume V^e.[15][16] This integral is typically evaluated numerically using Gauss-Legendre quadrature, which approximates the integrand at specific Gauss points \xi_i with weights w_i, ensuring exact integration for polynomials up to a certain degree matching the element's order. For a 2D quadrilateral element, a 2x2 Gauss quadrature rule (four points) is often sufficient for bilinear interpolation, transforming the integral to \mathbf{k}^e \approx \sum_{i=1}^{n_g} w_i \mathbf{B}(\xi_i)^T \mathbf{D} \mathbf{B}(\xi_i) \det(\mathbf{J}(\xi_i)), where n_g is the number of Gauss points and \mathbf{J} is the Jacobian of the coordinate transformation.[15][17][18] In isoparametric elements, which use the same shape functions for geometry and displacements, the physical coordinates are mapped from a reference (parent) element via \mathbf{x} = \sum N_i(\xi) \mathbf{x}_i, requiring the Jacobian matrix \mathbf{J} = \partial \mathbf{x}/\partial \xi to account for distortions in curved or irregular elements; its determinant scales the differential volume in the integral. This transformation ensures accurate representation of arbitrary shapes while maintaining computational efficiency.[19][20][21] While the above focuses on linear problems, extensions to material nonlinearity involve the tangent stiffness matrix, which updates \mathbf{D} to the current secant or tangent modulus at each iteration of a nonlinear solver, preserving the integral form but requiring re-evaluation per load step.[22][23] As a specific example, consider a bilinear quadrilateral element in 2D for the Poisson equation -\nabla^2 u = f, where the stiffness matrix simplifies to \mathbf{k}^e = \int_{\Omega^e} \nabla \mathbf{N}^T \nabla \mathbf{N} \, d\Omega with identity constitutive matrix. For a rectangular element of width a and height b, the shape function derivatives in local coordinates yield explicit entries after 2x2 Gauss quadrature; for instance, the (1,1) entry is \frac{1}{3} \left( \frac{a}{b} + \frac{b}{a} \right), but numerical evaluation at points like (\pm 1/\sqrt{3}, \pm 1/\sqrt{3}) with weights 1 handles general cases efficiently.[24]Global matrix assembly
The global stiffness matrix in the finite element method is constructed by superimposing the contributions from each element's stiffness matrix, based on the mesh connectivity, to form a single matrix that governs the entire discretized domain.[25] This assembly process, central to the direct stiffness method, ensures kinematic compatibility across element boundaries and force equilibrium at shared nodes.[26] Connectivity between elements is defined by node numbering in the mesh, which maps local degrees of freedom in each element to global indices. During assembly, a scatter-gather operation distributes the local stiffness matrix entries to their corresponding global positions and sums overlapping contributions; specifically, for an element's local entry K_e(m,n), the global update is K_\text{global}(i,j) += K_e(m,n), where i and j are the global node numbers linked to local indices m and n via the element's connectivity array.[7] The element stiffness matrices themselves are typically obtained from numerical integration over the element domain. The assembled global stiffness matrix exhibits a sparse structure, with non-zero entries limited to degrees of freedom connected through the mesh topology, reflecting the localized interactions in the physical system.[26] To enhance computational efficiency in storage and solution phases, node renumbering techniques, such as the Cuthill-McKee algorithm, are applied to minimize the matrix bandwidth—the maximum difference in node indices for non-zero off-diagonal terms—thereby reducing fill-in during factorization.[27] The global load vector is assembled analogously, by scattering and summing element-level force contributions (from distributed loads or body forces) to the appropriate global degrees of freedom based on connectivity.[28] Consider a simple one-dimensional example: a bar discretized into three equal-length elements with four nodes (1, 2, 3, 4), assuming uniform axial stiffness k = \frac{AE}{L} for each element, where A is the cross-sectional area, E is the modulus, and L is the element length. Each element's local stiffness matrix is K_e = k \begin{bmatrix} 1 & -1 \\ -1 & 1 \end{bmatrix}, with local nodes mapped as follows: element 1 (nodes 1–2), element 2 (nodes 2–3), element 3 (nodes 3–4).- For element 1, scatter to global rows/columns 1 and 2: adds k to K_\text{global}(1,1) and K_\text{global}(2,2), -k to K_\text{global}(1,2) and K_\text{global}(2,1).
- For element 2, scatter to global rows/columns 2 and 3: adds k to K_\text{global}(2,2) and K_\text{global}(3,3), -k to K_\text{global}(2,3) and K_\text{global}(3,2); note the summation at row/column 2.
- For element 3, scatter to global rows/columns 3 and 4: adds k to K_\text{global}(3,3) and K_\text{global}(4,4), -k to K_\text{global}(3,4) and K_\text{global}(4,3); summation occurs at row/column 3.