hedge.mesh

1 """Mesh topology representation.""" 2 3 from __future__ import division 4 5 __copyright__ = "Copyright (C) 2007 Andreas Kloeckner" 6 7 __license__ = """ 8 This program is free software: you can redistribute it and/or modify 9 it under the terms of the GNU General Public License as published by 10 the Free Software Foundation, either version 3 of the License, or 11 (at your option) any later version. 12 13 This program is distributed in the hope that it will be useful, 14 but WITHOUT ANY WARRANTY; without even the implied warranty of 15 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 16 GNU General Public License for more details. 17 18 You should have received a copy of the GNU General Public License 19 along with this program. If not, see U{http://www.gnu.org/licenses/}. 20 """ 21 22 23 24 25 import pytools 26 import numpy 27 28 # make sure AffineMap monkeypatch happens 29 import hedge.tools

30 31 32 33 34 -class TAG_NONE(object):

35 """A boundary or volume tag representing an empty boundary or volume.""" 36 pass

37 -class TAG_ALL(object):

38 """A boundary or volume tag representing the entire boundary or volume. 39 40 In the case of the boundary, TAG_ALL does not include rank boundaries, 41 or, more generally, anything tagged with TAG_NO_BOUNDARY.""" 42 pass

43 -class TAG_REALLY_ALL(object):

44 """A boundary tag representing the entire boundary. 45 46 Unlike L{TAG_ALL}, this includes rank boundaries, 47 or, more generally, everything tagged with TAG_NO_BOUNDARY.""" 48 pass

49 -class TAG_NO_BOUNDARY(object):

50 """A boundary tag indicating that this edge should not fall under TAG_ALL.""" 51 pass

52 -class TAG_RANK_BOUNDARY(object):

53 """A boundary tag indicating the boundary with a neighboring rank."""

54 - def __init__(self, rank):

55 self.rank = rank

56

57 - def __repr__(self):

58 return "TAG_RANK_BOUNDARY(%d)" % self.rank

59

60 - def __eq__(self, other):

61 return isinstance(other, TAG_RANK_BOUNDARY) and self.rank == other.rank

62

63 - def __ne__(self, other):

64 return not self.__eq__(other)

65

66 - def __hash__(self):

67 return 0xaffe ^ hash(self.rank)

68

69 70 71 72 73 74 -def make_element(class_, id, vertex_indices, all_vertices):

75 vertices = [all_vertices[v] for v in vertex_indices] 76 map = class_.get_map_unit_to_global(vertices) 77 new_vertex_indices = \ 78 class_._reorder_vertices(vertex_indices, vertices, map) 79 if new_vertex_indices: 80 vertex_indices = new_vertex_indices 81 vertices = [all_vertices[v] for v in vertex_indices] 82 map = class_.get_map_unit_to_global(vertices) 83 84 face_normals, face_jacobians = \ 85 class_.face_normals_and_jacobians(vertices, map) 86 87 return class_(id, vertex_indices, map, map.inverted(), 88 face_normals, face_jacobians)

89

90 91 92 93 -class Element(pytools.Record):

94 __slots__ = ["id", "vertex_indices", "map", "inverse_map", "face_normals", 95 "face_jacobians"] 96

97 - def __init__(self, id, vertex_indices, map, inverse_map, face_normals, 98 face_jacobians):

99 pytools.Record.__init__(self, locals())

100 101 @staticmethod

102 - def _reorder_vertices(vertex_indices, vertices):

103 return vertex_indices

104

105 - def bounding_box(self, vertices):

106 my_verts = numpy.array([vertices[vi] for vi in self.vertex_indices]) 107 return numpy.min(my_verts, axis=0), numpy.max(my_verts, axis=0)

108

109 - def centroid(self, vertices):

110 my_verts = numpy.array([vertices[vi] for vi in self.vertex_indices]) 111 return numpy.average(my_verts, axis=0)

112

113 114 115 116 117 118 -class SimplicialElement(Element):

119 __slots__ = [] 120 121 @property

122 - def faces(self):

123 return self.face_vertices(self.vertex_indices)

124 125 @classmethod

126 - def get_map_unit_to_global(cls, vertices):

127 """Return an affine map that maps the unit coordinates of the reference 128 element to a global element at a location given by its `vertices'. 129 """ 130 from hedge._internal import get_simplex_map_unit_to_global 131 return get_simplex_map_unit_to_global(cls.dimensions, vertices)

132

133 - def contains_point(self, x):

134 unit_coords = self.inverse_map(x) 135 for xi in unit_coords: 136 if xi < -1: 137 return False 138 return sum(unit_coords) <= -(self.dimensions-2)

139

140 141 142 143 -class Interval(SimplicialElement):

144 dimensions = 1 145 @staticmethod

146 - def face_vertices(vertices):

147 return [(vertices[0],), (vertices[1],) ]

148 149 @classmethod

150 - def _reorder_vertices(cls, vertex_indices, vertices, map):

151 vi = vertex_indices 152 if vertices[0][0] > vertices[1][0]: # make sure we're ordered left-right 153 return (vi[1], vi[0]) 154 else: 155 return None

156 157 @staticmethod

158 - def face_normals_and_jacobians(vertices, affine_map):

159 """Compute the normals and face jacobians of the unit element 160 transformed according to `affine_map'. 161 162 Returns a pair of lists [normals], [jacobians]. 163 """ 164 if affine_map.jacobian() < 0: 165 return [ 166 numpy.array([1], dtype=float), 167 numpy.array([-1], dtype=float) 168 ], [1, 1] 169 else: 170 return [ 171 numpy.array([-1], dtype=float), 172 numpy.array([1], dtype=float) 173 ], [1, 1]

174

175 176 177 178 -class Triangle(SimplicialElement):

179 dimensions = 2 180 181 __slots__ = [] 182 183 @staticmethod

184 - def face_vertices(vertices):

185 return [(vertices[0], vertices[1]), 186 (vertices[1], vertices[2]), 187 (vertices[0], vertices[2]) 188 ]

189 190 @classmethod

