## Finite Elements

The finite element method (FEM) or analysis (FEA) is a *numerical method* for solving many engineering domain problems. It's areas include structural analysis, heat transfer and fluid flow around an object. 

The method yields approximate values of the unknowns using a discrete number of points over the domain of the object in question. To solve the problem, it integrates by parts, dividing a large problem into smaller, simpler objects that are called *Finite Elements*. The simple equations that model the finite elements are then assembled into larger systems of equations to model the entire problem. 

### 1D problem

The temperature $T$ of a bar of length 1 at equilibrium satisfies:

$$
\partial_{xx}T + f(x)=0
$$

The heat source is given by $f(x)$. We'll hold the temperature fixed at the right edge by setting $T(1)=0$. We'll allow heat to conduct through the left edge by setting $\partial_xT(0)=0$.

We need a way of handling the second-derivative term, so we can introduce a weighting function $w(x)$ and integrate by parts to get:

$$
  \left[ \partial_x T(x) w(x) \right]_0^1 - \int_0^1 \text{d}x \, \partial_x T(x) \partial_x w(x) = - \int_0^1 \text{d}x \, f(x) w(x).
$$

and after factoring the boundary conditions we get:

$$
\int_0^1 \text{d}x\partial_xT(x)\partial_xw(x)=\int_0^1 \text{d}xw(x)f(x)
$$

This is considered the **weak form** of the problem.

### Step 1. A simple 3-node example

To get an approximate solution, we split the domain into intervals/elements. Let's begin with a simple example where we choose $[0, 1/2]$ and $[1/2, 1]$, which means there are nodes at $\{0, 1/2, 1\}$. To repesent this, we'll introduce some *shape functions* $N_A(x)$ and write:

$$
T(x)=\sum_A T_AN_A(x)
$$

here $A$ is a counter associated with the nodes, so $A=\{0,1,2\}$: the nodes are $x_A$. The shape function $N_A(x)=1$ at the associated node, and $0$ at all other nodes. In this case, let's choose the shape functions to be piecewise linear, which are shape functions look as:

$$
\begin{align}
  N_0(x) &= \begin{cases} 1 - 2 x, & 0 \le x \le 1/2 \\ 0, & 1/2 \le x \le 1, \end{cases}, \\
  N_1(x) &= \begin{cases} 2 x, & 0 \le x \le 1/2 \\ 2 - 2 x, & 1/2 \le x \le 1, \end{cases}, \\
  N_2(x) &= \begin{cases} 0, & 0 \le x \le 1/2 \\ 2 x - 1, & 1/2 \le x \le 1. \end{cases}
\end{align}
$$


```python
import matplotlib.pyplot as plt
%matplotlib inline

x = np.linspace(0.,1.,1000)
# we use np.where to select x <= 0.5, then replace true with top function, false with bottom function
N0 = np.where(x <= .5, 1-2*x, np.zeros_like(x))
N1 = np.where(x <= .5, 2*x, 2-2*x)
N2 = np.where(x <= .5, np.zeros_like(x), 2*x-1)
plt.plot(x, N0, label='$N_0$')
plt.plot(x, N1, label='$N_1$')
plt.plot(x, N2, label='$N_2$')
plt.legend()
plt.show()
```

We can now write out the weak form of this problem as:

$$
  \sum_A \sum_B T_A w_B \int_0^1 \text{d}x \, \partial_x N_A(x) \partial_x N_B(x) = \sum_B w_B \int_0^1 \text{d}x \, N_A(x) f(x).
$$

This has to be true for *any* choice of weight function, so for *any* coefficients $w_B$. This is true only if

$$
  \sum_A T_A \int_0^1 \text{d}x \, \partial_x N_A(x) \partial_x N_B(x) = \int_0^1 \text{d}x \, N_B(x) f(x).
$$

which can be written as a linear system matrix problem:

$$
KT=F
$$

The coefficients of the **stiffness matrix** $K$ are given as:

$$
  K_{AB} = \int_0^1 \text{d}x \, \partial_x N_A(x) \partial_x N_B(x).
$$

The coefficients of the **force vector** $F$ are given by

$$
  F_{A} = \int_0^1 \text{d}x \, N_A(x) f(x).
$$

In order to calculate the **stiffness** matrix and **force** vectors, we need to calculate the derivatives of the shape functions:

