(c) Juan Gomez 2019. Thanks to Universidad EAFIT for support. This material is part of the course Introduction to Finite Element Analysis

# A simple finite element code
## Preliminary-Discrete systems

The strategy behind a finite element algorithm is that of using the concept of discrete elements to render a continous system into a discrete problem. A continous system is governed by a set of partial differential equations and boundary conditions defining a boundary value problem, while a discrete problem is governed by a set of linear (or non-linear) algebraic equations. As it turns out, the idea behind a finite element algorithm is to combine mathematical principles and numerical methods to render the continuous BVP into a discrete system of equations. The following example of a set of masses connected by springs is a useful introductory problem to study finite element methods as the system is discrete in nature allowing us to follow most of the algorithmic aspects of a finite element code without the complex parafernalia of numerical methods or mathematical principles.

Here the problem consists of an assemblage of masses joined by different springs submitted to static loads. The springs will play the role of finite elements, while the masses would represent special points called nodes within the finite element jargon. One such a system is shown in the figure.

<center></center>


### Equilibrium equations for a typical spring (element).

Consider a typical spring (finite element) of stiffness coefficient $k$  like the one shown in the figure

<center></center>

and with displacements and forces at each end denoted by $u_i$ and $f_i$ respectively and where $i = 1,2$ depending on the specific end. If the spring is subject to different displacements of the nodal points one obtains the following forces written in matrix form like:

$$
    \begin{Bmatrix}
        f_1\\
        f_2
    \end{Bmatrix} =
    K\begin{bmatrix}
          1.0 & -1.0\\
        - 1.0 & 1.0
    \end{bmatrix}
    \begin{Bmatrix}
        u_1\\
        u_2
    \end{Bmatrix}
$$


### Equilibrium equations for a typical mass.

Consider now the equilibrium equation:


$$
f_2^i + f_1^{i + 1} + m_j\frac{dV_j}{dt} = P_j.
$$

for a typical mass $m_j$ as the one shown in the figure


<center></center>

Note that as indicated by the free body diagram, the mass is connected to springs $i$ and $i+1$ and is also under the action of an external load $P$. At the same time the dashed vector respresents the inertial load $m_j\frac{dV_j}{dt}$.



Letting the displacement of the $m_j$ mass be $u_j$ and expressing spring forces in terms of displacements in the equlibrium equation gives:

$$
(K^i+K^{i+1})u_j-K^iu_{j-1}-K^{i+1}u_{j+1}+m_j\frac{dV_j}{dt}=P_j
$$


The "finite element" equations for the complete spring-mass system written in general matrix form like:

$$
\left[ {{K_G}} \right]\left\{ {{U_G}} \right\} + \left[ M \right]\left\{ {{A_G}} \right\} = \left\{ {{F_G}} \right\}.
$$

is obtained after considering the equilibrium relations for each mass.

### Computer implementation.

To write the equilibrium equations in a systematic fashion suitable for a general finite element code consider the following 3-mass system.

<center></center>


Writing the equilibrium equations for the springs $i$ and $i+1$ in terms of displacements $u_{j - 1}$, $u_j$ and $u_{j + 1}$ we have:

$$
\left\{ {\begin{array}{*{20}{c}}
{f_1^i}\\
{f_2^i}
\end{array}} \right\} = \left[ {\begin{array}{*{20}{c}}
{k_{11}^i}&{k_{12}^i}\\
{k_{21}^i}&{k_{22}^i}
\end{array}} \right]\left\{ {\begin{array}{*{20}{c}}
{{u_{j - 1}}}\\
{{u_j}}
\end{array}} \right\}
$$

and

$$
\left\{ {\begin{array}{*{20}{c}}
{f_1^{i + 1}}\\
{f_2^{i + 1}}
\end{array}} \right\} = \left[ {\begin{array}{*{20}{c}}
{k_{11}^{i + 1}}&{k_{12}^{i + 1}}\\
{k_{21}^{i + 1}}&{k_{22}^{i + 1}}
\end{array}} \right]\left\{ {\begin{array}{*{20}{c}}
{{u_j}}\\
{{u_{j + 1}}}
\end{array}} \right\}
$$