191 - def _reorder_vertices(cls, vertex_indices, vertices, map):

192 vi = vertex_indices 193 if map.jacobian() > 0: 194 return (vi[0], vi[2], vi[1]) 195 else: 196 return None

197 198 @staticmethod

199 - def face_normals_and_jacobians(vertices, affine_map):

200 """Compute the normals and face jacobians of the unit element 201 transformed according to `affine_map'. 202 203 Returns a pair of lists [normals], [jacobians]. 204 """ 205 from hedge.tools import sign 206 207 m = affine_map.matrix 208 orient = sign(affine_map.jacobian()) 209 face1 = m[:,1] - m[:,0] 210 raw_normals = [ 211 orient*numpy.array([m[1,0], -m[0,0]]), 212 orient*numpy.array([face1[1], -face1[0]]), 213 orient*numpy.array([-m[1,1], m[0,1]]), 214 ] 215 216 face_lengths = [numpy.linalg.norm(fn) for fn in raw_normals] 217 return [n/fl for n, fl in zip(raw_normals, face_lengths)], \ 218 face_lengths

219

220 221 222 223 -class Tetrahedron(SimplicialElement):

224 dimensions = 3 225 226 __slots__ = [] 227

228 - def _face_vertices(vertices):

229 return [(vertices[0],vertices[1],vertices[2]), 230 (vertices[0],vertices[1],vertices[3]), 231 (vertices[0],vertices[2],vertices[3]), 232 (vertices[1],vertices[2],vertices[3]), 233 ]

234 face_vertices = staticmethod(_face_vertices) 235 236 face_vertex_numbers = _face_vertices([0,1,2,3]) 237 238 @classmethod

239 - def _reorder_vertices(cls, vertex_indices, vertices, map):

240 vi = vertex_indices 241 if map.jacobian() > 0: 242 return (vi[0], vi[1], vi[3], vi[2]) 243 else: 244 return None

245 246 @classmethod

247 - def face_normals_and_jacobians(cls, vertices, affine_map):

248 """Compute the normals and face jacobians of the unit element 249 transformed according to `affine_map'. 250 251 Returns a pair of lists [normals], [jacobians]. 252 """ 253 from hedge._internal import tetrahedron_fj_and_normal 254 from hedge.tools import sign 255 256 return tetrahedron_fj_and_normal( 257 sign(affine_map.jacobian()), 258 cls.face_vertex_numbers, 259 vertices)

260

261 262 263 264 -class Mesh(pytools.Record):

265 """Information about the geometry and connectivity of a finite 266 element mesh. (Note: no information about the discretization 267 is stored here.) 268 269 @ivar points: list of Pylinear vectors of node coordinates 270 @ivar elements: list of Element instances 271 @ivar interfaces: a list of pairs:: 272 273 ((element instance 1, face index 1), (element instance 2, face index 2)) 274 275 enumerating elements bordering one another. The relation "element 1 touches 276 element 2" is always reflexive, but this list will only contain one entry 277 per element pair. 278 @ivar tag_to_boundary: a mapping of the form:: 279 boundary_tag -> [(element instance, face index)]) 280 281 The boundary tag TAG_NONE always refers to an empty boundary. 282 The boundary tag TAG_ALL always refers to the entire boundary. 283 @ivar tag_to_elements: a mapping of the form 284 element_tag -> [element instances] 285 286 The boundary tag TAG_NONE always refers to an empty domain. 287 The boundary tag TAG_ALL always refers to the entire domain. 288 @ivar periodicity: A list of tuples (minus_tag, plus_tag) or None 289 indicating the tags of the boundaries to be matched together 290 as periodic. There is one tuple per axis, so that for example 291 a 3D mesh has three tuples. 292 @ivar periodic_opposite_faces: a mapping of the form:: 293 (face_vertex_indices) -> 294 (opposite_face_vertex_indices), axis 295 296 This maps a face to its periodicity-induced opposite. 297 298 @ivar periodic_opposite_vertices: a mapping of the form:: 299 vertex_index -> [(opposite_vertex_index, axis), ...] 300 301 This maps one vertex to a list of its periodicity-induced 302 opposites. 303 """ 304

305 - def both_interfaces(self):

306 for face1, face2 in self.interfaces: 307 yield face1, face2 308 yield face2, face1

309 310 @property

311 - def dimensions(self):

312 return self.points.shape[1]

313

314 - def bounding_box(self):

315 try: 316 return self._bounding_box 317 except AttributeError: 318 self._bounding_box = ( 319 numpy.min(self.points, axis=0), 320 numpy.max(self.points, axis=0), 321 ) 322 return self._bounding_box

323

324 - def element_adjacency_graph(self):

325 """Return a dictionary mapping each element id to a 326 list of adjacent element ids. 327 """ 328 adjacency = {} 329 for (e1, f1), (e2, f2) in self.interfaces: 330 adjacency.setdefault(e1.id, set()).add(e2.id) 331 adjacency.setdefault(e2.id, set()).add(e1.id) 332 return adjacency

333

334 335 336 337 338 -def _build_mesh_data_dict(points, elements, boundary_tagger, periodicity, is_rankbdry_face):

