Today, FEM is probably the most common method for solving a wide range of applied engineering problems. Historically, he emerged from mechanics, but was subsequently applied to all sorts of non-mechanical problems.

Today there is a wide variety of software packages, such as ANSYS , Abaqus , Patran , Cosmos , etc. These software packages allow you to solve the problems of structural mechanics, fluid mechanics, thermodynamics, electrodynamics and many others. The implementation of the method itself, as a rule, is considered quite complex and voluminous.

Here I want to show that at present, using modern tools, writing the simplest FEM calculator from scratch for a two-dimensional plane-stress state problem is not something very complicated and cumbersome. I chose this kind of task because it was the first successful example of applying the finite element method. And of course it is the easiest to implement. I am going to use a linear, three-nodal element, since this is the only flat element in which case numerical integration is not required, as will be shown below. For elements of a higher order, with the exception of the integration operation (which is not entirely trivial, but its implementation is quite interesting), the idea is absolutely the same.

Picture to attract attention:


Historically, the first practical application of FEM was in the field of mechanics, which significantly affected the terminology and the first interpretations of the method. In the simplest case, the essence of the method can be described as follows: the continuous medium is replaced by an equivalent hinge system, and the solution of statically indeterminate systems is a well-known and studied problem of mechanics. This simplified, engineering interpretation contributed to the widespread use of the method, however, strictly speaking, such an understanding of the method is incorrect. The exact mathematical justification of the method was given much later than the first successful applications of the method, but this allowed us to expand the scope for a wider range of problems, not only from the field of mechanics.

I am not going to describe the method in detail, there is a lot of literature about it, instead I will focus on the implementation of the method. Requires minimal knowledge of mechanics, or the same sopromat. I will also be glad, reviews of people not related to mechanics, whether at least the idea was clear. The method will be implemented in C ++, however, I will not use any very specific features of this language and I hope that the code will be understandable to people who do not know this language.

Since this is just an example of the implementation of the method, so as not to complicate and leave everything in the most simple form for understanding, I will prefer brevity to the detriment of many programming principles. This is not an example of writing good, maintainable, and reliable code; it is an example of implementing FEM. So there will be many simplifications to concentrate on the main goal.

Input data

Before starting the method itself, let's find out in what form we will present the input data. The object in question should be divided into a large number of small elements, in our case, triangles. Thus, we replace the continuous medium with a set of nodes and finite elements forming a grid. The figure below shows an example of a grid with numbered nodes and elements.



In the figure we see nine nodes and eight elements. To describe the grid, you need two lists:

• A list of nodes that defines the coordinates of each node.
• List of items.

In the list of elements, each element is determined by the set of indices of the nodes that the element forms. In our case, we only have triangular elements, so we will use only three indexes for each element.
As an example, for the grid shown above, we will have the following lists:

List of nodes:

0.000000    31.500000
15.516667   31.500000
0.000000    19.489759
18.804134   23.248226
0.000000    0.000000
20.479981   11.365822
27.516667   19.500000
27.516667   11.365822
27.516667   0.000000

List of items:

1 3 2
2 3 4
4 6 7
4 3 6
7 6 8
3 5 6
6 5 9
6 9 8

It is worth noting that there are several ways to specify the same element. Node indices can be listed clockwise or counterclockwise. The counterclockwise enumeration is usually used, so we will assume that all elements are specified in this way.

Let's create some input file with a full description of the task. To avoid confusion and not complicate, it is better to use indexing, which starts from zero, and not from one, since in C / C ++ arrays are indexed from scratch. For the first test input file, we will create the simplest grid possible:



Let the first line be a description of the material properties. For example, for steel, Poisson's ratio ν = 0.3 and Young's modulus E = 2000 MPa:

0.3 2000

Then follows a line with the number of nodes and then the list of nodes itself:

4
0.0 0.0
1.0 0.0
0.0 1.0
1.0 1.0

Then a line with the number of elements follows, then a list of elements:

2
0 1 2
1 3 2

Now, we need to set the binding, or as they say, the conditions on the border of the first kind. As such boundary conditions, we will fix the movements of the nodes. You can fix movements along the axes independently of each other, i.e. we can prohibit movement along the x-axis or the y-axis, or both at once. In the general case, you can specify a certain initial displacement, however, we restrict ourselves to only a zero initial displacement. Thus, we will have a list of nodes with a given type of restriction on movement. Type of restriction will be indicated by number:

• 1 - movement in the x axis direction is prohibited
• 2 - movement in the direction of the y axis is prohibited
• 3 - it is forbidden to move in both x and y directions

The first line defines the number of restrictions:

2
0 3
1 2

Then, we must set the load. We will work only with nodal forces. Strictly speaking, nodal forces should not be regarded as forces in the general sense of the word, I will talk about this later, but for now, let's think of them simply as forces acting on a node. We need a list with node indices and corresponding force vectors. The first line determines the amount of applied force:

