44#include "Teuchos_UnitTestHarness.hpp"
45#include "Teuchos_TestingHelpers.hpp"
46#include "Teuchos_UnitTestRepository.hpp"
47#include "Teuchos_GlobalMPISession.hpp"
51#ifdef HAVE_STOKHOS_FORUQTK
58 template <
typename OrdinalType,
typename ValueType>
71 int order =
setup.basis.order();
72 TEUCHOS_TEST_EQUALITY(order,
setup.p, out, success);
76 int size =
setup.basis.size();
77 TEUCHOS_TEST_EQUALITY(size,
setup.p+1, out, success);
81 const Teuchos::Array<double>& n1 =
setup.basis.norm_squared();
82 Teuchos::Array<double> n2(
setup.p+1);
83 for (
int i=0; i<=
setup.p; i++)
90 Teuchos::Array<double> n1(
setup.p+1);
91 Teuchos::Array<double> n2(
setup.p+1);
92 for (
int i=0; i<=
setup.p; i++) {
93 n1[i] =
setup.basis.norm_squared(i);
101 const Teuchos::Array<double>& n1 =
setup.basis.norm_squared();
102 Teuchos::Array<double> n2(
setup.p+1);
103 Teuchos::Array<double> x, w;
104 Teuchos::Array< Teuchos::Array<double> > v;
105 setup.basis.getQuadPoints(2*
setup.p, x, w, v);
107 for (
int i=0; i<=
setup.p; i++) {
109 for (
int j=0;
j<nqp;
j++)
110 n2[i] += w[
j]*v[
j][i]*v[
j][i];
118 const Teuchos::Array<double>& norms =
setup.basis.norm_squared();
119 Teuchos::Array<double> x, w;
120 Teuchos::Array< Teuchos::Array<double> > v;
121 setup.basis.getQuadPoints(2*
setup.p, x, w, v);
123 Teuchos::SerialDenseMatrix<int,double> mat(
setup.p+1,
setup.p+1);
124 for (
int i=0; i<=
setup.p; i++) {
126 for (
int k=0; k<nqp; k++)
127 mat(i,
j) += w[k]*v[k][i]*v[k][
j];
128 mat(i,
j) /= std::sqrt(norms[i]*norms[
j]);
132 success = mat.normInf() <
setup.atol;
134 out <<
"\n Error, mat.normInf() < atol = " << mat.normInf()
135 <<
" < " <<
setup.atol <<
": failed!\n";
136 out <<
"mat = " <<
printMat(mat) << std::endl;
141 Teuchos::RCP< Stokhos::Dense3Tensor<int, double> > Cijk =
142 setup.basis.computeTripleProductTensor();
144 Teuchos::Array<double> x, w;
145 Teuchos::Array< Teuchos::Array<double> > v;
146 setup.basis.getQuadPoints(3*
setup.p, x, w, v);
149 for (
int i=0; i<=
setup.p; i++) {
151 for (
int k=0; k<=
setup.p; k++) {
154 for (
int qp=0; qp<nqp; qp++)
155 c += w[qp]*v[qp][i]*v[qp][
j]*v[qp][k];
156 std::stringstream ss;
157 ss <<
"Cijk(" << i <<
"," <<
j <<
"," << k <<
")";
167 Teuchos::RCP< Teuchos::SerialDenseMatrix<int, double> > Bij =
168 setup.basis.computeDerivDoubleProductTensor();
170 Teuchos::Array<double> x, w;
171 Teuchos::Array< Teuchos::Array<double> > v,
val, deriv;
172 setup.basis.getQuadPoints(3*
setup.p, x, w, v);
176 for (
int i=0; i<nqp; i++) {
178 deriv[i].resize(
setup.p+1);
179 setup.basis.evaluateBasesAndDerivatives(x[i],
val[i], deriv[i]);
183 for (
int i=0; i<=
setup.p; i++) {
186 for (
int qp=0; qp<nqp; qp++)
187 b += w[qp]*deriv[qp][i]*
val[qp][
j];
188 std::stringstream ss;
189 ss <<
"Bij(" << i <<
"," <<
j <<
")";
237 Teuchos::Array<double> v1(
setup.p+1), v2(
setup.p+1);
238 setup.basis.evaluateBases(x, v1);
241 b2 = std::sqrt(1.0/3.0);
246 for (
int i=2; i<=
setup.p; i++) {
247 b1 = std::sqrt((i*i)/((2.0*i+1.0)*(2.0*i-1.0)));
248 v2[i] = (x*v2[i-1] - b2*v2[i-2]) / b1;
257 Teuchos::Array<double> val1(
setup.p+1), deriv1(
setup.p+1),
259 setup.basis.evaluateBasesAndDerivatives(x, val1, deriv1);
270 b2 = std::sqrt(1.0/3.0);
275 for (
int i=2; i<=
setup.p; i++) {
276 b1 = std::sqrt((i*i)/((2.0*i+1.0)*(2.0*i-1.0)));
277 val2[i] = (x*val2[i-1] - b2*val2[i-2])/b1;
281 b2 = std::sqrt(1.0/3.0);
286 for (
int i=2; i<=
setup.p; i++) {
287 b1 = std::sqrt((i*i)/((2.0*i+1.0)*(2.0*i-1.0)));
288 deriv2[i] = (val2[i-1] + x*deriv2[i-1] - b2*deriv2[i-2])/b1;
302 Teuchos::Array<double> v1(
setup.p+1), v2(
setup.p+1);
303 for (
int i=0; i<=
setup.p; i++)
304 v1[i] =
setup.basis.evaluate(x, i);
307 b2 = std::sqrt(1.0/3.0);
312 for (
int i=2; i<=
setup.p; i++) {
313 b1 = std::sqrt((i*i)/((2.0*i+1.0)*(2.0*i-1.0)));
314 v2[i] = (x*v2[i-1] - b2*v2[i-2])/b1;
326 setup.basis.getRecurrenceCoefficients(a1, b1, c1, d1);
329 a2[0] = 0.0; b2[0] = 1.0; c2[0] = 1.0; d2[0] = 1.0;
330 for (
int i=1; i<=
setup.p; i++) {
332 b2[i] = std::sqrt((i*i)/((2.0*i+1.0)*(2.0*i-1.0)));
352 setup.basis.getRecurrenceCoefficients(a1, b1, c1, d1);
354 Teuchos::Array<double> a2(
setup.p+1);
355 const Teuchos::Array<double>& n1 =
setup.basis.norm_squared();
356 Teuchos::Array<double> x, w;
357 Teuchos::Array< Teuchos::Array<double> > v;
358 setup.basis.getQuadPoints(2*
setup.p, x, w, v);
360 for (
int i=0; i<=
setup.p; i++) {
362 for (
int j=0;
j<nqp;
j++)
363 a2[i] += w[
j]*x[
j]*v[
j][i]*v[
j][i];
364 a2[i] = a2[i]*c1[i]/n1[i];
376 setup.basis.getRecurrenceCoefficients(a1, b1, c1, d1);
378 Teuchos::Array<double> b2(
setup.p+1);
379 const Teuchos::Array<double>& n1 =
setup.basis.norm_squared();
381 for (
int i=1; i<=
setup.p; i++) {
382 b2[i] = d1[i]*(c1[i]/c1[i-1])*(n1[i]/n1[i-1]);
388#ifdef HAVE_STOKHOS_DAKOTA
390 int n =
static_cast<int>(std::ceil((2*
setup.p+1)/2.0));
391 Teuchos::Array<double> x1, w1;
392 Teuchos::Array< Teuchos::Array<double> > v1;
393 setup.basis.getQuadPoints(2*
setup.p, x1, w1, v1);
395 Teuchos::Array<double> x2(n), w2(n);
396 Teuchos::Array< Teuchos::Array<double> > v2(n);
397 webbur::legendre_compute(n, &x2[0], &w2[0]);
399 for (
int i=0; i<n; i++) {
401 v2[i].resize(
setup.p+1);
402 setup.basis.evaluateBases(x2[i], v2[i]);
409 for (
int i=0; i<n; i++) {
410 std::stringstream ss1, ss2;
411 ss1 <<
"v1[" << i <<
"]";
412 ss2 <<
"v2[" << i <<
"]";
420#ifdef HAVE_STOKHOS_FORUQTK
422 int n =
static_cast<int>(std::ceil((
setup.p+1)/2.0));
423 Teuchos::Array<double> x1, w1;
424 Teuchos::Array< Teuchos::Array<double> > v1;
425 setup.basis.getQuadPoints(
setup.p, x1, w1, v1);
427 Teuchos::Array<double> x2(n), w2(n);
428 Teuchos::Array< Teuchos::Array<double> > v2(n);
431 double endpts[2] = {0.0, 0.0};
432 Teuchos::Array<double> b(n);
435 GAUSSQ_F77(&kind, &n, &alpha, &beta, &kpts, endpts, &b[0], &x2[0], &w2[0]);
437 for (
int i=0; i<n; i++) {
439 v2[i].resize(
setup.p+1);
440 setup.basis.evaluateBases(x2[i], v2[i]);
447 for (
int i=0; i<n; i++) {
448 std::stringstream ss1, ss2;
449 ss1 <<
"v1[" << i <<
"]";
450 ss2 <<
"v2[" << i <<
"]";
461 Teuchos::GlobalMPISession mpiSession(&argc, &
argv);
462 return Teuchos::UnitTestRepository::runUnitTestsFromMain(argc,
argv);
int main(int argc, char *argv[])
Legendre polynomial basis.
UnitTestSetup< int, double > setup
TEUCHOS_UNIT_TEST(Stokhos_LegendreBasis, Order)
bool compareValues(const ValueType &a1, const std::string &a1_name, const ValueType &a2, const std::string &a2_name, const ValueType &rel_tol, const ValueType &abs_tol, Teuchos::FancyOStream &out)
bool compareArrays(const Array1 &a1, const std::string &a1_name, const Array2 &a2, const std::string &a2_name, const ValueType &rel_tol, const ValueType &abs_tol, Teuchos::FancyOStream &out)
Stokhos::LegendreBasis< OrdinalType, ValueType > basis
void printMat(const char *name, Epetra_IntSerialDenseMatrix &matrix)