1 """This module defines evaluation algorithms for converting
2 converting UFL expressions to swiginac representation."""
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
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
41 - def __init__(self, formrep, itgrep, data, on_facet):
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
50
51 - def expr(self, o, *ops):
52 sfc_error("Evaluation not implemented for expr %s." % type(o).__name__)
53
55 sfc_error("Evaluation not implemented for terminal %s." % type(o).__name__)
56
57
58
60 return swiginac.numeric(0)
61
63 return swiginac.numeric(o.value())
64
66
67
68 if component:
69
70 c, = component
71 else:
72
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
86
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
94 c, = component
95 else:
96
97 c = 0
98
99 return self.formrep.n_sym[c]
100
101 - def argument(self, o, component=(), derivatives=()):
102
103
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)
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
119 return s
120
121 - def coefficient(self, o, component=(), derivatives=()):
122
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)
130 else:
131
132 s = self.formrep.w_sym(iarg, component)
133 e = self.formrep.w_expr(iarg, component, False, self.on_facet)
134
135 if e.nops() == 0:
136 return e
137 return s
138
139
140
142 return tuple(map(int, o))
143
154
167
168
169
170 - def power(self, o, a, b):
172
173 - def sum(self, o, *ops):
175
178
181
182 - def abs(self, o, a):
183 return swiginac.abs(a)
184
185
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
218 sfc_error("Should strip away variables before building graph.")
219
220
221
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)
226
227
228 self._formrep = formrep
229 self._use_symbols = use_symbols
230 self._on_facet = on_facet
231
232
233 self._current_basis_function = tuple(0 for i in range(formrep.rank))
234
235
236 self._components = Stack()
237 self._index2value = StackDict()
238
239
240
241 self._variable2symbol = {}
242 self._tokens = []
243
244
245 self.nsd = self._formrep.cell.nsd
246
248
249 t = self._tokens
250 self._tokens = []
251 return t
252
254 self._current_basis_function = tuple(iota)
255
257 "Return current component tuple."
258 if len(self._components):
259 return self._components.peek()
260 return ()
261
262
263
265 sfc_error("Missing ufl to swiginac handler for type %s" % str(type(x)))
266
268 sfc_error("Missing ufl to swiginac handler for terminal type %s" % str(type(x)))
269
270
271
273
274 return swiginac.numeric(0)
275
277
278 return swiginac.numeric(x._value)
279
281 c = self.component()
282 v = 1 if c[0] == c[1] else 0
283 return swiginac.numeric(v)
284
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
296 iarg = x.count()
297 c = self.component()
298 if self._use_symbols:
299 return self._formrep.w_sym(iarg, c)
300 else:
301
302 return self._formrep.w_expr(iarg, c, False, self._on_facet)
303
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
310 c, = self.component()
311 return self._formrep.x_sym[c]
312
313
314
316 return self.visit(x._expression)
317
319 c = self.component()
320 index_values = tuple(self._index2value[k] for k in x._expression.free_indices())
321
322
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
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
336
337 - def sum(self, x, *ops):
339
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
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
358
359 - def power(self, x, a, b):
361
362 - def abs(self, x, a):
363 return swiginac.abs(a)
364
365
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
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
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
399
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
407 self._components.push(())
408
409
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
417 r = self.visit(e)
418
419
420 self._components.pop()
421 return r
422
424
425 c = self.component()
426 c0, c1 = c[0], c[1:]
427 op = x.operands()[c0]
428
429 self._components.push(c1)
430 r = self.visit(op)
431 self._components.pop()
432 return r
433
435
436
437 expression, indices = x.operands()
438 sfc_assert(expression.shape() == (), "Expecting scalar base expression.")
439
440
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
448 result = self.visit(expression)
449
450
451 for i in range(len(comp)):
452 self._index2value.pop()
453 self._components.pop()
454 return result
455
456
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
466
467
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!")
473
474
475 c = self.component()
476 der = self.visit(ii)
477
478
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
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
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
509 sfc_error("Derivative shouldn't occur here, you must apply AD first!")
510
511
512
514 sfc_error("TODO: Restrictions not implemented!")
515 return y
516
518 sfc_error("TODO: Restrictions not implemented!")
519 return y
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547