Npte that we have used a row-column index notation for the stiffness coefficients in order to facilitate the computer implementation. The equations for the $m_j$ mass then reads:

$$
k_{21}^i{u_{j - 1}} + (k_{22}^i + k_{11}^{i + 1}){u_j} + k_{12}^{i + 1}{u_{j + 1}} + m_j\frac{dV_j}{dt} = {P_j}.
$$


Considering also the contributions from the springs $K^i$ and $K^{i+1}$ to the equilibrium of masses $m_{j-1}$ and $m_{j+1}$ respectively we have the following block from the complete system of equations.



$$
\left[ {\begin{array}{*{20}{c}}
{}&{}&{}&{}\\
{}&{k_{11}^i}&{k_{12}^i}&{}\\
{}&{k_{21}^i}&{k_{22}^i + k_{11}^{i + 1}}&{k_{12}^{i + 1}}\\
{}&{}&{k_{21}^{i + 1}}&{k_{22}^{i + 1}}
\end{array}} \right]
$$


Considering now the complete system of masses and springs leads to a system of linear equations of the form
$$
\left[ {{K_G}} \right]\left\{ {{U_G}} \right\} + \left[ M \right]\left\{ {{A_G}} \right\} = \left\{ {{F_G}} \right\}.
$$

where each equation represents the equilibrium of a given mass.

#### Assemblage.

The construction of the global matrices governing the equilibrium of each mass in the system may be achieved in a very systematic way after adding up the contribution from each spring to the global matrix. This process is called assembly of the global equilibrium equations. This assembly operation can be performed after establishing the connection between the global and local degrees of freedom. This can be done through an operator storing in each row the identifiers for the global degrees of freedom corresponding to each element. For instance, in the 3-mass system the springs (or elementts) $i$ and $i+1$ have end displacements $j-1$ and $j$ and $j$ and $j+1$ respectively. These indices provide all the required information to conduct the assembly process. The matrix storing the global indices for all the elements in the model is called here the **DME()** operator and given like:


$$
DME = \left[ {\begin{array}{*{20}{c}}
{}&{}\\
{j - 1}&j\\
j&{j + 1}\\
{}&{}
\end{array}} \right]
$$


With the **DME()** operator available the assembly proceeds as indicated next:


$$\begin{array}{l}
{K_{j - 1,j - 1}} \leftarrow {K_{j - 1,j - 1}} + k_{11}^i\\
{K_{j - 1,j}} \leftarrow {K_{j - 1,j}} + k_{12}^i\\
{K_{j,j - 1}} \leftarrow {K_{j,j - 1}} + k_{21}^i\\
{K_{j,j}} \leftarrow {K_{j,j}} + k_{22}^i
\end{array}
$$

and

$$
\begin{array}{l}
{K_{j,j}} \leftarrow {K_{j,j}} + k_{11}^{i + 1}\\
{K_{j,j + 1}} \leftarrow {K_{j,j + 1}} + k_{12}^{i + 1}\\
{K_{j + 1,j}} \leftarrow {K_{j + 1,j}} + k_{21}^{i + 1}\\
{K_{j + 1,j + 1}} \leftarrow {K_{j + 1,j + 1}} + k_{22}^{i + 1}
\end{array}
$$

Notice the connection between the local indices, here corresponding to $1$ and $2$ and the possitions in the global matrix, here corresponding to $j-1$, $j$ and $j+1$.

### Example


```python
%matplotlib inline        
import matplotlib.pyplot as plt
import numpy as np
import sympy as sym
```

Consider the following system

<center></center>

The required input files containing the input data for the masses (nodes), springs (elements), stiffness coefficients (materials) and loads are in the folder **files** of this REPO. The input files are read in the following piece of code:


