```python
from numpy import array
from sympy import symbols, Eq, Matrix, diff, Derivative, simplify, factor, expand, latex, init_printing, collect
init_printing()
from IPython.display import display, Math
```

## Define the interpolation functions


```python
x1, y1, x2, y2, x3, y3, x4, y4 = symbols('x_1, y_1, x_2, y_2, x_3, y_3, x_4, y_4')
r, s = symbols('r, s')

# Define the interpolation functions
h1 = factor(1/4*(1+r)*(1+s))
h2 = factor(1/4*(1-r)*(1+s))
h3 = factor(1/4*(1-r)*(1-s))
h4 = factor(1/4*(1+r)*(1-s))

display(Math('h_1 = ' + latex(h1)))
display(Math('h_2 = ' + latex(h2)))
display(Math('h_2 = ' + latex(h3)))
display(Math('h_2 = ' + latex(h4)))
```


$\displaystyle h_1 = 0.25 \left(r + 1\right) \left(s + 1\right)$



$\displaystyle h_2 = - 0.25 \left(r - 1\right) \left(s + 1\right)$



$\displaystyle h_2 = 0.25 \left(r - 1\right) \left(s - 1\right)$



$\displaystyle h_2 = - 0.25 \left(r + 1\right) \left(s - 1\right)$


# Relate the (x, y) Local Coordinate System to the (r, s) Natural Coordinate System


```python
# Relate the (x, y) coordinate system to the (r, s) coordinate system
x = h1*x1 + h2*x2 + h3*x3 + h4*x4
y = h1*y1 + h2*y2 + h3*y3 + h4*y4

display(Math('x = ' + latex(x)))
display(Math('y = ' + latex(y)))
```


$\displaystyle x = 0.25 x_{1} \left(r + 1\right) \left(s + 1\right) - 0.25 x_{2} \left(r - 1\right) \left(s + 1\right) + 0.25 x_{3} \left(r - 1\right) \left(s - 1\right) - 0.25 x_{4} \left(r + 1\right) \left(s - 1\right)$



$\displaystyle y = 0.25 y_{1} \left(r + 1\right) \left(s + 1\right) - 0.25 y_{2} \left(r - 1\right) \left(s + 1\right) + 0.25 y_{3} \left(r - 1\right) \left(s - 1\right) - 0.25 y_{4} \left(r + 1\right) \left(s - 1\right)$


# Find the Jacobian Matrix
The Jacobian matrix converts derivatives with respect to x and y into derivatives with respect to r and s

$J = \begin{bmatrix} \frac{dx}{dr} & \frac{dy}{dr} \\ \frac{dx}{ds} & \frac{dy}{ds} \end{bmatrix}$


```python
# Calculate the Jacobian matrix
J = Matrix([[diff(x, r), diff(y, r)],
            [diff(x, s), diff(y, s)]])

display(Math('J = ' + latex(factor(J))))
```


$\displaystyle J = \left[\begin{matrix}0.25 x_{1} \left(s + 1\right) - 0.25 x_{2} \left(s + 1\right) + 0.25 x_{3} \left(s - 1\right) - 0.25 x_{4} \left(s - 1\right) & 0.25 y_{1} \left(s + 1\right) - 0.25 y_{2} \left(s + 1\right) + 0.25 y_{3} \left(s - 1\right) - 0.25 y_{4} \left(s - 1\right)\\0.25 x_{1} \left(r + 1\right) - 0.25 x_{2} \left(r - 1\right) + 0.25 x_{3} \left(r - 1\right) - 0.25 x_{4} \left(r + 1\right) & 0.25 y_{1} \left(r + 1\right) - 0.25 y_{2} \left(r - 1\right) + 0.25 y_{3} \left(r - 1\right) - 0.25 y_{4} \left(r + 1\right)\end{matrix}\right]$



