Note
Go to the end to download the full example code.
MeshOptim0#
Optimization of a happy mesh with quality criteria.
ratio = 37.500 %
ratio = 57.143 %
ratio = 66.809 %
/home/docs/checkouts/readthedocs.org/user_builds/easyfea/checkouts/v1.5.4/examples/Meshes/MeshOptim0.py:94: MatplotlibDeprecationWarning: Passing the range parameter of hist() positionally is deprecated since Matplotlib 3.10; the parameter will become keyword-only in 3.12.
axHist.hist(qual_e, 11, (0, 1))
12 from EasyFEA import Display, Folder, plt, np, Mesher, ElemType, Mesh, PyVista
13 from EasyFEA.Geoms import Point, Circle, CircleArc, Contour
14 from EasyFEA.fem import Mesh_Optim
15
16 if __name__ == "__main__":
17 Display.Clear()
18
19 # ----------------------------------------------
20 # Configuration
21 # ----------------------------------------------
22
23 # outputs
24 folder = Folder.Join(Folder.RESULTS_DIR, "Meshes", "Optim2D", mkdir=True)
25
26 # geom
27 D = 1
28 r = D * 1 / 4
29 e = 0.1
30 b = D * 0.1
31
32 # criteria
33 criteria = "aspect"
34 quality = 0.8 # lower bound of the target quality
35 ratio = 0.6 # the ratio of mesh elements that must satisfy the target
36 iterMax = 20 # Maximum number of iterations
37 coef = 1 / 2 # Scaling coefficient for the optimization process
38
39 # ----------------------------------------------
40 # Mesh
41 # ----------------------------------------------
42 elemType = ElemType.TRI3
43 mS = e / 3 * 10
44
45 # face
46 circle = Circle(Point(), D, mS)
47
48 # eyes
49 theta = 45 * np.pi / 180
50 sin = np.sin(theta)
51 cos = np.cos(theta)
52 eye1 = Circle(Point(-r * sin, r * cos), e, mS)
53 eye2 = Circle(Point(r * sin, r * cos), e, mS)
54
55 # happy smile
56 s = (D * 0.6) / 2
57 s1 = CircleArc(Point(-s), Point(s), Point(), coef=-1, meshSize=mS)
58 s2 = CircleArc(Point(s), Point(s - e), Point(s - e / 2), coef=-1, meshSize=mS)
59 s3 = CircleArc(Point(s - e), Point(-s + e), Point(), meshSize=mS)
60 s4 = CircleArc(Point(-s + e), Point(-s), Point(-s + e / 2), coef=-1, meshSize=mS)
61 happy = Contour([s1, s2, s3, s4])
62
63 inclusions = [eye1, eye2, happy]
64
65 def DoMesh(refineGeom=None) -> Mesh:
66 """Function used for mesh generation"""
67 return Mesher().Mesh_2D(circle, inclusions, elemType, [], [refineGeom])
68
69 # Construct the initial mesh
70 mesh = DoMesh()
71 PyVista.Plot_Mesh(mesh).show()
72
73 mesh, ratio = Mesh_Optim(DoMesh, folder, criteria, quality, ratio)
74
75 # ----------------------------------------------
76 # Plot
77 # ----------------------------------------------
78
79 qual_e = mesh.Get_Quality(criteria, False)
80
81 PyVista.Plot_Mesh(mesh).show()
82 PyVista.Plot(
83 mesh,
84 qual_e,
85 nodeValues=False,
86 plotMesh=True,
87 clim=(0, quality),
88 cmap="viridis",
89 ).show()
90
91 axHist = Display.Init_Axes()
92
93 axHist.hist(qual_e, 11, (0, 1))
94 axHist.set_xlabel("quality")
95 axHist.set_ylabel("elements")
96 axHist.set_title(f"ratio = {ratio * 100:.3f} %")
97 axHist.vlines([quality], [0], [ratio * mesh.Ne], color="red")
98
99 plt.show()
Total running time of the script: (0 minutes 0.823 seconds)