44#include "Teuchos_UnitTestHarness.hpp"
45#include "Teuchos_TestingHelpers.hpp"
46#include "Teuchos_UnitTestRepository.hpp"
47#include "Teuchos_GlobalMPISession.hpp"
55 template <
typename OrdinalType,
typename ValueType>
60 Teuchos::RCP<const Stokhos::CompletePolynomialBasis<OrdinalType,ValueType> >
basis;
61 Teuchos::RCP<Stokhos::Sparse3Tensor<int,double> >
Cijk;
62 Teuchos::RCP<const Stokhos::Quadrature<OrdinalType,ValueType> >
quad;
63 Teuchos::RCP< Stokhos::ForUQTKOrthogPolyExpansion<OrdinalType,ValueType> >
exp;
64 Stokhos::OrthogPolyApprox<OrdinalType,ValueType> x,
y,
u,
u2,
cx,
cu,
cu2,
sx,
su,
su2;
73 const OrdinalType d = 2;
74 const OrdinalType p = 7;
77 Teuchos::Array< Teuchos::RCP<const Stokhos::OneDOrthogPolyBasis<OrdinalType,ValueType> > > bases(d);
78 for (OrdinalType i=0; i<d; i++)
85 Cijk =
basis->computeTripleProductTensor();
93 Teuchos::rcp(
new Stokhos::ForUQTKOrthogPolyExpansion<OrdinalType,ValueType>(
basis,
Cijk, Stokhos::ForUQTKOrthogPolyExpansion<OrdinalType,ValueType>::INTEGRATION));
109 for (OrdinalType i=0; i<d; i++) {
114 for (OrdinalType i=0; i<d; i++)
118 template <
class Func>
123 const Teuchos::Array<ValueType>& weights =
quad->getQuadWeights();
124 const Teuchos::Array< Teuchos::Array<ValueType> >& points =
125 quad->getQuadPoints();
126 const Teuchos::Array< Teuchos::Array<ValueType> >& values =
127 quad->getBasisAtQuadPoints();
128 OrdinalType nqp = weights.size();
131 for (OrdinalType i=0; i<c.
size(); i++)
136 for (OrdinalType k=0; k<nqp; k++) {
137 ValueType
val =
a.evaluate(points[k], values[k]);
139 for (
int i=0; i<c.
size(); i++)
140 c[i] += weights[k]*
val*values[k][i] /
basis->norm_squared(i);
144 template <
class Func>
150 const Teuchos::Array<ValueType>& weights =
quad->getQuadWeights();
151 const Teuchos::Array< Teuchos::Array<ValueType> >& points =
152 quad->getQuadPoints();
153 const Teuchos::Array< Teuchos::Array<ValueType> >& values =
154 quad->getBasisAtQuadPoints();
155 OrdinalType nqp = weights.size();
158 for (OrdinalType i=0; i<c.
size(); i++)
163 for (OrdinalType k=0; k<nqp; k++) {
164 ValueType val1 =
a.
evaluate(points[k], values[k]);
165 ValueType val2 = b.
evaluate(points[k], values[k]);
166 ValueType
val = func(val1, val2);
167 for (
int i=0; i<c.
size(); i++)
168 c[i] += weights[k]*
val*values[k][i] /
basis->norm_squared(i);
172 template <
class Func>
179 const Teuchos::Array<ValueType>& weights =
quad->getQuadWeights();
180 const Teuchos::Array< Teuchos::Array<ValueType> >& points =
181 quad->getQuadPoints();
182 const Teuchos::Array< Teuchos::Array<ValueType> >& values =
183 quad->getBasisAtQuadPoints();
184 OrdinalType nqp = weights.size();
187 for (OrdinalType i=0; i<c.
size(); i++)
192 for (OrdinalType k=0; k<nqp; k++) {
193 ValueType val2 = b.
evaluate(points[k], values[k]);
194 ValueType
val = func(
a, val2);
195 for (
int i=0; i<c.
size(); i++)
196 c[i] += weights[k]*
val*values[k][i] /
basis->norm_squared(i);
200 template <
class Func>
207 const Teuchos::Array<ValueType>& weights =
quad->getQuadWeights();
208 const Teuchos::Array< Teuchos::Array<ValueType> >& points =
209 quad->getQuadPoints();
210 const Teuchos::Array< Teuchos::Array<ValueType> >& values =
211 quad->getBasisAtQuadPoints();
212 OrdinalType nqp = weights.size();
215 for (OrdinalType i=0; i<c.
size(); i++)
220 for (OrdinalType k=0; k<nqp; k++) {
221 ValueType val1 =
a.
evaluate(points[k], values[k]);
222 ValueType
val = func(val1, b);
223 for (
int i=0; i<c.
size(); i++)
224 c[i] += weights[k]*
val*values[k][i] /
basis->norm_squared(i);
279 return std::log(a+std::sqrt(a*a+1.0));
284 return std::log(a+std::sqrt(a*a-1.0));
289 return 0.5*std::log((1.0+a)/(1.0-a));
294 double operator() (
double a,
double b)
const {
return a + b; }
297 double operator() (
double a,
double b)
const {
return a - b; }
300 double operator() (
double a,
double b)
const {
return a * b; }
303 double operator() (
double a,
double b)
const {
return a / b; }
306 double operator() (
double a,
double b)
const {
return std::pow(a,b); }
821 setup.exp->plus(ru, v, w);
937 setup.exp->minus(ru, v, w);
1053 setup.exp->times(ru, v, w);
1121 setup.exp->divide(ru, v, w);
1230 setup.exp->plusEqual(ru, v);
1283 setup.exp->minusEqual(ru, v);
1284 setup.exp->unaryMinus(v, v);
1382 Teuchos::GlobalMPISession mpiSession(&argc, &
argv);
1383 return Teuchos::UnitTestRepository::runUnitTestsFromMain(argc,
argv);
int main(int argc, char *argv[])
Multivariate orthogonal polynomial basis generated from a total-order complete-polynomial tensor prod...
Legendre polynomial basis.
Class to store coefficients of a projection onto an orthogonal polynomial basis.
void reset(const Teuchos::RCP< const Stokhos::OrthogPolyBasis< ordinal_type, value_type > > &new_basis, ordinal_type sz=0)
Reset to a new basis.
reference term(ordinal_type dimension, ordinal_type order)
Get coefficient term for given dimension and order.
ordinal_type size() const
Return size.
value_type evaluate(const Teuchos::Array< value_type > &point) const
Evaluate polynomial approximation at a point.
Defines quadrature for a tensor product basis by tensor products of 1-D quadrature rules.
ScalarType f(const Teuchos::Array< ScalarType > &x, double a, double b)
UnitTestSetup< int, double > setup
TEUCHOS_UNIT_TEST(Stokhos_ForUQTKExpansion, UMinus)
bool comparePCEs(const PCEType &a1, const std::string &a1_name, const Stokhos::OrthogPolyApprox< OrdinalType, ValueType > &a2, const std::string &a2_name, const ValueType &rel_tol, const ValueType &abs_tol, Teuchos::FancyOStream &out)
double operator()(double a) const
double operator()(double a) const
double operator()(double a) const
double operator()(double a) const
double operator()(double a) const
double operator()(double a) const
double operator()(double a) const
double operator()(double a) const
double operator()(double a) const
double operator()(double a, double b) const
double operator()(double a) const
double operator()(double a) const
double operator()(double a) const
double operator()(double a, double b) const
double operator()(double a, double b) const
double operator()(double a, double b) const
double operator()(double a) const
double operator()(double a) const
double operator()(double a) const
double operator()(double a) const
double operator()(double a) const
double operator()(double a, double b) const
double operator()(double a) const
Stokhos::OrthogPolyApprox< OrdinalType, ValueType > u2
Stokhos::OrthogPolyApprox< OrdinalType, ValueType > u
Stokhos::OrthogPolyApprox< OrdinalType, ValueType > y
void computePCE2RC(Stokhos::OrthogPolyApprox< OrdinalType, ValueType > &c, const Stokhos::OrthogPolyApprox< OrdinalType, ValueType > &a, ValueType b)
void computePCE1(Stokhos::OrthogPolyApprox< OrdinalType, ValueType > &c, const Stokhos::OrthogPolyApprox< OrdinalType, ValueType > &a)
Stokhos::OrthogPolyApprox< OrdinalType, ValueType > x
void computePCE2LC(Stokhos::OrthogPolyApprox< OrdinalType, ValueType > &c, ValueType a, const Stokhos::OrthogPolyApprox< OrdinalType, ValueType > &b)
Stokhos::OrthogPolyApprox< OrdinalType, ValueType > su
void computePCE2(Stokhos::OrthogPolyApprox< OrdinalType, ValueType > &c, const Stokhos::OrthogPolyApprox< OrdinalType, ValueType > &a, const Stokhos::OrthogPolyApprox< OrdinalType, ValueType > &b)
Teuchos::RCP< Stokhos::ForUQTKOrthogPolyExpansion< OrdinalType, ValueType > > exp
Stokhos::OrthogPolyApprox< OrdinalType, ValueType > cu2
Teuchos::RCP< const Stokhos::Quadrature< OrdinalType, ValueType > > quad
Teuchos::RCP< const Stokhos::CompletePolynomialBasis< OrdinalType, ValueType > > basis
Stokhos::OrthogPolyApprox< OrdinalType, ValueType > sx
Teuchos::RCP< Stokhos::Sparse3Tensor< int, double > > Cijk
Stokhos::OrthogPolyApprox< OrdinalType, ValueType > su2
Stokhos::OrthogPolyApprox< OrdinalType, ValueType > cx
Stokhos::OrthogPolyApprox< OrdinalType, ValueType > cu