Note
Go to the end to download the full example code.
MeshOptim3#
Mesh optimization using the ZZ1 criterion for a letter weigher.
error = 80.031 %
error = 63.003 %
error = 34.483 %
error = 19.716 %
error = 9.033 %
error = 2.937 %
error = 1.279 %
error = 0.694 %
Generate movie 1/8 (12.50 %) 807.71 ms
Generate movie 2/8 (25.00 %) 609.85 ms
Generate movie 3/8 (37.50 %) 504.05 ms
Generate movie 4/8 (50.00 %) 394.28 ms
Generate movie 5/8 (62.50 %) 304.59 ms
Generate movie 6/8 (75.00 %) 234.47 ms
Generate movie 7/8 (87.50 %) 136.68 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 ElemType,
22 Mesh,
23 Simulations,
24 Paraview,
25 PyVista,
26 )
27 from EasyFEA.Geoms import Point, Points
28 from EasyFEA.FEM import Calc_projector
29
30 if __name__ == "__main__":
31 Display.Clear()
32
33 # ----------------------------------------------
34 # Configuration
35 # ----------------------------------------------
36 dim = 2
37
38 # outputs
39 folder = Folder.Results_Dir()
40 plotProj = False
41 makeMovie = True
42 makeParaview = False
43
44 # geom
45 L = 80
46 h1 = L / 4
47 e1 = h1 * 0.1
48 h2 = (h1 - e1) * 0.95
49 e2 = (h1 - e1 - h2) / 2
50 r = h2 / 4
51 l = L / 2
52 b = h1
53
54 # load
55 P = 800 # N
56 lineLoad = P / h1 # N/mm
57 surfLoad = P / h1 / b # N/mm2
58
59 # criteria
60 threshold = (
61 1 / 100 if dim == 2 else 0.04
62 ) # Target error for the optimization process
63 iterMax = 20 # Maximum number of iterations
64 coef = 1 / 10 # Scaling coefficient for the optimization process
65
66 # ----------------------------------------------
67 # Mesh
68 # ----------------------------------------------
69 if dim == 2:
70 elemType = ElemType.TRI3
71 else:
72 elemType = ElemType.TETRA4
73 meshSize = r / 3
74
75 pt1 = Point()
76 pt2 = Point(L - h1)
77 pt3 = pt2 + [e1, -e1]
78 pt4 = Point(L, -e1)
79 pt5 = pt4 + [0, h1]
80 pt6 = Point(h1, h1 - e1)
81 pt7 = pt6 + [-e1, e1]
82 pt8 = pt1 + [0, h1]
83
84 points = Points([pt1, pt2, pt3, pt4, pt5, pt6, pt7, pt8], meshSize)
85
86 p1 = Point(L / 2 - l / 2, e2, r=r)
87 p2 = Point(L / 2 - l / 2 + 2 * r, e2, r=r)
88 p3 = Point(L / 2 - l / 2 + 2 * r, e2 + r)
89 p4 = p3.copy()
90 p4.Symmetry((L / 2, (h1 - e1) / 2), (1, 0))
91 p5 = p2.copy()
92 p5.Symmetry((L / 2, (h1 - e1) / 2), (1, 0))
93 p6 = p1.copy()
94 p6.Symmetry((L / 2, (h1 - e1) / 2), (1, 0))
95 p7 = p6.copy()
96 p7.Symmetry((L / 2, (h1 - e1) / 2), (0, 1))
97 p8 = p5.copy()
98 p8.Symmetry((L / 2, (h1 - e1) / 2), (0, 1))
99 p9 = p4.copy()
100 p9.Symmetry((L / 2, (h1 - e1) / 2), (0, 1))
101 p10 = p3.copy()
102 p10.Symmetry((L / 2, (h1 - e1) / 2), (0, 1))
103 p11 = p2.copy()
104 p11.Symmetry((L / 2, (h1 - e1) / 2), (0, 1))
105 p12 = p1.copy()
106 p12.Symmetry((L / 2, (h1 - e1) / 2), (0, 1))
107
108 inclusion = Points(
109 [p1, p2, p3, p4, p5, p6, p7, p8, p9, p10, p11, p12], meshSize, isFilled=False
110 )
111 inclusions = [inclusion]
112
113 PyVista.Plot_Geoms([points.Get_Contour(), inclusion.Get_Contour()]).show()
114
115 def DoMesh(refineGeom=None) -> Mesh:
116 """Function used to generate the mesh."""
117 if dim == 2:
118 return points.Mesh_2D(inclusions, elemType, [], [refineGeom])
119 else:
120 return points.Mesh_Extrude(
121 inclusions, [0, 0, b], [], elemType, [], [refineGeom]
122 )
123
124 # Construct the initial mesh
125 mesh = DoMesh()
126 PyVista.Plot_Mesh(mesh).show()
127
128 # ----------------------------------------------
129 # Material and Simulation
130 # ----------------------------------------------
131 material = Models.Elastic.Isotropic(dim, E=210000, v=0.3, thickness=b)
132 simu = Simulations.Elastic(mesh, material)
133 simu.rho = 8100 * 1e-9
134
135 def DoSimu(refineGeom: str):
136 simu.mesh = DoMesh(refineGeom)
137
138 # get the nodes
139 nodes_Fixed = simu.mesh.Nodes_Conditions(lambda x, y, z: x == 0)
140 nodes_Load = simu.mesh.Nodes_Conditions(lambda x, y, z: x == L)
141
142 # do the simulation
143 simu.Bc_Init()
144 simu.add_dirichlet(
145 nodes_Fixed, [0] * dim, simu.Get_unknowns(), description="Fixed"
146 )
147 simu.add_surfLoad(nodes_Load, [-surfLoad], ["y"])
148
149 simu.Solve()
150
151 simu.Save_Iter()
152
153 return simu
154
155 simu = Simulations.Mesh_Optim_ZZ1(DoSimu, folder, threshold, iterMax, 1 / 10)
156 PyVista.Plot_BoundaryConditions(simu).show()
157
158 # ----------------------------------------------
159 # Plot
160 # ----------------------------------------------
161 PyVista.Plot_Mesh(simu.mesh).show()
162 PyVista.Plot(simu, "ZZ1_e", nodeValues=False, nColors=11).show()
163
164 if plotProj:
165 simu.Set_Iter(0)
166 mesh0 = simu.mesh
167 u0 = np.reshape(simu.displacement, (mesh0.Nn, -1))
168
169 simu.Set_Iter(1)
170 mesh1 = simu.mesh
171
172 proj = Calc_projector(mesh0, mesh1)
173 uProj = np.zeros((mesh1.Nn, dim), dtype=float)
174 for d in range(dim):
175 uProj[:, d] = proj @ u0[:, d]
176
177 ax = Display.Plot_Result(
178 mesh0, np.linalg.norm(u0, axis=1), plotMesh=True, title="u0"
179 )
180 ax.plot(*mesh1.coord[:, :dim].T, ls="", marker="+", c="k", label="new nodes")
181 ax.legend()
182 Display.Plot_Result(
183 mesh1, np.linalg.norm(uProj, axis=1), plotMesh=True, title="uProj"
184 )
185
186 if makeParaview:
187 Paraview.Save_simu(simu, folder, nodeFields=["ZZ1_e"])
188
189 if makeMovie:
190
191 def func(plotter, n):
192 simu.Set_Iter(n)
193
194 PyVista.Plot_Mesh(simu, plotter=plotter)
195
196 zz1 = simu._Calc_ZZ1()[0]
197
198 plotter.add_title(f"ZZ1 = {zz1 * 100:.2f} %")
199
200 PyVista.Movie_func(func, len(simu.results), folder, "letterWeigher.gif")
201
202 Tic.Plot_History()
203 plt.show()
Total running time of the script: (0 minutes 9.391 seconds)