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