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