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