30 quadrature_option::family_type
f = opt.get_family();
34 if (
f == quadrature_option::equispaced) {
37 wx (
x(1/
T(3),1/
T(3),1/(3)),
T(1)/6);
40 T w = (1/
T(6))/
T(
int(n));
44 wx (
x(
T(
int(i))/r,
T(
int(j))/r,
T(
int(k))/r), w);
54 if (
f == quadrature_option::superconvergent) {
55 switch (opt.get_order()) {
59 f = quadrature_option::gauss;
62 error_macro (
"unsupported superconvergent("<<opt.get_order()<<
")");
69 if (
f == quadrature_option::gauss_lobatto) {
70 switch (opt.get_order()) {
74 wx(
x(
T(0),
T(0),
T(0)),
T(1)/24);
75 wx(
x(
T(1),
T(0),
T(0)),
T(1)/24);
76 wx(
x(
T(0),
T(1),
T(0)),
T(1)/24);
77 wx(
x(
T(0),
T(0),
T(1)),
T(1)/24);
83 error_macro (
"unsupported Gauss-Lobatto("<<opt.get_order()<<
")");
91 "unsupported quadrature family \"" << opt.get_family_name() <<
"\"");
93 switch (opt.get_order()) {
101 wx (
x(0.25, 0.25, 0.25),
T(1)/6);
107 T a = (5 - sqrt(
T(5)))/20;
108 T b = (5 + 3*sqrt(
T(5)))/20;
109 wx (
x(a, a, a),
T(1)/24);
110 wx (
x(a, a, b),
T(1)/24);
111 wx (
x(a, b, a),
T(1)/24);
112 wx (
x(b, a, a),
T(1)/24);
124 T b1 = (7 + sqrt(
T(15)))/34;
125 T b2 = (7 - sqrt(
T(15)))/34;
126 T c1 = (13 - 3*sqrt(
T(15)))/34;
127 T c2 = (13 + 3*sqrt(
T(15)))/34;
128 T w1 = (2665 - 14*sqrt(
T(15)))/226800;
129 T w2 = (2665 + 14*sqrt(
T(15)))/226800;
130 T d = (5 - sqrt(
T(15)))/20;
131 T e = (5 + sqrt(
T(15)))/20;
132 wx (
x(a, a, a),
T(8)/405);
133 wx (
x(b1, b1, b1), w1);
134 wx (
x(b1, b1, c1), w1);
135 wx (
x(b1, c1, b1), w1);
136 wx (
x(c1, b1, b1), w1);
137 wx (
x(b2, b2, b2), w2);
138 wx (
x(b2, b2, c2), w2);
139 wx (
x(b2, c2, b2), w2);
140 wx (
x(c2, b2, b2), w2);
141 wx (
x(
d,
d, e),
T(5)/567);
142 wx (
x(
d, e,
d),
T(5)/567);
143 wx (
x(e,
d,
d),
T(5)/567);
144 wx (
x(
d, e, e),
T(5)/567);
145 wx (
x(e,
d, e),
T(5)/567);
146 wx (
x(e, e,
d),
T(5)/567);
152 size_t r = opt.get_order();
156 vector<T> zeta0(n0), omega0(n0);
157 vector<T> zeta1(n1), omega1(n1);
158 vector<T> zeta2(n2), omega2(n2);
167 T eta_0 = (1+zeta0[i])*(1-zeta1[j])*(1-zeta2[k])/8;
168 T eta_1 = (1+zeta1[j])*(1-zeta2[k])/4;
169 T eta_2 = (1+zeta2[k])/2;
170 T J = (1-zeta1[j])*sqr(1-zeta2[k])/64;
171 wx (
x(eta_0,eta_1,eta_2), J*omega0[i]*omega1[j]*omega2[k]);
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
void gauss_jacobi(Size R, typename std::iterator_traits< OutputIterator1 >::value_type alpha, typename std::iterator_traits< OutputIterator1 >::value_type beta, OutputIterator1 zeta, OutputIterator2 omega)