```python
import numpy as np
import sympy as sym

x_1, x_2, x_3, x_4 = sym.symbols('x_1, x_2, x_3, x_4')
y_1, y_2, y_3, y_4 = sym.symbols('y_1, y_2, y_3, y_4')
XY = sym.Matrix([sym.symbols('x_1, x_2, x_3, x_4'), sym.symbols('y_1, y_2, y_3, y_4')]).transpose()
xi, eta = sym.symbols('xi, eta')
basis = sym.Matrix([xi, eta])
N1 = (1-xi)*(1-eta)/4
N2 = (1+xi)*(1-eta)/4
N3 = (1+xi)*(1+eta)/4
N4 = (1-xi)*(1+eta)/4
NN = sym.Matrix([N1, N2, N3, N4])

coordinate = (NN.transpose()*XY).transpose()

NN_diff = NN.jacobian(basis)
jacob = coordinate.jacobian(basis)

jacob_det = sym.det(jacob)
jacob_inv = sym.inv_quick(jacob)

NN_diff_global = (NN_diff*jacob_inv)
temp = NN_diff_global*NN_diff_global.transpose()/jacob_det

sq3 = np.sqrt(1/3)
```


```python
NN
```




$\displaystyle \left[\begin{matrix}\frac{\left(1 - \eta\right) \left(1 - \xi\right)}{4}\\\frac{\left(1 - \eta\right) \left(\xi + 1\right)}{4}\\\frac{\left(\eta + 1\right) \left(\xi + 1\right)}{4}\\\frac{\left(1 - \xi\right) \left(\eta + 1\right)}{4}\end{matrix}\right]$




```python
coordinate
```




$\displaystyle \left[\begin{matrix}\frac{x_{1} \left(1 - \eta\right) \left(1 - \xi\right)}{4} + \frac{x_{2} \left(1 - \eta\right) \left(\xi + 1\right)}{4} + \frac{x_{3} \left(\eta + 1\right) \left(\xi + 1\right)}{4} + \frac{x_{4} \left(1 - \xi\right) \left(\eta + 1\right)}{4}\\\frac{y_{1} \left(1 - \eta\right) \left(1 - \xi\right)}{4} + \frac{y_{2} \left(1 - \eta\right) \left(\xi + 1\right)}{4} + \frac{y_{3} \left(\eta + 1\right) \left(\xi + 1\right)}{4} + \frac{y_{4} \left(1 - \xi\right) \left(\eta + 1\right)}{4}\end{matrix}\right]$




```python
jacob
```




$\displaystyle \left[\begin{matrix}- \frac{x_{1} \left(1 - \eta\right)}{4} + \frac{x_{2} \left(1 - \eta\right)}{4} + \frac{x_{3} \left(\eta + 1\right)}{4} - \frac{x_{4} \left(\eta + 1\right)}{4} & - \frac{x_{1} \left(1 - \xi\right)}{4} - \frac{x_{2} \left(\xi + 1\right)}{4} + \frac{x_{3} \left(\xi + 1\right)}{4} + \frac{x_{4} \left(1 - \xi\right)}{4}\\- \frac{y_{1} \left(1 - \eta\right)}{4} + \frac{y_{2} \left(1 - \eta\right)}{4} + \frac{y_{3} \left(\eta + 1\right)}{4} - \frac{y_{4} \left(\eta + 1\right)}{4} & - \frac{y_{1} \left(1 - \xi\right)}{4} - \frac{y_{2} \left(\xi + 1\right)}{4} + \frac{y_{3} \left(\xi + 1\right)}{4} + \frac{y_{4} \left(1 - \xi\right)}{4}\end{matrix}\right]$




```python
jacob.subs(([x_1, -0.7], [x_2, 0.2], [x_3, 0.8], [x_4, 0], [y_1, 0.3], [y_2, -0.8], [y_3, 0], [y_4, 0.8]))
```




$\displaystyle \left[\begin{matrix}0.425 - 0.025 \eta & 0.325 - 0.025 \xi\\0.075 \eta - 0.475 & 0.075 \xi + 0.325\end{matrix}\right]$




```python
jacob.subs(([xi, sq3], [eta, -sq3], [x_1, -0.8], [x_2, 0], [x_3, 0.8], [x_4, 0], [y_1, 0], [y_2, -0.8], [y_3, 0], [y_4, 0.8]))
```




