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, points_name="quad_points", weigths_name="quad_weights"):
30 code = "" 31 code += "static const double %s[%d][%d] = \n" % (points_name, rule.num_points, rule.nsd) 32 code += " {\n" 33 code += ",\n".join( [ " { %s }" % (", ".join(map(str,p))) for p in rule.points ] ) 34 code += "\n };\n" 35 code += "\n" 36 code += "static const double %s[%d] =\n" % (weigths_name, rule.num_points) 37 code += " {\n" 38 code += " " + ", ".join( map(str, rule.weights) ) 39 code += "\n };\n" 40 return code
41
42 -def gen_quadrature_rule_definitions(rules, points_name="facet_quad_points", weigths_name="facet_quad_weights"):
43 """Generate static code for quadrature rule definitions for a list of rules.""" 44 nr = len(rules) 45 nq = rules[0].num_points 46 nsd = rules[0].nsd 47 48 code = "" 49 50 rule_snippets = [] 51 for r in range(nr): 52 rule = rules[r] 53 54 point_snippets = [] 55 for p in rule.points: 56 c = "{ " + ", ".join(map(str,p)) + " }" 57 point_snippets.append( indent(c) ) 58 59 cc = "{ // quadrature points for rule %d:\n" % r 60 cc += ",\n".join(point_snippets) + "\n" 61 cc += "}" 62 rule_snippets.append( indent(cc) ) 63 64 code += "// [number of rules][number of points][space dimension]\n" 65 code += "static const double %s[%d][%d][%d] = \n" % (points_name, nr, nq, nsd) 66 code += indent("{") + "\n" 67 code += indent(",\n".join(rule_snippets)) + "\n" 68 code += indent("};") + "\n\n" 69 70 rule_snippets= [] 71 for r in range(nr): 72 rule = rules[r] 73 c = "// weights for rule %d:\n" % r 74 c += "{ " 75 c += ", ".join(map(str, rule.weights)) 76 c += " }" 77 rule_snippets.append( indent(c) ) 78 79 code += "// [number of rules][number of points]\n" 80 code += "static const double %s[%d][%d] =\n" % (weigths_name, nr, nq) 81 code += indent("{") + "\n" 82 code += indent(",\n".join(rule_snippets)) + "\n" 83 code += indent("};") + "\n\n" 84 85 return code
86 87 88 if __name__ == '__main__': 89 inner_code = "// Compute A here\n" 90 91 order = 3 92 polygon = "triangle" 93 94 rule = find_quadrature_rule(polygon, order, "gauss") 95 96 code = gen_quadrature_rule_definition(rule) 97 print code 98 99 rules = [rule, rule, rule] 100 code = gen_quadrature_rule_definitions(rules) 101 print code 102