PaperlayerView exact revision on arXiv

Source-linked AI summary

Isogeometric analysis: an overview and computer implementation aspects

Vinh Phu Nguyen, Stéphane P. A. Bordas, Timon Rabczuk

arXiv:1205.2129v2math.NAcs.MS

TL;DR

IGA uses CAD-generated basis functions to discretize geometry and unknown fields, while addressing challenges including faceted contact surfaces and watertight NURBS geometries. The paper presents a Matlab implementation spanning multiple dimensions and mechanics problems, with examples demonstrating smooth contact surfaces, exact-solution capture, and B-spline/NURBS flexibility, while noting that further development is needed.

  • Problem

    Conventional geometry discretizations can produce faceted contact surfaces, while NURBS have difficulty producing watertight geometries.

  • Method

    The paper presents IGA using CAD-generated basis functions for geometry and unknown fields, with associated basis-function construction, refinement, and Matlab implementation details.

  • Results

    The implementation covers one-, two-, and three-dimensional isogeometric finite element analysis, including smooth contact surfaces and examples capturing exact solutions with B-splines/NURBS.

  • Takeaways & Limitations

    IGA provides a unified analysis framework using the same CAD-generated basis functions for geometry and unknown fields, with demonstrated flexibility across example problems.

  • Takeaways & Limitations

    Further development is still needed beyond the paper's preliminary concepts and implementation details, and NURBS remain unable to produce watertight geometries.

Abstract

from arXiv · show

Isogeometric analysis (IGA) represents a recently developed technology in computational mechanics that offers the possibility of integrating methods for analysis and Computer Aided Design (CAD) into a single, unified process. The implications to practical engineering design scenarios are profound, since the time taken from design to analysis is greatly reduced, leading to dramatic gains in efficiency. The tight coupling of CAD and analysis within IGA requires knowledge from both fields and it is one of the goals of the present paper to outline much of the commonly used notation. In this manuscript, through a clear and simple Matlab implementation, we present an introduction to IGA applied to the Finite Element (FE) method and related computer implementation aspects. Furthermore, implemen- tation of the extended IGA which incorporates enrichment functions through the partition of unity method (PUM) is also presented, where several examples for both two-dimensional and three-dimensional fracture are illustrated. The open source Matlab code which accompanies the present paper can be applied to one, two and three-dimensional problems for linear elasticity, linear elastic fracture mechanics, structural mechanics (beams/plates/shells including large displacements and rotations) and Poisson problems with or without enrichment. The Bezier extraction concept that allows FE analysis to be performed efficiently on T-spline geometries is also incorporated. The article includes a summary of recent trends and developments within the field of IGA.

1. Introduction

IGA unifies CAD geometry representation and finite-element analysis by using CAD-generated basis functions for both geometry and unknown fields. The paper surveys its applications, implementation strategies, alternatives, and remaining computational and geometric limitations.

  • Underlying concepts of isogeometric analysis: NURBS can represent conic and circular sections exactly, and IGA uses them as both geometry and analysis discretisations.This CAD-analysis coupling maintains geometry throughout analysis and avoids geometrical error at all stages.
  • Applications: IGA has been applied across structural, solid, fluid, contact, optimisation, and experimental mechanics, with applications also extending to meshfree methods.Reported examples include smooth contact surfaces, plate and shell formulations, fluid-structure interaction, and NURBS-based digital image correlation.
  • Applications: Smooth NURBS bases support rotation-free plate and shell elements, reduce shear-locking relative to standard FE elements, and enable high-order-continuity formulations for higher-order PDEs.Multi-patch rotation-free elements require special treatment at patch boundaries where basis functions are C0 continuous.
  • Applications: For structural vibration, k-refinement can provide more robust and accurate frequency spectra than higher-order FE p-methods by eliminating optical branches.With lumped mass matrices, however, accuracy is limited to second order for any basis order because high-order isogeometric lumped mass matrices are unavailable.
  • Shortcomings of NURBS and alternatives: NURBS have computational and geometric limitations, including globally coupled refinement, inefficient adaptivity, difficulty producing watertight geometries, and complications from tensor-product structure.T-splines address these deficiencies with single-patch watertight geometries and local refinement, although linear independence is not guaranteed on generic T-meshes and adaptive knot insertion can be complex in 3D.
  • Contributions and outline: The paper provides Matlab implementations for one-, two-, and three-dimensional IGA, structural and fracture problems, visualisation, boundary conditions, and extended IGA with discontinuity enrichment.Bezier extraction is included to support efficient finite-element analysis on T-spline geometries, while the authors state that substantial research remains before unified design and analysis is fully realised.