```python
B_kappa = Matrix([[0,      0,      -diff(h1, r), 0,      0,      -diff(h2, r), 0,      0,      -diff(h2, r), 0,      0,      -diff(h4, r)],
                  [0, diff(h1, s),       0,      0, diff(h2, s),       0,      0, diff(h3, s),       0,      0, diff(h4, s),       0     ],
                  [0, diff(h1, r), -diff(h1, s), 0, diff(h2, r), -diff(h2, s), 0, diff(h3, r), -diff(h3, s), 0, diff(h4, r), -diff(h4, s)]])

display(Math('B_\kappa = J^{-1}(1/4)' + latex(B_kappa*4)))
```


$\displaystyle B_\kappa = J^{-1}(1/4)\left[\begin{array}{cccccccccccc}0 & 0 & - 1.0 s - 1.0 & 0 & 0 & s + 1.0 & 0 & 0 & s + 1.0 & 0 & 0 & s - 1.0\\0 & r + 1.0 & 0 & 0 & 1.0 - 1.0 r & 0 & 0 & r - 1.0 & 0 & 0 & - 1.0 r - 1.0 & 0\\0 & s + 1.0 & - 1.0 r - 1.0 & 0 & - 1.0 s - 1.0 & r - 1.0 & 0 & s - 1.0 & 1.0 - 1.0 r & 0 & 1.0 - 1.0 s & r + 1.0\end{array}\right]$



```python
dH = Matrix([[diff(h1, r), diff(h2, r), diff(h3, r), diff(h4, r)],
             [diff(h1, s), diff(h2, s), diff(h3, s), diff(h4, s)]])

B_m = Matrix([[diff(h1, r),      0,      diff(h2, r),      0,      diff(h3, r),      0,      diff(h4, r),      0     ],
              [     0,      diff(h1, s),      0,      diff(h2, s),      0,      diff(h3, s),      0,      diff(h4, s)],
              [diff(h1, s), diff(h1, r), diff(h2, s), diff(h2, r), diff(h3, s), diff(h3, r), diff(h4, s), diff(h4, r)]])

display(Math('dH = (1/4)' + latex(dH*4)))
display(Math('B_m = J^{-1}(1/4)' + latex(B_m*4)))
print(B_m*4)
```


$\displaystyle dH = (1/4)\left[\begin{matrix}s + 1.0 & - 1.0 s - 1.0 & s - 1.0 & 1.0 - 1.0 s\\r + 1.0 & 1.0 - 1.0 r & r - 1.0 & - 1.0 r - 1.0\end{matrix}\right]$



$\displaystyle B_m = J^{-1}(1/4)\left[\begin{matrix}s + 1.0 & 0 & - 1.0 s - 1.0 & 0 & s - 1.0 & 0 & 1.0 - 1.0 s & 0\\0 & r + 1.0 & 0 & 1.0 - 1.0 r & 0 & r - 1.0 & 0 & - 1.0 r - 1.0\\r + 1.0 & s + 1.0 & 1.0 - 1.0 r & - 1.0 s - 1.0 & r - 1.0 & s - 1.0 & - 1.0 r - 1.0 & 1.0 - 1.0 s\end{matrix}\right]$


    Matrix([[s + 1.0, 0, -1.0*s - 1.0, 0, s - 1.0, 0, 1.0 - 1.0*s, 0], [0, r + 1.0, 0, 1.0 - 1.0*r, 0, r - 1.0, 0, -1.0*r - 1.0], [r + 1.0, s + 1.0, 1.0 - 1.0*r, -1.0*s - 1.0, r - 1.0, s - 1.0, -1.0*r - 1.0, 1.0 - 1.0*s]])



