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

Source Code for Module sfc.geometry.affine_map

  1  """This module contains prewritten code for the computation 
  2  of the geometry mapping to the reference cell.""" 
  3   
  4  # Copyright (C) 2007 Martin Sandve Alnes and Simula Resarch Laboratory 
  5  # 
  6  # This file is part of SyFi. 
  7  # 
  8  # SyFi 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 2 of the License, or 
 11  # (at your option) any later version. 
 12  # 
 13  # SyFi 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 SyFi. If not, see <http://www.gnu.org/licenses/>. 
 20  # 
 21  # First added:  2007-05-28 
 22  # Last changed: 2007-05-28 
 23   
 24  geometry_string1d = """\ 
 25  // coordinates 
 26  double x0 = c.coordinates[0][0]; 
 27  double x1 = c.coordinates[1][0]; 
 28   
 29  // affine map 
 30  double G00 = x1 - x0; 
 31  """ 
 32   
 33  detG_string1d = """\ 
 34  double detG_tmp = G00; 
 35  double detG = fabs(detG_tmp); 
 36  """ 
 37   
 38  GinvT_string1d = """\ 
 39  double GinvT00 =  1.0 / G00; 
 40  """ 
 41   
 42   
 43   
 44  geometry_string2d = """\ 
 45  // coordinates 
 46  double x0 = c.coordinates[0][0]; double y0 = c.coordinates[0][1]; 
 47  double x1 = c.coordinates[1][0]; double y1 = c.coordinates[1][1]; 
 48  double x2 = c.coordinates[2][0]; double y2 = c.coordinates[2][1]; 
 49   
 50  // affine map 
 51  double G00 = x1 - x0; 
 52  double G01 = x2 - x0; 
 53   
 54  double G10 = y1 - y0; 
 55  double G11 = y2 - y0; 
 56  """ 
 57   
 58  detG_string2d = """\ 
 59  double detG_tmp = G00*G11-G01*G10; 
 60  double detG = fabs(detG_tmp); 
 61  """ 
 62   
 63  GinvT_string2d = """\ 
 64  double GinvT00 =  G11 / detG_tmp; 
 65  double GinvT01 = -G10 / detG_tmp; 
 66  double GinvT10 = -G01 / detG_tmp; 
 67  double GinvT11 =  G00 / detG_tmp; 
 68  """ 
 69   
 70   
 71  geometry_string3d = """\ 
 72  // coordinates 
 73  double x0 = c.coordinates[0][0]; double y0 = c.coordinates[0][1]; double z0 = c.coordinates[0][2]; 
 74  double x1 = c.coordinates[1][0]; double y1 = c.coordinates[1][1]; double z1 = c.coordinates[1][2]; 
 75  double x2 = c.coordinates[2][0]; double y2 = c.coordinates[2][1]; double z2 = c.coordinates[2][2]; 
 76  double x3 = c.coordinates[3][0]; double y3 = c.coordinates[3][1]; double z3 = c.coordinates[3][2]; 
 77   
 78  // affine map 
 79  double G00 = x1 - x0; 
 80  double G10 = y1 - y0; 
 81  double G20 = z1 - z0; 
 82   
 83  double G01 = x2 - x0; 
 84  double G11 = y2 - y0; 
 85  double G21 = z2 - z0; 
 86   
 87  double G02 = x3 - x0; 
 88  double G12 = y3 - y0; 
 89  double G22 = z3 - z0; 
 90  """ 
 91   
 92  detG_string3d = """\ 
 93  double detG_tmp =    G00*(G11*G22-G21*G12) 
 94                     - G01*(G10*G22-G20*G12) 
 95                     + G02*(G10*G21-G20*G11); 
 96  double detG = fabs(detG_tmp); 
 97  """ 
 98   
 99  # TODO: this can be optimized further, f.ex. G11*G22 occurs both below and in detG 
100  # TODO: write this in symbolic form 
101  GinvT_string3d = """\ 
102  double GinvT00 = ( G11*G22-G12*G21) / detG_tmp; 
103  double GinvT01 = (-G22*G10+G12*G20) / detG_tmp; 
104  double GinvT02 = (-G11*G20+G10*G21) / detG_tmp; 
105  double GinvT10 = ( G02*G21-G22*G01) / detG_tmp; 
106  double GinvT11 = ( G22*G00-G02*G20) / detG_tmp; 
107  double GinvT12 = (-G21*G00+G20*G01) / detG_tmp; 
108  double GinvT20 = ( G12*G01-G11*G02) / detG_tmp; 
109  double GinvT21 = ( G10*G02-G12*G00) / detG_tmp; 
110  double GinvT22 = ( G11*G00-G10*G01) / detG_tmp; 
111  """ 
112   
113 -def gen_geometry_code(nsd, detG=True, GinvT=False):
114 if nsd != 1 and nsd != 2 and nsd != 3: 115 raise RuntimeError("gen_geometry_code not implemented for nsd != 1, 2 or 3") 116 code = eval("geometry_string%dd" % nsd) 117 if GinvT or detG: 118 code += "\n" + eval("detG_string%dd" % nsd) 119 if GinvT: 120 code += "\n" + eval("GinvT_string%dd" % nsd) # TODO: use GinvT_string(nsd) instead after verification 121 return code + "\n"
122 123
124 -def GinvT_string(nsd): # TODO: verify and use this above
125 from sfc.CodeFormatter import gen_token_assignments 126 from sfc.symbol_factory import symbolic_matrix 127 G_sym = symbolic_matrix(nsd, nsd, "G") 128 GinvT_sym = symbolic_matrix(nsd, nsd, "GinvT") 129 GinvT = G_sym.inverse().transpose() 130 code = gen_token_assignments( [( GinvT_sym[i], GinvT[i] ) for i in range(nsd**2) ] ) 131 return code 132 133 134 if __name__ == '__main__': 135 print GinvT_string(2) 136 print "\n\n" 137 print gen_geometry_code(2, True, True) 138