Note
Go to the end to download the full example code.
MeshOptim2#
Mesh optimization with the ZZ1 criterion for a bending part.
error = 90.444 %
error = 89.449 %
error = 41.121 %
error = 24.547 %
error = 9.902 %
error = 3.761 %
error = 1.372 %
error = 0.441 %
Generate movie 1/8 (12.50 %) 769.94 ms
Generate movie 2/8 (25.00 %) 546.36 ms
Generate movie 3/8 (37.50 %) 450.24 ms
Generate movie 4/8 (50.00 %) 359.86 ms
Generate movie 5/8 (62.50 %) 275.99 ms
Generate movie 6/8 (75.00 %) 189.18 ms
Generate movie 7/8 (87.50 %) 104.41 ms
Generate movie 8/8 (100.00 %) 0.00 µs
13 import matplotlib.pyplot as plt
14 import numpy as np
15
16 from EasyFEA import (
17 Display,
18 Folder,
19 Models,
20 Tic,
21 Mesher,
22 ElemType,
23 Mesh,
24 Simulations,
25 Paraview,
26 PyVista,
27 )
28 from EasyFEA.Geoms import Point, Points
29 from EasyFEA.FEM import Calc_projector
30
31 if __name__ == "__main__":
32 Display.Clear()
33
34 # ----------------------------------------------
35 # Configuration
36 # ----------------------------------------------
37 dim = 2
38
39 # geom
40 L = 120
41 h = L * 2 / 3
42 b = h
43 r = h / (2 + 1e-2)
44 e = (L - 2 * r) / 2
45
46 # outputs
47 folder = Folder.Results_Dir()
48 plotProj = False
49 makeMovie = True
50 makeParaview = False
51
52 # load
53 P = 800 # N
54 lineLoad = P / h # N/mm
55 surfLoad = P / h / b # N/mm2
56
57 # criteria
58 threshold = (
59 1 / 100 if dim == 2 else 0.04
60 ) # Target error for the optimization process
61 iterMax = 20 # Maximum number of iterations
62 coef = 1 / 10 # Scaling coefficient for the optimization process
63
64 # ----------------------------------------------
65 # Mesh
66 # ----------------------------------------------
67 if dim == 2:
68 elemType = ElemType.TRI3
69 else:
70 elemType = ElemType.TETRA4
71 meshSize = h / 10
72
73 pt1 = Point()
74 pt2 = Point(e, 0)
75 pt3 = Point(e, r, r=r)
76 pt4 = Point(L - e, r, r=r)
77 pt5 = Point(L - e, 0)
78 pt6 = Point(L, 0)
79
80 pt7 = Point(L, h)
81 pt8 = Point(L - e, h)
82 pt9 = Point(L - e, h - r, r=r)
83 pt10 = Point(e, h - r, r=r)
84 pt11 = Point(e, h)
85 pt12 = Point(0, h)
86
87 points = Points(
88 [pt1, pt2, pt3, pt4, pt5, pt6, pt7, pt8, pt9, pt10, pt11, pt12], meshSize
89 )
90 inclusions = []
91
92 PyVista.Plot_Geoms(points.Get_Contour()).show()
93
94 def DoMesh(refineGeom=None) -> Mesh:
95 """Function used to generate the mesh."""
96 if dim == 2:
97 return Mesher().Mesh_2D(points, inclusions, elemType, [], [refineGeom])
98 else:
99 return Mesher().Mesh_Extrude(
100 points, inclusions, [0, 0, b], [], elemType, [], [refineGeom]
101 )
102
103 # Construct the initial mesh
104 mesh = DoMesh()
105 PyVista.Plot_Mesh(mesh).show()
106
107 # ----------------------------------------------
108 # Material and Simulation
109 # ----------------------------------------------
110 material = Models.Elastic.Isotropic(dim, E=210000, v=0.3, thickness=b)
111 simu = Simulations.Elastic(mesh, material)
112 simu.rho = 8100 * 1e-9
113
114 def DoSimu(refineGeom: str):
115 simu.mesh = DoMesh(refineGeom)
116
117 # get the nodes
118 nodes_Fixed = simu.mesh.Nodes_Conditions(lambda x, y, z: x == 0)
119 nodes_Load = simu.mesh.Nodes_Conditions(lambda x, y, z: x == L)
120
121 # do the simulation
122 simu.Bc_Init()
123 simu.add_dirichlet(
124 nodes_Fixed, [0] * dim, simu.Get_unknowns(), description="Fixed"
125 )
126 simu.add_surfLoad(nodes_Load, [-surfLoad], ["y"])
127
128 simu.Solve()
129
130 simu.Save_Iter()
131
132 return simu
133
134 simu = Simulations.Mesh_Optim_ZZ1(DoSimu, folder, threshold, iterMax, 1 / 10)
135 PyVista.Plot_BoundaryConditions(simu).show()
136
137 # ----------------------------------------------
138 # Results
139 # ----------------------------------------------
140 PyVista.Plot_Mesh(simu.mesh).show()
141 PyVista.Plot(simu, "ZZ1_e", nodeValues=False, nColors=11).show()
142
143 if plotProj:
144 simu.Set_Iter(0)
145 mesh0 = simu.mesh
146 u0 = np.reshape(simu.displacement, (mesh0.Nn, -1))
147
148 simu.Set_Iter(1)
149 mesh1 = simu.mesh
150
151 proj = Calc_projector(mesh0, mesh1)
152 uProj = np.zeros((mesh1.Nn, dim), dtype=float)
153 for d in range(dim):
154 uProj[:, d] = proj @ u0[:, d]
155
156 ax = Display.Plot_Result(
157 mesh0, np.linalg.norm(u0, axis=1), plotMesh=True, title="u0"
158 )
159 ax.plot(*mesh1.coord[:, :dim].T, ls="", marker="+", c="k", label="new nodes")
160 ax.legend()
161 Display.Plot_Result(
162 mesh1, np.linalg.norm(uProj, axis=1), plotMesh=True, title="uProj"
163 )
164
165 if makeParaview:
166 Paraview.Save_simu(simu, folder, nodeFields=["ZZ1_e"])
167
168 if makeMovie:
169
170 def func(plotter, n):
171 simu.Set_Iter(n)
172
173 PyVista.Plot_Mesh(simu, plotter=plotter)
174
175 zz1 = simu._Calc_ZZ1()[0]
176
177 plotter.add_title(f"ZZ1 = {zz1 * 100:.2f} %")
178
179 PyVista.Movie_func(func, len(simu.results), folder, "lmt.gif")
180
181 Tic.Plot_History()
182 plt.show()
Total running time of the script: (0 minutes 4.866 seconds)