339 # create face_map, which is a mapping of 340 # (vertices on a face) -> 341 # [(element, face_idx) for elements bordering that face] 342 face_map = {} 343 for el in elements: 344 for fid, face_vertices in enumerate(el.faces): 345 face_map.setdefault(frozenset(face_vertices), []).append((el, fid)) 346 347 # build non-periodic connectivity structures 348 interfaces = [] 349 tag_to_boundary = { 350 TAG_NONE: [], 351 TAG_ALL: [], 352 TAG_REALLY_ALL: [], 353 } 354 355 all_tags = set([TAG_ALL, TAG_REALLY_ALL]) 356 357 for face_vertices, els_faces in face_map.iteritems(): 358 if len(els_faces) == 2: 359 interfaces.append(els_faces) 360 elif len(els_faces) == 1: 361 el_face = el, face = els_faces[0] 362 tags = set(boundary_tagger(face_vertices, el, face, points)) - all_tags 363 364 if isinstance(tags, str): 365 from warnings import warn 366 warn("Received string as tag list") 367 368 for btag in tags: 369 tag_to_boundary.setdefault(btag, []) \ 370 .append(el_face) 371 372 if TAG_NO_BOUNDARY not in tags: 373 # TAG_NO_BOUNDARY is used to mark rank interfaces 374 # as not being part of the boundary 375 tag_to_boundary[TAG_ALL].append(el_face) 376 377 tag_to_boundary[TAG_REALLY_ALL].append(el_face) 378 else: 379 raise RuntimeError, "face can at most border two elements" 380 381 # add periodicity-induced connectivity 382 from pytools import flatten, reverse_dictionary 383 384 periodic_opposite_faces = {} 385 periodic_opposite_vertices = {} 386 387 for tag_bdries in tag_to_boundary.itervalues(): 388 assert len(set(tag_bdries)) == len(tag_bdries) 389 390 for axis, axis_periodicity in enumerate(periodicity): 391 if axis_periodicity is not None: 392 # find faces on +-axis boundaries 393 minus_tag, plus_tag = axis_periodicity 394 minus_faces = tag_to_boundary.get(minus_tag, []) 395 plus_faces = tag_to_boundary.get(plus_tag, []) 396 397 # find vertex indices and points on these faces 398 minus_vertex_indices = list(set(flatten(el.faces[face] 399 for el, face in minus_faces))) 400 plus_vertex_indices = list(set(flatten(el.faces[face] 401 for el, face in plus_faces))) 402 403 minus_z_points = [points[pi] for pi in minus_vertex_indices] 404 plus_z_points = [points[pi] for pi in plus_vertex_indices] 405 406 # find a mapping from -axis to +axis vertices 407 from hedge.tools import find_matching_vertices_along_axis 408 409 minus_to_plus, not_found = find_matching_vertices_along_axis( 410 axis, minus_z_points, plus_z_points, 411 minus_vertex_indices, plus_vertex_indices) 412 plus_to_minus = reverse_dictionary(minus_to_plus) 413 414 for a, b in minus_to_plus.iteritems(): 415 periodic_opposite_vertices.setdefault(a, []).append((b, axis)) 416 periodic_opposite_vertices.setdefault(b, []).append((a, axis)) 417 418 # establish face connectivity 419 for minus_face in minus_faces: 420 minus_el, minus_fi = minus_face 421 minus_fvi = minus_el.faces[minus_fi] 422 423 try: 424 mapped_plus_fvi = tuple(minus_to_plus[i] for i in minus_fvi) 425 plus_faces = face_map[frozenset(mapped_plus_fvi)] 426 assert len(plus_faces) == 1 427 except KeyError: 428 # is our periodic counterpart is in a different mesh clump? 429 if is_rankbdry_face(minus_face): 430 # if so, cool. parallel handler will take care of it. 431 continue 432 else: 433 # if not, bad. 434 raise 435 436 plus_face = plus_faces[0] 437 interfaces.append([minus_face, plus_face]) 438 439 plus_el, plus_fi = plus_face 440 plus_fvi = plus_el.faces[plus_fi] 441 442 mapped_minus_fvi = tuple(plus_to_minus[i] for i in plus_fvi) 443 444 # periodic_opposite_faces maps face vertex tuples from 445 # one end of the periodic domain to the other, while 446 # correspondence between each entry 447 448 periodic_opposite_faces[minus_fvi] = mapped_plus_fvi, axis 449 periodic_opposite_faces[plus_fvi] = mapped_minus_fvi, axis 450 451 tag_to_boundary[TAG_ALL].remove(plus_face) 452 tag_to_boundary[TAG_ALL].remove(minus_face) 453 454 tag_to_boundary[TAG_REALLY_ALL].remove(plus_face) 455 tag_to_boundary[TAG_REALLY_ALL].remove(minus_face) 456 457 return { 458 "interfaces": interfaces, 459 "tag_to_boundary": tag_to_boundary, 460 "periodicity": periodicity, 461 "periodic_opposite_faces": periodic_opposite_faces, 462 "periodic_opposite_vertices": periodic_opposite_vertices, 463 }

464

465 466 467 468 -def make_conformal_mesh(points, elements, 469 boundary_tagger=(lambda fvi, el, fn, all_v: []), 470 element_tagger=(lambda el, all_v: []), 471 periodicity=None, 472 _is_rankbdry_face=(lambda (el, face): False), 473 ):

474 """Construct a simplical mesh. 475 476 Face indices follow the convention for the respective element, 477 such as Triangle or Tetrahedron, in this module. 478 479 @param points: an iterable of vertex coordinates, given as vectors. 480 @param elements: an iterable of tuples of indices into points, 481 giving element endpoints. 482 @param boundary_tagger: a function that takes the arguments 483 C{(set_of_face_vertex_indices, element, face_number, all_vertices)} 484 It returns a list of tags that apply to this surface. 485 @param element_tagger: a function that takes the arguments 486 (element, all_vertices) and returns the a list of tags that apply 487 to that element. 488 @param periodicity: either None or is a list of tuples 489 just like the one documented for the C{periodicity} 490 member of class L{Mesh}. 491 @param _is_rankbdry_face: an implementation detail, 492 should not be used from user code. It is a function 493 returning whether a given face identified by 494 C{(element instance, face_nr)} is cut by a parallel 495 mesh partition. 496 """ 497 if len(points) == 0: 498 raise ValueError, "mesh contains no points" 499 500 dim = len(points[0]) 501 if dim == 1: 502 el_class = Interval 503 elif dim == 2: 504 el_class = Triangle 505 elif dim == 3: 506 el_class = Tetrahedron 507 else: 508 raise ValueError, "%d-dimensional meshes are unsupported" % dim 509 510 # build points and elements 511 new_points = numpy.array(points, dtype=float, order="C") 512 513 element_objs = [make_element(el_class, id, vert_indices, new_points) 514 for id, vert_indices in enumerate(elements)] 515 516 # tag elements 517 tag_to_elements = {TAG_NONE: [], TAG_ALL: []} 518 for el in element_objs: 519 for el_tag in element_tagger(el, new_points): 520 tag_to_elements.setdefault(el_tag, []).append(el) 521 tag_to_elements[TAG_ALL].append(el) 522 523 # build connectivity 524 if periodicity is None: 525 periodicity = dim*[None] 526 assert len(periodicity) == dim 527 528 mdd = _build_mesh_data_dict( 529 new_points, element_objs, boundary_tagger, periodicity, _is_rankbdry_face) 530 mdd["tag_to_elements"] = tag_to_elements 531 return ConformalMesh(new_points, element_objs, **mdd)