2
2 0.0 1.0
3 0.0 1.0

It can be easily seen that the problem under consideration is a square-shaped body, the bottom of which is fixed, and tensile forces act on the upper face.



Entire input file:

0.3 2000
4
0.0 0.0
1.0 0.0
0.0 1.0
1.0 1.0
2
0 1 2
1 3 2
2
0 3
1 2
2
2 0.0 1.0
3 0.0 1.0

Eigen - math library

Before starting to write code, I would like to talk about the mathematical library - Eigen . It is a mature, reliable and efficient library that has a very clean and expressive API. It has many compile-time optimizations, the library is able to perform explicit vectorization and supports SSE 2/3/4, ARM and NEON instructions. As for me, this library is remarkable if only because it allows you to implement complex mathematical calculations in a concise and expressive way.
I would like to make a brief description of the part of the Eigen API that we will use later. Those readers who are familiar with this library may skip this section.

There are two types of matrices: dense and sparse. We will use both types. In the case of a dense matrix, all elements are in memory, in the order of the indices, both column-major (default) and raw-major placement are supported. This type of matrix is ​​good for relatively small matrices or matrices with a small number of zero elements. Sparse matrices are good for storing large matrices with a small number of nonzero elements. We will use a sparse matrix for the global stiffness matrix.

Dense Matrices

To use dense matrices we will need to include a header . There are two main types of dense matrices: fixed and dynamic. A fixed-size matrix is ​​a matrix whose size is known at compile time. In the case of a matrix of dynamic size, its size can be set at the stage of code execution; moreover, the size of the dynamic matrix can even be changed on the fly. Of course, you need to give preference to matrices with a fixed size wherever possible. Memory for dynamic matrices is allocated on the heap, and an unknown matrix size limits the optimizations that the compiler can do. Fixed matrices can be allocated on the stack, the compiler can deploy loops and much more. It is worth noting that a mixed type is also possible when the matrix has a fixed number of rows, but the number of columns is dynamic, or vice versa.

Fixed size matrices can be defined as follows:

Eigen::Matrix m;

Where m is a float matrix with a fixed size of 3 × 3. You can also use the following predefined types:

Eigen::Matrix2i
Eigen::Matrix3i
Eigen::Matrix4i
Eigen::Matrix2f
Eigen::Matrix3f
Eigen::Matrix4f
Eigen::Matrix2d
Eigen::Matrix3d
Eigen::Matrix4d

There are a few more predefined types; these are basic. The number indicates the dimension of the square matrix and the letter indicates the type of matrix element: i is an integer, f is a single-precision floating-point number, d is a double-precision floating-point number.

For vectors there is no separate type, vectors are the same matrices. The column vector (sometimes row vectors are found in the literature, we have to refine it), can be declared as follows:

Eigen::Matrix v;

Or you can use one of the predefined types (the notation is the same as for matrices):

Eigen::Vector2i
Eigen::Vector3i
Eigen::Vector4i
Eigen::Vector2f
Eigen::Vector3f
Eigen::Vector4f
Eigen::Vector2d
Eigen::Vector3d
Eigen::Vector4d

For quick reference, I wrote the following example:

#include 
int main()
{
    //Fixed size 3x3 matrix of floats
    Eigen::Matrix3f A;
    A   <<  1, 0, 1,
            2, 5, 0,
            1, 1, 2;
    std::cout   << "Matrix A:" 
                << std::endl
                << A 
                << std::endl;
    //Access to matrix element [1, 1]:
    std::cout   << "A[1, 1]:"
                << std::endl
                << A(1, 1) 
                << std::endl;
    //Fixed size 3 vector of floats
    Eigen::Vector3f B;
    B << 1, 0, 1;
    std::cout   << "Vector B:"
                << std::endl
                << B
                << std::endl;
    //Access to vector element [1]:
    std::cout   << "B[1]:"
                << std::endl
                << B(1)
                << std::endl;
    //Multiplication of vector and matrix
    Eigen::Vector3f C = A * B;
    std::cout   << "C = A * B:"
                << std::endl
                << C
                << std::endl;
    //Dynamic size matrix of floats
    Eigen::MatrixXf D;
    D.resize(3, 3);
    //Let matrix D be an identity matrix:
    D.setIdentity();
    std::cout   << "Matrix D:" 
                << std::endl
                << D
                << std::endl;
    //Multiplication of matrix and matrix
    Eigen::MatrixXf E = A * D;
    std::cout   << "E = A * D:" 
                << std::endl
                << E
                << std::endl;
    return 0;
}

Conclusion:

Matrix A:
1 0 1
2 5 0
1 1 2
A[1, 1]:
5
Vector B:
1
0
1
B[1]:
0
C = A * B:
2
2
3
Matrix D:
1 0 0
0 1 0
0 0 1
E = A * D:
1 0 1
2 5 0
1 1 2