$$
\begin{align}
  \partial_x N_0(x) &= \begin{cases} - 2, & 0 \le x \le 1/2 \\ 0, & 1/2 \le x \le 1, \end{cases}, \\
  \partial_x N_1(x) &= \begin{cases} 2, & 0 \le x \le 1/2 \\ - 2, & 1/2 \le x \le 1, \end{cases}, \\
  \partial_x N_2(x) &= \begin{cases} 0, & 0 \le x \le 1/2 \\ 2, & 1/2 \le x \le 1. \end{cases}
\end{align}
$$

Therefore the product $\partial_xN_A\partial_xN_B$ required for calculating the stiffness matrix is either $0$ or $\pm 4$:

$$
\begin{align}
  \partial_x N_0(x) \partial_x N_0(x) &= \begin{cases} 4, & 0 \le x \le 1/2 \\ 0, & 1/2 \le x \le 1, \end{cases}, \\
  \partial_x N_0(x) \partial_x N_1(x) &= \begin{cases} -4, & 0 \le x \le 1/2 \\ 0, & 1/2 \le x \le 1, \end{cases}, \\
  \partial_x N_1(x) \partial_x N_0(x) &= \begin{cases} -4, & 0 \le x \le 1/2 \\ 0, & 1/2 \le x \le 1, \end{cases}, \\
  \partial_x N_1(x) \partial_x N_1(x) &= \begin{cases} 0, & 0 \le x \le 1/2 \\ 4, & 1/2 \le x \le 1. \end{cases}
\end{align}
$$

This gives:

$$
K_{AB}= \begin{pmatrix} 2 & -2 \\ -2 & 4 \end{pmatrix}
$$

Meanwhile, the force vector is

$$
\begin{align}
  F_0 &= \int_0^{1/2} \text{d}x \, x - 2 x^2 \\ &= \frac{1}{8} - \frac{1}{12} = \frac{1}{24}, \\
  F_1 &= \int_0^{1/2} \text{d}x \, 2 x^2 + \int_{1/2}^1 \text{d}x \, 2x - 2 x^2 \\ &= \frac{1}{12} + 1 - \frac{2}{3} - \frac{1}{4} + \frac{1}{12} = \frac{1}{4}.
\end{align}
$$


### Step 2.

We will demonstrate this by computing analytically the **stiffness** matrix and **force vector** using *sympy*:


```python
import sympy as smp
smp.init_printing()
```


```python
x = smp.symbols("x")
N = [
    smp.Piecewise((1-2*x, x <=1/2), (0, x > 1/2)),
    smp.Piecewise((2*x, x<=1/2), (2-2*x, x>1/2)),
    smp.Piecewise((0, x<1/2), (2*x-1, x>1/2))
]
K = np.zeros((2,2))
for A in range(2):
    for B in range(2):
        K[A,B] = smp.integrate(smp.diff(N[A], x)*smp.diff(N[B], x), (x, 0 , 1))
print(K)
```

    [[ 2. -2.]
     [-2.  4.]]


And now for the force vector:


```python
f = lambda x : x
F = np.zeros(2)
for A in range(2):
    F[A] = smp.integrate(N[A] * f(x), (x, 0, 1))
print(F)
```

    [0.04166667 0.25      ]


And finally we will **numerically** solve $KT=F$.


```python
T_A = np.linalg.solve(K,F)

N_num = []
for A, N_A in enumerate(N):
    N_num.append(smp.lambdify((x,), N[A], "numpy"))

def T(x):
    Tvals = np.zeros_like(x)
    for A in range(2):
        Tvals += T_A[A] * N_num[A](x)
    return Tvals
```


```python
coords = np.linspace(0,1,1000)
plt.plot(coords, T(coords), label="Numerical")
plt.plot(coords, (1-coords**3)/6, label="Exact")
plt.legend()
plt.show()
```


```python

```

527ddafb13ec8429d0ec1c0bb66310701b06d5d0
ipynb
Jupyter Notebook
Extras/Finite-Elements-Task.ipynb
terminal-w/PythonCourse
22b4c5390cf0fd6639270adde2b5ccc10f10af94
MIT
2022-03-10T17:02:22.000Z
2022-03-19T14:49:30.000Z
Extras/Finite-Elements-Task.ipynb
terminal-w/PythonCourse
22b4c5390cf0fd6639270adde2b5ccc10f10af94
MIT
Extras/Finite-Elements-Task.ipynb
terminal-w/PythonCourse
22b4c5390cf0fd6639270adde2b5ccc10f10af94
MIT
2021-05-14T14:11:48.000Z
2021-05-14T14:11:48.000Z
Qwen/Qwen-72B
1. YES
2. YES