PeriDEM 0.2.0
PeriDEM -- Peridynamics-based high-fidelity model for granular media
Loading...
Searching...
No Matches
Functions
geom_util Namespace Reference

Functions

 rotate (p, theta, axis)
 
 get_ref_rect_points (center, Lx, Ly, add_center=False)
 
 get_ref_triangle_points (center, radius, add_center=False)
 
 get_ref_hex_points (center, radius, add_center=False)
 
 get_ref_drum_points (center, radius, width, add_center=False)
 
 does_intersect (p, particles, padding)
 
 does_intersect_rect (p, particles, padding, rect, is_3d=False)
 
 get_gmsh_entities ()
 
 gmsh_transform (m, offset_entity, offset_node, offset_element, tx, ty, tz)
 
 gmsh_transform_general (m, offset_entity, offset_node, offset_element, angle, axis=[0, 0, 1])
 
 gmsh_translate (xc)
 
 check_hanging_nodes ()
 

Function Documentation

◆ check_hanging_nodes()

geom_util.check_hanging_nodes ( ) 
Check for hanging/orphan nodes in the GMSH mesh that are not connected to any elements.

Returns:
-------
list
    List of node tags that are hanging (not connected to any elements)
bool
    True if hanging nodes were found, False otherwise

Definition at line 383 of file geom_util.py.

383def check_hanging_nodes():
384 """Check for hanging/orphan nodes in the GMSH mesh that are not connected to any elements.
385
386 Returns:
387 -------
388 list
389 List of node tags that are hanging (not connected to any elements)
390 bool
391 True if hanging nodes were found, False otherwise
392 """
393 # Get all nodes in the mesh
394 all_node_tags, _, _ = gmsh.model.mesh.getNodes()
395 all_node_tags = set(all_node_tags)
396
397 # Get all nodes that are part of elements
398 connected_nodes = set()
399
400 # Iterate through all entities in the model
401 for entity in gmsh.model.getEntities():
402 # Get elements for this entity
403 elem_types, elem_tags, elem_node_tags = gmsh.model.mesh.getElements(entity[0], entity[1])
404
405 # Add all nodes from elements to connected set
406 for nodes in elem_node_tags:
407 connected_nodes.update(nodes)
408
409 # Find hanging nodes (nodes that exist but aren't in any elements)
410 hanging_nodes = list(all_node_tags - connected_nodes)
411
412 has_hanging = len(hanging_nodes) > 0
413 if has_hanging:
414 print(f"Warning: Found {len(hanging_nodes)} hanging nodes with tags: {hanging_nodes}")
415
416 return hanging_nodes, has_hanging

Referenced by gmsh_particles.annulus_circle_mesh_symmetric(), gmsh_particles.annulus_rectangle_mesh(), gmsh_particles.circle_mesh_symmetric(), gmsh_particles.cuboid_mesh_symmetric(), gmsh_particles.cylinder_mesh_symmetric(), gmsh_particles.cylindrical2d_wall_mesh(), gmsh_particles.drum2d_mesh_symmetric(), gmsh_particles.ellipse_mesh_symmetric(), gmsh_particles.ellipsoid_mesh_symmetric(), gmsh_particles.hexagon_mesh_symmetric(), gmsh_particles.open_pipe_mesh(), gmsh_particles.open_rectangle_mesh(), gmsh_particles.polygon_mesh_symmetric(), gmsh_particles.rectangle_mesh_symmetric(), gmsh_particles.sphere_mesh_symmetric(), and gmsh_particles.triangle_mesh_symmetric().

Here is the caller graph for this function:

◆ does_intersect()

geom_util.does_intersect (   p,
  particles,
  padding 
) 
Returns true if particle p intersects or is near enough to existing particles

Parameters
----------
p : list
    Coordinates of center and radius of particle [x,y,z,r]
particles : list
    List of center + radius of multiple particles. E.g. particles[0] is a list containing coordinates of center and radius.
padding: float
    Minimum distance between circle boundaries such that if two circles

Returns
-------
bool
    True if particle intersects or is near enough to one of the particle in the list

Definition at line 225 of file geom_util.py.

