hedge.element

1 """Local function space 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 numpy 26 import numpy.linalg as la 27 import pyublas 28 from hedge.tools import AffineMap 29 import hedge._internal 30 from math import sqrt 31 from pytools import memoize_method 32 import hedge.mesh 33 34 35 36 37 __all__ = ["TriangularElement", "TetrahedralElement"]

38 39 40 41 42 -class WarpFactorCalculator:

43 """Calculator for Warburton's warp factor. 44 45 See T. Warburton, 46 "An explicit construction of interpolation nodes on the simplex" 47 Journal of Engineering Mathematics Vol 56, No 3, p. 247-262, 2006 48 """ 49

50 - def __init__(self, N):

51 from hedge.quadrature import legendre_gauss_lobatto_points 52 from hedge.interpolation import newton_interpolation_function 53 54 # Find lgl and equidistant interpolation points 55 r_lgl = legendre_gauss_lobatto_points(N) 56 r_eq = numpy.linspace(-1,1,N+1) 57 58 self.int_f = newton_interpolation_function(r_eq, r_lgl - r_eq)

59

60 - def __call__(self, x):

61 if abs(x) > 1-1e-10: 62 return 0 63 else: 64 return self.int_f(x)/(1-x**2)

65

66 67 68 69 -class FaceVertexMismatch(Exception):

70 - def __init__(self, value):

71 self.value = value

72 - def __str__(self):

73 return repr(self.value)

74

75 76 77 78 -class TriangleWarper:

79 - def __init__(self, alpha, order):

80 self.alpha = alpha 81 self.warp = WarpFactorCalculator(order) 82 83 cls = TriangularElement 84 85 from pytools import wandering_element 86 from tools import normalize 87 88 vertices = [cls.barycentric_to_equilateral(bary) 89 for bary in wandering_element(cls.dimensions+1)] 90 all_vertex_indices = range(cls.dimensions+1) 91 face_vertex_indices = cls.geometry \ 92 .face_vertices(all_vertex_indices) 93 faces_vertices = cls.geometry.face_vertices(vertices) 94 95 edgedirs = [normalize(v2-v1) for v1, v2 in faces_vertices] 96 opp_vertex_indices = [ 97 (set(all_vertex_indices) - set(fvi)).__iter__().next() 98 for fvi in face_vertex_indices] 99 100 self.loop_info = zip( 101 face_vertex_indices, 102 edgedirs, 103 opp_vertex_indices)

104

105 - def __call__(self, bp):

106 shifts = [] 107 108 from operator import add, mul 109 110 for fvi, edgedir, opp_vertex_index in self.loop_info: 111 blend = 4*reduce(mul, (bp[i] for i in fvi)) 112 warp_amount = blend*self.warp(bp[fvi[1]]-bp[fvi[0]]) \ 113 * (1 + (self.alpha*bp[opp_vertex_index])**2) 114 shifts.append(warp_amount*edgedir) 115 116 return reduce(add, shifts)

117 118 119 120 121 TriangleBasisFunction = hedge._internal.TriangleBasisFunction 122 GradTriangleBasisFunction = hedge._internal.GradTriangleBasisFunction 123 TetrahedronBasisFunction = hedge._internal.TetrahedronBasisFunction 124 GradTetrahedronBasisFunction = hedge._internal.GradTetrahedronBasisFunction

125 126 127 128 129 -class Element(object):

130 - def basis_functions(self):

131 """Get a sequence of functions that form a basis of the approximation space.""" 132 raise NotImplementedError

133

134 - def grad_basis_functions(self):

135 """Get the gradient functions of the basis_functions(), in the same order.""" 136 raise NotImplementedError

137 138 # matrices ---------------------------------------------------------------- 139 @memoize_method

140 - def vandermonde(self):

141 from hedge.polynomial import generic_vandermonde 142 143 return generic_vandermonde( 144 list(self.unit_nodes()), 145 list(self.basis_functions()))

146 147 @memoize_method

148 - def inverse_mass_matrix(self):

149 """Return the inverse of the mass matrix of the unit element 150 with respect to the nodal coefficients. Divide by the Jacobian 151 to obtain the global mass matrix. 152 """ 153 154 # see doc/hedge-notes.tm 155 v = self.vandermonde() 156 return numpy.dot(v, v.T)

157 158 @memoize_method

159 - def mass_matrix(self):

160 """Return the mass matrix of the unit element with respect 161 to the nodal coefficients. Multiply by the Jacobian to obtain 162 the global mass matrix. 163 """ 164 165 return numpy.asarray(la.inv(self.inverse_mass_matrix()), order="C")

166 167 @memoize_method

168 - def grad_vandermonde(self):

169 """Compute the Vandermonde matrices of the grad_basis_functions(). 170 Return a list of these matrices.""" 171 172 from hedge.polynomial import generic_multi_vandermonde 173 174 return generic_multi_vandermonde( 175 list(self.unit_nodes()), 176 list(self.grad_basis_functions()))

177 178 @memoize_method

