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

Source Code for Module sfc.representation.formrepresentation

  1  # -*- coding: utf-8 -*- 
  2  """ 
  3  This module contains a representation class for forms, storing information about all form aspects. 
  4  """ 
  5   
  6  # Copyright (C) 2008-2009 Martin Sandve Alnes and Simula Resarch Laboratory 
  7  # 
  8  # This file is part of SyFi. 
  9  # 
 10  # SyFi is free software: you can redistribute it and/or modify 
 11  # it under the terms of the GNU General Public License as published by 
 12  # the Free Software Foundation, either version 2 of the License, or 
 13  # (at your option) any later version. 
 14  # 
 15  # SyFi is distributed in the hope that it will be useful, 
 16  # but WITHOUT ANY WARRANTY; without even the implied warranty of 
 17  # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 
 18  # GNU General Public License for more details. 
 19  # 
 20  # You should have received a copy of the GNU General Public License 
 21  # along with SyFi. If not, see <http://www.gnu.org/licenses/>. 
 22  # 
 23  # First added:  2008-08-13 
 24  # Last changed: 2009-04-19 
 25   
 26  import copy 
 27  import math 
 28   
 29  import swiginac 
 30  import SyFi 
 31   
 32  import ufl 
 33  from ufl.algorithms import extract_max_quadrature_element_degree,\ 
 34      estimate_quadrature_degree, estimate_total_polynomial_degree 
 35   
 36  from sfc.common import sfc_assert 
 37  from sfc.common.names import finite_element_classname, dof_map_classname, \ 
 38                             integral_classname, form_classname 
 39  from sfc.common.utilities import matrix_to_list, list_to_vector, index_string, dot_product 
 40  from sfc.geometry import UFCCell, geometry_mapping 
 41  from sfc.quadrature import find_quadrature_rule 
 42  from sfc.symbolic_utils import symbolic_matrix, symbolic_vector, symbol, symbols, cross, sub, inner 
 43  from sfc.representation.elementrepresentation import create_syfi_polygon 
 44   
