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

Source Code for Module sfc.geometry.sfc_cell

  1  #!/usr/bin/env python 
  2  # -*- coding: utf-8 -*- 
  3  """ 
  4  This module contains a class SFCCell to represent the 
  5  properties of a cell in a easily accessible way. 
  6  """ 
  7   
  8  # Copyright (C) 2008 Martin Sandve Alnes and Simula Resarch Laboratory 
  9  # 
 10  # This file is part of SyFi. 
 11  # 
 12  # SyFi is free software: you can redistribute it and/or modify 
 13  # it under the terms of the GNU General Public License as published by 
 14  # the Free Software Foundation, either version 2 of the License, or 
 15  # (at your option) any later version. 
 16  # 
 17  # SyFi is distributed in the hope that it will be useful, 
 18  # but WITHOUT ANY WARRANTY; without even the implied warranty of 
 19  # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 
 20  # GNU General Public License for more details. 
 21  # 
 22  # You should have received a copy of the GNU General Public License 
 23  # along with SyFi. If not, see <http://www.gnu.org/licenses/>. 
 24  # 
 25  # First added:  2008-08-12 
 26  # Last changed: 2008-08-12 
 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   
36 -class SFCCellBase(object):
37 - def __init__(self, polygon):
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
48 - def find_entity(self, xi):
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
56 - def entity_check(self, d, i, p):
57 # this is a bit ugly, could benefit from a cleanup 58 59 # make p into a list 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 # check if we match a vertex exactly 66 if d == 0: 67 #eq = self.vertex_equations[i] 68 return bool(self.vertices[i] == matrix(self.nsd, 1, p)) 69 70 # get implicit equation for this entity 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 # check if implicit equation is zero in this point, which means p is on the entity 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
84 - def facet_check(self, i, p):
85 return self.entity_check(-1, i, p)
86
87 - def vertex_check(self, i, p):
88 return self.entity_check(0, i, p)
89
90 - def edge_check(self, i, p):
91 return self.entity_check(1, i, p)
92
93 - def face_check(self, i, p):
94 return self.entity_check(2, i, p)
95
96 - def __eq__(self, other):
97 if self.shape == other.shape: 98 return True 99 return False
100
101 - def __ne__(self, other):
102 if self.shape == other.shape: 103 return False 104 return True
105
106 - def __str__(self):
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
119 -class SFCInterval(SFCCellBase):
120 - def __init__(self, polygon):
121 SFCCellBase.__init__(self, polygon) 122 x = self.xi[0] 123 124 self.shape = "interval" 125 self.facet_shape = "point" 126 127 # ... dimensions: 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 # ... connectivity: 136 self.facet_vertices = [(0,), (1,)] 137 138 # ... geometry: 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 # ... volume: 144 self.reference_volume = 1.0 145 146 # ... normal vector: 147 n = matrix(1, 1, [1]) 148 self.facet_n = [-n, +n] # TODO: is this what we want to do for the interval "normal"? 149 150 # ... implicit equations: 151 self.facet_equations = [x, x-1]
152
153 -class SFCTriangle(SFCCellBase):
154 - def __init__(self, polygon):
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 # ... dimensions: 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 # ... connectivity: 171 # facet vertices in counterclockwise direction around reference cell: 172 self.facet_vertices = [ 173 (1, 2), 174 (2, 0), 175 (0, 1) 176 ] 177 # ... geometry: 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 # ... volume: 183 self.reference_volume = 0.5 184 185 # ... normal vector: 186 self.facet_n = [] 187 for facet in range(self.num_facets): 188 fromvert, tovert = self.facet_vertices[facet] 189 190 # FIXME: finish facet stuff 191 t = self.vertices[fromvert] - self.vertices[tovert] 192 n = matrix(2, 1, [t[1], -t[0]]) 193 n = n / sqrt(inner(n,n)) # TODO: Is this approach efficient? FIXME: Is this correct? Scaling to unit normal on reference domain, but this isn't necessarily unit length on global domain... 194 195 self.facet_n.append(n) 196 197 # ... implicit equations: 198 self.edge_equations = [x+y-1, x, y] 199 self.facet_equations = self.edge_equations
200
201 -class SFCTetrahedron(SFCCellBase):
202 - def __init__(self, polygon):
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 # ... dimensions: 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 # ... connectivity: 220 # facet vertices in counterclockwise direction around reference cell: 221 self.facet_vertices = [ 222 (1, 2, 3), 223 (0, 3, 2), # (0, 2, 3) 224 (0, 1, 3), 225 (0, 2, 1) # (0, 1, 2) 226 ] 227 228 # ... geometry: 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 # ... volume: 234 self.reference_volume = 1.0/6.0 # FIXME: Is this correct? 235 236 # ... normal vector: 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)) # TODO: is this approach efficient? FIXME: Is this correct? Scaling to unit normal on reference domain... 245 246 self.facet_n.append(n) 247 248 #self.local_facet_n = [ 249 # matrix(self.nsd, 1, (1.0, 1.0, 1.0)), 250 # matrix(self.nsd, 1, (1.0, 0.0, 0.0)), 251 # matrix(self.nsd, 1, (0.0, 1.0, 0.0)), 252 # matrix(self.nsd, 1, (0.0, 0.0, 1.0)), 253 # ] 254 255 # ... implicit equations: 256 self.face_equations = [ (x+y+z-1), x, y, z] 257 self.facet_equations = self.face_equations 258 259 # in a tetrahedron, no two edges are parallel WARNiNG FIXME: this edge code is not verified! 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
269 -class SFCQuadrilateral(SFCCellBase):
270 - def __init__(self, polygon):
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 # ... dimensions: 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 # ... connectivity: 287 # facet vertices in counterclockwise direction around reference cell: 288 self.facet_vertices = [ 289 (2, 3), 290 (1, 2), 291 (3, 0), 292 (0, 1) 293 ] 294 295 # ... geometry: 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 # ... volume: 301 self.reference_volume = 1.0 302 303 # ... normal vector: 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 # ... implicit equations: 314 self.edge_equations = [x-1, y-1, x, y] 315 self.facet_equations = self.edge_equations
316
317 -class SFCHexahedron(SFCCellBase):
318 - def __init__(self, polygon):
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 # ... dimensions: 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 # ... connectivity: 336 # facet vertices in counterclockwise direction around reference cell: 337 self.facet_vertices = [ 338 (4, 5, 6, 7), 339 (2, 3, 7, 6), # (2, 3, 6, 7) 340 (1, 2, 6, 5), # (1, 2, 5, 6) 341 (0, 4, 7, 3), # (0, 3, 4, 7) 342 (0, 1, 5, 4), # (0, 1, 4, 5) 343 (0, 3, 2, 1) # (0, 1, 2, 3) 344 ] 345 346 # Current Lagrange basis function order in SyFi: 347 #0 [[0, 0, 0], 0] 348 #1 [[0, 0, 1], 0] 349 #2 [[1, 0, 0], 0] 350 #3 [[1, 0, 1], 0] 351 #4 [[0, 1, 0], 0] 352 #5 [[0, 1, 1], 0] 353 #6 [[1, 1, 0], 0] 354 #7 [[1, 1, 1], 0] 355 # UFC vertex order for hexes requires reordering: TODO: don't need to do this? 356 #0->0 [[0, 0, 0], 0] 357 #2->1 [[1, 0, 0], 0] 358 #6->2 [[1, 1, 0], 0] 359 #4->3 [[0, 1, 0], 0] 360 #1->4 [[0, 0, 1], 0] 361 #3->5 [[1, 0, 1], 0] 362 #7->6 [[1, 1, 1], 0] 363 #5->7 [[0, 1, 1], 0] 364 365 # ... geometry: 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 # ... volume: 371 self.reference_volume = 1.0 372 373 # ... normal vector: 374 self.facet_n = [] 375 for facet in range(self.num_facets): 376 # counterclockwise ordering seen from outside of reference cell: 377 vert = self.facet_vertices[facet] 378 379 # computing normal as the cross product of the facet diagonals 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 #n = n / swiginac.sqrt(inner(n,n)) # TODO: is this approach efficient? FIXME: Is this correct? Scaling to unit normal on reference domain... 384 385 self.facet_n.append(n) 386 387 # ... implicit equations: 388 # all faces coincide with a cartesian plane 389 self.face_equations = [z-1, y-1, x-1, x, y, z] 390 self.facet_equations = self.face_equations 391 392 # points on an edge satisfy the equations of two faces 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 ) # TODO: a better equation?
410
411 -def SFCCell(polygon):
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 # TODO: Test in unit test suite instead, add more tests
427 -def test_ufc_cell():
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