179 - def differentiation_matrices(self):

180 """Return matrices that map the nodal values of a function 181 to the nodal values of its derivative in each of the unit 182 coordinate directions. 183 """ 184 185 from hedge.tools import leftsolve 186 # see doc/hedge-notes.tm 187 v = self.vandermonde() 188 return [leftsolve(v, vdiff) for vdiff in self.grad_vandermonde()]

189 190 @memoize_method

191 - def multi_face_mass_matrix(self):

192 """Return a matrix that combines the effect of multiple face 193 mass matrices applied to a vector of the shape:: 194 195 [face_1_dofs]+[face_2_dofs]+...+[face_n_dofs] 196 197 Observe that this automatically maps this vector to a volume 198 contribution. 199 """ 200 fnc = self.face_node_count() 201 202 result = numpy.zeros( 203 (self.node_count(), self.face_count()*fnc), 204 dtype=numpy.float) 205 206 fmm = self.face_mass_matrix() 207 208 for i_face, f_indices in enumerate(self.face_indices()): 209 for i_dof, f_index in enumerate(f_indices): 210 result[f_index, i_face*fnc:(i_face+1)*fnc] = fmm[i_dof] 211 212 return result

213 214 @memoize_method

215 - def lifting_matrix(self):

216 """Return a matrix that combines the effect of the inverse 217 mass matrix applied after the multi-face mass matrix to a vector 218 of the shape:: 219 220 [face_1_dofs]+[face_2_dofs]+...+[face_n_dofs] 221 222 Observe that this automatically maps this vector to a volume 223 contribution. 224 """ 225 return numpy.dot(self.inverse_mass_matrix(), 226 self.multi_face_mass_matrix())

227

228 - def find_diff_mat_permutation(self, target_idx):

229 """Find a permuation C{p} such that:: 230 231 diff_mats = self.differentiation_matrices() 232 diff_mats[0][p][:,p] == diff_mats[target_idx] 233 234 The permutation is returned as a numpy array of type intp. 235 """ 236 237 node_tups = self.node_tuples() 238 d = self.dimensions 239 240 from pytools import get_read_from_map_from_permutation 241 242 def transpose(tup, i, j): 243 l = list(tup) 244 l[i], l[j] = l[j], l[i] 245 return tuple(l)

246 247 p = numpy.array(get_read_from_map_from_permutation( 248 node_tups, 249 [transpose(nt, 0, target_idx) for nt in node_tups]), 250 dtype=numpy.intp) 251 252 dmats = self.differentiation_matrices() 253 assert la.norm(dmats[0][p][:,p]-dmats[target_idx]) < 1e-12 254 255 return p

256

257 258 259 260 261 262 -class SimplicialElement(Element):

263 # queries ----------------------------------------------------------------- 264 @property

265 - def has_facial_nodes(self):

266 return self.order > 0

267 268 # numbering --------------------------------------------------------------- 269 @memoize_method

270 - def node_count(self):

271 """Return the number of interpolation nodes in this element.""" 272 d = self.dimensions 273 o = self.order 274 from operator import mul 275 from pytools import factorial 276 return int(reduce(mul, (o+1+i for i in range(d)))/factorial(d))

277 278 @memoize_method

279 - def face_count(self):

280 return len(self.face_indices())

281 282 @memoize_method

283 - def face_node_count(self):

284 return len(self.face_indices()[0])

285 286 @memoize_method

287 - def vertex_indices(self):

288 """Return the list of the vertices' node indices.""" 289 from pytools import wandering_element 290 291 result = [] 292 293 node_tup_to_idx = dict( 294 (ituple, idx) 295 for idx, ituple in enumerate(self.node_tuples())) 296 297 vertex_tuples = [self.dimensions * (0,)] \ 298 + list(wandering_element(self.dimensions, wanderer=self.order)) 299 300 return [node_tup_to_idx[vt] for vt in vertex_tuples]

301 302 @memoize_method

303 - def face_indices(self):

304 """Return a list of face index lists. Each face index list contains 305 the local node numbers of the nodes on that face. 306 """ 307 308 node_tup_to_idx = dict( 309 (ituple, idx) 310 for idx, ituple in enumerate(self.node_tuples())) 311 312 from pytools import generate_nonnegative_integer_tuples_summing_to_at_most 313 314 enum_order_nodes_gen = generate_nonnegative_integer_tuples_summing_to_at_most( 315 self.order, self.dimensions) 316 317 faces = [[] for i in range(self.dimensions+1)] 318 319 for node_tup in enum_order_nodes_gen: 320 for face_idx in self.faces_for_node_tuple(node_tup): 321 faces[face_idx].append(node_tup_to_idx[node_tup]) 322 323 return [tuple(fi) for fi in faces]

324 325 # node wrangling ----------------------------------------------------------

326 - def equidistant_barycentric_nodes(self):

327 """Generate equidistant nodes in barycentric coordinates.""" 328 for indices in self.node_tuples(): 329 divided = tuple(i/self.order for i in indices) 330 yield (1-sum(divided),) + divided

