1
2
3 """
4 This module contains a class UFCCell to represent the properties of a cell in a easily accessible way.
5 """
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27 import swiginac
28 from swiginac import matrix, exp, sqrt
29
30 import SyFi
31
32 from sfc.symbolic_utils import cross, inner, symbols
33 from sfc.common.output import sfc_debug
34
35
36
39 sfc_debug("Entering UFCCell.__init__")
40
41 assert isinstance(polygon, SyFi.Polygon)
42 self.polygon = polygon
43 name = polygon.str()
44 x, y, z = symbols("xyz")
45
46 if isinstance(polygon, SyFi.ReferenceLine) or name == "ReferenceLine":
47 self.shape = "interval"
48 self.facet_shape = "point"
49
50
51 self.nsd = 1
52 self.num_vertices = 2
53 self.num_edges = 1
54 self.num_faces = 0
55 self.num_facets = self.num_vertices
56 self.num_entities = (self.num_vertices, self.num_edges)
57
58
59 self.facet_vertices = [
60 (0,),
61 (1,)
62 ]
63
64
65 self.vertices = [matrix(self.nsd, 1, polygon.vertex(i)) for i in range(self.num_vertices)]
66 self.facet_polygons = self.vertices
67
68
69 n = matrix(1, 1, [1])
70 self.facet_n = [-n, +n]
71
72
73 self.facet_equations = [x, x-1]
74
75 elif isinstance(polygon, SyFi.ReferenceTriangle) or name == "ReferenceTriangle":
76 self.shape = "triangle"
77 self.facet_shape = "interval"
78
79
80 self.nsd = 2
81 self.num_vertices = 3
82 self.num_edges = 3
83 self.num_faces = 1
84 self.num_facets = self.num_edges
85 self.num_entities = (self.num_vertices, self.num_edges, self.num_faces)
86
87
88
89 self.facet_vertices = [
90 (1, 2),
91 (2, 0),
92 (0, 1)
93 ]
94
95 self.vertices = [matrix(self.nsd, 1, polygon.vertex(i)) for i in range(self.num_vertices)]
96 self.facet_polygons = [polygon.line(i) for i in range(self.num_facets)]
97
98
99 self.facet_n = []
100 for facet in range(self.num_facets):
101 fromvert, tovert = self.facet_vertices[facet]
102
103
104 t = self.vertices[fromvert] - self.vertices[tovert]
105 n = matrix(2, 1, [t[1], -t[0]])
106 n = n / sqrt(inner(n,n))
107
108 self.facet_n.append(n)
109
110
111 self.edge_equations = [x+y-1, x, y]
112 self.facet_equations = self.edge_equations
113
114 elif isinstance(polygon, SyFi.ReferenceTetrahedron) or name == "ReferenceTetrahedron":
115 self.shape = "tetrahedron"
116 self.facet_shape = "triangle"
117
118
119 self.nsd = 3
120 self.num_vertices = 4
121 self.num_edges = 6
122 self.num_faces = 4
123 self.num_facets = self.num_faces
124 self.num_entities = (self.num_vertices, self.num_edges, self.num_faces, 1)
125
126
127
128 self.facet_vertices = [
129 (1, 2, 3),
130 (0, 3, 2),
131 (0, 1, 3),
132 (0, 2, 1)
133 ]
134
135
136 self.vertices = [matrix(self.nsd, 1, polygon.vertex(i)) for i in range(self.num_vertices)]
137 self.facet_polygons = [polygon.triangle(i) for i in range(self.num_facets)]
138
139
140 self.facet_n = []
141 for facet in range(self.num_facets):
142 vert = self.facet_vertices[facet]
143
144 t01 = self.vertices[vert[1]] - self.vertices[vert[0]]
145 t02 = self.vertices[vert[2]] - self.vertices[vert[0]]
146 n = cross(t01, t02)
147 n = n / sqrt(inner(n,n))
148
149 self.facet_n.append(n)
150
151
152 self.face_equations = [ (x+y+z-1), x, y, z]
153 self.facet_equations = self.face_equations
154
155
156 p = matrix(3, 1, [x, y, z])
157 v0, v1, v2, v3 = self.vertices
158 self.edge_equations = [ cross( (p - v3), (v2 - v3) ),
159 cross( (p - v3), (v1 - v3) ),
160 cross( (p - v2), (v1 - v2) ),
161 cross( (p - v3), (v0 - v3) ),
162 cross( (p - v2), (v0 - v2) ),
163 cross( (p - v1), (v0 - v1) ) ]
164
165 elif isinstance(polygon, SyFi.ReferenceRectangle) or name == "ReferenceRectangle":
166 self.shape = "quadrilateral"
167 self.facet_shape = "interval"
168
169
170 self.nsd = 2
171 self.num_vertices = 4
172 self.num_edges = 4
173 self.num_faces = 1
174 self.num_facets = self.num_edges
175 self.num_entities = (self.num_vertices, self.num_edges, self.num_faces)
176
177
178
179 self.facet_vertices = [
180 (2, 3),
181 (1, 2),
182 (3, 0),
183 (0, 1)
184 ]
185
186
187 self.vertices = [matrix(self.nsd, 1, polygon.vertex(i)) for i in range(self.num_vertices)]
188 self.facet_polygons = [polygon.line(i) for i in range(self.num_facets)]
189
190
191 self.facet_n = []
192 for facet in range(self.num_facets):
193 fromvert, tovert = self.facet_vertices[facet]
194
195 t = self.vertices[fromvert] - self.vertices[tovert]
196 n = matrix(2, 1, [t[1], -t[0]])
197
198 self.facet_n.append(n)
199
200
201 self.edge_equations = [x-1, y-1, x, y]
202 self.facet_equations = self.edge_equations
203
204 elif isinstance(polygon, SyFi.ReferenceBox) or name == "ReferenceBox":
205 self.shape = "hexahedron"
206 self.facet_shape = "quadrilateral"
207
208
209 self.nsd = 3
210 self.num_vertices = 8
211 self.num_edges = 12
212 self.num_faces = 6
213 self.num_facets = self.num_faces
214 self.num_entities = (self.num_vertices, self.num_edges, self.num_faces, 1)
215
216
217
218 self.facet_vertices = [
219 (4, 5, 6, 7),
220 (2, 3, 7, 6),
221 (1, 2, 6, 5),
222 (0, 4, 7, 3),
223 (0, 1, 5, 4),
224 (0, 3, 2, 1)
225 ]
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247 self.vertices = [matrix(self.nsd, 1, polygon.vertex(i)) for i in range(self.num_vertices)]
248 self.facet_polygons = [polygon.rectangle(i) for i in range(self.num_facets)]
249
250 self.facet_n = []
251 for facet in range(self.num_facets):
252
253 vert = self.facet_vertices[facet]
254
255
256 t02 = self.vertices[vert[0]] - self.vertices[vert[2]]
257 t13 = self.vertices[vert[1]] - self.vertices[vert[3]]
258 n = cross(t02, t13)
259
260
261 self.facet_n.append(n)
262
263
264
265 self.face_equations = [z-1, y-1, x-1, x, y, z]
266 self.facet_equations = self.face_equations
267
268
269 f = self.face_equations
270 e = [0]*12
271 e[0 ] = (f[0], f[1])
272 e[1 ] = (f[0], f[2])
273 e[2 ] = (f[0], f[3])
274 e[3 ] = (f[0], f[4])
275 e[4 ] = (f[1], f[3])
276 e[5 ] = (f[1], f[2])
277 e[6 ] = (f[1], f[5])
278 e[7 ] = (f[2], f[4])
279 e[8 ] = (f[2], f[5])
280 e[9 ] = (f[3], f[4])
281 e[10] = (f[3], f[5])
282 e[11] = (f[4], f[5])
283 self.edge_equations = []
284 for f1, f2 in e:
285 self.edge_equations.append( exp(f1)*exp(f2) - 1 )
286
287 else:
288 raise RuntimeError("Unknown polygon type %s." % name)
289
290 sfc_debug("Leaving UFCCell.__init__")
291
293 "Find which cell entity the coordinate xi lies on."
294 for i in range(self.nsd):
295 for j in range(self.num_entities[i]):
296 if self.entity_check(i, j, xi):
297 return (i, j)
298 return (self.nsd, 0)
299
301
302
303
304 if isinstance(p, swiginac.matrix):
305 p = [p[k] for k in range(len(p))]
306 elif isinstance(p, swiginac.basic):
307 p = [p]
308
309
310 if d == 0:
311
312 return bool(self.vertices[i] == matrix(self.nsd, 1, p))
313
314
315 if d == -1:
316 eq = self.facet_equations[i]
317 if d == 1:
318 eq = self.edge_equations[i]
319 if d == 2:
320 eq = self.face_equations[i]
321
322
323 x = symbols(["x", "y", "z"])
324 for j in range(len(p)):
325 eq = eq.subs(x[j] == p[j])
326 return inner(eq, eq).expand().is_zero()
327
330
333
336
339
341 if self.shape == other.shape:
342 return True
343 return False
344
346 if self.shape == other.shape:
347 return False
348 return True
349
351 s = "Cell\n"
352 s += " shape: %s\n" % self.shape
353 s += " nsd: %d\n" % self.nsd
354 s += " num_vertices: %d\n" % self.num_vertices
355 s += " num_edges: %d\n" % self.num_edges
356 s += " num_faces: %d\n" % self.num_faces
357 s += " num_facets: %d\n" % self.num_facets
358 s += " vertices: %s\n" % str(self.vertices)
359 s += " facet_equations: %s\n" % str(self.facet_equations)
360 return s
361
362
363 if __name__ == "__main__":
364 print ""
365 polygon = SyFi.ReferenceLine()
366 cell = UFCCell(polygon)
367 print cell
368
369 print ""
370 polygon = SyFi.ReferenceTriangle()
371 cell = UFCCell(polygon)
372 print cell
373
374 print ""
375 polygon = SyFi.ReferenceTetrahedron()
376 cell = UFCCell(polygon)
377 print cell
378