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.9.0/examples/LinearizedElasticity/MeshOptim0.py:98: 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))
13 import matplotlib.pyplot as plt
14 import numpy as np
15
16 from EasyFEA import Display, Folder, ElemType, Mesh, PyVista
17 from EasyFEA.Geoms import Point, Circle, CircleArc, Contour
18 from EasyFEA.FEM import Mesh_Optim
19
20 if __name__ == "__main__":
21 Display.Clear()
22
23 # ----------------------------------------------
24 # Configuration
25 # ----------------------------------------------
26
27 # outputs
28 folder = Folder.Results_Dir()
29
30 # geom
31 D = 1
32 r = D * 1 / 4
33 e = 0.1
34 b = D * 0.1
35
36 # criteria
37 criteria = "aspect"
38 quality = 0.8 # lower bound of the target quality
39 ratio = 0.6 # the ratio of mesh elements that must satisfy the target
40 iterMax = 20 # Maximum number of iterations
41 coef = 1 / 2 # Scaling coefficient for the optimization process
42
43 # ----------------------------------------------
44 # Mesh
45 # ----------------------------------------------
46 elemType = ElemType.TRI3
47 mS = e / 3 * 10
48
49 # face
50 circle = Circle(Point(), D, mS)
51
52 # eyes
53 theta = 45 * np.pi / 180
54 sin = np.sin(theta)
55 cos = np.cos(theta)
56 eye1 = Circle(Point(-r * sin, r * cos), e, mS)
57 eye2 = Circle(Point(r * sin, r * cos), e, mS)
58
59 # happy smile
60 s = (D * 0.6) / 2
61 s1 = CircleArc(Point(-s), Point(s), Point(), coef=-1, meshSize=mS)
62 s2 = CircleArc(Point(s), Point(s - e), Point(s - e / 2), coef=-1, meshSize=mS)
63 s3 = CircleArc(Point(s - e), Point(-s + e), Point(), meshSize=mS)
64 s4 = CircleArc(Point(-s + e), Point(-s), Point(-s + e / 2), coef=-1, meshSize=mS)
65 happy = Contour([s1, s2, s3, s4])
66
67 inclusions = [eye1, eye2, happy]
68
69 def DoMesh(refineGeom=None) -> Mesh:
70 """Function used for mesh generation"""
71 return circle.Mesh_2D(inclusions, elemType, [], [refineGeom])
72
73 # Construct the initial mesh
74 mesh = DoMesh()
75 PyVista.Plot_Mesh(mesh).show()
76
77 mesh, ratio = Mesh_Optim(DoMesh, folder, criteria, quality, ratio)
78
79 # ----------------------------------------------
80 # Plot
81 # ----------------------------------------------
82
83 qual_e = mesh.Get_Quality(criteria, False)
84
85 PyVista.Plot_Mesh(mesh).show()
86 PyVista.Plot(
87 mesh,
88 qual_e,
89 nodeValues=False,
90 plotMesh=True,
91 clim=(0, quality),
92 cmap="viridis",
93 ).show()
94
95 axHist = Display.Init_Axes()
96
97 axHist.hist(qual_e, 11, (0, 1))
98 axHist.set_xlabel("quality")
99 axHist.set_ylabel("elements")
100 axHist.set_title(f"ratio = {ratio * 100:.3f} %")
101 axHist.vlines([quality], [0], [ratio * mesh.Ne], color="red")
102
103 plt.show()
Total running time of the script: (0 minutes 0.824 seconds)