331

332 - def equidistant_equilateral_nodes(self):

333 """Generate equidistant nodes in equilateral coordinates.""" 334 335 for bary in self.equidistant_barycentric_nodes(): 336 yield self.barycentric_to_equilateral(bary)

337

338 - def equidistant_unit_nodes(self):

339 """Generate equidistant nodes in unit coordinates.""" 340 341 for bary in self.equidistant_barycentric_nodes(): 342 yield self.equilateral_to_unit(self.barycentric_to_equilateral(bary))

343 344 @memoize_method

345 - def unit_nodes(self):

346 """Generate the warped nodes in unit coordinates (r,s,...).""" 347 return [self.equilateral_to_unit(node) 348 for node in self.equilateral_nodes()]

349 350 # basis functions ---------------------------------------------------------

351 - def generate_mode_identifiers(self):

352 """Generate a hashable objects identifying each basis function, in order. 353 354 The output from this function is required to be in the same order 355 as that of L{basis_functions} and L{grad_basis_functions}, and thereby 356 also from L{vandermonde}. 357 """ 358 from pytools import generate_nonnegative_integer_tuples_summing_to_at_most 359 360 return generate_nonnegative_integer_tuples_summing_to_at_most( 361 self.order, self.dimensions)

362 363 # time step scaling -------------------------------------------------------

364 - def dt_non_geometric_factor(self):

365 unodes = self.unit_nodes() 366 vertex_indices = self.vertex_indices() 367 return 2/3* \ 368 min(min(min( 369 la.norm(unodes[face_node_index]-unodes[vertex_index]) 370 for vertex_index in vertex_indices 371 if vertex_index != face_node_index) 372 for face_node_index in face_indices) 373 for face_indices in self.face_indices())

374

375 376 377 378 -class IntervalElementBase(SimplicialElement):

379 dimensions = 1 380 has_local_jacobians = False 381 geometry = hedge.mesh.Interval 382 383 # numbering --------------------------------------------------------------- 384 @memoize_method

385 - def node_tuples(self):

386 """Generate tuples enumerating the node indices present 387 in this element. Each tuple has a length equal to the dimension 388 of the element. The tuples constituents are non-negative integers 389 whose sum is less than or equal to the order of the element. 390 391 The order in which these nodes are generated dictates the local 392 node numbering. 393 """ 394 return [(i,) for i in range(self.order+1)]

395

396 - def faces_for_node_tuple(self, node_idx):

397 """Return the list of face indices of faces on which the node 398 represented by C{node_idx} lies. 399 """ 400 401 if node_idx == (0,): 402 if self.order == 0: 403 return [0, 1] 404 else: 405 return [0] 406 elif node_idx == (self.order,): 407 return [1] 408 else: 409 return []

410 411 # node wrangling ---------------------------------------------------------- 412 @memoize_method

413 - def get_submesh_indices(self):

414 """Return a list of tuples of indices into the node list that 415 generate a tesselation of the reference element.""" 416 417 return [(i,i+1) for i in range(self.order)]

418 419 # face operations --------------------------------------------------------- 420 @memoize_method

421 - def face_mass_matrix(self):

422 return numpy.array([[1]], dtype=float)

423 424 @staticmethod

425 - def get_face_index_shuffle_to_match(face_1_vertices, face_2_vertices):

426 if set(face_1_vertices) != set(face_2_vertices): 427 raise FaceVertexMismatch("face vertices do not match") 428 429 class IntervalFaceIndexShuffle: 430 def __hash__(self): 431 return 0x3472477

432 433 def __eq__(self, other): 434 return True

435 436 def __call__(self, indices): 437 return indices 438 439 return IntervalFaceIndexShuffle() 440

441 442 443 444 445 -class IntervalElement(IntervalElementBase):

446 """An arbitrary-order polynomial finite interval element. 447 448 Coordinate systems used: 449 ======================== 450 451 unit coordinates (r):: 452 453 ---[--------0--------]---> 454 -1 1 455 """ 456

457 - def __init__(self, order):

458 self.order = order

459 460 # node wrangling ----------------------------------------------------------

461 - def nodes(self):

462 """Generate warped nodes in unit coordinates (r,).""" 463 464 if self.order == 0: 465 return [ numpy.array([0.5]) ] 466 else: 467 from hedge.quadrature import legendre_gauss_lobatto_points 468 return [numpy.array([x]) for x in legendre_gauss_lobatto_points(self.order)]

469 470 equilateral_nodes = nodes 471 unit_nodes = nodes 472 473 # basis functions --------------------------------------------------------- 474 @memoize_method

475 - def basis_functions(self):

476 """Get a sequence of functions that form a basis of the approximation space. 477 478 The approximation space is spanned by the polynomials::: 479 480 r**i for i <= N 481 """ 482 from hedge.polynomial import LegendreFunction 483 484 class VectorLF: 485 def __init__(self, n): 486 self.lf = LegendreFunction(n)

