00001 #include "SundanceQuadQuadrature.hpp"
00002 #include "SundanceOut.hpp"
00003 #include "SundanceGauss1D.hpp"
00004 #include "PlayaTabs.hpp"
00005
00006 using namespace Sundance;
00007 using namespace Teuchos;
00008
00009 void QuadQuadrature::getPoints(int order, Array<double>& wgt, Array<double>& x,
00010 Array<double>& y)
00011 {
00012 int p = order + 1;
00013 p = p + (p % 2);
00014 int nNodes = p / 2;
00015 Gauss1D rule(nNodes, 0.0, 1.0);
00016 Array<double> s = rule.nodes();
00017 Array<double> t = s;
00018 Array<double> w = rule.weights();
00019 int n = rule.nPoints();
00020
00021 wgt.resize(n * n);
00022 x.resize(n * n);
00023 y.resize(n * n);
00024
00025 int k = 0;
00026 for (int i = 0; i < n; i++)
00027 {
00028 double p = s[i];
00029 for (int j = 0; j < n; j++, k++)
00030 {
00031 double q = t[j];
00032 x[k] = p;
00033 y[k] = q;
00034 wgt[k] = w[i] * w[j];
00035 }
00036 }
00037 }
00038
00039 bool QuadQuadrature::test(int p)
00040 {
00041 Array<double> w;
00042 Array<double> x;
00043 Array<double> y;
00044
00045 getPoints(p, w, x, y);
00046 bool pass = true;
00047 for (int a = 0; a <= p; a++)
00048 {
00049 int bMax = p - a;
00050 for (int b = 0; b <= bMax; b++)
00051 {
00052 double sum = 0.0;
00053 for (int q = 0; q < w.length(); q++)
00054 {
00055 sum += w[q] * pow(x[q], (double) a) * pow(y[q], (double) b);
00056 }
00057 double err = fabs(sum - exact(a, b));
00058 bool localPass = err < 1.0e-14;
00059 pass = pass && localPass;
00060 if (!localPass)
00061 {
00062 fprintf(stderr,
00063 "order=%d m (%d, %d) q=%22.15g exact=%22.15g\n", p, a,
00064 b, sum, exact(a, b));
00065 std::cerr << "error = " << err << std::endl;
00066 }
00067 }
00068 }
00069 return pass;
00070 }
00071
00072 double QuadQuadrature::exact(int a, int b)
00073 {
00074 return 1.0 / (a + 1) / (b + 1);
00075 }
00076