Funciones de forma: el corazón de la interpolación en el MEF.

En el artículo anterior vimos cómo se construye y se ensambla la matriz de rigidez [K]. Pero esa matriz da por hecho algo que todavía no habíamos explicado: el ordenador no resuelve el comportamiento de una pieza en un continuo infinito de puntos, sino que calcula los desplazamientos únicamente en un número discreto de nodos. Los componentes reales, en cambio, se deforman de manera continua en todo su volumen. ¿Cómo sabe entonces el software qué ocurre en los infinitos puntos que hay entre un nodo y otro? La respuesta son las funciones de forma —también llamadas funciones de interpolación— y son, literalmente, el corazón matemático de cualquier elemento finito.

¿Qué son las funciones de forma?

Las funciones de forma son expresiones matemáticas —casi siempre polinomios— que interpolan los valores nodales para describir cómo varía una magnitud de campo, como el desplazamiento u(x,y) o la temperatura T(x,y), en cualquier punto interior del elemento.

Si llamamos u(x) al desplazamiento en el interior de un elemento, éste se expresa en función de las funciones de forma N_i y de los desplazamientos nodales u_i mediante:

u(x) = N₁u₁ + N₂u₂ + … + N_nu_n

donde n es el número de nodos del elemento. Para que una función de forma sea físicamente válida tiene que cumplir dos condiciones estrictas.

Propiedad de Kronecker. La función N_i debe valer exactamente 1 en su propio nodo i, y exactamente 0 en todos los demás nodos del elemento. Esto garantiza que, justo en la posición de un nodo, el desplazamiento calculado coincida con el valor nodal —ni más, ni menos.

Partición de la unidad. En cualquier punto interior del elemento, la suma de todas las funciones de forma debe valer 1:

N₁ + N₂ + … + N_n = 1

Esta propiedad es la que permite que el elemento se mueva como cuerpo rígido —un desplazamiento constante, sin deformación— sin generar tensiones artificiales en su interior. Un elemento cuyas funciones de forma no sumen 1 fallaría el test más básico de la física: desplazar la pieza entera sin tocarla no debería tensionar nada.

Fig. 1 — La función de forma N₂ de un elemento cuadrático de 3 nodos: vale exactamente 1 en su propio nodo y 0 en los otros dos —la propiedad de Kronecker, visible de un vistazo.

Jerarquía de las funciones de forma: lineales, parabólicas y cúbicas

El orden del polinomio elegido para representar el campo de desplazamientos define la familia y el orden del elemento finito. Y como la deformación se obtiene derivando ese polinomio, cada grado de desplazamiento produce un campo de deformación un grado inferior.

Orden	Polinomio de desplazamiento	Campo de deformación/tensión
1er grado	Lineal — u = a₁ + a₂x	Constante
2º grado	Cuadrático — u = a₁ + a₂x + a₃x²	Lineal
3er grado	Cúbico — u = a₁ + a₂x + a₃x² + a₄x³	Cuadrático

Funciones de forma lineales (primer orden)

Se basan en polinomios de primer grado y representan, geométricamente, líneas rectas en 1D o planos inclinados en 2D y 3D entre nodos. Solo tienen nodos en las esquinas o extremos del elemento. Son la familia del elemento de barra de 2 nodos, el triángulo de tensión constante de 3 nodos (CST), el cuadrilátero de 4 nodos (Q4) y el tetraedro de 4 nodos (TET4).

Funciones de forma parabólicas o cuadráticas (segundo orden)

Emplean polinomios de segundo grado —en 2D, del tipo u(x,y) = a₁ + a₂x + a₃y + a₄x² + a₅xy + a₆y²— y para resolver esos coeficientes adicionales necesitan incorporar nodos intermedios (midside nodes) en las aristas del elemento. Son la familia del triángulo de deformación lineal de 6 nodos (LST), el cuadrilátero de 8 nodos (Q8), el tetraedro de 10 nodos (TET10) y el hexaedro de 20 nodos (HEX20). Su ventaja geométrica es directa: sus aristas y caras pueden curvarse físicamente, adaptándose mucho mejor a geometrías orgánicas o redondeadas que sus equivalentes lineales.

Funciones de forma cúbicas (tercer orden)

Utilizan polinomios de tercer grado. Su aplicación más conocida está en los elementos de viga sometidos a flexión: la curva de deflexión v(x) se representa mediante un polinomio cúbico —la interpolación de Hermite—. Como el elemento tiene 4 grados de libertad (un desplazamiento transversal v y una rotación φ = dv/dx por nodo), las funciones de forma de Hermite garantizan la continuidad tanto del desplazamiento como de la pendiente entre elementos adyacentes. También existen elementos planos de orden superior, como el QST (triángulo de deformación cuadrática de 10 nodos) o el Q12, que emplean términos cúbicos completos de la serie de Pascal para representaciones hiperprecisas.