$\displaystyle \left[\begin{matrix}0.5 & 0.5\\-0.5 & 0.5\end{matrix}\right]$




```python
jaco_inv = sym.inv_quick(jaco)
jaco_inv
```




$\displaystyle \left[\begin{matrix}1.25 & -1.25\\1.25 & 1.25\end{matrix}\right]$




```python
NN_diff.subs(([xi, sq3], [eta, -sq3]))
```




$\displaystyle \left[\begin{matrix}-0.394337567297406 & -0.105662432702594\\0.394337567297406 & -0.394337567297406\\0.105662432702594 & 0.394337567297406\\-0.105662432702594 & 0.105662432702594\end{matrix}\right]$




```python
NN_diff
```




$\displaystyle \left[\begin{matrix}\frac{\eta}{4} - \frac{1}{4} & \frac{\xi}{4} - \frac{1}{4}\\\frac{1}{4} - \frac{\eta}{4} & - \frac{\xi}{4} - \frac{1}{4}\\\frac{\eta}{4} + \frac{1}{4} & \frac{\xi}{4} + \frac{1}{4}\\- \frac{\eta}{4} - \frac{1}{4} & \frac{1}{4} - \frac{\xi}{4}\end{matrix}\right]$




```python
jacob_det
```




$\displaystyle - \frac{\eta x_{1} y_{2}}{8} + \frac{\eta x_{1} y_{3}}{8} + \frac{\eta x_{2} y_{1}}{8} - \frac{\eta x_{2} y_{4}}{8} - \frac{\eta x_{3} y_{1}}{8} + \frac{\eta x_{3} y_{4}}{8} + \frac{\eta x_{4} y_{2}}{8} - \frac{\eta x_{4} y_{3}}{8} - \frac{x_{1} \xi y_{3}}{8} + \frac{x_{1} \xi y_{4}}{8} + \frac{x_{1} y_{2}}{8} - \frac{x_{1} y_{4}}{8} + \frac{x_{2} \xi y_{3}}{8} - \frac{x_{2} \xi y_{4}}{8} - \frac{x_{2} y_{1}}{8} + \frac{x_{2} y_{3}}{8} + \frac{x_{3} \xi y_{1}}{8} - \frac{x_{3} \xi y_{2}}{8} - \frac{x_{3} y_{2}}{8} + \frac{x_{3} y_{4}}{8} - \frac{x_{4} \xi y_{1}}{8} + \frac{x_{4} \xi y_{2}}{8} + \frac{x_{4} y_{1}}{8} - \frac{x_{4} y_{3}}{8}$




```python
jacob_inv
```