225def does_intersect(p, particles, padding):
226 """
227 Returns true if particle p intersects or is near enough to existing particles
228
229 Parameters
230 ----------
231 p : list
232 Coordinates of center and radius of particle [x,y,z,r]
233 particles : list
234 List of center + radius of multiple particles. E.g. particles[0] is a list containing coordinates of center and radius.
235 padding: float
236 Minimum distance between circle boundaries such that if two circles
237
238 Returns
239 -------
240 bool
241 True if particle intersects or is near enough to one of the particle in the list
242 """
243 if len(p) < 4: raise Exception('p = {} must have atleast 4 elements'.format(p))
244 if len(particles) == 0: raise Exception('particles = {} can not be empty'.format(particles))
245 if padding < 0.: raise Exception('padding = {} can not be negative'.format(padding))
246
247 for q in particles:
248 if len(q) < 4: raise Exception('q = {} in particles must have atleast 4 elements'.format(q))
249 pq = np.array([p[i] - q[i] for i in range(3)])
250 if np.linalg.norm(pq) <= p[3] + q[3] + padding:
251 return True
252 return False
253

◆ does_intersect_rect()

geom_util.does_intersect_rect (   p,
  particles,
  padding,
  rect,
  is_3d = False 
) 
Returns true if particle p is sufficiently close or outside the rectangle (in 2d) or cuboid (in 3d)

Parameters
----------
p : list
    Coordinates of center and radius of particle [x,y,z,r]
particles : list
    List of center + radius of multiple particles. E.g. particles[0] is a list containing coordinates of center and radius.
padding: float
    Minimum distance between circle boundaries such that if two circles
rect: list
    Coordinates of left-bottom and right-top corner points of rectangle (2d) or cuboid (3d). E.g. [x1 y1, z1, x2, y2, z2]
is_3d: bool
    True if we are dealing with cuboid

Returns
-------
bool
    True if particle intersects or is near enough to the rectangle

Definition at line 254 of file geom_util.py.

254def does_intersect_rect(p, particles, padding, rect, is_3d = False):
255 """
256 Returns true if particle p is sufficiently close or outside the rectangle (in 2d) or cuboid (in 3d)
257
258 Parameters
259 ----------
260 p : list
261 Coordinates of center and radius of particle [x,y,z,r]
262 particles : list
263 List of center + radius of multiple particles. E.g. particles[0] is a list containing coordinates of center and radius.
264 padding: float
265 Minimum distance between circle boundaries such that if two circles
266 rect: list
267 Coordinates of left-bottom and right-top corner points of rectangle (2d) or cuboid (3d). E.g. [x1 y1, z1, x2, y2, z2]
268 is_3d: bool
269 True if we are dealing with cuboid
270
271 Returns
272 -------
273 bool
274 True if particle intersects or is near enough to the rectangle
275 """
276 if len(p) < 4: raise Exception('p = {} must have atleast 4 elements'.format(p))
277 if len(particles) == 0: raise Exception('particles = {} can not be empty'.format(particles))
278 if padding < 0.: raise Exception('padding = {} can not be negative'.format(padding))
279 if len(rect) < 6: raise Exception('rect = {} must have 6 elements'.format(rect))
280
281 pr = [p[0] - p[3], p[1] - p[3], p[2], p[0] + p[3], p[1] + p[3], p[2]]
282 if is_3d:
283 pr[2] -= p[3]
284 pr[5] += p[3]
285
286 if pr[0] < rect[0] + padding or pr[1] < rect[1] + padding or pr[3] > rect[3] - padding or pr[4] > rect[4] - padding:
287 if is_3d:
288 if pr[2] < rect[2] + padding or pr[5] > rect[5] - padding:
289 return True
290 else:
291 return True
292 return False
293
294

◆ get_gmsh_entities()

geom_util.get_gmsh_entities ( )

Definition at line 295 of file geom_util.py.

295def get_gmsh_entities():
296 m = {}
297 for e in gmsh.model.getEntities():
298 bnd = gmsh.model.getBoundary([e])
299 nod = gmsh.model.mesh.getNodes(e[0], e[1])
300 ele = gmsh.model.mesh.getElements(e[0], e[1])
301 m[e] = (bnd, nod, ele)
302
303 return m
304
305# transform the mesh and create new discrete entities to store it
306# source - https://gitlab.onelab.info/gmsh/gmsh/-/blob/master/examples/api/mirror_mesh.py

