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 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 shape = "" 70 if polygon.str().find("Triangle") >= 0: 71 shape = "triangle" 72 if polygon.str().find("Line") >= 0: 73 shape = "interval" 74 if polygon.str().find("Tetrahedron") >= 0: 75 shape = "tetrahedron" 76 table = find_quadrature_rule(shape, order, family="gauss") 77 78 x, y, z = symbols(["x", "y", "z"]) 79 xi0, xi1, xi2 = symbols(["xi0", "xi1", "xi2"]) 80 81 G, x0 = geometry_mapping(polygon) # FIXME: passing Rectangle (or Box) here will yield an affine mapping only correct for straight rectangles 82 xi = swiginac.matrix(nsd, 1, [xi0, xi1, xi2]) 83 global_point = G*xi + x0 84 ss = function.subs(x == global_point[0,0]) 85 if global_point.rows() > 1: ss = ss.subs(y == global_point[1,0]) 86 if global_point.rows() > 2: ss = ss.subs(z == global_point[2,0]) 87 88 sum = 0.0 89 for k in range(table.num_points): 90 p = table.points[k] 91 w = table.weights[k] 92 s = ss.subs(xi0 == p[0]) 93 if len(p) > 1: s = s.subs(xi1 == p[1]) 94 if len(p) > 2: s = s.subs(xi2 == p[2]) 95 sum += s*w 96 sum *= G.determinant() 97 return sum
98
99 -def integral(polygon, function, integration_mode=True, integration_order=None):
100 """ 101 A small wrapper on top of the functions 102 symbolic_integral and numeric_integral. 103 """ 104 # TODO: Support Taylor series expansion as well 105 if integration_mode == "symbolic": 106 return symbolic_integral(polygon, function) 107 elif integration_mode == "quadrature": 108 return numeric_integral(polygon, function, integration_order) 109 sfc_error("Not implemented integration mode %s" % integration_mode)
110