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