Referenced by gmsh_particles.annulus_circle_mesh_symmetric(), gmsh_particles.annulus_rectangle_mesh(), gmsh_particles.circle_mesh_symmetric(), gmsh_particles.cuboid_mesh_symmetric(), gmsh_particles.cylinder_mesh_symmetric(), gmsh_particles.drum2d_mesh_symmetric(), gmsh_particles.ellipse_mesh_symmetric(), gmsh_particles.ellipsoid_mesh_symmetric(), gmsh_particles.hexagon_mesh_symmetric(), gmsh_particles.open_rectangle_mesh(), gmsh_particles.polygon_mesh_symmetric(), gmsh_particles.rectangle_mesh_symmetric(), gmsh_particles.sphere_mesh_symmetric(), and gmsh_particles.triangle_mesh_symmetric().

Here is the caller graph for this function:

◆ get_ref_drum_points()

geom_util.get_ref_drum_points (   center,
  radius,
  width,
  add_center = False 
) 
Returns size points on reference concave polygon (we refer to it as drum)

Reference concave polygon:

             v3                                v2
               + -------------------------------- +
                \                               /
                   \                         /
                     +         o           +
                    /v4        x          v1 \
                 /                              \ 
               + -------------------------------- +
               v5                               v6

Radius: x-v2
Width: v1-v4
Axis vector: x-v1
Rotation axis: [0, 0, 1]

Parameters
----------
center: list
    Coordinates of center of polygon
radius : float
    Radius of polygon (distance between center x and vertex v2)
width : float
    Width of neck, i.e. distance between vertex v2 and v4
add_center : bool
    True if we include the center to the returned list (first point in the returned list will be the center if the flag is true)

Returns
-------
list
    Coordinates of points

Definition at line 161 of file geom_util.py.

161def get_ref_drum_points(center, radius, width, add_center = False):
162 """
163 Returns size points on reference concave polygon (we refer to it as drum)
164
165 Reference concave polygon:
166
167 v3 v2
168 + -------------------------------- +
169 \ /
170 \ /
171 + o +
172 /v4 x v1 \
173 / \
174 + -------------------------------- +
175 v5 v6
176
177 Radius: x-v2
178 Width: v1-v4
179 Axis vector: x-v1
180 Rotation axis: [0, 0, 1]
181
182 Parameters
183 ----------
184 center: list
185 Coordinates of center of polygon
186 radius : float
187 Radius of polygon (distance between center x and vertex v2)
188 width : float
189 Width of neck, i.e. distance between vertex v2 and v4
190 add_center : bool
191 True if we include the center to the returned list (first point in the returned list will be the center if the flag is true)
192
193 Returns
194 -------
195 list
196 Coordinates of points
197 """
198 axis = [1., 0., 0.]
199 rotate_axis = [0., 0., 1.]
200
201 x1 = rotate(axis, np.pi/3., rotate_axis)
202 x2 = rotate(axis, -np.pi/3., rotate_axis)
203
204 points = []
205 if add_center:
206 points.append(center)
207
208 # v1
209 points.append([center[i] + width*0.5*axis[i] for i in range(3)])
210 # v2
211 points.append([center[i] + radius*x1[i] for i in range(3)])
212 # v3
213 points.append([center[i] + radius*x1[i] - radius*axis[i] for i in range(3)])
214 # v4
215 points.append([center[i] - width*0.5*axis[i] for i in range(3)])
216 # v5
217 v6 = [center[i] + radius*x2[i] for i in range(3)]
218 points.append([v6[i] - radius*axis[i] for i in range(3)])
219 # v6
220 points.append(v6)
221
222 return points
223
224

References rotate().

Referenced by gmsh_particles.drum2d_mesh_symmetric().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ get_ref_hex_points()

geom_util.get_ref_hex_points (   center,
  radius,
  add_center = False 
) 
Returns size points on reference hexagon

Reference hexagon:

                        v3           v2
                         +           +


                      +         o           +
                     v4         x            v1

                         +           +
                        v5           v6

