1
2 import swiginac
3 import SyFi
4
5 from sfc.symbolic_utils import symbolic_matrix, symbolic_vector, symbol, symbols, cross, sub, inner
6 from sfc.common.utilities import matrix_to_list, list_to_vector, index_string, dot_product
7 from sfc.geometry import SFCCell, geometry_mapping
8
9
10
11
13 if cell == "interval":
14 return SyFi.ReferenceLine()
15 elif cell == "triangle":
16 return SyFi.ReferenceTriangle()
17 elif cell == "tetrahedron":
18 return SyFi.ReferenceTetrahedron()
19 elif cell == "quadrilateral":
20 return SyFi.ReferenceRectangle()
21 elif cell == "hexahedron":
22 return SyFi.ReferenceBox()
23 raise ValueError("Unknown element cell '%s'." % cell)
24
25
27 if shape == "interval":
28 polygon = SyFi.Line(*p)
29 elif shape == "triangle":
30 polygon = SyFi.Triangle(*p)
31 elif shape == "tetrahedron":
32 polygon = SyFi.Tetrahedron(*p)
33 elif shape == "quadrilateral":
34
35
36 polygon = SyFi.Rectangle(*p)
37 elif shape == "hexahedron":
38
39
40 polygon = SyFi.Box(*p)
41 return polygon
42
43
48
50 s = ""
51 s += " Geometry tokens:\n"
52 s += " self.sfc_cell = %s\n" % self.sfc_cell
53
54 s += " self.vx_sym = %s\n" % self.vx_sym
55 s += " self.vx_expr = %s\n" % self.vx_expr
56 s += " self.xi_sym = %s\n" % self.xi_sym
57 s += " self.x_sym = %s\n" % self.x_sym
58 s += " self.x_expr = %s\n" % self.x_expr
59 s += " self.G_sym = %s\n" % self.G_sym
60 s += " self.G_expr, = %s\n" % self.G_expr,
61 s += " self.detGtmp_sym = %s\n" % self.detGtmp_sym
62 s += " self.detGtmp_expr = %s\n" % self.detGtmp_expr
63 s += " self.detG_sym = %s\n" % self.detG_sym
64 s += " self.detG_expr = %s\n" % self.detG_expr
65 s += " self.Ginv_sym = %s\n" % self.Ginv_sym
66 s += " self.Ginv_expr = %s\n" % self.Ginv_expr
67 s += " self.quad_weight_sym = %s\n" % self.quad_weight_sym
68 return s
69
71 "Build tokens for variables derived from cell."
72
73
74 self.syfi_polygon = create_syfi_polygon(self.polygon_name)
75 self.sfc_cell = SFCCell(self.syfi_polygon)
76
77 self.cell = self.sfc_cell
78
79 sfc_cell = self.sfc_cell
80 nsd = sfc_cell.nsd
81
82
83 self.vx_sym = []
84 self.vx_expr = []
85 for i in range(sfc_cell.num_vertices):
86 s = swiginac.matrix(nsd, 1, symbols("vx%d_%d" % (i,j) for j in range(nsd)))
87 e = swiginac.matrix(nsd, 1, symbols("c.coordinates[%d][%d]" % (i,j)
88 for j in range(nsd)))
89 self.vx_sym.append(s)
90 self.vx_expr.append(e)
91
92
93
94
95
96 p = [matrix_to_list(vx) for vx in self.vx_sym]
97 self.global_polygon = build_global_polygon(sfc_cell.shape, p)
98
99
100
101 self.G_sym = symbolic_matrix(nsd, nsd, "G")
102
103 self.x0_sym = self.vx_sym[0]
104
105 self.G_expr, self.x0_expr = geometry_mapping(self.global_polygon)
106
107
108 self.detGtmp_sym = symbol("detGtmp")
109 self.detGtmp_expr = swiginac.determinant(self.G_sym)
110 self.detG_sym = symbol("detG")
111
112 self.detG_expr = swiginac.abs(self.detGtmp_sym)
113
114
115 self.detG_sign_sym = symbol("detG_sign")
116 self.detG_sign_expr = self.detG_sym / swiginac.abs(self.detG_sym)
117
118
119 self.Ginv_sym = symbolic_matrix(nsd, nsd, "Ginv")
120 self.Ginv_expr = ((self.detGtmp_expr * swiginac.inverse(self.G_sym))
121 / self.detGtmp_sym)
122
123
124
125
126 self.xi_sym = swiginac.matrix(nsd, 1, list(sfc_cell.xi[:nsd]))
127
128 x_sym = symbols(("x%d" % i for i in range(nsd)))
129 x_expr = [self.x0_sym[i] + sum(self.G_sym[i,j]*self.xi_sym[j] for j in range(nsd))
130 for i in range(nsd)]
131 self.x_sym = swiginac.matrix(nsd, 1, x_sym)
132 self.x_expr = swiginac.matrix(nsd, 1, x_expr)
133
134
135 self.quad_weight_sym = symbol("quad_weight")
136 self.D_sym = symbol("D")
137 self.D_expr = self.quad_weight_sym * self.detG_sym
138
139
140
141
142
143
144 self.n_tmp_sym = symbolic_vector(nsd, "ntmp")
145 self.n_tmp_expr = [[sum(self.G_sym[j,i]*n[j] for j in range(nsd))
146 for i in range(nsd)]
147 for n in sfc_cell.facet_n]
148
149 self.nnorm_sym = symbol('nnorm')
150 self.nnorm_expr = swiginac.sqrt(sum(self.n_tmp_sym[i]**2 for i in range(nsd)))
151
152 self.n_sym = symbolic_vector(nsd, "n")
153 self.n_expr = [(self.n_tmp_sym[i]/self.nnorm_sym).evalm() for i in range(nsd)]
154
155 self.facet_D_sym = symbol("facet_D")
156 self.facet_D_expr = build_facet_D(sfc_cell, self.G_sym)
157
171
173 one = swiginac.numeric(1)
174 return [one]*len(facet_polygons)
175
177 facet_D_expr = []
178 for facet, facet_polygon in enumerate(facet_polygons):
179
180 v0 = G_sym * list_to_vector(facet_polygon.vertex(0))
181 v1 = G_sym * list_to_vector(facet_polygon.vertex(1))
182 v = sub(v1, v0)
183 D = swiginac.sqrt( inner(v, v) )
184 if facet == 0:
185 D = D / swiginac.sqrt(2)
186
187
188 facet_D_expr.append(D)
189 return facet_D_expr
190
192 facet_D_expr = []
193 for facet, facet_polygon in enumerate(facet_polygons):
194
195 v0 = facet_polygon.vertex(0)
196 v1 = facet_polygon.vertex(1)
197 v2 = facet_polygon.vertex(2)
198
199 v0 = G_sym * list_to_vector(v0)
200 v1 = G_sym * list_to_vector(v1)
201 v2 = G_sym * list_to_vector(v2)
202
203
204
205
206
207
208
209 c = cross(sub(v1, v0), sub(v2, v0))
210 D = swiginac.sqrt( inner(c, c) )
211 if facet == 0:
212 D = D / swiginac.sqrt(3)
213
214
215 facet_D_expr.append(D)
216 return facet_D_expr
217
219 facet_D_expr = []
220 for facet, facet_polygon in enumerate(facet_polygons):
221
222 v0 = G_sym * list_to_vector(facet_polygon.vertex(0))
223 v1 = G_sym * list_to_vector(facet_polygon.vertex(1))
224 v = sub(v1, v0)
225 D = swiginac.sqrt(inner(v,v))
226
227
228 facet_D_expr.append(D)
229 return facet_D_expr
230
232 facet_D_expr = []
233 for facet, facet_polygon in enumerate(facet_polygons):
234
235 v0 = G_sym * list_to_vector(facet_polygon.vertex(0))
236 v1 = G_sym * list_to_vector(facet_polygon.vertex(1))
237 v2 = G_sym * list_to_vector(facet_polygon.vertex(2))
238 v3 = G_sym * list_to_vector(facet_polygon.vertex(3))
239
240 m0 = (v0 + v1)/2
241 m1 = (v1 + v2)/2
242 m2 = (v2 + v3)/2
243 m3 = (v3 + v0)/2
244
245 c = cross(sub(m2, m0), sub(m1, m3))
246 D = swiginac.sqrt( inner(c, c) )
247
248
249 facet_D_expr.append(D)
250 return facet_D_expr
251