2. A brief introduction to B-splines/NURBS

B-splines and NURBS represent geometry through parametric functions and provide the basis for IGA discretization. Their knot vectors, continuity properties, and refinement capabilities support CAD-compatible analysis and boundary-condition treatment.

  • Parametric representation: Parametric functions map parameter-space coordinates to physical geometry, allowing curves, surfaces, and solids to be represented and geometrically sampled.B-splines and NURBS use this parametric representation to sweep out the desired shape.
  • Knot vectors: A knot vector is an ordered set of parameter values that partitions parametric space into knot spans and defines basis-function order and multiplicity.Open knot vectors commonly repeat the first and last knots p + 1 times.
  • Basis-function properties: Open-knot B-splines interpolate their start and end points, which facilitates imposing boundary conditions in analysis.Geometrically, the curve coincides with its start and end control points.
  • Basis-function properties: B-spline basis functions form a partition of unity, remain non-negative, are linearly independent, and have compact support over p + 1 knot spans.Continuity across a knot ξi is Cp−mi, where mi is its multiplicity; general basis functions are not interpolatory at control points.
  • Basis-function properties: Knot multiplicity controls continuity: quadratic functions are C1 elsewhere but attain C0 continuity at the repeated knot ξ = 4.This continuity control is useful for modeling discontinuities such as cracks or material interfaces.
  • Refinement algorithms: B-spline discretizations support knot insertion, degree elevation, and other refinement algorithms for modeling intricate geometries and capturing rapid field variations.Knot insertion is analogous to h-refinement, while degree elevation is analogous to p-refinement in conventional FEM.

3. NURBS as a basis for analysis: isogeometric finite element formulation

IGA uses CAD-generated B-spline and NURBS bases to represent both geometry and unknown fields, connecting CAD and analysis while retaining exact geometry. The formulation organizes index, parametric, physical, and parent spaces for finite-element analysis.

  • 3.1.1. Index space: Index space is generated from knot vectors, but this implementation retains only elements with non-zero parametric area and therefore does not require index space.
  • 3.1.2. Parametric space: Parametric space consists of non-zero knot intervals that define element domains and can be normalized to a unit interval, square, or cube.
  • 3.1.3. Physical space: B-spline and NURBS mappings transform parametric coordinates into physical coordinates, with control grids defining connectivity and interior control points generally non-interpolatory.
  • 3.1.4. Parent space: Parent space is introduced as [−1, 1]^dp so numerical integration routines can operate on a standard domain.
  • 3.2. Isogeometric formulation: IGA uses CAD-generated basis functions to discretize both geometry and unknown fields, combining CAD and analysis into a unified process.
  • 3.2. Isogeometric formulation: Unlike conventional discretizations, IGA maintains exact CAD geometry throughout analysis rather than introducing geometrical approximation error.
  • 3.3. Isogeometric discretisation: The field approximation uses NURBS basis functions and control variables, whose coefficients are generally non-interpolatory at control points.
  • 3.3. Isogeometric discretisation: Element-level implementation uses parent-coordinate mappings, local-to-global connectivity, and nen = (p + 1)^dp non-zero basis functions per element.

3.5. One-dimensional IGA formulation

