Package sfc :: Package representation :: Module geometryrepresentation
[hide private]
[frames] | no frames]

Source Code for Module sfc.representation.geometryrepresentation

  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  # TODO: Replace nsd with topological and geometric dimensions 
 10  # everywhere for clarity and future flexibility 
 11   
12 -def create_syfi_polygon(cell):
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
26 -def build_global_polygon(shape, p):
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 #SyFi.Triangle(p[0], p[1], p[3]) 35 #SyFi.Rectangle(p[0], p[2]) # TODO: Better affine mapping? 36 polygon = SyFi.Rectangle(*p) 37 elif shape == "hexahedron": 38 # SyFi.Tetrahedron(p[0], p[1], p[3], p[4]) 39 # SyFi.Box(p[0], p[6]) # TODO: Better affine mapping? 40 polygon = SyFi.Box(*p) 41 return polygon
42 43
44 -class GeometryRepresentation(object):
45 - def __init__(self, ufl_cell_name):
46 self.polygon_name = ufl_cell_name 47 self._build_geometry_tokens()
48
49 - def __str__(self):
50 s = "" 51 s += " Geometry tokens:\n" 52 s += " self.sfc_cell = %s\n" % self.sfc_cell 53 #s += " self.global_cell = %s\n" % self.global_cell 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
70 - def _build_geometry_tokens(self):
71 "Build tokens for variables derived from cell." 72 73 # FIXME: Use unambiguous names: ufl_cell, sfc_cell, syfi_polygon, global_syfi_polygon 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 # FIXME: Rename all references to this member to sfc_cell 78 79 sfc_cell = self.sfc_cell 80 nsd = sfc_cell.nsd 81 82 # vx[i][j] = component j of coordinate of vertex i in cell 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 # Build a global polygon from the vertex coordinates: 93 # (using Triangle and Tetrahedron for the hypercube affine mappings 94 # is correct for rectangular and trapezoidal shapes, while using the polygons 95 # Box and Rectangle from SyFi would only support straight rectangular shapes) 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 # Geometry mapping 100 # FIXME: Make sure we don't mess up transpose here again: 101 self.G_sym = symbolic_matrix(nsd, nsd, "G") 102 # FIXME: Make sure this is consistent: 103 self.x0_sym = self.vx_sym[0] #symbolic_vector(nsd, "x0") 104 # TODO: Non-affine mappings: 105 self.G_expr, self.x0_expr = geometry_mapping(self.global_polygon) 106 107 # FIXME: Rename so that detG includes the sign and detGabs is the absolute, detGtmp/detG is unclear 108 self.detGtmp_sym = symbol("detGtmp") 109 self.detGtmp_expr = swiginac.determinant(self.G_sym) 110 self.detG_sym = symbol("detG") 111 #self.detG_expr = self.detGtmp_sym 112 self.detG_expr = swiginac.abs(self.detGtmp_sym) 113 114 # Sign of det G, negative if cell is inverted 115 self.detG_sign_sym = symbol("detG_sign") 116 self.detG_sign_expr = self.detG_sym / swiginac.abs(self.detG_sym) 117 118 # Inverse of geometry mapping (expression simplification using cofactor) 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 #self.Ginv_expr = swiginac.inverse(self.G_sym) 123 124 # Local and global coordinates (getting local xi names from cell) 125 # TODO: Change in cell to "xi[%d]", then we can change x to "x[%d]". 126 self.xi_sym = swiginac.matrix(nsd, 1, list(sfc_cell.xi[:nsd])) 127 #self.x_sym = symbols(("x[%d]" % i for i in range(nsd))) 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)] # FIXME: Transpose? 131 self.x_sym = swiginac.matrix(nsd, 1, x_sym) 132 self.x_expr = swiginac.matrix(nsd, 1, x_expr) 133 134 # Quadrature variables (local coordinates are self.xi_sym) 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 # Facet mapping 140 #self.facet_G_sym = FIXME 141 #self.facet_G_expr[facet] = FIXME 142 143 # Facet normal expressions with sign scaling for inverted cells 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
158 -def build_facet_D(sfc_cell, G_sym):
159 sh = sfc_cell.shape 160 fp = sfc_cell.facet_polygons 161 if sh == "interval": 162 return build_facet_D_interval(fp, G_sym) 163 elif sh == "triangle": 164 return build_facet_D_triangle(fp, G_sym) 165 elif sh == "tetrahedron": 166 return build_facet_D_tetrahedron(fp, G_sym) 167 elif sh == "quadrilateral": 168 return build_facet_D_quadrilateral(fp, G_sym) 169 elif sh == "hexahedron": 170 return build_facet_D_hexahedron(fp, G_sym)
171
172 -def build_facet_D_interval(facet_polygons, G_sym):
173 one = swiginac.numeric(1) 174 return [one]*len(facet_polygons)
175
176 -def build_facet_D_triangle(facet_polygons, G_sym):
177 facet_D_expr = [] 178 for facet, facet_polygon in enumerate(facet_polygons): 179 # facet_polygon is a line 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) ) # |v| 184 if facet == 0: # diagonal facet 185 D = D / swiginac.sqrt(2) # FIXME: Is this right? 186 # TODO: Should in general use determinant of mapping 187 # instead of this linear approximation 188 facet_D_expr.append(D) 189 return facet_D_expr
190
191 -def build_facet_D_tetrahedron(facet_polygons, G_sym):
192 facet_D_expr = [] 193 for facet, facet_polygon in enumerate(facet_polygons): 194 # facet_polygon is a triangle 195 v0 = facet_polygon.vertex(0) 196 v1 = facet_polygon.vertex(1) 197 v2 = facet_polygon.vertex(2) 198 # map local reference coordinates to get global cell coordinates 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 # vx[i][j] = component j of coordinate of vertex i in cell 204 #vx_sym[i][j] 205 206 # A = |a x b|/2, D = A/(1/2) = |a x b|, 207 # with A = global triangle area, and 208 # D being A over the reference area 1/2 209 c = cross(sub(v1, v0), sub(v2, v0)) 210 D = swiginac.sqrt( inner(c, c) ) # |c| 211 if facet == 0: # diagonal facet 212 D = D / swiginac.sqrt(3) # FIXME: Is this right? 213 # TODO: Should in general use determinant of mapping 214 # instead of this linear approximation 215 facet_D_expr.append(D) 216 return facet_D_expr
217
218 -def build_facet_D_quadrilateral(facet_polygons, G_sym):
219 facet_D_expr = [] 220 for facet, facet_polygon in enumerate(facet_polygons): 221 # facet_polygon is a line 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 # TODO: Should in general use determinant of mapping 227 # instead of this linear approximation 228 facet_D_expr.append(D) 229 return facet_D_expr
230
231 -def build_facet_D_hexahedron(facet_polygons, G_sym):
232 facet_D_expr = [] 233 for facet, facet_polygon in enumerate(facet_polygons): 234 # facet_polygon is a quadrilateral 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 # compute midpoints 240 m0 = (v0 + v1)/2 241 m1 = (v1 + v2)/2 242 m2 = (v2 + v3)/2 243 m3 = (v3 + v0)/2 244 # area is length of cross product of vectors between opposing midpoints 245 c = cross(sub(m2, m0), sub(m1, m3)) 246 D = swiginac.sqrt( inner(c, c) ) 247 # TODO: Should in general use determinant of mapping 248 # instead of this linear approximation 249 facet_D_expr.append(D) 250 return facet_D_expr
251