```python
def readin():
    nodes = np.loadtxt('files/' + 'sprnodes.txt', ndmin=2)
    mats = np.loadtxt('files/' + 'sprmater.txt', ndmin=2)
    elements = np.loadtxt('files/' + 'spreles.txt', ndmin=2, dtype=np.int)
    loads = np.loadtxt('files/' + 'sprloads.txt', ndmin=2)

    return nodes, mats, elements, loads
```

In the nodes file, storing the information from each mass, there is a $-1$ or a $0$ value indicateing if a given mass is restrained or free. Such data allows the code to assign an equation number to each free mass as done in the subroutine **eqcounter**


```python
def eqcounter(nodes):

    nn = nodes.shape[0]
    IBC = np.zeros([nn, 1], dtype=np.integer)
    neq = 0
    for i in range(nn):
        IBC[i] = int(nodes[i, 2])
        if IBC[i] == 0:
            IBC[i] = neq
            neq = neq + 1

    return neq, IBC
```

The equation number assigned to each mass is now used to create the **DME()** operator. Note that each row contains the identifiers for the end displacements in the current spring.


```python
def DME(nodes, elements):

    nels = elements.shape[0]
    DME = np.zeros([nels, 2], dtype=np.integer)

    neq, IBC = eqcounter(nodes)
    ndof = 2
    nnodes = 2
    for i in range(nels):
        for j in range(nnodes):
            kk = elements[i, j+3]
            DME[i, j] = IBC[kk]
    return DME, IBC, neq
```

Using the **DME()** operator it is now possible to assemble the global matrix of stiffness coefficients in terms of equations of the type:

$$\begin{array}{l}
{K_{j - 1,j - 1}} \leftarrow {K_{j - 1,j - 1}} + k_{11}^i\\
\end{array}
$$



```python
def assembly(elements, mats, nodes, neq, DME, uel=None):

    IELCON = np.zeros([2], dtype=np.integer)
    KG = np.zeros((neq, neq))
    nels = elements.shape[0]
    nnodes = 2
    ndof = 2
    for el in range(nels):
        elcoor = np.zeros([nnodes])
        im = np.int(elements[el, 2])
        par0 = mats[im]
        for j in range(nnodes):
            IELCON[j] = elements[el, j+3]
            elcoor[j] = nodes[IELCON[j], 1]
        kloc = uelspring(par0)
        dme = DME[el, :ndof]
        for row in range(ndof):
            glob_row = dme[row]
            if glob_row != -1:
                for col in range(ndof):
                    glob_col = dme[col]
                    if glob_col != -1:
                        KG[glob_row, glob_col] = KG[glob_row, glob_col] +\
                            kloc[row, col]

    return KG
```

**Question: What is the function of the following subroutines?**


```python
def uelspring(kcof):
    """1D-2-noded Spring element

    Kcof : float
      Stiffness coefficient (>0).

    Returns
    -------
    kl : ndarray
      Local stiffness matrix for the element (2, 2).


    """
    kl = np.zeros([2, 2])
    kl[0, 0] = kcof
    kl[0, 1] = -kcof
    kl[1, 0] = -kcof
    kl[1, 1] = kcof

    return kl
```


```python
def loadasem(loads, IBC, neq, nl):
    """Assembles the global Right Hand Side Vector RHSG

    Parameters
    ----------
    loads : ndarray
      Array with the loads imposed in the system.
    IBC : ndarray (int)
      Array that maps the nodes with number of equations.
    neq : int
      Number of equations in the system after removing the nodes
      with imposed displacements.
    nl : int
      Number of loads.

    Returns
    -------
    RHSG : ndarray
      Array with the right hand side vector.

    """
    RHSG = np.zeros([neq])
    for i in range(nl):
        il = int(loads[i, 0])
        ilx = IBC[il]
        if ilx != -1:
            RHSG[ilx] = loads[i, 1]

    return RHSG
```

The main program is then compossed of the following steps:
* Read the model
* Build the DME() operator
* Assembly the global system of equations
* Solve for the global displacements $UG$


