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  from sfc.geometry.UFCCell import UFCCell 
 33   
 34   
35 -def geometry_mapping(polygon):
36 """ 37 This function computes the geometry mapping from 38 the corresponding reference polygon to polygon. 39 It returns the tuple (G,x0) in the 40 mapping x = G xi + x0 41 """ 42 43 vlen = SyFi.cvar.nsd 44 45 # If we get a point, return this trivial mapping 46 if isinstance(polygon, swiginac.matrix): 47 # x = 1 * xi + x0 = x0 because xi is always zero at the point 48 assert vlen == 1 49 assert polygon.rows() == 1 50 assert polygon.cols() == 1 51 G = swiginac.matrix(1, 1, [1]) 52 x0 = polygon 53 return G, x0 54 55 Gl = [] 56 v0 = polygon.vertex(0) 57 v1 = polygon.vertex(1) 58 59 assert vlen == len(v0) 60 61 if isinstance(polygon, SyFi.Line): 62 Gl.extend(v1[k] - v0[k] for k in range(vlen)) 63 64 elif isinstance(polygon, SyFi.Triangle): 65 if vlen == 2 or vlen == 3: 66 v2 = polygon.vertex(2) 67 Gl.extend(v1[k] - v0[k] for k in range(vlen)) 68 Gl.extend(v2[k] - v0[k] for k in range(vlen)) 69 else: 70 raise RuntimeError("Unsupported vertex size.") 71 72 elif isinstance(polygon, SyFi.Rectangle): 73 if vlen == 2 or vlen == 3: # TODO: improve mapping 74 v2 = polygon.vertex(2) 75 v3 = polygon.vertex(3) 76 Gl.extend(v1[k] - v0[k] for k in range(vlen)) 77 Gl.extend(v3[k] - v0[k] for k in range(vlen)) 78 else: 79 raise RuntimeError("Unsupported vertex size.") 80 81 elif isinstance(polygon, SyFi.Tetrahedron): 82 if vlen != 3: 83 raise RuntimeError("Unsupported vertex size.") 84 v2 = polygon.vertex(2) 85 v3 = polygon.vertex(3) 86 Gl.extend(v1[k] - v0[k] for k in range(vlen)) 87 Gl.extend(v2[k] - v0[k] for k in range(vlen)) 88 Gl.extend(v3[k] - v0[k] for k in range(vlen)) 89 90 elif isinstance(polygon, SyFi.Box): 91 if vlen != 3: 92 raise RuntimeError("Unsupported vertex size.") 93 v2 = polygon.vertex(2) # TODO: improve mapping 94 v3 = polygon.vertex(3) 95 v4 = polygon.vertex(4) 96 Gl.extend(v1[k] - v0[k] for k in range(vlen)) 97 Gl.extend(v3[k] - v0[k] for k in range(vlen)) 98 Gl.extend(v4[k] - v0[k] for k in range(vlen)) 99 100 else: 101 raise RuntimeError("Unsupported polygon type.") 102 103 G = swiginac.matrix(vlen, vlen, Gl).transpose() 104 x0 = swiginac.matrix(vlen, 1, v0) 105 106 #print "G =", G 107 108 return G, x0
109