Note
Go to the end to download the full example code.
Modal1#
Modal analysis of a wall structure.
12 from EasyFEA import (
13 Display,
14 Models,
15 np,
16 ElemType,
17 Simulations,
18 PyVista,
19 )
20 from EasyFEA.Geoms import Domain
21
22 from scipy.sparse import linalg, eye
23
24 if __name__ == "__main__":
25 Display.Clear()
26
27 # ----------------------------------------------
28 # Configuration
29 # ----------------------------------------------
30
31 dim = 3
32
33 # outputs
34 isFixed = True
35 Nmode = 3
36
37 # ----------------------------------------------
38 # Mesh
39 # ----------------------------------------------
40 contour = Domain((0, 0), (1, 1), 1 / 10)
41 thickness = 1 / 10
42
43 if dim == 2:
44 mesh = contour.Mesh_2D([], ElemType.QUAD9, isOrganised=True)
45 else:
46 mesh = contour.Mesh_Extrude(
47 [], [0, 0, -thickness], [2], ElemType.HEXA27, isOrganised=True
48 )
49 nodesY0 = mesh.Nodes_Conditions(lambda x, y, z: y == 0)
50 nodesSupY0 = mesh.Nodes_Conditions(lambda x, y, z: y > 0)
51
52 PyVista.Plot_Mesh(mesh).show()
53
54 # ----------------------------------------------
55 # Simulation
56 # ----------------------------------------------
57 material = Models.ElasIsot(dim, planeStress=True, thickness=thickness)
58
59 simu = Simulations.ElasticSimu(mesh, material)
60
61 simu.Solver_Set_Hyperbolic_Algorithm(0.1)
62
63 K, C, M, F = simu.Get_K_C_M_F()
64
65 if isFixed:
66 simu.add_dirichlet(nodesY0, [0] * dim, simu.Get_unknowns())
67 known, unknown = simu.Bc_dofs_known_unknown(simu.problemType)
68 K_t = K[unknown, :].tocsc()[:, unknown].tocsr()
69 M_t = M[unknown, :].tocsc()[:, unknown].tocsr()
70
71 else:
72 K_t = K + K.min() * eye(K.shape[0]) * 1e-12
73 M_t = M
74
75 # eigenValues, eigenVectors = linalg.eigs(K_t, Nmode, M_t, which="SR")
76 eigenValues, eigenVectors = linalg.eigs(K_t, Nmode, M_t, sigma=0, which="LR")
77
78 eigenValues = eigenValues.real
79 eigenVectors = eigenVectors.real
80 freq_t = np.sqrt(eigenValues) / 2 / np.pi
81
82 # ----------------------------------------------
83 # Plot modes
84 # ----------------------------------------------
85 for n in range(eigenValues.size):
86 if isFixed:
87 mode = np.zeros((mesh.Nn, dim))
88 mode[nodesSupY0, :] = np.reshape(eigenVectors[:, n], (-1, dim))
89 else:
90 mode = np.reshape(eigenVectors[:, n], (-1, dim))
91
92 simu._Set_solutions(simu.problemType, mode.ravel())
93 simu.Save_Iter()
94
95 sol = np.linalg.norm(mode, axis=1)
96 deformFactor = 1 / 5 / np.abs(sol).max()
97
98 plotter = PyVista.Plot(simu, alpha=0.5)
99 PyVista.Plot(simu, None, deformFactor, alpha=0.8, color="r", plotter=plotter)
100 plotter.add_title(f"mode {n + 1}")
101 plotter.show()
102
103 axModes = Display.Init_Axes()
104 axModes.plot(np.arange(eigenValues.size), freq_t, ls="", marker="+")
105 axModes.set_xlabel("modes")
106 axModes.set_ylabel("freq [Hz]")
107 axModes.grid()
108
109 Display.plt.show()
Total running time of the script: (0 minutes 2.852 seconds)