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

Source Code for Module sfc.representation.integralrep_regular

  1  # -*- coding: utf-8 -*- 
  2  """ 
  3  This module contains representation classes for integrals. 
  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  # Modified by Kent-Andre Mardal, 2010. 
 24  # 
 25  # First added:  2008-08-13 
 26  # Last changed: 2009-03-12 
 27   
 28  from itertools import izip 
 29   
 30  #import SyFi 
 31  import swiginac 
 32   
 33  #import ufl 
 34  from ufl.permutation import compute_indices 
 35  from ufl.algorithms import (Graph, partition, expand_indices2, expand_indices, 
 36                              strip_variables, tree_format) 
 37  from ufl.classes import Expr, Terminal, UtilityType, SpatialDerivative, Indexed 
 38   
 39  from sfc.common import sfc_assert, sfc_warning, sfc_debug, sfc_error 
 40  from sfc.common.utilities import indices_subset 
 41  from sfc.symbolic_utils import symbols, symbol 
 42  from sfc.codegeneration.codeformatting import indent, row_major_index_string 
 43  from sfc.representation.swiginac_eval import SwiginacEvaluator, EvaluateAsSwiginac 
 44   
 45  from sfc.representation.integralrepresentationbase import IntegralRepresentationBase 
 46   
47 -def find_locals(Vin, Vout, P):
48 # Count the number of vertices depending on each vertex 49 Vcount = [len(vi) for vi in Vin] 50 # Count down all dependencies internally in the partition, O(|E| log|P|) 51 Pset = set(P) 52 for i in P: # v_i is in our partition 53 for j in Vout[i]: # v_i depends on v_j 54 if j in Pset: # is v_j in our partition? 55 Vcount[j] -= 1 # don't count with internal dependencies in partition 56 # Return mapping (v_i -> is_local_to_P) for all v_i in P 57 return dict((i, (Vcount[i] == 0)) for i in P)
58
59 -def is_simple(e):
60 return (e.nops() == 0)
61 # TODO: Try something like this: 62 #return (e.nops() == 0) or \ 63 # (isinstance(e, (swiginac.add, swiginac.mul)) \ 64 # and all(o.nops() for o in e.ops())) 65
66 -class IntegralData(object):
67 - def __init__(self, itgrep, integral):
68 self.itgrep = itgrep 69 self.integral = integral 70 71 # TODO: Override these and other options with metadata 72 metadata = integral.measure().metadata() 73 self.integration_method = itgrep.options.integration_method 74 self.integration_order = itgrep.options.integration_order 75 76 self.quad_rule = None 77 self.facet_quad_rules = None 78 79 self.G = None 80 self.Vdeps = None 81 self.Vin_count = None 82 self.partitions = {} 83 self.evaluator = None 84 self.evaluate = None 85 86 self.V_symbols = {} 87 self.vertex_data_set = {}
88
89 - def store(self, value, i, basis_functions):
90 #print "STORING AT", i, basis_functions, str(value) 91 vd = self.vertex_data_set.get(basis_functions) 92 if vd is None: 93 vd = {} 94 self.vertex_data_set[basis_functions] = vd 95 vd[i] = value
96
97 - def get_vd(self, j, basis_functions):
98 jkeep = [(("v%d" % k) in self.Vdeps[j]) 99 for k in range(self.itgrep.formrep.rank)] 100 iota, = indices_subset([basis_functions], jkeep) 101 if all(i is None for i in iota): 102 iota = () 103 vd = self.vertex_data_set[iota] 104 return vd
105
106 - def fetch_storage(self, j, basis_functions):
107 vd = self.get_vd(j, basis_functions) 108 return vd[j]
109
110 - def free_storage(self, j, basis_functions):
111 "Delete stored data for j." 112 vd = self.get_vd(j, basis_functions)
113 #del vd[j] # FIXME: Enable when counts are corrected 114
115 -class IntegralRepresentation(IntegralRepresentationBase):
116 - def __init__(self, integrals, formrep, on_facet):
117 IntegralRepresentationBase.__init__(self, integrals, formrep, on_facet) 118 sfc_debug("Entering IntegralRepresentation.__init__") 119 120 # Shape of element tensor 121 A_shape = [] 122 for i in range(self.formrep.rank): 123 element = self.formrep.formdata.elements[i] 124 rep = self.formrep.element_reps[element] 125 A_shape.append(rep.local_dimension) 126 self.A_shape = tuple(A_shape) 127 128 # All permutations of element tensor indices 129 # TODO: Not same for interior_facet_integrals, move to subclasses 130 self.indices = compute_indices(self.A_shape) 131 132 # Symbols for output tensor entry A[iota] 133 Asym = symbols("A[%s]" % row_major_index_string(i, self.A_shape) 134 for i in self.indices) 135 self.A_sym = dict(izip(self.indices, Asym)) 136 137 # Data structures to hold integrals for each measure 138 # TODO: Consider integral metadata here to mix and match integration modes 139 if self.options.integration_method == "symbolic": 140 # Only a single symbolic integral is allowed: 141 assert len(self.integrals) == 1 142 self.symbolic_integral = self.integrals[0] 143 self.quadrature_integrals = [] 144 elif self.options.integration_method == "quadrature": 145 # Multiple quadrature integrals are possible in principle, each with 146 # different quadrature rules, but this is not fully implemented: 147 self.symbolic_integral = None 148 self.quadrature_integrals = self.integrals 149 150 # Data structures for token symbols 151 self._free_symbols = [] 152 self._symbol_counter = 0 153 self._build_integral_data() 154 155 sfc_debug("Leaving IntegralRepresentation.__init__")
156
157 - def _build_integral_data(self):
158 # Mapping from measure to IntegralData object 159 self.integral_data = {} 160 161 # Compute data about integrals for each measure 162 fd = self.formrep.formdata 163 for integral in self.integrals: 164 # TODO: Split IntegralData classes into one for each integration mode 165 if self.options.integration_method == "symbolic": 166 data = self._compute_symbolic_integral_data(integral) 167 elif self.options.integration_method == "quadrature": 168 data = self._compute_quadrature_integral_data(integral) 169 self.integral_data[integral.measure()] = data
170
171 - def _compute_symbolic_integral_data(self, integral):
172 data = IntegralData(self, integral) 173 data.evaluator = SwiginacEvaluator(self.formrep, use_symbols=False, 174 on_facet=self._on_facet) 175 return data
176
177 - def _compute_quadrature_integral_data(self, integral):
178 data = IntegralData(self, integral) 179 180 # FIXME: Build a separate quadrature rule for each integral, if necessary 181 data.quad_rule = self.formrep.quad_rule 182 data.facet_quad_rules = self.formrep.facet_quad_rules 183 184 if self.options.safemode: 185 # Produce naive code, good for debugging and verification 186 data.evaluator = SwiginacEvaluator(self.formrep, use_symbols=True, 187 on_facet=self._on_facet) 188 189 else: 190 # Create an evaluation multifunction 191 data.evaluate = EvaluateAsSwiginac(self.formrep, self, data, 192 on_facet=self._on_facet) 193 194 # Expand all indices before building graph 195 integrand = integral.integrand() 196 integrand = strip_variables(integrand) 197 if self.options.use_expand_indices2: 198 integrand = expand_indices2(integrand) 199 else: 200 integrand = expand_indices(integrand) 201 202 # Build linearized computational graph 203 data.G = Graph(integrand) 204 205 # Count dependencies for each node 206 V = data.G.V() 207 n = len(V) 208 data.Vin_count = [0]*n 209 for i, vs in enumerate(data.G.Vin()): 210 for j in vs: 211 vj = V[j] 212 # ... FIXME: Need to store one counter for each iota relevant for i, and to add "size(iota_space_j)" instead of 1 213 data.Vin_count[i] += 1 214 #data.Vin_count[-1] += ... # FIXME: Increment further for integrand root vertex 215 216 # Partition graph based on dependencies 217 #basis_function_deps = [set(("v%d" % i, "x")) \ 218 # for (i, bf) in enumerate(fd.basis_functions)] 219 #function_deps = [set(("w", "c", "x",)) \ 220 # for (i, bf) in enumerate(fd.functions)] 221 # TODO: Input actual dependencies of functions and 222 # basis functions and their derivatives here 223 data.partitions, data.Vdeps = partition(data.G) 224 return data
225 226 # --- Printing --- 227
228 - def __str__(self):
229 s = "" 230 s += "IntegralRepresentation:\n" 231 s += " classname: %s\n" % self.classname 232 s += " A_shape: %s\n" % self.A_shape 233 s += " indices: %s\n" % self.indices 234 s += " A_sym: %s\n" % self.A_sym 235 236 # Not computing whole A at any point because of memory usage 237 #s += " A:\n" 238 #for ii in self.indices: 239 # s += " %s = %s\n" % (self.A_sym[ii], self.A[ii]) 240 241 s += " UFL Integrals:\n" 242 for integral in self.integrals: 243 s += indent(str(integral)) + "\n\n" 244 return s.strip("\n")
245 246 # --- Utilities for allocation of symbols and association of data with graph vertices --- 247
248 - def allocate(self, data, j):
249 """Allocate symbol(s) for vertex j. Either gets symbol from the free 250 symbol set or creates a new symbol and increases the symbol counter.""" 251 if self._free_symbols: 252 return self._free_symbols.pop() 253 s = symbol("s[%d]" % self._symbol_counter) 254 self._symbol_counter += 1 255 data.V_symbols[j] = s 256 return s
257
258 - def free_symbols(self, data, j):
259 "Delete stored data for j and make its symbols available again." 260 s = data.V_symbols.get(j) 261 if s is None: 262 sfc_debug("Trying to deallocate symbols that are not allocated!") 263 return 264 self._free_symbols.append(s) 265 del data.V_symbols[j]
266 267 # --- Generators for tokens of various kinds --- 268
269 - def iter_partition(self, data, deps, basis_functions=()):
270 sfc_debug("Entering IntegralRepresentation.iter_partition") 271 272 deps = frozenset(deps) 273 274 P = data.partitions.get(deps) 275 if not P: 276 sfc_debug("Leaving IntegralRepresentation.iter_partition, empty") 277 return 278 279 data.evaluate.current_basis_function = basis_functions 280 281 # Get graph components 282 G = data.G 283 V = G.V() 284 E = G.E() 285 Vin = G.Vin() 286 Vout = G.Vout() 287 288 # Figure out which variables are only needed within the partition 289 is_local = find_locals(Vin, Vout, P) 290 291 # For all vertices in partition P 292 for i in P: 293 v = V[i] 294 295 if isinstance(v, UtilityType): 296 # Skip labels and indices 297 continue 298 299 if v.shape(): 300 # Skip expressions with shape: they will be evaluated when indexed. 301 if not isinstance(v, (Terminal, SpatialDerivative)): 302 print "="*30 303 print "type:", type(v) 304 print "str:", str(v) 305 print "child types:", [type(V[j]) for j in Vout[i]] 306 print "child str:" 307 print "\n".join( " vertex %d: %s" % (j, str(V[j])) for j in Vout[i] ) 308 print "number of parents:", len(Vin[i]) 309 if len(Vin[i]) < 5: 310 print "parent types:", [type(V[j]) for j in Vin[i]] 311 #print "parent str:" 312 #print "\n".join( " vertex %d: %s" % (j, str(V[j])) for j in Vin[i] ) 313 sfc_error("Expecting all indexing to have been propagated to terminals?") 314 continue 315 316 if Vin[i] and all(isinstance(V[j], SpatialDerivative) for j in Vin[i]): 317 # Skipping expressions that only occur later in 318 # differentiated form: we don't need their non-differentiated 319 # expression so they will be evaluated when indexed. 320 if not isinstance(v, (Terminal, SpatialDerivative, Indexed)): 321 print "="*30 322 print type(v) 323 print str(v) 324 sfc_error("Expecting all indexing to have been propagated to terminals?") 325 continue 326 327 if isinstance(v, (Indexed, SpatialDerivative)): 328 ops = v.operands() 329 330 # Evaluate vertex v with given mapped ops 331 if not all(isinstance(o, (Expr, swiginac.basic)) for o in ops): 332 print ";"*80 333 print tree_format(v) 334 print str(v) 335 print type(ops) 336 print str(ops) 337 print repr(ops) 338 print "types:" 339 print "\n".join(str(type(o)) for o in ops) 340 print ";"*80 341 342 e = data.evaluate(v, *ops) 343 344 else: 345 # Get already computed operands 346 # (if a vertex isn't already computed 347 # at this point, that is a bug) 348 ops = [] 349 for j in Vout[i]: 350 try: 351 # Fetch expression or symbol for vertex j 352 e = data.fetch_storage(j, basis_functions) 353 except: 354 print "Failed to fetch expression for vertex %d," % j 355 print " V[%d] = %s" % (j, repr(V[j])) 356 print " parent V[%d] = %s" % (i, repr(V[i])) 357 raise RuntimeError 358 ops.append(e) 359 ops = tuple(ops) 360 e = data.evaluate(v, *ops) 361 362 # TODO: Make sure that dependencies of skipped expressions 363 # are counted down when evaluated later! 364 # Or do we need that? We don't make symbols for them anyway. 365 366 # Count down number of uses of v's dependencies 367 for j in Vout[i]: # For each vertex j that depend on vertex i 368 data.Vin_count[j] -= 1 369 if False: # data.Vin_count[j] == 0: # FIXME: This doesn't work yet, counters are wrong 370 # Make the symbols for j available again 371 self.free_symbols(data, j) 372 data.free_storage(j, basis_functions) 373 374 # Store some result for vertex i and yield token based on some heuristic 375 if is_simple(e): 376 # Store simple expression associated with 377 # vertex j, no code generation necessary 378 data.store(e, i, basis_functions) 379 else: 380 if is_local[i]: 381 pass # TODO: Use this information to make a local set of symbols "sl[%d]" for each partition 382 383 # Allocate a symbol and remember it, just throw away the expression 384 s = self.allocate(data, i) 385 data.store(s, i, basis_functions) 386 # Only yield when we have a variable to generate code for! 387 yield (s, e) 388 389 sfc_debug("Leaving IntegralRepresentation.iter_partition")
390 391 # Some stuff from the old code: 392 #bf = self._current_basis_function # FIXME: Filter depending on deps! 393 #key = tuple(chain((count,), component, index_values, bf)) 394 #compstr = "_".join("%d" % k for k in key) 395 #vname = "s_%s" % compstr 396 397 # Register token (s, e) with variable v in evaluator, 398 # such that it can return it from evaluator.variable 399 #self.evaluator._variable2symbol[key] = vsym 400
401 - def iter_member_quad_tokens(self, data):
402 """Return an iterator over member tokens dependent of spatial variables. 403 404 Overload in subclasses!""" 405 # FIXME: Make this into precomputation of a static array of constants 406 407 assert data.integration_method == "quadrature" 408 409 # TODO: Precompute basis function values symbolically here 410 # instead of generating a loop in the constructor 411 412 # TODO: Skip what's not needed! 413 414 # Precompute all basis functions 415 fr = self.formrep 416 fd = fr.formdata 417 generated = set() 418 for iarg in range(fr.rank + fr.num_coefficients): 419 elm = fd.elements[iarg] 420 rep = fr.element_reps[elm] 421 for i in range(rep.local_dimension): 422 for component in rep.value_components: 423 # Yield basis function itself 424 s = fr.v_sym(iarg, i, component, self._on_facet) 425 if not (s == 0 or s in generated): 426 e = fr.v_expr(iarg, i, component) 427 t = (s, e) 428 yield t 429 generated.add(s) 430 # Yield its derivatives w.r.t. local coordinates 431 for d in range(fr.cell.nsd): 432 der = (d,) 433 s = fr.dv_sym(iarg, i, component, der, self._on_facet) 434 if not (s == 0 or s in generated): 435 e = fr.dv_expr(iarg, i, component, der) 436 t = (s, e) 437 yield t 438 generated.add(s)
439
440 - def iter_geometry_tokens(self):
441 """Return an iterator over geometry tokens independent of spatial variables.""" 442 443 fr = self.formrep 444 445 # TODO: Skip what's not needed! 446 447 # vx 448 for (ss,ee) in zip(fr.vx_sym, fr.vx_expr): 449 for i in range(ss.nops()): 450 yield (ss.op(i), ee.op(i)) 451 452 # x0 is just == vx0 453 #(ss,ee) = (fr.x0_sym, fr.x0_expr) 454 #for i in range(ss.nops()): 455 # yield (ss.op(i), ee.op(i)) 456 457 # G 458 (ss,ee) = (fr.G_sym, fr.G_expr) 459 for i in range(ss.nops()): 460 yield (ss.op(i), ee.op(i)) 461 462 # detG 463 yield (fr.detGtmp_sym, fr.detGtmp_expr) 464 yield (fr.detG_sym, fr.detG_expr) 465 466 # Ginv 467 (ss,ee) = (fr.Ginv_sym, fr.Ginv_expr) 468 for i in range(ss.nops()): 469 yield (ss.op(i), ee.op(i)) 470 471 if self._on_facet: 472 # Needed for normal vector TODO: Skip if not needed 473 yield (fr.detG_sign_sym, fr.detG_sign_expr) 474 else: 475 if self.symbolic_integral is not None: 476 # Scaling by cell volume factor, determinant of coordinate mapping 477 yield (fr.D_sym, fr.detG_sym)
478
479 - def iter_runtime_quad_tokens(self, data):
480 "Return an iterator over runtime tokens dependent of spatial variables. Overload in subclasses!" 481 assert data.integration_method == "quadrature" 482 483 fr = self.formrep 484 gr = fr.geomrep 485 fd = fr.formdata 486 nsd = gr.sfc_cell.nsd 487 488 # TODO: yield geometry tokens like G etc here instead if they depend on x,y,z 489 490 # TODO: Skip what's not needed! 491 492 # global coordinates 493 (ss,ee) = (gr.x_sym, gr.x_expr) 494 for i in range(ss.nops()): 495 yield (ss.op(i), ee.op(i)) 496 497 # Generate all basis function derivatives 498 generated = set() 499 for iarg in range(fr.rank + fr.num_coefficients): 500 elm = fd.elements[iarg] 501 rep = fr.element_reps[elm] 502 for i in range(rep.local_dimension): 503 for component in rep.value_components: 504 # Yield first order derivatives w.r.t. global coordinates 505 for d in range(fr.cell.nsd): 506 der = (d,) 507 s = fr.Dv_sym(iarg, i, component, der, self._on_facet) 508 if not (s == 0 or s in generated): 509 e = fr.Dv_expr(iarg, i, component, der, True, self._on_facet) 510 t = (s, e) 511 yield t 512 generated.add(s) 513 514 # Currently placing all w in here: 515 generated = set() 516 for iarg in range(fr.num_coefficients): 517 elm = fd.elements[fr.rank+iarg] 518 rep = fr.element_reps[elm] 519 for component in rep.value_components: 520 # Yield coefficient function itself 521 s = fr.w_sym(iarg, component) 522 if not s in generated: 523 e = fr.w_expr(iarg, component, True, self._on_facet) 524 t = (s, e) 525 yield t 526 generated.add(s) 527 # Yield first order derivatives w.r.t. global coordinates 528 for d in range(fr.cell.nsd): 529 der = (d,) 530 s = fr.Dw_sym(iarg, component, der) 531 if not (s == 0 or s in generated): 532 e = fr.Dw_expr(iarg, component, der, True, self._on_facet) 533 t = (s, e) 534 yield t 535 generated.add(s) 536 537 # Scaling factor 538 if self._on_facet: 539 # TODO: yield tokens that depend on both x and facet here 540 541 # Scaling by facet area factor 542 D_expr = fr.quad_weight_sym*fr.facet_D_sym 543 else: 544 # Scaling by cell volume factor, determinant of coordinate mapping 545 D_expr = fr.quad_weight_sym*fr.detG_sym 546 yield (fr.D_sym, D_expr)
547 548 # --- Element tensor producers 549
550 - def iter_A_tokens(self, data, facet=None):
551 "Iterate over all A[iota] tokens." 552 for iota in self.indices: 553 A_sym = self.A_sym[iota] 554 A_expr = self.compute_A(data, iota, facet) 555 yield (A_sym, A_expr)
556
557 - def compute_A(self, data, iota, facet=None):
558 "Compute expression for A[iota]. Overload in subclasses!" 559 raise NotImplementedError
560 561
562 -class QuadratureCellIntegralRepresentation(IntegralRepresentation):
563 - def __init__(self, integrals, formrep):
564 IntegralRepresentation.__init__(self, integrals, formrep, False)
565
566 - def compute_A(self, data, iota, facet=None):
567 "Compute expression for A[iota]. Overload in subclasses!" 568 if self.options.safemode: 569 integrand = data.integral.integrand() 570 data.evaluator.update(iota) 571 integrand = data.evaluator.visit(integrand) 572 else: 573 n = len(data.G.V()) 574 integrand = data.vertex_data_set[iota][n-1] 575 576 D = self.formrep.D_sym 577 A = integrand * D 578 579 if self.formrep.options.output.enable_debug_prints: 580 sfc_debug("In compute_A(", iota, "):") 581 sfc_debug(" data.integral.integrand() = ", data.integral.integrand()) 582 sfc_debug(" integrand = ", integrand) 583 sfc_debug(" A = ", A) 584 585 return A
586
587 -class SymbolicCellIntegralRepresentation(IntegralRepresentation):
588 - def __init__(self, integrals, formrep):
589 IntegralRepresentation.__init__(self, integrals, formrep, False)
590
591 - def compute_A(self, data, iota, facet=None):
592 "Compute expression for A[iota]. Overload in subclasses!" 593 594 integrand = data.integral.integrand() 595 data.evaluator.update(iota) 596 integrand = data.evaluator.visit(integrand) 597 detG = self.formrep.detG_sym 598 polygon = self.formrep.cell.polygon 599 A = polygon.integrate(integrand) * detG 600 601 if self.formrep.options.output.enable_debug_prints: 602 sfc_debug("In compute_A(", iota, "):") 603 sfc_debug(" data.integral.integrand() = ", data.integral.integrand()) 604 sfc_debug(" integrand = ", integrand) 605 sfc_debug(" A = ", A) 606 607 return A
608 609
610 -class ExteriorFacetIntegralRepresentationBase(IntegralRepresentation):
611 - def __init__(self, integrals, formrep):
612 IntegralRepresentation.__init__(self, integrals, formrep, True)
613
614 - def iter_facet_tokens(self, facet):
615 fr = self.formrep 616 gr = fr.geomrep 617 nsd = fr.cell.nsd 618 619 # Facet mapping FIXME (needed anywere? quadrature on facets?) 620 #for (s,e) in zip(fr.facet_G_sym, fr.facet_G_expr[facet]): 621 # yield (s,e) 622 623 s = gr.facet_D_sym 624 e = gr.facet_D_expr[facet] 625 yield (s, e) 626 627 # Facet normal 628 for i in range(nsd): 629 yield (gr.n_tmp_sym[i], gr.n_tmp_expr[facet][i]) 630 yield (gr.nnorm_sym, gr.nnorm_expr) 631 for i in range(nsd): 632 yield (gr.n_sym[i], gr.n_expr[i]) 633 634 # Scaling factor 635 if self.symbolic_integral is not None: 636 yield (gr.D_sym, gr.facet_D_sym)
637 638
639 -class QuadratureExteriorFacetIntegralRepresentation(ExteriorFacetIntegralRepresentationBase):
640 - def __init__(self, integrals, formrep):
641 ExteriorFacetIntegralRepresentationBase.__init__(self, integrals, formrep)
642
643 - def compute_A(self, data, iota, facet=None):
644 "Compute expression for A[iota]. Overload in subclasses!" 645 646 integrand = data.integral.integrand() 647 648 sfc_assert(facet is None, "Not expecting facet number.") 649 650 if self.options.safemode: 651 data.evaluator.update(iota) 652 integrand = data.evaluator.visit(integrand) 653 else: 654 n = len(data.G.V()) 655 try: 656 integrand = data.vertex_data_set[iota][n-1] 657 except: 658 print data.vertex_data_set 659 raise RuntimeError 660 661 D = self.formrep.D_sym 662 A = integrand * D 663 664 if self.formrep.options.output.enable_debug_prints: 665 sfc_debug("In compute_A(", iota, "):") 666 sfc_debug(" data.integral.integrand() = ", data.integral.integrand()) 667 sfc_debug(" integrand = ", integrand) 668 sfc_debug(" A = ", A) 669 670 return A
671 672
673 -class SymbolicExteriorFacetIntegralRepresentation(ExteriorFacetIntegralRepresentationBase):
674 - def __init__(self, integrals, formrep):
675 ExteriorFacetIntegralRepresentationBase.__init__(self, integrals, formrep)
676
677 - def compute_A(self, data, iota, facet=None):
678 "Compute expression for A[iota]. Overload in subclasses!" 679 680 integrand = data.integral.integrand() 681 682 sfc_assert(isinstance(facet, int), "Expecting facet number.") 683 684 data.evaluator.update(iota) 685 integrand = data.evaluator.visit(integrand) 686 687 D = self.formrep.D_sym 688 polygon = self.formrep.cell.facet_polygons[facet] 689 if isinstance(polygon, swiginac.matrix): 690 # 1D 691 repmap = swiginac.exmap() 692 repmap[self.formrep.xi_sym[0]] = polygon[0] 693 A = integrand.subs(repmap) * D 694 else: 695 A = polygon.integrate(integrand) * D 696 697 if self.formrep.options.output.enable_debug_prints: 698 sfc_debug("In compute_A(", iota, "):") 699 sfc_debug(" data.integral.integrand() = ", data.integral.integrand()) 700 sfc_debug(" integrand = ", integrand) 701 sfc_debug(" A = ", A) 702 703 return A
704 705
706 -class InteriorFacetIntegralRepresentation(IntegralRepresentation):
707 - def __init__(self, integrals, formrep):
708 IntegralRepresentation.__init__(self, integrals, formrep, True)
709 710
711 -def CellIntegralRepresentation(integrals, formrep):
712 "Integral representation selection switch." 713 m = formrep.options.code.integral.integration_method 714 if m == "quadrature": 715 return QuadratureCellIntegralRepresentation(integrals, formrep) 716 elif m == "symbolic": 717 return SymbolicCellIntegralRepresentation(integrals, formrep)
718
719 -def ExteriorFacetIntegralRepresentation(integrals, formrep):
720 "Integral representation selection switch." 721 m = formrep.options.code.integral.integration_method 722 if m == "quadrature": 723 return QuadratureExteriorFacetIntegralRepresentation(integrals, formrep) 724 elif m == "symbolic": 725 return SymbolicExteriorFacetIntegralRepresentation(integrals, formrep)
726