Перейдём к вычислению матрицы жёсткости для трёхмерного случая. В трёхмерном пространстве градиент имеет вид $\nabla \upsilon = \left( \frac{\displaystyle \partial \upsilon}{\displaystyle \partial x}, \frac{\displaystyle \partial \upsilon}{\displaystyle \partial y}, \frac{\displaystyle \partial \upsilon}{\displaystyle \partial z} \right)$, а скалярное произведение градиентов соответственно $\nabla \upsilon \cdot \nabla \upsilon = \left( \frac{\displaystyle \partial \upsilon}{\displaystyle \partial x} \right)^2 + \left( \frac{\displaystyle \partial \upsilon}{\displaystyle \partial y} \right)^2 + \left( \frac{\displaystyle \partial \upsilon}{\displaystyle \partial z} \right)^2$. Учитывая, что область $M$ представляет собой координатное пространство, разбитое на симплексы-тетраэдры, исследуемую часть функционала для одного тетраэдра с вершинами

можно записать как

Функция $\upsilon(x, y, z) = \sum_{i=1}^N \upsilon_i(x, y, z)$. Пробная функция на тетраэдре имеет вид $\upsilon_{(i)(i+3)}(x, y, z) = q_i \cdot \phi_i(x, y, z) + q_{i+1} \cdot \phi_{i+1}(x, y, z) + q_{i+2} \cdot \phi_{i+2}(x, y, z) + q_{i+3} \cdot \phi_{i+3}(x, y, z)$. Аналогично двумерному случаю, функции <<крышек>> для тетраэдра имеют линейный вид

Вычислим частные производные пробной функции

Заметим, что производные не зависят от $x$, $y$ и $z$ и являются константами на тетраэдре. Подставим (6.15$$\begin{cases} \frac{\displaystyle \partial \upsilon_{(i)(i+3)}(x, y, z)}{\displaystyle \partial x} = q_i \cdot b_i + q_{i+1} \cdot b_{i+1} + q_{i+2} \cdot b_{i+2} + q_{i+3} \cdot b_{i+3}\\ \frac{\displaystyle \partial \upsilon_{(i)(i+3)}(x, y, z)}{\displaystyle \partial y} = q_i \cdot c_i + q_{i+1} \cdot c_{i+1} + q_{i+2} \cdot c_{i+2} + q_{i+3} \cdot c_{i+3}\\ \frac{\displaystyle \partial \upsilon_{(i)(i+3)}(x, y, z)}{\displaystyle \partial z} = q_i \cdot d_i + q_{i+1} \cdot d_{i+1} + q_{i+2} \cdot d_{i+2} + q_{i+3} \cdot d_{i+3} \end{cases}$$) в (6.13$$\int_{\text{тет}} \left[ \left( \frac{\displaystyle \partial \upsilon}{\displaystyle \partial x} \right)^2 + \left( \frac{\displaystyle \partial \upsilon}{\displaystyle \partial y} \right)^2 + \left( \frac{\displaystyle \partial \upsilon}{\displaystyle \partial z} \right)^2 \right] \,dV.$$)

где $V_{\text{тет}}$ — объём тетраэдра, который вычисляется по формуле

где $\Delta$ — определитель матрицы

Раскроем квадраты и перегруппируем члены

Коэффициенты $b_k$, $c_k$ и $d_k$ для каждой вершины $k$ тетраэдра вычисляются через миноры определителя $\Delta$. Для вершины с индексом $k$ коэффициенты имеют вид

где $\Delta_k^{(x)}$, $\Delta_k^{(y)}$ и $\Delta_k^{(z)}$ — миноры, получаемые вычёркиванием $k$-й строки и соответствующего столбца ($x$, $y$ или $z$) из матрицы (6.17$$\begin{aligned}\Delta = \begin{vmatrix} x_i & y_i & z_i & 1\\ x_{i+1} & y_{i+1} & z_{i+1} & 1\\ x_{i+2} & y_{i+2} & z_{i+2} & 1\\ x_{i+3} & y_{i+3} & z_{i+3} & 1 \end{vmatrix}.\end{aligned}$$).

Введём обозначения для элементов локальной матрицы жёсткости тетраэдра

Таким образом, локальная матрица жёсткости для тетраэдрального элемента имеет вид

Локальная матрица жёсткости является симметричной, то есть $k_{mn}^{\text{лок}} = k_{nm}^{\text{лок}}$. Глобальная матрица жёсткости $\mathbf{K}$ получается путём суммирования вкладов от всех тетраэдральных элементов сетки методом сборки: элементы локальных матриц добавляются к соответствующим элементам глобальной матрицы согласно глобальной нумерации узлов. Размерность глобальной матрицы жёсткости равна $N \times N$, где $N$ — общее количество узлов сетки.