487 488 def __call__(self, x): 489 return self.lf(x[0])

490 491 return [VectorLF(idx[0]) for idx in self.generate_mode_identifiers()] 492

493 - def grad_basis_functions(self):

494 """Get the gradient functions of the basis_functions(), in the same order.""" 495 from hedge.polynomial import DiffLegendreFunction 496 497 class DiffVectorLF: 498 def __init__(self, n): 499 self.dlf = DiffLegendreFunction(n)

500 501 def __call__(self, x): 502 return numpy.array([self.dlf(x[0])]) 503 504 return [DiffVectorLF(idx[0]) 505 for idx in self.generate_mode_identifiers()] 506 507 # time step scaling -------------------------------------------------------

508 - def dt_non_geometric_factor(self):

509 if self.order == 0: 510 return 1 511 else: 512 unodes = self.unit_nodes() 513 return la.norm(unodes[0] - unodes[1]) * 0.85

514

515 - def dt_geometric_factor(self, vertices, el):

516 return abs(el.map.jacobian())

517

518 519 520 521 -class TriangularElement(SimplicialElement):

522 """An arbitrary-order triangular finite element. 523 524 Coordinate systems used: 525 ======================== 526 527 unit coordinates (r,s):: 528 529 C 530 |\\ 531 | \\ 532 | O 533 | \\ 534 | \\ 535 A-----B 536 537 Points in unit coordinates:: 538 539 O = (0,0) 540 A = (-1,-1) 541 B = (1,-1) 542 C = (-1,1) 543 544 equilateral coordinates (x,y):: 545 546 C 547 / \\ 548 / \\ 549 / \\ 550 / O \\ 551 / \\ 552 A-----------B 553 554 Points in equilateral coordinates:: 555 556 O = (0,0) 557 A = (-1,-1/sqrt(3)) 558 B = (1,-1/sqrt(3)) 559 C = (0,2/sqrt(3)) 560 561 When global vertices are passed in, they are mapped to the 562 reference vertices A, B, C in order. 563 564 Faces are always ordered AB, BC, AC. 565 """ 566 567 # In case you were wondering: the double backslashes in the docstring 568 # are required because single backslashes only escape their subsequent 569 # newlines, and thus end up not yielding a correct docstring. 570 571 dimensions = 2 572 has_local_jacobians = False 573 geometry = hedge.mesh.Triangle 574

575 - def __init__(self, order, fancy_node_ordering=True):

576 self.order = order 577 self.fancy_node_ordering = fancy_node_ordering

578 579 # numbering --------------------------------------------------------------- 580 @memoize_method

581 - def node_tuples(self):

582 """Generate tuples enumerating the node indices present 583 in this element. Each tuple has a length equal to the dimension 584 of the element. The tuples constituents are non-negative integers 585 whose sum is less than or equal to the order of the element. 586 587 The order in which these nodes are generated dictates the local 588 node numbering. 589 """ 590 from pytools import generate_nonnegative_integer_tuples_summing_to_at_most 591 node_tups = list(generate_nonnegative_integer_tuples_summing_to_at_most( 592 self.order, self.dimensions)) 593 594 if not self.fancy_node_ordering: 595 return node_tups 596 597 faces_to_nodes = {} 598 for node_tup in node_tups: 599 faces_to_nodes.setdefault( 600 frozenset(self.faces_for_node_tuple(node_tup)), 601 []).append(node_tup) 602 603 result = [] 604 def add_face_nodes(faces): 605 result.extend(faces_to_nodes.get(frozenset(faces), []))

606 607 add_face_nodes([]) 608 add_face_nodes([0]) 609 add_face_nodes([0,1]) 610 add_face_nodes([1]) 611 add_face_nodes([1,2]) 612 add_face_nodes([2]) 613 add_face_nodes([0,2]) 614 615 assert set(result) == set(node_tups) 616 assert len(result) == len(node_tups) 617 618 return result

619

620 - def faces_for_node_tuple(self, node_tuple):

621 """Return the list of face indices of faces on which the node 622 represented by C{node_tuple} lies. 623 """ 624 m, n = node_tuple 625 626 result = [] 627 if n == 0: 628 result.append(0) 629 if n+m == self.order: 630 result.append(1) 631 if m == 0: 632 result.append(2) 633 return result

634 635 # node wrangling ---------------------------------------------------------- 636 @staticmethod

637 - def barycentric_to_equilateral((lambda1, lambda2, lambda3)):

638 """Return the equilateral (x,y) coordinate corresponding 639 to the barycentric coordinates (lambda1..lambdaN).""" 640 641 # reflects vertices in equilateral coordinates 642 return numpy.array([ 643 (-lambda1 +lambda2 ), 644 (-lambda1 -lambda2 +2*lambda3)/sqrt(3.0)])

645 646 # see doc/hedge-notes.tm 647 equilateral_to_unit = AffineMap( 648 numpy.array([[1,-1/sqrt(3)], [0,2/sqrt(3)]]), 649 numpy.array([-1/3,-1/3])) 650