$\displaystyle \left[\begin{matrix}\frac{- \frac{y_{1} \left(1 - \xi\right)}{4} - \frac{y_{2} \left(\xi + 1\right)}{4} + \frac{y_{3} \left(\xi + 1\right)}{4} + \frac{y_{4} \left(1 - \xi\right)}{4}}{\left(- \frac{x_{1} \left(1 - \eta\right)}{4} + \frac{x_{2} \left(1 - \eta\right)}{4} + \frac{x_{3} \left(\eta + 1\right)}{4} - \frac{x_{4} \left(\eta + 1\right)}{4}\right) \left(- \frac{y_{1} \left(1 - \xi\right)}{4} - \frac{y_{2} \left(\xi + 1\right)}{4} + \frac{y_{3} \left(\xi + 1\right)}{4} + \frac{y_{4} \left(1 - \xi\right)}{4}\right) - \left(- \frac{x_{1} \left(1 - \xi\right)}{4} - \frac{x_{2} \left(\xi + 1\right)}{4} + \frac{x_{3} \left(\xi + 1\right)}{4} + \frac{x_{4} \left(1 - \xi\right)}{4}\right) \left(- \frac{y_{1} \left(1 - \eta\right)}{4} + \frac{y_{2} \left(1 - \eta\right)}{4} + \frac{y_{3} \left(\eta + 1\right)}{4} - \frac{y_{4} \left(\eta + 1\right)}{4}\right)} & \frac{\frac{x_{1} \left(1 - \xi\right)}{4} + \frac{x_{2} \left(\xi + 1\right)}{4} - \frac{x_{3} \left(\xi + 1\right)}{4} - \frac{x_{4} \left(1 - \xi\right)}{4}}{\left(- \frac{x_{1} \left(1 - \eta\right)}{4} + \frac{x_{2} \left(1 - \eta\right)}{4} + \frac{x_{3} \left(\eta + 1\right)}{4} - \frac{x_{4} \left(\eta + 1\right)}{4}\right) \left(- \frac{y_{1} \left(1 - \xi\right)}{4} - \frac{y_{2} \left(\xi + 1\right)}{4} + \frac{y_{3} \left(\xi + 1\right)}{4} + \frac{y_{4} \left(1 - \xi\right)}{4}\right) - \left(- \frac{x_{1} \left(1 - \xi\right)}{4} - \frac{x_{2} \left(\xi + 1\right)}{4} + \frac{x_{3} \left(\xi + 1\right)}{4} + \frac{x_{4} \left(1 - \xi\right)}{4}\right) \left(- \frac{y_{1} \left(1 - \eta\right)}{4} + \frac{y_{2} \left(1 - \eta\right)}{4} + \frac{y_{3} \left(\eta + 1\right)}{4} - \frac{y_{4} \left(\eta + 1\right)}{4}\right)}\\\frac{\frac{y_{1} \left(1 - \eta\right)}{4} - \frac{y_{2} \left(1 - \eta\right)}{4} - \frac{y_{3} \left(\eta + 1\right)}{4} + \frac{y_{4} \left(\eta + 1\right)}{4}}{\left(- \frac{x_{1} \left(1 - \eta\right)}{4} + \frac{x_{2} \left(1 - \eta\right)}{4} + \frac{x_{3} \left(\eta + 1\right)}{4} - \frac{x_{4} \left(\eta + 1\right)}{4}\right) \left(- \frac{y_{1} \left(1 - \xi\right)}{4} - \frac{y_{2} \left(\xi + 1\right)}{4} + \frac{y_{3} \left(\xi + 1\right)}{4} + \frac{y_{4} \left(1 - \xi\right)}{4}\right) - \left(- \frac{x_{1} \left(1 - \xi\right)}{4} - \frac{x_{2} \left(\xi + 1\right)}{4} + \frac{x_{3} \left(\xi + 1\right)}{4} + \frac{x_{4} \left(1 - \xi\right)}{4}\right) \left(- \frac{y_{1} \left(1 - \eta\right)}{4} + \frac{y_{2} \left(1 - \eta\right)}{4} + \frac{y_{3} \left(\eta + 1\right)}{4} - \frac{y_{4} \left(\eta + 1\right)}{4}\right)} & \frac{- \frac{x_{1} \left(1 - \eta\right)}{4} + \frac{x_{2} \left(1 - \eta\right)}{4} + \frac{x_{3} \left(\eta + 1\right)}{4} - \frac{x_{4} \left(\eta + 1\right)}{4}}{\left(- \frac{x_{1} \left(1 - \eta\right)}{4} + \frac{x_{2} \left(1 - \eta\right)}{4} + \frac{x_{3} \left(\eta + 1\right)}{4} - \frac{x_{4} \left(\eta + 1\right)}{4}\right) \left(- \frac{y_{1} \left(1 - \xi\right)}{4} - \frac{y_{2} \left(\xi + 1\right)}{4} + \frac{y_{3} \left(\xi + 1\right)}{4} + \frac{y_{4} \left(1 - \xi\right)}{4}\right) - \left(- \frac{x_{1} \left(1 - \xi\right)}{4} - \frac{x_{2} \left(\xi + 1\right)}{4} + \frac{x_{3} \left(\xi + 1\right)}{4} + \frac{x_{4} \left(1 - \xi\right)}{4}\right) \left(- \frac{y_{1} \left(1 - \eta\right)}{4} + \frac{y_{2} \left(1 - \eta\right)}{4} + \frac{y_{3} \left(\eta + 1\right)}{4} - \frac{y_{4} \left(\eta + 1\right)}{4}\right)}\end{matrix}\right]$




```python
D
```




$\displaystyle \left[\begin{matrix}1 & 1\\1 & 1\end{matrix}\right]$




