Теперь вычислим вектор нагрузки для трёхмерного случая. Вектор нагрузки связан с интегралом произведения функции источника на пробную функцию. Рассмотрим интеграл для одного тетраэдра с вершинами

Пробная функция на тетраэдре имеет вид $\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)$. Для вычисления интеграла необходимо интерполировать функцию источника $f(x, y, z)$ на тетраэдре линейной функцией. Представим $f(x, y, z)$ в виде $\widetilde{f}(x, y, z) = e_1 + e_2 \cdot x + e_3 \cdot y + e_4 \cdot z$. Коэффициенты $e_1, e_2, e_3, e_4$ определяются из условий интерполяции в узлах

Поскольку функции <<крышки>> $\phi_n$ равны единице в своём узле и нулю в остальных, линейный интерполянт функции источника, выраженный через них, принимает вид

где $f_i, f_{i+1}, f_{i+2}, f_{i+3}$ — значения функции источника в вершинах тетраэдра.

Подставим интерполированную функцию источника и пробную функцию в (6.48$$\int_{\text{тет}} f(x, y, z) \cdot \upsilon \,dV.$$). Примем во внимание (6.32$$\begin{cases} \phi_i(x, y, z) = a_i + b_i \cdot x + c_i \cdot y + d_i \cdot z\\ \phi_{i+1}(x, y, z) = a_{i+1} + b_{i+1} \cdot x + c_{i+1} \cdot y + d_{i+1} \cdot z\\ \phi_{i+2}(x, y, z) = a_{i+2} + b_{i+2} \cdot x + c_{i+2} \cdot y + d_{i+2} \cdot z\\ \phi_{i+3}(x, y, z) = a_{i+3} + b_{i+3} \cdot x + c_{i+3} \cdot y + d_{i+3} \cdot z \end{cases}$$) из раздела о функциях <<крышек>>

где индексы $m$ и $n$ пробегают вершины тетраэдра $\{i, i+1, i+2, i+3\}$.

Воспользуемся соотношениями (6.33$$\begin{split} &\int_{\text{тет}} \phi_m \cdot \phi_n \,dV = \begin{cases} \frac{\displaystyle V_{\text{тет}}}{\displaystyle 10}, & m = n\\ \frac{\displaystyle V_{\text{тет}}}{\displaystyle 20}, & m \neq n \end{cases} \end{split}$$) для интегралов произведений функций <<крышек>>, выведенными при вычислении матрицы демпфирования, где $V_{\text{тет}}$ — объём тетраэдра, который вычисляется по формуле (6.16$$V_{\text{тет}} = \frac{\displaystyle |\Delta|}{\displaystyle 6},$$). После приведения подобных получаем

Введём обозначения для элементов локального вектора нагрузки тетраэдра

Таким образом, локальный вектор нагрузки для трёхмерного тетраэдрального элемента имеет вид

Заметим, что полученный вектор нагрузки совпадает с произведением локальной матрицы демпфирования (6.35$$\begin{aligned}\mathbf{C}_{\text{тет}} = \begin{bmatrix} c_{(i)(i)} & c_{(i)(i+1)} & c_{(i)(i+2)} & c_{(i)(i+3)}\\ c_{(i+1)(i)} & c_{(i+1)(i+1)} & c_{(i+1)(i+2)} & c_{(i+1)(i+3)}\\ c_{(i+2)(i)} & c_{(i+2)(i+1)} & c_{(i+2)(i+2)} & c_{(i+2)(i+3)}\\ c_{(i+3)(i)} & c_{(i+3)(i+1)} & c_{(i+3)(i+2)} & c_{(i+3)(i+3)} \end{bmatrix} = \frac{\displaystyle V_{\text{тет}}}{\displaystyle 20} \begin{bmatrix} 2 & 1 & 1 & 1\\ 1 & 2 & 1 & 1\\ 1 & 1 & 2 & 1\\ 1 & 1 & 1 & 2 \end{bmatrix}.\end{aligned}$$) на вектор узловых значений источника, то есть $\mathbf{R}_{\text{тет}} = \mathbf{C}_{\text{тет}} \cdot \mathbf{f}$.

Глобальный вектор нагрузки $\mathbf{R}$ получается путём суммирования вкладов от всех тетраэдральных элементов сетки методом сборки: элементы локальных векторов добавляются к соответствующим элементам глобального вектора согласно глобальной нумерации узлов. Размерность глобального вектора нагрузки равна $N$, где $N$ — общее количество узлов сетки. Важно отметить, что для внутренних узлов происходит суммирование вкладов от всех смежных тетраэдральных элементов, содержащих данный узел. Для граничных узлов вектор нагрузки определяется граничными условиями задачи.