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

Source Code for Module sfc.geometry.UFCCell

  1  #!/usr/bin/env python 
  2  # -*- coding: utf-8 -*- 
  3  """ 
  4  This module contains a class UFCCell to represent the properties of a cell in a easily accessible way. 
  5  """ 
  6   
  7  # Copyright (C) 2008 Martin Sandve Alnes and Simula Resarch Laboratory 
  8  # 
  9  # This file is part of SyFi. 
 10  # 
 11  # SyFi is free software: you can redistribute it and/or modify 
 12  # it under the terms of the GNU General Public License as published by 
 13  # the Free Software Foundation, either version 2 of the License, or 
 14  # (at your option) any later version. 
 15  # 
 16  # SyFi is distributed in the hope that it will be useful, 
 17  # but WITHOUT ANY WARRANTY; without even the implied warranty of 
 18  # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 
 19  # GNU General Public License for more details. 
 20  # 
 21  # You should have received a copy of the GNU General Public License 
 22  # along with SyFi. If not, see <http://www.gnu.org/licenses/>. 
 23  # 
 24  # First added:  2008-08-12 
 25  # Last changed: 2008-08-12 
 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  # TODO: this code could do with some tests! 
 36   
37 -class UFCCell:
38 - def __init__(self, polygon):
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 # ... dimensions: 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 # ... connectivity: 59 self.facet_vertices = [ 60 (0,), 61 (1,) 62 ] 63 64 # ... geometry: 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 # ... normal vector: 69 n = matrix(1, 1, [1]) 70 self.facet_n = [-n, +n] # TODO: is this what we want to do for the interval "normal"? 71 72 # ... implicit equations: 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 # ... dimensions: 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 # ... connectivity: 88 # facet vertices in counterclockwise direction around reference cell: 89 self.facet_vertices = [ 90 (1, 2), 91 (2, 0), 92 (0, 1) 93 ] 94 # ... geometry: 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 # ... normal vector: 99 self.facet_n = [] 100 for facet in range(self.num_facets): 101 fromvert, tovert = self.facet_vertices[facet] 102 103 # FIXME: finish facet stuff 104 t = self.vertices[fromvert] - self.vertices[tovert] 105 n = matrix(2, 1, [t[1], -t[0]]) 106 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... 107 108 self.facet_n.append(n) 109 110 # ... implicit equations: 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 # ... dimensions: 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 # ... connectivity: 127 # facet vertices in counterclockwise direction around reference cell: 128 self.facet_vertices = [ 129 (1, 2, 3), 130 (0, 3, 2), # (0, 2, 3) 131 (0, 1, 3), 132 (0, 2, 1) # (0, 1, 2) 133 ] 134 135 # ... geometry: 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 # ... normal vector: 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)) # TODO: is this approach efficient? FIXME: Is this correct? Scaling to unit normal on reference domain... 148 149 self.facet_n.append(n) 150 151 # ... implicit equations: 152 self.face_equations = [ (x+y+z-1), x, y, z] 153 self.facet_equations = self.face_equations 154 155 # in a tetrahedron, no two edges are parallel WARNiNG FIXME: this edge code is not verified! 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 # ... dimensions: 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 # ... connectivity: 178 # facet vertices in counterclockwise direction around reference cell: 179 self.facet_vertices = [ 180 (2, 3), 181 (1, 2), 182 (3, 0), 183 (0, 1) 184 ] 185 186 # ... geometry: 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 # ... normal vector: 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 # ... implicit equations: 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 # ... dimensions: 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 # ... connectivity: 217 # facet vertices in counterclockwise direction around reference cell: 218 self.facet_vertices = [ 219 (4, 5, 6, 7), 220 (2, 3, 7, 6), # (2, 3, 6, 7) 221 (1, 2, 6, 5), # (1, 2, 5, 6) 222 (0, 4, 7, 3), # (0, 3, 4, 7) 223 (0, 1, 5, 4), # (0, 1, 4, 5) 224 (0, 3, 2, 1) # (0, 1, 2, 3) 225 ] 226 227 # Current Lagrange basis function order in SyFi: 228 #0 [[0, 0, 0], 0] 229 #1 [[0, 0, 1], 0] 230 #2 [[1, 0, 0], 0] 231 #3 [[1, 0, 1], 0] 232 #4 [[0, 1, 0], 0] 233 #5 [[0, 1, 1], 0] 234 #6 [[1, 1, 0], 0] 235 #7 [[1, 1, 1], 0] 236 # UFC vertex order for hexes requires reordering: TODO: don't need to do this? 237 #0->0 [[0, 0, 0], 0] 238 #2->1 [[1, 0, 0], 0] 239 #6->2 [[1, 1, 0], 0] 240 #4->3 [[0, 1, 0], 0] 241 #1->4 [[0, 0, 1], 0] 242 #3->5 [[1, 0, 1], 0] 243 #7->6 [[1, 1, 1], 0] 244 #5->7 [[0, 1, 1], 0] 245 246 # ... geometry: 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 # counterclockwise ordering seen from outside of reference cell: 253 vert = self.facet_vertices[facet] 254 255 # computing normal as the cross product of the facet diagonals 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 #n = n / swiginac.sqrt(inner(n,n)) # TODO: is this approach efficient? FIXME: Is this correct? Scaling to unit normal on reference domain... 260 261 self.facet_n.append(n) 262 263 # ... implicit equations: 264 # all faces coincide with a cartesian plane 265 self.face_equations = [z-1, y-1, x-1, x, y, z] 266 self.facet_equations = self.face_equations 267 268 # points on an edge satisfy the equations of two faces 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 ) # TODO: a better equation? 286 287 else: 288 raise RuntimeError("Unknown polygon type %s." % name) 289 290 sfc_debug("Leaving UFCCell.__init__")
291
292 - def find_entity(self, xi):
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
300 - def entity_check(self, d, i, p):
301 # this is a bit ugly, could benefit from a cleanup 302 303 # make p into a list 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 # check if we match a vertex exactly 310 if d == 0: 311 #eq = self.vertex_equations[i] 312 return bool(self.vertices[i] == matrix(self.nsd, 1, p)) 313 314 # get implicit equation for this entity 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 # check if implicit equation is zero in this point, which means p is on the entity 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
328 - def facet_check(self, i, p):
329 return self.entity_check(-1, i, p)
330
331 - def vertex_check(self, i, p):
332 return self.entity_check(0, i, p)
333
334 - def edge_check(self, i, p):
335 return self.entity_check(1, i, p)
336
337 - def face_check(self, i, p):
338 return self.entity_check(2, i, p)
339
340 - def __eq__(self, other):
341 if self.shape == other.shape: 342 return True 343 return False
344
345 - def __ne__(self, other):
346 if self.shape == other.shape: 347 return False 348 return True
349
350 - def __str__(self):
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