651 - def equilateral_nodes(self):

652 """Generate warped nodes in equilateral coordinates (x,y).""" 653 654 # port of Hesthaven/Warburton's Nodes2D routine 655 # note that the order of the barycentric coordinates is changed 656 # match the order of the equilateral vertices 657 658 # Not much is left of the original routine--it was very redundant. 659 # The test suite still contains the original code and verifies this 660 # one against it. 661 662 # Set optimized parameter alpha, depending on order N 663 alpha_opt = [0.0000, 0.0000, 1.4152, 0.1001, 0.2751, 0.9800, 1.0999, 664 1.2832, 1.3648, 1.4773, 1.4959, 1.5743, 1.5770, 1.6223, 1.6258] 665 666 try: 667 alpha = alpha_opt[self.order-1] 668 except IndexError: 669 alpha = 5/3 670 671 warp = TriangleWarper(alpha, self.order) 672 673 for bp in self.equidistant_barycentric_nodes(): 674 yield self.barycentric_to_equilateral(bp) + warp(bp)

675 676 @memoize_method

677 - def get_submesh_indices(self):

678 """Return a list of tuples of indices into the node list that 679 generate a tesselation of the reference element.""" 680 681 node_dict = dict( 682 (ituple, idx) 683 for idx, ituple in enumerate(self.node_tuples())) 684 685 result = [] 686 for i, j in self.node_tuples(): 687 if i+j < self.order: 688 result.append( 689 (node_dict[i,j], node_dict[i+1,j], node_dict[i,j+1])) 690 if i+j < self.order-1: 691 result.append( 692 (node_dict[i+1,j+1], node_dict[i,j+1], node_dict[i+1,j])) 693 return result

694 695 # basis functions --------------------------------------------------------- 696 @memoize_method

697 - def basis_functions(self):

698 """Get a sequence of functions that form a basis of the approximation space. 699 700 The approximation space is spanned by the polynomials::: 701 702 r**i * s**j for i+j <= N 703 """ 704 return [TriangleBasisFunction(*idx) for idx in 705 self.generate_mode_identifiers()]

706

707 - def grad_basis_functions(self):

708 """Get the gradient functions of the basis_functions(), in the same order.""" 709 return [GradTriangleBasisFunction(*idx) for idx in 710 self.generate_mode_identifiers()]

711 712 # face operations --------------------------------------------------------- 713 @memoize_method

714 - def face_mass_matrix(self):

715 from hedge.polynomial import legendre_vandermonde 716 unodes = self.unit_nodes() 717 face_vandermonde = legendre_vandermonde( 718 [unodes[i][0] for i in self.face_indices()[0]], 719 self.order) 720 721 return numpy.asarray( 722 la.inv( 723 numpy.dot(face_vandermonde, face_vandermonde.T)), 724 order="C")

725 726 @staticmethod

727 - def get_face_index_shuffle_to_match(face_1_vertices, face_2_vertices):

728 if set(face_1_vertices) != set(face_2_vertices): 729 raise FaceVertexMismatch("face vertices do not match") 730 731 class TriangleFaceIndexShuffle: 732 def __init__(self, operations): 733 self.operations = operations

734 735 def __hash__(self): 736 return hash(self.operations) 737 738 def __eq__(self, other): 739 return self.operations == other.operations 740 741 def __call__(self, indices): 742 for op in self.operations: 743 if op == "flip": 744 indices = indices[::-1] 745 else: 746 raise RuntimeError, "invalid operation" 747 return indices 748 749 750 if face_1_vertices != face_2_vertices: 751 assert face_1_vertices[::-1] == face_2_vertices 752 return TriangleFaceIndexShuffle(("flip",)) 753 else: 754 return TriangleFaceIndexShuffle(()) 755 756 # time step scaling -------------------------------------------------------

757 - def dt_geometric_factor(self, vertices, el):

758 area = abs(2*el.map.jacobian()) 759 semiperimeter = sum(la.norm(vertices[vi1]-vertices[vi2]) 760 for vi1, vi2 in [(0,1), (1,2), (2,0)])/2 761 return area/semiperimeter

762

763 764 765 766 -class TetrahedralElement(SimplicialElement):