The one-dimensional formulation constructs a Galerkin finite-dimensional problem with NURBS trial and test spaces, assembles element contributions, and verifies the implementation against exact solutions. The strong-gradient example shows that continuity can be adapted locally to capture sharp features.

  • 3.5. One-dimensional IGA formulation: The one-dimensional problem partitions the boundary into Dirichlet and Neumann parts, while explicitly excluding Robin boundary conditions.
  • 3.5. One-dimensional IGA formulation: The Galerkin formulation uses NURBS-generated finite-dimensional trial and test spaces, with the solution and weighting functions constructed from the same basis.
  • 3.5. One-dimensional IGA formulation: Dirichlet data is incorporated through the finite-dimensional boundary term, and the discrete weak form becomes a matrix system assembled from element stiffness matrices and force vectors.
  • 3.5.1. Assembly of system of equations: The implementation assembles global stiffness and force arrays by looping over non-zero knot intervals, integrating at Gauss points, and inserting element entries through connectivity.
  • 3.5.1. Assembly of system of equations: The computed quadratic NURBS solutions for two and four elements converge to the exact solution, verifying the presented IGA implementation.
  • 3.5.1. Assembly of system of equations: For a sharp peak at x = 0.5, smooth high-order bases are unsuitable, while linear elements with C0 continuity at the peak perform better.
  • 3.5.1. Assembly of system of equations: Inserting knots at x = 0.5 and nearby locations increases local continuity control and captures the sharp peak more accurately.
  • 3.5.1. Assembly of system of equations: With 16 cubic elements and 25 control points, the exact strong-gradient solution is well captured, demonstrating flexible control of order and continuity.

4. Elasticity: two-dimensional implementation

The two-dimensional elasticity implementation applies the Galerkin method with NURBS shape functions for displacement and its variations. It defines element strains, connectivity, and data structures for assembling the elasticity system.

  • 4. Elasticity: two-dimensional implementation: The resulting procedure highlights differences from conventional finite-element discretizations while retaining familiar elasticity notation and assembly structure.
  • 4. Elasticity: two-dimensional implementation: The implementation formulates two-dimensional linear elasticity on a domain whose boundary separates prescribed displacements from prescribed tractions.
  • 4. Elasticity: two-dimensional implementation: The Galerkin discretization uses the same NURBS shape functions for displacement and displacement variations.
  • 4. Elasticity: two-dimensional implementation: Each two-dimensional control point carries two unknown displacement components, uA = {uxA, uyA}, with analogous modifications for three dimensions.
  • 4. Elasticity: two-dimensional implementation: The strain-displacement matrix is constructed from physical derivatives of the NURBS shape functions, and the resulting equations produce the elasticity system.
  • 4. Elasticity: two-dimensional implementation: For the example surface, basis orders p = 1 and q = 2 with 8 control points define the two-dimensional discretization.
  • 4. Elasticity: two-dimensional implementation: Element basis functions and global indices are specified explicitly, while control-point and connectivity matrices store geometry and assembly data.

B11 B21 B12 B22 B13 B23 B14 B24

The implementation covers element integration, boundary-force assembly, direct and least-squares Dirichlet enforcement, and the associated limitations of simple control-point treatment. NURBS boundary structure supports localized enforcement, while general boundary data requires specialized methods.

  • Assembly and elasticity implementation: Element stiffness matrices are computed from basis-function derivatives, the strain-displacement matrix, material matrix, and Jacobian factors at Gauss points.
  • Assembly and elasticity implementation: The global stiffness matrix is assembled using element connectivity, with each control point contributing two positions in the global displacement vector.
  • Boundary-force assembly: A boundary mesh enables external traction integration using the one-dimensional assembly procedure.
  • Dirichlet boundary conditions: Homogeneous Dirichlet conditions are imposed by setting corresponding control variables to zero, while uniform inhomogeneous conditions use prescribed values.
  • Dirichlet boundary conditions: For open knot vectors, corner control-point values can directly impose Dirichlet conditions because the NURBS basis satisfies the Kronecker delta property there.
  • Dirichlet boundary conditions: Interior prescribed displacements and non-uniform boundary conditions require specialized treatments rather than simple direct control-point assignment.
  • Least-squares boundary enforcement: The least-squares method determines boundary control-point parameters by fitting prescribed displacements at collocation points and solving Aq = b.
  • Least-squares boundary enforcement: The Matlab procedure distributes collocation points uniformly in parameter space, and the paper notes that their influence on accuracy is discussed elsewhere.