Radius: x-v1
Axis vector: x-v1
Rotation axis: [0, 0, 1]

Parameters
----------
center: list
    Coordinates of center of hexagon
radius : float
    Radius of hexagon
add_center : bool
    True if we include the center to the returned list (first point in the returned list will be the center if the flag is true)

Returns
-------
list
    Coordinates of points

Definition at line 115 of file geom_util.py.

115def get_ref_hex_points(center, radius, add_center = False):
116 """
117 Returns size points on reference hexagon
118
119 Reference hexagon:
120
121 v3 v2
122 + +
123
124
125 + o +
126 v4 x v1
127
128 + +
129 v5 v6
130
131 Radius: x-v1
132 Axis vector: x-v1
133 Rotation axis: [0, 0, 1]
134
135 Parameters
136 ----------
137 center: list
138 Coordinates of center of hexagon
139 radius : float
140 Radius of hexagon
141 add_center : bool
142 True if we include the center to the returned list (first point in the returned list will be the center if the flag is true)
143
144 Returns
145 -------
146 list
147 Coordinates of points
148 """
149 axis = [1., 0., 0.]
150 rotate_axis = [0., 0., 1.]
151 points = []
152 if add_center:
153 points.append(center)
154
155 for i in range(6):
156 xi = rotate(axis, i*np.pi/3., rotate_axis)
157 points.append([center[i] + radius * xi[i] for i in range(3)])
158
159 return points
160

References rotate().

Here is the call graph for this function:

◆ get_ref_rect_points()

geom_util.get_ref_rect_points (   center,
  Lx,
  Ly,
  add_center = False 
) 
Returns size points on reference rectangle

Reference rectangle:

+-------------------+
|                   |
|                   |
|                   |
+-------------------+

Parameters
----------
center: list
    Coordinates of center of rectangle
Lx : float
    Length of rectangle in x-direction
Ly : float
    Length of rectangle in y-direction
add_center : bool
    True if we include the center to the returned list (first point in the returned list will be the center if the flag is true)

Returns
-------
list
    Coordinates of points

Definition at line 32 of file geom_util.py.

32def get_ref_rect_points(center, Lx, Ly, add_center = False):
33 """
34 Returns size points on reference rectangle
35
36 Reference rectangle:
37
38 +-------------------+
39 | |
40 | |
41 | |
42 +-------------------+
43
44 Parameters
45 ----------
46 center: list
47 Coordinates of center of rectangle
48 Lx : float
49 Length of rectangle in x-direction
50 Ly : float
51 Length of rectangle in y-direction
52 add_center : bool
53 True if we include the center to the returned list (first point in the returned list will be the center if the flag is true)
54
55 Returns
56 -------
57 list
58 Coordinates of points
59 """
60
61 p1 = [center[0] - Lx/2., center[1] - Ly/2., center[2]]
62 p2 = [center[0] + Lx/2., center[1] - Ly/2., center[2]]
63 p3 = [center[0] + Lx/2., center[1] + Ly/2., center[2]]
64 p4 = [center[0] - Lx/2., center[1] + Ly/2., center[2]]
65
66 return [p1, p2, p3, p4] if not add_center else [center, p1, p2, p3, p4]
67

◆ get_ref_triangle_points()

geom_util.get_ref_triangle_points (   center,
  radius,
  add_center = False 
) 
Returns size points on reference triangle

Reference triangle:
            
#    + v2
#    | \
#    |      \
#    |            \  
#    |    o           + v1
#    |            /
#    |      /
#    | /
#    + v3

Radius: o-v1 
Axis vector: o-v1
Rotation axis: [0, 0, 1]

Parameters
----------
center: list
    Coordinates of center of triangle
radius : float
    Radius of triangle
add_center : bool
    True if we include the center to the returned list (first point in the returned list will be the center if the flag is true)

Returns
-------
list
    Coordinates of points

Definition at line 68 of file geom_util.py.

