1 """This module contains prewritten code for the computation
2 of the geometry mapping to the reference cell."""
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
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
100
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
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)
121 return code + "\n"
122
123
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