532

533 534 535 536 -class ConformalMesh(Mesh):

537 """A mesh whose elements' faces exactly match up with one another. 538 539 See the Mesh class for data members provided by this class. 540 """ 541

542 - def __init__(self, points, elements, interfaces, tag_to_boundary, tag_to_elements, 543 periodicity, periodic_opposite_faces, periodic_opposite_vertices):

544 """This constructor is for internal use only. Use L{make_conformal_mesh} instead. 545 """ 546 Mesh.__init__(self, locals())

547

548 - def get_reorder_oldnumbers(self, method):

549 if method == "cuthill": 550 from hedge.tools import cuthill_mckee 551 return cuthill_mckee(self.element_adjacency_graph()) 552 else: 553 raise ValueError, "invalid mesh reorder method"

554

555 - def reordered_by(self, method):

556 old_numbers = self.get_reorder_oldnumbers(method) 557 return self.reordered(old_numbers)

558

559 - def reordered(self, old_numbers):

560 """Return a copy of this C{Mesh} whose elements are 561 reordered using such that for each element C{i}, 562 C{old_numbers[i]} gives the previous number of that 563 element. 564 """ 565 566 elements = [self.elements[old_numbers[i]].copy(id=i) 567 for i in range(len(self.elements))] 568 569 old2new_el = dict( 570 (self.elements[old_numbers[i]], new_el) 571 for i, new_el in enumerate(elements) 572 ) 573 574 # sort interfaces by element id -- this is actually the most important part 575 def face_cmp(face1, face2): 576 (face1_el1, _), (face1_el2, _) = face1 577 (face2_el1, _), (face2_el2, _) = face2 578 579 return cmp( 580 min(face1_el1.id, face1_el2.id), 581 min(face2_el1.id, face2_el2.id))

582 583 interfaces = [ 584 ((old2new_el[e1], f1), (old2new_el[e2], f2)) 585 for (e1,f1), (e2,f2) in self.interfaces] 586 interfaces.sort(face_cmp) 587 588 tag_to_boundary = dict( 589 (tag, [(old2new_el[old_el], fnr) for old_el, fnr in elfaces]) 590 for tag, elfaces in self.tag_to_boundary.iteritems()) 591 592 tag_to_elements = dict( 593 (tag, [old2new_el[old_el] for old_el in tag_els]) 594 for tag, tag_els in self.tag_to_elements.iteritems()) 595 596 return ConformalMesh( 597 self.points, elements, interfaces, 598 tag_to_boundary, tag_to_elements, self.periodicity, 599 self.periodic_opposite_faces, self.periodic_opposite_vertices 600 )

601

602 603 604 605 -def check_bc_coverage(mesh, bc_tags, incomplete_ok=False):

606 """Verify boundary condition coverage. 607 608 Given a list of boundary tags as C{bc_tags}, this function verifies 609 that 610 1. the union of all these boundaries gives the complete boundary, 611 2. all these boundaries are disjoint. 612 613 @arg incomplete_ok: Do not report an error if some faces are not covered 614 by the boundary conditions. 615 """ 616 617 bdry_to_tag = {} 618 all_bdry_faces = mesh.tag_to_boundary[TAG_ALL] 619 bdry_face_countdown = len(all_bdry_faces) 620 621 for tag in bc_tags: 622 for el_face in mesh.tag_to_boundary.get(tag, []): 623 if el_face in bdry_to_tag: 624 raise RuntimeError, "Duplicate BCs %s found on (el=%d,face=%d)" % ( 625 [bdry_to_tag[el_face], tag],el_face[0].id, el_face[1]) 626 else: 627 bdry_to_tag[el_face] = tag 628 bdry_face_countdown -= 1 629 630 if bdry_face_countdown > 0 and not incomplete_ok: 631 raise RuntimeError, "No BCs on faces %s" % ( 632 set(all_bdry_faces)-set(bdry_to_tag.keys())) 633 elif bdry_face_countdown < 0: 634 raise RuntimeError, "More BCs were assigned than boundary faces are present " \ 635 "(did something screw up your periodicity?)"

636

637 638 639 640 -class MeshPyFaceMarkerLookup:

641 - def __init__(self, meshpy_output):

642 self.fvi2fm = dict((frozenset(fvi), marker) for fvi, marker in 643 zip(meshpy_output.facets, meshpy_output.facet_markers))

644

645 - def __call__(self, fvi):

646 return self.fvi2fm[frozenset(fvi)]

647

648 649 650 651 # mesh producers for simple geometries ---------------------------------------- 652 -def make_1d_mesh(points, left_tag=None, right_tag=None, periodic=False, 653 boundary_tagger=None, element_tagger=None):

654 def force_array(pt): 655 if not isinstance(pt, numpy.ndarray): 656 return numpy.array([pt]) 657 else: 658 return pt

659 660 def my_boundary_tagger(fvi, el, fn, all_v): 661 if el.face_normals[fn][0] < 0: 662 return [left_tag] 663 else: 664 return [right_tag] 665 666 kwargs = {} 667 if element_tagger is not None: 668 kwargs["element_tagger"] = element_tagger 669 670 if periodic: 671 left_tag = "x_minus" 672 right_tag = "x_plus" 673 return make_conformal_mesh( 674 [force_array(pt) for pt in points], 675 [(i,i+1) for i in range(len(points)-1)], 676 periodicity=[("x_minus", "x_plus")], 677 boundary_tagger=my_boundary_tagger, 678 **kwargs) 679 else: 680 return make_conformal_mesh( 681 [force_array(pt) for pt in points], 682 [(i,i+1) for i in range(len(points)-1)], 683 boundary_tagger=boundary_tagger or my_boundary_tagger, 684 **kwargs) 685

