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

Source Code for Module sfc.representation.swiginac_eval

  1  """This module defines evaluation algorithms for converting 
  2  converting UFL expressions to swiginac representation.""" 
  3   
  4  # Copyright (C) 2008-2009 Martin Sandve Alnes and Simula Resarch Laboratory 
  5  # 
  6  # This file is part of SyFi. 
  7  # 
  8  # SyFi is free software: you can redistribute it and/or modify 
  9  # it under the terms of the GNU General Public License as published by 
 10  # the Free Software Foundation, either version 2 of the License, or 
 11  # (at your option) any later version. 
 12  # 
 13  # SyFi is distributed in the hope that it will be useful, 
 14  # but WITHOUT ANY WARRANTY; without even the implied warranty of 
 15  # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 
 16  # GNU General Public License for more details. 
 17  # 
 18  # You should have received a copy of the GNU General Public License 
 19  # along with SyFi. If not, see <http://www.gnu.org/licenses/>. 
 20  # 
 21  # Modified by Kent-Andre Mardal, 2010. 
 22  # 
 23  # First added:  2008-08-22 
 24  # Last changed: 2009-03-19 
 25   
 26  from collections import defaultdict 
 27  from itertools import izip, chain 
 28   
 29  import swiginac 
 30   
 31  from ufl import * 
 32  from ufl.classes import * 
 33  from ufl.common import some_key, product, Stack, StackDict 
 34  from ufl.algorithms.transformations import Transformer, MultiFunction 
 35  from ufl.permutation import compute_indices 
 36   
 37  from sfc.common import sfc_assert, sfc_error, sfc_warning 
 38  from sfc.symbolic_utils import symbol, symbols 
 39   
