00001 #include "SundanceFeketeTriangleQuadrature.hpp"
00002 #include "PlayaExceptions.hpp"
00003 #include "SundanceOut.hpp"
00004 #include "SundancePoint.hpp"
00005 #include "PlayaTabs.hpp"
00006
00007 using namespace Sundance;
00008 using namespace Teuchos;
00009
00010 extern "C"
00011 {
00012
00013 void dgemv_(char *trans, int *m, int *n, double *alpha, double *a, int *lda,
00014 double *x, int *incx, double *beta, double *y, int *incy);
00015
00016
00017 void dgetrf_(const int* M, const int* N, double* A, const int* lda,
00018 const int* iPiv, int* info);
00019
00020
00021 void dgetri_(const int* n, double* a, const int* lda, const int* iPiv,
00022 double* work, const int* lwork, int* info);
00023 }
00024
00025
00026
00027
00028
00029
00030 void FeketeTriangleQuadrature::getPoints(int order, Array<double>& wgt, Array<
00031 double>& x, Array<double>& y)
00032 {
00033 Array<double> w;
00034 Array<int> multiplicity;
00035 Array<Array<double> > q;
00036
00037 if (order == 1)
00038 {
00039
00040 multiplicity = tuple(3);
00041 q.resize(1);
00042 w = tuple(1.0 / 3.0);
00043 q[0] = tuple(1.0, 0.0, 0.0);
00044 }
00045 else if (order == 2)
00046 {
00047
00048 multiplicity = tuple(3, 3);
00049 q.resize(2);
00050 w = tuple(0.0, 1.0 / 3.0);
00051 q[0] = tuple(1.0, 0.0, 0.0);
00052 q[1] = tuple(0.0, 0.5, 0.5);
00053 }
00054 else if (order == 3)
00055 {
00056 multiplicity = tuple(1, 3, 6);
00057 q.resize(3);
00058 w = tuple(0.45, 1.0 / 60.0, 1.0 / 12.0);
00059 q[0] = tuple(1.0 / 3.0, 1.0 / 3.0, 1.0 / 3.0);
00060 q[1] = tuple(1.0, 0.0, 0.0);
00061 q[2] = tuple(0.000000000000000, 0.276393202250021, 0.723606797749979);
00062 }
00063 else if (order == 4)
00064 {
00065 multiplicity = tuple(3, 3, 3, 6);
00066 q.resize(4);
00067 w = tuple(-0.002443433685378, 0.040684938686721, 0.200360939516545,
00068 0.047365444407723);
00069 q[0] = tuple(1.0, 0.0, 0.0);
00070 q[1] = tuple(0.0, 0.5, 0.5);
00071 q[2] = tuple(0.551583507555305, 0.224208246222347, 0.224208246222347);
00072 q[3] = tuple(0.000000000000000, 0.172673164646012, 0.827326835353988);
00073 }
00074 else if (order == 5)
00075 {
00076 multiplicity = tuple(3, 3, 3, 6, 6);
00077 q.resize(5);
00078 w = tuple(0.005695992854379, 0.125094845896272, 0.116460019007601,
00079 0.012270695896559, 0.030770541890982);
00080 q[0] = tuple(1.0, 0.0, 0.0);
00081 q[1] = tuple(0.684472514501908, 0.157763742749046, 0.157763742749046);
00082 q[2] = tuple(0.171245477332074, 0.414377261333963, 0.414377261333963);
00083 q[3] = tuple(0.000000000000000, 0.117472338035268, 0.882527661964732);
00084 q[4] = tuple(0.000000000000000, 0.357384241759678, 0.642615758240322);
00085 }
00086 else if (order == 6)
00087 {
00088 multiplicity = tuple(1, 3, 3, 3, 6, 6, 6);
00089 q.resize(7);
00090 w = tuple(0.078877036160935, -0.001814501354491, 0.021979435329947,
00091 0.055816419204784, 0.013312944349338, 0.013197985920130,
00092 0.089018887113590);
00093 q[0] = tuple(1.0 / 3.0, 1.0 / 3.0, 1.0 / 3.0);
00094 q[1] = tuple(1.0, 0.0, 0.0);
00095 q[2] = tuple(0.0, 0.5, 0.5);
00096 q[3] = tuple(0.769520861253251, 0.115239569373374, 0.115239569373375);
00097 q[4] = tuple(0.000000000000000, 0.084888051860717, 0.915111948139283);
00098 q[5] = tuple(0.000000000000000, 0.265575603264643, 0.734424396735357);
00099 q[6] = tuple(0.127944523240232, 0.317617605129902, 0.554437871629866);
00100 }
00101 else if (order == 9)
00102 {
00103 multiplicity = tuple(1, 3, 3, 3, 3, 6, 6, 6, 6, 6, 6, 6);
00104 q.resize(12);
00105 w = tuple(0.054830037550851, 0.001827879484114, 0.019001036234018,
00106 0.036420799722133, 0.030554508665561, -0.000151089387317,
00107 0.005590374483628, 0.004964725734335, 0.006760863782205,
00108 0.020422905832759, 0.035099032364350, 0.040939402211986);
00109 q[0] = tuple(1.0 / 3.0, 1.0 / 3.0, 1.0 / 3.0);
00110 q[1] = tuple(1.0, 0.0, 0.0);
00111 q[2] = tuple(0.888511356254118, 0.055744321872941, 0.055744321872941);
00112 q[3] = tuple(0.643273196352075, 0.178363401823962, 0.178363401823962);
00113 q[4] = tuple(0.067157296821421, 0.466421351589290, 0.466421351589290);
00114 q[5] = tuple(0.000000000000000, 0.040233045916771, 0.959766954083229);
00115 q[6] = tuple(0.000000000000000, 0.130613067447248, 0.869386932552752);
00116 q[7] = tuple(0.000000000000000, 0.261037525094778, 0.738962474905222);
00117 q[8] = tuple(0.000000000000000, 0.417360521166807, 0.582639478833193);
00118 q[9] = tuple(0.063288094133481, 0.162084689378845, 0.774627216487674);
00119 q[10] = tuple(0.066372381175690, 0.304692583561054, 0.628935035263256);
00120 q[11] = tuple(0.185067879630320, 0.326702931348863, 0.488229189020817);
00121 }
00122
00123
00124
00125
00126
00127
00128
00129
00130
00131
00132
00133
00134
00135
00136
00137
00138
00139
00140
00141
00142
00143
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154
00155
00156
00157
00158
00159
00160 else
00161 {
00162 #ifndef TRILINOS_7
00163 SUNDANCE_ERROR("symmetric Fekete quadrature rule order "
00164 << order <<
00165 " for triangles not available");
00166 #else
00167 SUNDANCE_ERROR7("symmetric Fekete quadrature rule order "
00168 << order <<
00169 " for triangles not available");
00170 #endif
00171 }
00172
00173 for (int i = 0; i < q.length(); i++)
00174 {
00175 Array<Array<double> > qPerm;
00176 permute(multiplicity[i], q[i], qPerm);
00177 for (int j = 0; j < multiplicity[i]; j++)
00178 {
00179 x.append(qPerm[j][0]);
00180 y.append(qPerm[j][1]);
00181 wgt.append(w[i]);
00182 }
00183 }
00184
00185 }
00186
00187 bool FeketeTriangleQuadrature::supportsOrder(int order)
00188 {
00189 if ((order >= 1 && order <= 6) || order == 9)
00190 return true;
00191 return false;
00192 }
00193
00194
00195
00196
00197
00198
00199 void FeketeTriangleQuadrature::computeBasisCoeffs(const int order, Array<double>& basisCoeffs)
00200 {
00201
00202 Array<Point> feketePts;
00203
00204 Array<double> x;
00205 Array<double> y;
00206 Array<double> w;
00207 getPoints(order, w, x, y);
00208 int nFeketePts = w.length();
00209 feketePts.resize(nFeketePts);
00210 for (int i = 0; i < nFeketePts; i++)
00211 feketePts[i] = Point(x[i], y[i]);
00212
00213
00214
00215
00216 basisCoeffs.resize(nFeketePts * nFeketePts);
00217
00218
00219
00220 for (int n = 0; n < nFeketePts; n++)
00221 {
00222
00223
00224 double* start = &(basisCoeffs[n * nFeketePts]);
00225 evalPKDpolynomials(order, feketePts[n][0], feketePts[n][1], start);
00226 }
00227
00228
00229
00230 int lapack_err = 0;
00231 Array<int> pivot;
00232 pivot.resize(nFeketePts);
00233 Array<double> work;
00234 work.resize(1);
00235 int lwork = -1;
00236
00237
00238 ::dgetrf_(&nFeketePts, &nFeketePts, &(basisCoeffs[0]), &nFeketePts,
00239 &(pivot[0]), &lapack_err);
00240
00241 TEUCHOS_TEST_FOR_EXCEPTION(
00242 lapack_err != 0,
00243 std::runtime_error,
00244 "FeketeTriangleQuadrature::computeBasisCoeffs(): factorization of generalized Vandermonde matrix failed");
00245
00246
00247 ::dgetri_(&nFeketePts, &(basisCoeffs[0]), &nFeketePts, &(pivot[0]),
00248 &(work[0]), &lwork, &lapack_err);
00249 lwork = (int) work[0];
00250 work.resize(lwork);
00251 ::dgetri_(&nFeketePts, &(basisCoeffs[0]), &nFeketePts, &(pivot[0]),
00252 &(work[0]), &lwork, &lapack_err);
00253
00254 TEUCHOS_TEST_FOR_EXCEPTION(
00255 lapack_err != 0,
00256 std::runtime_error,
00257 "FeketeTriangleQuadrature::computeBasisCoeffs(): inversion of generalized Vandermonde matrix failed");
00258 }
00259
00260
00261
00262
00263
00264
00265 void FeketeTriangleQuadrature::evalPKDpolynomials(int order, double x,
00266 double y, double* resultPtr)
00267 {
00268 int i = 0;
00269 double Legendre = 1.0;
00270 double Legendre1 = 0.0;
00271 for (int k = 0; k <= order; k++)
00272 {
00273 double Jacobi = 1.0;
00274 double Jacobi1 = 0.0;
00275
00276
00277 for (int l = 0; l <= order - k; l++)
00278 {
00279
00280 resultPtr[i++] = Legendre * pow(x + y, k) * Jacobi;
00281
00282
00283 double Jacobi2 = Jacobi1;
00284 Jacobi1 = Jacobi;
00285 int c1 = 2 * (l + 1) * (l + 2 * k + 2) * (2 * l + 2 * k + 1);
00286 int c2 = (2 * l + 2 * k + 2) * (2 * k + 1) * (2 * k + 1);
00287 int c3 = (2 * l + 2 * k + 1) * (2 * l + 2 * k + 2) * (2 * l + 2 * k
00288 + 3);
00289 int c4 = 2 * (l + 2 * k + 1) * l * (2 * l + 2 * k + 3);
00290 Jacobi = ((c2 + c3 * (1.0 - 2.0 * x - 2.0 * y)) * Jacobi1 - c4
00291 * Jacobi2) / c1;
00292 }
00293
00294
00295 double Legendre2 = Legendre1;
00296 Legendre1 = Legendre;
00297 double xy = x + y;
00298 if (::fabs(xy) > 0.0)
00299 {
00300 xy = (2.0 * k + 1.0) * ((x - y) / xy) * Legendre1;
00301 }
00302 Legendre = (xy - k * Legendre2) / (k + 1);
00303 }
00304 }
00305
00306 void FeketeTriangleQuadrature::permute(int m, const Array<double>& q, Array<
00307 Array<double> >& qPerm)
00308 {
00309 qPerm.resize(m);
00310 if (m == 1)
00311 {
00312 qPerm[0] = q;
00313 }
00314 else if (m == 3)
00315 {
00316 qPerm[0] = tuple(q[0], q[1], q[2]);
00317 qPerm[1] = tuple(q[1], q[0], q[2]);
00318 qPerm[2] = tuple(q[2], q[1], q[0]);
00319 }
00320 else if (m == 6)
00321 {
00322 qPerm[0] = tuple(q[0], q[1], q[2]);
00323 qPerm[1] = tuple(q[0], q[2], q[1]);
00324 qPerm[2] = tuple(q[1], q[0], q[2]);
00325 qPerm[3] = tuple(q[1], q[2], q[0]);
00326 qPerm[4] = tuple(q[2], q[1], q[0]);
00327 qPerm[5] = tuple(q[2], q[0], q[1]);
00328 }
00329 else
00330 {
00331 #ifndef TRILINOS_7
00332 SUNDANCE_ERROR("invalid multiplicity "
00333 << m <<
00334 " in FeketeTriangleQuadrature::permute()");
00335 #else
00336 SUNDANCE_ERROR7("invalid multiplicity "
00337 << m <<
00338 " in FeketeTriangleQuadrature::permute()");
00339 #endif
00340 }
00341 }
00342
00343 bool FeketeTriangleQuadrature::test(int p)
00344 {
00345 Array<double> w;
00346 Array<double> x;
00347 Array<double> y;
00348
00349 getPoints(p, w, x, y);
00350 bool pass = true;
00351
00352 for (int a = 0; a <= p; a++)
00353 {
00354 for (int b = 0; b < p - a; b++)
00355 {
00356 int cMax = p - a - b;
00357 for (int c = 0; c <= cMax; c++)
00358 {
00359 double sum = 0.0;
00360 for (int q = 0; q < w.length(); q++)
00361 {
00362 sum += 0.5 * w[q] * pow(x[q], (double) a) * pow(y[q],
00363 (double) b) * pow(1.0 - x[q] - y[q], (double) c);
00364 }
00365 double err = fabs(sum - exact(a, b, c));
00366 bool localPass = err < 1.0e-14;
00367 pass = pass && localPass;
00368 if (!localPass)
00369 {
00370 fprintf(
00371 stderr,
00372 "order=%d m (%d, %d, %d) q=%22.15g exact=%22.15g\n",
00373 p, a, b, c, sum, exact(a, b, c));
00374 std::cerr << "error = " << err << std::endl;
00375 }
00376 }
00377 }
00378 }
00379 return pass;
00380 }
00381
00382 double FeketeTriangleQuadrature::exact(int a, int b, int c)
00383 {
00384 return fact(a) * fact(b) * fact(c) / fact(a + b + c + 2);
00385 }
00386
00387 double FeketeTriangleQuadrature::fact(int x)
00388 {
00389 if (x == 0)
00390 return 1.0;
00391 return x * fact(x - 1);
00392 }
00393