```python
a = NN_diff_global.subs(([xi, -sq3], [eta, -sq3], [x_1, -0.8], [x_2, 0], [x_3, 0.8], [x_4, 0], [y_1, 0], [y_2, -0.8], [y_3, 0], [y_4, 0.8]))
```


```python
b = NN_diff_global.subs(([xi, sq3], [eta, -sq3], [x_1, -0.8], [x_2, 0], [x_3, 0.8], [x_4, 0], [y_1, 0], [y_2, -0.8], [y_3, 0], [y_4, 0.8]))
```


```python
c = NN_diff_global.subs(([xi, sq3], [eta, +sq3], [x_1, -0.8], [x_2, 0], [x_3, 0.8], [x_4, 0], [y_1, 0], [y_2, -0.8], [y_3, 0], [y_4, 0.8]))
```


```python
d = NN_diff_global.subs(([xi, -sq3], [eta, +sq3], [x_1, -0.8], [x_2, 0], [x_3, 0.8], [x_4, 0], [y_1, 0], [y_2, -0.8], [y_3, 0], [y_4, 0.8]))
```


```python
a
```




$\displaystyle \left[\begin{matrix}-0.985843918243516 & 0\\0.360843918243516 & -0.625\\0.264156081756484 & 0\\0.360843918243516 & 0.625\end{matrix}\right]$




```python
aa = np.array(a, dtype=np.float32)
bb = np.array(b, dtype=np.float32)
cc = np.array(c, dtype=np.float32)
dd = np.array(d, dtype=np.float32)
```


```python
aa
```




    array([[-0.9858439 ,  0.        ],
           [ 0.36084393, -0.625     ],
           [ 0.26415607,  0.        ],
           [ 0.36084393,  0.625     ]], dtype=float32)




```python
np.einsum('mi,ni->imn', aa, aa)*0.32
```




    array([[[ 0.31100422, -0.11383545, -0.08333333, -0.11383545],
            [-0.11383545,  0.04166667,  0.03050212,  0.04166667],
            [-0.08333333,  0.03050212,  0.0223291 ,  0.03050212],
            [-0.11383545,  0.04166667,  0.03050212,  0.04166667]],
    
           [[ 0.        ,  0.        ,  0.        ,  0.        ],
            [ 0.        ,  0.125     ,  0.        , -0.125     ],
            [ 0.        ,  0.        ,  0.        ,  0.        ],
            [ 0.        , -0.125     ,  0.        ,  0.125     ]]],
          dtype=float32)




```python
(np.einsum('mi,ni->imn', aa, aa)+np.einsum('mi,ni->imn', bb, bb)+\
np.einsum('mi,ni->imn', cc, cc)+np.einsum('mi,ni->imn', dd, dd))*0.32
```




    array([[[ 0.5833333 , -0.08333333, -0.41666663, -0.08333333],
            [-0.08333333,  0.08333334, -0.08333333,  0.08333334],
            [-0.41666663, -0.08333333,  0.5833333 , -0.08333333],
            [-0.08333333,  0.08333334, -0.08333333,  0.08333334]],
    
           [[ 0.08333334, -0.08333333,  0.08333334, -0.08333333],
            [-0.08333333,  0.5833333 , -0.08333333, -0.41666663],
            [ 0.08333334, -0.08333333,  0.08333334, -0.08333333],
            [-0.08333333, -0.41666663, -0.08333333,  0.5833333 ]]],
          dtype=float32)




```python

```

e11011c15250383c3c63214d0a31e762f294cf01
ipynb
Jupyter Notebook
equations4twoDimensionalElement.ipynb
AndrewWangJZ/pyfem
8e7df6aa69c1c761bb8ec67302847e30a83190b4
MIT
2022-03-10T17:22:53.000Z
2022-03-10T17:22:53.000Z
equations4twoDimensionalElement.ipynb
AndrewWangJZ/pyfem
8e7df6aa69c1c761bb8ec67302847e30a83190b4
MIT
equations4twoDimensionalElement.ipynb
AndrewWangJZ/pyfem
8e7df6aa69c1c761bb8ec67302847e30a83190b4
MIT
2022-03-10T12:47:34.000Z
2022-03-10T13:25:18.000Z
Qwen/Qwen-72B
1. YES
2. YES