El vínculo matemático entre desplazamiento, deformación y tensión

Para entender por qué el orden de la función de forma importa tanto, hay que mirar la cadena de cálculo interna del solver. El desplazamiento se determina resolviendo:

u(x) = [N]{d}

La deformación, según la teoría de la elasticidad, es la derivada espacial del desplazamiento:

{ε} = du/dx = [dN/dx]{d} = [B]{d}

donde [B] es la matriz de gradientes o matriz deformación-desplazamiento —la misma [B] que ya apareció al construir la matriz de rigidez: [K] = ∫[B]^T[D][B] dV—. Y la tensión se calcula multiplicando la deformación por la matriz constitutiva del material [D] (ley de Hooke):

{σ} = [D]{ε} = [D][B]{d}

Como la deformación y la tensión se obtienen derivando la función de forma, el polinomio que describe la tensión en el interior del elemento es siempre un grado menor que el que describe el desplazamiento. Esa reducción de un grado —y no otra cosa— es la que separa el comportamiento de un elemento lineal del de uno cuadrático.

Por qué los elementos de orden superior capturan mejor los gradientes de tensión

El problema de los elementos lineales. Si u(x) es lineal, su derivada du/dx es una constante. En un elemento lineal —un CST de 3 nodos, un TET4 de 4 nodos— la deformación y la tensión son, por tanto, rigurosamente constantes en todo el volumen del elemento. Colócalo en una zona donde la tensión cambia deprisa —un agujero, un redondeo, una sección en flexión— y el elemento no podrá mostrar esa variación: solo reportará un único valor promedio. Para capturar un gradiente alto con elementos lineales hace falta una malla extremadamente fina, con docenas de elementos apilados aproximando la curva real en forma de escalones.

La potencia de los elementos cuadráticos. Si u(x) es cuadrático, su derivada du/dx es lineal. En un elemento cuadrático —un LST de 6 nodos, un TET10, un HEX20— la tensión varía linealmente a través del propio elemento. Un solo elemento cuadrático puede capturar con fidelidad cómo la tensión sube desde un valor bajo en un extremo hasta un pico en el otro, lo que permite predecir concentraciones de tensión y fenómenos de flexión con mallas mucho más gruesas y un margen de error notablemente menor.

Fig. 2 — El mismo gradiente de tensión real (línea discontinua), aproximado por elementos lineales (escalones constantes) frente a elementos cuadráticos (rampas lineales que enlazan sin discontinuidad). Menos elementos cuadráticos, más fidelidad.

Eliminación del bloqueo por cortante (shear locking)

Cuando elementos lineales cuadriláteros o sólidos de integración completa se someten a flexión, sus bordes rectos no pueden curvarse. Al intentar deformarse, desarrollan tensiones cortantes parásitas que absorben energía e inmovilizan numéricamente el modelo —el bloqueo por cortante—, haciendo que la estructura parezca más rígida de lo que es en realidad. Las funciones de forma cuadráticas permiten que los bordes del elemento se curven físicamente, eliminando por completo este problema.

La tensión no vive en el elemento por derecho propio: es lo que queda después de derivar la función de forma. Un elemento lineal no puede mentir de forma sutil —solo puede promediar; uno cuadrático puede seguir la curva real.

Criterios de compatibilidad y completitud (C⁰ frente a C¹)

Para que los resultados converjan hacia la solución real a medida que se refina la malla, las funciones de forma deben cumplir tres criterios.

Completitud. La función de forma debe poder representar los movimientos de cuerpo rígido y los estados de deformación constante —por eso no se pueden omitir términos polinómicos de orden inferior en favor de términos de orden superior dentro de la serie de Pascal.

Continuidad C⁰. Exige que la variable de campo —el desplazamiento— sea continua a lo largo de las fronteras entre elementos, evitando huecos, grietas o solapamientos al deformarse. Los elementos de barra, membrana, sólido y cáscara son elementos C⁰.

Continuidad C¹. Exige que tanto la variable de campo como su primera derivada —las pendientes o rotaciones— sean continuas a través de los límites del elemento. Es un requisito indispensable en los elementos de flexión de vigas y placas, donde una discontinuidad de pendiente entre elementos vecinos carecería de sentido físico.

Resumen comparativo para la selección en ingeniería

