Package meshpy :: Module geometry
[hide private]
[frames] | no frames]

Source Code for Module meshpy.geometry

  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  # {{{ geometry building 
 40   
41 -def bounding_box(points):
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
47 -def is_multi_polygon(facets):
48 if not len(facets): 49 return False 50 51 try: 52 facets[0][0][0] # facet 0, poly 0, point 0 53 except TypeError: 54 # pure python raises this 55 return False 56 except IndexError: 57 # numpy raises this 58 return False 59 else: 60 return True
61 62
63 -def offset_point_indices(facets, offset):
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
72 -class GeometryBuilder(object):
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 """
81 - def __init__(self):
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
137 - def dimensions(self):
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
152 - def mesher_module(self):
153 dim = self.dimensions() 154 if dim == 2: 155 import meshpy.triangle 156 return meshpy.triangle 157 elif dim == 3: 158 import meshpy.tet 159 return meshpy.tet 160 else: 161 raise ValueError("unsupported dimensionality %d" % dim)
162
163 - def bounding_box(self):
164 return bounding_box(self.points)
165
166 - def center(self):
167 a, b = bounding_box(self.points) 168 return (a+b)/2
169
170 - def wrap_in_box(self, distance, subdivisions=None):
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
182 - def apply_transform(self, f):
183 self.points = [f(x) for x in self.points]
184 185 # }}} 186 187 188 # {{{ actual geometries 189
190 -class Marker:
191 MINUS_X = 1 192 PLUS_X = 2 193 MINUS_Y = 3 194 PLUS_Y = 4 195 MINUS_Z = 5 196 PLUS_Z = 6 197 SHELL = 100 198 199 FIRST_USER_MARKER = 1000
200 201
202 -def make_box(a, b, subdivisions=None):
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 # CAUTION: Do not change point or facet order here. 216 # Other code depends on this staying the way it is. 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 # 7--------6 233 # /| /| 234 # 4--------5 | z 235 # | | | | ^ 236 # | 3------|-2 | y 237 # |/ |/ |/ 238 # 0--------1 +--->x 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
280 -def make_circle(r, center=(0, 0), subdivisions=40, marker=Marker.SHELL):
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
298 -def make_ball(r, subdivisions=10):
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 # {{{ extrusions 333
334 -def _is_same_float(a, b, threshold=1e-10):
335 if abs(a) > abs(b): 336 a, b = b, a 337 338 # now abs(a) <= abs(b) always 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 # need to generate fresh ring, continue 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 # we're moving purely outward--this doesn't need fans, only plane 450 # surfaces. Special casing this leads to more freedom for TetGen 451 # and hence better meshes. 452 453 if r1 == 0: 454 # make opening surface 455 if r2 != 0: 456 add_polygons([get_ring(ring2_idx)], marker=marker) 457 elif r2 == 0: 458 # make closing surface 459 add_polygons([get_ring(ring1_idx)], marker=marker) 460 else: 461 # make single-surface interface with hole 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 # make opening fan 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 # make closing fan 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 # make quad strip 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 # pre-populate ring dict with ring_point_indices 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
533 -def generate_surface_of_revolution(rz_points, 534 closure=EXT_OPEN, radial_subdiv=16, 535 point_idx_offset=0, ring_point_indices=None, 536 ring_markers=None, rz_closure_marker=0):
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 # vim: foldmethod=marker 551