45 -class FormRepresentation(object):
46 - def __init__(self, formdata, element_reps, options):
47 48 # UFL data about the form 49 self.formdata = formdata 50 51 # Mapping : ufl.FiniteElement -> sfc.ElementRepresentation 52 self.element_reps = element_reps 53 54 # SFC options 55 self.options = options 56 57 # Fetch some common variables from formdata for shorter names in rest of code 58 self.form = self.formdata.preprocessed_form 59 self.rank = self.formdata.rank 60 self.num_coefficients = self.formdata.num_coefficients 61 self.domain = self.formdata.cell.domain() 62 63 # Some strings 64 self.signature = repr(self.form) 65 self._define_classnames() 66 67 # Build some symbols and expressions 68 self._build_geometry_tokens() 69 self._build_argument_tokens() 70 71 self._define_quadrature_rule() 72 73 # Build running count of all unique ufl elements 74 # mapping : ufl_element -> running count of unique elements 75 self.element_number = {} 76 for e in self.formdata.sub_elements: 77 sfc_assert(e in self.element_reps, "Missing ElementRepresentation!") 78 if not e in self.element_number: 79 self.element_number[e] = len(self.element_number)
80
81 - def _define_quadrature_rule(self):
82 # Get integration order from options or elements 83 d = self.options.code.integral.integration_order 84 if d is None: 85 d = extract_max_quadrature_element_degree(self.form) 86 if d is None: 87 d = estimate_quadrature_degree(self.form) # TODO: Option for which estimation method to choose! 88 if d is None: 89 d = estimate_total_polynomial_degree(self.form) 90 sfc_assert(d is not None, "All quadrature degree estimation rules failed, this shouldn't happen!") 91 self.quad_order = d 92 93 # Determine quadrature rule from options or elements 94 if self.options.code.integral.integration_method == "quadrature": 95 self.quad_rule = find_quadrature_rule(self.domain, self.quad_order) 96 self._build_facet_quadrature_rules() 97 else: 98 self.quad_rule = None 99 self.facet_quad_rules = None
100
102 # pick a rule 103 cell = self.cell 104 reference_quad_rule = find_quadrature_rule(cell.facet_shape, self.quad_order, family="gauss") 105 106 self.facet_quad_rules = [] 107 for facet in range(cell.num_facets): 108 # facet polygon on the reference cell 109 facet_polygon = cell.facet_polygons[facet] 110 111 # TODO: That dimensions get right here is rather magical, 112 # I suspect that initSyFi plays a major role here... 113 # affine mapping between reference 114 facet_G, facet_x0 = geometry_mapping(facet_polygon) 115 116 # scaling factor for diagonal facets, 117 # for scaling quadrature rule to a 118 # different reference domain, this is NOT the 119 # same as the global/local scaling factor facet_D! 120 D = 1.0 121 if facet == 0: 122 if cell.shape == "triangle": 123 D = math.sqrt(2) 124 elif cell.shape == "tetrahedron": 125 D = math.sqrt(3) 126 127 # scale weights by scaling factor 128 weights = [ float(w)*D for w in reference_quad_rule.weights ] 129 130 # apply affine mapping to quadrature points 131 points = [] 132 for p in reference_quad_rule.points: 133 # extend (n-1)D point p into nD point x: 134 x = swiginac.matrix(len(facet_x0), 1, list(p)+[0.0]) 135 # apply mapping in nD space: 136 newx = facet_G * x + facet_x0 137 # store as a list: 138 xlist = [float(newx[i].evalf()) for i in range(len(newx))] 139 points.append(xlist) 140 141 # make a copy of reference quadrature rule and modify it: 142 facet_quad_rule = copy.copy(reference_quad_rule) 143 144 if cell.shape != "interval": 145 facet_quad_rule.nsd += 1 146 assert facet_quad_rule.nsd == cell.nsd 147 facet_quad_rule.comment += "\nMapped to polygon in higher dimension.\n" 148 facet_quad_rule.weights = weights 149 facet_quad_rule.points = points 150 self.facet_quad_rules.append(facet_quad_rule)
151
152 - def _define_classnames(self):
153 # Generate names for all classes 154 fd = self.formdata 155 self.classname = form_classname(self.form, self.options) 156 self.fe_names = [finite_element_classname(e) for e in fd.elements] 157 self.dm_names = [dof_map_classname(e) for e in fd.elements] 158 159 self.itg_names = dict((itg, integral_classname(itg, self.classname)) \ 160 for itg in self.form.integrals()) 161 f = self.form 162 self.citg_names = dict((itg.measure().domain_id(), self.itg_names[itg]) for itg in f.cell_integrals()) 163 self.eitg_names = dict((itg.measure().domain_id(), self.itg_names[itg]) for itg in f.exterior_facet_integrals()) 164 self.iitg_names = dict((itg.measure().domain_id(), self.itg_names[itg]) for itg in f.interior_facet_integrals())
165
166 - def _build_geometry_tokens(self):
167 "Build tokens for variables derived from cell." 168 169 self.polygon = create_syfi_polygon(self.domain) 170 self.cell = UFCCell(self.polygon) 171 172 cell = self.cell 173 nsd = cell.nsd 174 175 # vx[i][j] = component j of coordinate of vertex i in cell 176 self.vx_sym = [] 177 self.vx_expr = [] 178 for i in range(cell.num_vertices): 179 s = swiginac.matrix(nsd, 1, symbols("vx%d_%d" % (i,j) for j in range(nsd))) 180 e = swiginac.matrix(nsd, 1, symbols("c.coordinates[%d][%d]" % (i,j) for j in range(nsd))) 181 self.vx_sym.append(s) 182 self.vx_expr.append(e) 183 184 # Build a global polygon from the vertex coordinates: 185 # (using Triangle and Tetrahedron for the hypercube affine mappings 186 # is correct for rectangular and trapezoidal shapes, while using the polygons 187 # Box and Rectangle from SyFi would only support straight rectangular shapes) 188 p = [matrix_to_list(vx) for vx in self.vx_sym] 189 if cell.shape == "interval": polygon = SyFi.Line(*p) 190 elif cell.shape == "triangle": polygon = SyFi.Triangle(*p) 191 elif cell.shape == "tetrahedron": polygon = SyFi.Tetrahedron(*p) 192 elif cell.shape == "quadrilateral": polygon = SyFi.Rectangle(*p) # SyFi.Triangle(p[0], p[1], p[3]) #SyFi.Rectangle(p[0], p[2]) # TODO: Better affine mapping? 193 elif cell.shape == "hexahedron": polygon = SyFi.Box(*p) # SyFi.Tetrahedron(p[0], p[1], p[3], p[4]) # SyFi.Box(p[0], p[6]) # TODO: Better affine mapping? 194 self.global_polygon = polygon 195 #self.global_cell = UFCCell(polygon) 196 197 # Geometry mapping 198 self.G_sym = symbolic_matrix(nsd, nsd, "G") # FIXME: Make sure we don't mess up transpose here again... 199 self.x0_sym = self.vx_sym[0] #symbolic_vector(nsd, "x0") # FIXME: Make sure this is consistent. 200 self.G_expr, self.x0_expr = geometry_mapping(self.global_polygon) # TODO: Non-affine mappings? 201 202 self.detGtmp_sym = symbol("detGtmp") 203 self.detGtmp_expr = swiginac.determinant(self.G_sym) 204 205 self.detG_sym = symbol("detG") 206 self.detG_expr = swiginac.abs(self.detGtmp_sym) 207 208 # Sign of det G, negative if cell is inverted 209 self.detG_sign_sym = symbol("detG_sign") 210 self.detG_sign_expr = self.detG_sym / swiginac.abs(self.detG_sym) 211 212 # Inverse of geometry mapping (expression simplification using cofactor) 213 self.Ginv_sym = symbolic_matrix(nsd, nsd, "Ginv") 214 self.Ginv_expr = (self.detGtmp_expr*swiginac.inverse(self.G_sym)) / self.detGtmp_sym 215 #self.Ginv_expr = swiginac.inverse(self.G_sym) 216 217 # Local and global coordinates 218 self.xi_sym = symbols(("x","y","z")[:nsd]) 219 # TODO: Change to xiN, but x,y,z is used as local coordinate symbols in rest of SyFi. 220 #self.xi_sym = symbols(("xi0", "xi1", "xi2")[:nsd]) 221 self.x_sym = symbols(("x0","x1","x2")[:nsd]) 222 x_expr = [sum(self.G_sym[i,j]*self.xi_sym[j] for j in range(nsd)) for i in range(nsd)] # FIXME: Transpose? 223 self.x_expr = swiginac.matrix(nsd, 1, x_expr) 224 225 # Quadrature variables (coordinates are self.xi_sym?) 226 self.quad_weight_sym = symbol("quad_weight") 227 self.D_sym = symbol("D") 228 229 # Facet normal expressions with sign scaling for inverted cells 230 self.n_sym = symbolic_vector(nsd, "n") 231 self.n_expr = [self.detG_sign_sym*n for n in self.cell.facet_n] 232 233 # Facet mapping 234 #self.facet_G_sym = FIXME 235 #self.facet_G_expr[facet] = FIXME 236 237 self.facet_D_sym = symbol("facet_D") 238 self.facet_D_expr = [] 239 240 # --- Compute facet_D 241 sqrt = swiginac.sqrt 242 if cell.shape == "interval": 243 self.facet_D_expr = [swiginac.numeric(1) for facet in range(cell.num_facets)] 244 245 elif cell.shape == "triangle": 246 for facet in range(cell.num_facets): 247 facet_polygon = cell.facet_polygons[facet] 248 # facet_polygon is a line 249 v0 = self.G_sym * list_to_vector(facet_polygon.vertex(0)) 250 v1 = self.G_sym * list_to_vector(facet_polygon.vertex(1)) 251 v = sub(v1, v0) 252 D = sqrt( inner(v, v) ) # |v| 253 if facet == 0: # diagonal facet 254 D = D / sqrt(2) # FIXME: Is this right? 255 # TODO: Should in general use determinant of mapping instead of this linear approximation 256 self.facet_D_expr.append(D) 257 258 elif cell.shape == "tetrahedron": 259 for facet in range(cell.num_facets): 260 facet_polygon = cell.facet_polygons[facet] 261 # facet_polygon is a triangle 262 v0 = facet_polygon.vertex(0) 263 v1 = facet_polygon.vertex(1) 264 v2 = facet_polygon.vertex(2) 265 # map local reference coordinates to get global cell coordinates 266 v0 = self.G_sym * list_to_vector(v0) 267 v1 = self.G_sym * list_to_vector(v1) 268 v2 = self.G_sym * list_to_vector(v2) 269 270 # vx[i][j] = component j of coordinate of vertex i in cell 271 #self.vx_sym[i][j] 272 273 # A = |a x b|/2, D = A/(1/2) = |a x b|, 274 # with A = global triangle area, and 275 # D being A over the reference area 1/2 276 c = cross(sub(v1, v0), sub(v2, v0)) 277 D = sqrt( inner(c, c) ) # |c| 278 if facet == 0: # diagonal facet 279 D = D / sqrt(3) # FIXME: Is this right? 280 # TODO: Should in general use determinant of mapping instead of this linear approximation 281 self.facet_D_expr.append(D) 282 283 elif cell.shape == "quadrilateral": 284 for facet in range(cell.num_facets): 285 facet_polygon = cell.facet_polygons[facet] 286 # facet_polygon is a line 287 v0 = self.G_sym * list_to_vector(facet_polygon.vertex(0)) 288 v1 = self.G_sym * list_to_vector(facet_polygon.vertex(1)) 289 v = sub(v1, v0) 290 D = sqrt(inner(v,v)) 291 # TODO: Should in general use determinant of mapping instead of this linear approximation 292 self.facet_D_expr.append(D) 293 294 elif cell.shape == "hexahedron": 295 for facet in range(cell.num_facets): 296 facet_polygon = cell.facet_polygons[facet] 297 # facet_polygon is a quadrilateral 298 v0 = self.G_sym * list_to_vector(facet_polygon.vertex(0)) 299 v1 = self.G_sym * list_to_vector(facet_polygon.vertex(1)) 300 v2 = self.G_sym * list_to_vector(facet_polygon.vertex(2)) 301 v3 = self.G_sym * list_to_vector(facet_polygon.vertex(3)) 302 # compute midpoints 303 m0 = (v0 + v1)/2 304 m1 = (v1 + v2)/2 305 m2 = (v2 + v3)/2 306 m3 = (v3 + v0)/2 307 # area is length of cross product of vectors between opposing midpoints 308 c = cross(sub(m2, m0), sub(m1, m3)) 309 D = sqrt( inner(c, c) ) 310 # TODO: Should in general use determinant of mapping instead of this linear approximation 311 self.facet_D_expr.append(D)
312
313 - def _build_argument_tokens(self):
314 # --- Coefficient dofs 315 self.w_dofs = [] 316 for i in range(self.num_coefficients): 317 rep = self.element_reps[self.formdata.elements[self.rank + i]] 318 wdofs = symbols("w[%d][%d]" % (i, j) for j in range(rep.local_dimension)) 319 self.w_dofs.append(wdofs)
320
321 - def __str__(self):
322 s = "" 323 s += "FormRepresentation:\n" 324 s += "\n" # TODO: Add more here 325 s += " Geometry tokens:\n" 326 s += " self.cell = %s\n" % self.cell 327 #s += " self.global_cell = %s\n" % self.global_cell 328 s += " self.vx_sym = %s\n" % self.vx_sym 329 s += " self.vx_expr = %s\n" % self.vx_expr 330 s += " self.xi_sym = %s\n" % self.xi_sym 331 s += " self.x_sym = %s\n" % self.x_sym 332 s += " self.x_expr = %s\n" % self.x_expr 333 s += " self.G_sym = %s\n" % self.G_sym 334 s += " self.G_expr, = %s\n" % self.G_expr, 335 s += " self.detGtmp_sym = %s\n" % self.detGtmp_sym 336 s += " self.detGtmp_expr = %s\n" % self.detGtmp_expr 337 s += " self.detG_sym = %s\n" % self.detG_sym 338 s += " self.detG_expr = %s\n" % self.detG_expr 339 s += " self.Ginv_sym = %s\n" % self.Ginv_sym 340 s += " self.Ginv_expr = %s\n" % self.Ginv_expr 341 s += " self.quad_weight_sym = %s\n" % self.quad_weight_sym 342 s += " Argument tokens:\n" 343 s += " self.w_dofs = %s\n" % self.w_dofs 344 return s
345 346 #def element_representation(self, iarg): 347 # return self.element_reps[self.formdata.elements[iarg]] 348 349 #def unique_element_number(self, iarg): 350 # return self.element_number[self.formdata.elements[iarg]] 351 352 #=============================================================================== 353 # def _build_argument_cache(self): 354 # use_symbols = self.integration_method == "quadrature" 355 # 356 # self.v_cache = [] 357 # self.Dv_cache = [] 358 # for iarg in range(self.rank + self.num_coefficients): 359 # elm = self.formdata.elements[iarg] 360 # rep = self.element_reps[elm] 361 # self.v_cache.append([]) 362 # self.Dv_cache.append([]) 363 # for i in range(rep.local_dimension): 364 # for component in rep.value_components: 365 # ve = self.v_expr(iarg, i, component) 366 # vs = self.v_sym(iarg, i, component) 367 # self.v_cache[iarg][i][component] = (vs, ve) 368 # 369 # for derivatives in [(d,) for d in range(self.cell.nsd)]: 370 # Dve = self.Dv_expr(iarg, i, component, derivatives, use_symbols) 371 # Dvs = self.Dv_sym(iarg, i, component, derivatives) 372 # self.Dv_cache[iarg][i][component][derivatives] = (Dvs, Dve) 373 #=============================================================================== 374 375 # --- Basis function and coefficient data for arguments 376
377 - def v_expr(self, iarg, i, component):
378 elm = self.formdata.elements[iarg] 379 #component, elm = elm.extract_component(component) 380 381 rep = self.element_reps[elm] 382 v = rep.basis_function(i, component) 383 return v
384
385 - def v_sym(self, iarg, i, component, on_facet):
386 elm = self.formdata.elements[iarg] 387 388 e = self.v_expr(iarg, i, component) 389 e2 = e.evalf() 390 if e2.nops() == 0: 391 if e2 == 0: # FIXME: Remove this line after fixing code generation, this leads to "double 0;" statements 392 return e2 393 394 scomponent, selm = elm.extract_component(component) 395 num = self.element_number[selm] 396 397 prefix = "qt[facet][iq]." if on_facet else "qt[iq]." 398 suffix = index_string((num, i) + scomponent) 399 name = "%sv_%s" % (prefix, suffix) 400 s = symbol(name) 401 return s
402
403 - def w_expr(self, iarg, component, use_symbols, on_facet):
404 elm = self.formdata.elements[self.rank+iarg] 405 rep = self.element_reps[elm] 406 dim = rep.local_dimension 407 if use_symbols: 408 vs = [self.v_sym(iarg+self.rank, i, component, on_facet=on_facet) for i in range(dim)] 409 else: 410 vs = [self.v_expr(iarg+self.rank, i, component) for i in range(dim)] 411 w = dot_product(self.w_dofs[iarg], vs) 412 return w
413
414 - def w_sym(self, iarg, component):
415 name = "w_%s" % index_string((iarg,) + component) 416 s = symbol(name) 417 return s
418 419 # --- Basis function and coefficient data for derivatives of arguments 420
421 - def ddxi(self, f, i):
422 return swiginac.diff(f, self.xi_sym[i])
423
424 - def ddx(self, f, i):
425 return sum(self.Ginv_sym[j,i]*self.ddxi(f, j) for j in range(self.cell.nsd)) # TODO: Verify this line (in particular transposed or not!)
426
427 - def dv_expr(self, iarg, i, component, d):
428 """Return expression for dv/dxi, with v being a particular 429 component of basis function i in argument space iarg, 430 and d is a tuple (i.e. multiindex) of xi directions (local coordinates).""" 431 sfc_assert(len(d) == 1, "Higher order derivatives not implemented.") 432 d = tuple(sorted(d)) 433 v = self.v_expr(iarg, i, component) 434 dv = v 435 for k in d: 436 dv = self.ddxi(dv, k) 437 return dv
438
439 - def dv_sym(self, iarg, i, component, d, on_facet):
440 sfc_assert(len(d) == 1, "Higher order derivatives not implemented.") 441 d = tuple(sorted(d)) 442 443 e = self.dv_expr(iarg, i, component, d) 444 e2 = e.evalf() 445 if e2.nops() == 0: 446 if e2 == 0: # FIXME: Remove this line after fixing code generation, this leads to "double 0;" statements 447 return e2 448 449 prefix = "qt[facet][iq]." if on_facet else "qt[iq]." 450 suffix = index_string((iarg, i) + component + d) 451 name = "%sdv_dxi_%s" % (prefix, suffix) 452 return symbol(name)
453
454 - def Dv_expr(self, iarg, i, component, d, use_symbols, on_facet):
455 """Return expression for dv/dx, with v being a particular 456 component of basis function i in argument space iarg, 457 and d is a tuple (i.e. multiindex) of x directions (global coordinates).""" 458 sfc_assert(len(d) == 1, "Higher order derivatives not implemented.") 459 d = tuple(sorted(d)) 460 # Get local derivatives in all directions 461 if use_symbols: 462 dv = [self.dv_sym(iarg, i, component, (j,), on_facet=on_facet) for j in range(self.cell.nsd)] 463 else: 464 dv = [self.dv_expr(iarg, i, component, (j,)) for j in range(self.cell.nsd)] 465 # Apply mapping to dv 466 j, = d 467 Dv = sum(self.Ginv_sym[k, j]*dv[k] for k in range(self.cell.nsd)) # TODO: Verify this line (in particular transposed or not!) 468 return Dv
469
470 - def Dv_sym(self, iarg, i, component, d, on_facet):
471 sfc_assert(len(d) == 1, "Higher order derivatives not implemented.") 472 d = tuple(sorted(d)) 473 474 e = self.Dv_expr(iarg, i, component, d, False, on_facet=on_facet) 475 e2 = e.evalf() 476 if e2.nops() == 0: 477 if e2 == 0: # FIXME: Remove this line after fixing code generation, this leads to "double 0;" statements 478 return e2 479 480 name = "Dv_%s" % index_string((iarg, i) + component + d) 481 return symbol(name)
482
483 - def Dw_expr(self, iarg, component, d, use_symbols, on_facet):
484 sfc_assert(len(d) == 1, "Higher order derivatives not implemented.") 485 d = tuple(sorted(d)) 486 elm = self.formdata.elements[self.rank+iarg] 487 rep = self.element_reps[elm] 488 dim = rep.local_dimension 489 if use_symbols: 490 vs = [self.Dv_sym(self.rank+iarg, i, component, d, on_facet=on_facet) for i in range(dim)] 491 else: 492 vs = [self.Dv_expr(self.rank+iarg, i, component, d, use_symbols=False, on_facet=on_facet) for i in range(dim)] 493 w = dot_product(self.w_dofs[iarg], vs) 494 return w
495
496 - def Dw_sym(self, iarg, component, d):
497 sfc_assert(len(d) == 1, "Higher order derivatives not implemented.") 498 d = tuple(sorted(d)) 499 name = "Dw_%s" % index_string((iarg,) + component + d) 500 return symbol(name)
501