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