1 from __future__ import division
2 import numpy as np
3
4 __doc__ = """
5
6 Geometry builder
7 ----------------
8
9 .. autoclass:: GeometryBuilder
10
11 Geometries
12 ----------
13
14 These functions are designed so that their output can be splat-passed to
15 :meth:`GeometryBuilder.add_geometry`::
16
17 builder = GeometryBuilder()
18 builder.add_geometry(*make_ball(10))
19
20 .. autoclass:: Marker
21 :members:
22 :undoc-members:
23
24 .. autofunction:: make_box
25 .. autofunction:: make_circle
26 .. autofunction:: make_ball
27 .. autofunction:: make_cylinder
28
29 Extrusions and surfaces of revolution
30 -------------------------------------
31
32 .. data:: EXT_OPEN
33 .. data:: EXT_CLOSED_IN_RZ
34
35 .. autofunction:: generate_extrusion
36 .. autofunction:: generate_surface_of_revolution
37 """
38
39
40
42 return (
43 np.asarray(np.min(points, axis=0), dtype=np.float64),
44 np.asarray(np.max(points, axis=0), dtype=np.float64))
45
46
48 if not len(facets):
49 return False
50
51 try:
52 facets[0][0][0]
53 except TypeError:
54
55 return False
56 except IndexError:
57
58 return False
59 else:
60 return True
61
62
64 if is_multi_polygon(facets):
65 return [[tuple(p_i+offset for p_i in poly)
66 for poly in facet]
67 for facet in facets]
68 else:
69 return [tuple(p_i+offset for p_i in facet) for facet in facets]
70
71
73 """
74 .. automethod:: add_geometry
75 .. automethod:: set
76 .. automethod:: wrap_in_box
77 .. automethod:: bounding_box
78 .. automethod:: center
79 .. automethod:: apply_transform
80 """
82 self.points = []
83 self.facets = []
84 self.facet_hole_starts = None
85 self.facet_markers = None
86 self.point_markers = None
87
88 - def add_geometry(self, points, facets, facet_hole_starts=None,
89 facet_markers=None, point_markers=None):
90 if isinstance(facet_markers, int):
91 facet_markers = len(facets) * [facet_markers]
92
93 if facet_hole_starts and not self.facet_hole_starts:
94 self.facet_hole_starts = len(self.facets) * []
95 if facet_markers and not self.facet_markers:
96 self.facet_markers = len(self.facets) * [0]
97 if point_markers and not self.point_markers:
98 self.point_markers = len(self.points) * [0]
99
100 if not facet_hole_starts and self.facet_hole_starts:
101 facet_hole_starts = len(facets) * [[]]
102 if not facet_markers and self.facet_markers:
103 facet_markers = len(facets) * [0]
104 if not point_markers and self.point_markers:
105 point_markers = len(points) * [0]
106
107 if is_multi_polygon(facets) and not is_multi_polygon(self.facets):
108 self.facets = [[facet] for facet in self.facets]
109
110 if not is_multi_polygon(facets) and is_multi_polygon(self.facets):
111 facets = [[facet] for facet in facets]
112
113 self.facets.extend(offset_point_indices(facets, len(self.points)))
114 self.points.extend(points)
115
116 if facet_markers:
117 self.facet_markers.extend(facet_markers)
118 assert len(facets) == len(facet_markers)
119 if facet_hole_starts:
120 self.facet_hole_starts.extend(facet_hole_starts)
121 assert len(facets) == len(facet_hole_starts)
122 if point_markers:
123 self.point_markers.extend(point_markers)
124 assert len(points) == len(point_markers)
125
126 - def add_cycle(self, points, facet_markers=None, point_markers=None):
127 def make_facets():
128 end = len(points)-1
129 for i in range(end):
130 yield i, i+1
131 yield end, 0
132
133 self.add_geometry(points, list(make_facets()),
134 facet_markers=facet_markers,
135 point_markers=point_markers)
136
138 return len(self.points[0])
139
140 - def set(self, mesh_info):
141 """Transfer the built geometry into a :class:`meshpy.triangle.MeshInfo`
142 or a :class:`meshpy.tet.MeshInfo`.
143 """
144
145 mesh_info.set_points(self.points, self.point_markers)
146 if self.facet_hole_starts or is_multi_polygon(self.facets):
147 mesh_info.set_facets_ex(self.facets,
148 self.facet_hole_starts, self.facet_markers)
149 else:
150 mesh_info.set_facets(self.facets, self.facet_markers)
151
162
165
169
171 """
172 :param subdivisions: is a tuple of integers specifying
173 the number of subdivisions along each axis.
174 """
175
176 a, b = bounding_box(self.points)
177 points, facets, _, facet_markers = \
178 make_box(a-distance, b+distance, subdivisions)
179
180 self.add_geometry(points, facets, facet_markers=facet_markers)
181
184
185
186
187
188
189
200
201
203 """
204 :param subdivisions: is a tuple of integers specifying
205 the number of subdivisions along each axis.
206 """
207
208 a = [float(ai) for ai in a]
209 b = [float(bi) for bi in b]
210
211 assert len(a) == len(b)
212
213 dimensions = len(a)
214 if dimensions == 2:
215
216
217
218 points = [
219 (a[0], a[1]),
220 (b[0], a[1]),
221 (b[0], b[1]),
222 (a[0], b[1]),
223 ]
224
225 facets = [(0, 1), (1, 2), (2, 3), (3, 0)]
226
227 facet_markers = [
228 Marker.MINUS_Y, Marker.PLUS_X,
229 Marker.PLUS_Y, Marker.MINUS_X]
230
231 elif dimensions == 3:
232
233
234
235
236
237
238
239
240 points = [
241 (a[0], a[1], a[2]),
242 (b[0], a[1], a[2]),
243 (b[0], b[1], a[2]),
244 (a[0], b[1], a[2]),
245 (a[0], a[1], b[2]),
246 (b[0], a[1], b[2]),
247 (b[0], b[1], b[2]),
248 (a[0], b[1], b[2]),
249 ]
250
251 facets = [
252 (0, 1, 2, 3),
253 (0, 1, 5, 4),
254 (1, 2, 6, 5),
255 (7, 6, 2, 3),
256 (7, 3, 0, 4),
257 (4, 5, 6, 7)
258 ]
259
260 facet_markers = [Marker.MINUS_Z, Marker.MINUS_Y, Marker.PLUS_X,
261 Marker.PLUS_Y, Marker.MINUS_X, Marker.PLUS_Z]
262 else:
263 raise ValueError("unsupported dimension count: %d" % len(a))
264
265 if subdivisions is not None:
266 if dimensions != 2:
267 raise NotImplementedError(
268 "subdivision not implemented for any "
269 "dimension count other than 2")
270
271 from meshpy.triangle import subdivide_facets
272 points, facets, facet_markers = subdivide_facets(
273 [subdivisions[0], subdivisions[1],
274 subdivisions[0], subdivisions[1]],
275 points, facets, facet_markers)
276
277 return points, facets, None, facet_markers
278
279
281 def round_trip_connect(seq):
282 result = []
283 for i in range(len(seq)):
284 result.append((i, (i+1) % len(seq)))
285 return result
286
287 phi = np.linspace(0, 2*np.pi, num=subdivisions, endpoint=False)
288 cx, cy = center
289 x = r*np.cos(phi) + cx
290 y = r*np.sin(phi) + cy
291
292 return ([np.array(pt) for pt in zip(x, y)],
293 round_trip_connect(range(subdivisions)),
294 None,
295 subdivisions*[marker])
296
297
299 from math import pi, cos, sin
300
301 dphi = pi/subdivisions
302
303 def truncate(my_r):
304 if abs(my_r) < 1e-9*r:
305 return 0
306 else:
307 return my_r
308
309 rz = [(truncate(r*sin(i*dphi)), r*cos(i*dphi)) for i in range(subdivisions+1)]
310
311 return generate_surface_of_revolution(
312 rz, closure=EXT_OPEN, radial_subdiv=subdivisions)
313
314
315 -def make_cylinder(radius, height, radial_subdivisions=10,
316 height_subdivisions=1):
317 dz = height/height_subdivisions
318 rz = [(0, 0)] \
319 + [(radius, i*dz) for i in range(height_subdivisions+1)] \
320 + [(0, height)]
321 ring_markers = [Marker.MINUS_Z] \
322 + ((height_subdivisions)*[Marker.SHELL]) \
323 + [Marker.PLUS_Z]
324
325 return generate_surface_of_revolution(rz,
326 closure=EXT_OPEN, radial_subdiv=radial_subdivisions,
327 ring_markers=ring_markers)
328
329
330
331
332
333
335 if abs(a) > abs(b):
336 a, b = b, a
337
338
339 return abs(b) < threshold or abs(a-b) < threshold*abs(b)
340
341
342 EXT_OPEN = 0
343 EXT_CLOSED_IN_RZ = 1
344
345
346 -def generate_extrusion(rz_points, base_shape, closure=EXT_OPEN,
347 point_idx_offset=0, ring_point_indices=None,
348 ring_markers=None, rz_closure_marker=0):
349 """Extrude a given connected *base_shape* (a list of (x,y) points)
350 along the z axis. For each step in the extrusion, the base shape
351 is multiplied by a radius and shifted in the z direction. Radius
352 and z offset are given by *rz_points*, which is a list of
353 (r, z) tuples.
354
355 Returns ``(points, facets, facet_holestarts, markers)``, where *points* is
356 a list of (3D) points and facets is a list of polygons. Each polygon is, in
357 turn, represented by a tuple of indices into *points*. If
358 *point_idx_offset* is not zero, these indices start at that number.
359 *markers* is a list equal in length to *facets*, each specifying the facet
360 marker of that facet. *facet_holestarts* is also equal in length to
361 *facets*, each element is a list of hole starting points for the
362 corresponding facet.
363
364 Use L{MeshInfo.set_facets_ex} to add the extrusion to a L{MeshInfo}
365 structure.
366
367 The extrusion proceeds by generating quadrilaterals connecting each
368 ring. If any given radius in *rz_points* is 0, triangle fans are
369 produced instead of quads to provide non-degenerate closure.
370
371 If *closure* is :data:`EXT_OPEN`, no efforts are made to put end caps on the
372 extrusion.
373
374 If *closure* is :data:`EXT_CLOSED_IN_RZ`, then a torus-like structure
375 is assumed and the last ring is just connected to the first.
376
377 If *ring_markers* is not None, it is an list of markers added to each
378 ring. There should be len(rz_points)-1 entries in this list.
379 If rings are added because of closure options, they receive the
380 corresponding *XXX_closure_marker*. If *facet_markers* is given, this function
381 returns (points, facets, markers), where markers is is a list containing
382 a marker for each generated facet. Unspecified markers generally
383 default to 0.
384
385 If *ring_point_indices* is given, it must be a list of the same
386 length as *rz_points*. Each entry in the list may either be None,
387 or a list of point indices. This list must contain the same number
388 of points as the *base_shape*; it is taken as the indices of
389 pre-existing points that are to be used for the given ring, instead
390 of generating new points.
391 """
392
393 assert len(rz_points) > 0
394
395 if ring_markers is not None:
396 assert len(rz_points) == len(ring_markers)+1
397
398 def get_ring(ring_idx):
399 try:
400 return rings[ring_idx]
401 except KeyError:
402
403 pass
404
405 p_indices = None
406 if ring_point_indices is not None:
407 p_indices = ring_point_indices[ring_idx]
408
409 first_idx = point_idx_offset+len(points)
410
411 r, z = rz_points[ring_idx]
412
413 if r == 0:
414 p_indices = (first_idx,)
415 points.append((0, 0, z))
416 else:
417 p_indices = tuple(xrange(first_idx, first_idx+len(base_shape)))
418 points.extend([(x*r, y*r, z) for (x, y) in base_shape])
419
420 rings[ring_idx] = p_indices
421 return p_indices
422
423 def pair_with_successor(l):
424 n = len(l)
425 return [(l[i], l[(i+1) % n]) for i in range(n)]
426
427 def add_polygons(new_polys, marker):
428 """Add several new facets, each polygon in new_polys corresponding
429 to a new facet.
430 """
431 facets.extend([poly] for poly in new_polys)
432 markers.extend(len(new_polys)*[marker])
433 holelists.extend(len(new_polys)*[[]])
434
435 def add_facet(facet_polygons, holestarts, marker):
436 """Add a single facet, with each polygon in *facet_polygons*
437 belonging to a single facet.
438 """
439 facets.append(facet_polygons)
440 markers.append(marker)
441 holelists.append(holestarts)
442
443 def connect_ring(ring1_idx, ring2_idx, marker):
444 r1, z1 = rz_points[ring1_idx]
445 r2, z2 = rz_points[ring2_idx]
446
447 if _is_same_float(z2, z1):
448 assert not _is_same_float(r1, r2)
449
450
451
452
453 if r1 == 0:
454
455 if r2 != 0:
456 add_polygons([get_ring(ring2_idx)], marker=marker)
457 elif r2 == 0:
458
459 add_polygons([get_ring(ring1_idx)], marker=marker)
460 else:
461
462 add_facet([
463 get_ring(ring1_idx),
464 get_ring(ring2_idx),
465 ],
466 holestarts=[(0, 0, z1)], marker=marker)
467 else:
468 ring1 = get_ring(ring1_idx)
469 ring2 = get_ring(ring2_idx)
470 if r1 == 0:
471
472 assert len(ring1) == 1
473 start_pt = ring1[0]
474
475 if r2 != 0:
476 add_polygons(
477 [(start_pt, succ, pt)
478 for pt, succ in pair_with_successor(ring2)],
479 marker=marker)
480 elif r2 == 0:
481
482 assert len(ring2) == 1
483 end_pt = ring2[0]
484 add_polygons(
485 [(pt, succ, end_pt)
486 for pt, succ in pair_with_successor(ring1)],
487 marker=marker)
488 else:
489
490 pairs1 = pair_with_successor(ring1)
491 pairs2 = pair_with_successor(ring2)
492 add_polygons(
493 [(a, b, c, d) for ((a, b), (d, c)) in zip(pairs1, pairs2)],
494 marker=marker)
495
496 points = []
497 facets = []
498 markers = []
499 holelists = []
500
501 rings = {}
502
503
504 if ring_point_indices is not None:
505 for i, ring_points in enumerate(ring_point_indices):
506 if ring_points is not None:
507 assert isinstance(ring_points, tuple)
508
509 if rz_points[i][0] == 0:
510 assert len(ring_points) == 1
511 else:
512 assert len(ring_points) == len(base_shape), \
513 ("Ring points length (%d) does not "
514 "match base shape length (%d)"
515 % (len(ring_points), len(base_shape)))
516
517 rings[i] = ring_points
518
519 for i in range(len(rz_points)-1):
520 if ring_markers is not None:
521 ring_marker = ring_markers[i]
522 else:
523 ring_marker = 0
524
525 connect_ring(i, i+1, ring_marker)
526
527 if closure == EXT_CLOSED_IN_RZ:
528 connect_ring(len(rz_points)-1, 0, rz_closure_marker)
529
530 return points, facets, holelists, markers
531
532
537 from math import sin, cos, pi
538
539 dphi = 2*pi/radial_subdiv
540 base_shape = [(cos(dphi*i), sin(dphi*i)) for i in range(radial_subdiv)]
541 return generate_extrusion(
542 rz_points, base_shape, closure=closure,
543 point_idx_offset=point_idx_offset,
544 ring_point_indices=ring_point_indices,
545 ring_markers=ring_markers, rz_closure_marker=rz_closure_marker,
546 )
547
548
549
550
551