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