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

Source Code for Module sfc.representation.exteriorfacetintegralrepresentation

  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  # First added:  2008-08-13 
 24  # Last changed: 2009-03-16 
 25   
 26  import swiginac 
 27  from sfc.common import sfc_assert, sfc_debug 
 28  from sfc.representation.integralrepresentation import IntegralRepresentation 
 29   
30 -class ExteriorFacetIntegralRepresentation(IntegralRepresentation):
31 - def __init__(self, integrals, formrep):
32 IntegralRepresentation.__init__(self, integrals, formrep, True)
33
34 - def iter_facet_tokens(self, facet):
35 fr = self.formrep 36 nsd = fr.cell.nsd 37 38 # Facet mapping FIXME (needed anywere? quadrature on facets?) 39 #for (s,e) in zip(fr.facet_G_sym, fr.facet_G_expr[facet]): 40 # yield (s,e) 41 42 s = fr.facet_D_sym 43 e = fr.facet_D_expr[facet] 44 yield (s, e) 45 46 # Facet normal 47 for i in range(nsd): 48 yield (fr.n_sym[i], fr.n_expr[facet][i]) 49 50 # Scaling factor 51 if self.symbolic_integral is not None: 52 yield (fr.D_sym, fr.facet_D_sym)
53
54 - def compute_A(self, data, iota, facet=None):
55 "Compute expression for A[iota]. Overload in subclasses!" 56 57 integrand = data.integral.integrand() 58 59 if data.integration_method == "quadrature": 60 sfc_assert(facet is None, "Not expecting facet number.") 61 62 if self.options.safemode: 63 data.evaluator.update(iota) 64 integrand = data.evaluator.visit(integrand) 65 else: 66 n = len(data.G.V()) 67 try: 68 integrand = data.vertex_data_set[iota][n-1] 69 except: 70 print data.vertex_data_set 71 raise RuntimeError 72 73 D = self.formrep.D_sym 74 A = integrand * D 75 76 if self.formrep.options.output.enable_debug_prints: 77 sfc_debug("In compute_A(", iota, "):") 78 sfc_debug(" data.integral.integrand() = ", data.integral.integrand()) 79 sfc_debug(" integrand = ", integrand) 80 sfc_debug(" A = ", A) 81 82 elif data.integration_method == "symbolic": 83 sfc_assert(isinstance(facet, int), "Expecting facet number.") 84 85 data.evaluator.update(iota) 86 integrand = data.evaluator.visit(integrand) 87 88 D = self.formrep.D_sym 89 polygon = self.formrep.cell.facet_polygons[facet] 90 if isinstance(polygon, swiginac.matrix): 91 # 1D 92 repmap = swiginac.exmap() 93 repmap[self.formrep.xi_sym[0]] = polygon[0] 94 A = integrand.subs(repmap) * D 95 else: 96 A = polygon.integrate(integrand) * D 97 98 if self.formrep.options.output.enable_debug_prints: 99 sfc_debug("In compute_A(", iota, "):") 100 sfc_debug(" data.integral.integrand() = ", data.integral.integrand()) 101 sfc_debug(" integrand = ", integrand) 102 sfc_debug(" A = ", A) 103 104 return A
105