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   
  5   
  6   
  7  # geometry building ----------------------------------------------------------- 
8 -def bounding_box(points):
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
16 -def is_multi_polygon(facets):
17 if not len(facets): 18 return False 19 20 try: 21 facets[0][0][0] # facet 0, poly 0, point 0 22 except TypeError: 23 # pure python raises this 24 return False 25 except IndexError: 26 # numpy raises this 27 return False 28 else: 29 return True
30 31 32 33
34 -def offset_point_indices(facets, offset):
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
45 -class GeometryBuilder(object):
46 - def __init__(self):
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
102 - def dimensions(self):
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
113 - def mesher_module(self):
114 dim = self.dimensions() 115 if dim == 2: 116 import meshpy.triangle 117 return meshpy.triangle 118 elif dim == 3: 119 import meshpy.tet 120 return meshpy.tet 121 else: 122 raise ValueError, "unsupported dimensionality %d" % dim
123
124 - def bounding_box(self):
125 return bounding_box(self.points)
126
127 - def center(self):
128 a, b = bounding_box(self.points) 129 return (a+b)/2
130
131 - def wrap_in_box(self, distance, subdivisions=None):
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
143 - def apply_transform(self, f):
144 self.points = [f(x) for x in self.points]
145 146 147 148 149 # actual geometries -----------------------------------------------------------
150 -class Marker:
151 MINUS_X = 1 152 PLUS_X = 2 153 MINUS_Y = 3 154 PLUS_Y = 4 155 MINUS_Z = 5 156 PLUS_Z = 6 157 SHELL = 100 158 159 FIRST_USER_MARKER = 1000
160 161 162 163
164 -def make_box(a, b, subdivisions=None):
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 # CAUTION: Do not change point or facet order here. 178 # Other code depends on this staying the way it is. 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 # 7--------6 195 # /| /| 196 # 4--------5 | z 197 # | | | | ^ 198 # | 3------|-2 | y 199 # |/ |/ |/ 200 # 0--------1 +--->x 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
244 -def make_circle(r, center=(0,0), subdivisions=40, marker=Marker.SHELL):
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
265 -def make_ball(r, subdivisions=10):
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 # extrusions ------------------------------------------------------------------
304 -def _is_same_float(a, b, threshold=1e-10):
305 if abs(a) > abs(b): 306 a,b = b,a 307 308 # now abs(a) <= abs(b) always 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 # need to generate fresh ring, continue 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 # we're moving purely outward--this doesn't need fans, only plane 422 # surfaces. Special casing this leads to more freedom for TetGen 423 # and hence better meshes. 424 425 if r1 == 0: 426 # make opening surface 427 if r2 != 0: 428 add_polygons([get_ring(ring2_idx)], marker=marker) 429 elif r2 == 0: 430 # make closing surface 431 add_polygons([get_ring(ring1_idx)], marker=marker) 432 else: 433 # make single-surface interface with hole 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 # make opening fan 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 # make closing fan 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 # make quad strip 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 # pre-populate ring dict with ring_point_indices 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
506 -def generate_surface_of_revolution(rz_points, 507 closure=EXT_OPEN, radial_subdiv=16, 508 point_idx_offset=0, ring_point_indices=None, 509 ring_markers=None, rz_closure_marker=0):
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