# Linear 2D solution

## Init symbols for *sympy*


```python
from sympy import *
from geom_util import *
from sympy.vector import CoordSys3D
import matplotlib.pyplot as plt
import sys
sys.path.append("../")

%matplotlib inline

%reload_ext autoreload
%autoreload 2
%aimport geom_util
```


```python
# Any tweaks that normally go in .matplotlibrc, etc., should explicitly go here
%config InlineBackend.figure_format='retina'
plt.rcParams['figure.figsize'] = (12, 12)

plt.rc('text', usetex=True)
    
plt.rc('font', family='serif')

init_printing()
```


```python
N = CoordSys3D('N')
alpha1, alpha2, alpha3 = symbols("alpha_1 alpha_2 alpha_3", real = True, positive=True)
A,K,rho = symbols("A K rho")
```


```python
B=Matrix([[0, 1/(A*(K*alpha3 + 1)), 0, 0, 0, 0, 0, 0, K/(K*alpha3 + 1), 0, 0, 0], [0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 1/(A*(K*alpha3 + 1)), 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0], [-K/(K*alpha3 + 1), 0, 0, 0, 0, 0, 0, 0, 0, 1/(A*(K*alpha3 + 1)), 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1]])
B
```


```python
E=zeros(6,9)
E[0,0]=1
E[1,4]=1
E[2,8]=1
E[3,1]=1
E[3,3]=1
E[4,2]=1
E[4,6]=1
E[5,5]=1
E[5,7]=1
E
```


```python
mu = Symbol('mu')
la = Symbol('lambda')
C_tensor = getIsotropicStiffnessTensor(mu, la)
C = convertStiffnessTensorToMatrix(C_tensor)
C
```


```python
S=B.T*E.T*C*E*B*A*(1+alpha3*K)
S
```


```python
M=Matrix([[rho, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, rho, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, rho, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]])
M=M*A*(1+alpha3*K)
M
```

## Cartesian coordinates


```python
import fem.geometry as g
import fem.model as m
import fem.material as mat
import fem.solver as s
import fem.mesh as me
import plot

stiffness_matrix_func = lambdify([A, K, mu, la, alpha3], S, "numpy")
mass_matrix_func = lambdify([A, K, rho, alpha3], M, "numpy")


def stiffness_matrix(material, geometry, x1, x2, x3):
    A,K = geometry.get_A_and_K(x1,x2,x3)
    return stiffness_matrix_func(A, K, material.mu(), material.lam(), x3)

def mass_matrix(material, geometry, x1, x2, x3):
    A,K = geometry.get_A_and_K(x1,x2,x3)
    return mass_matrix_func(A, K, material.rho, x3)



def generate_layers(thickness, layers_count, material):
    layer_top = thickness / 2
    layer_thickness = thickness / layers_count
    layers = set()
    for i in range(layers_count):
        layer = m.Layer(layer_top - layer_thickness, layer_top, material, i)
        layers.add(layer)
        layer_top -= layer_thickness
    return layers


def solve(geometry, thickness, linear, N_width, N_height):
    layers_count = 1
    layers = generate_layers(thickness, layers_count, mat.IsotropicMaterial.steel())
    model = m.Model(geometry, layers, m.Model.FIXED_BOTTOM_LEFT_RIGHT_POINTS)
    mesh = me.Mesh.generate(width, layers, N_width, N_height, m.Model.FIXED_BOTTOM_LEFT_RIGHT_POINTS)
    lam, vec = s.solve(model, mesh, stiffness_matrix, mass_matrix)
    return lam, vec, mesh, geometry


width = 2
curvature = 0.8
thickness = 0.05

corrugation_amplitude = 0.05
corrugation_frequency = 20

geometry = g.General(width, curvature, corrugation_amplitude, corrugation_frequency)

N_width = 50
N_height = 4


lam, vec, mesh, geometry = solve(geometry, thickness, False, N_width, N_height)
results = s.convert_to_results(lam, vec, mesh, geometry)

results_index = 0
    
plot.plot_init_and_deformed_geometry_in_cartesian(results[results_index], 0, width, -thickness / 2, thickness / 2, 0, geometry.to_cartesian_coordinates)    
    
to_print = 20
if (len(results) < to_print):
    to_print = len(results)
    
for i in range(to_print):
    print(results[i].rad_per_sec_to_Hz(results[i].freq))

```

fdd3f8ec9d344fd859ca9a9ed2f20c1af985ba2c
ipynb
Jupyter Notebook
py/notebooks/.ipynb_checkpoints/LinearSol2DFEM-checkpoint.ipynb
tarashor/vibrations
fcbadd545e2a8f1ef4b4bff399b396ee5f16924b
MIT
2018-02-05T18:05:08.000Z
2018-02-05T18:05:08.000Z
py/notebooks/.ipynb_checkpoints/LinearSol2DFEM-checkpoint.ipynb
tarashor/vibrations
fcbadd545e2a8f1ef4b4bff399b396ee5f16924b
MIT
py/notebooks/.ipynb_checkpoints/LinearSol2DFEM-checkpoint.ipynb
tarashor/vibrations
fcbadd545e2a8f1ef4b4bff399b396ee5f16924b
MIT
Qwen/Qwen-72B
1. YES
2. YES