58 template <
typename OrdinalType,
typename ValueType>
68 UnitTestSetup<int,double>
setup;
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);
84 for (
int i=1; i<=
setup.p; i++)
91 Teuchos::Array<double> n1(
setup.p+1);
92 Teuchos::Array<double> n2(
setup.p+1);
93 n1[0] =
setup.basis.norm_squared(0);
95 for (
int i=1; i<=
setup.p; i++) {
96 n1[i] =
setup.basis.norm_squared(i);
104 const Teuchos::Array<double>& n1 =
setup.basis.norm_squared();
105 Teuchos::Array<double> n2(
setup.p+1);
106 Teuchos::Array<double> x, w;
107 Teuchos::Array< Teuchos::Array<double> > v;
108 setup.basis.getQuadPoints(2*
setup.p, x, w, v);
110 for (
int i=0; i<=
setup.p; i++) {
112 for (
int j=0;
j<nqp;
j++)
113 n2[i] += w[
j]*v[
j][i]*v[
j][i];
121 const Teuchos::Array<double>& norms =
setup.basis.norm_squared();
122 Teuchos::Array<double> x, w;
123 Teuchos::Array< Teuchos::Array<double> > v;
124 setup.basis.getQuadPoints(2*
setup.p, x, w, v);
126 Teuchos::SerialDenseMatrix<int,double> mat(
setup.p+1,
setup.p+1);
127 for (
int i=0; i<=
setup.p; i++) {
129 for (
int k=0; k<nqp; k++)
130 mat(i,
j) += w[k]*v[k][i]*v[k][
j];
131 mat(i,
j) /= std::sqrt(norms[i]*norms[
j]);
135 success = mat.normInf() <
setup.atol;
137 out <<
"\n Error, mat.normInf() < atol = " << mat.normInf()
138 <<
" < " <<
setup.atol <<
": failed!\n";
139 out <<
"mat = " <<
printMat(mat) << std::endl;
144 Teuchos::RCP< Stokhos::Dense3Tensor<int, double> > Cijk =
145 setup.basis.computeTripleProductTensor();
147 Teuchos::Array<double> x, w;
148 Teuchos::Array< Teuchos::Array<double> > v;
149 setup.basis.getQuadPoints(3*
setup.p, x, w, v);
152 for (
int i=0; i<=
setup.p; i++) {
154 for (
int k=0; k<=
setup.p; k++) {
157 for (
int qp=0; qp<nqp; qp++)
158 c += w[qp]*v[qp][i]*v[qp][
j]*v[qp][k];
159 std::stringstream ss;
160 ss <<
"Cijk(" << i <<
"," <<
j <<
"," << k <<
")";
170 Teuchos::RCP< Teuchos::SerialDenseMatrix<int, double> > Bij =
171 setup.basis.computeDerivDoubleProductTensor();
173 Teuchos::Array<double> x, w;
174 Teuchos::Array< Teuchos::Array<double> > v,
val, deriv;
175 setup.basis.getQuadPoints(2*
setup.p, x, w, v);
179 for (
int i=0; i<nqp; i++) {
181 deriv[i].resize(
setup.p+1);
182 setup.basis.evaluateBasesAndDerivatives(x[i],
val[i], deriv[i]);
186 for (
int i=0; i<=
setup.p; i++) {
189 for (
int qp=0; qp<nqp; qp++)
190 b += w[qp]*deriv[qp][i]*
val[qp][
j];
191 std::stringstream ss;
192 ss <<
"Bij(" << i <<
"," <<
j <<
")";
201 Teuchos::RCP< Teuchos::SerialDenseMatrix<int, double> > Bij =
202 setup.basis.computeDerivDoubleProductTensor();
203 const Teuchos::Array<double>& n =
setup.basis.norm_squared();
206 for (
int i=0; i<=
setup.p; i++) {
211 std::stringstream ss;
212 ss <<
"Bij(" << i <<
"," <<
j <<
")";
222 Teuchos::Array<double> v1(
setup.p+1), v2(
setup.p+1);
223 setup.basis.evaluateBases(x, v1);
232 for (
int i=2; i<=
setup.p; i++)
233 v2[i] = x*v2[i-1] - (i-1.0)*v2[i-2];
240 Teuchos::Array<double> val1(
setup.p+1), deriv1(
setup.p+1),
242 setup.basis.evaluateBasesAndDerivatives(x, val1, deriv1);
251 for (
int i=2; i<=
setup.p; i++)
252 val2[i] = x*val2[i-1] - (i-1.0)*val2[i-2];
257 for (
int i=2; i<=
setup.p; i++)
258 deriv2[i] = i*val2[i-1];
270 Teuchos::Array<double> v1(
setup.p+1), v2(
setup.p+1);
271 for (
int i=0; i<=
setup.p; i++)
272 v1[i] =
setup.basis.evaluate(x, i);
281 for (
int i=2; i<=
setup.p; i++)
282 v2[i] = x*v2[i-1] - (i-1.0)*v2[i-2];
292 setup.basis.getRecurrenceCoefficients(a1, b1, c1, d1);
295 a2[0] = 0.0; b2[0] = 1.0; c2[0] = 1.0; d2[0] = 1.0;
296 for (
int i=1; i<=
setup.p; i++) {
318 setup.basis.getRecurrenceCoefficients(a1, b1, c1, d1);
320 Teuchos::Array<double> a2(
setup.p+1);
321 const Teuchos::Array<double>& n1 =
setup.basis.norm_squared();
322 Teuchos::Array<double> x, w;
323 Teuchos::Array< Teuchos::Array<double> > v;
324 setup.basis.getQuadPoints(2*
setup.p, x, w, v);
326 for (
int i=0; i<=
setup.p; i++) {
328 for (
int j=0;
j<nqp;
j++)
329 a2[i] += w[
j]*x[
j]*v[
j][i]*v[
j][i];
330 a2[i] = a2[i]*c1[i]/n1[i];
342 setup.basis.getRecurrenceCoefficients(a1, b1, c1, d1);
344 Teuchos::Array<double> b2(
setup.p+1);
345 const Teuchos::Array<double>& n1 =
setup.basis.norm_squared();
347 for (
int i=1; i<=
setup.p; i++) {
348 b2[i] = d1[i]*(c1[i]/c1[i-1])*(n1[i]/n1[i-1]);
354#ifdef HAVE_STOKHOS_DAKOTA
356 int n =
static_cast<int>(std::ceil((
setup.p+1)/2.0));
357 Teuchos::Array<double> x1, w1;
358 Teuchos::Array< Teuchos::Array<double> > v1;
359 setup.basis.getQuadPoints(
setup.p, x1, w1, v1);
361 Teuchos::Array<double> x2(n), w2(n);
362 Teuchos::Array< Teuchos::Array<double> > v2(n);
363 webbur::hermite_compute(n, &x2[0], &w2[0]);
365 for (
int i=0; i<n; i++) {
366 w2[i] *= 0.5/std::sqrt(std::atan(1.0));
367 x2[i] *= std::sqrt(2.0);
368 v2[i].resize(
setup.p+1);
369 setup.basis.evaluateBases(x2[i], v2[i]);
376 for (
int i=0; i<n; i++) {
377 std::stringstream ss1, ss2;
378 ss1 <<
"v1[" << i <<
"]";
379 ss2 <<
"v2[" << i <<
"]";
387#ifdef HAVE_STOKHOS_FORUQTK
389 int n =
static_cast<int>(std::ceil((
setup.p+1)/2.0));
390 Teuchos::Array<double> x1, w1;
391 Teuchos::Array< Teuchos::Array<double> > v1;
392 setup.basis.getQuadPoints(
setup.p, x1, w1, v1);
394 Teuchos::Array<double> x2(n), w2(n);
395 Teuchos::Array< Teuchos::Array<double> > v2(n);
398 double endpts[2] = {0.0, 0.0};
399 Teuchos::Array<double> b(n);
402 GAUSSQ_F77(&kind, &n, &alpha, &beta, &kpts, endpts, &b[0], &x2[0], &w2[0]);
404 for (
int i=0; i<n; i++) {
405 w2[i] *= 0.5/std::sqrt(std::atan(1.0));
406 x2[i] *= std::sqrt(2.0);
407 v2[i].resize(
setup.p+1);
408 setup.basis.evaluateBases(x2[i], v2[i]);
415 for (
int i=0; i<n; i++) {
416 std::stringstream ss1, ss2;
417 ss1 <<
"v1[" << i <<
"]";
418 ss2 <<
"v2[" << i <<
"]";
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)
void printMat(const char *name, Epetra_IntSerialDenseMatrix &matrix)