1 Introduction to the Finite Element Method

Engineering analysis of a continuum often begins with a mathematical model. Temperature, displacement, velocity, concentration, and electromagnetic potential are treated as fields defined over a body or region of space. Balance laws constrain how these fields can vary, constitutive equations describe the response of a particular material, and boundary and initial conditions connect the model to its surroundings. Together, these ingredients produce partial differential equations (PDEs) whose solutions represent the idealized response of the continuum.

Closed-form solutions are uncommon once geometry, material properties, or loading become realistic. A numerical method is therefore needed to replace the continuum problem by a finite calculation. The finite element method (FEM) makes this replacement through a sequence of mathematically distinct steps: it recasts the governing equations in variational form, restricts the admissible fields to a finite-dimensional space constructed over a mesh, assembles the resulting algebraic equations, and solves those equations on a computer.

The variational structure of FEM connects physical principles, approximation theory, and numerical linear algebra. It allows us to ask whether a formulation is well posed, whether a discretization is stable, how the approximation changes under mesh refinement, and what evidence supports a computed result. These questions are part of the analysis, alongside the calculation of the numerical solution.

This chapter introduces these stages without relying on the technical details developed later. It distinguishes the physical model, the PDE, the variational problem, the finite element approximation, and the algebraic system. The following chapters examine each stage in greater detail. The notation used throughout the book is collected in the Notation and Conventions reference and defined in context as it enters the analysis.

1.1 PDE-governed engineering problems

A continuum model begins by choosing the fields and mechanisms needed to answer a particular question about a physical system. Consider heat transfer in a solid body occupying a domain $\Om\subset\R^d$. The primary unknown is the temperature $T(\bx,t)$. Conservation of energy relates its rate of change to heat flux and internal heat generation. Fourier’s law, the constitutive relation in this model, relates the heat flux $\bq$ to the temperature gradient. For an isotropic material, the resulting equation has the form

\[ \rho c\,\dot T-\divg(k\grad T)=r \qquad\text{in }\Om, \tag{1.1}\]

where $\rho$ is density, $c$ is specific heat, $k$ is thermal conductivity, and $r$ is a volumetric heat source. A complete problem must also prescribe suitable temperature or heat-flux conditions on the boundary and an initial temperature.

Now consider a deformable solid. The primary unknown is the displacement field $\bu(\bx)$, which describes how a material point moves from its original position. For small quasistatic deformations, the balance of linear momentum may be written schematically as

\[ -\divg\bsig=\bb \qquad\text{in }\Om, \tag{1.2}\]

where $\bsig$ is the Cauchy stress tensor and $\bb$ is the body force per unit volume. Kinematics relate the displacement gradient to strain. For example, small-deformation theory uses $\beps = (\nabla \bu + \nabla \bu^T)/2$. A constitutive law then relates strain to stress; linear elasticity, for example, uses the linear relation $\bsig = \mathbb{C} \beps$. Displacements may be prescribed on one portion of the boundary and tractions on another.

Equation 1.1 and Equation 1.2 describe different physics and have unknowns of different types: temperature is a scalar field, whereas displacement is a vector-valued field. Nevertheless, the two problems share an important structure:

a field is sought over a domain;
a balance law governs the field in the interior;
a constitutive law models material behavior;
boundary and possibly initial data complete the problem;
quantities such as heat flux, stress, or reaction force are derived from the solution.

This common structure allows one mathematical framework to treat a broad class of continuum problems. It does not make the models interchangeable. Choosing a field, balance law, constitutive response, and boundary data remains a physical modeling decision, and no discretization can correct an inappropriate model.

1.1.1 Models answer specific questions

A physical system admits many mathematical models. A thermal analysis may neglect radiation, treat conductivity as constant, or assume a steady state. A structural analysis may assume infinitesimal strain, quasistatic loading, and a linearly elastic material. Each assumption limits the conditions under which the model is appropriate. For example, a model based on linearized strain and small deformations is generally unsuitable for rubbers and other soft materials undergoing large deformations.

Suppose that a metal support is heated along one surface and fixed mechanically along another. A thermal model may predict temperature and heat flux while saying nothing about thermal stress. A linear elastic model may predict small displacements under moderate loading but become inappropriate after large rotations or plastic deformation. These are limitations of the model, and mesh refinement cannot remove them. If observations include stresses at several points, a mechanical model that omits thermal effects may fail to represent the data even when its finite element discretization is highly accurate.

For this reason, a finite element calculation should begin with an explicit statement of:

the domain and relevant geometric idealizations;
the primary unknown fields;
the balance laws and constitutive assumptions;
boundary and initial conditions;
parameters and source terms;
the output quantities, or quantities of interest (QoIs), needed from the solution.

This statement separates model construction from numerical solution and gives a basis for interpreting the result later.

1.2 From a physical model to a numerical solution

The book organizes each model problem through the following sequence:

\[ \begin{array}{ccccccc} \boxed{\substack{\text{physical}\\\text{model}}} &\longrightarrow& \boxed{\substack{\text{balance laws}\\\text{and PDEs}}} &\longrightarrow& \boxed{\substack{\text{boundary and}\\\text{initial data}}} &\longrightarrow& \boxed{\substack{\text{variational}\\\text{problem}}} \\ &&&&&&\Big\downarrow \\ \boxed{\substack{\text{verification and}\\\text{interpretation}}} &\longleftarrow& \boxed{\substack{\text{numerical}\\\text{solution}}} &\longleftarrow& \boxed{\substack{\text{algebraic}\\\text{system}}} &\longleftarrow& \boxed{\substack{\text{finite element}\\\text{problem}}}. \end{array} \tag{1.3}\]

