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<const Stokhos::Quadrature<OrdinalType,ValueType> >
quad;
63 Teuchos::RCP< Stokhos::QuadOrthogPolyExpansion<OrdinalType,ValueType> >
exp,
exp_linear;
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++)
89 Cijk =
basis->computeTripleProductTensor();
112 for (OrdinalType i=0; i<d; i++) {
117 for (OrdinalType i=0; i<d; i++)
121 template <
class Func>
126 const Teuchos::Array<ValueType>& weights =
quad->getQuadWeights();
127 const Teuchos::Array< Teuchos::Array<ValueType> >& points =
128 quad->getQuadPoints();
129 const Teuchos::Array< Teuchos::Array<ValueType> >& values =
130 quad->getBasisAtQuadPoints();
131 OrdinalType nqp = weights.size();
134 for (OrdinalType i=0; i<c.
size(); i++)
139 for (OrdinalType k=0; k<nqp; k++) {
140 ValueType
val =
a.evaluate(points[k], values[k]);
142 for (
int i=0; i<c.
size(); i++)
143 c[i] += weights[k]*
val*values[k][i] /
basis->norm_squared(i);
147 template <
class Func>
153 const Teuchos::Array<ValueType>& weights =
quad->getQuadWeights();
154 const Teuchos::Array< Teuchos::Array<ValueType> >& points =
155 quad->getQuadPoints();
156 const Teuchos::Array< Teuchos::Array<ValueType> >& values =
157 quad->getBasisAtQuadPoints();
158 OrdinalType nqp = weights.size();
161 for (OrdinalType i=0; i<c.
size(); i++)
166 for (OrdinalType k=0; k<nqp; k++) {
167 ValueType val1 =
a.evaluate(points[k], values[k]);
168 ValueType val2 = b.
evaluate(points[k], values[k]);
169 ValueType
val = func(val1, val2);
170 for (
int i=0; i<c.
size(); i++)
171 c[i] += weights[k]*
val*values[k][i] /
basis->norm_squared(i);
175 template <
class Func>
182 const Teuchos::Array<ValueType>& weights =
quad->getQuadWeights();
183 const Teuchos::Array< Teuchos::Array<ValueType> >& points =
184 quad->getQuadPoints();
185 const Teuchos::Array< Teuchos::Array<ValueType> >& values =
186 quad->getBasisAtQuadPoints();
187 OrdinalType nqp = weights.size();
190 for (OrdinalType i=0; i<c.
size(); i++)
195 for (OrdinalType k=0; k<nqp; k++) {
196 ValueType val2 = b.
evaluate(points[k], values[k]);
197 ValueType
val = func(
a, val2);
198 for (
int i=0; i<c.
size(); i++)
199 c[i] += weights[k]*
val*values[k][i] /
basis->norm_squared(i);
203 template <
class Func>
210 const Teuchos::Array<ValueType>& weights =
quad->getQuadWeights();
211 const Teuchos::Array< Teuchos::Array<ValueType> >& points =
212 quad->getQuadPoints();
213 const Teuchos::Array< Teuchos::Array<ValueType> >& values =
214 quad->getBasisAtQuadPoints();
215 OrdinalType nqp = weights.size();
218 for (OrdinalType i=0; i<c.
size(); i++)
223 for (OrdinalType k=0; k<nqp; k++) {
224 ValueType val1 =
a.evaluate(points[k], values[k]);
225 ValueType
val = func(val1, b);
226 for (
int i=0; i<c.
size(); i++)
227 c[i] += weights[k]*
val*values[k][i] /
basis->norm_squared(i);
282 return std::log(a+std::sqrt(a*a+1.0));
287 return std::log(a+std::sqrt(a*a-1.0));
292 return 0.5*std::log((1.0+a)/(1.0-a));
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 a / b; }
309 double operator() (
double a,
double b)
const {
return std::pow(a,b); }
699 setup.exp->plus(ru, v, w);
815 setup.exp->minus(ru, v, w);
931 setup.exp->times(ru, v, w);
1063 setup.exp->divide(ru, v, w);
1248 setup.exp->plusEqual(ru, v);
1301 setup.exp->minusEqual(ru, v);
1302 setup.exp->unaryMinus(v, v);
1449 Teuchos::GlobalMPISession mpiSession(&argc, &
argv);
1450 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.
Orthogonal polynomial expansions based on numerical quadrature.
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)
TEUCHOS_UNIT_TEST(Stokhos_QuadExpansion, UMinus)
UnitTestSetup< int, double > setup
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
void computePCE2LC(Stokhos::OrthogPolyApprox< OrdinalType, ValueType > &c, ValueType a, const Stokhos::OrthogPolyApprox< OrdinalType, ValueType > &b)
void computePCE2RC(Stokhos::OrthogPolyApprox< OrdinalType, ValueType > &c, const Stokhos::OrthogPolyApprox< OrdinalType, ValueType > &a, ValueType b)
Stokhos::OrthogPolyApprox< OrdinalType, ValueType > x
void computePCE2(Stokhos::OrthogPolyApprox< OrdinalType, ValueType > &c, const Stokhos::OrthogPolyApprox< OrdinalType, ValueType > &a, const Stokhos::OrthogPolyApprox< OrdinalType, ValueType > &b)
Stokhos::OrthogPolyApprox< OrdinalType, ValueType > cx
Stokhos::OrthogPolyApprox< OrdinalType, ValueType > sx
Teuchos::RCP< Stokhos::QuadOrthogPolyExpansion< OrdinalType, ValueType > > exp
Stokhos::OrthogPolyApprox< OrdinalType, ValueType > su
Teuchos::RCP< const Stokhos::Quadrature< OrdinalType, ValueType > > quad
Teuchos::RCP< Stokhos::QuadOrthogPolyExpansion< OrdinalType, ValueType > > exp_linear
Stokhos::OrthogPolyApprox< OrdinalType, ValueType > cu2
Stokhos::OrthogPolyApprox< OrdinalType, ValueType > y
Stokhos::OrthogPolyApprox< OrdinalType, ValueType > su2
Teuchos::RCP< Stokhos::Sparse3Tensor< int, double > > Cijk
void computePCE1(Stokhos::OrthogPolyApprox< OrdinalType, ValueType > &c, const Stokhos::OrthogPolyApprox< OrdinalType, ValueType > &a)
Teuchos::RCP< Stokhos::Sparse3Tensor< int, double > > Cijk_linear
Stokhos::OrthogPolyApprox< OrdinalType, ValueType > u2
Stokhos::OrthogPolyApprox< OrdinalType, ValueType > u
Stokhos::OrthogPolyApprox< OrdinalType, ValueType > cu
Teuchos::RCP< const Stokhos::CompletePolynomialBasis< OrdinalType, ValueType > > basis