767 """An arbitrary-order tetrahedral finite element. 768 769 Coordinate systems used: 770 ======================== 771 772 unit coordinates (r,s,t):: 773 774 ^ s 775 | 776 C 777 /|\\ 778 / | \\ 779 / | \\ 780 / | \\ 781 / O| \\ 782 / __A-----B---> r 783 /_--^ ___--^^ 784 ,D--^^^ 785 t L 786 787 (squint, and it might start making sense...) 788 789 Points in unit coordinates:: 790 791 O=( 0, 0, 0) 792 A=(-1,-1,-1) 793 B=(+1,-1,-1) 794 C=(-1,+1,-1) 795 D=(-1,-1,+1) 796 797 Points in equilateral coordinates (x,y,z):: 798 799 O = (0,0) 800 A = (-1,-1/sqrt(3),-1/sqrt(6)) 801 B = ( 1,-1/sqrt(3),-1/sqrt(6)) 802 C = ( 0, 2/sqrt(3),-1/sqrt(6)) 803 D = ( 0, 0, 3/sqrt(6)) 804 805 When global vertices are passed in, they are mapped to the 806 reference vertices A, B, C, D in order. 807 808 Faces are always ordered ABC, ABD, ACD, BCD. 809 """ 810 811 # In case you were wondering: the double backslashes in the docstring 812 # above are required because single backslashes only escape their subsequent 813 # newlines, and thus end up not yielding a correct docstring. 814 815 dimensions = 3 816 has_local_jacobians = False 817 geometry = hedge.mesh.Tetrahedron 818

819 - def __init__(self, order):

820 self.order = order

821 822 # numbering --------------------------------------------------------------- 823 @memoize_method

824 - def node_tuples(self):

825 """Generate tuples enumerating the node indices present 826 in this element. Each tuple has a length equal to the dimension 827 of the element. The tuples constituents are non-negative integers 828 whose sum is less than or equal to the order of the element. 829 830 The order in which these nodes are generated dictates the local 831 node numbering. 832 """ 833 from pytools import generate_nonnegative_integer_tuples_summing_to_at_most 834 node_tups = list(generate_nonnegative_integer_tuples_summing_to_at_most( 835 self.order, self.dimensions)) 836 837 if False: 838 # hand-tuned node order 839 faces_to_nodes = {} 840 for node_tup in node_tups: 841 faces_to_nodes.setdefault( 842 frozenset(self.faces_for_node_tuple(node_tup)), 843 []).append(node_tup) 844 845 result = [] 846 def add_face_nodes(faces): 847 result.extend(faces_to_nodes.get(frozenset(faces), []))

848 849 add_face_nodes([0,3]) 850 add_face_nodes([0]) 851 add_face_nodes([0,2]) 852 add_face_nodes([0,1]) 853 add_face_nodes([0,1,2]) 854 add_face_nodes([0,1,3]) 855 add_face_nodes([0,2,3]) 856 add_face_nodes([1]) 857 add_face_nodes([1,2]) 858 add_face_nodes([1,2,3]) 859 add_face_nodes([1,3]) 860 add_face_nodes([2]) 861 add_face_nodes([2,3]) 862 add_face_nodes([3]) 863 add_face_nodes([]) 864 865 assert set(result) == set(node_tups) 866 assert len(result) == len(node_tups) 867 868 if True: 869 # average-sort heuristic node order 870 from pytools import average 871 872 def order_number_for_node_tuple(nt): 873 faces = self.faces_for_node_tuple(nt) 874 if not faces: 875 return -1 876 elif len(faces) >= 3: 877 return 1000 878 else: 879 return average(faces)

880 881 def cmp_node_tuples(nt1, nt2): 882 return cmp( 883 order_number_for_node_tuple(nt1), 884 order_number_for_node_tuple(nt2)) 885 886 result = node_tups 887 #result.sort(cmp_node_tuples) 888 889 for i, nt in enumerate(result): 890 fnt = self.faces_for_node_tuple(nt) 891 #print i, nt, fnt 892 893 return result 894

895 - def faces_for_node_tuple(self, node_tuple):

896 """Return the list of face indices of faces on which the node 897 represented by C{node_tuple} lies. 898 """ 899 m,n,o = node_tuple 900 result = [] 901 902 if o == 0: 903 result.append(0) 904 if n == 0: 905 result.append(1) 906 if m == 0: 907 result.append(2) 908 if n+m+o == self.order: 909 result.append(3) 910 911 return result

912 913 # node wrangling ---------------------------------------------------------- 914 @staticmethod

915 - def barycentric_to_equilateral((lambda1, lambda2, lambda3, lambda4)):

916 """Return the equilateral (x,y) coordinate corresponding 917 to the barycentric coordinates (lambda1..lambdaN).""" 918 919 # reflects vertices in equilateral coordinates 920 return numpy.array([ 921 (-lambda1 +lambda2 ), 922 (-lambda1 -lambda2 +2*lambda3 )/sqrt(3.0), 923 (-lambda1 -lambda2 - lambda3 +3*lambda4)/sqrt(6.0), 924 ])

925 926 # see doc/hedge-notes.tm 927 equilateral_to_unit = AffineMap( 928 numpy.array([ 929 [1,-1/sqrt(3),-1/sqrt(6)], 930 [0, 2/sqrt(3),-1/sqrt(6)], 931 [0, 0, sqrt(6)/2] 932 ]), 933 numpy.array([-1/2,-1/2,-1/2])) 934

935 - def equilateral_nodes(self):

