1
2
3 """
4 This module handles finding quadrature rules based on certain properties.
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.quadrature.QuadRule import QuadRule
28 from sfc.quadrature.quad_tables import load_rule
29 from ufl.geometry import domain2dim
30
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
48
49 nsd = domain2dim[polygon]
50 if nsd == 1:
51 return rule
52
53
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