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


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

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

Gallery generated by Sphinx-Gallery