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<Teuchos::SerialDenseMatrix<int,double> >
Bij;
62 Teuchos::RCP<Stokhos::Sparse3Tensor<int,double> >
Cijk;
63 Teuchos::RCP<Stokhos::Dense3Tensor<int,double> >
Dijk;
64 Teuchos::RCP<const Stokhos::Quadrature<OrdinalType,ValueType> >
quad;
65 Teuchos::RCP< Stokhos::DerivOrthogPolyExpansion<OrdinalType,ValueType> >
exp;
66 Stokhos::OrthogPolyApprox<OrdinalType,ValueType> x,
y,
u,
u2,
cx,
cu,
cu2,
sx,
su,
su2;
75 const OrdinalType d = 1;
76 const OrdinalType p = 7;
79 Teuchos::Array< Teuchos::RCP<const Stokhos::OneDOrthogPolyBasis<OrdinalType,ValueType> > > bases(d);
80 for (OrdinalType i=0; i<d; i++)
91 Bij =
basis->computeDerivDoubleProductTensor();
92 Cijk =
basis->computeTripleProductTensor();
113 for (OrdinalType i=0; i<d; i++) {
118 for (OrdinalType i=0; i<d; i++)
122 template <
class Func>
127 const Teuchos::Array<ValueType>& weights =
quad->getQuadWeights();
128 const Teuchos::Array< Teuchos::Array<ValueType> >& points =
129 quad->getQuadPoints();
130 const Teuchos::Array< Teuchos::Array<ValueType> >& values =
131 quad->getBasisAtQuadPoints();
132 OrdinalType nqp = weights.size();
135 for (OrdinalType i=0; i<c.
size(); i++)
140 for (OrdinalType k=0; k<nqp; k++) {
141 ValueType
val =
a.evaluate(points[k], values[k]);
143 for (
int i=0; i<c.
size(); i++)
144 c[i] += weights[k]*
val*values[k][i] /
basis->norm_squared(i);
148 template <
class Func>
154 const Teuchos::Array<ValueType>& weights =
quad->getQuadWeights();
155 const Teuchos::Array< Teuchos::Array<ValueType> >& points =
156 quad->getQuadPoints();
157 const Teuchos::Array< Teuchos::Array<ValueType> >& values =
158 quad->getBasisAtQuadPoints();
159 OrdinalType nqp = weights.size();
162 for (OrdinalType i=0; i<c.
size(); i++)
167 for (OrdinalType k=0; k<nqp; k++) {
168 ValueType val1 =
a.
evaluate(points[k], values[k]);
169 ValueType val2 = b.
evaluate(points[k], values[k]);
170 ValueType
val = func(val1, val2);
171 for (
int i=0; i<c.
size(); i++)
172 c[i] += weights[k]*
val*values[k][i] /
basis->norm_squared(i);
176 template <
class Func>
183 const Teuchos::Array<ValueType>& weights =
quad->getQuadWeights();
184 const Teuchos::Array< Teuchos::Array<ValueType> >& points =
185 quad->getQuadPoints();
186 const Teuchos::Array< Teuchos::Array<ValueType> >& values =
187 quad->getBasisAtQuadPoints();
188 OrdinalType nqp = weights.size();
191 for (OrdinalType i=0; i<c.
size(); i++)
196 for (OrdinalType k=0; k<nqp; k++) {
197 ValueType val2 = b.
evaluate(points[k], values[k]);
198 ValueType
val = func(
a, val2);
199 for (
int i=0; i<c.
size(); i++)
200 c[i] += weights[k]*
val*values[k][i] /
basis->norm_squared(i);
204 template <
class Func>
211 const Teuchos::Array<ValueType>& weights =
quad->getQuadWeights();
212 const Teuchos::Array< Teuchos::Array<ValueType> >& points =
213 quad->getQuadPoints();
214 const Teuchos::Array< Teuchos::Array<ValueType> >& values =
215 quad->getBasisAtQuadPoints();
216 OrdinalType nqp = weights.size();
219 for (OrdinalType i=0; i<c.
size(); i++)
224 for (OrdinalType k=0; k<nqp; k++) {
225 ValueType val1 =
a.
evaluate(points[k], values[k]);
226 ValueType
val = func(val1, b);
227 for (
int i=0; i<c.
size(); i++)
228 c[i] += weights[k]*
val*values[k][i] /
basis->norm_squared(i);
283 return std::log(a+std::sqrt(a*a+1.0));
288 return std::log(a+std::sqrt(a*a-1.0));
293 return 0.5*std::log((1.0+a)/(1.0-a));
298 double operator() (
double a,
double b)
const {
return a + b; }
301 double operator() (
double a,
double b)
const {
return a - b; }
304 double operator() (
double a,
double b)
const {
return a * b; }
307 double operator() (
double a,
double b)
const {
return a / b; }
310 double operator() (
double a,
double b)
const {
return std::pow(a,b); }
700 setup.exp->plus(ru, v, w);
816 setup.exp->minus(ru, v, w);
932 setup.exp->times(ru, v, w);
1048 setup.exp->divide(ru, v, w);
1233 setup.exp->plusEqual(ru, v);
1286 setup.exp->minusEqual(ru, v);
1287 setup.exp->unaryMinus(v, v);
1427 Teuchos::GlobalMPISession mpiSession(&argc, &
argv);
1428 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...
Othogonal polynomial expansions based on derivative calculations.
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_DerivExpansion, 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
void computePCE2LC(Stokhos::OrthogPolyApprox< OrdinalType, ValueType > &c, ValueType a, const Stokhos::OrthogPolyApprox< OrdinalType, ValueType > &b)
Teuchos::RCP< Stokhos::Dense3Tensor< int, double > > Dijk
Stokhos::OrthogPolyApprox< OrdinalType, ValueType > u
Teuchos::RCP< const Stokhos::Quadrature< OrdinalType, ValueType > > quad
Stokhos::OrthogPolyApprox< OrdinalType, ValueType > cu2
Stokhos::OrthogPolyApprox< OrdinalType, ValueType > cx
Teuchos::RCP< Stokhos::Sparse3Tensor< int, double > > Cijk
Stokhos::OrthogPolyApprox< OrdinalType, ValueType > cu
Stokhos::OrthogPolyApprox< OrdinalType, ValueType > sx
Teuchos::RCP< Stokhos::DerivOrthogPolyExpansion< OrdinalType, ValueType > > exp
Stokhos::OrthogPolyApprox< OrdinalType, ValueType > u2
Stokhos::OrthogPolyApprox< OrdinalType, ValueType > su2
Stokhos::OrthogPolyApprox< OrdinalType, ValueType > x
Stokhos::OrthogPolyApprox< OrdinalType, ValueType > su
void computePCE2(Stokhos::OrthogPolyApprox< OrdinalType, ValueType > &c, const Stokhos::OrthogPolyApprox< OrdinalType, ValueType > &a, const Stokhos::OrthogPolyApprox< OrdinalType, ValueType > &b)
Stokhos::OrthogPolyApprox< OrdinalType, ValueType > y
void computePCE1(Stokhos::OrthogPolyApprox< OrdinalType, ValueType > &c, const Stokhos::OrthogPolyApprox< OrdinalType, ValueType > &a)
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, double > > Bij
Teuchos::RCP< const Stokhos::CompletePolynomialBasis< OrdinalType, ValueType > > basis
void computePCE2RC(Stokhos::OrthogPolyApprox< OrdinalType, ValueType > &c, const Stokhos::OrthogPolyApprox< OrdinalType, ValueType > &a, ValueType b)