```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)))
```

# 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)))
```

# 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))))
```


```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)))
```


```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)
```

# Calculate the Tensor Shear Strain Components
The tensor shear strain components can be calculated using Reference 1, equations 5.102:

$\gamma_{xz}=\gamma_{rz}sin{\beta}-\gamma_{sz}sin{\alpha}$

$\gamma_{yz}=-\gamma_{rz}cos{\beta}+\gamma_{sz}cos{\alpha}$

Where the valus for $\gamma_{rz}$ and $\gamma_{sz}$ can be found from equations 5.103:

$\gamma_{rz}=\frac{\sqrt{(C_x+rB_x)^2+(C_y+rB_y)^2}}{8|J|}([1+s][\frac{w_1-w_2}{2}+\frac{x_1-x_2}{4}(\theta_y^1+\theta_y^2)-\frac{y_1-y_2}{4}(\theta_x^1+\theta_x^2)]+[1-s][\frac{w_4-w_3}{2}+\frac{x_4-x_3}{4}(\theta_y^4+\theta_y^3)-\frac{y_4-y_3}{4}(\theta_x^4+\theta_x^3)])$

$\gamma_{sz}=\frac{\sqrt{(A_x+sB_x)^2+(A_y+sB_y)^2}}{8|J|}([1+r][\frac{w_1-w_4}{2}+\frac{x_1-x_4}{4}(\theta_y^1+\theta_y^4)-\frac{y_1-y_4}{4}(\theta_x^1+\theta_x^4)]+[1-r][\frac{w_2-w_3}{2}+\frac{x_2-x_3}{4}(\theta_y^2+\theta_y^3)-\frac{y_2-y_3}{4}(\theta_x^2+\theta_x^3)])$

In order to express the tensor shear strains in matrix form we'll need to collect the terms associated with the displacement matrix [d]. The coefficients associated with these terms form the [B] matrix for shear.

$\gamma = \begin{bmatrix} \gamma_{xz} \\ \gamma_{yz} \end{bmatrix} = \begin{bmatrix} B_{\gamma_{xz}} \\ B_{\gamma_{yz}} \end{bmatrix} \begin{bmatrix} w_1 \\ \theta_x^1 \\ \theta_y^1 \\ w_2 \\ \theta_x^2 \\ \theta_y^2 \\ w_3 \\ \theta_x^3 \\ \theta_y^3 \\ w_4 \\ \theta_x^4 \\ \theta_y^4 \end{bmatrix}$

It's easier to convert $\gamma_{rz}$ and $\gamma_{sz}$ (as opposed to $\gamma_{xz}$ and $\gamma_{yz}$) into matrix format by collecting [d] terms, and then apply equations 5.102 afterward:


```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_terms = (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_terms = (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_terms*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_terms*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

c3520b4717f50ba519b53b54a22eca15ea21698d
ipynb
Jupyter Notebook
Derivations/MITC4 Quad Element.ipynb
Simple-Structures/PyNite
42417c57200d55fced6a89154823dd68e7cd15af
MIT
2021-03-19T12:02:40.000Z
2021-03-19T12:02:40.000Z
Derivations/MITC4 Quad Element.ipynb
Simple-Structures/PyNite
42417c57200d55fced6a89154823dd68e7cd15af
MIT
Derivations/MITC4 Quad Element.ipynb
Simple-Structures/PyNite
42417c57200d55fced6a89154823dd68e7cd15af
MIT
Qwen/Qwen-72B
1. YES
2. YES