Característica	Elementos lineales (1er orden)	Elementos cuadráticos (2º orden)
Grado del polinomio (u)	1er grado	2º grado
Variación de tensión (σ)	Constante en el elemento	Lineal en el elemento
Ubicación de nodos	Solo en las esquinas	Esquinas y puntos medios
Adaptación a curvas	Bordes rectos	Bordes curvos reales
Sensibilidad a la distorsión	Muy alta	Baja
Casos recomendados	Mallas masivas, contactos	Gradientes de tensión, flexión, TET10

La matriz de rigidez, el mallado y las funciones de forma no son tres temas separados: son la misma idea vista desde tres ángulos. Entender qué polinomio hay detrás de cada elemento —y qué grado tiene su derivada— es lo que separa elegir una malla por instinto de elegirla con criterio. Si trabajas con tetraedros en zonas críticas con gradientes de tensión o flexión, evita el TET4 lineal y usa TET10 cuadrático: un número menor de elementos cuadráticos entrega una solución más precisa y fiable que un modelo saturado con millones de elementos lineales.

Queda un último eslabón en esta cadena: cómo integra el software estas funciones de forma sobre el volumen del elemento —la cuadratura de Gauss— y qué ocurre cuando se elige integración reducida en lugar de completa, el origen del temido hourglassing. Lo veremos cuando hablemos de cuándo confiar en un mallado.

FIN

In the previous article we saw how the stiffness matrix [K] is built and assembled. But that matrix takes something for granted that we hadn't explained yet: the computer doesn't solve a part's behaviour over an infinite continuum of points — it only computes displacements at a discrete number of nodes. Real components, on the other hand, deform continuously throughout their entire volume. So how does the software know what happens at the infinite points between one node and the next? The answer is shape functions — also called interpolation functions — and they are, quite literally, the mathematical heart of every finite element.

What are shape functions?

Shape functions are mathematical expressions — almost always polynomials — that interpolate nodal values to describe how a field quantity, such as displacement u(x,y) or temperature T(x,y), varies at any interior point of the element.

If we call u(x) the displacement inside an element, it is expressed in terms of the shape functions N_i and the nodal displacements u_i as:

u(x) = N₁u₁ + N₂u₂ + … + N_nu_n

where n is the number of nodes in the element. For a shape function to be physically valid it must satisfy two strict conditions.

Kronecker property. N_i must equal exactly 1 at its own node i, and exactly 0 at every other node in the element. This guarantees that, right at a node's location, the computed displacement matches the nodal value — no more, no less.

Partition of unity. At any interior point of the element, the sum of all shape functions must equal 1:

N₁ + N₂ + … + N_n = 1

This property is what lets the element undergo rigid-body motion — a constant displacement, with no deformation — without generating artificial internal stresses. An element whose shape functions don't sum to 1 would fail the most basic physical test there is: moving the whole part without touching it shouldn't stress anything.

Fig. 1 — The shape function N₂ of a 3-node quadratic element: it equals exactly 1 at its own node and 0 at the other two — the Kronecker property, visible at a glance.

Hierarchy of shape functions: linear, parabolic and cubic

The order of the polynomial chosen to represent the displacement field defines the finite element's family and order. And since strain comes from differentiating that polynomial, every degree of displacement produces a strain field one degree lower.

Order	Displacement polynomial	Strain/stress field
1st degree	Linear — u = a₁ + a₂x	Constant
2nd degree	Quadratic — u = a₁ + a₂x + a₃x²	Linear
3rd degree	Cubic — u = a₁ + a₂x + a₃x² + a₄x³	Quadratic

Linear shape functions (first order)

These are based on first-degree polynomials and represent, geometrically, straight lines in 1D or tilted planes in 2D and 3D between nodes. They only have nodes at the corners or ends of the element. This is the family of the 2-node bar element, the 3-node constant-strain triangle (CST), the 4-node quadrilateral (Q4) and the 4-node tetrahedron (TET4).

Parabolic or quadratic shape functions (second order)

These use second-degree polynomials — in 2D, of the form u(x,y) = a₁ + a₂x + a₃y + a₄x² + a₅xy + a₆y² — and solving for those extra coefficients requires adding midside nodes along the element's edges. This is the family of the 6-node linear-strain triangle (LST), the 8-node quadrilateral (Q8), the 10-node tetrahedron (TET10) and the 20-node hexahedron (HEX20). Their geometric advantage is direct: their edges and faces can physically curve, adapting far better to organic or rounded geometry than their linear equivalents.

Cubic shape functions (third order)