686 687 688 689 690 -def make_uniform_1d_mesh(a, b, el_count, left_tag=None, right_tag=None, periodic=False, 691 boundary_tagger=None):

692 dx = (b-a)/el_count 693 return make_1d_mesh( 694 [a+dx*i for i in range(el_count+1)], 695 left_tag=left_tag, 696 right_tag=right_tag, 697 periodic=periodic, 698 boundary_tagger=boundary_tagger)

699

700 701 702 703 -def make_single_element_mesh(a=-0.5, b=0.5, 704 boundary_tagger=(lambda vertices, face_indices: [])):

705 n = 2 706 node_dict = {} 707 points = [] 708 points_1d = numpy.linspace(a, b, n) 709 for j in range(n): 710 for i in range(n): 711 node_dict[i,j] = len(points) 712 points.append(numpy.array([points_1d[i], points_1d[j]])) 713 714 elements = [( 715 node_dict[1,1], 716 node_dict[0,1], 717 node_dict[1,0], 718 )] 719 720 boundary_faces = [(3,1), (1,2), (2,3)] 721 722 boundary_tags = dict( 723 (frozenset(seg), 724 boundary_tagger(points, seg)) 725 for seg in boundary_faces) 726 727 return make_conformal_mesh( 728 points, 729 elements, 730 boundary_tags)

731

732 733 734 735 -def make_regular_rect_mesh(a=(0,0), b=(1,1), n=(5,5), periodicity=None, 736 boundary_tagger=(lambda fvi, el, fn, all_v: [])):

737 """Create a semi-structured rectangular mesh. 738 739 @arg a: the lower left hand point of the rectangle 740 @arg b: the upper right hand point of the rectangle 741 @arg n: a tuple of integers indicating the total number of points 742 on [a,b]. 743 @arg periodicity: either None, or a tuple of bools specifying whether 744 the mesh is to be periodic in x and y. 745 """ 746 if min(n) < 2: 747 raise ValueError("need at least two points in each direction") 748 749 node_dict = {} 750 points = [] 751 points_1d = [numpy.linspace(a_i, b_i, n_i) 752 for a_i, b_i, n_i in zip(a, b, n)] 753 754 for j in range(n[1]): 755 for i in range(n[0]): 756 node_dict[i,j] = len(points) 757 points.append(numpy.array([points_1d[0][i], points_1d[1][j]])) 758 759 elements = [] 760 761 if periodicity is None: 762 periodicity = (False, False) 763 764 axes = ["x", "y"] 765 mesh_periodicity = [] 766 periodic_tags = set() 767 for i, axis in enumerate(axes): 768 if periodicity[i]: 769 minus_tag = "minus_"+axis 770 plus_tag = "plus_"+axis 771 mesh_periodicity.append((minus_tag, plus_tag)) 772 periodic_tags.add(minus_tag) 773 periodic_tags.add(plus_tag) 774 else: 775 mesh_periodicity.append(None) 776 777 fvi2fm = {} 778 779 for i in range(n[0]-1): 780 for j in range(n[1]-1): 781 782 # c--d 783 # | | 784 # a--b 785 786 a = node_dict[i,j] 787 b = node_dict[i+1,j] 788 c = node_dict[i,j+1] 789 d = node_dict[i+1,j+1] 790 791 elements.append((a,b,c)) 792 elements.append((d,c,b)) 793 794 if i == 0: fvi2fm[frozenset((a,c))] = "minus_x" 795 if i == n[0]-2: fvi2fm[frozenset((b,d))] = "plus_x" 796 if j == 0: fvi2fm[frozenset((a,b))] = "minus_y" 797 if j == n[1]-2: fvi2fm[frozenset((c,d))] = "plus_y" 798 799 def wrapped_boundary_tagger(fvi, el, fn, all_v): 800 btag = fvi2fm[frozenset(fvi)] 801 if btag in periodic_tags: 802 return [btag] 803 else: 804 return [btag] + boundary_tagger(fvi, el, fn, all_v)

805 806 return make_conformal_mesh(points, elements, wrapped_boundary_tagger, 807 periodicity=mesh_periodicity) 808

809 810 811 812 -def make_centered_regular_rect_mesh(a=(0,0), b=(1,1), n=(5,5), periodicity=None, 813 post_refine_factor=1, boundary_tagger=(lambda fvi, el, fn, all_v: [])):