Each box in Equation 1.3 denotes a distinct stage of modeling or computation.

1.2.1 Physical model

The physical model identifies the idealized system, its state variables, and the mechanisms retained in the analysis. Material response, source terms, and loading assumptions enter at this level. At this stage, the analyst specifies what is being modeled and the assumptions under which the model will be used. The governing differential equations follow from those choices.

1.2.2 Balance laws, constitutive equations, and PDEs

Balance laws express conservation of quantities such as energy, mass, or linear momentum. Constitutive equations close the system by relating fluxes or stresses to the primary fields. Combining these statements produces governing PDEs such as Equation 1.1 and Equation 1.2. The resulting equations are usually called the strong form when they are required to hold pointwise and all derivatives are interpreted classically.

1.2.3 Boundary and initial data

A differential equation alone is not enough to obtain a unique solution. Boundary conditions describe interaction with the environment: a temperature may be prescribed, a heat flux imposed, a displacement constrained, or a traction applied. Time-dependent problems also require initial data. These conditions are part of the mathematical model, not implementation details to be added after discretization.

1.2.4 Variational problem

The variational problem asks that the governing equations hold through their action on an admissible class of test functions. Derivatives are redistributed using integration by parts, and boundary terms become part of the formulation. The unknown is still a function drawn from an infinite-dimensional space. No finite-dimensional approximation has been introduced at this stage. Under appropriate assumptions, the variational problem is equivalent to the strong form, but it admits solutions with weaker differentiability. Such a solution is called a weak solution to distinguish it from a classical, or strong, solution of the strong form.

1.2.5 Finite element problem

The finite element problem restricts the trial and test functions to finite-dimensional spaces constructed over a mesh. This is the approximation step. Mesh geometry, polynomial order, interelement continuity, and the treatment of essential boundary data determine the discrete spaces.

1.2.6 Algebraic system

After finite element discretization, the variational problem becomes a system of algebraic equations for the coefficients associated with the finite element degrees of freedom. A linear scalar PDE may produce

\[ \mathbf K\mathbf U=\mathbf F, \tag{1.4}\]

where $\mathbf U$ contains coefficients of the finite element field. Nonlinear, mixed, constrained, or transient formulations lead to different algebraic structures. The coefficient vector is a representation of the discrete field; it is not the continuum solution itself.

1.2.7 Numerical solution

An algebraic solver computes an approximation to the coefficient vector. Direct and iterative solvers, preconditioners, stopping tolerances, and floating-point arithmetic enter at this stage. Solving the algebraic equations inaccurately can obscure the behavior of an otherwise stable discretization.

1.2.8 Verification and interpretation

The computed field must be checked and interpreted. Verification asks whether the equations were discretized and solved correctly. Interpretation asks whether derived quantities are meaningful for the original question. Keeping these levels separate prevents common reasoning errors. A solver failure is not automatically a defect in the PDE, a mesh-refinement trend does not validate a constitutive model, and agreement with the assembled equations does not prove agreement with the continuum problem.

1.3 Why the variational formulation matters

The finite element method is naturally built from a variational formulation rather than directly from pointwise differential equations. To see why, consider the representative scalar problem

\[ -\divg(k\grad u)=f \qquad\text{in }\Om. \tag{1.5}\]

The differential operator contains second derivatives of $u$ when $k$ is smooth. A classical solution must possess enough pointwise differentiability for those derivatives to exist and satisfy the equation. Such regularity can fail near corners, across material interfaces, or where boundary conditions change.

In a variational formulation, the equation is multiplied by a test function and integrated over the domain. Integration by parts transfers one derivative from the unknown to the test function and exposes a boundary term. The resulting statement has the abstract form

\[ \text{find }u\in\Ug\text{ such that} \qquad a(u,v)=\ell(v) \quad\text{for every }v\in\Vzero. \tag{1.6}\]

Here $\Ug$ is a trial space containing fields that satisfy prescribed values on the Dirichlet boundary, and $\Vzero$ is the associated test space with homogeneous boundary values. The form $a(u,v)$ contains the differential operator after integration by parts, while $\ell(v)$ contains sources and prescribed fluxes or tractions. Chapter 3 constructs this statement carefully.

Several important consequences are already visible.

1.3.1 Weaker differentiability requirements

The variational problem can admit fields whose derivatives exist only in the weak sense. Weak derivatives and Sobolev spaces provide the precise language for this enlargement of the admissible class. Chapter 2 develops that language. This broader solution space is important because piecewise polynomial finite element functions are generally not globally smooth enough to satisfy a second-order PDE pointwise, but they can satisfy the corresponding variational equation.

1.3.2 Boundary conditions enter in different ways

Prescribed values are built into the admissible trial space and are called essential boundary conditions. Fluxes and tractions naturally arise in boundary integrals and are called natural boundary conditions. This distinction follows from the variational structure.

1.3.3 Physical principles appear naturally

In mechanics, the test function can be interpreted as a virtual displacement, and the variational equation expresses virtual work. In symmetric conservative problems, the same solution may make an energy functional stationary or minimize it over an admissible set. A variational formulation and an energy principle are closely related in those cases, but they are not synonymous: many useful variational problems do not arise from minimization.

1.3.4 Properties of the variational problem

The form $a(u,v)=\ell(v)$ can be examined for properties such as symmetry, continuity, and coercivity. For a conforming Galerkin method, appropriate continuous properties are inherited by the finite element problem. These properties govern existence, uniqueness, and stability and later help determine the quality of the finite element approximation.