```python
r, s = symbols('r, s')
h = symbols('h')
x1, x2, x3, x4 = symbols('x1, x2, x3, x4')
y1, y2, y3, y4 = symbols('y1, y2, y3, y4')
w1, w2, w3, w4 = symbols('w1, w2, w3, w4')
theta_x1, theta_x2, theta_x3, theta_x4 = symbols('theta_x1, theta_x2, theta_x3, theta_x4')
theta_y1, theta_y2, theta_y3, theta_y4 = symbols('theta_y1, theta_y2, theta_y3, theta_y4')

gamma_rz = (1 + s)*((w1 - w2)/2 + (x1 - x2)/4*(theta_y1 + theta_y2) - (y1 - y2)/4*(theta_x1 + theta_x2)) \
         + (1 - s)*((w4 - w3)/2 + (x4 - x3)/4*(theta_y4 + theta_y3) - (y4 - y3)/4*(theta_x4 + theta_x3))

gamma_sz = (1 + r)*((w1 - w4)/2 + (x1 - x4)/4*(theta_y1 + theta_y4) - (y1 - y4)/4*(theta_x1 + theta_x4)) \
         + (1 - r)*((w2 - w3)/2 + (x2 - x3)/4*(theta_y2 + theta_y3) - (y2 - y3)/4*(theta_x2 + theta_x3))

display(1/4*collect(expand(gamma_rz*4), [w1, theta_x1, theta_y1, w2, theta_x2, theta_y2, w3, theta_x3, theta_y3, w4, theta_x4, theta_y4]))
print('')
display(1/4*collect(expand(gamma_sz*4), [w1, theta_x1, theta_y1, w2, theta_x2, theta_y2, w3, theta_x3, theta_y3, w4, theta_x4, theta_y4]))
```

# References

The following references were used in the formulation of this element. The primary reference was the first one because it had a section directly relating to isoparametric general plate bending elements, but all three were used. The first reference is a free download from MIT's website.

1. "Finite Element Procedures, 2nd Edition", Klaus-Jurgen Bathe
2. "Finite Element Analysis Fundamentals", Richard H. Gallagher
3. "A First Course in the Finite Element Method, 4th Edition", Daryl L. Logan

# Write the Equations that Map the (x, y) Coordinates to the (r, s) Coordinates

For any invertible matrix $M = \begin{bmatrix} a & b \\ c & d \end{bmatrix}$, the inverse can be found as follows:


```python
a, b, c, d = symbols('a, b, c, d')

M = Matrix([[a, b],
           [c, d]])

display(Math('M^{-1} = ' + latex(M.inv())))
```

Recognizing that the denominator in each term of this inverse is the determinant of $M$, the inverse can be rewritten as:

$M^{-1} = \frac{1}{|M|}\begin{bmatrix} d & -b \\ -c & a \end{bmatrix}$

Applying this to the Jacobian, we can express its the inverse as:

$J^{-1} = \frac{1}{|J|} \begin{bmatrix} \frac{dy}{ds} & \frac{-dy}{dr} \\ \frac{-dx}{ds} & \frac{dx}{dr} \end{bmatrix}$

We can relate the derivatives of the (r, s) coordinate system to the derivatives of the (x, y) coordinate system using the Jacobian operator:

$\begin{bmatrix} \frac{d}{dr} \\ \frac{d}{ds} \end{bmatrix} = \begin{bmatrix} \frac{dx}{dr} & \frac{dy}{dr} \\ \frac{dx}{ds} & \frac{dy}{ds}\end{bmatrix} \begin{bmatrix} \frac{d}{dx} \\ \frac{d}{dy} \end{bmatrix} = J \begin{bmatrix} \frac{d}{dx} \\ \frac{d}{dy} \end{bmatrix}$

Rearranging the equation gives us the formula for any derivative of the (x, y) coordinate system in terms of (r, s):

$\begin{bmatrix} \frac{d}{dx} \\ \frac{d}{dy} \end{bmatrix} = \frac{1}{|J|} \begin{bmatrix} \frac{dy}{ds} & \frac{-dy}{dr} \\ \frac{-dx}{ds} & \frac{dx}{dr} \end{bmatrix} \begin{bmatrix} \frac{d}{dr} \\ \frac{d}{ds} \end{bmatrix}$

Finally, multiplying through, we can obtain equations for the converting derivatives of any function with respect to x or y to derivatives with respect to r and s, respectively:

$\frac{d}{dx} = \frac{1}{|J|} (\frac{dy}{ds} \frac{d}{dr} + \frac{-dy}{dr} \frac{d}{ds})$

$\frac{d}{dy} = \frac{1}{|J|} (\frac{-dx}{ds} \frac{d}{dr} + \frac{dx}{dr} \frac{d}{ds})$

The bending strains vary through the thickness of the plate. Using the curvatures of the plate they can be expressed as:

