Package sfc :: Package geometry :: Module mappings
[hide private]
[frames] | no frames]

Source Code for Module sfc.geometry.mappings

  1      # -*- coding: utf-8 -*- 
  2  """ 
  3  This module contains functions for computing 
  4  integrals both symbolically and numerically. 
  5  """ 
  6   
  7  # Copyright (C) 2007-2009 Martin Sandve Alnes and Simula Resarch Laboratory 
  8  # 
  9  # This file is part of SyFi. 
 10  # 
 11  # SyFi is free software: you can redistribute it and/or modify 
 12  # it under the terms of the GNU General Public License as published by 
 13  # the Free Software Foundation, either version 2 of the License, or 
 14  # (at your option) any later version. 
 15  # 
 16  # SyFi is distributed in the hope that it will be useful, 
 17  # but WITHOUT ANY WARRANTY; without even the implied warranty of 
 18  # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 
 19  # GNU General Public License for more details. 
 20  # 
 21  # You should have received a copy of the GNU General Public License 
 22  # along with SyFi. If not, see <http://www.gnu.org/licenses/>. 
 23  # 
 24  # First added:  2007-11-07 
 25  # Last changed: 2009-03-04 
 26   
 27  import swiginac 
 28  import SyFi 
 29   
 30  from sfc.quadrature import find_quadrature_rule 
 31  from sfc.symbolic_utils import symbol, symbols 
 32   
33 -def geometry_mapping(polygon):
34 """ 35 This function computes the geometry mapping from 36 the corresponding reference polygon to polygon. 37 It returns the tuple (G,x0) in the 38 mapping x = G xi + x0 39 """ 40 41 vlen = SyFi.cvar.nsd 42 43 # If we get a point, return this trivial mapping 44 if isinstance(polygon, swiginac.matrix): 45 # x = 1 * xi + x0 = x0 because xi is always zero at the point 46 assert vlen == 1 47 assert polygon.rows() == 1 48 assert polygon.cols() == 1 49 G = swiginac.matrix(1, 1, [1]) 50 x0 = polygon 51 return G, x0 52 53 Gl = [] 54 v0 = polygon.vertex(0) 55 v1 = polygon.vertex(1) 56 57 if not vlen == len(v0): 58 print 59 print vlen 60 print v0 61 print v1 62 print polygon 63 print 64 65 assert vlen == len(v0) 66 67 if isinstance(polygon, SyFi.Line): 68 Gl.extend(v1[k] - v0[k] for k in range(vlen)) 69 70 elif isinstance(polygon, SyFi.Triangle): 71 if vlen == 2 or vlen == 3: 72 v2 = polygon.vertex(2) 73 Gl.extend(v1[k] - v0[k] for k in range(vlen)) 74 Gl.extend(v2[k] - v0[k] for k in range(vlen)) 75 else: 76 raise RuntimeError("Unsupported vertex size.") 77 78 elif isinstance(polygon, SyFi.Rectangle): 79 if vlen == 2 or vlen == 3: # TODO: improve mapping 80 v2 = polygon.vertex(2) 81 v3 = polygon.vertex(3) 82 Gl.extend(v1[k] - v0[k] for k in range(vlen)) 83 Gl.extend(v3[k] - v0[k] for k in range(vlen)) 84 else: 85 raise RuntimeError("Unsupported vertex size.") 86 87 elif isinstance(polygon, SyFi.Tetrahedron): 88 if vlen != 3: 89 raise RuntimeError("Unsupported vertex size.") 90 v2 = polygon.vertex(2) 91 v3 = polygon.vertex(3) 92 Gl.extend(v1[k] - v0[k] for k in range(vlen)) 93 Gl.extend(v2[k] - v0[k] for k in range(vlen)) 94 Gl.extend(v3[k] - v0[k] for k in range(vlen)) 95 96 elif isinstance(polygon, SyFi.Box): 97 if vlen != 3: 98 raise RuntimeError("Unsupported vertex size.") 99 v2 = polygon.vertex(2) # TODO: improve mapping 100 v3 = polygon.vertex(3) 101 v4 = polygon.vertex(4) 102 Gl.extend(v1[k] - v0[k] for k in range(vlen)) 103 Gl.extend(v3[k] - v0[k] for k in range(vlen)) 104 Gl.extend(v4[k] - v0[k] for k in range(vlen)) 105 106 else: 107 raise RuntimeError("Unsupported polygon type.") 108 109 G = swiginac.matrix(vlen, vlen, Gl).transpose() 110 x0 = swiginac.matrix(vlen, 1, v0) 111 112 #print "G =", G 113 114 return G, x0
115