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

Source Code for Module meshpy.geometry

  1  from __future__ import division 
  2  import numpy 
  3   
  4   
  5   
  6   
  7  # geometry building ----------------------------------------------------------- 
8 -def bounding_box(points):
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
16 -def is_multi_polygon(facets):
17 if not facets: 18 return False 19 20 try: 21 facets[0][0][0] # facet 0, poly 0, point 0 22 except TypeError: 23 return False 24 else: 25 return True
26 27 28 29
30 -def offset_point_indices(facets, offset):
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
41 -class GeometryBuilder(object):
42 - def __init__(self):
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
98 - def dimensions(self):
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
109 - def mesher_module(self):
110 dim = self.dimensions() 111 if dim == 2: 112 import meshpy.triangle 113 return meshpy.triangle 114 elif dim == 3: 115 import meshpy.tet 116 return meshpy.tet 117 else: 118 raise ValueError, "unsupported dimensionality %d" % dim
119
120 - def bounding_box(self):
121 return bounding_box(self.points)
122
123 - def center(self):
124 a, b = bounding_box(self.points) 125 return (a+b)/2
126
127 - def wrap_in_box(self, distance, subdivisions=None):
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
139 - def apply_transform(self, f):
140 self.points = [f(x) for x in self.points]
141 142 143 144 145 # actual geometries -----------------------------------------------------------
146 -class Marker:
147 MINUS_X = 1 148 PLUS_X = 2 149 MINUS_Y = 3 150 PLUS_Y = 4 151 MINUS_Z = 5 152 PLUS_Z = 6 153 SHELL = 100 154 155 FIRST_USER_MARKER = 1000
156 157 158 159
160 -def make_box(a, b, subdivisions=None):
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 # CAUTION: Do not change point or facet order here. 171 # Other code depends on this staying the way it is. 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 # 7--------6 188 # /| /| 189 # 4--------5 | z 190 # | | | | ^ 191 # | 3------|-2 | y 192 # |/ |/ |/ 193 # 0--------1 +--->x 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
237 -def make_circle(r, center=(0,0), subdivisions=40, marker=Marker.SHELL):
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
258 -def make_ball(r, subdivisions=10):
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 # extrusions ------------------------------------------------------------------
297 -def _is_same_float(a, b, threshold=1e-10):
298 if abs(a) > abs(b): 299 a,b = b,a 300 301 # now abs(a) <= abs(b) always 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 # need to generate fresh ring, continue 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 # we're moving purely outward--this doesn't need fans, only plane 415 # surfaces. Special casing this leads to more freedom for TetGen 416 # and hence better meshes. 417 418 if r1 == 0: 419 # make opening surface 420 if r2 != 0: 421 add_polygons([get_ring(ring2_idx)], marker=marker) 422 elif r2 == 0: 423 # make closing surface 424 add_polygons([get_ring(ring1_idx)], marker=marker) 425 else: 426 # make single-surface interface with hole 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 # make opening fan 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 # make closing fan 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 # make quad strip 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 # pre-populate ring dict with ring_point_indices 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
499 -def generate_surface_of_revolution(rz_points, 500 closure=EXT_OPEN, radial_subdiv=16, 501 point_idx_offset=0, ring_point_indices=None, 502 ring_markers=None, rz_closure_marker=0):
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