$\begin{bmatrix} \epsilon_{xx} \\ \epsilon_{yy} \\ \gamma_{xy} \end{bmatrix} = -z \begin{bmatrix} \frac{d \beta_{x}}{dx} \\ \frac{d \beta_{y}}{dy} \\ \frac{d \beta_{x}}{dy} + \frac{d\beta_{y}}{dx} \end{bmatrix} = \frac{-z}{|J|} \begin{bmatrix} \frac{dy}{ds} \frac{d\beta_x}{dr} + \frac{-dy}{dr} \frac{d\beta_x}{ds} \\ \frac{-dx}{ds} \frac{d\beta_y}{dr} + \frac{dx}{dr} \frac{d\beta_y}{ds} \\ (\frac{-dx}{ds} \frac{d\beta_x}{dr} + \frac{dx}{dr} \frac{d\beta_x}{ds}) + (\frac{dy}{ds} \frac{d\beta_y}{dr} + \frac{-dy}{dr} \frac{d\beta_y}{ds})\end{bmatrix}$


```python

```


```python
w_zi, theta_xi, theta_yi = symbols('w_{zi}, \\theta_{xi}, \\theta_{yi}')
w_zj, theta_xj, theta_yj = symbols('w_{zj}, \\theta_{xj}, \\theta_{yj}')
w_zm, theta_xm, theta_ym = symbols('w_{zm}, \\theta_{xm}, \\theta_{ym}')
w_zn, theta_xn, theta_yn = symbols('w_{zn}, \\theta_{xn}, \\theta_{yn}')

d = Matrix([w_zi, theta_xi, theta_yi,
            w_zj, theta_xj, theta_yj,
            w_zm, theta_xm, theta_ym,
            w_zn, theta_xn, theta_yn])
d
```


```python
w = Matrix([h1, 0, 0, h2, 0, 0, h3, 0, 0, h4, 0, 0]).T*d
beta_x = -1*Matrix([0, 0, h1, 0, 0, h2, 0, 0, h3, 0, 0, h4]).T*d
beta_y = Matrix([0, h1, 0, 0, h2, 0, 0, h3, 0, 0, h4, 0]).T*d
display(w, beta_x, beta_y)
```


```python
beta_x_dx = 1/J.det()*(diff(y, s)*diff(beta_x, r) - diff(y, r)*diff(beta_x, s))
beta_x_dy = 1/J.det()*(-diff(x, s)*diff(beta_x, r) + diff(x, r)*diff(beta_x, s))
beta_y_dx = 1/J.det()*(diff(y, s)*diff(beta_y, r) - diff(y, r)*diff(beta_y, s))
beta_y_dy = 1/J.det()*(-diff(x, s)*diff(beta_y, r) + diff(x, r)*diff(beta_y, s))
```


```python
B_kappa = Matrix([[0, 0,           diff(-h1, r), 0, 0,           diff(-h2, r), 0, 0,           diff(-h3, r), 0, 0,           diff(-h4, r)],
                  [0, diff(h1, s), 0,            0, diff(h2, s), 0,            0, diff(h3, s), 0,            0, diff(h4, s), 0],
                  [0, diff(h1, r), diff(-h1, s), 0, diff(h2, r), diff(-h2, s), 0, diff(h3, r), diff(-h3, s), 0, diff(h4, r), diff(-h4, s)]])
B_kappa
```


```python
x = 
gr = diff(x, r)
```

a12ed697a2b94e82b13392791cf8cca9f8e91f0c
ipynb
Jupyter Notebook
Derivations/Isoparametric Quad Element.ipynb
AmirHosseinNamadchi/PyNite
8cc1fe3262e1efe029c6860394d2436601272e33
MIT
Derivations/Isoparametric Quad Element.ipynb
AmirHosseinNamadchi/PyNite
8cc1fe3262e1efe029c6860394d2436601272e33
MIT
Derivations/Isoparametric Quad Element.ipynb
AmirHosseinNamadchi/PyNite
8cc1fe3262e1efe029c6860394d2436601272e33
MIT
Qwen/Qwen-72B
1. YES
2. YES