1
2
3 """
4 This module contains a class SFCCell to represent the
5 properties of a cell in a easily accessible way.
6 """
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28 import swiginac
29 from swiginac import matrix, exp, sqrt
30
31 import SyFi
32
33 from sfc.symbolic_utils import cross, inner, symbols
34 from sfc.common.output import sfc_debug, sfc_error
35
38 assert isinstance(polygon, SyFi.Polygon)
39 self.polygon = polygon
40 if 1:
41 x, y, z = symbols("xyz")
42 self.xi = (x,y,z)
43 else:
44 d = 3
45 x, y, z = symbols("xi[%d]" % i for i in range(d))
46 self.xi = (x,y,z)
47
49 "Find which cell entity the coordinate xi lies on."
50 for i in range(self.nsd):
51 for j in range(self.num_entities[i]):
52 if self.entity_check(i, j, xi):
53 return (i, j)
54 return (self.nsd, 0)
55
57
58
59
60 if isinstance(p, swiginac.matrix):
61 p = [p[k] for k in range(len(p))]
62 elif isinstance(p, swiginac.basic):
63 p = [p]
64
65
66 if d == 0:
67
68 return bool(self.vertices[i] == matrix(self.nsd, 1, p))
69
70
71 if d == -1:
72 eq = self.facet_equations[i]
73 if d == 1:
74 eq = self.edge_equations[i]
75 if d == 2:
76 eq = self.face_equations[i]
77
78
79 x = symbols(["x", "y", "z"])
80 for j in range(len(p)):
81 eq = eq.subs(x[j] == p[j])
82 return inner(eq, eq).expand().is_zero()
83
86
89
92
95
97 if self.shape == other.shape:
98 return True
99 return False
100
102 if self.shape == other.shape:
103 return False
104 return True
105
107 s = "Cell\n"
108 s += " shape: %s\n" % self.shape
109 s += " nsd: %d\n" % self.nsd
110 s += " num_vertices: %d\n" % self.num_vertices
111 s += " num_edges: %d\n" % self.num_edges
112 s += " num_faces: %d\n" % self.num_faces
113 s += " num_facets: %d\n" % self.num_facets
114 s += " reference_volume: %d\n" % self.reference_volume
115 s += " vertices: %s\n" % str(self.vertices)
116 s += " facet_equations: %s\n" % str(self.facet_equations)
117 return s
118
121 SFCCellBase.__init__(self, polygon)
122 x = self.xi[0]
123
124 self.shape = "interval"
125 self.facet_shape = "point"
126
127
128 self.nsd = 1
129 self.num_vertices = 2
130 self.num_edges = 1
131 self.num_faces = 0
132 self.num_facets = self.num_vertices
133 self.num_entities = (self.num_vertices, self.num_edges)
134
135
136 self.facet_vertices = [(0,), (1,)]
137
138
139 self.vertices = [matrix(self.nsd, 1, self.polygon.vertex(i))
140 for i in range(self.num_vertices)]
141 self.facet_polygons = self.vertices
142
143
144 self.reference_volume = 1.0
145
146
147 n = matrix(1, 1, [1])
148 self.facet_n = [-n, +n]
149
150
151 self.facet_equations = [x, x-1]
152
155 SFCCellBase.__init__(self, polygon)
156 x = self.xi[0]
157 y = self.xi[1]
158
159 self.shape = "triangle"
160 self.facet_shape = "interval"
161
162
163 self.nsd = 2
164 self.num_vertices = 3
165 self.num_edges = 3
166 self.num_faces = 1
167 self.num_facets = self.num_edges
168 self.num_entities = (self.num_vertices, self.num_edges, self.num_faces)
169
170
171
172 self.facet_vertices = [
173 (1, 2),
174 (2, 0),
175 (0, 1)
176 ]
177
178 self.vertices = [matrix(self.nsd, 1, self.polygon.vertex(i))
179 for i in range(self.num_vertices)]
180 self.facet_polygons = [self.polygon.line(i) for i in range(self.num_facets)]
181
182
183 self.reference_volume = 0.5
184
185
186 self.facet_n = []
187 for facet in range(self.num_facets):
188 fromvert, tovert = self.facet_vertices[facet]
189
190
191 t = self.vertices[fromvert] - self.vertices[tovert]
192 n = matrix(2, 1, [t[1], -t[0]])
193 n = n / sqrt(inner(n,n))
194
195 self.facet_n.append(n)
196
197
198 self.edge_equations = [x+y-1, x, y]
199 self.facet_equations = self.edge_equations
200
203 SFCCellBase.__init__(self, polygon)
204 x = self.xi[0]
205 y = self.xi[1]
206 z = self.xi[2]
207
208 self.shape = "tetrahedron"
209 self.facet_shape = "triangle"
210
211
212 self.nsd = 3
213 self.num_vertices = 4
214 self.num_edges = 6
215 self.num_faces = 4
216 self.num_facets = self.num_faces
217 self.num_entities = (self.num_vertices, self.num_edges, self.num_faces, 1)
218
219
220
221 self.facet_vertices = [
222 (1, 2, 3),
223 (0, 3, 2),
224 (0, 1, 3),
225 (0, 2, 1)
226 ]
227
228
229 self.vertices = [matrix(self.nsd, 1, self.polygon.vertex(i))
230 for i in range(self.num_vertices)]
231 self.facet_polygons = [self.polygon.triangle(i) for i in range(self.num_facets)]
232
233
234 self.reference_volume = 1.0/6.0
235
236
237 self.facet_n = []
238 for facet in range(self.num_facets):
239 vert = self.facet_vertices[facet]
240
241 t01 = self.vertices[vert[1]] - self.vertices[vert[0]]
242 t02 = self.vertices[vert[2]] - self.vertices[vert[0]]
243 n = cross(t01, t02)
244 n = n / sqrt(inner(n,n))
245
246 self.facet_n.append(n)
247
248
249
250
251
252
253
254
255
256 self.face_equations = [ (x+y+z-1), x, y, z]
257 self.facet_equations = self.face_equations
258
259
260 p = matrix(3, 1, [x, y, z])
261 v0, v1, v2, v3 = self.vertices
262 self.edge_equations = [ cross( (p - v3), (v2 - v3) ),
263 cross( (p - v3), (v1 - v3) ),
264 cross( (p - v2), (v1 - v2) ),
265 cross( (p - v3), (v0 - v3) ),
266 cross( (p - v2), (v0 - v2) ),
267 cross( (p - v1), (v0 - v1) ) ]
268
271 SFCCellBase.__init__(self, polygon)
272 x = self.xi[0]
273 y = self.xi[1]
274
275 self.shape = "quadrilateral"
276 self.facet_shape = "interval"
277
278
279 self.nsd = 2
280 self.num_vertices = 4
281 self.num_edges = 4
282 self.num_faces = 1
283 self.num_facets = self.num_edges
284 self.num_entities = (self.num_vertices, self.num_edges, self.num_faces)
285
286
287
288 self.facet_vertices = [
289 (2, 3),
290 (1, 2),
291 (3, 0),
292 (0, 1)
293 ]
294
295
296 self.vertices = [matrix(self.nsd, 1, self.polygon.vertex(i))
297 for i in range(self.num_vertices)]
298 self.facet_polygons = [self.polygon.line(i) for i in range(self.num_facets)]
299
300
301 self.reference_volume = 1.0
302
303
304 self.facet_n = []
305 for facet in range(self.num_facets):
306 fromvert, tovert = self.facet_vertices[facet]
307
308 t = self.vertices[fromvert] - self.vertices[tovert]
309 n = matrix(2, 1, [t[1], -t[0]])
310
311 self.facet_n.append(n)
312
313
314 self.edge_equations = [x-1, y-1, x, y]
315 self.facet_equations = self.edge_equations
316
319 SFCCellBase.__init__(self, polygon)
320 x = self.xi[0]
321 y = self.xi[1]
322 z = self.xi[2]
323
324 self.shape = "hexahedron"
325 self.facet_shape = "quadrilateral"
326
327
328 self.nsd = 3
329 self.num_vertices = 8
330 self.num_edges = 12
331 self.num_faces = 6
332 self.num_facets = self.num_faces
333 self.num_entities = (self.num_vertices, self.num_edges, self.num_faces, 1)
334
335
336
337 self.facet_vertices = [
338 (4, 5, 6, 7),
339 (2, 3, 7, 6),
340 (1, 2, 6, 5),
341 (0, 4, 7, 3),
342 (0, 1, 5, 4),
343 (0, 3, 2, 1)
344 ]
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366 self.vertices = [matrix(self.nsd, 1, self.polygon.vertex(i))
367 for i in range(self.num_vertices)]
368 self.facet_polygons = [self.polygon.rectangle(i) for i in range(self.num_facets)]
369
370
371 self.reference_volume = 1.0
372
373
374 self.facet_n = []
375 for facet in range(self.num_facets):
376
377 vert = self.facet_vertices[facet]
378
379
380 t02 = self.vertices[vert[0]] - self.vertices[vert[2]]
381 t13 = self.vertices[vert[1]] - self.vertices[vert[3]]
382 n = cross(t02, t13)
383
384
385 self.facet_n.append(n)
386
387
388
389 self.face_equations = [z-1, y-1, x-1, x, y, z]
390 self.facet_equations = self.face_equations
391
392
393 f = self.face_equations
394 e = [0]*12
395 e[0 ] = (f[0], f[1])
396 e[1 ] = (f[0], f[2])
397 e[2 ] = (f[0], f[3])
398 e[3 ] = (f[0], f[4])
399 e[4 ] = (f[1], f[3])
400 e[5 ] = (f[1], f[2])
401 e[6 ] = (f[1], f[5])
402 e[7 ] = (f[2], f[4])
403 e[8 ] = (f[2], f[5])
404 e[9 ] = (f[3], f[4])
405 e[10] = (f[3], f[5])
406 e[11] = (f[4], f[5])
407 self.edge_equations = []
408 for f1, f2 in e:
409 self.edge_equations.append( exp(f1)*exp(f2) - 1 )
410
412 if isinstance(polygon, SyFi.ReferenceLine) or polygon.str() == "ReferenceLine":
413 return SFCInterval(polygon)
414 elif isinstance(polygon, SyFi.ReferenceTriangle) or polygon.str() == "ReferenceTriangle":
415 return SFCTriangle(polygon)
416 elif isinstance(polygon, SyFi.ReferenceTetrahedron) or polygon.str() == "ReferenceTetrahedron":
417 return SFCTetrahedron(polygon)
418 elif isinstance(polygon, SyFi.ReferenceRectangle) or polygon.str() == "ReferenceRectangle":
419 return SFCQuadrilateral(polygon)
420 elif isinstance(polygon, SyFi.ReferenceBox) or polygon.str() == "ReferenceBox":
421 return SFCQuadrilateral(polygon)
422 else:
423 sfc_error("Unknown polygon type %s." % polygon.str())
424
425
426
428 print ""
429 polygon = SyFi.ReferenceLine()
430 cell = SFCCell(polygon)
431 print cell
432
433 print ""
434 polygon = SyFi.ReferenceTriangle()
435 cell = SFCCell(polygon)
436 print cell
437
438 print ""
439 polygon = SyFi.ReferenceTetrahedron()
440 cell = SFCCell(polygon)
441 print cell
442
443 if __name__ == "__main__":
444 test_ufc_cell()
445