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::AlgebraicOrthogPolyExpansion<OrdinalType,ValueType> >
exp,
exp_linear;
64 Teuchos::RCP< Stokhos::QuadOrthogPolyExpansion<OrdinalType,ValueType> >
qexp;
65 Stokhos::OrthogPolyApprox<OrdinalType,ValueType> x,
y,
u,
u2,
cx,
cu,
cu2,
sx,
su,
su2;
74 const OrdinalType d = 2;
75 const OrdinalType p = 7;
78 Teuchos::Array< Teuchos::RCP<const Stokhos::OneDOrthogPolyBasis<OrdinalType,ValueType> > > bases(d);
79 for (OrdinalType i=0; i<d; i++)
90 Cijk =
basis->computeTripleProductTensor();
117 for (OrdinalType i=0; i<d; i++) {
122 for (OrdinalType i=0; i<d; i++)
126 template <
class Func>
131 const Teuchos::Array<ValueType>& weights =
quad->getQuadWeights();
132 const Teuchos::Array< Teuchos::Array<ValueType> >& points =
133 quad->getQuadPoints();
134 const Teuchos::Array< Teuchos::Array<ValueType> >& values =
135 quad->getBasisAtQuadPoints();
136 OrdinalType nqp = weights.size();
139 for (OrdinalType i=0; i<c.
size(); i++)
144 for (OrdinalType k=0; k<nqp; k++) {
145 ValueType
val =
a.evaluate(points[k], values[k]);
147 for (
int i=0; i<c.
size(); i++)
148 c[i] += weights[k]*
val*values[k][i] /
basis->norm_squared(i);
152 template <
class Func>
158 const Teuchos::Array<ValueType>& weights =
quad->getQuadWeights();
159 const Teuchos::Array< Teuchos::Array<ValueType> >& points =
160 quad->getQuadPoints();
161 const Teuchos::Array< Teuchos::Array<ValueType> >& values =
162 quad->getBasisAtQuadPoints();
163 OrdinalType nqp = weights.size();
166 for (OrdinalType i=0; i<c.
size(); i++)
171 for (OrdinalType k=0; k<nqp; k++) {
172 ValueType val1 =
a.evaluate(points[k], values[k]);
173 ValueType val2 = b.
evaluate(points[k], values[k]);
174 ValueType
val = func(val1, val2);
175 for (
int i=0; i<c.
size(); i++)
176 c[i] += weights[k]*
val*values[k][i] /
basis->norm_squared(i);
180 template <
class Func>
187 const Teuchos::Array<ValueType>& weights =
quad->getQuadWeights();
188 const Teuchos::Array< Teuchos::Array<ValueType> >& points =
189 quad->getQuadPoints();
190 const Teuchos::Array< Teuchos::Array<ValueType> >& values =
191 quad->getBasisAtQuadPoints();
192 OrdinalType nqp = weights.size();
195 for (OrdinalType i=0; i<c.
size(); i++)
200 for (OrdinalType k=0; k<nqp; k++) {
201 ValueType val2 = b.
evaluate(points[k], values[k]);
202 ValueType
val = func(
a, val2);
203 for (
int i=0; i<c.
size(); i++)
204 c[i] += weights[k]*
val*values[k][i] /
basis->norm_squared(i);
208 template <
class Func>
215 const Teuchos::Array<ValueType>& weights =
quad->getQuadWeights();
216 const Teuchos::Array< Teuchos::Array<ValueType> >& points =
217 quad->getQuadPoints();
218 const Teuchos::Array< Teuchos::Array<ValueType> >& values =
219 quad->getBasisAtQuadPoints();
220 OrdinalType nqp = weights.size();
223 for (OrdinalType i=0; i<c.
size(); i++)
228 for (OrdinalType k=0; k<nqp; k++) {
229 ValueType val1 =
a.evaluate(points[k], values[k]);
230 ValueType
val = func(val1, b);
231 for (
int i=0; i<c.
size(); i++)
232 c[i] += weights[k]*
val*values[k][i] /
basis->norm_squared(i);
246 return std::log(a+std::sqrt(a*a+1.0));
251 return std::log(a+std::sqrt(a*a-1.0));
256 return 0.5*std::log((1.0+a)/(1.0-a));
261 double operator() (
double a,
double b)
const {
return a + b; }
264 double operator() (
double a,
double b)
const {
return a - b; }
267 double operator() (
double a,
double b)
const {
return a * b; }
270 double operator() (
double a,
double b)
const {
return a / b; }
439 setup.exp->plus(ru, v, w);
555 setup.exp->minus(ru, v, w);
671 setup.exp->times(ru, v, w);
830 setup.exp->plusEqual(ru, v);
883 setup.exp->minusEqual(ru, v);
884 setup.exp->unaryMinus(v, v);
1006 Teuchos::GlobalMPISession mpiSession(&argc, &
argv);
1007 return Teuchos::UnitTestRepository::runUnitTestsFromMain(argc,
argv);
int main(int argc, char *argv[])
Orthogonal polynomial expansions limited to algebraic operations.
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)
UnitTestSetup< int, double > setup
TEUCHOS_UNIT_TEST(Stokhos_AlgebraicExpansion, 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, double b) 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
Teuchos::RCP< Stokhos::AlgebraicOrthogPolyExpansion< OrdinalType, ValueType > > exp
void computePCE2RC(Stokhos::OrthogPolyApprox< OrdinalType, ValueType > &c, const Stokhos::OrthogPolyApprox< OrdinalType, ValueType > &a, ValueType b)
Stokhos::OrthogPolyApprox< OrdinalType, ValueType > u2
Stokhos::OrthogPolyApprox< OrdinalType, ValueType > x
Stokhos::OrthogPolyApprox< OrdinalType, ValueType > u
Stokhos::OrthogPolyApprox< OrdinalType, ValueType > su
Stokhos::OrthogPolyApprox< OrdinalType, ValueType > cu2
Stokhos::OrthogPolyApprox< OrdinalType, ValueType > cu
Stokhos::OrthogPolyApprox< OrdinalType, ValueType > cx
Teuchos::RCP< const Stokhos::CompletePolynomialBasis< OrdinalType, ValueType > > basis
Stokhos::OrthogPolyApprox< OrdinalType, ValueType > y
Stokhos::OrthogPolyApprox< OrdinalType, ValueType > su2
Teuchos::RCP< const Stokhos::Quadrature< OrdinalType, ValueType > > quad
void computePCE2(Stokhos::OrthogPolyApprox< OrdinalType, ValueType > &c, const Stokhos::OrthogPolyApprox< OrdinalType, ValueType > &a, const Stokhos::OrthogPolyApprox< OrdinalType, ValueType > &b)
void computePCE2LC(Stokhos::OrthogPolyApprox< OrdinalType, ValueType > &c, ValueType a, const Stokhos::OrthogPolyApprox< OrdinalType, ValueType > &b)
Teuchos::RCP< Stokhos::Sparse3Tensor< int, double > > Cijk
void computePCE1(Stokhos::OrthogPolyApprox< OrdinalType, ValueType > &c, const Stokhos::OrthogPolyApprox< OrdinalType, ValueType > &a)
Stokhos::OrthogPolyApprox< OrdinalType, ValueType > sx
Teuchos::RCP< Stokhos::Sparse3Tensor< int, double > > Cijk_linear
Teuchos::RCP< Stokhos::AlgebraicOrthogPolyExpansion< OrdinalType, ValueType > > exp_linear
Teuchos::RCP< Stokhos::QuadOrthogPolyExpansion< OrdinalType, ValueType > > qexp