Package sfc :: Package symbolic_utils :: Module integration
[hide private]
[frames] | no frames]

Source Code for Module sfc.symbolic_utils.integration

  1  # -*- coding: utf-8 -*- 
  2  """ 
  3  This module contains functions for computing  
  4  integrals both symbolically and numerically. 
  5  """ 
  6   
  7  # Copyright (C) 2007-2009 Martin Sandve Alnes and Simula Resarch Laboratory 
  8  # 
  9  # This file is part of SyFi. 
 10  # 
 11  # SyFi is free software: you can redistribute it and/or modify 
 12  # it under the terms of the GNU General Public License as published by 
 13  # the Free Software Foundation, either version 2 of the License, or 
 14  # (at your option) any later version. 
 15  # 
 16  # SyFi is distributed in the hope that it will be useful, 
 17  # but WITHOUT ANY WARRANTY; without even the implied warranty of 
 18  # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 
 19  # GNU General Public License for more details. 
 20  # 
 21  # You should have received a copy of the GNU General Public License 
 22  # along with SyFi. If not, see <http://www.gnu.org/licenses/>. 
 23  # 
 24  # First added:  2007-11-07 
 25  # Last changed: 2009-03-19 
 26   
 27  import swiginac 
 28  import SyFi 
 29   
 30  from sfc.common import sfc_error 
 31  from sfc.quadrature import find_quadrature_rule 
 32  from sfc.geometry import UFCCell, geometry_mapping 
 33  from sfc.symbolic_utils.symbolic_utils import symbol, symbols 
 34   
35 -def symbolic_integral(polygon, function):
36 """ 37 This function computes the integral of the specified function over 38 the specified polygon symbolically. 39 The function may take as arguments the global coordinates (x,y,z) 40 as well as coordinates on the reference polygon (xi0, xi1, xi2). 41 """ 42 nsd = polygon.no_space_dim() 43 G, x0 = geometry_mapping(polygon) # FIXME: passing Rectangle (or Box) here will yield an affine mapping only correct for straight rectangles 44 45 # TODO: Get symbols as input 46 x, y, z = symbols(["x", "y", "z"]) 47 xi0, xi1, xi2 = symbols(["xi0", "xi1", "xi2"]) 48 p = swiginac.matrix(nsd, 1, [x,y,z]) 49 Ginv = G.inverse() 50 51 xi = Ginv*(p-x0) 52 ss = function.subs(xi0 == xi[0,0]) 53 if xi.rows() > 1: ss = ss.subs(xi1 == xi[1,0]) 54 if xi.rows() > 2: ss = ss.subs(xi2 == xi[2,0]) 55 56 return polygon.integrate(ss).evalf()
57
58 -def numeric_integral(polygon, function, order):
59 """ 60 This function computes the integral of the specified function over 61 the specified polygon numerically. 62 The function may take as arguments the global coordinates (x,y,z) 63 as well as coordinates on the reference polygon (xi0, xi1, xi2). 64 """ 65 66 # TODO: Get symbols as input 67 68 nsd = polygon.no_space_dim() 69 # UFCCell is only defined on reference domains ? 70 # OLD code 71 # cell = UFCCell(polygon) 72 # table = find_quadrature_rule(cell.shape, order, family="gauss") 73 shape = "" 74 if polygon.str().find("Triangle") >= 0: 75 shape = "triangle" 76 if polygon.str().find("Line") >= 0: 77 shape = "interval" 78 if polygon.str().find("Tetrahedron") >= 0: 79 shape = "tetrahedron" 80 table = find_quadrature_rule(shape, order, family="gauss") 81 82 x, y, z = symbols(["x", "y", "z"]) 83 xi0, xi1, xi2 = symbols(["xi0", "xi1", "xi2"]) 84 85 G, x0 = geometry_mapping(polygon) # FIXME: passing Rectangle (or Box) here will yield an affine mapping only correct for straight rectangles 86 xi = swiginac.matrix(nsd, 1, [xi0, xi1, xi2]) 87 global_point = G*xi + x0 88 ss = function.subs(x == global_point[0,0]) 89 if global_point.rows() > 1: ss = ss.subs(y == global_point[1,0]) 90 if global_point.rows() > 2: ss = ss.subs(z == global_point[2,0]) 91 92 sum = 0.0 93 for k in range(table.num_points): 94 p = table.points[k] 95 w = table.weights[k] 96 s = ss.subs(xi0 == p[0]) 97 if len(p) > 1: s = s.subs(xi1 == p[1]) 98 if len(p) > 2: s = s.subs(xi2 == p[2]) 99 sum += s*w 100 sum *= G.determinant() 101 return sum
102
103 -def integral(polygon, function, integration_mode=True, integration_order=None):
104 """ 105 A small wrapper on top of the functions 106 symbolic_integral and numeric_integral. 107 """ 108 # TODO: Support Taylor series expansion as well 109 if integration_mode == "symbolic": 110 return symbolic_integral(polygon, function) 111 elif integration_mode == "quadrature": 112 return numeric_integral(polygon, function, integration_order) 113 sfc_error("Not implemented integration mode %s" % integration_mode)
114