1.4 Approximation and discretization

The variational problem remains infinite dimensional: its unknown is a function, and the equation must hold for every admissible test function. Computation requires finite-dimensional trial and test spaces. FEM creates those spaces from a mesh and local approximation rules.

1.4.1 Meshes and local approximation

A mesh $\Th$ partitions the domain into nonoverlapping elements $K$. Typical elements include intervals in one dimension, triangles and quadrilaterals in two dimensions, and tetrahedra and hexahedra in three dimensions. On each element, the unknown is represented by a simple local function, usually a polynomial. Neighboring representations are connected so that the assembled field has the interelement continuity required by the variational problem. Discontinuous Galerkin methods relax this continuity requirement and instead couple neighboring elements through interface terms (Arnold et al. 2002). This book focuses on conforming finite element methods; discontinuous Galerkin methods are outside its present scope.

Each geometric member of the mesh is called a cell or element. A finite element is more precisely the combination of an element, a local function space, and a set of degrees of freedom.

1.4.2 Finite-dimensional spaces

Let $\Vh$ denote a finite-dimensional test space constructed over $\Th$. A conforming method chooses $\Vh\subset\Vzero$, so every discrete test function is also admissible for the continuous variational problem. The corresponding Galerkin statement has the schematic form

\[ \text{find }\uh\in\Uh\text{ such that} \qquad a(\uh,\vh)=\ell(\vh) \quad\text{for every }\vh\in\Vh. \tag{1.7}\]

The subscript $h$ indicates dependence on the mesh. It does not mean that the discrete field is simply the exact field evaluated at grid points. The finite element solution is determined by the discrete variational equation over an entire function space.

1.4.3 Local-to-global construction

Finite element basis functions have local support: each is nonzero only on a small neighborhood of elements. Integrals can therefore be evaluated element by element and accumulated into global matrices and vectors. Two basis functions that do not overlap make no joint contribution, which produces sparse algebraic systems.

This local-to-global construction is one of the defining computational features of FEM. It permits complex domains to be assembled from standard element types, allows material data to vary by region, and makes refinement possible in selected parts of the mesh. It also permits distributed computing, in which different processors assemble contributions from the portions of the mesh they own.

1.4.4 Reference elements and geometric flexibility

Calculations on many physical elements are organized by mapping them from a small collection of reference elements. Shape functions, quadrature rules, and derivatives can be defined once on a reference element and transformed to each physical element. When the geometry and the unknown field are represented using the same type of basis, the construction is called isoparametric.

These ideas explain how FEM accommodates curved boundaries and complicated geometries while retaining a systematic implementation. Chapter 5 develops the finite-element triple, Lagrange spaces, element mappings, quadrature, assembly, and discrete boundary conditions in detail.

1.5 Assessing the reliability of a computed result

A completed solve establishes that an algebraic procedure returned a result. Further checks are needed to determine whether that result approximates the intended physical system with acceptable accuracy. The following sources of discrepancy should be considered separately.

Source	Question
Modeling error	Do the physical assumptions represent the system of interest?
Data uncertainty	Are geometry, material parameters, loads, and boundary data known adequately?
Discretization error	Does the finite-dimensional space approximate the continuum solution sufficiently well?
Quadrature error	Were element and boundary integrals evaluated accurately enough?
Algebraic error	Were the discrete equations solved to an adequate tolerance?
Implementation error	Does the code represent the intended formulation?
Interpretation error	Are reported outputs consistent with the model and discretization?

The distinction between verification and validation is useful. Verification asks whether the mathematical model has been solved correctly. Validation asks whether that model is an adequate representation of the physical system for the intended use. Mesh refinement, exact solutions, manufactured solutions, conservation checks, and independent implementations support verification. Comparison with experimental or trusted observational data supports validation.

1.5.2 Derived quantities require separate attention

The primary unknown may be temperature or displacement, while the engineering quantity of interest is heat flux, stress, reaction force, average temperature, or compliance. Such outputs involve derivatives, boundary reactions, or integrals and need not inherit the same accuracy as the primary field. A credible analysis identifies the quantity of interest before choosing the mesh and verification strategy.

1.6 Computational framework used in this book

The computational material follows the same order as the mathematics: model, variational form, discrete spaces, algebraic problem, solver, and verification.

1.6.1 Python and Jupyter

Python provides the common language for numerical examples, preprocessing, postprocessing, and verification. Jupyter allows executable computations, figures, and explanatory text to remain together. Short examples appear directly in the chapters; longer investigations belong in notebooks, and reusable implementations belong in source modules or examples.

1.6.2 FEniCSx and UFL

FEniCSx represents meshes, finite element spaces, functions, boundary conditions, and variational forms as computational objects. The Unified Form Language (UFL) expresses integrals using notation close to their mathematical form. From those forms, the software generates element-level computations and assembles global algebraic objects.

1.6.3 PETSc and scalable algebra

The assembled matrices and vectors are managed through PETSc. PETSc supplies linear and nonlinear solvers, preconditioners, distributed data structures, and parallel algorithms.

1.7 Scope and organization of the book

The book begins with a focused review and then develops the analytical and computational tools required for advanced formulations.

1.7.1 Review and foundations

Chapter 2 introduces vector and tensor operations, function spaces, weak derivatives, Sobolev spaces, norms, bilinear forms, and well-posedness. Chapter 3 derives variational formulations from boundary-value problems and explains trial and test spaces, essential and natural boundary conditions, and energy stationarity. Chapter 4 develops heat conduction, linear elasticity, and the finite-deformation quantities needed for hyperelasticity. Chapter 5 constructs finite element spaces and connects meshes, basis functions, mappings, quadrature, assembly, boundary conditions, and FEniCSx.