40 -class EvaluateAsSwiginac(MultiFunction):
41 - def __init__(self, formrep, itgrep, data, on_facet): # TODO: Remove unused arguments after implementing:
42 MultiFunction.__init__(self) 43 self.formrep = formrep 44 self.itgrep = itgrep 45 self.data = data 46 self.on_facet = on_facet 47 self.current_basis_function = (None,)*formrep.rank
48 49 ### Fallback handlers: 50
51 - def expr(self, o, *ops):
52 sfc_error("Evaluation not implemented for expr %s." % type(o).__name__)
53
54 - def terminal(self, o, *ops):
55 sfc_error("Evaluation not implemented for terminal %s." % type(o).__name__)
56 57 ### Terminals: 58
59 - def zero(self, o):
60 return swiginac.numeric(0)
61
62 - def scalar_value(self, o):
63 return swiginac.numeric(o.value())
64
65 - def spatial_coordinate(self, o, component=(), derivatives=()):
66 # Called by indexed 67 68 if component: 69 # 2D, 3D 70 c, = component 71 else: 72 # 1D 73 c = 0 74 75 if derivatives: 76 if len(derivatives) > 1: 77 return swiginac.numeric(0) 78 d, = derivatives 79 if d == c: 80 return swiginac.numeric(1) 81 return swiginac.numeric(0) 82 83 return self.formrep.x_sym[c]
84
85 - def facet_normal(self, o, component=(), derivatives=()):
86 # Called by indexed 87 sfc_assert(self.on_facet, "Expecting to be on a facet in facet_normal.") 88 89 if derivatives: 90 return swiginac.numeric(0) 91 92 if component: 93 # 2D, 3D 94 c, = component 95 else: 96 # 1D 97 c = 0 98 99 return self.formrep.n_sym[c]
100
101 - def argument(self, o, component=(), derivatives=()):
102 103 # Assuming renumbered arguments! 104 iarg = o.count() 105 sfc_assert(iarg >= 0, "Argument count shouldn't be negative.") 106 sfc_assert(isinstance(component, tuple), "Expecting tuple for component.") 107 108 j = self.current_basis_function[iarg] 109 110 if derivatives: 111 s = self.formrep.Dv_sym(iarg, j, component, derivatives, self.on_facet) 112 e = self.formrep.Dv_expr(iarg, j, component, derivatives, False, self.on_facet) # FIXME: use_symbols = False -> can do better 113 else: 114 s = self.formrep.v_sym(iarg, j, component, self.on_facet) 115 e = self.formrep.v_expr(iarg, j, component) 116 117 if e.nops() == 0: 118 return e # FIXME: Avoid generating code for s when not using it 119 return s
120
121 - def coefficient(self, o, component=(), derivatives=()):
122 # Assuming renumbered arguments! 123 iarg = o.count() 124 sfc_assert(iarg >= 0, "Coefficient count shouldn't be negative.") 125 sfc_assert(isinstance(component, tuple), "Expecting tuple for component.") 126 127 if derivatives: 128 s = self.formrep.Dw_sym(iarg, component, derivatives) 129 e = self.formrep.Dw_expr(iarg, component, derivatives, False, self.on_facet) # FIXME: use_symbols = False -> can do better 130 else: 131 # w^i_h(x) = \sum_j w[i][j] * phi^i_j(x) 132 s = self.formrep.w_sym(iarg, component) 133 e = self.formrep.w_expr(iarg, component, False, self.on_facet) # FIXME: use_symbols = False -> can do better 134 135 if e.nops() == 0: 136 return e # FIXME: Avoid generating code for s when not using it 137 return s
138 139 ### Indexing and derivatives: 140
141 - def multi_index(self, o):
142 return tuple(map(int, o))
143
144 - def spatial_derivative(self, o, f, i, component=(), derivatives=()):
145 derivatives = sorted(derivatives + self.multi_index(i)) 146 if isinstance(f, Indexed): 147 # Since expand_indices moves Indexed in to the terminals, 148 # SpatialDerivative can be outside an Indexed 149 sfc_assert(component == (), "Expecting no component outside of Indexed!") 150 A, ii = f.operands() 151 component = self.multi_index(ii) 152 return self(A, component, derivatives) 153 return self(f, component, derivatives)
154
155 - def indexed(self, o, A, ii):
156 # Passes on control to one of: 157 #def argument(self, o, component): 158 #def coefficient(self, o, component): 159 #def facet_normal(self, o, component): 160 #def spatial_coordinate(self, o, component): 161 #def spatial_derivative(self, o, component): 162 component = self.multi_index(ii) 163 if isinstance(A, SpatialDerivative): 164 f, i = A.operands() 165 return self.spatial_derivative(A, f, i, component) 166 return self(A, component)
167 168 ### Algebraic operators: 169
170 - def power(self, o, a, b):
171 return a**b
172
173 - def sum(self, o, *ops):
174 return sum(ops)
175
176 - def product(self, o, *ops):
177 return product(ops)
178
179 - def division(self, o, a, b):
180 return a / b
181
182 - def abs(self, o, a):
183 return swiginac.abs(a)
184 185 ### Basic math functions: 186
187 - def sqrt(self, o, a):
188 return swiginac.sqrt(a)
189
190 - def exp(self, o, a):
191 return swiginac.exp(a)
192
193 - def ln(self, o, a):
194 return swiginac.log(a)
195
196 - def cos(self, o, a):
197 return swiginac.cos(a)
198
199 - def sin(self, o, a):
200 return swiginac.sin(a)
201
202 - def tan(self, o, a):
203 return swiginac.tan(a)
204
205 - def acos(self, o, a):
206 return swiginac.acos(a)
207
208 - def asin(self, o, a):
209 return swiginac.asin(a)
210
211 - def atan(self, o, a):
212 return swiginac.atan(a)
213 214 215 216
217 - def variable(self, o):
218 sfc_error("Should strip away variables before building graph.")
219 220 # FIXME: Implement all missing operators here 221
222 -class SwiginacEvaluator(Transformer):
223 "Algorithm for evaluation of an UFL expression as a swiginac expression."
224 - def __init__(self, formrep, use_symbols, on_facet):
225 Transformer.__init__(self)#, variable_cache) 226 227 # input 228 self._formrep = formrep 229 self._use_symbols = use_symbols 230 self._on_facet = on_facet 231 232 # current basis function configuration 233 self._current_basis_function = tuple(0 for i in range(formrep.rank)) 234 235 # current indexing status 236 self._components = Stack() 237 self._index2value = StackDict() 238 239 # code and cache structures 240 # FIXME: Need pre-initialized self._variable2symbol and self._tokens 241 self._variable2symbol = {} 242 self._tokens = [] 243 244 # convenience variables 245 self.nsd = self._formrep.cell.nsd
246
247 - def pop_tokens(self):
248 # TODO: Make a generator approach to this? Allow handlers to trigger a "token yield"? 249 t = self._tokens 250 self._tokens = [] 251 return t
252
253 - def update(self, iota):
254 self._current_basis_function = tuple(iota)
255
256 - def component(self):
257 "Return current component tuple." 258 if len(self._components): 259 return self._components.peek() 260 return ()
261 262 ### Fallback handlers: 263
264 - def expr(self, x):
265 sfc_error("Missing ufl to swiginac handler for type %s" % str(type(x)))
266
267 - def terminal(self, x):
268 sfc_error("Missing ufl to swiginac handler for terminal type %s" % str(type(x)))
269 270 ### Handlers for basic terminal objects: 271
272 - def zero(self, x):
273 #sfc_assert(len(self.component()) == len(x.shape()), "Index component length mismatch in zero tensor!") 274 return swiginac.numeric(0)
275
276 - def scalar_value(self, x):
277 #sfc_assert(self.component() == (), "Shouldn't have any component at this point.") 278 return swiginac.numeric(x._value)
279
280 - def identity(self, x):
281 c = self.component() 282 v = 1 if c[0] == c[1] else 0 283 return swiginac.numeric(v)
284
285 - def argument(self, x):
286 iarg = x.count() 287 sfc_assert(iarg >= 0, "Argument count shouldn't be negative.") 288 j = self._current_basis_function[iarg] 289 c = self.component() 290 if self._use_symbols: 291 return self._formrep.v_sym(iarg, j, c, self._on_facet) 292 else: 293 return self._formrep.v_expr(iarg, j, c)
294
295 - def coefficient(self, x):
296 iarg = x.count() 297 c = self.component() 298 if self._use_symbols: 299 return self._formrep.w_sym(iarg, c) 300 else: 301 # w^i_h(x) = \sum_j w[i][j] * phi^i_j(x) 302 return self._formrep.w_expr(iarg, c, False, self._on_facet)
303
304 - def facet_normal(self, x):
305 sfc_assert(self._on_facet, "Expecting to be on a facet in facet_normal.") 306 c, = self.component() 307 return self._formrep.n_sym[c]
308
309 - def spatial_coordinate(self, x):
310 c, = self.component() 311 return self._formrep.x_sym[c]
312 313 ### Handler for variables: 314
315 - def variable(self, x):
316 return self.visit(x._expression)
317
318 - def garbage(self, x): # TODO: Maybe some of the ideas here can be used in code generation
319 c = self.component() 320 index_values = tuple(self._index2value[k] for k in x._expression.free_indices()) 321 322 # TODO: Doesn't always depend on _current_basis_function, this is crap: 323 key = (x.count(), c, index_values, self._current_basis_function) 324 vsym = self._variable2symbol.get(key) 325 326 if vsym is None: 327 expr = self.visit(x._expression) 328 # TODO: Doesn't always depend on _current_basis_function, this is crap: 329 compstr = "_".join("%d" % k for k in chain(c, index_values, self._current_basis_function)) 330 vname = "_".join(("t_%d" % x.count(), compstr)) 331 vsym = symbol(vname) 332 self._variable2symbol[key] = vsym 333 self._tokens.append((vsym, expr))
334 335 ### Handlers for basic algebra: 336
337 - def sum(self, x, *ops):
338 return sum(ops)
339
340 - def index_sum(self, x):
341 ops = [] 342 summand, multiindex = x.operands() 343 index, = multiindex 344 for i in range(x.dimension()): 345 self._index2value.push(index, i) 346 ops.append(self.visit(summand)) 347 self._index2value.pop() 348 return sum(ops)
349
350 - def product(self, x):
351 sfc_assert(not self.component(), "Non-empty indexing component in product!") 352 ops = [self.visit(o) for o in x.operands()] 353 return product(ops)
354 # ... 355
356 - def division(self, x, a, b):
357 return a / b
358
359 - def power(self, x, a, b):
360 return a ** b
361
362 - def abs(self, x, a):
363 return swiginac.abs(a)
364 365 ### Basic math functions:
366 - def sqrt(self, x, y):
367 return swiginac.sqrt(y)
368
369 - def exp(self, x, y):
370 return swiginac.exp(y)
371
372 - def ln(self, x, y):
373 return swiginac.log(y)
374
375 - def cos(self, x, y):
376 return swiginac.cos(y)
377
378 - def sin(self, x, y):
379 return swiginac.sin(y)
380 381 ### Index handling:
382 - def multi_index(self, x):
383 subcomp = [] 384 for i in x: 385 if isinstance(i, FixedIndex): 386 subcomp.append(i._value) 387 elif isinstance(i, Index): 388 subcomp.append(self._index2value[i]) 389 return tuple(subcomp)
390
391 - def indexed(self, x):
392 A, ii = x.operands() 393 self._components.push(self.visit(ii)) 394 result = self.visit(A) 395 self._components.pop() 396 return result
397 398 ### Container handling: 399
400 - def old_list_tensor(self, x): # doesn't support e.g. building a matrix from vector rows
401 component = self.component() 402 sfc_assert(len(component) > 0 and \ 403 all(isinstance(i, int) for i in component), 404 "Can't index tensor with %s." % repr(component)) 405 406 # Hide indexing when evaluating subexpression 407 self._components.push(()) 408 409 # Get scalar UFL subexpression from tensor 410 e = x 411 for i in component: 412 e = e._expressions[i] 413 sfc_assert(e.shape() == (), "Expecting scalar expression "\ 414 "after extracting component from tensor.") 415 416 # Apply conversion to scalar subexpression 417 r = self.visit(e) 418 419 # Return to previous component state 420 self._components.pop() 421 return r 422
423 - def list_tensor(self, x):
424 # Pick the right subtensor and subcomponent 425 c = self.component() 426 c0, c1 = c[0], c[1:] 427 op = x.operands()[c0] 428 # Evaluate subtensor with this subcomponent 429 self._components.push(c1) 430 r = self.visit(op) 431 self._components.pop() 432 return r
433
434 - def component_tensor(self, x):
435 # this function evaluates the tensor expression 436 # with indices equal to the current component tuple 437 expression, indices = x.operands() 438 sfc_assert(expression.shape() == (), "Expecting scalar base expression.") 439 440 # update index map with component tuple values 441 comp = self.component() 442 sfc_assert(len(indices) == len(comp), "Index/component mismatch.") 443 for i, v in izip(indices._indices, comp): 444 self._index2value.push(i, v) 445 self._components.push(()) 446 447 # evaluate with these indices 448 result = self.visit(expression) 449 450 # revert index map 451 for i in range(len(comp)): 452 self._index2value.pop() 453 self._components.pop() 454 return result
455 456 ### Differentiation: 457
458 - def _ddx(self, f, i):
459 """Differentiate swiginac expression f w.r.t. x_i, using 460 df/dx_i = df/dxi_j dxi_j/dx_i.""" 461 Ginv = self._formrep.Ginv_sym 462 xi = self._formrep.xi_sym 463 return sum(Ginv[j, i] * swiginac.diff(f, xi[j]) for j in range(self.nsd))
464
465 - def spatial_derivative(self, x):
466 # Assuming that AD has been applied, so 467 # the expression to differentiate is always a Terminal. 468 469 f, ii = x.operands() 470 471 sfc_assert(isinstance(f, Terminal), \ 472 "Expecting to differentiate a Terminal object, you must apply AD first!") # The exception is higher order derivatives, ignoring for now 473 474 # Get component and derivative directions 475 c = self.component() 476 der = self.visit(ii) 477 478 # --- Handle derivatives of basis functions 479 if isinstance(f, Argument): 480 iarg = f.count() 481 i = self._current_basis_function[iarg] 482 if self._use_symbols: 483 return self._formrep.Dv_sym(iarg, i, c, der, self._on_facet) 484 else: 485 return self._formrep.Dv_expr(iarg, i, c, der, False, self._on_facet) 486 487 # --- Handle derivatives of coefficient functions 488 if isinstance(f, Coefficient): 489 iarg = f.count() 490 if self._use_symbols: 491 return self._formrep.Dw_sym(iarg, c, der) 492 else: 493 return self._formrep.Dw_expr(iarg, c, der, False, self._on_facet) 494 495 # --- Handle derivatives of geometry objects 496 if isinstance(f, FacetNormal): 497 return swiginac.numeric(0.0) 498 499 if isinstance(f, SpatialCoordinate): 500 c, = c 501 if der[0] == c: 502 return swiginac.numeric(1.0) 503 else: 504 return swiginac.numeric(0.0) 505 506 sfc_error("Eh?")
507
508 - def derivative(self, x):
509 sfc_error("Derivative shouldn't occur here, you must apply AD first!")
510 511 ### Interior facet stuff: 512
513 - def positive_restricted(self, x, y):
514 sfc_error("TODO: Restrictions not implemented!") 515 return y
516
517 - def negative_restricted(self, x, y):
518 sfc_error("TODO: Restrictions not implemented!") 519 return y
520 521 522 ### These require code structure and thus shouldn't occur in SwiginacEvaluator 523 # (i.e. any conditionals should be handled externally) 524 #d[EQ] = 525 #d[NE] = 526 #d[LE] = 527 #d[GE] = 528 #d[LT] = 529 #d[GT] = 530 531 ### These are replaced by expand_compounds, so we skip them here: 532 #d[Identity] = 533 #d[Transposed] = 534 #d[Outer] = 535 #d[Inner] = 536 #d[Dot] = 537 #d[Cross] = 538 #d[Trace] = 539 #d[Determinant] = 540 #d[Inverse] = 541 #d[Deviatoric] = 542 #d[Cofactor] = 543 #d[Grad] = 544 #d[Div] = 545 #d[Curl] = 546 #d[Rot] = 547