```python
nodes, mats, elements, loads = readin()
DME, IBC, neq = DME(nodes, elements)
KG = assembly(elements, mats, nodes, neq, DME)
RHSG = loadasem(loads, IBC, neq, 3)
UG = np.linalg.solve(KG, RHSG)
print(UG)
```

    [0.002  0.0025 0.0045]


### References

* Bathe, Klaus-Jürgen. (2006) Finite element procedures. Klaus-Jurgen Bathe. Prentice Hall International.

* Juan Gómez, Nicolás Guarín-Zapata (2018). SolidsPy: 2D-Finite Element Analysis with Python, <https://github.com/AppliedMechanics-EAFIT/SolidsPy>.


```python
from IPython.core.display import HTML
def css_styling():
    styles = open('./nb_style.css', 'r').read()
    return HTML(styles)
css_styling()
```





<link href='http://fonts.googleapis.com/css?family=Fenix' rel='stylesheet' type='text/css'>
<link href='http://fonts.googleapis.com/css?family=Alegreya+Sans:100,300,400,500,700,800,900,100italic,300italic,400italic,500italic,700italic,800italic,900italic' rel='stylesheet' type='text/css'>
<link href='http://fonts.googleapis.com/css?family=Source+Code+Pro:300,400' rel='stylesheet' type='text/css'>

<style>

/*
Template for Notebooks for ModelaciÃ³n computacional.

Based on Lorena Barba template available at:

    https://github.com/barbagroup/AeroPython/blob/master/styles/custom.css
*/

/* Fonts */
@font-face {
font-family: "Computer Modern";
src: url('http://mirrors.ctan.org/fonts/cm-unicode/fonts/otf/cmunss.otf');
}

/* Text */
div.cell{
width:800px;
margin-left:16% !important;
margin-right:auto;
}
h1 {
font-family: 'Alegreya Sans', sans-serif;
}
h2 {
font-family: 'Fenix', serif;
}
h3{
font-family: 'Fenix', serif;
margin-top:12px;
margin-bottom: 3px;
}
h4{
font-family: 'Fenix', serif;
}
h5 {
font-family: 'Alegreya Sans', sans-serif;
}	
div.text_cell_render{
font-family: 'Alegreya Sans',Computer Modern, "Helvetica Neue", Arial, Helvetica, Geneva, sans-serif;
line-height: 135%;
font-size: 120%;
width:600px;
margin-left:auto;
margin-right:auto;
}
.CodeMirror{
font-family: "Source Code Pro";
font-size: 90%;
}
/* .prompt{
display: None;
}*/
.text_cell_render h1 {
font-weight: 200;
font-size: 50pt;
line-height: 100%;
color:#CD2305;
margin-bottom: 0.5em;
margin-top: 0.5em;
display: block;
}	
.text_cell_render h5 {
font-weight: 300;
font-size: 16pt;
color: #CD2305;
font-style: italic;
margin-bottom: .5em;
margin-top: 0.5em;
display: block;
}
.warning{
color: rgb( 240, 20, 20 )
}
</style>









```python

```

03bc82ec5301f55d1dc5aa0faac63713e09bd024
ipynb
Jupyter Notebook
notebooks/00b_springs.ipynb
AppliedMechanics-EAFIT/Introductory-Finite-Elements
a4b44d8bf29bcd40185e51ee036f38102f9c6a72
MIT
2019-11-26T13:28:30.000Z
2022-02-16T17:57:11.000Z
notebooks/00b_springs.ipynb
jgomezc1/Introductory-Finite-Elements
a4b44d8bf29bcd40185e51ee036f38102f9c6a72
MIT
notebooks/00b_springs.ipynb
jgomezc1/Introductory-Finite-Elements
a4b44d8bf29bcd40185e51ee036f38102f9c6a72
MIT
2020-02-17T07:24:59.000Z
2022-03-02T07:54:28.000Z
Qwen/Qwen-72B
1. YES
2. YES