1.7.2 Error analysis and adaptivity

A priori analysis relates stability and approximation properties to convergence. A posteriori analysis seeks computable information about the error in a completed solution. Adaptive FEM uses that information to refine or enrich the discretization where additional resolution is most valuable.

1.7.3 Constraints and mixed formulations

Constraints arise in contact, incompressibility, and other admissibility conditions. Mixed methods introduce multiple primary fields, such as displacement and pressure or flux and temperature. Their stability depends on compatibility between spaces and cannot be inferred from the behavior of a single scalar coercive problem.

1.7.4 Nonlinear and transient problems

Finite-deformation mechanics produces nonlinear residual equations. Their solution requires consistent linearization, tangent operators, and iterative methods. Transient problems add time discretization and require separate control of spatial and temporal errors.

1.7.5 PDE-constrained optimization

Optimization treats a PDE solution as a constraint while adjusting design or control variables to improve an objective. Adjoint equations make sensitivities computationally tractable and bring together variational formulations, linearized operators, solver design, and quantities of interest.

Across all these topics, the same questions recur:

What is the physical and mathematical problem?
Which trial and test spaces express its structure?
Is the continuous problem well posed?
Is the discrete formulation stable and consistent?
What controls the approximation error?
How should the algebraic equations be solved?
What evidence supports the reported result?

The chapters return to these questions when introducing each formulation and its computational implementation.

1.8 Key ideas

A continuum model, its strong PDE, its variational formulation, its finite element approximation, and its algebraic system are distinct mathematical objects.
Boundary and initial conditions are part of the model, not postprocessing or software details.
Variational formulations reduce differentiability requirements, expose natural boundary terms, and provide the structure used by FEM.
Finite element spaces combine local approximation on mesh elements with a global continuity structure.
Local support leads to elementwise computation and sparse global systems.
A computed solution is credible only when modeling, discretization, algebra, implementation, and interpretation errors are considered separately.
Computation should follow the mathematical formulation and include explicit verification evidence.

1.9 Exercises

A rectangular plate $\Om=(0,L)\times(0,H)$ has constant thermal conductivity $k$ and a volumetric heat source $r$. The temperature is prescribed as $T_0$ on the left edge, the normal heat flux is prescribed as $\bar q$ on the right edge, and the top and bottom edges are insulated. Without solving the problem, state:
1. the primary unknown and its type (scalar, vector, or tensor);
2. Fourier’s constitutive law for $\bq$;
3. the steady governing PDE;
4. the boundary condition on each edge, classified as essential or natural;
5. three distinct QoIs that could be computed from the solution.
A two-dimensional cantilever of length $L$ and height $H$ is clamped at $x=0$ and carries a prescribed downward traction at $x=L$. It is modeled using small-strain linear elasticity. Prepare a table containing two modeling assumptions, two discretization choices, and two verification checks. For each entry, state whether mesh refinement can address it and justify the answer in one sentence.
Classify each item below as belonging primarily to the physical model, variational problem, finite element approximation, algebraic system, or numerical solver. Give a one-sentence justification for each classification:
1. a prescribed traction;
2. a test function;
3. a quadratic triangular element;
4. a Krylov-method stopping tolerance;
5. a coefficient vector $\mathbf U$;
6. a constitutive tensor $\mathbb C$.
Let $u_h$ be continuous and piecewise linear on a conforming triangular mesh. State where its classical second derivatives exist and explain why $u_h$ generally cannot satisfy Equation 1.5 pointwise across element boundaries. Then explain, using the integration-by-parts argument from this chapter, why $u_h$ can still be admissible in the corresponding variational problem.
A finite element analysis of the cantilever in Exercise 2 produces a smooth displacement plot, but the reported QoI is the average normal stress over a specified region near the clamp. Design a mesh-refinement study for this QoI. State the sequence of meshes to be compared, the value to record on each mesh, and a numerical criterion for deciding that the QoI is sufficiently resolved. Explain why visual smoothness of the displacement field is not that criterion.
A tensile coupon is tested in a laboratory, and the same test is represented by a structural finite element model. Give two verification activities and two validation activities for this study. For each activity, identify the reference result or data needed and the discrepancy that would be measured.
For each item below, identify the box in Equation 1.3 where it first enters and state one error or uncertainty it may introduce:
1. calibration of a material parameter from test data;
2. selection of a numerical quadrature rule;
3. local mesh refinement;
4. selection of a linear solver and stopping tolerance;
5. comparison of a computed QoI with experimental measurements.

1.10 Further reading

Broad treatments of the mathematical and variational foundations of FEM are available in (Oden and Reddy 1976; Süli 2012). Mechanics-oriented developments and computational constructions can be found in (Jog 1978; Becker et al. 1981), while (Larson and Bengzon 2010; Arnold 2011) provide further perspectives on approximation, implementation, stability, and error assessment. Oden and Demkowicz provide the functional-analysis background used in this book (Oden and Demkowicz 2018). Discontinuous Galerkin methods are outside the present scope; Arnold et al. provide a unified analysis of these methods for elliptic problems (Arnold et al. 2002).