These use third-degree polynomials. Their best-known application is in beam elements under bending: the deflection curve v(x) is represented by a cubic polynomial — Hermite interpolation. Since the element has 4 degrees of freedom (a transverse displacement v and a rotation φ = dv/dx per node), Hermite shape functions guarantee continuity of both displacement and slope between adjacent elements. Higher-order planar elements also exist, such as the QST (10-node quadratic-strain triangle) or the Q12, which use complete cubic terms of Pascal's triangle for hyper-accurate representations.

The mathematical link between displacement, strain and stress

To understand why the order of the shape function matters so much, look at the solver's internal calculation chain. Displacement is determined by solving:

u(x) = [N]{d}

Strain, from elasticity theory, is the spatial derivative of displacement:

{ε} = du/dx = [dN/dx]{d} = [B]{d}

where [B] is the gradient, or strain-displacement, matrix — the very same [B] that already showed up when building the stiffness matrix: [K] = ∫[B]^T[D][B] dV. And stress is obtained by multiplying strain by the material's constitutive matrix [D] (Hooke's law):

{σ} = [D]{ε} = [D][B]{d}

Because strain and stress come from differentiating the shape function, the polynomial describing stress inside the element is always one degree lower than the one describing displacement. That one-degree reduction — and nothing else — is what separates a linear element's behaviour from a quadratic one's.

Why higher-order elements capture stress gradients better

The problem with linear elements. If u(x) is linear, its derivative du/dx is a constant. In a linear element — a 3-node CST, a 4-node TET4 — strain and stress are therefore rigorously constant throughout the element's volume. Place it in a region where stress changes quickly — a hole, a fillet, a section under bending — and the element simply cannot show that variation: it can only report a single averaged value. Capturing a steep gradient with linear elements requires an extremely fine mesh, with dozens of stacked elements approximating the real curve as a staircase.

The power of quadratic elements. If u(x) is quadratic, its derivative du/dx is linear. In a quadratic element — a 6-node LST, a TET10, a HEX20 — stress varies linearly across the element itself. A single quadratic element can faithfully capture how stress climbs from a low value at one end to a peak at the other, letting you predict stress concentrations and bending behaviour with much coarser meshes and a noticeably smaller margin of error.

Fig. 2 — The same real stress gradient (dashed line), approximated by linear elements (constant steps) versus quadratic elements (linear ramps that join without a break). Fewer quadratic elements, more fidelity.

Eliminating shear locking

When fully-integrated linear quadrilateral or solid elements are subjected to bending, their straight edges cannot curve. As they try to deform, they develop parasitic shear stresses that absorb energy and numerically stiffen the model — shear locking — making the structure appear stiffer than it actually is. Quadratic shape functions let the element's edges physically curve, eliminating this problem entirely.

Stress doesn't live in the element in its own right: it's what's left after differentiating the shape function. A linear element can't lie subtly — it can only average; a quadratic one can follow the real curve.

Compatibility and completeness criteria (C⁰ vs. C¹)

For results to converge to the real solution as the mesh is refined, shape functions must satisfy three criteria.

Completeness. The shape function must be able to represent rigid-body motion and states of constant strain — which is why lower-order polynomial terms can't be dropped in favour of higher-order ones within Pascal's triangle.

C⁰ continuity. Requires the field variable — displacement — to be continuous across element boundaries, avoiding gaps, cracks or overlaps as the mesh deforms. Bar, membrane, solid and shell elements are all C⁰ elements.

C¹ continuity. Requires both the field variable and its first derivative — slopes or rotations — to be continuous across element boundaries. This is an indispensable requirement for beam and plate bending elements, where a slope discontinuity between neighbouring elements would make no physical sense.

Comparative summary for engineering selection

Feature	Linear elements (1st order)	Quadratic elements (2nd order)
Polynomial degree (u)	1st degree	2nd degree
Stress variation (σ)	Constant within element	Linear within element
Node location	Corners only	Corners and midsides
Curve adaptation	Straight edges	Real curved edges
Distortion sensitivity	Very high	Low
Recommended cases	Massive meshes, contact	Stress gradients, bending, TET10

The stiffness matrix, the mesh and the shape functions aren't three separate topics: they're the same idea seen from three angles. Understanding which polynomial sits behind each element — and what degree its derivative has — is what separates choosing a mesh by instinct from choosing it with judgement. If you work with tetrahedra in critical zones with stress gradients or bending, avoid the linear TET4 and use a quadratic TET10 instead: fewer quadratic elements deliver a more accurate, more reliable solution than a model saturated with millions of linear ones.

One link remains in this chain: how the software integrates these shape functions over the element's volume — Gauss quadrature — and what happens when reduced integration is chosen instead of full integration, the origin of the dreaded hourglassing. We'll get to that when we talk about when to trust a mesh.

END