How Magnetite functions.
Magnetite covers meshing, solving, and post-processing—all in one program. Let's break it down.
If you haven't already, checkout the main readme; this details how you can submit geometry into Magnetite.
The goal of meshing is to turn our 2D model into Nodes and Elements. Nodes are effectively just points (with some extra metadata), and elements are triangles made from three Nodes. Thats not bad, is it?
Magnetite supports both CSV inputs and SVG inputs. In both cases, the desired end result is a vector of vectors of Nodes. We need a vector of vectors because we need a way to differentiate between external geometry and internal geometry. Because we know that there's only ever going to be one external geometry, we can reserve the
Once we've established the nodes in the model, we can translate them into a .geo file. This will serve as the input to Gmsh. There is a Gmsh SDK for Rust, but this method had the benefit of debugging the .geo file at will. At some point, I may switch to this better method.
With the .geo file in hand, we can run gmsh geom.geo -2 -o geom.msh to get our geom.msh. Simple enough!
Next, we open up the geom.msh file and start parsing it. There are a few different sections in a .msh file, but we are concerned with the $Elements section and the $Nodes section. First, we'll register the nodes, then the elements. There's some tricks we follow to ignore elements from other dimensions, and, thanks to the .msh file format being so developer-friendly, this is no issue.
Once we have established the nodes and elements in our model, we're ready to apply the boundary conditions described in the input json. We'll loop though the nodes and check its position against all of the registered boundary rules. If the number boundary conditions stay constant, the time complexity remains at
Once we've finally parsed all of the nodes and elements in our model, we're ready to solve. But wait, how do we solve something like this?
The math in this solver is based heavily on this University of New Mexico Paper, plus some other resources. For a comprehensive list of citations, view the bottom of the main readme.
My personal favorite law in physics is Hooke's law. It says:
Ok, so what's the big deal? Spring's alone aren't that interesting, but what's cool is that everything is a spring. We're taught that the reason why we don't fall through the floor is due to "the normal force" — and this is true! But what is the normal force? It's just Hooke's law! When you walk down the road, you compress it, just enough for it to respond with the force to match your weight. As with everything, there's cavitates to this, but this is tue for linear-elastic systems.
So, if everything is a spring, why don't we represent our model as one big spring? Well, if we can, we do! In engineering statics, we say
But, what happens when the geometry gets too complicated? You'd be hard-pressed to derive an equation for your three-story building using simple stress-strain relationships. This is where the finite-element method (FEM) comes into play. If we break down complex geometries that we are unable to solve into thousands of elements that we can solve, we are able to create a system of equations that governs the reveals the behavior of the system as a whole. This looks like:
Here,
Before we get too deep into the math, let's go back to Hooke's law again.
Consider what happens when we compress the string over a distance
If we were to consider
... where
Cool, but what's the significance? If we know the potential energy for an element in our system, we can leverage the The Principle of Minimum Potential Energy (MPE) to arrive at a stiffness matrix. This principle states:
"For conservative structural systems, of all the kinematically admissible deformations, those corresponding to the equilibrium state extremize (i.e., minimize or maximize) the total potential energy. If the extremum is a minimum, the equilibrium state is stable."
In other terms, the deformations that minimize the potential energy are the equilibrium deformations; because we're building a static mechanical simulator, we are interested in the static (i.e., stable) solution, which is found at the MPE.
Note: when referring to potential energy, we are referring to total potential energy including work potential—i.e., the opposite of work done on the system by external systems. For this reason, potential energy, for a 1D spring system, is represented as
$KE=KE_{sp}+WP$ — where$WP$ is work potential.
For instance, consider the following system of springs (from Indian Institute of Science):
This image is missing
$F_3$ ; it is applied to the lowest arrow attached to the [3] block.
We are interested in solving for the equilibrium displacements for
We need to also consider the work done on the system by external forces. Work is force times distance (assuming force is constant), so we can say that the two forces do the combined work on this system:
Remember, work potential is the opposite of the work done on the system by external forces. Therefore, the total potential energy is:
As stated in the principle of MPE, the potential energy for each
We now have a system we're able to solve. This derivation is for a 1D spring system, but we can follow the same steps to arrive at a stiffness matrix for planar elements.
- Obtain an equation potential energy
- Apply the Principle of MPE
- Factor out displacements to create a stiffness matrix
- Setup the equation in the form
${F}=[K]{u}$
Before we get started working with 2D triangular elements, we'll need to consider some new properties at arise with the increased dimension. Let's look at 1D Hooke's law again, in it's stress-strain form:
where:
-
$\sigma$ is stress -
$\epsilon$ is strain -
$E$ is the material elasticity (Young's Modulus)
Let's see how this changes when we move to 2D.
First, let's look at a stress element in 2D:
The symbols represent:
-
$\sigma$ – normal stress -
$\tau$ – shear stress
We'll also need to look at the poisson ratio (
When the green block is stretched, it becomes thinner in its other axes; the magnitude of this effect is described by poisson ratio.
Now, we have the following definitions for the properties that define a 2D stress element:
- Axial Stress (
$\sigma$ ) - Axial Strain (
$\epsilon$ ) - Shear Stress (
$\tau$ ) - Shear Strain (
$\gamma$ ) (also not shown, described here) - Young's Modulus (
$E$ ) - Poisson Ratio (
$v$ )
2D square elements are simple enough, but we're dealing with triangles—that makes things a bit more tricky. Let's examine how a triangle behaves when it's stressed.
Triangles have three vertices, which we call nodes. Each node can deform in the
Now, let's consider what happens when the triangle deforms. There's a few ways to approach this problem, but the simplest way is to interpolate stress over the area of the triangle. This approach sacrifices accuracy, but its a reasonable approximation for small displacements. To apply this theory, we split the element into three sub-triangles with a common vertex at the interior point of the element:
From University of New Mexico
Considering that the area of the element is the sum of the areas
Because
With these definitions, we can estimate the displacements in the element with:
And, we can also solve for the position of interior point:
This will allow us to do some cool stuff in a moment—hang in there.
Finding the potential energy in a spring element is pretty simple; however, the process becomes a lot more complex for 2D elements. We use the same idea of integrating the force over displacement, as we did for the spring, for the plane as well. Instead of looking the force, we'll look at the stress
Remember, we're working in 2D, so
We want our model to be defined in two dimensions. If we assume a constant thickness
Our goal is to define stress and strain so that we can apply the Principle of MPE to create a stiffness matrix for our problem. Once we have the stiffness matrix, we solve the system:
where ...
-
${F}$ is a column vector of nodal forces -
$[K]$ is a square, sparse stiffness matrix -
${U}$ is a column vector of nodal displacements
For small stresses in an isotropic material (i.e., a material whose properties are consistent regardless of direction), the 2D version of Hooke's law states:
The derivations of these equations are very complex. Feel free to read more about the derivation here
Let's use the equations above to write stress as a matrix, in terms of strain. This will be helpful for solving our big equation in a moment. If we solve the following equation from above for
... we get:
This can be substituted into another equation from above ...
Now, it is possible to solve for
Lastly, we can rearrange ...
... to solve for
Phew! That's a lot. Let's clean this up into one nice vector:
This allows us to simply state that
Let's look back at our problem:
We've been able to write an expression for stress in terms of strain, so we can re-write this equation as:
We're close, but we now need an expression for strain. We know that strain can be found from displacement:
Subscript e referes to the element as a whole, not a node.
Recall the equations we derived for a triangular element:
We have expressions for
This can be expressed in the following system:
We can flip this equation to get expressions for
Now, we can take the partial derivatives of the shape functions we derived for triangular elements. This yields the following:
Substituting these values into the matrix above, we can say:
We're almost there! Let's get rid of the inverse on the second matrix. Notice that this matrix is actually the Jacobian of the multivariable displacement function. For this reason, we'll let...
Now, our equation becomes:
There's one more trick we can apply here. It turns out that the area of a triangle can be defined as
If we multiply to solve for
Again, recall how we derived these equations:
We can take the partial derivatives of these to solve for
Substituting this into the expression before yields:
Great! Now we do everything again for
Finally, our original strain matrix becomes:
Which, (thankfully) simplifies down to:
Where
$x_{ij} = (x_i-x_j)$ and$y_{ij} = (y_i-y_j)$ for simplicity
We'll label these two matrices as...
...such that
Finally, we have an expression for the strain matrix that we can substitute back into our potential energy equation:
The equation we have created for potential energy is a double-integral over some region A. However, all of the stuff in the integral is constant; therefore, we can pull it out of the integral:
Finally, we can say that:
Notice how similar this expression is to
Finally, the last step is to sum the element stiffness matrices into a total stiffness matrix that represents the entire system; I'll talk about how this is built in the following chapter
There's a lot of math in the last step that would make it appear as if a lot is going on, but there honestly isn't; the solver needs to take the following steps:
- Determine the element stiffness matrix of each element
- Compile these element stiffness matrices into a total stiffness matrix
- Set up the finite-element equation
- Solve for displacements and forces.
All the hard work for this step was discussed in the previous section. The equations I referenced were initially derived from people much smarter than myself. In the code, all we have to do is plug-and-chug.
The element stiffness matrix, for each element, is given as:
Thus, we need to find
Where
$x_{ij} = (x_i-x_j)$ and$y_{ij} = (y_i-y_j)$ for simplicity
Furthermore, we assume that thickness
Using nalgebra, we can create a stiffness matrix for each element.
Each element stiffness matrix is square. The rows and columns correspond to pairs of nodal displacements:
Notice how there is a stiffness between every two nodes?
The example stiffness matrix above would correspond to the stiffness matrix for an element with nodes 1, 2, and 3 because it defines a stiffness at each of these nodes.
The total stiffness matrix will look similar, just with more nodes. For each element stiffness matrix, we match each row-column node pair from the local stiffness matrix to the total stiffness matrix. This is easier explained by showing the code:
for (stiffness_mat, element) in std::iter::zip(element_stiffness_matrices, elements) {
for (local_row, node_row) in element.nodes.iter().enumerate() {
for (local_col, node_col) in element.nodes.iter().enumerate() {
let global_row = node_row * 2;
let global_col = node_col * 2;
let local_row = local_row * 2;
let local_col = local_col * 2;
// Add RowX ColX
total_stiffness_matrix[(global_row, global_col)] +=
stiffness_mat[(local_row, local_col)];
// Add RowX ColY
total_stiffness_matrix[(global_row, global_col + 1)] +=
stiffness_mat[(local_row, local_col + 1)];
// Add RowY ColX
total_stiffness_matrix[(global_row + 1, global_col)] +=
stiffness_mat[(local_row + 1, local_col)];
// Add RowY ColY
total_stiffness_matrix[(global_row + 1, global_col + 1)] +=
stiffness_mat[(local_row + 1, local_col + 1)];
}
}
}We iterate through each element element and it's corresponding local stiffness matrix stiffness_mat. Then, we iterate over the array of nodes that belong to the element; we do this twice to get every pair. Next, we calculate global_row and global_col; this will be location of the stiffness for the two nodes in the total stiffness matrix. We do the same thing for local_row and local_col to find the location of the stiffness in the local stiffness matrix.
In the derivation section, you might have realized how the vector
Notice how each node appears twice? It's
Most of the work has been done by the mesher by parsing the mesh into elements and nodes with boundary conditions. We set up our equation to be:
However, we will almost always have unknown forces and unknown displacements. This makes things a bit tricker to solve. It turns out that, in a properly constrained model, whenever a force is known, the displacement is unknown; the same is true vice-versa. We can leverage this artifact by initially ignoring all unknown forces. For the purpose of this example, imagine
Now, all the forces on the left are known, but there's still a mixture of known and unknown displacements on the right hand side. When multiplying across, the first three columns of the stiffness matrix will always correspond to unknown displacements, and the last three columns will always correspond to known displacements. Therefore, if we partition the matrix depending on the number of unknown displacements, then subtract the known partition from the nodal forces, we can solve the matrix for the unknown displacements.
Once we have found the displacements, we rewrite the equation again—this time, leaving out known forces:
Because all the nodal displacements are now known, we can quickly multiply across to solve for the nodal forces.
Once we've solved for nodal displacements, it's easy to go back and solve for stresses. In the derivation section, we found that stress can be represented by strain:
... and that strain can be represented by
Because we know the nodal displacements of each element, it's easy to plug the solved values into this equation to solve for stress.
A full citation list can be found at the bottom of the main readme