5. Extended isogeometric finite element method

Extended IGA uses partition-of-unity enrichment to represent discontinuities such as cracks, holes, and inclusions independently of the discretization. XIGA combines NURBS bases with Heaviside and crack-tip branch enrichments, while retaining a standard Galerkin system and requiring only modest implementation changes.

  • PUM-based XIGA models internal boundaries independently of the discretization, allowing crack growth without remeshing.
  • The extended displacement approximation adds enriched degrees of freedom to standard NURBS degrees of freedom for traction-free crack modeling.
  • Control points are enriched according to support topology: S_c contains supports cut by the crack, while S_f contains supports containing the crack tip.
  • The formulation uses a Heaviside function for the crack discontinuity and branch functions spanning the crack-tip displacement field.
  • The enriched Galerkin system retains K u = f, with an enlarged strain-displacement matrix combining standard and enriched contributions.
  • The implementation supports extended IGA for cracks and material interfaces, while knot insertion is mainly suitable for cracks and can become tedious with unknown trajectories.

6. MIGFEM- A Matlab IGA (X)FEM code

MIGFEM is an open-source Matlab implementation of IGA and extended IGA spanning multiple dimensions, physical models, refinements, and structural formulations. Its mesh and assembly procedures adapt CAD data into element connectivity and otherwise remain close to conventional FEM workflows.

  • MIGFEM supports one-, two-, and three-dimensional linear elasticity, extended crack and interface modeling, geometrically nonlinear solids, and structural mechanics.
  • The code provides global h-, p-, and k-refinement, 2D/3D stationary traction-free cracks and interfaces, visualization, and several treatments of inhomogeneous Dirichlet conditions.
  • Additional capabilities include compatible multi-patch 2D formulations, Bézier extraction for T-splines, structural elements, time integration schemes, and Rhino3d surface extraction.
  • Unlike conventional FEM inputs, IGA receives CAD data consisting of knot vectors, control points, and basis-function orders, from which element connectivity is constructed.
  • The 2D mesh generator derives element ranges and connectivity from knot vectors and basis orders, then maps each element to its parametric coordinates.
  • NURBS basis routines compute first and second derivatives in 1D, 2D, and 3D, with MEX implementations used to improve performance.
  • IGA assembly is nearly identical to FEM, differing mainly through parameter-space element integration and the additional parent-to-parametric mapping.
  • The code reuses finite-element visualization techniques for IGA results and includes a simpler implementation approach than the referenced alternative.

7. Structural mechanics

NURBS continuity enables direct rotation-free Kirchhoff plate and related structural formulations. The implementation uses transverse control-point deflections, second derivatives, and paired control-point constraints to impose rotations and symmetry.

  • High-order NURBS/T-spline continuity makes rotation-free thin beam, plate, and shell elements direct to implement.
  • In the Kirchhoff plate formulation, both geometry and deflection use NURBS, with one transverse displacement unknown at each control point.
  • The presented plate model assumes isotropic elasticity and uses a displacement-curvature matrix containing second derivatives of the shape functions.
  • Clamped rotations are imposed by fixing two rows of boundary control-point deflections because those control points define the boundary tangent.
  • For symmetry lines, zero rotation is enforced by constraining the deflections of two control-point rows to have the same value.
  • The implementation computes spatial first and second derivatives before forming and assembling the plate element stiffness matrix.

8. Verification examples