For more detailed information, it is better to refer to the documentation .

Sparse Matrices

Sparse matrix is ​​a slightly more complex type of matrix. The basic idea is not to store the entire matrix in memory as it is, but to store only nonzero elements. Quite often in practical problems there are large matrices with a small number of nonzero elements. When assembling the global stiffness matrix in the FEM model, the number of nonzero elements can be less than 0.1% of the total number of elements.

For modern FEM packages, it is not a problem to solve a problem with hundreds of thousands of nodes, moreover, on a quite ordinary modern PC. If we try to allocate a place to store the global stiffness matrix, we will encounter the following problem:



The memory size required to store the matrix is ​​huge! A sparse matrix will require 10,000 times less memory.

Sparse matrices use memory more efficiently, since they store only nonzero elements. Sparse matrices can be represented in two ways: compressed and uncompressed. In an uncompressed form, it is easy to add or remove an element from the matrix, but this type of representation is not optimal in terms of memory usage. Eigen, when working with a sparse matrix, uses a kind of compressed representation, slightly more optimized for the dynamic addition / removal of elements. Eigen also knows how to convert such a matrix to a compressed form, moreover, it does it transparently, you do not need to do this explicitly. Most algorithms require a compressed form of the matrix, and the use of any of these algorithms automatically converts the matrix into a compressed form. And vice versa,

How can we define a matrix? A good way to do this is to use the so-called Triplets . This is a data structure, and this is a template class, which is one matrix element in combination with indexes that determine its position in the matrix. We can specify a sparse matrix directly from the triplets vector .

For example, we have the following sparse matrix:

0   3   0   0   0
22  0   0   0   0
0   5   0   1   0
0   0   0   0   0
0   0   14  0   8

The first thing we should do is include the appropriate Eigen library header: . Then, we declare an empty sparse 5x5 matrix. Next, we define the triplets vector and fill it with values. The triplet constructor accepts the following arguments: i-index, j-index, value.

#include 
#include 
int main()
{
    Eigen::SparseMatrix A(5, 5);
    std::vector< Eigen::Triplet > triplets;
    triplets.push_back(Eigen::Triplet(0, 1, 3));
    triplets.push_back(Eigen::Triplet(1, 0, 22));
    triplets.push_back(Eigen::Triplet(2, 1, 5));
    triplets.push_back(Eigen::Triplet(2, 3, 1));
    triplets.push_back(Eigen::Triplet(4, 2, 14));
    triplets.push_back(Eigen::Triplet(4, 4, 8));
    A.setFromTriplets(triplets.begin(), triplets.end());
    A.insert(0, 0);
    std::cout << A;
    A.makeCompressed();
    std::cout << std::endl << A;
}

The conclusion will be as follows:

Nonzero entries:
(0,0) (22,1) (_,_) (3,0) (5,2) (_,_) (_,_) (14,4) (_,_) (_,_) (1,2) (_,_) (_,_) (8,4) (_,_) (_,_)
Outer pointers:
0 3 7 10 13  $
Inner non zeros:
2 2 1 1 1  $
0 3 0 0 0
22 0 0 0 0
0 5 0 1 0
0 0 0 0 0
0 0 14 0 8
Nonzero entries:
(0,0) (22,1) (3,0) (5,2) (14,4) (1,2) (8,4)
Outer pointers:
0 2 4 5 6  $
0 3 0 0 0
22 0 0 0 0
0 5 0 1 0
0 0 0 0 0
0 0 14 0 8

Data structures

To store the data that we are going to read from the input file, as well as intermediate data, we must declare our structures. The simplest data structure describing the finite element is shown below. It consists of an array of three elements that store the indices of the nodes forming the finite element, as well as the matrix B. This matrix is ​​called the gradient matrix and we will return to it later. The CalculateStiffnessMatrix method will also be described later.

struct Element
{
    void CalculateStiffnessMatrix(const Eigen::Matrix3f& D, std::vector >& triplets);
    Eigen::Matrix B;
    int nodesIds[3];
};

Another simple structure that we will need is a structure for describing bindings. This is a simple structure consisting of an enumerated type that defines the type of restriction, and the index of the node on which the restriction is imposed.

struct Constraint
{
    enum Type
    {
        UX = 1 << 0,
        UY = 1 << 1,
        UXY = UX | UY
    };
    int node;
    Type type;
};

To keep things simple, let's define some global variables. Creating any global objects is always a bad idea, but for this example it’s completely acceptable. We will need the following global variables:

• Number of nodes
• Vector with x-coordinate nodes
• Vector with y-coordinate of nodes
• Elements Vector
• Pins Vector
• Load Vector

In the code, we define them as follows:

