PeriDEM 0.2.0
PeriDEM -- Peridynamics-based high-fidelity model for granular media
Loading...
Searching...
No Matches
gmsh_particles.py
Go to the documentation of this file.
1import os
2import sys
3import numpy as np
4import gmsh
5import math
6
7from util import *
8from geom_util import *
9
10
11
12def circle_mesh_symmetric(xc = [0., 0., 0.], r = 1., h = 0.1, filename = 'mesh', vtk_out = False, symmetric_mesh = True):
13 """
14 Create a circular mesh, either symmetric or full.
15 xc - center point for symmetric mesh, or bottom-left corner for non-symmetric mesh
16 r - radius of the circle
17 h - mesh size
18 symmetric_mesh - if True, creates 1/4 mesh and mirrors it. If False, creates full circle
19 """
20
21 gmsh.initialize()
22 gmsh.option.setNumber("Mesh.MshFileVersion", 2.2)
23
24 if symmetric_mesh:
25 # we first assume circle center is at origin and then translate the symmetric mesh to the desired location
26 xc_mesh = [0., 0., 0.]
27 p1 = gmsh.model.geo.addPoint(xc_mesh[0], xc_mesh[1], xc_mesh[2], h)
28 p2 = gmsh.model.geo.addPoint(xc_mesh[0] + r, xc_mesh[1], xc_mesh[2], h)
29 p3 = gmsh.model.geo.addPoint(xc_mesh[0], xc_mesh[1] + r, xc_mesh[2], h)
30
31 l1 = gmsh.model.geo.addCircleArc(p2, p1, p3)
32 l2 = gmsh.model.geo.addLine(p1, p2)
33 l3 = gmsh.model.geo.addLine(p3, p1)
34
35 c1 = gmsh.model.geo.addCurveLoop([l2, l1, l3])
36
37 p1 = gmsh.model.geo.addPlaneSurface([c1])
38
39 #gmsh.model.occ.addBox(0,0,0, 1,0.5,0.5)
40 gmsh.model.geo.synchronize()
41 # gmsh.model.mesh.setSize(gmsh.model.getEntities(0), 0.1)
42 # gmsh.model.mesh.setSize([(0, 2)], 0.01)
43 gmsh.model.mesh.generate(3)
44 # gmsh.write('mesh.vtk')
45
46 # get the mesh data
47 m = get_gmsh_entities()
48
49 gmsh_transform(m, 1000, 1000000, 1000000, -1, 1, 1) # x-axis
50 gmsh_transform(m, 2000, 2000000, 2000000, 1, -1, 1) # y-axis
51 gmsh_transform(m, 3000, 3000000, 3000000, -1, -1, 1) # z-axis
52
53 # remove the duplicate nodes that will have been created on the internal
54 # boundaries
55 gmsh.model.mesh.removeDuplicateNodes()
56
57 # translate the mesh to the specified center coordinates
58 gmsh_translate(xc)
59
60 else:
61 # Create circle using built-in gmsh circle
62 c = gmsh.model.occ.addCircle(xc[0], xc[1], xc[2], r)
63
64 # Create curve loop from circle
65 cl = gmsh.model.occ.addCurveLoop([c])
66
67 # Create surface from curve loop
68 s = gmsh.model.occ.addPlaneSurface([cl])
69
70 # Add center point
71 p = gmsh.model.occ.addPoint(xc[0], xc[1], xc[2], h)
72
73 # Synchronize
74 gmsh.model.occ.synchronize()
75
76 # Embed center point in surface
77 gmsh.model.mesh.embed(0, [p], 2, s)
78
79 # generate mesh
80 gmsh.model.mesh.generate(3)
81
82 # Check for hanging nodes before writing
83 check_hanging_nodes()
84
85 # write
86 gmsh.write(filename + '.msh')
87 if vtk_out:
88 gmsh.write(filename + '.vtk')
89
90 gmsh.finalize()
91
92def ellipse_mesh_symmetric(xc = [0., 0., 0.], rx = 1., ry = 0.5, h = 0.1, filename = 'mesh', vtk_out = False, symmetric_mesh = True):
93 """
94 Create an elliptical mesh, either symmetric or full.
95 xc - center point for symmetric mesh, or bottom-left corner for non-symmetric mesh
96 rx - radius in x direction
97 ry - radius in y direction
98 h - mesh size
99 symmetric_mesh - if True, creates 1/4 mesh and mirrors it. If False, creates full ellipse
100 """
101
102 gmsh.initialize()
103 gmsh.option.setNumber("Mesh.MshFileVersion", 2.2)
104
105 if symmetric_mesh:
106 # we first assume ellipse center is at origin and then translate the symmetric mesh to the desired location
107 xc_mesh = [0., 0., 0.]
108 p1 = gmsh.model.geo.addPoint(xc_mesh[0], xc_mesh[1], xc_mesh[2], h) # Center
109 p2 = gmsh.model.geo.addPoint(xc_mesh[0] + rx, xc_mesh[1], xc_mesh[2], h) # Right
110 p3 = gmsh.model.geo.addPoint(xc_mesh[0], xc_mesh[1] + ry, xc_mesh[2], h) # Top
111
112 # Create elliptical arc using 3 points
113 l1 = gmsh.model.geo.addEllipseArc(p2, p1, p2, p3) # Right to top
114 l2 = gmsh.model.geo.addLine(p1, p2) # Center to right
115 l3 = gmsh.model.geo.addLine(p3, p1) # Top to center
116
117 # Create curve loop and surface
118 c1 = gmsh.model.geo.addCurveLoop([l2, l1, l3])
119 s1 = gmsh.model.geo.addPlaneSurface([c1])
120
121 gmsh.model.geo.synchronize()
122 gmsh.model.mesh.generate(3)
123
124 # get the mesh data
125 m = get_gmsh_entities()
126
127 # Mirror the mesh to create full ellipse
128 gmsh_transform(m, 1000, 1000000, 1000000, -1, 1, 1) # x-axis
129 gmsh_transform(m, 2000, 2000000, 2000000, 1, -1, 1) # y-axis
130 gmsh_transform(m, 3000, 3000000, 3000000, -1, -1, 1) # z-axis
131
132 # Remove duplicate nodes on internal boundaries
133 gmsh.model.mesh.removeDuplicateNodes()
134
135 # Translate to specified center coordinates
136 gmsh_translate(xc)
137
138 else:
139 # Create points for full ellipse
140 p1 = gmsh.model.geo.addPoint(xc[0], xc[1], xc[2], h) # Center
141 p2 = gmsh.model.geo.addPoint(xc[0] + rx, xc[1], xc[2], h) # Right
142 p3 = gmsh.model.geo.addPoint(xc[0], xc[1] + ry, xc[2], h) # Top
143 p4 = gmsh.model.geo.addPoint(xc[0] - rx, xc[1], xc[2], h) # Left
144 p5 = gmsh.model.geo.addPoint(xc[0], xc[1] - ry, xc[2], h) # Bottom
145
146 # Create elliptical arcs
147 l1 = gmsh.model.geo.addEllipseArc(p2, p1, p2, p3) # Right to top
148 l2 = gmsh.model.geo.addEllipseArc(p3, p1, p3, p4) # Top to left
149 l3 = gmsh.model.geo.addEllipseArc(p4, p1, p4, p5) # Left to bottom
150 l4 = gmsh.model.geo.addEllipseArc(p5, p1, p5, p2) # Bottom to right
151
152 # Create curve loop and surface
153 c1 = gmsh.model.geo.addCurveLoop([l1, l2, l3, l4])
154 s1 = gmsh.model.geo.addPlaneSurface([c1])
155
156 # synchronize
157 gmsh.model.geo.synchronize()
158
159 # embed p0 (dim 0) in surface 1 (dim 2)
160 gmsh.model.mesh.embed(0, [p1], 2, s1)
161
162 # generate mesh
163 gmsh.model.mesh.generate(3)
164
165 # Check for hanging nodes before writing
166 check_hanging_nodes()
167
168 # Write output files
169 gmsh.write(filename + '.msh')
170 if vtk_out:
171 gmsh.write(filename + '.vtk')
172
173 gmsh.finalize()
174
175def sphere_mesh_symmetric(xc = [0., 0., 0.], r = 1., h = 0.1, filename = 'mesh', vtk_out = False, symmetric_mesh = True):
176 """
177 Create a spherical mesh, either symmetric or full.
178 xc - center point coordinates [x, y, z]
179 r - radius of the sphere
180 h - mesh size
181 symmetric_mesh - if True, creates 1/8 mesh and mirrors it. If False, creates full sphere
182 """
183 gmsh.initialize()
184 gmsh.option.setNumber("Mesh.MshFileVersion", 2.2)
185
186 if symmetric_mesh:
187 # Create 1/8th of sphere at origin, then mirror it
188 xc_mesh = [0., 0., 0.]
189
190 # Create first octant of sphere using OpenCASCADE
191 sphere = gmsh.model.occ.addSphere(xc_mesh[0], xc_mesh[1], xc_mesh[2], r)
192
193 # cubes to cut the sphere
194 cut_left = gmsh.model.occ.addBox(xc_mesh[0]-r, xc_mesh[1]-r, xc_mesh[2]-r, r, 2*r, 2*r)
195 cut_back = gmsh.model.occ.addBox(xc_mesh[0]-r, xc_mesh[1]-r, xc_mesh[2]-r, 2*r, r, 2*r)
196 cut_bottom = gmsh.model.occ.addBox(xc_mesh[0]-r, xc_mesh[1]-r, xc_mesh[2]-r, 2*r, 2*r, r)
197
198 # Cut sphere to get first octant using 3 planes
199 gmsh.model.occ.cut([(3,sphere)], [(3,cut_left), (3,cut_back), (3,cut_bottom)])
200
201 gmsh.model.occ.synchronize()
202
203 # Set mesh size
204 gmsh.model.mesh.setSize(gmsh.model.getEntities(0), h)
205
206 # Generate 3D mesh
207 gmsh.model.mesh.generate(3)
208
209 # gmsh.fltk.run()
210
211 # Get mesh data for mirroring
212 m = get_gmsh_entities()
213
214 # Mirror mesh in x direction
215 gmsh_transform(m, 1000, 10000, 100000, -1, 1, 1)
216 # Mirror in y direction
217 gmsh_transform(m, 2000, 20000, 200000, 1, -1, 1)
218 # Mirror in z direction
219 gmsh_transform(m, 3000, 30000, 300000, 1, 1, -1)
220 # Mirror in xy plane
221 gmsh_transform(m, 4000, 40000, 400000, -1, -1, 1)
222 # Mirror in yz plane
223 gmsh_transform(m, 5000, 50000, 500000, 1, -1, -1)
224 # Mirror in xz plane
225 gmsh_transform(m, 6000, 60000, 600000, -1, 1, -1)
226 # Mirror in xyz
227 gmsh_transform(m, 7000, 70000, 700000, -1, -1, -1)
228
229 # Translate to specified center if not at origin
230 gmsh_translate(xc)
231
232 else:
233 # Create points for full sphere
234 # Create full sphere using OpenCASCADE
235 sphere = gmsh.model.occ.addSphere(xc[0], xc[1], xc[2], r)
236
237 # Synchronize the OpenCASCADE geometry
238 gmsh.model.occ.synchronize()
239
240 # Set mesh size
241 gmsh.model.mesh.setSize(gmsh.model.getEntities(0), h)
242
243 # Synchronize
244 gmsh.model.geo.synchronize()
245
246 # Embed center point in volume
247 p1 = gmsh.model.geo.addPoint(xc[0], xc[1], xc[2], h)
248 gmsh.model.mesh.embed(0, [p1], 3, sphere)
249
250 # Generate 3D mesh
251 gmsh.model.mesh.generate(3)
252
253 # Check for hanging nodes before writing
254 check_hanging_nodes()
255
256 # Write output files
257 gmsh.write(filename + '.msh')
258 if vtk_out:
259 gmsh.write(filename + '.vtk')
260
261 gmsh.finalize()
262
263def rectangle_mesh_symmetric(xc = [0., 0., 0.], Lx = 1., Ly = 1., h = 0.1, filename = 'mesh', vtk_out = False, symmetric_mesh = True):
264 """
265 Create a rectangular mesh, either symmetric or full.
266 xc - center point for symmetric mesh, or bottom-left corner for non-symmetric mesh
267 Lx - length in x direction
268 Ly - length in y direction
269 h - mesh size
270 symmetric_mesh - if True, creates 1/4 mesh and mirrors it. If False, creates full rectangle
271 """
272 # Initialize gmsh
273 gmsh.initialize()
274 gmsh.option.setNumber("Mesh.MshFileVersion", 2.2)
275
276 if symmetric_mesh:
277 # center at origin
278 xc_mesh = [0., 0., 0.]
279
280 # Create points for 1/4 rectangle
281 p1 = gmsh.model.geo.addPoint(xc_mesh[0], xc_mesh[1], xc_mesh[2], h)
282 p2 = gmsh.model.geo.addPoint(xc_mesh[0]+0.5*Lx, xc_mesh[1], xc_mesh[2], h)
283 p3 = gmsh.model.geo.addPoint(xc_mesh[0]+0.5*Lx, xc_mesh[1]+0.5*Ly, xc_mesh[2], h)
284 p4 = gmsh.model.geo.addPoint(xc_mesh[0], xc_mesh[1]+0.5*Ly, xc_mesh[2], h)
285
286 # Create lines
287 l1 = gmsh.model.geo.addLine(p1, p2)
288 l2 = gmsh.model.geo.addLine(p2, p3)
289 l3 = gmsh.model.geo.addLine(p3, p4)
290 l4 = gmsh.model.geo.addLine(p4, p1)
291
292 # Create curve loop and surface
293 c1 = gmsh.model.geo.addCurveLoop([l1, l2, l3, l4])
294 s1 = gmsh.model.geo.addPlaneSurface([c1])
295
296 # Synchronize and generate mesh
297 gmsh.model.geo.synchronize()
298 gmsh.model.mesh.generate(3)
299
300 # Get the mesh data for mirroring
301 m = get_gmsh_entities()
302
303 # Mirror the mesh in all quadrants
304 gmsh_transform(m, 1000, 1000000, 1000000, -1, 1, 1)
305 gmsh_transform(m, 2000, 2000000, 2000000, 1, -1, 1)
306 gmsh_transform(m, 3000, 3000000, 3000000, -1, -1, 1)
307
308 # Remove duplicate nodes
309 gmsh.model.mesh.removeDuplicateNodes()
310
311 # translate the mesh to the specified center coordinates
312 gmsh_translate(xc)
313
314 else:
315 # Create points for full rectangle
316 p0 = gmsh.model.geo.addPoint(xc[0], xc[1], xc[2], h) # Center
317 p1 = gmsh.model.geo.addPoint(xc[0] - 0.5*Lx, xc[1] - 0.5*Ly, xc[2], h) # Bottom left
318 p2 = gmsh.model.geo.addPoint(xc[0] + 0.5*Lx, xc[1] - 0.5*Ly, xc[2], h) # Bottom right
319 p3 = gmsh.model.geo.addPoint(xc[0] + 0.5*Lx, xc[1] + 0.5*Ly, xc[2], h) # Top right
320 p4 = gmsh.model.geo.addPoint(xc[0] - 0.5*Lx, xc[1] + 0.5*Ly, xc[2], h) # Top left
321
322 # Create lines
323 l1 = gmsh.model.geo.addLine(p1, p2) # Bottom
324 l2 = gmsh.model.geo.addLine(p2, p3) # Right
325 l3 = gmsh.model.geo.addLine(p3, p4) # Top
326 l4 = gmsh.model.geo.addLine(p4, p1) # Left
327
328 # Create curve loop and surface
329 c1 = gmsh.model.geo.addCurveLoop([l1, l2, l3, l4])
330 s1 = gmsh.model.geo.addPlaneSurface([c1])
331
332 # synchronize
333 gmsh.model.geo.synchronize()
334
335 # embed p0 (dim 0) in surface 1 (dim 2)
336 gmsh.model.mesh.embed(0, [p0], 2, s1)
337
338 # generate mesh
339 gmsh.model.mesh.generate(3)
340
341 # Check for hanging nodes before writing
342 check_hanging_nodes()
343
344 # write
345 gmsh.write(filename + '.msh')
346 if vtk_out:
347 gmsh.write(filename + '.vtk')
348
349 gmsh.finalize()
350
351def hexagon_mesh_symmetric(xc = [0., 0., 0.], r = 1., h = 0.1, filename = 'mesh', vtk_out = False, symmetric_mesh = True):
352 """
353 Create a hexagon mesh, either symmetric or full.
354 xc - center point for symmetric mesh, or bottom-left corner for non-symmetric mesh
355 r - radius of the hexagon
356 h - mesh size
357 symmetric_mesh - if True, creates 1/6 mesh and mirrors it. If False, creates full hexagon
358 """
359 gmsh.initialize()
360 gmsh.option.setNumber("Mesh.MshFileVersion", 2.2)
361
362 if symmetric_mesh:
363 # Create 1/6 hexagon mesh at origin first, then rotate and translate
364 xc_mesh = [0., 0., 0.]
365 p1 = gmsh.model.geo.addPoint(xc_mesh[0], xc_mesh[1], xc_mesh[2], h) # Center (origin)
366 p2 = gmsh.model.geo.addPoint(xc_mesh[0] + r, xc_mesh[1], xc_mesh[2], h) # Point on x-axis (p1)
367 p3 = gmsh.model.geo.addPoint(xc_mesh[0] + r*np.cos(np.pi/3), xc_mesh[1] + r*np.sin(np.pi/3), xc_mesh[2], h) # p2
368
369 # Create lines for 1/6th of hexagon
370 l1 = gmsh.model.geo.addLine(p1, p2) # Center to p1
371 l2 = gmsh.model.geo.addLine(p2, p3) # p1 to p2
372 l3 = gmsh.model.geo.addLine(p3, p1) # p2 to center
373
374 # Create curve loop and surface
375 c1 = gmsh.model.geo.addCurveLoop([l1, l2, l3])
376 s1 = gmsh.model.geo.addPlaneSurface([c1])
377
378 gmsh.model.geo.synchronize()
379 gmsh.model.mesh.generate(3)
380
381 # Get mesh data for rotation
382 m = get_gmsh_entities()
383
384 # Rotate the mesh 5 times by 60 degrees
385 for i in range(5):
386 angle = (i + 1) * np.pi/3
387 gmsh_transform_general(m, 1000*(i+1), 1000000*(i+1), 1000000*(i+1), angle)
388
389 # Remove duplicate nodes
390 gmsh.model.mesh.removeDuplicateNodes()
391
392 # Translate to specified center
393 gmsh_translate(xc)
394
395 else:
396 # For non-symmetric mesh, create full hexagon directly
397 points = []
398 for i in range(6):
399 angle = i * np.pi/3
400 px = xc[0] + r * np.cos(angle)
401 py = xc[1] + r * np.sin(angle)
402 points.append(gmsh.model.geo.addPoint(px, py, xc[2], h))
403
404 points.append(gmsh.model.geo.addPoint(xc[0], xc[1], xc[2], h))
405
406 # Create lines connecting points
407 lines = []
408 for i in range(5):
409 lines.append(gmsh.model.geo.addLine(points[i], points[i+1]))
410 lines.append(gmsh.model.geo.addLine(points[5], points[0]))
411
412 # Create curve loop and surface
413 c1 = gmsh.model.geo.addCurveLoop(lines)
414 s1 = gmsh.model.geo.addPlaneSurface([c1])
415
416 # synchronize
417 gmsh.model.geo.synchronize()
418
419 # embed p1 (dim 0) in surface 1 (dim 2)
420 gmsh.model.mesh.embed(0, [points[-1]], 2, s1)
421
422 # generate mesh
423 gmsh.model.mesh.generate(3)
424
425 # Check for hanging nodes before writing
426 check_hanging_nodes()
427
428 # Write output files
429 gmsh.write(filename + '.msh')
430 if vtk_out:
431 gmsh.write(filename + '.vtk')
432
433 gmsh.finalize()
434
435def drum2d_mesh_symmetric(xc = [0., 0., 0.], r = 1., width = 1., h = 0.1, filename = 'mesh', vtk_out = False, symmetric_mesh = True):
436 """
437 Create a drum2d mesh, either symmetric or full.
438 xc - center point for symmetric mesh, or bottom-left corner for non-symmetric mesh
439 r - radius of the drum2d
440 width - width of the drum2d
441 h - mesh size
442 symmetric_mesh - if True, creates 1/4 mesh and mirrors it. If False, creates full drum2d
443 """
444 gmsh.initialize()
445 gmsh.option.setNumber("Mesh.MshFileVersion", 2.2)
446
447 if symmetric_mesh:
448 # Create 1/4 drum2d mesh at origin first, then rotate and translate
449 xc_mesh = [0., 0., 0.]
450 x_drum = get_ref_drum_points(xc_mesh, r, width)
451 v1, v2, v3 = x_drum[0], x_drum[1], x_drum[2]
452 v23 = [0.5*(v2[0] + v3[0]), 0.5*(v2[1] + v3[1]), 0.5*(v2[2] + v3[2])]
453
454 # 1/4th mesh is a polygon with 4 points (o, v1, v2, (v2+v3)/2)
455 # Create points for 1/4th mesh
456 p1 = gmsh.model.geo.addPoint(xc_mesh[0], xc_mesh[1], xc_mesh[2], h) # Center point
457 p2 = gmsh.model.geo.addPoint(v1[0], v1[1], v1[2], h)
458 p3 = gmsh.model.geo.addPoint(v2[0], v2[1], v2[2], h)
459 p4 = gmsh.model.geo.addPoint(v23[0], v23[1], v23[2], h)
460
461 # Create lines
462 lines = []
463 points = [p1, p2, p3, p4, p1] # Add p1 again at end to complete loop
464 for i in range(4):
465 lines.append(gmsh.model.geo.addLine(points[i], points[i+1]))
466
467 # Create curve loop and surface
468 c1 = gmsh.model.geo.addCurveLoop(lines)
469 s1 = gmsh.model.geo.addPlaneSurface([c1])
470
471 gmsh.model.geo.synchronize()
472 gmsh.model.mesh.generate(3)
473
474 # Get mesh data for mirroring
475 m = get_gmsh_entities()
476
477 # Mirror mesh to create full drum2d
478 gmsh_transform(m, 1000, 1000000, 1000000, -1, 1, 1) # x-axis
479 gmsh_transform(m, 2000, 2000000, 2000000, 1, -1, 1) # y-axis
480 gmsh_transform(m, 3000, 3000000, 3000000, -1, -1, 1) # z-axis
481
482 # Remove duplicate nodes on internal boundaries
483 gmsh.model.mesh.removeDuplicateNodes()
484
485 # Translate to specified center
486 gmsh_translate(xc)
487
488 else:
489 # For non-symmetric mesh, create full drum2d directly
490 x_drum = get_ref_drum_points(xc, r, width)
491 # add points
492 points = []
493 points.append(gmsh.model.geo.addPoint(xc[0], xc[1], xc[2], h))
494 for i in range(len(x_drum)):
495 points.append(gmsh.model.geo.addPoint(x_drum[i][0], x_drum[i][1], x_drum[i][2], h))
496
497 # Create lines connecting vertices
498 lines = []
499 points.append(points[1]) # to close the loop
500 for i in range(len(x_drum)):
501 lines.append(gmsh.model.geo.addLine(points[i+1], points[i+2]))
502
503 # Create curve loop and surface
504 c1 = gmsh.model.geo.addCurveLoop(lines)
505 s1 = gmsh.model.geo.addPlaneSurface([c1])
506
507 # synchronize
508 gmsh.model.geo.synchronize()
509
510 # embed p1 (dim 0) in surface 1 (dim 2)
511 gmsh.model.mesh.embed(0, [points[0]], 2, s1)
512
513 # generate mesh
514 gmsh.model.mesh.generate(3)
515
516 # Check for hanging nodes before writing
517 check_hanging_nodes()
518
519 # Write output files
520 gmsh.write(filename + '.msh')
521 if vtk_out:
522 gmsh.write(filename + '.vtk')
523
524 gmsh.finalize()
525
526def polygon_mesh_symmetric(points, theta, xc=[0., 0., 0.], h=0.1, filename='mesh', vtk_out=False, symmetric_mesh = True):
527 """
528 Create a symmetric mesh by rotating a polygon segment n times, where n = 360/theta.
529
530 Parameters:
531 points - List of [x,y,z] coordinates defining the polygon segment. First point must be [0,0,0]
532 and first two edges must form angle theta at origin
533 theta - Angle in radians between first two edges at origin. Must divide 2pi evenly.
534 xc - Center point for final translation
535 h - Mesh size
536 filename - Output filename
537 vtk_out - Whether to output VTK file
538 """
539 # Validate inputs
540 if not points or len(points) < 3:
541 raise ValueError("Need at least 3 points to define a polygon segment")
542
543 if not np.allclose(points[0], [0., 0., 0.]):
544 raise ValueError("First point must be at origin [0,0,0]")
545
546 # Check if theta divides 360 evenly
547 n = int(round(2*np.pi/theta))
548 if not np.isclose(2*np.pi, n * theta):
549 raise ValueError(f"Angle {theta} degrees must divide 360 evenly")
550
551 gmsh.initialize()
552 gmsh.option.setNumber("Mesh.MshFileVersion", 2.2)
553
554 # Create points for the polygon segment
555 gmsh_points = []
556 for pt in points:
557 gmsh_points.append(gmsh.model.geo.addPoint(pt[0], pt[1], pt[2], h))
558
559 # Create lines connecting points
560 lines = []
561 for i in range(len(gmsh_points)-1):
562 lines.append(gmsh.model.geo.addLine(gmsh_points[i], gmsh_points[i+1]))
563 # Close the polygon segment
564 lines.append(gmsh.model.geo.addLine(gmsh_points[-1], gmsh_points[0]))
565
566 # Create curve loop and surface
567 c1 = gmsh.model.geo.addCurveLoop(lines)
568 s1 = gmsh.model.geo.addPlaneSurface([c1])
569
570 # Synchronize and generate initial mesh
571 gmsh.model.geo.synchronize()
572 gmsh.model.mesh.generate(3)
573
574 # Get mesh data for rotation
575 m = get_gmsh_entities()
576
577 # Rotate n-1 times by theta
578 for i in range(1, n):
579 angle = i * theta # Convert to radians
580 gmsh_transform_general(m, 1000*i, 1000000*i, 1000000*i, angle)
581
582 # Remove duplicate nodes
583 gmsh.model.mesh.removeDuplicateNodes()
584
585 # Translate to final position if needed
586 gmsh_translate(xc)
587
588 # Check for hanging nodes before writing
589 check_hanging_nodes()
590
591 # Write output files
592 gmsh.write(filename + '.msh')
593 if vtk_out:
594 gmsh.write(filename + '.vtk')
595
596 gmsh.finalize()
597
598def cylindrical2d_wall_mesh(center=[0., 0., 0.], outer_radius=1.0, inner_radius=0.8,
599 bar_width=0.2, bar_length=0.3, h=0.1, filename='mesh', vtk_out=False):
600 """
601 Create a cylindrical wall mesh with bars.
602
603 Parameters:
604 center - Center point coordinates [x, y, z]
605 outer_radius - Outer radius of the wall
606 inner_radius - Inner radius of the wall
607 bar_width - Width of the bars
608 bar_length - Length of the bars
609 h - Mesh size
610 filename - Output filename
611 vtk_out - Whether to output VTK file
612 """
613 gmsh.initialize()
614 gmsh.option.setNumber("Mesh.MshFileVersion", 2.2)
615
616 # Create points
617 # Center point
618 p1 = gmsh.model.geo.addPoint(center[0], center[1], center[2], h)
619
620 # Outer circle points
621 p2 = gmsh.model.geo.addPoint(center[0] + outer_radius, center[1], center[2], h)
622 p3 = gmsh.model.geo.addPoint(center[0], center[1] + outer_radius, center[2], h)
623 p4 = gmsh.model.geo.addPoint(center[0] - outer_radius, center[1], center[2], h)
624 p5 = gmsh.model.geo.addPoint(center[0], center[1] - outer_radius, center[2], h)
625
626 # Inner points with bar
627 p6 = gmsh.model.geo.addPoint(center[0] + inner_radius, center[1] + 0.5*bar_width, center[2], h)
628 p7 = gmsh.model.geo.addPoint(center[0] + inner_radius, center[1] - 0.5*bar_width, center[2], h)
629 p8 = gmsh.model.geo.addPoint(center[0], center[1] + inner_radius, center[2], h)
630 p9 = gmsh.model.geo.addPoint(center[0] - inner_radius, center[1], center[2], h)
631 p10 = gmsh.model.geo.addPoint(center[0], center[1] - inner_radius, center[2], h)
632
633 # Bar end points
634 p11 = gmsh.model.geo.addPoint(center[0] + inner_radius - bar_length, center[1] + 0.5*bar_width, center[2], h)
635 p12 = gmsh.model.geo.addPoint(center[0] + inner_radius - bar_length, center[1] - 0.5*bar_width, center[2], h)
636
637 # Create circular arcs for outer circle
638 c1 = gmsh.model.geo.addCircleArc(p2, p1, p3)
639 c2 = gmsh.model.geo.addCircleArc(p3, p1, p4)
640 c3 = gmsh.model.geo.addCircleArc(p4, p1, p5)
641 c4 = gmsh.model.geo.addCircleArc(p5, p1, p2)
642
643 # Create circular arcs for inner circle
644 c5 = gmsh.model.geo.addCircleArc(p6, p1, p8)
645 c6 = gmsh.model.geo.addCircleArc(p8, p1, p9)
646 c7 = gmsh.model.geo.addCircleArc(p9, p1, p10)
647 c8 = gmsh.model.geo.addCircleArc(p10, p1, p7)
648
649 # Create lines for the bar
650 l1 = gmsh.model.geo.addLine(p6, p11)
651 l2 = gmsh.model.geo.addLine(p11, p12)
652 l3 = gmsh.model.geo.addLine(p12, p7)
653
654 # Create curve loops
655 outer_loop = gmsh.model.geo.addCurveLoop([c1, c2, c3, c4])
656 inner_loop = gmsh.model.geo.addCurveLoop([c5, c6, c7, c8, -l3, -l2, -l1])
657
658 # Create surface with hole
659 s1 = gmsh.model.geo.addPlaneSurface([outer_loop, inner_loop])
660
661 # Synchronize geometry before adding physical groups
662 gmsh.model.geo.synchronize()
663
664 # Add physical group for the surface
665 gmsh.model.addPhysicalGroup(2, [s1], 1)
666
667 # Generate mesh
668 gmsh.model.mesh.generate(2)
669
670 # Check for hanging nodes before writing
671 check_hanging_nodes()
672
673 # Write output files
674 gmsh.write(filename + '.msh')
675 if vtk_out:
676 gmsh.write(filename + '.vtk')
677
678 gmsh.finalize()
679
680def triangle_mesh_symmetric(xc=[0., 0., 0.], r=1., h=0.1, filename='mesh', vtk_out=False, symmetric_mesh=True):
681 """
682 Create a triangle mesh, either symmetric or full.
683
684 Parameters:
685 xc - center point coordinates [x, y, z]
686 r - radius of the circumscribed circle
687 h - mesh size
688 filename - output filename
689 vtk_out - whether to output VTK file
690 symmetric_mesh - if True, creates 1/3 mesh and rotates it twice. If False, creates full triangle
691 """
692 gmsh.initialize()
693 gmsh.option.setNumber("Mesh.MshFileVersion", 2.2)
694
695 if symmetric_mesh:
696 # Create 1/3 triangle mesh at origin first, then rotate and translate
697 xc_mesh = [0., 0., 0.]
698 # Get reference triangle points
699 x_tri = get_ref_triangle_points(xc_mesh, r)
700 v1, v2, v3 = x_tri[0], x_tri[1], x_tri[2]
701
702 # Create points for 1/3rd mesh (triangle from center to first two vertices)
703 p1 = gmsh.model.geo.addPoint(xc_mesh[0], xc_mesh[1], xc_mesh[2], h) # Center
704 p2 = gmsh.model.geo.addPoint(v1[0], v1[1], v1[2], h) # First vertex
705 p3 = gmsh.model.geo.addPoint(v2[0], v2[1], v2[2], h) # Second vertex
706
707 # Create lines
708 l1 = gmsh.model.geo.addLine(p1, p2) # Center to first vertex
709 l2 = gmsh.model.geo.addLine(p2, p3) # First to second vertex
710 l3 = gmsh.model.geo.addLine(p3, p1) # Second vertex back to center
711
712 # Create curve loop and surface
713 c1 = gmsh.model.geo.addCurveLoop([l1, l2, l3])
714 s1 = gmsh.model.geo.addPlaneSurface([c1])
715
716 gmsh.model.geo.synchronize()
717 gmsh.model.mesh.generate(3)
718
719 # Get mesh data for rotation
720 m = get_gmsh_entities()
721
722 # Rotate twice by 120 degrees to complete the triangle
723 for i in range(2):
724 angle = (i + 1) * 2*np.pi/3 # 120 degrees in radians
725 gmsh_transform_general(m, 1000*(i+1), 1000000*(i+1), 1000000*(i+1), angle)
726
727 # Remove duplicate nodes
728 gmsh.model.mesh.removeDuplicateNodes()
729
730 # Translate to specified center
731 gmsh_translate(xc)
732
733 else:
734 # For non-symmetric mesh, create full triangle directly
735 x_tri = get_ref_triangle_points(xc, r)
736
737 # Create points for the full triangle
738 points = []
739 points.append(gmsh.model.geo.addPoint(xc[0], xc[1], xc[2], h)) # Center point
740 for v in x_tri:
741 points.append(gmsh.model.geo.addPoint(v[0], v[1], v[2], h))
742
743 # Create lines connecting vertices
744 lines = []
745 for i in range(2):
746 lines.append(gmsh.model.geo.addLine(points[i+1], points[i+2]))
747 lines.append(gmsh.model.geo.addLine(points[3], points[1]))
748
749 # Create curve loop and surface
750 c1 = gmsh.model.geo.addCurveLoop(lines)
751 s1 = gmsh.model.geo.addPlaneSurface([c1])
752
753 # Synchronize
754 gmsh.model.geo.synchronize()
755
756 # Embed center point in surface
757 gmsh.model.mesh.embed(0, [points[0]], 2, s1)
758
759 # Generate mesh
760 gmsh.model.mesh.generate(3)
761
762 # Check for hanging nodes before writing
763 check_hanging_nodes()
764
765 # Write output files
766 gmsh.write(filename + '.msh')
767 if vtk_out:
768 gmsh.write(filename + '.vtk')
769
770 gmsh.finalize()
771
772def cuboid_mesh_symmetric(xc = [0., 0., 0.], Lx = 1., Ly = 1., Lz = 1., h = 0.1, filename = 'mesh', vtk_out = False, symmetric_mesh = True):
773 """
774 Create a cuboid mesh, either symmetric or full.
775 xc - center point for symmetric mesh, or bottom-left-back corner for non-symmetric mesh
776 Lx - length in x direction
777 Ly - length in y direction
778 Lz - length in z direction
779 h - mesh size
780 symmetric_mesh - if True, creates 1/8 mesh and mirrors it. If False, creates full cuboid
781 """
782 gmsh.initialize()
783 gmsh.option.setNumber("Mesh.MshFileVersion", 2.2)
784
785 if symmetric_mesh:
786 # Create 1/8th of cuboid at origin, then mirror it
787 xc_mesh = [0., 0., 0.]
788
789 # Create box for first octant
790 box = gmsh.model.occ.addBox(xc_mesh[0], xc_mesh[1], xc_mesh[2],
791 0.5*Lx, 0.5*Ly, 0.5*Lz)
792
793 gmsh.model.occ.synchronize()
794
795 # Set mesh size
796 gmsh.model.mesh.setSize(gmsh.model.getEntities(0), h)
797
798 # Generate 3D mesh
799 gmsh.model.mesh.generate(3)
800
801 # Get mesh data for mirroring
802 m = get_gmsh_entities()
803
804 # Mirror the mesh to create full cuboid
805 # First create the other octants in the positive z hemisphere
806 gmsh_transform(m, 1000, 1000000, 1000000, -1, 1, 1) # Mirror across yz plane
807 gmsh_transform(m, 2000, 2000000, 2000000, 1, -1, 1) # Mirror across xz plane
808 gmsh_transform(m, 3000, 3000000, 3000000, -1, -1, 1) # Mirror across z axis
809
810 # Then mirror the entire upper hemisphere to create lower hemisphere
811 gmsh_transform(m, 4000, 4000000, 4000000, 1, 1, -1) # Mirror across xy plane
812 gmsh_transform(m, 5000, 5000000, 5000000, -1, 1, -1) # Mirror across yz plane
813 gmsh_transform(m, 6000, 6000000, 6000000, 1, -1, -1) # Mirror across xz plane
814 gmsh_transform(m, 7000, 7000000, 7000000, -1, -1, -1) # Mirror across z axis
815
816 # Remove duplicate nodes
817 gmsh.model.mesh.removeDuplicateNodes()
818
819 # Translate to specified center coordinates
820 gmsh_translate(xc)
821
822 else:
823 # Create full cuboid directly using OpenCASCADE
824 # Create box centered at xc
825 box = gmsh.model.occ.addBox(xc[0] - 0.5*Lx, xc[1] - 0.5*Ly, xc[2] - 0.5*Lz,
826 Lx, Ly, Lz)
827
828 # Add center point
829 p1 = gmsh.model.occ.addPoint(xc[0], xc[1], xc[2], h)
830
831 # Synchronize
832 gmsh.model.occ.synchronize()
833
834 # Embed center point in volume
835 gmsh.model.mesh.embed(0, [p1], 3, box)
836
837 # Set mesh size
838 gmsh.model.mesh.setSize(gmsh.model.getEntities(0), h)
839
840 # Generate 3D mesh
841 gmsh.model.mesh.generate(3)
842
843 # Check for hanging nodes before writing
844 check_hanging_nodes()
845
846 # Write output files
847 gmsh.write(filename + '.msh')
848 if vtk_out:
849 gmsh.write(filename + '.vtk')
850
851 gmsh.finalize()
852
853def ellipsoid_mesh_symmetric(xc = [0., 0., 0.], rx = 1., ry = 0.5, rz = 0.3, h = 0.1, filename = 'mesh', vtk_out = False, symmetric_mesh = True):
854 """
855 Create an ellipsoid mesh, either symmetric or full.
856 xc - center point coordinates [x, y, z]
857 rx - radius in x direction
858 ry - radius in y direction
859 rz - radius in z direction
860 h - mesh size
861 symmetric_mesh - if True, creates 1/8 mesh and mirrors it. If False, creates full ellipsoid
862 """
863 gmsh.initialize()
864 gmsh.option.setNumber("Mesh.MshFileVersion", 2.2)
865
866 if symmetric_mesh:
867 # Create 1/8th of ellipsoid at origin, then mirror it
868 xc_mesh = [0., 0., 0.]
869
870 # Create sphere of radius 1
871 sphere = gmsh.model.occ.addSphere(xc_mesh[0], xc_mesh[1], xc_mesh[2], 1.0)
872
873 # Scale it to create ellipsoid
874 gmsh.model.occ.dilate([(3, sphere)], xc_mesh[0], xc_mesh[1], xc_mesh[2], rx, ry, rz)
875
876 # Create cutting boxes for first octant
877 cut_left = gmsh.model.occ.addBox(xc_mesh[0]-rx, xc_mesh[1]-ry, xc_mesh[2]-rz, rx, 2*ry, 2*rz)
878 cut_back = gmsh.model.occ.addBox(xc_mesh[0]-rx, xc_mesh[1]-ry, xc_mesh[2]-rz, 2*rx, ry, 2*rz)
879 cut_bottom = gmsh.model.occ.addBox(xc_mesh[0]-rx, xc_mesh[1]-ry, xc_mesh[2]-rz, 2*rx, 2*ry, rz)
880
881 # Cut sphere to get first octant
882 gmsh.model.occ.cut([(3,sphere)], [(3,cut_left), (3,cut_back), (3,cut_bottom)])
883
884 gmsh.model.occ.synchronize()
885
886 # Set mesh size
887 gmsh.model.mesh.setSize(gmsh.model.getEntities(0), h)
888
889 # Generate 3D mesh
890 gmsh.model.mesh.generate(3)
891
892 # Get mesh data for mirroring
893 m = get_gmsh_entities()
894
895 # Mirror the mesh to create full ellipsoid
896 # First create the other octants in the positive z hemisphere
897 gmsh_transform(m, 1000, 1000000, 1000000, -1, 1, 1) # Mirror across yz plane
898 gmsh_transform(m, 2000, 2000000, 2000000, 1, -1, 1) # Mirror across xz plane
899 gmsh_transform(m, 3000, 3000000, 3000000, -1, -1, 1) # Mirror across z axis
900
901 # Then mirror the entire upper hemisphere to create lower hemisphere
902 gmsh_transform(m, 4000, 4000000, 4000000, 1, 1, -1) # Mirror across xy plane
903 gmsh_transform(m, 5000, 5000000, 5000000, -1, 1, -1) # Mirror across yz plane
904 gmsh_transform(m, 6000, 6000000, 6000000, 1, -1, -1) # Mirror across xz plane
905 gmsh_transform(m, 7000, 7000000, 7000000, -1, -1, -1) # Mirror across z axis
906
907 # Remove duplicate nodes
908 gmsh.model.mesh.removeDuplicateNodes()
909
910 # Translate to specified center coordinates
911 gmsh_translate(xc)
912
913 else:
914 # Create full ellipsoid directly using OpenCASCADE
915 # First create sphere of radius 1
916 sphere = gmsh.model.occ.addSphere(xc[0], xc[1], xc[2], 1.0)
917
918 # Scale it to create ellipsoid
919 gmsh.model.occ.dilate([(3, sphere)], xc[0], xc[1], xc[2], rx, ry, rz)
920
921 # Synchronize
922 gmsh.model.occ.synchronize()
923
924 # Set mesh size
925 gmsh.model.mesh.setSize(gmsh.model.getEntities(0), h)
926
927 # Generate 3D mesh
928 gmsh.model.mesh.generate(3)
929
930 # Check for hanging nodes before writing
931 check_hanging_nodes()
932
933 # Write output files
934 gmsh.write(filename + '.msh')
935 if vtk_out:
936 gmsh.write(filename + '.vtk')
937
938 gmsh.finalize()
939
940def cylinder_mesh_symmetric(xc = [0., 0., 0.], r = 1., h = 1., mesh_size = 0.1, filename = 'mesh', vtk_out = False, symmetric_mesh = True):
941 """
942 Create a cylinder mesh, either symmetric or full.
943 xc - center point coordinates [x, y, z] (center of bottom face)
944 r - radius of cylinder
945 h - height of cylinder
946 mesh_size - mesh element size
947 symmetric_mesh - if True, creates 1/8 mesh and mirrors it. If False, creates full cylinder
948 """
949 gmsh.initialize()
950 gmsh.option.setNumber("Mesh.MshFileVersion", 2.2)
951
952 if symmetric_mesh:
953 # Create 1/8th of cylinder at origin, then mirror it
954 xc_mesh = [0., 0., 0.]
955
956 # Create full cylinder first
957 cyl = gmsh.model.occ.addCylinder(xc_mesh[0], xc_mesh[1], xc_mesh[2],
958 0, 0, h, # height vector in z direction
959 r)
960
961 # Create cutting boxes for first octant
962 cut_left = gmsh.model.occ.addBox(xc_mesh[0]-r, xc_mesh[1]-r, xc_mesh[2], r, 2*r, h)
963 cut_back = gmsh.model.occ.addBox(xc_mesh[0]-r, xc_mesh[1]-r, xc_mesh[2], 2*r, r, h)
964
965 # Cut cylinder to get first quarter
966 gmsh.model.occ.cut([(3,cyl)], [(3,cut_left), (3,cut_back)])
967
968 gmsh.model.occ.synchronize()
969
970 # Set mesh size
971 gmsh.model.mesh.setSize(gmsh.model.getEntities(0), mesh_size)
972
973 # Generate 3D mesh
974 gmsh.model.mesh.generate(3)
975
976 # Get mesh data for mirroring
977 m = get_gmsh_entities()
978
979 # Mirror the mesh to create full cylinder
980 # First create the other quarters
981 gmsh_transform(m, 1000, 1000000, 1000000, -1, 1, 1) # Mirror across yz plane
982 gmsh_transform(m, 2000, 2000000, 2000000, 1, -1, 1) # Mirror across xz plane
983 gmsh_transform(m, 3000, 3000000, 3000000, -1, -1, 1) # Mirror across z axis
984
985 # Remove duplicate nodes
986 gmsh.model.mesh.removeDuplicateNodes()
987
988 # Translate to specified center coordinates
989 gmsh_translate(xc)
990
991 else:
992 # Create full cylinder directly using OpenCASCADE
993 cyl = gmsh.model.occ.addCylinder(xc[0], xc[1], xc[2], # base center
994 0, 0, h, # height vector in z direction
995 r) # radius
996
997 # Add center points at top and bottom faces
998 p1 = gmsh.model.occ.addPoint(xc[0], xc[1], xc[2], mesh_size) # bottom center
999 p2 = gmsh.model.occ.addPoint(xc[0], xc[1], xc[2] + h, mesh_size) # top center
1000
1001 # Synchronize
1002 gmsh.model.occ.synchronize()
1003
1004 # Embed center points in respective faces
1005 # First get all surfaces
1006 surfaces = gmsh.model.getEntities(2)
1007 # Find top and bottom surfaces (they are parallel to xy-plane)
1008 for s in surfaces:
1009 com = gmsh.model.occ.getCenterOfMass(s[0], s[1])
1010 if abs(com[2] - xc[2]) < 1e-6: # bottom face
1011 gmsh.model.mesh.embed(0, [p1], 2, s[1])
1012 elif abs(com[2] - (xc[2] + h)) < 1e-6: # top face
1013 gmsh.model.mesh.embed(0, [p2], 2, s[1])
1014
1015 # Set mesh size
1016 gmsh.model.mesh.setSize(gmsh.model.getEntities(0), mesh_size)
1017
1018 # Generate 3D mesh
1019 gmsh.model.mesh.generate(3)
1020
1021 # Check for hanging nodes before writing
1022 check_hanging_nodes()
1023
1024 # Write output files
1025 gmsh.write(filename + '.msh')
1026 if vtk_out:
1027 gmsh.write(filename + '.vtk')
1028
1029 gmsh.finalize()
1030
1031def annulus_circle_mesh_symmetric(xc = [0., 0., 0.], r_outer = 1., r_inner = 0.5, h = 0.1, filename = 'mesh', vtk_out = False, symmetric_mesh = True):
1032 """
1033 Create an annulus (ring) mesh, either symmetric or full.
1034 xc - center point coordinates [x, y, z]
1035 r_outer - radius of outer circle
1036 r_inner - radius of inner circle
1037 h - mesh size
1038 symmetric_mesh - if True, creates 1/4 mesh and mirrors it. If False, creates full annulus
1039 """
1040 if r_inner >= r_outer:
1041 raise ValueError("Inner radius must be smaller than outer radius")
1042
1043 gmsh.initialize()
1044 gmsh.option.setNumber("Mesh.MshFileVersion", 2.2)
1045
1046 if symmetric_mesh:
1047 xc_mesh = [0., 0., 0.]
1048
1049 # add points for the quadrant of annulus circle
1050 p1 = gmsh.model.geo.addPoint(xc_mesh[0], xc_mesh[1], xc_mesh[2], h)
1051 p2 = gmsh.model.geo.addPoint(xc_mesh[0] + r_outer, xc_mesh[1], xc_mesh[2], h)
1052 p3 = gmsh.model.geo.addPoint(xc_mesh[0], xc_mesh[1] + r_outer, xc_mesh[2], h)
1053 p4 = gmsh.model.geo.addPoint(xc_mesh[0] + r_inner, xc_mesh[1], xc_mesh[2], h)
1054 p5 = gmsh.model.geo.addPoint(xc_mesh[0], xc_mesh[1] + r_inner, xc_mesh[2], h)
1055
1056 l1 = gmsh.model.geo.addCircleArc(p2, p1, p3)
1057 l2 = gmsh.model.geo.addCircleArc(p4, p1, p5)
1058 l3 = gmsh.model.geo.addLine(p4, p2)
1059 l4 = gmsh.model.geo.addLine(p5, p3)
1060
1061 c1 = gmsh.model.geo.addCurveLoop([l3, l1, -l4, -l2])
1062 s1 = gmsh.model.geo.addPlaneSurface([c1])
1063
1064 gmsh.model.geo.synchronize()
1065 gmsh.model.mesh.generate(3)
1066
1067 m = get_gmsh_entities()
1068
1069 gmsh_transform(m, 1000, 1000000, 1000000, -1, 1, 1) # x-axis
1070 gmsh_transform(m, 2000, 2000000, 2000000, 1, -1, 1) # y-axis
1071 gmsh_transform(m, 3000, 3000000, 3000000, -1, -1, 1) # z-axis
1072
1073 # remove the duplicate nodes that will have been created on the internal
1074 # boundaries
1075 gmsh.model.mesh.removeDuplicateNodes()
1076
1077 # translate the mesh to the specified center coordinates
1078 gmsh_translate(xc)
1079
1080 else:
1081
1082 # Outer circle points
1083 p1 = gmsh.model.geo.addPoint(xc[0], xc[1], xc[2], h)
1084 p2 = gmsh.model.geo.addPoint(xc[0] + r_outer, xc[1], xc[2], h)
1085 p3 = gmsh.model.geo.addPoint(xc[0], xc[1] + r_outer, xc[2], h)
1086 p4 = gmsh.model.geo.addPoint(xc[0] - r_outer, xc[1], xc[2], h)
1087 p5 = gmsh.model.geo.addPoint(xc[0], xc[1] - r_outer, xc[2], h)
1088
1089 # Inner circle points
1090 p6 = gmsh.model.geo.addPoint(xc[0] + r_inner, xc[1], xc[2], h)
1091 p7 = gmsh.model.geo.addPoint(xc[0], xc[1] + r_inner, xc[2], h)
1092 p8 = gmsh.model.geo.addPoint(xc[0] - r_inner, xc[1], xc[2], h)
1093 p9 = gmsh.model.geo.addPoint(xc[0], xc[1] - r_inner, xc[2], h)
1094
1095 # Create circular arcs for outer circle
1096 c1 = gmsh.model.geo.addCircleArc(p2, p1, p3)
1097 c2 = gmsh.model.geo.addCircleArc(p3, p1, p4)
1098 c3 = gmsh.model.geo.addCircleArc(p4, p1, p5)
1099 c4 = gmsh.model.geo.addCircleArc(p5, p1, p2)
1100
1101 # Create circular arcs for inner circle
1102 c5 = gmsh.model.geo.addCircleArc(p6, p1, p7)
1103 c6 = gmsh.model.geo.addCircleArc(p7, p1, p8)
1104 c7 = gmsh.model.geo.addCircleArc(p8, p1, p9)
1105 c8 = gmsh.model.geo.addCircleArc(p9, p1, p6)
1106
1107 # Create curve loops
1108 outer_loop = gmsh.model.geo.addCurveLoop([c1, c2, c3, c4])
1109 inner_loop = gmsh.model.geo.addCurveLoop([c5, c6, c7, c8])
1110
1111 #gmsh.fltk.run()
1112
1113 # surface
1114 s1 = gmsh.model.geo.addPlaneSurface([outer_loop, inner_loop])
1115
1116 # Synchronize geometry before adding physical groups
1117 gmsh.model.geo.synchronize()
1118
1119 # Add physical group for the surface
1120 gmsh.model.addPhysicalGroup(2, [s1], 1)
1121
1122 # Generate mesh
1123 gmsh.model.mesh.generate(2)
1124
1125 # Check for hanging nodes before writing
1126 check_hanging_nodes()
1127
1128 # Write output files
1129 gmsh.write(filename + '.msh')
1130 if vtk_out:
1131 gmsh.write(filename + '.vtk')
1132
1133 gmsh.finalize()
1134
1135def annulus_rectangle_mesh(xc=[0., 0., 0.], Lx=1., Ly=1., hole_Lx=0.3, hole_Ly=0.3, h=0.1, filename='mesh', vtk_out=False, symmetric_mesh=True):
1136 """
1137 Create a rectangular mesh with a rectangular hole in the center (annulus rectangle).
1138
1139 Parameters
1140 ----------
1141 xc : list
1142 Center coordinates [x, y, z]
1143 Lx : float
1144 Length of outer rectangle in x direction
1145 Ly : float
1146 Length of outer rectangle in y direction
1147 hole_Lx : float
1148 Length of inner hole in x direction
1149 hole_Ly : float
1150 Length of inner hole in y direction
1151 h : float
1152 Mesh size
1153 filename : str
1154 Output filename without extension
1155 vtk_out : bool
1156 If True, also writes a VTK file
1157 symmetric_mesh : bool
1158 If True, creates 1/4 mesh and mirrors it. If False, creates full mesh
1159 """
1160 # Initialize gmsh
1161 gmsh.initialize()
1162 gmsh.option.setNumber("Mesh.MshFileVersion", 2.2)
1163
1164 if symmetric_mesh:
1165 # Create 1/4 mesh at origin first, then mirror and translate
1166 xc_mesh = [0., 0., 0.]
1167
1168 # Create points for outer rectangle (1/4)
1169 p1 = gmsh.model.geo.addPoint(xc_mesh[0] + 0.5*Lx, xc_mesh[1], xc_mesh[2], h)
1170 p2 = gmsh.model.geo.addPoint(xc_mesh[0] + 0.5*Lx, xc_mesh[1] + 0.5*Ly, xc_mesh[2], h)
1171 p3 = gmsh.model.geo.addPoint(xc_mesh[0], xc_mesh[1] + 0.5*Ly, xc_mesh[2], h)
1172 p4 = gmsh.model.geo.addPoint(xc_mesh[0] + 0.5*hole_Lx, xc_mesh[1], xc_mesh[2], h)
1173 p5 = gmsh.model.geo.addPoint(xc_mesh[0] + 0.5*hole_Lx, xc_mesh[1] + 0.5*hole_Ly, xc_mesh[2], h)
1174 p6 = gmsh.model.geo.addPoint(xc_mesh[0], xc_mesh[1] + 0.5*hole_Ly, xc_mesh[2], h)
1175
1176 # Create lines
1177 l1 = gmsh.model.geo.addLine(p1, p2)
1178 l2 = gmsh.model.geo.addLine(p2, p3)
1179 l3 = gmsh.model.geo.addLine(p3, p6)
1180 l4 = gmsh.model.geo.addLine(p6, p5)
1181 l5 = gmsh.model.geo.addLine(p5, p4)
1182 l6 = gmsh.model.geo.addLine(p4, p1)
1183
1184 # Create curve loops
1185 cl = gmsh.model.geo.addCurveLoop([l1, l2, l3, l4, l5, l6])
1186
1187 # Create plane surface with hole
1188 s = gmsh.model.geo.addPlaneSurface([cl])
1189
1190 # Synchronize and generate mesh
1191 gmsh.model.geo.synchronize()
1192 gmsh.model.mesh.generate(2)
1193
1194 # Get mesh data for mirroring
1195 m = get_gmsh_entities()
1196
1197 # Mirror the mesh in all quadrants
1198 gmsh_transform(m, 1000, 1000000, 1000000, -1, 1, 1) # Mirror across y-axis
1199 gmsh_transform(m, 2000, 2000000, 2000000, 1, -1, 1) # Mirror across x-axis
1200 gmsh_transform(m, 3000, 3000000, 3000000, -1, -1, 1) # Mirror across origin
1201
1202 # Remove duplicate nodes
1203 gmsh.model.mesh.removeDuplicateNodes()
1204
1205 # Translate to specified center coordinates
1206 gmsh_translate(xc)
1207
1208 else:
1209 # Create points for outer rectangle
1210 p1 = gmsh.model.geo.addPoint(xc[0] - 0.5*Lx, xc[1] - 0.5*Ly, xc[2], h) # Bottom left
1211 p2 = gmsh.model.geo.addPoint(xc[0] + 0.5*Lx, xc[1] - 0.5*Ly, xc[2], h) # Bottom right
1212 p3 = gmsh.model.geo.addPoint(xc[0] + 0.5*Lx, xc[1] + 0.5*Ly, xc[2], h) # Top right
1213 p4 = gmsh.model.geo.addPoint(xc[0] - 0.5*Lx, xc[1] + 0.5*Ly, xc[2], h) # Top left
1214
1215 # Create points for inner hole
1216 p5 = gmsh.model.geo.addPoint(xc[0] - 0.5*hole_Lx, xc[1] - 0.5*hole_Ly, xc[2], h) # Bottom left
1217 p6 = gmsh.model.geo.addPoint(xc[0] + 0.5*hole_Lx, xc[1] - 0.5*hole_Ly, xc[2], h) # Bottom right
1218 p7 = gmsh.model.geo.addPoint(xc[0] + 0.5*hole_Lx, xc[1] + 0.5*hole_Ly, xc[2], h) # Top right
1219 p8 = gmsh.model.geo.addPoint(xc[0] - 0.5*hole_Lx, xc[1] + 0.5*hole_Ly, xc[2], h) # Top left
1220
1221 # Create lines for outer rectangle
1222 l1 = gmsh.model.geo.addLine(p1, p2) # Bottom
1223 l2 = gmsh.model.geo.addLine(p2, p3) # Right
1224 l3 = gmsh.model.geo.addLine(p3, p4) # Top
1225 l4 = gmsh.model.geo.addLine(p4, p1) # Left
1226
1227 # Create lines for inner hole
1228 l5 = gmsh.model.geo.addLine(p5, p6) # Bottom
1229 l6 = gmsh.model.geo.addLine(p6, p7) # Right
1230 l7 = gmsh.model.geo.addLine(p7, p8) # Top
1231 l8 = gmsh.model.geo.addLine(p8, p5) # Left
1232
1233 # Create curve loops
1234 outer_loop = gmsh.model.geo.addCurveLoop([l1, l2, l3, l4])
1235 inner_loop = gmsh.model.geo.addCurveLoop([l5, l6, l7, l8])
1236
1237 # Create plane surface with hole
1238 s1 = gmsh.model.geo.addPlaneSurface([outer_loop, inner_loop])
1239
1240 # Synchronize and generate mesh
1241 gmsh.model.geo.synchronize()
1242 gmsh.model.mesh.generate(2)
1243
1244 # Check for hanging nodes before writing
1245 check_hanging_nodes()
1246
1247 # Write mesh files
1248 gmsh.write(filename + '.msh')
1249 if vtk_out:
1250 gmsh.write(filename + '.vtk')
1251
1252 # Finalize gmsh
1253 gmsh.finalize()
1254
1255def open_rectangle_mesh(xc = [0., 0., 0.], Lx = 1., Ly = 1., hole_Lx = 0.5, hole_Ly = 0.5, h = 0.1, filename = 'mesh', vtk_out = False):
1256 """
1257 Create a rectangular mesh with a rectangular hole in the center (annulus rectangle).
1258
1259 Parameters
1260 ----------
1261 xc : list
1262 Center coordinates [x, y, z]
1263 Lx : float
1264 Length of outer rectangle in x direction
1265 Ly : float
1266 Length of outer rectangle in y direction
1267 hole_Lx : float
1268 Length of inner hole in x direction
1269 hole_Ly : float
1270 Length of inner hole in y direction
1271 h : float
1272 Mesh size
1273 filename : str
1274 Output filename without extension
1275 vtk_out : bool
1276 If True, also writes a VTK file
1277 """
1278 # Initialize gmsh
1279 gmsh.initialize()
1280 gmsh.option.setNumber("Mesh.MshFileVersion", 2.2)
1281
1282 # Annulus rectangle with top removed
1283 # xc is the center of the cocentric rectangles
1284 # we mesh the right half of the annulus rectangle
1285 xc_mesh = [0., 0., 0.]
1286 points = []
1287 points.append(gmsh.model.geo.addPoint(xc_mesh[0] + 0.5*Lx, xc_mesh[1] - 0.5*Ly, xc_mesh[2], h))
1288 points.append(gmsh.model.geo.addPoint(xc_mesh[0] + 0.5*Lx, xc_mesh[1] + 0.5*Ly, xc_mesh[2], h))
1289 points.append(gmsh.model.geo.addPoint(xc_mesh[0] + 0.5*hole_Lx, xc_mesh[1] + 0.5*Ly, xc_mesh[2], h))
1290 points.append(gmsh.model.geo.addPoint(xc_mesh[0] + 0.5*hole_Lx, xc_mesh[1] - 0.5*hole_Ly, xc_mesh[2], h))
1291 points.append(gmsh.model.geo.addPoint(xc_mesh[0], xc_mesh[1] - 0.5*hole_Ly, xc_mesh[2], h))
1292 points.append(gmsh.model.geo.addPoint(xc_mesh[0], xc_mesh[1] - 0.5*Ly, xc_mesh[2], h))
1293
1294 lines = []
1295 for i in range(len(points) - 1):
1296 lines.append(gmsh.model.geo.addLine(points[i], points[i + 1]))
1297 lines.append(gmsh.model.geo.addLine(points[-1], points[0]))
1298 cl = gmsh.model.geo.addCurveLoop(lines)
1299 s = gmsh.model.geo.addPlaneSurface([cl])
1300 gmsh.model.geo.synchronize()
1301 gmsh.model.mesh.generate(2)
1302
1303 m = get_gmsh_entities()
1304
1305 # mirror the mesh across the x-axis
1306 gmsh_transform(m, 1000, 1000000, 1000000, -1, 1, 1)
1307 gmsh.model.mesh.removeDuplicateNodes()
1308 gmsh_translate(xc)
1309
1310 # Check for hanging nodes before writing
1311 check_hanging_nodes()
1312
1313 gmsh.write(filename + '.msh')
1314 if vtk_out:
1315 gmsh.write(filename + '.vtk')
1316
1317 gmsh.finalize()
1318
1319def open_pipe_mesh(xc=[0., 0., 0.], axis=[0., 0., 1.], length=2., outer_radius=1., wall_thickness=0.1, h=0.1, filename='mesh', vtk_out=False):
1320 """
1321 Create a 3D pipe mesh with specified axis, closed bottom and open top.
1322 The pipe is created by:
1323 1. Creating an outer cylinder
1324 2. Subtracting an inner cylinder to create walls
1325 3. Subtracting the top surface to create the opening
1326
1327 Parameters
1328 ----------
1329 xc : list
1330 Center coordinates [x, y, z] of the base center
1331 axis : list
1332 Axis vector defining pipe orientation [ax, ay, az]
1333 length : float
1334 Length of the pipe along axis
1335 outer_radius : float
1336 Outer radius of the pipe
1337 wall_thickness : float
1338 Thickness of the pipe wall and bottom
1339 h : float
1340 Mesh size
1341 filename : str
1342 Output filename without extension
1343 vtk_out : bool
1344 If True, also writes a VTK file
1345 """
1346 gmsh.initialize()
1347 gmsh.option.setNumber("Mesh.MshFileVersion", 2.2)
1348
1349 # Normalize axis vector
1350 axis_norm = np.sqrt(axis[0]**2 + axis[1]**2 + axis[2]**2)
1351 if axis_norm == 0:
1352 raise ValueError("Axis vector cannot be zero")
1353 axis = [x/axis_norm for x in axis]
1354
1355 # Create outer cylinder
1356 l = length + wall_thickness
1357 outer_cylinder = gmsh.model.occ.addCylinder(xc[0], xc[1], xc[2],
1358 axis[0]*l, axis[1]*l, axis[2]*l,
1359 outer_radius)
1360
1361 # Create inner cylinder (to be subtracted)
1362 inner_cylinder = gmsh.model.occ.addCylinder(xc[0]+wall_thickness*axis[0], xc[1]+wall_thickness*axis[1], xc[2]+wall_thickness*axis[2],
1363 axis[0]*l, axis[1]*l, axis[2]*l,
1364 outer_radius - wall_thickness)
1365
1366 # Create a box slightly larger than the cylinder at the top to cut off the top surface
1367 # First create a box aligned with coordinate axes
1368 dx = 2 * outer_radius
1369 box = gmsh.model.occ.addBox(xc[0] - dx, xc[1] - dx, xc[2] + l - wall_thickness,
1370 2*dx, 2*dx, 2*wall_thickness)
1371
1372 # If axis is not aligned with z, rotate the box to align with axis
1373 if not (abs(axis[0]) < 1e-10 and abs(axis[1]) < 1e-10):
1374 # Calculate rotation angle and axis
1375 z_axis = [0, 0, 1]
1376 rotation_axis = np.cross(z_axis, axis)
1377 rotation_angle = np.arccos(np.dot(z_axis, axis))
1378
1379 # Rotate the box around center point
1380 gmsh.model.occ.rotate([(3, box)],
1381 xc[0], xc[1], xc[2],
1382 rotation_axis[0], rotation_axis[1], rotation_axis[2],
1383 rotation_angle)
1384
1385 # First subtract inner cylinder from outer cylinder
1386 pipe = gmsh.model.occ.cut([(3, outer_cylinder)], [(3, inner_cylinder), (3, box)])
1387
1388 # Synchronize before meshing
1389 gmsh.model.occ.synchronize()
1390
1391 # Set mesh size
1392 gmsh.model.mesh.setSize(gmsh.model.getEntities(0), h)
1393
1394 # Generate 3D mesh
1395 gmsh.model.mesh.generate(3)
1396
1397 # Check for hanging nodes before writing
1398 check_hanging_nodes()
1399
1400 # Write mesh files
1401 gmsh.write(filename + '.msh')
1402 if vtk_out:
1403 gmsh.write(filename + '.vtk')
1404
1405 gmsh.finalize()
1406
1407
1409if __name__ == "__main__":
1410 inp_dir = './'
1411
1412 test_meshes = ['circle', 'ellipse', 'sphere', 'cuboid', 'ellipsoid', 'rectangle', \
1413 'hexagon', 'drum2d', 'triangle', 'polygon', 'cylindrical2d_wall', \
1414 'cylinder', 'annulus_circle', 'annulus_rectangle', 'open_rectangle', 'open_pipe']
1415 test_meshes = ['open_pipe']
1416 symm_flag = 0 # 0 - both, 1 - symmetric, 2 - non-symmetric
1417
1418 for mesh in test_meshes:
1419 if mesh == 'circle':
1420 if symm_flag == 0 or symm_flag == 1:
1421 circle_mesh_symmetric(xc = [-3., -3., 0.], r = 1., h = 0.1, filename = './' + mesh + '_sym', vtk_out = True, symmetric_mesh = True)
1422 if symm_flag == 0 or symm_flag == 2:
1423 circle_mesh_symmetric(xc = [-3., -3., 0.], r = 1., h = 0.1, filename = './' + mesh + '_non_sym', vtk_out = True, symmetric_mesh = False)
1424 elif mesh == 'ellipse':
1425 if symm_flag == 0 or symm_flag == 1:
1426 ellipse_mesh_symmetric(xc = [-1., -1., 0.], rx = 1., ry = 0.5, h = 0.1, filename = './' + mesh + '_sym', vtk_out = True, symmetric_mesh = True)
1427 if symm_flag == 0 or symm_flag == 2:
1428 ellipse_mesh_symmetric(xc = [-1., -1., 0.], rx = 1., ry = 0.5, h = 0.1, filename = './' + mesh + '_non_sym', vtk_out = True, symmetric_mesh = False)
1429 elif mesh == 'sphere':
1430 if symm_flag == 0 or symm_flag == 1:
1431 sphere_mesh_symmetric(xc = [-5., -5., 0.], r = 1., h = 0.1, filename = './' + mesh + '_sym', vtk_out = True, symmetric_mesh = True)
1432 if symm_flag == 0 or symm_flag == 2:
1433 sphere_mesh_symmetric(xc = [-5., -5., 0.], r = 1., h = 0.1, filename = './' + mesh + '_non_sym', vtk_out = True, symmetric_mesh = False)
1434 elif mesh == 'rectangle':
1435 if symm_flag == 0 or symm_flag == 1:
1436 rectangle_mesh_symmetric(xc = [1., 1., 0.], Lx = 1., Ly = 1., h = 0.1, filename = './' + mesh + '_sym', vtk_out = True, symmetric_mesh = True)
1437 if symm_flag == 0 or symm_flag == 2:
1438 rectangle_mesh_symmetric(xc = [1., 1., 0.], Lx = 1., Ly = 1., h = 0.1, filename = './' + mesh + '_non_sym', vtk_out = True, symmetric_mesh = False)
1439 elif mesh == 'hexagon':
1440 if symm_flag == 0 or symm_flag == 1:
1441 hexagon_mesh_symmetric(xc = [3., 3., 0.], r = 1., h = 0.1, filename = './' + mesh + '_sym', vtk_out = True, symmetric_mesh = True)
1442 if symm_flag == 0 or symm_flag == 2:
1443 hexagon_mesh_symmetric(xc = [3., 3., 0.], r = 1., h = 0.1, filename = './' + mesh + '_non_sym', vtk_out = True, symmetric_mesh = False)
1444 elif mesh == 'drum2d':
1445 if symm_flag == 0 or symm_flag == 1:
1446 drum2d_mesh_symmetric(xc = [5., 5., 0.], r = 1., width = 0.5, h = 0.1, filename = './' + mesh + '_sym', vtk_out = True, symmetric_mesh = True)
1447 if symm_flag == 0 or symm_flag == 2:
1448 drum2d_mesh_symmetric(xc = [5., 5., 0.], r = 1., width = 0.5, h = 0.1, filename = './' + mesh + '_non_sym', vtk_out = True, symmetric_mesh = False)
1449 elif mesh == 'triangle':
1450 if symm_flag == 0 or symm_flag == 1:
1451 triangle_mesh_symmetric(xc = [7., 7., 0.], r = 1., h = 0.1, filename = './' + mesh + '_sym', vtk_out = True, symmetric_mesh = True)
1452 if symm_flag == 0 or symm_flag == 2:
1453 triangle_mesh_symmetric(xc = [7., 7., 0.], r = 1., h = 0.1, filename = './' + mesh + '_non_sym', vtk_out = True, symmetric_mesh = False)
1454 elif mesh == 'polygon':
1455 theta = np.pi/6.
1456 R, a = 1., 0.25
1457 v1 = [0., 0., 0.]
1458 v2 = [R*np.cos(0.5*theta), -R*np.sin(0.5*theta), 0.]
1459 v4 = [R*np.cos(0.5*theta), R*np.sin(0.5*theta), 0.]
1460 v3 = [R + a, 0., 0.]
1461 if symm_flag == 0 or symm_flag == 1:
1462 polygon_mesh_symmetric(points = [v1, v2, v3, v4], theta = theta, xc = [9., 9., 9.], h = 0.1, filename = './' + mesh + '_sym', vtk_out = True, symmetric_mesh = True)
1463 if symm_flag == 0 or symm_flag == 2:
1464 polygon_mesh_symmetric(points = [v1, v2, v3, v4], theta = theta, xc = [9., 9., 9.], h = 0.1, filename = './' + mesh + '_non_sym', vtk_out = True, symmetric_mesh = False)
1465 elif mesh == 'cylindrical2d_wall':
1466 cylindrical2d_wall_mesh(
1467 center=[0., 0., 0.],
1468 outer_radius=1.0,
1469 inner_radius=0.8,
1470 bar_width=0.2,
1471 bar_length=0.3,
1472 h=0.05,
1473 filename='./' + mesh + '_sym',
1474 vtk_out=True
1475 )
1476 elif mesh == 'cuboid':
1477 if symm_flag == 0 or symm_flag == 1:
1478 cuboid_mesh_symmetric(xc = [0., 0., 0.], Lx = 1., Ly = 0.5, Lz = 0.3, h = 0.1, filename = './' + mesh + '_sym', vtk_out = True, symmetric_mesh = True)
1479 if symm_flag == 0 or symm_flag == 2:
1480 cuboid_mesh_symmetric(xc = [0., 0., 0.], Lx = 1., Ly = 0.5, Lz = 0.3, h = 0.1, filename = './' + mesh + '_non_sym', vtk_out = True, symmetric_mesh = False)
1481 elif mesh == 'ellipsoid':
1482 if symm_flag == 0 or symm_flag == 1:
1483 ellipsoid_mesh_symmetric(xc = [0., 0., 0.], rx = 1., ry = 0.5, rz = 0.3, h = 0.1, filename = './' + mesh + '_sym', vtk_out = True, symmetric_mesh = True)
1484 if symm_flag == 0 or symm_flag == 2:
1485 ellipsoid_mesh_symmetric(xc = [0., 0., 0.], rx = 1., ry = 0.5, rz = 0.3, h = 0.1, filename = './' + mesh + '_non_sym', vtk_out = True, symmetric_mesh = False)
1486 elif mesh == 'cylinder':
1487 if symm_flag == 0 or symm_flag == 1:
1488 cylinder_mesh_symmetric(xc = [0., 0., 0.], r = 1., h = 2., mesh_size = 0.1, filename = './' + mesh + '_sym', vtk_out = True, symmetric_mesh = True)
1489 if symm_flag == 0 or symm_flag == 2:
1490 cylinder_mesh_symmetric(xc = [0., 0., 0.], r = 1., h = 2., mesh_size = 0.1, filename = './' + mesh + '_non_sym', vtk_out = True, symmetric_mesh = False)
1491 elif mesh == 'annulus_circle':
1492 if symm_flag == 0 or symm_flag == 1:
1493 annulus_circle_mesh_symmetric(xc = [0., 0., 0.], r_outer = 1., r_inner = 0.5, h = 0.1, filename = './' + mesh + '_sym', vtk_out = True, symmetric_mesh = True)
1494 if symm_flag == 0 or symm_flag == 2:
1495 annulus_circle_mesh_symmetric(xc = [0., 0., 0.], r_outer = 1., r_inner = 0.5, h = 0.1, filename = './' + mesh + '_non_sym', vtk_out = True, symmetric_mesh = False)
1496 elif mesh == 'annulus_rectangle':
1497 if symm_flag == 0 or symm_flag == 1:
1498 annulus_rectangle_mesh(xc = [0., 0., 0.], Lx = 1., Ly = 0.8, hole_Lx = 0.2, hole_Ly = 0.4, h = 0.1, filename = './' + mesh + '_sym', vtk_out = True, symmetric_mesh = True)
1499 if symm_flag == 0 or symm_flag == 2:
1500 annulus_rectangle_mesh(xc = [0., 0., 0.], Lx = 1., Ly = 0.8, hole_Lx = 0.2, hole_Ly = 0.4, h = 0.1, filename = './' + mesh + '_non_sym', vtk_out = True, symmetric_mesh = False)
1501 elif mesh == 'open_rectangle':
1502 open_rectangle_mesh(xc = [0., 0., 0.], Lx = 1.2, Ly = 0.8, hole_Lx = 1., hole_Ly = 0.6, h = 0.1, filename = './' + mesh, vtk_out = True)
1503 elif mesh == 'open_pipe':
1504 open_pipe_mesh(xc=[0., 0., 0.], axis=[1., 1., 0.], length=2., outer_radius=0.5,
1505 wall_thickness=0.2, h=0.1, filename='./' + mesh, vtk_out=True)
1506 else:
1507 print(f"Mesh {mesh} not found")
1508
geom_util.gmsh_translate
gmsh_translate(xc)
Definition geom_util.py:361
geom_util.get_ref_drum_points
get_ref_drum_points(center, radius, width, add_center=False)
Definition geom_util.py:161
geom_util.check_hanging_nodes
check_hanging_nodes()
Definition geom_util.py:383
geom_util.get_gmsh_entities
get_gmsh_entities()
Definition geom_util.py:295
geom_util.gmsh_transform
gmsh_transform(m, offset_entity, offset_node, offset_element, tx, ty, tz)
Definition geom_util.py:307
geom_util.get_ref_triangle_points
get_ref_triangle_points(center, radius, add_center=False)
Definition geom_util.py:68
geom_util.gmsh_transform_general
gmsh_transform_general(m, offset_entity, offset_node, offset_element, angle, axis=[0, 0, 1])
Definition geom_util.py:328
gmsh_particles.open_pipe_mesh
open_pipe_mesh(xc=[0., 0., 0.], axis=[0., 0., 1.], length=2., outer_radius=1., wall_thickness=0.1, h=0.1, filename='mesh', vtk_out=False)
Definition gmsh_particles.py:1319
gmsh_particles.triangle_mesh_symmetric
triangle_mesh_symmetric(xc=[0., 0., 0.], r=1., h=0.1, filename='mesh', vtk_out=False, symmetric_mesh=True)
Definition gmsh_particles.py:680
gmsh_particles.sphere_mesh_symmetric
sphere_mesh_symmetric(xc=[0., 0., 0.], r=1., h=0.1, filename='mesh', vtk_out=False, symmetric_mesh=True)
Definition gmsh_particles.py:175
gmsh_particles.open_rectangle_mesh
open_rectangle_mesh(xc=[0., 0., 0.], Lx=1., Ly=1., hole_Lx=0.5, hole_Ly=0.5, h=0.1, filename='mesh', vtk_out=False)
Definition gmsh_particles.py:1255
gmsh_particles.cylinder_mesh_symmetric
cylinder_mesh_symmetric(xc=[0., 0., 0.], r=1., h=1., mesh_size=0.1, filename='mesh', vtk_out=False, symmetric_mesh=True)
Definition gmsh_particles.py:940
gmsh_particles.annulus_circle_mesh_symmetric
annulus_circle_mesh_symmetric(xc=[0., 0., 0.], r_outer=1., r_inner=0.5, h=0.1, filename='mesh', vtk_out=False, symmetric_mesh=True)
Definition gmsh_particles.py:1031
gmsh_particles.hexagon_mesh_symmetric
hexagon_mesh_symmetric(xc=[0., 0., 0.], r=1., h=0.1, filename='mesh', vtk_out=False, symmetric_mesh=True)
Definition gmsh_particles.py:351
gmsh_particles.circle_mesh_symmetric
circle_mesh_symmetric(xc=[0., 0., 0.], r=1., h=0.1, filename='mesh', vtk_out=False, symmetric_mesh=True)
Definition gmsh_particles.py:12
gmsh_particles.polygon_mesh_symmetric
polygon_mesh_symmetric(points, theta, xc=[0., 0., 0.], h=0.1, filename='mesh', vtk_out=False, symmetric_mesh=True)
Definition gmsh_particles.py:526
gmsh_particles.cylindrical2d_wall_mesh
cylindrical2d_wall_mesh(center=[0., 0., 0.], outer_radius=1.0, inner_radius=0.8, bar_width=0.2, bar_length=0.3, h=0.1, filename='mesh', vtk_out=False)
Definition gmsh_particles.py:599
gmsh_particles.drum2d_mesh_symmetric
drum2d_mesh_symmetric(xc=[0., 0., 0.], r=1., width=1., h=0.1, filename='mesh', vtk_out=False, symmetric_mesh=True)
Definition gmsh_particles.py:435
gmsh_particles.annulus_rectangle_mesh
annulus_rectangle_mesh(xc=[0., 0., 0.], Lx=1., Ly=1., hole_Lx=0.3, hole_Ly=0.3, h=0.1, filename='mesh', vtk_out=False, symmetric_mesh=True)
Definition gmsh_particles.py:1135
gmsh_particles.cuboid_mesh_symmetric
cuboid_mesh_symmetric(xc=[0., 0., 0.], Lx=1., Ly=1., Lz=1., h=0.1, filename='mesh', vtk_out=False, symmetric_mesh=True)
Definition gmsh_particles.py:772
gmsh_particles.rectangle_mesh_symmetric
rectangle_mesh_symmetric(xc=[0., 0., 0.], Lx=1., Ly=1., h=0.1, filename='mesh', vtk_out=False, symmetric_mesh=True)
Definition gmsh_particles.py:263
gmsh_particles.ellipsoid_mesh_symmetric
ellipsoid_mesh_symmetric(xc=[0., 0., 0.], rx=1., ry=0.5, rz=0.3, h=0.1, filename='mesh', vtk_out=False, symmetric_mesh=True)
Definition gmsh_particles.py:853
gmsh_particles.ellipse_mesh_symmetric
ellipse_mesh_symmetric(xc=[0., 0., 0.], rx=1., ry=0.5, h=0.1, filename='mesh', vtk_out=False, symmetric_mesh=True)
Definition gmsh_particles.py:92