936 """Generate warped nodes in equilateral coordinates (x,y).""" 937 938 # port of Hesthaven/Warburton's Nodes3D routine 939 940 # Set optimized parameter alpha, depending on order N 941 alpha_opt = [0,0,0,0.1002, 1.1332,1.5608,1.3413,1.2577,1.1603, 942 1.10153,0.6080,0.4523,0.8856,0.8717,0.9655] 943 944 try: 945 alpha = alpha_opt[self.order-1] 946 except IndexError: 947 alpha = 1 948 949 from pytools import wandering_element 950 951 vertices = [self.barycentric_to_equilateral(bary) 952 for bary in wandering_element(self.dimensions+1)] 953 all_vertex_indices = range(self.dimensions+1) 954 face_vertex_indices = self.geometry \ 955 .face_vertices(all_vertex_indices) 956 faces_vertices = self.geometry \ 957 .face_vertices(vertices) 958 959 bary_points = list(self.equidistant_barycentric_nodes()) 960 equi_points = [self.barycentric_to_equilateral(bp) 961 for bp in bary_points] 962 963 from tools import normalize 964 from operator import add, mul 965 966 tri_warp = TriangleWarper(alpha, self.order) 967 968 for fvi, (v1, v2, v3) in zip(face_vertex_indices, faces_vertices): 969 # find directions spanning the face: "base" and "altitude" 970 directions = [normalize(v2-v1), normalize((v3)-(v1+v2)/2)] 971 972 # the two should be orthogonal 973 assert abs(numpy.dot(directions[0],directions[1])) < 1e-16 974 975 # find the vertex opposite to the current face 976 opp_vertex_index = (set(all_vertex_indices) - set(fvi)).__iter__().next() 977 978 shifted = [] 979 for bp, ep in zip(bary_points, equi_points): 980 face_bp = [bp[i] for i in fvi] 981 982 blend = reduce(mul, face_bp) * (1+alpha*bp[opp_vertex_index])**2 983 984 for i in fvi: 985 denom = bp[i] + 0.5*bp[opp_vertex_index] 986 if abs(denom) > 1e-12: 987 blend /= denom 988 else: 989 blend = 0.5 # each edge gets shifted twice 990 break 991 992 shifted.append(ep + blend*reduce(add, 993 (tw*dir for tw, dir in zip(tri_warp(face_bp), directions)))) 994 995 equi_points = shifted 996 997 return equi_points

998 999 @memoize_method

1000 - def get_submesh_indices(self):

1001 """Return a list of tuples of indices into the node list that 1002 generate a tesselation of the reference element.""" 1003 1004 node_dict = dict( 1005 (ituple, idx) 1006 for idx, ituple in enumerate(self.node_tuples())) 1007 1008 def add_tuples(a, b): 1009 return tuple(ac+bc for ac, bc in zip(a,b))

1010 1011 def try_add_tet(d1, d2, d3, d4): 1012 try: 1013 result.append(( 1014 node_dict[add_tuples(current, d1)], 1015 node_dict[add_tuples(current, d2)], 1016 node_dict[add_tuples(current, d3)], 1017 node_dict[add_tuples(current, d4)], 1018 )) 1019 except KeyError, e: 1020 pass 1021 1022 result = [] 1023 for current in self.node_tuples(): 1024 # this is a tesselation of a cube into six tets. 1025 # subtets that fall outside of the master tet are simply not added. 1026 1027 # positively oriented 1028 try_add_tet((0,0,0), (1,0,0), (0,1,0), (0,0,1)) 1029 try_add_tet((1,0,1), (1,0,0), (0,0,1), (0,1,0)) 1030 try_add_tet((1,0,1), (0,1,1), (0,1,0), (0,0,1)) 1031 1032 try_add_tet((1,0,0), (0,1,0), (1,0,1), (1,1,0)) 1033 try_add_tet((0,1,1), (0,1,0), (1,1,0), (1,0,1)) 1034 try_add_tet((0,1,1), (1,1,1), (1,0,1), (1,1,0)) 1035 1036 # negatively oriented 1037 #try_add_tet((0,0,0), (1,0,0), (0,0,1), (0,1,0)) 1038 #try_add_tet((1,0,1), (1,0,0), (0,1,0), (0,0,1)) 1039 #try_add_tet((1,0,1), (0,1,1), (0,0,1), (0,1,0)) 1040 1041 #try_add_tet((1,0,0), (0,1,0), (1,1,0), (1,0,1)) 1042 #try_add_tet((0,1,1), (0,1,0), (1,0,1), (1,1,0)) 1043 #try_add_tet((0,1,1), (1,1,1), (1,1,0), (1,0,1)) 1044 return result 1045 1046 # basis functions --------------------------------------------------------- 1047 @memoize_method

1048 - def basis_functions(self):

1049 """Get a sequence of functions that form a basis of the approximation space. 1050 1051 The approximation space is spanned by the polynomials:: 1052 1053 r**i * s**j * t**k for i+j+k <= order 1054 """ 1055 return [TetrahedronBasisFunction(*idx) for idx in 1056 self.generate_mode_identifiers()]

1057

1058 - def grad_basis_functions(self):