814 """Create a semi-structured rectangular mesh. 815 816 @arg a: the lower left hand point of the rectangle 817 @arg b: the upper right hand point of the rectangle 818 @arg n: a tuple of integers indicating the total number of points 819 on [a,b]. 820 @arg periodicity: either None, or a tuple of bools specifying whether 821 the mesh is to be periodic in x and y. 822 """ 823 if min(n) < 2: 824 raise ValueError("need at least two points in each direction") 825 826 node_dict = {} 827 centered_node_dict = {} 828 points = [] 829 points_1d = [numpy.linspace(a_i, b_i, n_i) 830 for a_i, b_i, n_i in zip(a, b, n)] 831 dx = (numpy.array(b, dtype=numpy.float64) 832 - numpy.array(a, dtype=numpy.float64)) / (numpy.array(n)-1) 833 half_dx = dx/2 834 835 for j in range(n[1]): 836 for i in range(n[0]): 837 node_dict[i,j] = len(points) 838 points.append(numpy.array([points_1d[0][i], points_1d[1][j]])) 839 840 centered_points = [] 841 for j in range(n[1]-1): 842 for i in range(n[0]-1): 843 centered_node_dict[i,j] = len(points) 844 points.append(numpy.array([points_1d[0][i], points_1d[1][j]]) + half_dx) 845 846 elements = [] 847 848 if periodicity is None: 849 periodicity = (False, False) 850 851 axes = ["x", "y"] 852 mesh_periodicity = [] 853 periodic_tags = set() 854 for i, axis in enumerate(axes): 855 if periodicity[i]: 856 minus_tag = "minus_"+axis 857 plus_tag = "plus_"+axis 858 mesh_periodicity.append((minus_tag, plus_tag)) 859 periodic_tags.add(minus_tag) 860 periodic_tags.add(plus_tag) 861 else: 862 mesh_periodicity.append(None) 863 864 fvi2fm = {} 865 866 for i in range(n[0]-1): 867 for j in range(n[1]-1): 868 869 # c---d 870 # |\ /| 871 # | m | 872 # |/ \| 873 # a---b 874 875 a = node_dict[i,j] 876 b = node_dict[i+1,j] 877 c = node_dict[i,j+1] 878 d = node_dict[i+1,j+1] 879 880 m = centered_node_dict[i,j] 881 882 elements.append((a,b,m)) 883 elements.append((b,d,m)) 884 elements.append((d,c,m)) 885 elements.append((c,a,m)) 886 887 if i == 0: fvi2fm[frozenset((a,c))] = "minus_x" 888 if i == n[0]-2: fvi2fm[frozenset((b,d))] = "plus_x" 889 if j == 0: fvi2fm[frozenset((a,b))] = "minus_y" 890 if j == n[1]-2: fvi2fm[frozenset((c,d))] = "plus_y" 891 892 def wrapped_boundary_tagger(fvi, el, fn, all_v): 893 btag = fvi2fm[frozenset(fvi)] 894 if btag in periodic_tags: 895 return [btag] 896 else: 897 return [btag] + boundary_tagger(fvi, el, fn, all_v)

898 899 if post_refine_factor > 1: 900 from meshpy.tools import uniform_refine_triangles 901 points, elements, of2nf = uniform_refine_triangles( 902 points, elements, post_refine_factor) 903 old_fvi2fm = fvi2fm 904 fvi2fm = {} 905 906 for fvi, fm in old_fvi2fm.iteritems(): 907 for new_fvi in of2nf[fvi]: 908 fvi2fm[frozenset(new_fvi)] = fm 909 910 return make_conformal_mesh(points, elements, wrapped_boundary_tagger, 911 periodicity=mesh_periodicity) 912

913 914 915 916 917 -def make_regular_square_mesh(a=-0.5, b=0.5, n=5, periodicity=None, 918 boundary_tagger=(lambda fvi, el, fn, all_v: [])):

919 """Create a semi-structured square mesh. 920 921 @arg a: the lower x and y coordinate of the square 922 @arg b: the upper x and y coordinate of the square 923 @arg n: integer indicating the total number of points on [a,b]. 924 @arg periodicity: either None, or a tuple of bools specifying whether 925 the mesh is to be periodic in x and y. 926 """ 927 return make_regular_rect_mesh( 928 (a,a), (b,b), (n,n), periodicity, boundary_tagger)

929

930 931 932 933 -def finish_2d_rect_mesh(points, facets, facet_markers, marker2tag, refine_func, 934 periodicity, boundary_tagger):

935 """Semi-internal bottom-half routine for generation of rectangular 2D meshes.""" 936 import meshpy.triangle as triangle 937 938 mesh_info = triangle.MeshInfo() 939 mesh_info.set_points(points) 940 mesh_info.set_facets(facets, facet_markers) 941 942 #triangle.write_gnuplot_mesh("mesh.dat", mesh_info, True) 943 944 if periodicity is None: 945 periodicity = (False, False) 946 947 axes = ["x", "y"] 948 mesh_periodicity = [] 949 periodic_tags = set() 950 for i, axis in enumerate(axes): 951 if periodicity[i]: 952 minus_tag = "minus_"+axis 953 plus_tag = "plus_"+axis 954 mesh_periodicity.append((minus_tag, plus_tag)) 955 periodic_tags.add(minus_tag) 956 periodic_tags.add(plus_tag) 957 else: 958 mesh_periodicity.append(None) 959 960 generated_mesh = triangle.build(mesh_info, 961 refinement_func=refine_func, 962 allow_boundary_steiner=not (periodicity[0] or periodicity[1])) 963 964 fmlookup = MeshPyFaceMarkerLookup(generated_mesh) 965 966 def wrapped_boundary_tagger(fvi, el, fn, all_v): 967 btag = marker2tag[fmlookup(fvi)] 968 if btag in periodic_tags: 969 return [btag] 970 else: 971 return [btag] + boundary_tagger(fvi, el, fn, all_v)

972 973 return make_conformal_mesh( 974 generated_mesh.points, 975 generated_mesh.elements, 976 wrapped_boundary_tagger, 977 periodicity=mesh_periodicity) 978

979 980 981 982 -def _round_trip_connect(start, end):

983 for i in range(start, end): 984 yield i, i+1 985 yield end, start

986

987 988 989 990 -def make_rect_mesh(a=(0,0), b=(1,1), max_area=None, 991 boundary_tagger=(lambda fvi, el, fn, all_v: []), 992 periodicity=None, subdivisions=None, 993 refine_func=None):

994 """Create an unstructured rectangular mesh. 995 996 @arg a: the lower left hand point of the rectangle 997 @arg b: the upper right hand point of the rectangle 998 @arg max_area: maximum area of each triangle. 999 @arg periodicity: either None, or a tuple of bools specifying whether 1000 the mesh is to be periodic in x and y. 1001 @arg subdivisions: If not C{None}, this is a 2-tuple specifying 1002 the number of facet subdivisions in X and Y. 1003 @arg refine_func: A refinement function as taken by C{meshpy.triangle.build}. 1004 """ 1005 import meshpy.triangle as triangle 1006 1007 if max_area is not None: 1008 if refine_func is not None: 1009 raise ValueError, "cannot specify both refine_func and max_area" 1010 def refine_func(vertices, area): 1011 return area > max_area

