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