68def get_ref_triangle_points(center, radius, add_center = False):
69 """
70 Returns size points on reference triangle
71
72 Reference triangle:
73
74 # + v2
75 # | \
76 # | \
77 # | \
78 # | o + v1
79 # | /
80 # | /
81 # | /
82 # + v3
83
84 Radius: o-v1
85 Axis vector: o-v1
86 Rotation axis: [0, 0, 1]
87
88 Parameters
89 ----------
90 center: list
91 Coordinates of center of triangle
92 radius : float
93 Radius of triangle
94 add_center : bool
95 True if we include the center to the returned list (first point in the returned list will be the center if the flag is true)
96
97 Returns
98 -------
99 list
100 Coordinates of points
101 """
102 axis = [1., 0., 0.]
103 rotate_axis = [0., 0., 1.]
104
105 points = []
106 if add_center:
107 points.append(center)
108
109 for i in range(3):
110 xi = rotate(axis, i*2*np.pi/3., rotate_axis)
111 points.append([center[i] + radius * xi[i] for i in range(3)])
112
113 return points
114

References rotate().

Referenced by gmsh_particles.triangle_mesh_symmetric().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ gmsh_transform()

geom_util.gmsh_transform (   m,
  offset_entity,
  offset_node,
  offset_element,
  tx,
  ty,
  tz 
)

Definition at line 307 of file geom_util.py.

307def gmsh_transform(m, offset_entity, offset_node, offset_element, tx, ty, tz):
308 for e in sorted(m):
309 gmsh.model.addDiscreteEntity(
310 e[0], e[1] + offset_entity,
311 [(abs(b[1]) + offset_entity) * int(math.copysign(1, b[1])) for b in m[e][0]])
312 coord = []
313 for i in range(0, len(m[e][1][1]), 3):
314 x = m[e][1][1][i] * tx
315 y = m[e][1][1][i + 1] * ty
316 z = m[e][1][1][i + 2] * tz
317 coord.append(x)
318 coord.append(y)
319 coord.append(z)
320 gmsh.model.mesh.addNodes(e[0], e[1] + offset_entity,
321 [n + offset_node for n in m[e][1][0]], coord)
322 gmsh.model.mesh.addElements(e[0], e[1] + offset_entity, m[e][2][0],
323 [[t + offset_element for t in typ] for typ in m[e][2][1]],
324 [[n + offset_node for n in typ] for typ in m[e][2][2]])
325 if (tx * ty * tz) < 0: # reverse the orientation
326 gmsh.model.mesh.reverse([(e[0], e[1] + offset_entity)])
327

Referenced by gmsh_particles.annulus_circle_mesh_symmetric(), gmsh_particles.annulus_rectangle_mesh(), gmsh_particles.circle_mesh_symmetric(), gmsh_particles.cuboid_mesh_symmetric(), gmsh_particles.cylinder_mesh_symmetric(), gmsh_particles.drum2d_mesh_symmetric(), gmsh_particles.ellipse_mesh_symmetric(), gmsh_particles.ellipsoid_mesh_symmetric(), gmsh_particles.open_rectangle_mesh(), gmsh_particles.rectangle_mesh_symmetric(), and gmsh_particles.sphere_mesh_symmetric().

Here is the caller graph for this function:

◆ gmsh_transform_general()

geom_util.gmsh_transform_general (   m,
  offset_entity,
  offset_node,
  offset_element,
  angle,
  axis = [0, 0, 1] 
) 
Rotate mesh by angle (in radians) around specified axis

Parameters:
m - mesh data dictionary
offset_entity - offset for entity numbers
offset_node - offset for node numbers 
offset_element - offset for element numbers
angle - rotation angle in radians
axis - rotation axis vector (default is z-axis [0,0,1])

Definition at line 328 of file geom_util.py.