1012 1013 marker2tag = { 1014 1: "minus_x", 1015 2: "minus_y", 1016 3: "plus_x", 1017 4: "plus_y", 1018 } 1019 1020 points = [a, (b[0],a[1]), b, (a[0],b[1])] 1021 facets = list(_round_trip_connect(0, 3)) 1022 facet_markers = [2,3,4,1] 1023 1024 if subdivisions is not None: 1025 points, facets, facet_markers = triangle.subdivide_facets( 1026 [subdivisions[0], subdivisions[1], 1027 subdivisions[0], subdivisions[1]], 1028 points, facets, facet_markers) 1029 1030 return finish_2d_rect_mesh(points, facets, facet_markers, marker2tag, 1031 refine_func, periodicity, boundary_tagger) 1032

1033 1034 1035 1036 1037 -def make_rect_mesh_with_corner(a=(0,0), b=(1,1), max_area=None, 1038 boundary_tagger=(lambda fvi, el, fn, all_v: []), 1039 corner_fraction=(0.3, 0.3), 1040 refine_func=None):

1041 """Create an unstructured rectangular mesh with a reentrant 1042 corner at (-x, -y). 1043 1044 @arg a: the lower left hand point of the rectangle 1045 @arg b: the upper right hand point of the rectangle 1046 @arg max_area: maximum area of each triangle. 1047 @arg refine_func: A refinement function as taken by C{meshpy.triangle.build}. 1048 @arg corner_fraction: Tuple of fraction of the width taken up by 1049 the rentrant corner. 1050 """ 1051 import meshpy.triangle as triangle 1052 1053 if max_area is not None: 1054 if refine_func is not None: 1055 raise ValueError, "cannot specify both refine_func and max_area" 1056 def refine_func(vertices, area): 1057 return area > max_area

1058 1059 marker2tag = { 1060 1: "minus_x", 1061 2: "minus_y", 1062 3: "plus_x", 1063 4: "plus_y", 1064 4: "plus_y", 1065 5: "corner_plus_y", 1066 6: "corner_plus_x", 1067 } 1068 1069 a = numpy.asarray(a) 1070 b = numpy.asarray(b) 1071 diag = b-a 1072 w = diag.copy(); w[1] = 0 1073 h = diag.copy(); h[0] = 0 1074 1075 points = [ 1076 a+h*corner_fraction[1], 1077 a+h*corner_fraction[1]+w*corner_fraction[0], 1078 a+w*corner_fraction[0], 1079 a+w, 1080 a+w+h, 1081 a+h, 1082 ] 1083 facets = list(_round_trip_connect(0, 5)) 1084 facet_markers = [5,6,2,3,4,1] 1085 1086 import meshpy.triangle as triangle 1087 mesh_info = triangle.MeshInfo() 1088 mesh_info.set_points(points) 1089 mesh_info.set_facets(facets, facet_markers) 1090 1091 generated_mesh = triangle.build(mesh_info, 1092 refinement_func=refine_func) 1093 1094 return make_conformal_mesh( 1095 generated_mesh.points, 1096 generated_mesh.elements, 1097 boundary_tagger) 1098

1099 1100 1101 1102 1103 -def make_square_mesh(a=-0.5, b=0.5, max_area=4e-3, 1104 boundary_tagger=(lambda fvi, el, fn, all_v: [])):

1105 """Create an unstructured square mesh. 1106 1107 @arg a: the lower x and y coordinate of the square 1108 @arg b: the upper x and y coordinate of the square 1109 @arg max_area: maximum area of each triangle 1110 """ 1111 return make_rect_mesh((a,a), (b,b), max_area, boundary_tagger)

1112

1113 1114 1115 1116 -def make_disk_mesh(r=0.5, faces=50, max_area=4e-3, 1117 boundary_tagger=(lambda fvi, el, fn, all_v: [])):

1118 from math import cos, sin, pi 1119 1120 def needs_refinement(vertices, area): 1121 return area > max_area

1122 1123 points = [(r*cos(angle), r*sin(angle)) 1124 for angle in numpy.linspace(0, 2*pi, faces, endpoint=False)] 1125 1126 import meshpy.triangle as triangle 1127 1128 mesh_info = triangle.MeshInfo() 1129 mesh_info.set_points(points) 1130 mesh_info.set_facets( 1131 list(_round_trip_connect(0, faces-1)), 1132 faces*[1] 1133 ) 1134 1135 generated_mesh = triangle.build(mesh_info, refinement_func=needs_refinement) 1136 1137 return make_conformal_mesh( 1138 generated_mesh.points, 1139 generated_mesh.elements, 1140 boundary_tagger) 1141

1142 1143 1144 1145 -def make_ball_mesh(r=0.5, subdivisions=10, max_volume=None, 1146 boundary_tagger=(lambda fvi, el, fn, all_v: [])):

1147 from meshpy.tet import MeshInfo, build 1148 from meshpy.geometry import make_ball 1149 1150 points, facets, facet_holestarts, facet_markers = \ 1151 make_ball(r, subdivisions) 1152 1153 mesh_info = MeshInfo() 1154 1155 mesh_info.set_points(points) 1156 mesh_info.set_facets_ex(facets, facet_holestarts, facet_markers) 1157 generated_mesh = build(mesh_info, max_volume=max_volume) 1158 1159 return make_conformal_mesh( 1160 generated_mesh.points, 1161 generated_mesh.elements, 1162 boundary_tagger)

1163

1164 1165 1166 1167 1168 1169 1170 1171 -def _make_z_periodic_mesh(points, facets, facet_holestarts, facet_markers, height, 1172 max_volume, boundary_tagger):