1059 """Get the (r,s,...) gradient functions of the basis_functions(), 1060 in the same order. 1061 """ 1062 return [GradTetrahedronBasisFunction(*idx) for idx in 1063 self.generate_mode_identifiers()]

1064 1065 # face operations --------------------------------------------------------- 1066 @memoize_method

1067 - def face_mass_matrix(self):

1068 from hedge.polynomial import generic_vandermonde 1069 1070 from pytools import generate_nonnegative_integer_tuples_summing_to_at_most 1071 1072 basis = [TriangleBasisFunction(*node_tup) 1073 for node_tup in 1074 generate_nonnegative_integer_tuples_summing_to_at_most( 1075 self.order, self.dimensions-1)] 1076 1077 unodes = self.unit_nodes() 1078 face_indices = self.face_indices() 1079 1080 face_vandermonde = generic_vandermonde( 1081 [unodes[i][:2] for i in face_indices[0]], 1082 basis) 1083 1084 return numpy.asarray( 1085 la.inv( 1086 numpy.dot(face_vandermonde, face_vandermonde.T)), 1087 order="C")

1088

1089 - def get_face_index_shuffle_to_match(self, face_1_vertices, face_2_vertices):

1090 (a,b,c) = face_1_vertices 1091 f2_tuple = face_2_vertices 1092 1093 try: 1094 idx_map = self._shuffle_face_idx_map 1095 except AttributeError: 1096 idx_map = self._shuffle_face_idx_map = {} 1097 idx = 0 1098 for j in range(0, self.order+1): 1099 for i in range(0, self.order+1-j): 1100 self._shuffle_face_idx_map[i,j] = idx 1101 idx += 1 1102 1103 order = self.order 1104 1105 class TetrahedronFaceIndexShuffle: 1106 def __init__(self, operations): 1107 self.operations = operations

1108 1109 def __hash__(self): 1110 return hash(self.operations) 1111 1112 def __eq__(self, other): 1113 return self.operations == other.operations 1114 1115 def __call__(self, indices): 1116 for op in self.operations: 1117 if op == "flip": 1118 indices = self.flip(indices) 1119 elif op == "shift_left": 1120 indices = self.shift_left(indices) 1121 else: 1122 raise RuntimeError, "invalid operation" 1123 return indices 1124 1125 # flip and shift_left generate S_3 1126 def flip(self, indices): 1127 """Flip the indices along the unit hypotenuse.""" 1128 result = [] 1129 for j in range(0, order+1): 1130 for i in range(0, order+1-j): 1131 result.append(indices[idx_map[j,i]]) 1132 return tuple(result) 1133 1134 def shift_left(self, indices): 1135 """Rotate all edges to the left.""" 1136 result = len(indices)*[0] 1137 idx = 0 1138 for j in range(0, order+1): 1139 for i in range(0, order+1-j): 1140 result[idx_map[j, order-i-j]] = indices[idx] 1141 idx += 1 1142 return tuple(result) 1143 1144 # yay, enumerate S_3 by hand 1145 if f2_tuple == (a,b,c): 1146 #return face_2_indices 1147 return TetrahedronFaceIndexShuffle(()) 1148 elif f2_tuple == (a,c,b): 1149 #return flip(face_2_indices) 1150 return TetrahedronFaceIndexShuffle(("flip",)) 1151 elif f2_tuple == (b,c,a): 1152 # (b,c,a) -sl-> (c,a,b) -sl-> (a,b,c) 1153 #return shift_left(shift_left(face_2_indices)) 1154 return TetrahedronFaceIndexShuffle(("shift_left", "shift_left")) 1155 elif f2_tuple == (b,a,c): 1156 # (b,a,c) -sl-> (a,c,b) -fl-> (a,b,c) 1157 #return flip(shift_left(face_2_indices)) 1158 return TetrahedronFaceIndexShuffle(("shift_left", "flip")) 1159 elif f2_tuple == (c,a,b): 1160 # (c,a,b) -sl-> (a,b,c) 1161 #return shift_left(face_2_indices) 1162 return TetrahedronFaceIndexShuffle(("shift_left",)) 1163 elif f2_tuple == (c,b,a): 1164 # (c,b,a) -fl-> (c,a,b) -sl-> (a,b,c) 1165 #return shift_left(flip(face_2_indices)) 1166 return TetrahedronFaceIndexShuffle(("flip", "shift_left")) 1167 else: 1168 raise FaceVertexMismatch("face vertices do not match") 1169 1170 # time step scaling -------------------------------------------------------

1171 - def dt_geometric_factor(self, vertices, el):

1172 result = abs(el.map.jacobian())/max(abs(fj) for fj in el.face_jacobians) 1173 if self.order in [1, 2]: 1174 from warnings import warn 1175 warn("cowardly halving timestep for order 1 and 2 tets to avoid CFL issues") 1176 result /= 2 1177 1178 return result

1179 1180 1181 1182 1183 ELEMENTS = [IntervalElement, TriangularElement, TetrahedralElement] 1184

Source Code for Module hedge.element