Finite Element Method by Example in Qt/C++
Author
Krzysztof Napiontek
Last Updated
11 years ago
License
Creative Commons CC BY 4.0
Abstract
Learn FEM by example in a few steps.
\documentclass[a4paper]{article}
\usepackage[english]{babel}
\usepackage{amsmath}
\usepackage{tikz}
\usepackage{listings}
\usepackage{pxfonts}
\usepackage[margin=1in]{geometry}
\usepackage[nottoc]{tocbibind}
\usepackage{xcolor}
\usepackage{float}
\setcounter{MaxMatrixCols}{12}
\lstset{
belowcaptionskip=1\baselineskip,
breaklines=true,
xleftmargin=\parindent,
language=C++,
showstringspaces=false,
basicstyle=\ttfamily\small,
keywordstyle=\bfseries\color{green!40!black},
commentstyle=\itshape\color{purple!40!black},
identifierstyle=\color{blue},
stringstyle=\color{orange},
tabsize=2,
}
\title{The Finite Element Method by Example in Qt/C++}
\author{Krzysztof Napiontek\protect\\knapiontek@gmail.com}
\date{Version 1.0\protect\\\today}
\begin{document}
\maketitle
\begin{abstract}
The aim of this article is help with understanding The Finite Element Method in mechanics. The theoretical part is reduced to its minimum and most of the work is focused on diagrams, graphics and examples. The relevant part of source code calculating tetrahedron truss is presented at the end of each section along with its visualized output.
\end{abstract}
\tableofcontents
\pagebreak
\section{Introduction}
There are many works presenting FEM, most of them very complex and comprehensive with purpose of full explanation of all details of the method. A reader who needs explanation of general idea behind finds it difficult to understand. This work tries to extract the most important aspects in a way helping to understand all knowledge almost at the first read. The relevant source code added at the end of each section calculates and visualizes the simplest 3D truss object - tetrahedron. The Hook's Law is the only equation borrowed from mechanics. There is also a need to understand trigonometry and basic matrix operations before reading.
\section{The Concrete Example}
For the purpose of this article the complete concrete example was implemented. It includes data type definitions, construction of a matrix equation, its solution and in final the presentation in graphics format. A number of dumps of the equation data have been generated for better understanding parts of implementation. The source code presented in this work was tested with Qt 5.1 in Windows environment. It does not use any platform specific code and can be easily ported to other systems.
\bigskip
The final result of the example is presented on the figure \ref{fig:final:graphics}. It contains an original geometry of the truss, also geometry deformation by an external force and reactions in fixed points. Detailed explanation will be given in following chapters.
\begin{figure}[H]
\centering
\begin{tikzpicture}
\input{figure-final}
\node at (1.5,.5){$p_0$};
\node at (.8,4){$p_1$};
\node at (6.1,-.2){$p_2$};
\node at (.5,-.7){$p_3$};
\draw (4,5) -- (5,5) node[right]{original geometry};
\draw[dashed] (4,4.5) -- (5,4.5) node[right]{geometry deformation};
\draw[->,red] (4,4) -- (5,4) node[right]{applied external force};
\draw[->,red] (4,3.5) -- (5,3.5) node[right]{reactions in fixed points};
\draw[red,fill](4,3.5) circle[radius=.1];
\end{tikzpicture}
\caption{The Graphics Generated by the Example}
\label{fig:final:graphics}
\end{figure}
\subsection{Global Data Definition in the Example}
\lstinputlisting{data.cpp}
\section{The Finite Element Method}
The FEM analyse in this chapter finds relation between forces and displacements in global coordinates. At first it breaks down the truss into separate elements, finds relation of them in local coordinates using the Hook's Law and in next step it transforms it to global coordinate system for each element. Finally it aggregates all elements in one global matrix of stiffness. The lower case bold font is reserved for local and the upper case bold font for global matrices and vectors.
\subsection{Tetrahedron Truss Breakdown}
The first step of analyses requires the truss breakdown. Each element will be analysed separately and finally all of them will be aggregated in global matrix. Separation of elements is shown on the figure \ref{fig:breakdown}.
\begin{figure}[H]
\centering
\begin{tikzpicture}[scale=0.7]
\input{figure-breakdown}
\node at (1.3,.3){$p_0$};
\node at (1.2,5.2){$p_1$};
\node at (7.2,0.1){$p_2$};
\node at (-0.1,-1.1){$p_3$};
\node[red] at (1.5,2.0){$e_0$};
\node[red] at (3.7,3.2){$e_1$};
\node[red] at (3.5,0.5){$e_2$};
\node[red] at (0.6,-0.1){$e_3$};
\node[red] at (-0.1,2.0){$e_4$};
\node[red] at (3.5,-0.5){$e_5$};
\end{tikzpicture}
\caption{Truss Breakdown}
\label{fig:breakdown}
\end{figure}
\subsection{The Stiffness Matrix of an Element in Local Coordinates}
Relation between axial forces $f_0, f_1$ and axial displacements $\Delta p_0, \Delta p_1$ in local coordinates in respect of the Hook's Law as depicted on the figure \ref{fig:displacement:and:force} is described by equations \ref{eq:hook0} and \ref{eq:hook1} . Despite of 3D analyses most of the figures use 2D coordinates for simplicity.
\begin{equation}
f_0 = \frac {EA} L (\Delta p_0 - \Delta p_1)
\label{eq:hook0}
\end{equation}
\begin{equation}
f_1 = \frac {EA} L (\Delta p_1 - \Delta p_0)
\label{eq:hook1}
\end{equation}
\begin{figure}[H]
\centering
\begin{tikzpicture}[scale=.7]
%%% displacements
\node at (2,9){a) Displacements};
% global arrows
\draw[<->](0,8) node[right] {$y$} -- (0,0) -- (9,0) node[below] {$x$};
% auxyliary lines and local arrows
\draw[dashed](1,3) -- (4,3) -- (4,1);
\draw[<->](1,3) -- node[left] {$\Delta y_0$} (1,1) -- node[above] {$\Delta x_0$} (4,1);
\draw[->](1,1) -- node[above] {$\Delta p_0$} (4,3);
\draw[dashed](4,7.5) -- (8,7.5) -- (8,6);
\draw[<->](4,7.5) -- node[left] {$\Delta y_1$} (4,6) -- node[above] {$\Delta x_1$} (8,6);
\draw[->](4,6) -- node[above] {$\Delta p_1$} (8,7.5);
% original and displaced elements
\draw[thick,red](1,1) -- node[above left]{$L$} (4,6);
\draw[fill](1,1) circle[radius=.1];
\draw[fill](4,6) circle[radius=.1];
\draw[thick,red,dashed](4,3) -- node[below right]{$L'$} (8,7.5);
\draw[fill](4,3) circle[radius=.1];
\draw[fill](8,7.5) circle[radius=.1];
% labels of points
\node at (1.2,.5) {$(x_0,y_0)$};
\node at (4.9,5.5) {$(x_1,y_1)$};
%%% forces
\node at (13,9){b) Forces};
% global arrows
\draw[<->](12,8) node[right] {$y$} -- (12,0) -- (21,0) node[below] {$x$};
% auxyliary lines and local arrows
\draw[dashed](13,3) -- (13,1) -- (15,1);
\draw[<->](13,3) -- node[above] {$f_{x_0}$} (15,3) -- node[right] {$f_{y_0}$} (15,1);
\draw[->](15,3) -- node[below] {$f_0$} (13,1);
\draw[dashed](18,8) -- (20,8) -- (20,6);
\draw[<->](18,8) -- node[left] {$f_{y_1}$} (18,6) -- node[above] {$f_{x_1}$} (20,6);
\draw[->](18,6) -- node[above] {$f_1$} (20,8);
% element
\draw[thick,red,dashed](15,3) -- node[above left]{$L'$} (18,6);
\draw[fill](15,3) circle[radius=.1];
\draw[fill](18,6) circle[radius=.1];
% labels of points
\node at (16,2.8) {$(x_0,y_0)$};
\node at (18.6,5.5) {$(x_1,y_1)$};
\end{tikzpicture}
\caption{Displacements (a) and Stress Forces (b) in an Element}
\label{fig:displacement:and:force}
\end{figure}
In fact the method provides only approximated solution since it assumes that angles between elements stay the same after applying external forces. In reality such approach is accepted as long as deformations are small.
\bigskip
Using matrix notation we isolate local stiffness matrix \textbf{k}.
\begin{equation}
\begin{bmatrix}
f_0 \\
f_1
\end{bmatrix}
= \frac {EA} L
\begin{bmatrix}
1 & -1 \\
-1 & 1
\end{bmatrix}
\begin{bmatrix}
\Delta p_0 \\
\Delta p_1
\end{bmatrix}
\end{equation}
\begin{equation}
\mathbf k = \frac {EA} L
\begin{bmatrix}
1 & -1 \\
-1 & 1
\end{bmatrix}
\end{equation}
The stiffness matrix represents geometry in an algebraic form. Final form of matrix equation in local coordinates shown below. Lower case bold letters reserved for local coordinates.
\begin{equation}
\mathbf{f} = \mathbf{k} \Delta \mathbf{p}
\end{equation}
\subsection{Coordinate Transformation}
In previous chapter we described relations in local coordinates. In order to aggregate all elements into global system a number of transformations have to take place.
\subsubsection{Global and Local Coordinates Relations}
The figure \ref{fig:displacement:and:force}a can be used in order to describe local and global geometrical relations.
\begin{equation}
L = \sqrt { (x_1 - x_0)^2 + (y_1 - y_0)^2 + (z_1 - z_0)^2 }
\end{equation}
\begin{equation}
c_x = \cos x = \frac {x_1 - x_0} L
\end{equation}
\begin{equation}
c_y = \cos y = \frac {y_1 - y_0} L
\end{equation}
\begin{equation}
c_z = \cos z = \frac {z_1 - z_0} L
\end{equation}
\subsubsection{The Point Displacement Relations}
The displacement from the figure \ref{fig:displacement:and:force}a for the point $p_0$ is broken down on the figure \ref{fig:displacement:details}.
\begin{figure}[H]
\centering
\begin{tikzpicture}[scale=.5]
% global arrows
\draw[<->](0,8) node[right] {$y$} -- (0,0) -- (10,0) node[below] {$x$};
% auxyliary lines and local arrows
\draw[dashed](2,6) -- (8,6) -- (8,2);
\draw[<->](2,6) -- node[left] {$\Delta y_0$} (2,2) -- node[above] {$\Delta x_0$} (8,2);
% displacement
\draw[->,red](2,2) -- (8,6);
\draw[red,dashed](8,2) -- (6.15,4.8);
% labels
\node at (3,2.3) {$\alpha$};
\node at (7.1,2.3) {$\beta$};
\node at (7.7,3.2) {$\alpha$};
\node at (7.7,5.1) {$\beta$};
\node at (5.9,4.0) {$90^o$};
\node at (4.9,4.5) {$\Delta p_0$};
\draw[fill](2,2) circle[radius=.1] node[below]{$(x_0,y_0)$};
\end{tikzpicture}
\caption{Displacement of the Point $p_0$ in Details}
\label{fig:displacement:details}
\end{figure}
From congruence of triangles as on the figure \ref{fig:displacement:details} we obtain
\begin{equation}
\Delta p_0 = \Delta x_0 \cos x + \Delta y_0 \cos y + \Delta z_0 \cos z
\end{equation}
\begin{equation}
\Delta p_1 = \Delta x_1 \cos x + \Delta y_1 \cos y + \Delta z_1 \cos z
\end{equation}
In other words
\begin{equation}
\begin{bmatrix}
\Delta p_0 \\
\Delta p_1
\end{bmatrix}
=
\begin{bmatrix}
c_x & c_y & c_z & 0 & 0 & 0 \\
0 & 0 & 0 & c_x & c_y & c_z
\end{bmatrix}
\begin{bmatrix}
\Delta x_0 \\
\Delta y_0 \\
\Delta z_0 \\
\Delta x_1 \\
\Delta y_1 \\
\Delta z_1
\end{bmatrix}
\end{equation}
or simpler
\begin{equation}
\Delta \mathbf{p} = \mathbf{T} \Delta \mathbf{P}
\label{eq:displacement:local:vs:global}
\end{equation}
The equation \ref{eq:displacement:local:vs:global} describes displacement relation between local and global coordinates.
\subsubsection{Global and Local Forces Relations}
Forces given directly from geometry on the figure \ref{fig:displacement:and:force}b
\begin{equation}
f_{x_0} = c_x f_0 \\
\end{equation}
\begin{equation}
f_{y_0} = c_y f_0 \\
\end{equation}
\begin{equation}
f_{z_0} = c_z f_0 \\
\end{equation}
\begin{equation}
f_{x_1} = c_x f_1 \\
\end{equation}
\begin{equation}
f_{y_1} = c_y f_1 \\
\end{equation}
\begin{equation}
f_{z_1} = c_z f_1
\end{equation}
or in a matrix form
\begin{equation}
\begin{bmatrix}
f_{x_0} \\
f_{y_0} \\
f_{z_0} \\
f_{x_1} \\
f_{y_1} \\
f_{z_1}
\end{bmatrix}
=
\begin{bmatrix}
c_x & 0 \\
c_y & 0 \\
c_z & 0 \\
0 & c_x \\
0 & c_y \\
0 & c_z
\end{bmatrix}
\begin{bmatrix}
f_0 \\
f_1
\end{bmatrix}
\end{equation}
\begin{equation}
\mathbf{F = T^T f}
\label{eq:forces:local:vs:global}
\end{equation}
The equation \ref{eq:forces:local:vs:global} describes forces relation between local and global coordinates.
\subsection{The Stiffness Matrix of an Element in Global Coordinates}
Equations obtained in previous chapter: Hook's law, displacement and force relations between local and global coordinates are repeated below.
\begin{equation}
\mathbf{f} = \mathbf{k} \Delta \mathbf{p}
\end{equation}
\begin{equation}
\Delta \mathbf{p} = \mathbf{T} \Delta \mathbf{P}
\end{equation}
\begin{equation}
\mathbf{F = T^T f}
\end{equation}
It can be transformed into
\begin{equation}
\mathbf{f} = \mathbf{k} \Delta \mathbf{p}
\end{equation}
\begin{equation}
\mathbf{f} = \mathbf{k T} \Delta \mathbf{P}
\end{equation}
\begin{equation}
\mathbf{T^T f} = \mathbf{T^T k T} \Delta \mathbf{P}
\end{equation}
\begin{equation}
\mathbf{F} = \mathbf{T^T k T} \Delta \mathbf{P}
\end{equation}
Finally the extracted version of global stiffness matrix \textbf{K} in global coordinates takes form
\begin{equation}
\mathbf{K} = \mathbf{T^T k T} = \frac {EA} L
\begin{bmatrix}
c^2_x & c_xc_y & c_xc_z & -c^2_x & -c_xc_y & -c_xc_z \\
c_xc_y & c^2_y & c_yc_z & -c_xc_y & -c^2_y & -c_yc_z \\
c_xc_z & c_yc_z & c^2_z & -c_xc_z & -c_yc_z & -c^2_z \\
-c^2_x & -c_xc_y & -c_xc_z & c^2_x & c_xc_y & c_xc_z \\
-c_xc_y & -c^2_y & -c_yc_z & c_xc_y & c^2_y & c_yc_z \\
-c_xc_z & -c_yc_z & -c^2_z & c_xc_z & c_yc_z & c^2_z
\end{bmatrix}
\end{equation}
so the main equation of this work is formed as
\begin{equation}
\mathbf{F} = \mathbf{K} \Delta \mathbf{P}
\label{eq:main}
\end{equation}
and describes relation of global forces and deformation of whole truss geometry.
\subsection{Global Stiffness Matrix Aggregation}
Let's form equation \ref{eq:main} for an element $e_2(p_0,p_2)$ in Global Stiffness Matrix. The element index for cosines is omitted for simplicity ($c_x$ instead of $c_{2_x}$ etc.).
\begin{equation}
\frac {EA} L
\begin{bmatrix}
c^2_x & c_xc_y & c_xc_z & . & . & . & -c^2_x & -c_xc_y & -c_xc_z & . & . & . \\
c_xc_y & c^2_y & c_yc_z & . & . & . & -c_xc_y & -c^2_y & -c_yc_z & . & . & . \\
c_xc_z & c_yc_z & c^2_z & . & . & . & -c_xc_z & -c_yc_z & -c^2_z & . & . & . \\
. & . & . & . & . & . & . & . & . & . & . & . \\
. & . & . & . & . & . & . & . & . & . & . & . \\
. & . & . & . & . & . & . & . & . & . & . & . \\
-c^2_x & -c_xc_y & -c_xc_z & . & . & . & c^2_x & c_xc_y & c_xc_z & . & . & . \\
-c_xc_y & -c^2_y & -c_yc_z & . & . & . & c_xc_y & c^2_y & c_yc_z & . & . & . \\
-c_xc_z & -c_yc_z & -c^2_z & . & . & . & c_xc_z & c_yc_z & c^2_z & . & . & . \\
. & . & . & . & . & . & . & . & . & . & . & . \\
. & . & . & . & . & . & . & . & . & . & . & . \\
. & . & . & . & . & . & . & . & . & . & . & .
\end{bmatrix}
\begin{bmatrix}
\Delta x_0 \\
\Delta y_0 \\
\Delta z_0 \\
. \\
. \\
. \\
\Delta x_2 \\
\Delta y_2 \\
\Delta z_2 \\
. \\
. \\
.
\end{bmatrix}
=
\begin{bmatrix}
f_{x_0} \\
f_{y_0} \\
f_{z_0} \\
. \\
. \\
. \\
f_{x_2} \\
f_{y_2} \\
f_{z_2} \\
. \\
. \\
.
\end{bmatrix}
\label{eq:detail:K:element}
\end{equation}
By substituting parts of above equation by
\begin{equation}
k_2 = \frac {EA} L
\begin{bmatrix}
c^2_x & c_xc_y & c_xc_z \\
c_xc_y & c^2_y & c_yc_z \\
c_xc_z & c_yc_z & c^2_z
\end{bmatrix}
\;\;\;\;\;
\Delta p_i =
\begin{bmatrix}
\Delta x_i \\
\Delta y_i \\
\Delta z_i
\end{bmatrix}
\;\;\;\;\;
f_i =
\begin{bmatrix}
f_{x_i} \\
f_{y_i} \\
f_{z_i}
\end{bmatrix}
\end{equation}
we obtain simpler form ready for further analysis
\begin{equation}
\begin{bmatrix}
k_2 & . & -k_2 & . \\
. & . & . & . \\
-k_2 & . & k_2 & . \\
. & . & . & .
\end{bmatrix}
\begin{bmatrix}
\Delta p_0 \\
. \\
\Delta p_2 \\
.
\end{bmatrix}
=
\begin{bmatrix}
f_0 \\
. \\
f_2 \\
.
\end{bmatrix}
\end{equation}
\begin{equation}
\mathbf{K_2} \Delta \mathbf{P} = \mathbf{F_2}
\end{equation}
The sum of all directional forces based on the condition of equilibrium in a point is equal zero
\begin{equation}
\sum \mathbf{F_x} = 0\;\;\;\;\;\sum \mathbf{F_y} = 0\;\;\;\;\;\sum \mathbf{F_z} = 0
\end{equation}
\begin{figure}[H]
\centering
\begin{tikzpicture}
\input{figure-point0-forces}
\node at (3.3,1.7) {$p_0$};
\end{tikzpicture}
\caption{Force Equilibrium in the Point $p_0$ (reactions and stresses)}
\end{figure}
and displacement in the point is the same for all elements based on it. So we can write down
\begin{equation}
\sum\limits_{e=1}^n \mathbf{K_e} \Delta \mathbf{P} = \sum\limits_{e=1}^n \mathbf{F_e}\;\;\;\;\;\mbox{e - index of an element}
\end{equation}
Finally the simplified form of \textbf{K} for all elements (aggregation with collocation)
\begin{equation}
\begin{bmatrix}
k_0 + \textcolor{green}{k_2} + \textcolor{blue}{k_3} & -k_0 & \textcolor{green}{-k_2} & \textcolor{blue}{-k_3} \\
-k_0 & k_0 + \textcolor{red}{k_1} + \textcolor{cyan}{k_4} & \textcolor{red}{-k_1} & \textcolor{cyan}{-k_4} \\
\textcolor{green}{-k_2} & \textcolor{red}{-k_1} & \textcolor{red}{k_1} + \textcolor{green}{k_2} + k_5 & -k_5 \\
\textcolor{blue}{-k_3} & \textcolor{cyan}{-k_4} & -k_5 & \textcolor{blue}{k_3} + \textcolor{cyan}{k_4} + k_5
\end{bmatrix}
\begin{bmatrix}
\Delta p_0 \\
\Delta p_1 \\
\Delta p_2 \\
\Delta p_3
\end{bmatrix}
=
\begin{bmatrix}
f_0 \\
f_1 \\
f_2 \\
f_3
\end{bmatrix}
\end{equation}
The aggregated version of equation \ref{eq:detail:K:element} is omitted here due to a poor readability.
\bigskip
The numerical version of the stiffness matrix with no forces applied as generated by the concrete example
\begin{equation}
\resizebox{.93 \textwidth}{!} {$\input{equation-nil}$}
\end{equation}
\subsection{Fixing a Point Displacement}
Determinant of matrix \textbf{K} constructed in previous chapter is equal zero.
\begin{equation}
det(\mathbf{K}) = 0
\end{equation}
It means there is no unique solution and the truss is statically unstable. Applying external force does not cause any reaction. In order to stabilise it some coordinates in some points of geometry have to be locked. It can be understood as fixing the truss to the wall. The procedure is as follows:
\begin{figure}[H]
\centering
\begin{tikzpicture}
\input{figure-fix}
\node at (1.5,.5){$p_0$};
\node at (.8,4){$p_1$};
\node at (6.1,-.2){$p_2$};
\node at (.5,-.7){$p_3$};
\draw[red,fill](4,3.5) node[right]{Fixed Point Displacement} circle[radius=.1];
\end{tikzpicture}
\caption{Fixing Tetrahedron to the Wall}
\end{figure}
Consider linear equations with exactly one solution
\begin{equation}
a_1 x_0 + b_1 x_1 + c_1 x_2 = d_1
\end{equation}
\begin{equation}
a_2 x_0 + b_2 x_1 + c_2 x_2 = d_2
\end{equation}
\begin{equation}
a_3 x_0 + b_3 x_1 + c_3 x_2 = d_3
\end{equation}
or
\begin{equation}
\begin{bmatrix}
a_1 & b_1 & c_1 \\
a_2 & b_2 & c_2 \\
a_3 & b_3 & c_3
\end{bmatrix}
\begin{bmatrix}
x_0 \\
x_1 \\
x_2
\end{bmatrix}
=
\begin{bmatrix}
d_1 \\
d_2 \\
d_3
\end{bmatrix}
\end{equation}
Fixing a point displacement requires setting selected variable to zero (no point displacement in this direction). Let's choose $x_1=0$. In order to preserve equality the constant $d_2$ must be freed. Disregarding zeros and moving new variable to the left side we obtain:
\begin{equation}
a_1 x_0 + c_1 x_2 = d_1
\end{equation}
\begin{equation}
a_2 x_0 + c_2 x_2 - d_2 = 0
\end{equation}
\begin{equation}
a_3 x_0 + c_3 x_2 = d_3
\end{equation}
\begin{equation}
\begin{bmatrix}
a_1 & 0 & c_1 \\
a_2 & -1 & c_2 \\
a_3 & 0 & c_3
\end{bmatrix}
\begin{bmatrix}
x_0 \\
d_2 \\
x_2
\end{bmatrix}
=
\begin{bmatrix}
d_1 \\
0 \\
d_3
\end{bmatrix}
\end{equation}
Concluding - in order to fix some degrees of freedom the corresponding column of the stiffness matrix has to be zeroed and a fixed variable replaced with -1. Corresponding cell in solution vector will hold \textit{freed} constant.
\begin{equation}
\input{equation-fix}
\end{equation}
\subsection{Implementation}
\lstinputlisting{populate.cpp}
\section{Solving Linear Equations by LU Decomposition Method}
\subsection{Explanation}
The method takes advantage of the fact that solving equation where matrix is triangular is straightforward. All cells above (or below) diagonal are equal zero and solution can be achieved in one iteration. The method at first breaks down the matrix into a multiplication of 2 triangular ones and then calculates sub-equations as follows.
\bigskip
Consider linear equations with exactly one solution
\begin{equation}
a_1 x_0 + b_1 x_1 + c_1 x_2 = d_1
\end{equation}
\begin{equation}
a_2 x_0 + b_2 x_1 + c_2 x_2 = d_2
\end{equation}
\begin{equation}
a_3 x_0 + b_3 x_1 + c_3 x_2 = d_3
\end{equation}
or
\begin{equation}
\begin{bmatrix}
a_1 & b_1 & c_1 \\
a_2 & b_2 & c_2 \\
a_3 & b_3 & c_3
\end{bmatrix}
\begin{bmatrix}
x_0 \\
x_1 \\
x_2
\end{bmatrix}
=
\begin{bmatrix}
d_1 \\
d_2 \\
d_3
\end{bmatrix}
\end{equation}
\begin{equation}
\mathbf{K} \Delta \mathbf{P} = \mathbf{F}
\end{equation}
LU Decomposition Method is one of the simplest, though memory hungry solver. Let's decompose matrix \textbf{A} into 2 triangular matrices LU (lower and upper).
\begin{equation}
\begin{bmatrix}
a_1 & b_1 & c_1 \\
a_2 & b_2 & c_2 \\
a_3 & b_3 & c_3
\end{bmatrix}
=
\begin{bmatrix}
1 & 0 & 0 \\
\textcolor{green}{l_{21}} & 1 & 0 \\
\textcolor{green}{l_{31}} & \textcolor{green}{l_{32}} & 1
\end{bmatrix}
\begin{bmatrix}
\textcolor{blue}{u_{11}} & \textcolor{blue}{u_{12}} & \textcolor{blue}{u_{13}} \\
0 & \textcolor{blue}{u_{22}} & \textcolor{blue}{u_{23}} \\
0 & 0 & \textcolor{blue}{u_{33}}
\end{bmatrix}
\end{equation}
The equation can be rewritten as
\begin{equation}
\mathbf{LU} \Delta \mathbf{P} = \mathbf{F}
\end{equation}
When L and U are calculated the auxiliary vector Z can be simply solved
\begin{equation}
\mathbf{LZ} = \mathbf{F}
\end{equation}
Finally the vector $\Delta \textbf{P}$ can be solved in similar way
\begin{equation}
\mathbf{U} \Delta \mathbf{P} = \mathbf{Z}
\end{equation}
\subsection{Solving an Example}
Let's solve a simple example.
\begin{equation}
\begin{bmatrix}
1 & 2 & 3 \\
4 & 14 & 19 \\
5 & 58 & 80
\end{bmatrix}
\begin{bmatrix}
x_0 \\
x_1 \\
x_2
\end{bmatrix}
=
\begin{bmatrix}
10 \\
59 \\
211
\end{bmatrix}
\end{equation}
LU Decomposition
\begin{equation}
\begin{bmatrix}
1 & 2 & 3 \\
4 & 14 & 19 \\
5 & 58 & 80
\end{bmatrix}
=
\begin{bmatrix}
1 & 0 & 0 \\
l_{21} & 1 & 0 \\
l_{31} & l_{32} & 1
\end{bmatrix}
\begin{bmatrix}
u_{11} & u_{12} & u_{13} \\
0 & u_{22} & u_{23} \\
0 & 0 & u_{33}
\end{bmatrix}
\end{equation}
It can be noticed that general direction of solving cells of matrices L and U is as depicted on figure \ref{fig:lu:direction}.
\begin{equation}
\begin{array}{lllll}
1 = 1 * u_{11} & \Rightarrow & & \Rightarrow & \textcolor{blue}{u_{11}} = 1 \\
2 = 1 * u_{12} & \Rightarrow & & \Rightarrow & \textcolor{blue}{u_{12}} = 2 \\
3 = 1 * u_{13} & \Rightarrow & & \Rightarrow & \textcolor{blue}{u_{13}} = 3 \\
4 = l_{21}u_{11} & \Rightarrow & 4 = l_{21} * 1 & \Rightarrow & \textcolor{green}{l_{21}} = 4 \\
5 = l_{31}u_{11} & \Rightarrow & 5 = l_{31} * 1 & \Rightarrow & \textcolor{green}{l_{31}} = 5 \\
14 = l_{21}u_{12} + u_{22} & \Rightarrow & 14 = 4 * 2 + u_{22} & \Rightarrow & \textcolor{blue}{u_{22}} = 6 \\
19 = l_{21}u_{13} + u_{23} & \Rightarrow & 19 = 4 * 3 + u_{23} & \Rightarrow & \textcolor{blue}{u_{23}} = 7 \\
58 = l_{31}u_{12} + l_{32}u_{22} & \Rightarrow & 58 = 5 * 2 + l_{32} * 6 & \Rightarrow & \textcolor{green}{l_{32}} = 8 \\
80 = l_{31}u_{13} + l_{32}u_{23} + u_{33} & \Rightarrow & 80 = 5 * 3 + 8 * 7 + u_{33} & \Rightarrow & \textcolor{blue}{u_{33}} = 9
\end{array}
\end{equation}
\begin{figure}[H]
\centering
\begin{tikzpicture}[scale=.4]
\draw[thick](0,-.5) -- (-.5,-.5) -- (-.5,5.5) -- (0,5.5);
\draw[->,dashed](0,5) -- (5,0);
\draw[->,dashed](0,4.5) -- (0,0);
\draw[->,dashed](2.3,2.5) -- (2.3,0);
\draw[thick](5,-.5) -- (5.5,-.5) -- (5.5,5.5) -- (5,5.5);
\draw[thick](7,-.5) -- (6.5,-.5) -- (6.5,5.5) -- (7,5.5);
\draw[->,dashed](7,5) -- (12,0);
\draw[->,dashed](7.5,5) -- (12,5);
\draw[->,dashed](9.5,2.7) -- (12,2.7);
\draw[thick](12,-.5) -- (12.5,-.5) -- (12.5,5.5) -- (12,5.5);
\end{tikzpicture}
\caption{Direction of the Decomposition in Matrices L and U}
\label{fig:lu:direction}
\end{figure}
Verification of LU Decomposition
\begin{equation}
\begin{bmatrix}
1 & 2 & 3 \\
4 & 14 & 19 \\
5 & 58 & 80
\end{bmatrix}
=
\begin{bmatrix}
1 & 0 & 0 \\
4 & 1 & 0 \\
5 & 8 & 1
\end{bmatrix}
\begin{bmatrix}
1 & 2 & 3 \\
0 & 6 & 7 \\
0 & 0 & 9
\end{bmatrix}
\end{equation}
Solving the auxiliary Z vector
\begin{equation}
\begin{bmatrix}
1 & 0 & 0 \\
4 & 1 & 0 \\
5 & 8 & 1
\end{bmatrix}
\begin{bmatrix}
z_0 \\
z_1 \\
z_2
\end{bmatrix}
=
\begin{bmatrix}
10 \\
59 \\
211
\end{bmatrix}
\end{equation}
\begin{equation}
\begin{array}{llllll}
1 * z_0 = 10 & \Rightarrow & & \Rightarrow & z_0 = 10 \\
4 * z_0 + 1 * z_1 = 59 & \Rightarrow & 4 * 10 + z_1 = 59 & \Rightarrow & z_1 = 19 \\
5 * z_0 + 8 * z_1 + 1 * z_2 = 211 & \Rightarrow & 5 * 10 + 8 * 19 + 1 * z_2 = 211 & \Rightarrow & z_2 = 9
\end{array}
\end{equation}
Verification
\begin{equation}
\begin{bmatrix}
1 & 0 & 0 \\
4 & 1 & 0 \\
5 & 8 & 1
\end{bmatrix}
\begin{bmatrix}
10 \\
19 \\
9
\end{bmatrix}
=
\begin{bmatrix}
10 \\
59 \\
211
\end{bmatrix}
\end{equation}
The final solution
\begin{equation}
\begin{bmatrix}
1 & 2 & 3 \\
0 & 6 & 7 \\
0 & 0 & 9
\end{bmatrix}
\begin{bmatrix}
x_0 \\
x_1 \\
x_2
\end{bmatrix}
=
\begin{bmatrix}
10 \\
19 \\
9
\end{bmatrix}
\end{equation}
\begin{equation}
\begin{array}{llllll}
9 * x_2 = 9 & \Rightarrow & & \Rightarrow & x_2 = 1 \\
6 * x_1 + 7 * x_2 = 19 & \Rightarrow & 6 * x_1 + 7 * 1 = 19 & \Rightarrow & x_1 = 2 \\
1 * x_0 + 2 * x_1 + 3 * x_2 = 10 & \Rightarrow & 1 * x_0 + 2 * 2 + 3 * 1 = 10 & \Rightarrow & x_0 = 3
\end{array}
\end{equation}
Verification
\begin{equation}
\begin{bmatrix}
1 & 2 & 3 \\
0 & 6 & 7 \\
0 & 0 & 9
\end{bmatrix}
\begin{bmatrix}
3 \\
2 \\
1
\end{bmatrix}
=
\begin{bmatrix}
10 \\
19 \\
9
\end{bmatrix}
\end{equation}
Final verification of the original equation
\begin{equation}
\begin{bmatrix}
1 & 2 & 3 \\
4 & 14 & 19 \\
5 & 58 & 80
\end{bmatrix}
\begin{bmatrix}
3 \\
2 \\
1
\end{bmatrix}
=
\begin{bmatrix}
10 \\
59 \\
211
\end{bmatrix}
\end{equation}
\subsection{Pitfalls of Direct Methods}
The present example uses LU Decomposition for solving matrix equations. The main downside of such approach is that data blocks are composed mostly with zeros. It is very inefficient management of memory and can be used only for simple geometries. Bigger examples may not fit into memory of modern computers.
More realistic model should be solved by one of sparse matrix methods. The most promising of them is Conjugate Gradient Method. Further work on next version of this document will be focused on it.
\subsection{The Concrete Example Result Presentation}
The below matrix equation presents the stiffness matrix, the displacement vector and external forces applied in selected points. Some cells where $-1$ is shown are used for fixed displacements and were chosen to compute reaction forces in corresponding displacement vector. The simplicity dictates such approach for price of loosing symmetry of matrix.
\bigskip
The data in the form
\begin{equation}
\mathbf{K} \Delta \mathbf{P} = \mathbf{F}
\end{equation}
is generated below as a result of computation of the main example
\begin{equation}
\input{equation-full}
\end{equation}
\subsection{Implementation of LU Decomposition Method}
\lstinputlisting{calculate.cpp}
\section{The 3D Data Presentation in the Documentation}
This chapter tries to explain how to present three-dimensional data in the document. Using only $x,y$ and removing $z$ coordinate is not a good solution since we loose part of the information and presented figures are not clear. Rotation is much better as all of original coordinates can be used and after transformation the new coordinate $z'$ can be discarded without loosing general understanding.
\subsection{Point Rotation}
Let's consider rotation of a point in coordinates x, y
\begin{figure}[H]
\centering
\begin{tikzpicture}[scale=.5]
% global arrows
\draw[<->](0,10) node[right] {$y$} -- (0,0) -- (10,0) node[below] {$x$};
% auxyliary lines
\draw[dashed](8,0) -- (8,4) -- (0,4);
\draw[dashed](4,0) -- (4,8) -- (0,8);
% angles
\draw(2,0) arc(0:25:2);
\node at (10:2.6){$\alpha_0$};
\draw[->](3,1.5) arc(30:70:3);
\node at (45:2.7){$\alpha_1$};
% original and rotated lines
\draw[thick,red](0,0) -- node[above left]{$r$} (8,4) node[right]{$(x_0,y_0)$};
\draw[fill](8,4) circle[radius=.1];
\draw[thick,red,dashed](0,0) -- node[above left]{$r$} (4,8) node[right]{$(x_1,y_1)$};
\draw[fill](4,8) circle[radius=.1];
\end{tikzpicture}
\caption{Point Rotation}
\label{fig:point:rotation}
\end{figure}
From geometry on the figure \ref{fig:point:rotation} we can obtain
\begin{equation}
x_0 = r \cos \alpha_0
\end{equation}
\begin{equation}
y_0 = r \sin \alpha_0
\end{equation}
\begin{equation}
x_1 = r \cos (\alpha_0 + \alpha_1) = r \cos \alpha_0 \cos \alpha_1 - r \sin \alpha_0 \sin \alpha_1
\end{equation}
\begin{equation}
y_1 = r \sin (\alpha_0 + \alpha_1) = r \sin \alpha_0 \cos \alpha_1 + r \cos \alpha_0 \sin \alpha_1
\end{equation}
hence
\begin{equation}
x_1 = x_0 \cos \alpha_1 - y_0 \sin \alpha_1
\end{equation}
\begin{equation}
y_1 = x_0 \sin \alpha_1 + y_0 \cos \alpha_1
\end{equation}
Changes on the axis z can be expressed simply by $z_1 = z_0$. Constructing three-dimensional matrix we obtain rotation operator around the axis z
\begin{equation}
\begin{bmatrix}
x_1 \\
y_1 \\
z_1
\end{bmatrix}
=
\begin{bmatrix}
\cos \alpha_1 & - \sin \alpha_1 & 0 \\
\sin \alpha_1 & \cos \alpha_1 & 0 \\
0 & 0 & 1
\end{bmatrix}
\begin{bmatrix}
x_0 \\
y_0 \\
z_0
\end{bmatrix}
\end{equation}
\begin{equation}
\mathbf{P_1} = \mathbf{R_z P_0}
\end{equation}
\subsection{An Example of Rotated Tetrahedron}
The example rotates all points of geometry around the axis y and then x. It can be expressed as
\begin{equation}
\mathbf{P_1} = \mathbf{R_x R_y P_0}
\end{equation}
In details
\begin{equation}
\begin{bmatrix}
x_1 \\
y_1 \\
z_1
\end{bmatrix}
=
\begin{bmatrix}
1 & 0 & 0 \\
0 & c_x & s_x \\
0 & -s_x & c_x
\end{bmatrix}
\begin{bmatrix}
c_y & 0 & s_y \\
0 & 1 & 0 \\
-s_y & 0 & c_y
\end{bmatrix}
\end{equation}
\begin{equation}
\begin{bmatrix}
x_1 \\
y_1 \\
z_1
\end{bmatrix}
=
\begin{bmatrix}
c_y & 0 & s_y \\
-s_xs_y & c_x & s_xc_y \\
-c_xs_y & -s_x & c_xc_y
\end{bmatrix}
\begin{bmatrix}
x_0 \\
y_0 \\
z_0
\end{bmatrix}
\end{equation}
Finally the coordinate $z_1$ is cut off as it does not have representation in 2D. Visualized results are presented on the figure \ref{fig:object:rotation}.
\begin{figure}[H]
\centering
\begin{tikzpicture}
\input{figure-rotate}
\node at (-6,-.4){$p_0,p_3$};
\node at (-6.3,4){$p_1$};
\node at (-1,-.4){$p_2$};
\node at (1.5,.5){$p_0'$};
\node at (.8,4){$p_1'$};
\node at (6.1,-.2){$p_2'$};
\node at (.6,-.7){$p_3'$};
\node at (-4,4){a) Original};
\node at (3,4){b) Rotated};
\draw[<->](-2,3.5) node[left] {$y$} -- (-2,2.5) node[below] {$z=0$} -- (-1,2.5) node[below] {$x$};
\end{tikzpicture}
\caption{Original and Rotated Tetrahedron Example}
\label{fig:object:rotation}
\end{figure}
\subsection{Implementation of the Rotation Procedure}
\lstinputlisting{rotate.cpp}
\section{Summary}
Finally the article described all practical aspects of the finite method presentation. Starting with geometry definition, theoretical method explanation, solving matrix equations and 3D data presentation. Reader could confront mathematical aspects with practical implementation for better understanding.
\end{document}