int                         nodesCount;
int                         noadsCount;
Eigen::VectorXf             nodesX;
Eigen::VectorXf             nodesY;
Eigen::VectorXf             loads;
std::vector< Element >      elements;
std::vector< Constraint >   constraints;

Reading input file

Before we read something, we must know where to read and where to write the output. At the beginning of the main function, we check the number of input arguments, note that the first argument is always the path to the executable file. Thus, we should have three arguments, let the second be the path to the input file, and the third the path to the file with the output. For work with input / output, for a specific case it is most convenient to use file streams from the standard library.

The first thing we do is create threads for I / O:

int main(int argc, char *argv[])
{
    if ( argc != 3 )
    {
        std::cout << "usage: " << argv[0] << " 
\n";
        return 1;
    }
    std::ifstream infile(argv[1]);
    std::ofstream outfile(argv[2]);

The first line of the input file contains material properties, read them:

float poissonRatio, youngModulus;
infile >> poissonRatio >> youngModulus;

This is enough to build an elastic matrix of an isotropic material for a plane-stressed state, which is defined as follows:





Where did these expressions come from? They are easily obtained from Hooke's law in general, in fact, we can find an expression for the D matrix from the following relationships:



It is worth noting that the plane-stress state means that σ _{the Z} is equal to zero, but not e _{the Z} . A plane-stressed model is well suited for solving a number of engineering problems, where plane structures are considered and all forces act in a plane. At a time when calculating volumetric problems was too expensive, many tasks were reduced to flat, sacrificing accuracy.

In a plane-stressed state, nothing prevents the body from deforming in the direction normal to its plane, so the strain ε _Z is generally not equal to zero. It does not appear in the equations above, but can be easily obtained from the following relation, given that the stress σ _Z is zero:


We have all the initial data to set the elasticity matrix, let's do this:

Eigen::Matrix3f D;
D <<
    1.0f,           poissonRatio,   0.0f,
    poissonRatio,   1.0,            0.0f,
    0.0f,           0.0f,           (1.0f - poissonRatio) / 2.0f;
D *= youngModulus / (1.0f - pow(poissonRatio, 2.0f));

Next, we read the list with the coordinates of the nodes. First we read the number of nodes, then we set the size of the dynamic vectors x and y. Next, we just read the coordinates of the nodes in the loop, line by line.

infile >> nodesCount;
nodesX.resize(nodesCount);
nodesY.resize(nodesCount);
for (int i = 0; i < nodesCount; ++i)
{
    infile >> nodesX[i] >> nodesY[i];
}

Then, we read the list of items. Everything is the same, we read the number of elements, and then the node indices for each element:

int elementCount;
infile >> elementCount;
for (int i = 0; i < elementCount; ++i)
{
    Element element;
    infile >> element.nodesIds[0] >> element.nodesIds[1] >> element.nodesIds[2];
    elements.push_back(element);
}

Next, read the list of pins. All the same:

int constraintCount;
infile >> constraintCount;
for (int i = 0; i < constraintCount; ++i)
{
    Constraint constraint;
    int type;
    infile >> constraint.node >> type;
    constraint.type = static_cast(type);
    constraints.push_back(constraint);
}

You might notice static_cast , it is needed to convert an int type into a previously enumerated type for bindings that we declared.

Finally, you need to read the list of loads. There is one peculiarity associated with the load, because of which we will present the acting loads in the form of a dimension vector of a double number of nodes. The reason why we do so will be explained a bit later. The figure below illustrates this load vector:



Thus, for each node, we have two elements in the load vector that represent the x and y components of the force. Thus, the x-component of the force in a particular node will be stored in the element with the index id = 2 * node_id + 0, and the y-component of the force for this node will be stored in the element with the index id = 2 * node_id + 1.

First, set the size of the vector the applied effort is twice as large as the number of nodes, then zero the vector. We read the number of external forces and read their values ​​line by line:

int loadsCount;
loads.resize(2 * nodesCount);
loads.setZero();
infile >> loadsCount;
for (int i = 0; i < loadsCount; ++i)
{
    int node;
    float x, y;
    infile >> node >> x >> y;
    loads[2 * node + 0] = x;
    loads[2 * node + 1] = y;
}

Calculation of the global stiffness matrix

We consider a geometrically linear system with infinitesimal displacements. Moreover, we assume that the deformation occurs elastically, i.e. deformations are a linear function of stresses (Hooke's law). From the foundations of building mechanics, it can be shown that the movement of each node is a linear function of the applied forces. Thus, we can say that we have the following system of linear equations:



Where: K - stiffness matrix; δ is the displacement vector; R is the vector of loads, that is, the vector of external forces. This system of linear equations is often called an ensemble, since it represents a superposition of the stiffness matrices of each element, as will be shown below.

The displacement vector for two-dimensional problems can be defined as follows:





Where: u _i and v_i is the u-component and the v-component of the displacement of the i-th node.

Vector of external forces:





Where: R _xi and R _yi are the x-component and the y-component of the external force applied to the i-th node.

As we see, each element of the displacement vector and the load vector is a two-dimensional vector. Instead, we can define these vectors as follows:



That is, in fact, the same, but it simplifies the representation of these vectors in the code. This is an explanation of why we previously set the load vector in such a peculiar way.

How to build a stiffness matrix? The essence of the global stiffness matrix is ​​the superposition of the stiffness matrices of each element. If we consider one element separately from the rest, then we can determine the same stiffness matrix, but only for this element. This matrix determines the relationship between the movements of nodes and nodal forces. For example, for a 3-node element:





Where: [k] ^e is the stiffness matrix of the e-element; [δ] ^e is the vector of displacements of the nodes of the e-element; [F] ^e is the vector of nodal forces of the e-element; i, j, m are the indices of the nodes of the element.

What happens if one node is split between two elements? First of all, since this is still the same node, the movement, respectively, will not depend on the composition of which element we are considering this node. The second important consequence is that if we summarize all the nodal forces from each element in this node, the result will be equal to the external force, in other words, the load.

The sum of all the nodal forces in each node is equal to the sum of the external forces in this node, this follows from the principle of equilibrium. Despite the fact that this statement may seem quite logical and fair, in fact, everything is a little more complicated. Such a formulation is not accurate, because the nodal forces are not forces in the general sense of the word, and the fulfillment of the equilibrium conditions is not an obvious condition for them. Despite the fact that the statement is still true, it uses a mechanical principle, which limits the application of the method. There is a more rigorous, mathematical justification for the FEM from the standpoint of minimizing potential energy, which makes it possible to extend the method to a wider range of problems.

In my opinion, it is better to think about nodal forces as some kind of generalized forces, in the representation of Lagrange mechanics. The fact is that moving a node is actually a generalized moving, since it uniquely characterizes the field of displacements inside an element. Each component of the node’s movement can be interpreted as a generalized coordinate, therefore, the forces act not at the node’s point, but along the element’s boundary (or in volume for the volume forces), and they act precisely on the movement that characterizes the given node.

If you sum up all the equations for each node, the inter-element nodal forces will go away. On the right side of the equation, only external forces remain - loads. Thus, given this fact, we can write: The



figure below illustrates the above expression:



To build a global stiffness matrix, we need a triplets vector . In the loop, we will go through each element and fill this vector with the values ​​of stiffness matrices obtained from each element:

std::vector > triplets;
for (std::vector::iterator it = elements.begin(); it != elements.end(); ++it)
{
    it->CalculateStiffnessMatrix(D, triplets);
}

Where std :: vector > - vector triplets ; CalculateStiffnessMatrix is an element class method that fills the triplets vector .

As mentioned earlier, we can build a global sparse matrix directly from the triplets vector :

Eigen::SparseMatrix globalK(2 * nodesCount, 2 * nodesCount);
globalK.setFromTriplets(triplets.begin(), triplets.end());

Calculation of the stiffness matrix of individual elements

We figured out how to assemble a global stiffness matrix (ensemble) from matrices of elements. It remains to find out how to get the stiffness matrix of an element.

Form functions

Let's start by interpolating the movements within the element based on the movements of its nodes. If the displacements of the nodes are specified, then the displacement at any point of the element can be represented as:



Where [N] is a matrix of coordinate functions (x, y) . These functions are called form functions. Each displacement component, u and v, can be interpolated independently, then the equation above can be rewritten in the following form:



Or it can be written in a separate form:




How to obtain an inerpolating function? For simplicity, let's look at a scalar field. In the case of a three-node linear element, the interpolation function is linear. We will look for the interpolation expression in the following form:



If the values ​​of fin the nodes are known, then we can specify a system of three equations:



Or in matrix form:



From this system of equations, we can find an unknown vector of coefficients for the interpolating expression:



We introduce the notation



Finally, we get the interpolation expression:



Returning to the displacements, we can say that :





Then, the form functions will look like:

Deformations and gradient matrix

Knowing the displacement field, we can find the deformation field, according to the relations from the theory of elasticity:



Now we can replace u and v with function of the form:






Or we can write the same in combined form:



The matrix on the right side, usually denoted as matrix B , the so-called gradient matrix.



To find the matrix B , we must find all the partial derivatives of the form functions:





In the case of a linear element, we see that the partial derivatives of the form functions are actually constants, which will save us a lot of effort, as will be shown below. Multiplying the row vector with constants by the inverse matrix C we get:





Now we can calculate the B matrix. To construct the matrix C, we first define the coordinate vectors of the nodes:

Eigen::Vector3f x, y;
x << nodesX[nodesIds[0]], nodesX[nodesIds[1]], nodesX[nodesIds[2]];
y << nodesY[nodesIds[0]], nodesY[nodesIds[1]], nodesY[nodesIds[2]];

Then, we can build matrix C from three rows:

Eigen::Matrix3f C;
C << Eigen::Vector3f(1.0f, 1.0f, 1.0f), x, y;

Then we calculate the inverse matrix C and collect the matrix B:

Eigen::Matrix3f IC = C.inverse();
for (int i = 0; i < 3; i++)
{
	B(0, 2 * i + 0) = IC(1, i);
	B(0, 2 * i + 1) = 0.0f;
	B(1, 2 * i + 0) = 0.0f;
	B(1, 2 * i + 1) = IC(2, i);
	B(2, 2 * i + 0) = IC(2, i);
	B(2, 2 * i + 1) = IC(1, i);
}

Transition to stresses

As we mentioned earlier, the relationship between stress and strain can be written using the elasticity matrix D:



Thus, substituting the expression above into the relationship for strains, we get:



Here we see that stresses, like strains, are constants inside the element. This is a feature of linear elements. However, for elements of higher orders, continuity of stresses is also not satisfied. With an increase in the number of elements, these gaps should decrease, which is a criterion of convergence. In practice, the average strain value for each node is usually calculated, however, in our case, we will leave everything as it is.

The principle of virtual work

In order to move to the stiffness matrix, we must turn to virtual work. Let's say we have some virtual displacements in the nodes of the element: d [δ] ^e

Virtual displacements should be understood as any infinitely small displacements that can occur. We know that in the case of our task, there is only one solution and only one set of movements will take place, but here we look at everything from a slightly different angle. In other words, we take all infinitesimal displacements that are kinematically allowed, and then we find the only ones that satisfy the physical law.

For these virtual movements, the virtual work of the nodal forces will be:



Virtual work of internal stresses per unit volume for the same virtual movements:



Or:



Substituting the equation for the stresses, we get:



From the law of conservation of energy, the virtual work of external forces should be equal to the sum of the work of internal stresses over the volume of the element:



Since this ratio should be true for any virtual displacement, we can divide both sides of the equation into virtual displacement :



Nodal displacements can be taken out of the integral, since they are constant inside the element. Now we see that the integral is a proportionality coefficient between nodal displacements and nodal forces, which means that the integral is a stiffness matrix:



For a linear element, as we received earlier, the matrix B is a constant. If the material properties are also constant, then the integral becomes trivial:



Where A is the area of ​​the element, t is the thickness of the element. From the properties of a triangle, its area can be obtained as half the determinant of the matrix C:



In the end, we can write the following code to calculate the stiffness matrix:

Eigen::Matrix K = B.transpose() * D * B * C.determinant() / 2.0f;

Here I lowered the thickness, we assume that it is single with us. Usually, in practice, use the following system of units: MPa, mm ² . In this case, if we set the elastic modulus in megapascals, and the dimensions in millimeters, then the thickness will be one millimeter. If a different thickness is needed, then it can easily be introduced into the elastic modulus. In a more general case, it is best to set the thickness elementwise or even for each node.

The entire CalculateStiffnessMatrix method :

void Element::CalculateStiffnessMatrix(const Eigen::Matrix3f& D, std::vector >& triplets)
{
	Eigen::Vector3f x, y;
	x << nodesX[nodesIds[0]], nodesX[nodesIds[1]], nodesX[nodesIds[2]];
	y << nodesY[nodesIds[0]], nodesY[nodesIds[1]], nodesY[nodesIds[2]];
	Eigen::Matrix3f C;
	C << Eigen::Vector3f(1.0f, 1.0f, 1.0f), x, y;
	Eigen::Matrix3f IC = C.inverse();
	for (int i = 0; i < 3; i++)
	{
		B(0, 2 * i + 0) = IC(1, i);
		B(0, 2 * i + 1) = 0.0f;
		B(1, 2 * i + 0) = 0.0f;
		B(1, 2 * i + 1) = IC(2, i);
		B(2, 2 * i + 0) = IC(2, i);
		B(2, 2 * i + 1) = IC(1, i);
	}
	Eigen::Matrix K = B.transpose() * D * B * C.determinant() / 2.0f;
	for (int i = 0; i < 3; i++)
	{
		for (int j = 0; j < 3; j++)
		{
			Eigen::Triplet trplt11(2 * nodesIds[i] + 0, 2 * nodesIds[j] + 0, K(2 * i + 0, 2 * j + 0));
			Eigen::Triplet trplt12(2 * nodesIds[i] + 0, 2 * nodesIds[j] + 1, K(2 * i + 0, 2 * j + 1));
			Eigen::Triplet trplt21(2 * nodesIds[i] + 1, 2 * nodesIds[j] + 0, K(2 * i + 1, 2 * j + 0));
			Eigen::Triplet trplt22(2 * nodesIds[i] + 1, 2 * nodesIds[j] + 1, K(2 * i + 1, 2 * j + 1));
			triplets.push_back(trplt11);
			triplets.push_back(trplt12);
			triplets.push_back(trplt21);
			triplets.push_back(trplt22);
		}
	}
}

At the end of the method, after calculating the stiffness matrix, the Triplets vector is filled . Each Triplet is created using an array of element node indexes to determine its position in the global stiffness matrix. Just in case, I emphasize that we have two sets of indices, one local for the element matrix, the other global for the ensemble.

Joining Jobs

The resulting system of linear equations cannot be solved until we fix. Purely from a mechanical point of view, an unsecured system can perform any large displacements, and under the influence of external forces it will completely go into motion. Mathematically, this will lead to the fact that the obtained SLAE is not linearly independent, and therefore does not have a unique solution.

The movements of some nodes must be set to zero or to some predetermined value. In the simplest case, let's only consider fixing nodes, without initial offsets. Actually, this means that we are removing some equations from the system. But changing the amount of system equations on the fly is not a trivial task. Instead, you can use the following trick.

Say we have the following system of equations:



To fix a node, the corresponding matrix element must be set to 1, and all elements in the row and column containing this element must be set to zero. There should also be no external forces acting on the fixed unit in the direction in which the restriction acts. An equation with this node will obviously give a zero offset for this node, and zeros in the corresponding column will exclude this offset from other equations. Let's say we want to fix the zero node in the x-axis direction, and fix the first node in the y-axis direction:



To perform this operation, we first determine the indices of the equations that we are going to exclude from SLAE. To do this, we will go through the list of pins in a cycle and collect the index vector. Then, sorting through all nonzero elements of the stiffness matrix, we call the SetConstraints function . The following is the ApplyConstraints function :

void ApplyConstraints(Eigen::SparseMatrix& K, const std::vector& constraints)
{
	std::vector indicesToConstraint;
	for (std::vector::const_iterator it = constraints.begin(); it != constraints.end(); ++it)
	{
		if (it->type & Constraint::UX)
		{
			indicesToConstraint.push_back(2 * it->node + 0);
		}
		if (it->type & Constraint::UY)
		{
			indicesToConstraint.push_back(2 * it->node + 1);
		}
	}
	for (int k = 0; k < K.outerSize(); ++k)
	{
		for (Eigen::SparseMatrix::InnerIterator it(K, k); it; ++it)
		{
			for (std::vector::iterator idit = indicesToConstraint.begin(); idit != indicesToConstraint.end(); ++idit)
			{
				SetConstraints(it, *idit);
			}
		}
	}
}

The SetConstraints function checks whether the stiffness matrix element is part of the excluded equation, and if so, it sets it to zero or one, depending on whether the element is on the diagonal or not:

void SetConstraints(Eigen::SparseMatrix::InnerIterator& it, int index)
{
	if (it.row() == index || it.col() == index)
	{
		it.valueRef() = it.row() == it.col() ? 1.0f : 0.0f;
	}
}

Next, we simply call ApplyConstraints and pass the global stiffness matrix and the binding vector as arguments to it:

ApplyConstraints(globalK, constraints);

SLAU solution and output

The Eigen library has various sparse linear equation solvers available, we will use SimplicialLDLT , which is a fast direct solver. For demonstration purposes, we will derive the original stiffness matrix and the load vector, and then output the result of the solution, the displacement vector.

std::cout << "Global stiffness matrix:\n";
std::cout << static_cast >& > (globalK) << std::endl;
std::cout << "Loads vector:" << std::endl << loads << std::endl << std::endl;
Eigen::SimplicialLDLT > solver(globalK);
Eigen::VectorXf displacements = solver.solve(loads);
std::cout << "Displacements vector:" << std::endl << displacements << std::endl;
outfile << displacements << std::endl;

Of course, the matrix inference makes sense only for our test case, which contains only two elements.

Analysis of the result

As we saw earlier, we can move from displacements to strains and further to stresses:



However, we obtain a stress vector, which in elasticity theory generally speaking is a symmetric tensor of the second rank:



As you know, the tensor elements directly depend on the coordinate system, however the physical state itself is of course independent. Therefore, for analysis, it is better to go over as some kind of invariant quantity, it is best to scalar. Mises stresses are such an invariant :



This scalar quantity allows us to estimate stresses in a complex stress state, as if we were dealing with a uniaxial stress state. This criterion is not ideal, but in most cases, for static analysis it gives a more or less plausible picture, therefore it is most often used in practice.

Now we will go through all the elements in the cycle, and calculate the Mises stresses for each element:

std::cout << "Stresses:" << std::endl;
for (std::vector::iterator it = elements.begin(); it != elements.end(); ++it)
{
	Eigen::Matrix delta;
	delta << displacements.segment<2>(2 * it->nodesIds[0]),
	         displacements.segment<2>(2 * it->nodesIds[1]),
	         displacements.segment<2>(2 * it->nodesIds[2]);
	Eigen::Vector3f sigma = D * it->B * delta;
	float sigma_mises = sqrt(sigma[0] * sigma[0] - sigma[0] * sigma[1] + sigma[1] * sigma[1] + 3.0f * sigma[2] * sigma[2]);
	std::cout << sigma_mises << std::endl;
	outfile << sigma_mises << std::endl;
}

On this, the writing of our simplest FEM calculator can be considered complete. Let's see what conclusion we get for our test task:

Global stiffness matrix:
1 0 0 0 0 0 0 0
0 1 0 0 0 0 0 0
0 0 1483.52 0 0 714.286 -384.615 -384.615
0 0 0 1 0 0 0 0
0 0 0 0 1483.52 0 -1098.9 -329.67
0 0 714.286 0 0 1483.52 -384.615 -384.615
0 0 -384.615 0 -1098.9 -384.615 1483.52 714.286
0 0 -384.615 0 -329.67 -384.615 714.286 1483.52
Loads vector:
0
0
0
0
0
1
0
1
Deformations vector:
            0
            0
      -0.0003
            0
-5.27106e-011
        0.001
      -0.0003
        0.001
Stresses:
2
2

We see that we set the fixations correctly, the deformations of the fixed nodes in the corresponding directions are equal to zero. The deformations of the upper nodes are directed upward, in accordance with the tensile force. Deformations in the x-direction for the right pair of nodes are 0.0003, which is 0.3 part of the deformations of the upper nodes in the y-direction, which is the Poisson effect . The calculation result is fully consistent with expectations, so you can proceed to the tasks more interesting.



Run the cloc utility :

-------------------------------------------------------------------------------
Language                     files          blank        comment           code
-------------------------------------------------------------------------------
C++                              1             37              0            171
-------------------------------------------------------------------------------

As you can see, nothing at all - 171 lines of code.

Some pictures

To visualize the calculation results, a python script was written that draws an undeformed and deformed mesh with a stress field. Unfortunately, I don’t know how to do gradient fill with PIL , so the stresses inside the element are displayed as a solid constant fill. Abaqus can do visualization in such a way that, although it looks less beautiful, it’s closer to the truth.

Test task:



To get more complicated grids, a free student version of Abaqus was used . Another script was written to convert Abaqus ' input file . True, fixing and external loads will still have to be set manually.

Here we see a strip with a hole, the bottom of which is fixed, and a tensile force is applied to the upper edge:


Abstract design element, just for beauty. The bottom of the structure is fixed, and a concentrated force is applied to the upper protruding part, which I drew already on top of the picture:


Another example with a strip with a hole, this time only a quarter of the structure is considered due to symmetry. The left border is fixed along the x axis, the bottom - along the y axis. Maximum Received Voltage: 3.05152. This value, despite the coarseness of the grid (and maybe due to it), is quite close to the theoretical value - 3.


The same problem, but solved in Abaqus. As you can see, the maximum voltage is 3.052, which completely coincides with our result. Such an exact coincidence is, in principle, not surprising, since in practice it is difficult for a triangular element to do something in its own way, so that the result is different. For elements of a higher order, unfortunately, this 100% coincidence did not work out, which can be explained by a different implementation of numerical integration.


All source code, including examples, is available on github .

What remains behind the scenes

Honestly, behind the scenes a lot is left. But the article was already inflated, so I think you can stop there.
Not included (but wanted):

• Задание внешних нагрузок. Как было сказано, в исходном файле мы задаем узловые силы. Но к ним еще нужно перейти от внешних нагрузок. В случае сконцентрированной силы тут все тривиально, но в случае приложения давления (как в 2-м и 4-м примере), здесь все не так просто. В случае элементов более высоких порядков, результат даже немного неожиданный (с точки зрения инженерной интуиции).
• Не коснулся я важных свойств функций формы.
• Хотелось поговорить о элементах более высоких порядков, на примере изопараметрических элементов.
• Хотелось так-же показать, как выполнить численное интегрирование для элементов более высоких порядков. Это отдельная, весьма любопытная тема.
• С минимальными изменениями в коде, можно выполнить динамический анализ, когда учитываются силы инерции и система приходит в движение.
• It was possible, along with the consideration of the equilibrium of nodal forces, to give a more general justification of the method, excluding these nodal forces in general.

I think that many can rightly notice that in modern mathematical packages, such as MATLAB, there are tools for working with FEM, and in principle there is no need to write anything in C ++. However, here I still wanted to show a minimal implementation, with a minimum of dependencies, i.e. without using third-party tools. With the exception of the mathematical library, of course, but I think that all the functions used for readers should be extremely transparent.

PS The article is big, I'm sure there are a lot of typos, please in PM.