Package sfc :: Package quadrature :: Module quadrature
[hide private]
[frames] | no frames]

Source Code for Module sfc.quadrature.quadrature

 1  #!/usr/bin/env python 
 2  # -*- coding: utf-8 -*- 
 3  """ 
 4  This module handles finding quadrature rules based on certain properties. 
 5  """ 
 6   
 7  # Copyright (C) 2008 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:  2008-08-12 
25  # Last changed: 2008-08-12 
26   
27  from sfc.quadrature.QuadRule import QuadRule 
28  from sfc.quadrature.quad_tables import load_rule 
29  from ufl.geometry import domain2dim 
30   
31 -def find_quadrature_rule(polygon, order, family="gauss"):
32 if polygon == "point": 33 return QuadRule(polygon, 1, 0, [(0.0,)], [1.0], "Single point evaluation rule.") 34 35 if family == "gauss": 36 if polygon in ("triangle", "tetrahedron"): 37 family = "economical_gauss" 38 else: 39 family = "composite_gauss" 40 41 if family == "economical_gauss": 42 return load_rule(polygon, order, family) 43 44 if family == "composite_gauss": 45 46 rule = load_rule("interval", order, "gauss") 47 #print "rule = ", rule 48 49 nsd = domain2dim[polygon] 50 if nsd == 1: 51 return rule 52 53 # construct composite rule from 1D formula 54 assert polygon in ("quadrilateral", "hexahedron") 55 points = [] 56 weights = [] 57 rr = range(rule.num_points) 58 if nsd == 2: 59 for i in rr: 60 for j in rr: 61 points.append( (rule.points[i][0], rule.points[j][0]) ) 62 weights.append( rule.weights[i]*rule.weights[j] ) 63 elif nsd == 3: 64 for i in rr: 65 for j in rr: 66 for k in rr: 67 points.append( (rule.points[i][0], rule.points[j][0], rule.points[k][0]) ) 68 weights.append( rule.weights[i]*rule.weights[j]*rule.weights[k] ) 69 comment = "Computed composite_gauss rule on %s of order %d." % (polygon, order) 70 return QuadRule(polygon, nsd, order, points, weights, comment) 71 72 raise RuntimeError("Found no rule for a %s with order %d in family %s" % (polygon, order, family))
73 74 75 if __name__ == '__main__': 76 print "No test here." 77