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

Source Code for Module sfc.quadrature.quadrature_code

  1  #!/usr/bin/env python 
  2  # -*- coding: utf-8 -*- 
  3  """ 
  4  This module contains code generation tools for quadrature rules. 
  5  """ 
  6   
  7  # Copyright (C) 2008-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:  2008-08-13 
 25  # Last changed: 2009-03-19 
 26   
 27  from sfc.codegeneration.codeformatting import indent 
 28   
29 -def gen_quadrature_rule_definition(rule, 30 points_name="quadrature_points", 31 weigths_name="quadrature_weights"):
32 code = "" 33 code += "static const double %s[%d][%d] = \n" % (points_name, rule.num_points, 34 rule.nsd) 35 code += " {\n" 36 code += ",\n".join( [ " { %s }" % (", ".join(map(str,p))) 37 for p in rule.points ] ) 38 code += "\n };\n" 39 code += "\n" 40 code += "static const double %s[%d] =\n" % (weigths_name, rule.num_points) 41 code += " {\n" 42 code += " " + ", ".join( map(str, rule.weights) ) 43 code += "\n };\n" 44 return code
45
46 -def gen_quadrature_rule_definitions(rules, 47 points_name="quadrature_points", 48 weigths_name="quadrature_weights"):
49 """Generate static code for quadrature rule definitions for a list of rules.""" 50 nr = len(rules) 51 nq = rules[0].num_points 52 nsd = rules[0].nsd 53 54 code = "" 55 56 rule_snippets = [] 57 for r in range(nr): 58 rule = rules[r] 59 60 point_snippets = [] 61 for p in rule.points: 62 c = "{ " + ", ".join(map(str,p)) + " }" 63 point_snippets.append( indent(c) ) 64 65 cc = "{ // quadrature points for rule %d:\n" % r 66 cc += ",\n".join(point_snippets) + "\n" 67 cc += "}" 68 rule_snippets.append( indent(cc) ) 69 70 code += "// [number of rules][number of points][space dimension]\n" 71 code += "static const double %s[%d][%d][%d] = \n" % (points_name, nr, nq, nsd) 72 code += indent("{") + "\n" 73 code += indent(",\n".join(rule_snippets)) + "\n" 74 code += indent("};") + "\n\n" 75 76 rule_snippets= [] 77 for r in range(nr): 78 rule = rules[r] 79 c = "// weights for rule %d:\n" % r 80 c += "{ " 81 c += ", ".join(map(str, rule.weights)) 82 c += " }" 83 rule_snippets.append( indent(c) ) 84 85 code += "// [number of rules][number of points]\n" 86 code += "static const double %s[%d][%d] =\n" % (weigths_name, nr, nq) 87 code += indent("{") + "\n" 88 code += indent(",\n".join(rule_snippets)) + "\n" 89 code += indent("};") + "\n\n" 90 91 return code
92 93 94 # TODO: Move to unit tests
95 -def test_quadrature_code():
96 inner_code = "// Compute A here\n" 97 98 order = 3 99 polygon = "triangle" 100 101 rule = find_quadrature_rule(polygon, order, "gauss") 102 103 code = gen_quadrature_rule_definition(rule) 104 print code 105 106 rules = [rule, rule, rule] 107 code = gen_quadrature_rule_definitions(rules) 108 print code
109 110 if __name__ == '__main__': 111 test_quadrature_code() 112