1173 from meshpy.tet import MeshInfo, build 1174 from meshpy.geometry import Marker 1175 1176 mesh_info = MeshInfo() 1177 mesh_info.set_points(points) 1178 mesh_info.set_facets_ex(facets, facet_holestarts, facet_markers) 1179 1180 mesh_info.pbc_groups.resize(1) 1181 pbcg = mesh_info.pbc_groups[0] 1182 1183 pbcg.facet_marker_1 = Marker.MINUS_Z 1184 pbcg.facet_marker_2 = Marker.PLUS_Z 1185 1186 pbcg.set_transform(translation=[0,0,height]) 1187 1188 def zper_boundary_tagger(fvi, el, fn, all_v): 1189 # we only ask about *boundaries* 1190 # we should not try to have the user tag 1191 # the (periodicity-induced) interior faces 1192 1193 face_marker = fvi2fm[frozenset(fvi)] 1194 1195 if face_marker == Marker.MINUS_Z: 1196 return ["minus_z"] 1197 if face_marker == Marker.PLUS_Z: 1198 return ["plus_z"] 1199 1200 result = boundary_tagger(fvi, el, fn, all_v) 1201 if face_marker == Marker.SHELL: 1202 result.append("shell") 1203 return result

1204 1205 generated_mesh = build(mesh_info, max_volume=max_volume) 1206 fvi2fm = generated_mesh.face_vertex_indices_to_face_marker 1207 1208 return make_conformal_mesh( 1209 generated_mesh.points, 1210 generated_mesh.elements, 1211 zper_boundary_tagger, 1212 periodicity=[None, None, ("minus_z", "plus_z")]) 1213

1214 1215 1216 1217 -def make_cylinder_mesh(radius=0.5, height=1, radial_subdivisions=10, 1218 height_subdivisions=1, max_volume=None, periodic=False, 1219 boundary_tagger=(lambda fvi, el, fn, all_v: [])):

1220 from meshpy.tet import MeshInfo, build 1221 from meshpy.geometry import make_cylinder 1222 1223 points, facets, facet_holestarts, facet_markers = \ 1224 make_cylinder(radius, height, radial_subdivisions, 1225 height_subdivisions) 1226 1227 assert len(facets) == len(facet_markers) 1228 1229 if periodic: 1230 return _make_z_periodic_mesh( 1231 points, facets, facet_holestarts, facet_markers, 1232 height=height, 1233 max_volume=max_volume, 1234 boundary_tagger=boundary_tagger) 1235 else: 1236 mesh_info = MeshInfo() 1237 mesh_info.set_points(points) 1238 mesh_info.set_facets_ex(facets, facet_holestarts, facet_markers) 1239 1240 generated_mesh = build(mesh_info, max_volume=max_volume) 1241 1242 return make_conformal_mesh( 1243 generated_mesh.points, 1244 generated_mesh.elements, 1245 boundary_tagger)

1246

1247 1248 1249 1250 -def make_box_mesh(a=(0,0,0),b=(1,1,1), 1251 max_volume=None, periodicity=None, 1252 boundary_tagger=(lambda fvi, el, fn, all_v: []), 1253 return_meshpy_mesh=False):

1254 """Return a mesh for a brick from the origin to `dimensions'. 1255 1256 `max_volume' specifies the maximum volume for each tetrahedron. 1257 `periodicity' is either None, or a triple of bools, indicating 1258 whether periodic BCs are to be applied along that axis. 1259 See ConformalMesh.__init__ for the meaning of boundary_tagger. 1260 1261 A few stock boundary tags are provided for easy application 1262 of boundary conditions, namely plus_[xyz] and minus_[xyz] tag 1263 the appropriate faces of the brick. 1264 """ 1265 1266 def count(iterable): 1267 result = 0 1268 for i in iterable: 1269 result += 1 1270 return result

1271 1272 from meshpy.tet import MeshInfo, build 1273 from meshpy.geometry import make_box 1274 1275 points, facets, facet_markers = make_box(a, b) 1276 1277 mesh_info = MeshInfo() 1278 mesh_info.set_points(points) 1279 mesh_info.set_facets(facets, facet_markers) 1280 1281 if periodicity is None: 1282 periodicity = (False, False, False) 1283 1284 axes = ["x", "y", "z"] 1285 1286 per_count = count(p for p in periodicity if p) 1287 mesh_info.pbc_groups.resize(per_count) 1288 pbc_group_number = 0 1289 1290 marker_to_tag = {} 1291 mesh_periodicity = [] 1292 periodic_tags = set() 1293 1294 for axis, axis_per in enumerate(periodicity): 1295 minus_marker = 1+2*axis 1296 plus_marker = 2+2*axis 1297 1298 minus_tag = "minus_"+axes[axis] 1299 plus_tag = "plus_"+axes[axis] 1300 1301 marker_to_tag[minus_marker] = minus_tag 1302 marker_to_tag[plus_marker] = plus_tag 1303 1304 if axis_per: 1305 pbcg = mesh_info.pbc_groups[pbc_group_number] 1306 pbc_group_number +=1 1307 1308 pbcg.facet_marker_1 = minus_marker 1309 pbcg.facet_marker_2 = plus_marker 1310 1311 translation = [0,0,0] 1312 translation[axis] = b[axis]-a[axis] 1313 pbcg.set_transform(translation=translation) 1314 1315 mesh_periodicity.append((minus_tag, plus_tag)) 1316 periodic_tags.add(minus_tag) 1317 periodic_tags.add(plus_tag) 1318 else: 1319 mesh_periodicity.append(None) 1320 1321 generated_mesh = build(mesh_info, max_volume=max_volume) 1322 1323 fvi2fm = generated_mesh.face_vertex_indices_to_face_marker 1324 1325 def wrapped_boundary_tagger(fvi, el, fn, all_v): 1326 face_tag = marker_to_tag[fvi2fm[frozenset(fvi)]] 1327 1328 if face_tag in periodic_tags: 1329 return [face_tag] 1330 else: 1331 return [face_tag] + boundary_tagger(fvi, el, fn, all_v) 1332 1333 result = make_conformal_mesh( 1334 generated_mesh.points, 1335 generated_mesh.elements, 1336 wrapped_boundary_tagger, 1337 periodicity=mesh_periodicity) 1338 1339 if return_meshpy_mesh: 1340 return result, generated_mesh 1341 else: 1342 return result 1343

Source Code for Module hedge.mesh