Note

Go to the end to download the full example code.

Mesh8#

Meshing of a grooved 3D part with calculation of element quality.

Mesh8
Element type: TETRA4
Ne = 12305, Nn = 3996
volume = 658.400


 13 import numpy as np
 14 import gmsh
 15
 16 from EasyFEA import Display, Mesher, ElemType, PyVista
 17 from EasyFEA.Geoms import Point, Circle, Points
 18
 19 if __name__ == "__main__":
 20     Display.Clear()
 21
 22     # ----------------------------------------------
 23     # Geom
 24     # ----------------------------------------------
 25     R = 10
 26     e = 2
 27     r = R - e
 28     h = R * 2 / 3
 29
 30     meshSize = e / 3
 31
 32     center = Point()
 33
 34     circle_ext = Circle(center, R * 2, meshSize)
 35     circle_int = Circle(center, r * 2, meshSize)
 36
 37     useFillet = not circle_int.isFilled
 38     addCylinder = True
 39     addBox = True
 40     addRevolve = True
 41
 42     # ----------------------------------------------
 43     # Mesh
 44     # ----------------------------------------------
 45
 46     mesher = Mesher(False, False)
 47     # dim, elemType = 2, ElemType.TRI3
 48     dim, elemType = 3, ElemType.TETRA4
 49
 50     factory = mesher._factory
 51
 52     surfaces = mesher._Surfaces(circle_ext, [circle_int])[0]
 53     mesher._Extrude(surfaces, [0, 0, R], elemType)
 54     vol1 = factory.getEntities(3)
 55
 56     mesher._Synchronize()
 57
 58     if useFillet:
 59         surfs = gmsh.model.getBoundary(vol1)
 60         lines = gmsh.model.getBoundary(surfs, False)
 61         gmsh.model.occ.fillet([vol1[0][1]], [abs(line[1]) for line in lines], [e / 3])
 62
 63     if addCylinder:
 64         cylinder = factory.addCylinder(R + e, 0, h, -2 * R - e, 0, 0, e)
 65         factory.cut(vol1, [(3, cylinder)])
 66
 67     if addBox:
 68         box = factory.addBox(-R / 4, -R - e, R / 2, R / 2, 2 * R + e, R / 2)
 69         factory.cut(factory.getEntities(3), [(3, box)])
 70
 71     if addRevolve:
 72         p1 = Point(R - e / 2, 0, e, r=e / 4)
 73         p2 = Point(R, 0, e)
 74         p3 = Point(R, 0, e * 4)
 75         p4 = Point(R - e / 2, 0, e * 4, r=e / 4)
 76         contour = Points([p1, p2, p3, p4])
 77         surf = mesher._Surfaces(contour, [])[0][0]
 78
 79         rev1 = factory.revolve([(2, surf)], 0, 0, 0, 0, 0, R, np.pi)
 80         rev2 = factory.revolve([(2, surf)], 0, 0, 0, 0, 0, R, -np.pi)
 81         factory.cut(factory.getEntities(3), [rev1[1], rev2[1]])
 82
 83     mesher._Synchronize()
 84
 85     if meshSize > 0:
 86         mesher.Set_meshSize(meshSize)
 87
 88     mesher._Set_PhysicalGroups(
 89         setPoints=False, setLines=True, setSurfaces=True, setVolumes=False
 90     )
 91
 92     mesher._Mesh_Generate(dim, elemType)
 93
 94     mesh = mesher._Mesh_Get_Mesh()
 95
 96     # ----------------------------------------------
 97     # Plot
 98     # ----------------------------------------------
 99
100     print(mesh)
101
102     if dim == 3:
103         print(f"volume = {mesh.volume:.3f}")
104
105     plotter = PyVista._Plotter(shape=(1, 2))
106
107     PyVista.Plot_Mesh(mesh, plotter=plotter)
108
109     plotter.subplot(0, 1)
110     plotter.add_title("aspect ratio")
111     qual = mesh.Get_Quality()
112     PyVista.Plot_Elements(mesh, dimElem=1, plotter=plotter, color="k")
113     PyVista.Plot(
114         mesh,
115         qual,
116         nodeValues=False,
117         cmap="viridis",
118         clim=(0, 1),
119         plotMesh=True,
120         plotter=plotter,
121     ).show()

Total running time of the script: (0 minutes 1.073 seconds)

Gallery generated by Sphinx-Gallery