The verification suite applies IGA and extended IGA to elasticity, shells, and 2D/3D fracture examples. Results show accurate stress concentration capture, comparisons with reference or conventional formulations, and qualitative agreement where a full 3D SIF evaluation was not performed.

  • The verification examples cover 2D and 3D linear elasticity and fracture, including a plate with a hole, pinched cylinder, cracked plate, and 3D mode I crack.
  • For the plate with a circular hole, the numerical σxx field captures the stress concentration at (R, 3π/2) and remains smooth throughout.
  • The pinched-cylinder example uses a tri-quadratic NURBS solid mesh and analyzes one-eighth of the geometry through symmetry.
  • The pinched-cylinder solid-NURBS example illustrates 3D analysis, although the authors state that shell elements would model this shell-like structure more accurately.
  • The edge-cracked plate compares XFEM and XIGA using an interaction-integral computation of the mode I stress intensity factor.
  • The crack study evaluates linear and cubic B-spline bases on three meshes against a reference SIF of I = 1.6118.
  • The interaction-integral weight function equals unity near the crack tip and vanishes on an outer contour, using bilinear Lagrange shape functions.
  • The 3D fracture example demonstrates extended IGA implementation and surface-based Dirichlet boundary enforcement, but no 3D SIF computation was performed.

9. Conclusion

The paper presents a Matlab implementation of isogeometric finite element analysis across dimensions and mechanics problems, including enrichment and nonlinear capabilities. It also identifies unresolved challenges in developing IGA.

  • The Matlab implementation covers one-, two-, and three-dimensional isogeometric finite element analysis.
  • The implementation is intended to teach IGA concepts and differences from conventional finite element implementations through explicit source code.
  • The extended isogeometric formulation supports two- and three-dimensional linear elastic fracture analysis directly from CAD data.
  • The code includes geometrical nonlinearities for solid elements, multi-patch analysis, PUM enrichment for holes and inclusions, and least-squares enforcement of essential boundary conditions.
  • Mass matrices and time-integration schemes support transient analysis on CAD objects, and the open-source code is available for Linux and Mac OS machines.
  • Remaining challenges include constructing volume discretisations from CAD boundary representations, efficient integration schemes, and suitable error estimators.

Appendix A. Knot vector conventions

The appendix defines open knot-vector conventions and their simplified interval representation for B-spline and T-spline algorithms.

  • An open knot vector repeats its first and last knot values p + 1 times.
  • For open knot vectors, redundant endpoint entries can be removed, leaving endpoint values repeated p times in the simplified notation.
  • Knot interval vectors store successive differences as ∆ξ_i = ξ_i+1 − ξ_i.
  • T-splines are based entirely on knot interval vector notation.

Appendix B. Ck and Gk continuity

The appendix distinguishes parameter-based Ck continuity from geometric Gk continuity and describes their implications for curve derivatives and reparameterisation.

  • G1 continuity preserves tangent direction but not necessarily magnitude, and G2 continuity preserves acceleration direction but not necessarily magnitude.
  • Ck continuity requires derivatives through order k to be continuous and single-valued at a parameter value.
  • C1 continuity requires equal tangent magnitude and direction, while C2 additionally requires equal acceleration magnitude.
  • G continuity is independent of the parameter, allowing multiple parameter values while preserving geometric continuity.
  • A Gk curve can be made Ck-continuous through an appropriate reparameterisation.

Appendix C. Homogenous and non-homogenous coordinates

The appendix explains how homogeneous coordinates support rational basis-function algorithms and how to recover physical coordinates from them.

  • Homogeneous coordinates allow B-spline algorithms to be applied directly to rational NURBS basis functions.
  • A physical coordinate P_A = (x_A, y_A, z_A) with weight w_A is represented by a corresponding four-dimensional homogeneous coordinate.
  • Physical coordinates are recovered by dividing the homogeneous coordinate components by w_A.

Appendix D. B´ezier extraction

Bézier extraction expresses B-spline bases through Bernstein polynomials, allowing NURBS and T-splines to be incorporated into existing finite-element codes. Extraction operators are computed during preprocessing so standard FE solvers can use these geometries directly.

  • Bézier extraction facilitates incorporating NURBS and T-splines into finite-element codes.
  • B-spline basis functions can be written as linear combinations of Bernstein polynomials.
  • For element e, the shape functions use an elemental Bézier extraction operator together with Bernstein polynomials on the parent element.
  • Index space and knot vectors are embedded in Bézier extractors computed during preprocessing.
  • Existing finite-element solvers can therefore use NURBS and T-splines straightforwardly.
  • Examples of Bézier extraction are provided in the bezier-extraction folder.
Loading 1205.2129v2…