328def gmsh_transform_general(m, offset_entity, offset_node, offset_element, angle, axis=[0, 0, 1]):
329 """Rotate mesh by angle (in radians) around specified axis
330
331 Parameters:
332 m - mesh data dictionary
333 offset_entity - offset for entity numbers
334 offset_node - offset for node numbers
335 offset_element - offset for element numbers
336 angle - rotation angle in radians
337 axis - rotation axis vector (default is z-axis [0,0,1])
338 """
339
340 for e in sorted(m):
341 gmsh.model.addDiscreteEntity(
342 e[0], e[1] + offset_entity,
343 [(abs(b[1]) + offset_entity) * int(math.copysign(1, b[1])) for b in m[e][0]])
344 coord = []
345 # Apply rotation to each point
346 for i in range(0, len(m[e][1][1]), 3):
347 xi = np.array([m[e][1][1][i], m[e][1][1][i + 1], m[e][1][1][i + 2]])
348
349 xi_rot = rotate(xi, angle, axis)
350
351 coord.append(xi_rot[0])
352 coord.append(xi_rot[1])
353 coord.append(xi_rot[2])
354
355 gmsh.model.mesh.addNodes(e[0], e[1] + offset_entity,
356 [n + offset_node for n in m[e][1][0]], coord)
357 gmsh.model.mesh.addElements(e[0], e[1] + offset_entity, m[e][2][0],
358 [[t + offset_element for t in typ] for typ in m[e][2][1]],
359 [[n + offset_node for n in typ] for typ in m[e][2][2]])
360

References rotate().

Referenced by gmsh_particles.hexagon_mesh_symmetric(), gmsh_particles.polygon_mesh_symmetric(), and gmsh_particles.triangle_mesh_symmetric().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ gmsh_translate()

geom_util.gmsh_translate (   xc) 
Translate mesh by vector xc

Parameters:
xc - translation vector [x, y, z]

Definition at line 361 of file geom_util.py.

361def gmsh_translate(xc):
362 """Translate mesh by vector xc
363
364 Parameters:
365 xc - translation vector [x, y, z]
366 """
367 # Get the nodes
368 nodeTags, coord, paramCoord = gmsh.model.mesh.getNodes()
369 numNodes = len(nodeTags)
370
371 # Translation vector
372 translation_vector = np.array([xc[0], xc[1], xc[2]])
373
374 # Update each node individually
375 for i in range(numNodes):
376 nodeTag = nodeTags[i]
377 old_coord = coord[3*i:3*(i+1)]
378 new_coord = old_coord + translation_vector
379 # Handle paramCoord properly - use empty list if paramCoord is None or empty
380 param = paramCoord[3*i:3*(i+1)] if paramCoord is not None and len(paramCoord) > 0 else []
381 gmsh.model.mesh.setNode(nodeTag, new_coord, param)
382

Referenced by gmsh_particles.annulus_circle_mesh_symmetric(), gmsh_particles.annulus_rectangle_mesh(), gmsh_particles.circle_mesh_symmetric(), gmsh_particles.cuboid_mesh_symmetric(), gmsh_particles.cylinder_mesh_symmetric(), gmsh_particles.drum2d_mesh_symmetric(), gmsh_particles.ellipse_mesh_symmetric(), gmsh_particles.ellipsoid_mesh_symmetric(), gmsh_particles.hexagon_mesh_symmetric(), gmsh_particles.open_rectangle_mesh(), gmsh_particles.polygon_mesh_symmetric(), gmsh_particles.rectangle_mesh_symmetric(), gmsh_particles.sphere_mesh_symmetric(), and gmsh_particles.triangle_mesh_symmetric().

Here is the caller graph for this function:

◆ rotate()

geom_util.rotate (   p,
  theta,
  axis 
) 
Returns rotation of vector about specified axis by specified angle

Parameters
----------
p: list
    Coordinates of vector
theta : float
    Angle of rotation
axis : list
    Axis of rotation

Returns
-------
list
    Coordinates of rotated vector

Definition at line 7 of file geom_util.py.

7def rotate(p, theta, axis):
8 """
9 Returns rotation of vector about specified axis by specified angle
10
11 Parameters
12 ----------
13 p: list
14 Coordinates of vector
15 theta : float
16 Angle of rotation
17 axis : list
18 Axis of rotation
19
20 Returns
21 -------
22 list
23 Coordinates of rotated vector
24 """
25 p_np, axis_np = np.array(p), np.array(axis)
26 ct, st = np.cos(theta), np.sin(theta)
27 p_dot_n = np.dot(p_np,axis_np)
28 n_cross_p = np.cross(axis_np, p_np)
29
30 return (1. - ct) * p_dot_n * axis_np + ct * p_np + st * n_cross_p
31

Referenced by get_ref_drum_points(), get_ref_hex_points(), get_ref_triangle_points(), and gmsh_transform_general().

Here is the caller graph for this function: