1
2
3 """
4 This module contains code generation tools for quadrature rules.
5 """
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27 from sfc.codegeneration.codeformatting import indent
28
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
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