# Introduction to the Finite Element Method {#sec-introduction} ```{=latex} \chaptershorttitle{Introduction} ``` Engineering analysis of a continuum often begins with a mathematical model. Temperature, displacement, velocity, concentration, and electromagnetic potential are treated as fields defined over a body or region of space. Balance laws constrain how these fields can vary, constitutive equations describe the response of a particular material, and boundary and initial conditions connect the model to its surroundings. Together, these ingredients produce partial differential equations (PDEs) whose solutions represent the idealized response of the continuum. Closed-form solutions are uncommon once geometry, material properties, or loading become realistic. A numerical method is therefore needed to replace the continuum problem by a finite calculation. The finite element method (FEM) makes this replacement through a sequence of mathematically distinct steps: it recasts the governing equations in variational form, restricts the admissible fields to a finite-dimensional space constructed over a mesh, assembles the resulting algebraic equations, and solves those equations on a computer. The variational structure of FEM connects physical principles, approximation theory, and numerical linear algebra. It allows us to ask whether a formulation is well posed, whether a discretization is stable, how the approximation changes under mesh refinement, and what evidence supports a computed result. These questions are part of the analysis, alongside the calculation of the numerical solution. This chapter introduces these stages without relying on the technical details developed later. It distinguishes the physical model, the PDE, the variational problem, the finite element approximation, and the algebraic system. The following chapters examine each stage in greater detail. The notation used throughout the book is collected in the [Notation and Conventions](../notation.qmd) reference and defined in context as it enters the analysis. ## PDE-governed engineering problems {#sec-intro-pde-problems} A continuum model begins by choosing the fields and mechanisms needed to answer a particular question about a physical system. Consider heat transfer in a solid body occupying a domain $\Om\subset\R^d$. The primary unknown is the temperature $T(\bx,t)$. Conservation of energy relates its rate of change to heat flux and internal heat generation. Fourier's law, the constitutive relation in this model, relates the heat flux $\bq$ to the temperature gradient. For an isotropic material, the resulting equation has the form $$ \rho c\,\dot T-\divg(k\grad T)=r \qquad\text{in }\Om, $$ {#eq-intro-heat} where $\rho$ is density, $c$ is specific heat, $k$ is thermal conductivity, and $r$ is a volumetric heat source. A complete problem must also prescribe suitable temperature or heat-flux conditions on the boundary and an initial temperature. Now consider a deformable solid. The primary unknown is the displacement field $\bu(\bx)$, which describes how a material point moves from its original position. For small quasistatic deformations, the balance of linear momentum may be written schematically as $$ -\divg\bsig=\bb \qquad\text{in }\Om, $$ {#eq-intro-elasticity} where $\bsig$ is the Cauchy stress tensor and $\bb$ is the body force per unit volume. Kinematics relate the displacement gradient to strain. For example, small-deformation theory uses $\beps = (\nabla \bu + \nabla \bu^T)/2$. A constitutive law then relates strain to stress; linear elasticity, for example, uses the linear relation $\bsig = \mathbb{C} \beps$. Displacements may be prescribed on one portion of the boundary and tractions on another. @eq-intro-heat and @eq-intro-elasticity describe different physics and have unknowns of different types: temperature is a scalar field, whereas displacement is a vector-valued field. Nevertheless, the two problems share an important structure: - a field is sought over a domain; - a balance law governs the field in the interior; - a constitutive law models material behavior; - boundary and possibly initial data complete the problem; - quantities such as heat flux, stress, or reaction force are derived from the solution. This common structure allows one mathematical framework to treat a broad class of continuum problems. It does not make the models interchangeable. Choosing a field, balance law, constitutive response, and boundary data remains a physical modeling decision, and no discretization can correct an inappropriate model. ### Models answer specific questions A physical system admits many mathematical models. A thermal analysis may neglect radiation, treat conductivity as constant, or assume a steady state. A structural analysis may assume infinitesimal strain, quasistatic loading, and a linearly elastic material. Each assumption limits the conditions under which the model is appropriate. For example, a model based on linearized strain and small deformations is generally unsuitable for rubbers and other soft materials undergoing large deformations. Suppose that a metal support is heated along one surface and fixed mechanically along another. A thermal model may predict temperature and heat flux while saying nothing about thermal stress. A linear elastic model may predict small displacements under moderate loading but become inappropriate after large rotations or plastic deformation. These are limitations of the model, and mesh refinement cannot remove them. If observations include stresses at several points, a mechanical model that omits thermal effects may fail to represent the data even when its finite element discretization is highly accurate. For this reason, a finite element calculation should begin with an explicit statement of: 1. the domain and relevant geometric idealizations; 2. the primary unknown fields; 3. the balance laws and constitutive assumptions; 4. boundary and initial conditions; 5. parameters and source terms; 6. the output quantities, or quantities of interest (QoIs), needed from the solution. This statement separates model construction from numerical solution and gives a basis for interpreting the result later. ## From a physical model to a numerical solution {#sec-intro-model-to-solution} The book organizes each model problem through the following sequence: $$ \begin{array}{ccccccc} \boxed{\substack{\text{physical}\\\text{model}}} &\longrightarrow& \boxed{\substack{\text{balance laws}\\\text{and PDEs}}} &\longrightarrow& \boxed{\substack{\text{boundary and}\\\text{initial data}}} &\longrightarrow& \boxed{\substack{\text{variational}\\\text{problem}}} \\ &&&&&&\Big\downarrow \\ \boxed{\substack{\text{verification and}\\\text{interpretation}}} &\longleftarrow& \boxed{\substack{\text{numerical}\\\text{solution}}} &\longleftarrow& \boxed{\substack{\text{algebraic}\\\text{system}}} &\longleftarrow& \boxed{\substack{\text{finite element}\\\text{problem}}}. \end{array} $$ {#eq-intro-workflow} Each box in @eq-intro-workflow denotes a distinct stage of modeling or computation. ### Physical model The physical model identifies the idealized system, its state variables, and the mechanisms retained in the analysis. Material response, source terms, and loading assumptions enter at this level. At this stage, the analyst specifies what is being modeled and the assumptions under which the model will be used. The governing differential equations follow from those choices. ### Balance laws, constitutive equations, and PDEs Balance laws express conservation of quantities such as energy, mass, or linear momentum. Constitutive equations close the system by relating fluxes or stresses to the primary fields. Combining these statements produces governing PDEs such as @eq-intro-heat and @eq-intro-elasticity. The resulting equations are usually called the *strong form* when they are required to hold pointwise and all derivatives are interpreted classically. ### Boundary and initial data A differential equation alone is not enough to obtain a unique solution. Boundary conditions describe interaction with the environment: a temperature may be prescribed, a heat flux imposed, a displacement constrained, or a traction applied. Time-dependent problems also require initial data. These conditions are part of the mathematical model, not implementation details to be added after discretization. ### Variational problem The variational problem asks that the governing equations hold through their action on an admissible class of test functions. Derivatives are redistributed using integration by parts, and boundary terms become part of the formulation. The unknown is still a function drawn from an infinite-dimensional space. No finite-dimensional approximation has been introduced at this stage. Under appropriate assumptions, the variational problem is equivalent to the strong form, but it admits solutions with weaker differentiability. Such a solution is called a *weak solution* to distinguish it from a classical, or strong, solution of the strong form. ### Finite element problem The finite element problem restricts the trial and test functions to finite-dimensional spaces constructed over a mesh. This is the approximation step. Mesh geometry, polynomial order, interelement continuity, and the treatment of essential boundary data determine the discrete spaces. ### Algebraic system After finite element discretization, the variational problem becomes a system of algebraic equations for the coefficients associated with the finite element degrees of freedom. A linear scalar PDE may produce $$ \mathbf K\mathbf U=\mathbf F, $$ {#eq-intro-algebraic} where $\mathbf U$ contains coefficients of the finite element field. Nonlinear, mixed, constrained, or transient formulations lead to different algebraic structures. The coefficient vector is a representation of the discrete field; it is not the continuum solution itself. ### Numerical solution An algebraic solver computes an approximation to the coefficient vector. Direct and iterative solvers, preconditioners, stopping tolerances, and floating-point arithmetic enter at this stage. Solving the algebraic equations inaccurately can obscure the behavior of an otherwise stable discretization. ### Verification and interpretation The computed field must be checked and interpreted. Verification asks whether the equations were discretized and solved correctly. Interpretation asks whether derived quantities are meaningful for the original question. Keeping these levels separate prevents common reasoning errors. A solver failure is not automatically a defect in the PDE, a mesh-refinement trend does not validate a constitutive model, and agreement with the assembled equations does not prove agreement with the continuum problem. ## Why the variational formulation matters {#sec-intro-why-variational} The finite element method is naturally built from a variational formulation rather than directly from pointwise differential equations. To see why, consider the representative scalar problem $$ -\divg(k\grad u)=f \qquad\text{in }\Om. $$ {#eq-intro-diffusion-preview} The differential operator contains second derivatives of $u$ when $k$ is smooth. A classical solution must possess enough pointwise differentiability for those derivatives to exist and satisfy the equation. Such regularity can fail near corners, across material interfaces, or where boundary conditions change. In a variational formulation, the equation is multiplied by a test function and integrated over the domain. Integration by parts transfers one derivative from the unknown to the test function and exposes a boundary term. The resulting statement has the abstract form $$ \text{find }u\in\Ug\text{ such that} \qquad a(u,v)=\ell(v) \quad\text{for every }v\in\Vzero. $$ {#eq-intro-variational-preview} Here $\Ug$ is a trial space containing fields that satisfy prescribed values on the Dirichlet boundary, and $\Vzero$ is the associated test space with homogeneous boundary values. The form $a(u,v)$ contains the differential operator after integration by parts, while $\ell(v)$ contains sources and prescribed fluxes or tractions. Chapter 3 constructs this statement carefully. Several important consequences are already visible. ### Weaker differentiability requirements The variational problem can admit fields whose derivatives exist only in the weak sense. Weak derivatives and Sobolev spaces provide the precise language for this enlargement of the admissible class. Chapter 2 develops that language. This broader solution space is important because piecewise polynomial finite element functions are generally not globally smooth enough to satisfy a second-order PDE pointwise, but they can satisfy the corresponding variational equation. ### Boundary conditions enter in different ways Prescribed values are built into the admissible trial space and are called *essential* boundary conditions. Fluxes and tractions naturally arise in boundary integrals and are called *natural* boundary conditions. This distinction follows from the variational structure. ### Physical principles appear naturally In mechanics, the test function can be interpreted as a virtual displacement, and the variational equation expresses virtual work. In symmetric conservative problems, the same solution may make an energy functional stationary or minimize it over an admissible set. A variational formulation and an energy principle are closely related in those cases, but they are not synonymous: many useful variational problems do not arise from minimization. ### Properties of the variational problem The form $a(u,v)=\ell(v)$ can be examined for properties such as symmetry, continuity, and coercivity. For a conforming Galerkin method, appropriate continuous properties are inherited by the finite element problem. These properties govern existence, uniqueness, and stability and later help determine the quality of the finite element approximation. ## Approximation and discretization {#sec-intro-approximation} The variational problem remains infinite dimensional: its unknown is a function, and the equation must hold for every admissible test function. Computation requires finite-dimensional trial and test spaces. FEM creates those spaces from a mesh and local approximation rules. ### Meshes and local approximation A mesh $\Th$ partitions the domain into nonoverlapping elements $K$. Typical elements include intervals in one dimension, triangles and quadrilaterals in two dimensions, and tetrahedra and hexahedra in three dimensions. On each element, the unknown is represented by a simple local function, usually a polynomial. Neighboring representations are connected so that the assembled field has the interelement continuity required by the variational problem. Discontinuous Galerkin methods relax this continuity requirement and instead couple neighboring elements through interface terms [@arnold-brezzi-cockburn-marini2002]. This book focuses on conforming finite element methods; discontinuous Galerkin methods are outside its present scope. Each geometric member of the mesh is called a *cell* or *element*. A finite element is more precisely the combination of an element, a local function space, and a set of degrees of freedom. ### Finite-dimensional spaces Let $\Vh$ denote a finite-dimensional test space constructed over $\Th$. A conforming method chooses $\Vh\subset\Vzero$, so every discrete test function is also admissible for the continuous variational problem. The corresponding Galerkin statement has the schematic form $$ \text{find }\uh\in\Uh\text{ such that} \qquad a(\uh,\vh)=\ell(\vh) \quad\text{for every }\vh\in\Vh. $$ {#eq-intro-discrete-preview} The subscript $h$ indicates dependence on the mesh. It does not mean that the discrete field is simply the exact field evaluated at grid points. The finite element solution is determined by the discrete variational equation over an entire function space. ### Local-to-global construction Finite element basis functions have local support: each is nonzero only on a small neighborhood of elements. Integrals can therefore be evaluated element by element and accumulated into global matrices and vectors. Two basis functions that do not overlap make no joint contribution, which produces sparse algebraic systems. This local-to-global construction is one of the defining computational features of FEM. It permits complex domains to be assembled from standard element types, allows material data to vary by region, and makes refinement possible in selected parts of the mesh. It also permits distributed computing, in which different processors assemble contributions from the portions of the mesh they own. ### Reference elements and geometric flexibility Calculations on many physical elements are organized by mapping them from a small collection of reference elements. Shape functions, quadrature rules, and derivatives can be defined once on a reference element and transformed to each physical element. When the geometry and the unknown field are represented using the same type of basis, the construction is called *isoparametric*. These ideas explain how FEM accommodates curved boundaries and complicated geometries while retaining a systematic implementation. Chapter 5 develops the finite-element triple, Lagrange spaces, element mappings, quadrature, assembly, and discrete boundary conditions in detail. ## Assessing the reliability of a computed result {#sec-intro-reliability} A completed solve establishes that an algebraic procedure returned a result. Further checks are needed to determine whether that result approximates the intended physical system with acceptable accuracy. The following sources of discrepancy should be considered separately. | Source | Question | |---|---| | Modeling error | Do the physical assumptions represent the system of interest? | | Data uncertainty | Are geometry, material parameters, loads, and boundary data known adequately? | | Discretization error | Does the finite-dimensional space approximate the continuum solution sufficiently well? | | Quadrature error | Were element and boundary integrals evaluated accurately enough? | | Algebraic error | Were the discrete equations solved to an adequate tolerance? | | Implementation error | Does the code represent the intended formulation? | | Interpretation error | Are reported outputs consistent with the model and discretization? | The distinction between *verification* and *validation* is useful. Verification asks whether the mathematical model has been solved correctly. Validation asks whether that model is an adequate representation of the physical system for the intended use. Mesh refinement, exact solutions, manufactured solutions, conservation checks, and independent implementations support verification. Comparison with experimental or trusted observational data supports validation. ### Interpreting refinement studies Under the assumptions required for convergence, refining the mesh or increasing the polynomial order should bring a consistent and stable discretization closer to the continuum solution. A refinement study provides useful numerical evidence, but its results require careful interpretation. A poor model can converge, an incorrectly implemented method can show an accidental plateau, and some quantities converge more slowly than the primary field. Later chapters make refinement studies quantitative. A priori analysis relates the approximation error to mesh size, polynomial order, solution regularity, and stability. A posteriori analysis uses the computed solution and problem data to estimate where the discretization error is concentrated. Adaptive algorithms then use this information to distribute computational effort. ### Derived quantities require separate attention The primary unknown may be temperature or displacement, while the engineering quantity of interest is heat flux, stress, reaction force, average temperature, or compliance. Such outputs involve derivatives, boundary reactions, or integrals and need not inherit the same accuracy as the primary field. A credible analysis identifies the quantity of interest before choosing the mesh and verification strategy. ## Computational framework used in this book {#sec-intro-computational-framework} The computational material follows the same order as the mathematics: model, variational form, discrete spaces, algebraic problem, solver, and verification. ### Python and Jupyter Python provides the common language for numerical examples, preprocessing, postprocessing, and verification. Jupyter allows executable computations, figures, and explanatory text to remain together. Short examples appear directly in the chapters; longer investigations belong in notebooks, and reusable implementations belong in source modules or examples. ### FEniCSx and UFL FEniCSx represents meshes, finite element spaces, functions, boundary conditions, and variational forms as computational objects. The Unified Form Language (UFL) expresses integrals using notation close to their mathematical form. From those forms, the software generates element-level computations and assembles global algebraic objects. ### PETSc and scalable algebra The assembled matrices and vectors are managed through PETSc. PETSc supplies linear and nonlinear solvers, preconditioners, distributed data structures, and parallel algorithms. ## Scope and organization of the book {#sec-intro-organization} The book begins with a focused review and then develops the analytical and computational tools required for advanced formulations. ### Review and foundations Chapter 2 introduces vector and tensor operations, function spaces, weak derivatives, Sobolev spaces, norms, bilinear forms, and well-posedness. Chapter 3 derives variational formulations from boundary-value problems and explains trial and test spaces, essential and natural boundary conditions, and energy stationarity. Chapter 4 develops heat conduction, linear elasticity, and the finite-deformation quantities needed for hyperelasticity. Chapter 5 constructs finite element spaces and connects meshes, basis functions, mappings, quadrature, assembly, boundary conditions, and FEniCSx. ### Error analysis and adaptivity A priori analysis relates stability and approximation properties to convergence. A posteriori analysis seeks computable information about the error in a completed solution. Adaptive FEM uses that information to refine or enrich the discretization where additional resolution is most valuable. ### Constraints and mixed formulations Constraints arise in contact, incompressibility, and other admissibility conditions. Mixed methods introduce multiple primary fields, such as displacement and pressure or flux and temperature. Their stability depends on compatibility between spaces and cannot be inferred from the behavior of a single scalar coercive problem. ### Nonlinear and transient problems Finite-deformation mechanics produces nonlinear residual equations. Their solution requires consistent linearization, tangent operators, and iterative methods. Transient problems add time discretization and require separate control of spatial and temporal errors. ### PDE-constrained optimization Optimization treats a PDE solution as a constraint while adjusting design or control variables to improve an objective. Adjoint equations make sensitivities computationally tractable and bring together variational formulations, linearized operators, solver design, and quantities of interest. Across all these topics, the same questions recur: 1. What is the physical and mathematical problem? 2. Which trial and test spaces express its structure? 3. Is the continuous problem well posed? 4. Is the discrete formulation stable and consistent? 5. What controls the approximation error? 6. How should the algebraic equations be solved? 7. What evidence supports the reported result? The chapters return to these questions when introducing each formulation and its computational implementation. ## Key ideas {#sec-intro-key-ideas} - A continuum model, its strong PDE, its variational formulation, its finite element approximation, and its algebraic system are distinct mathematical objects. - Boundary and initial conditions are part of the model, not postprocessing or software details. - Variational formulations reduce differentiability requirements, expose natural boundary terms, and provide the structure used by FEM. - Finite element spaces combine local approximation on mesh elements with a global continuity structure. - Local support leads to elementwise computation and sparse global systems. - A computed solution is credible only when modeling, discretization, algebra, implementation, and interpretation errors are considered separately. - Computation should follow the mathematical formulation and include explicit verification evidence. ## Exercises {#sec-intro-exercises} 1. A rectangular plate $\Om=(0,L)\times(0,H)$ has constant thermal conductivity $k$ and a volumetric heat source $r$. The temperature is prescribed as $T_0$ on the left edge, the normal heat flux is prescribed as $\bar q$ on the right edge, and the top and bottom edges are insulated. Without solving the problem, state: a. the primary unknown and its type (scalar, vector, or tensor); b. Fourier's constitutive law for $\bq$; c. the steady governing PDE; d. the boundary condition on each edge, classified as essential or natural; e. three distinct QoIs that could be computed from the solution. 2. A two-dimensional cantilever of length $L$ and height $H$ is clamped at $x=0$ and carries a prescribed downward traction at $x=L$. It is modeled using small-strain linear elasticity. Prepare a table containing two modeling assumptions, two discretization choices, and two verification checks. For each entry, state whether mesh refinement can address it and justify the answer in one sentence. 3. Classify each item below as belonging primarily to the physical model, variational problem, finite element approximation, algebraic system, or numerical solver. Give a one-sentence justification for each classification: a. a prescribed traction; b. a test function; c. a quadratic triangular element; d. a Krylov-method stopping tolerance; e. a coefficient vector $\mathbf U$; f. a constitutive tensor $\mathbb C$. 4. Let $u_h$ be continuous and piecewise linear on a conforming triangular mesh. State where its classical second derivatives exist and explain why $u_h$ generally cannot satisfy @eq-intro-diffusion-preview pointwise across element boundaries. Then explain, using the integration-by-parts argument from this chapter, why $u_h$ can still be admissible in the corresponding variational problem. 5. A finite element analysis of the cantilever in Exercise 2 produces a smooth displacement plot, but the reported QoI is the average normal stress over a specified region near the clamp. Design a mesh-refinement study for this QoI. State the sequence of meshes to be compared, the value to record on each mesh, and a numerical criterion for deciding that the QoI is sufficiently resolved. Explain why visual smoothness of the displacement field is not that criterion. 6. A tensile coupon is tested in a laboratory, and the same test is represented by a structural finite element model. Give two verification activities and two validation activities for this study. For each activity, identify the reference result or data needed and the discrepancy that would be measured. 7. For each item below, identify the box in @eq-intro-workflow where it first enters and state one error or uncertainty it may introduce: a. calibration of a material parameter from test data; b. selection of a numerical quadrature rule; c. local mesh refinement; d. selection of a linear solver and stopping tolerance; e. comparison of a computed QoI with experimental measurements. ## Further reading {#sec-intro-further-reading} Broad treatments of the mathematical and variational foundations of FEM are available in [@oden-reddy-fem; @suli2012]. Mechanics-oriented developments and computational constructions can be found in [@jog; @becker-carey-oden1981], while [@larson-bengzon2010; @arnold2018] provide further perspectives on approximation, implementation, stability, and error assessment. Oden and Demkowicz provide the functional-analysis background used in this book [@oden-demkowicz]. Discontinuous Galerkin methods are outside the present scope; Arnold et al. provide a unified analysis of these methods for elliptic problems [@arnold-brezzi-cockburn-marini2002].