44#include "Teuchos_UnitTestHarness.hpp"
45#include "Teuchos_TestingHelpers.hpp"
46#include "Teuchos_UnitTestRepository.hpp"
47#include "Teuchos_GlobalMPISession.hpp"
69 mutable Teuchos::Array<double>
vec;
77 double operator() (
const double& a,
const double& b)
const {
84 mutable Teuchos::Array<double>
vec;
88template <
typename OrdinalType,
typename ValueType>
92 Teuchos::RCP<const Stokhos::CompletePolynomialBasis<OrdinalType,ValueType> >
basis;
93 Teuchos::RCP< Stokhos::QuadOrthogPolyExpansion<OrdinalType,ValueType> >
exp;
94 Teuchos::RCP<const Stokhos::GramSchmidtBasis<OrdinalType,ValueType> >
gs_basis;
95 Teuchos::RCP<const Stokhos::Quadrature<OrdinalType,ValueType> >
gs_quad;
102 const OrdinalType d = 3;
103 const OrdinalType p = 5;
106 Teuchos::Array< Teuchos::RCP<const Stokhos::OneDOrthogPolyBasis<OrdinalType,ValueType> > > bases(d);
107 for (OrdinalType i=0; i<d; i++)
116 for (OrdinalType i=0; i<d; i++)
120 Teuchos::RCP<const Stokhos::Quadrature<OrdinalType,ValueType> > quad =
124 Teuchos::RCP<Stokhos::Sparse3Tensor<int,double> > Cijk =
125 basis->computeTripleProductTensor();
128 Teuchos::RCP<Teuchos::ParameterList> exp_params =
129 Teuchos::rcp(
new Teuchos::ParameterList);
130 exp_params->set(
"Use Quadrature for Times",
true);
142 Teuchos::Array< Teuchos::RCP<const Stokhos::OneDOrthogPolyBasis<OrdinalType,ValueType> > > st_bases(2);
145 Teuchos::RCP<const Stokhos::OrthogPolyBasis<OrdinalType,ValueType> > st_basis =
148 u_st.term(0, 0) =
u.
mean();
149 u_st.term(0, 1) = 1.0;
151 v_st.
term(1, 1) = 1.0;
154 Teuchos::Array<ValueType> st_points_0;
155 Teuchos::Array<ValueType> st_weights_0;
156 Teuchos::Array< Teuchos::Array<ValueType> > st_values_0;
157 st_bases[0]->getQuadPoints(p+1, st_points_0, st_weights_0, st_values_0);
158 Teuchos::Array<ValueType> st_points_1;
159 Teuchos::Array<ValueType> st_weights_1;
160 Teuchos::Array< Teuchos::Array<ValueType> > st_values_1;
161 st_bases[1]->getQuadPoints(p+1, st_points_1, st_weights_1, st_values_1);
162 Teuchos::RCP< Teuchos::Array< Teuchos::Array<ValueType> > > st_points =
163 Teuchos::rcp(
new Teuchos::Array< Teuchos::Array<ValueType> >(st_points_0.size()));
164 for (
int i=0; i<st_points_0.size(); i++) {
165 (*st_points)[i].resize(2);
166 (*st_points)[i][0] = st_points_0[i];
167 (*st_points)[i][1] = st_points_1[i];
169 Teuchos::RCP< Teuchos::Array<ValueType> > st_weights =
170 Teuchos::rcp(
new Teuchos::Array<ValueType>);
171 *st_weights = st_weights_0;
177 Teuchos::RCP<Stokhos::Sparse3Tensor<int,double> > gs_Cijk =
178 gs_basis->computeTripleProductTensor();
182 Teuchos::RCP< const Teuchos::Array< Teuchos::Array<ValueType> > > points =
184 Teuchos::RCP< const Teuchos::Array<ValueType> > weights = st_weights;
215 success = comparePCEs(
setup.u,
"u", u2,
"u2",
setup.rtol,
setup.atol, out);
220 const Teuchos::Array<double>& norms =
setup.gs_basis->norm_squared();
221 const Teuchos::Array<double>& weights =
setup.gs_quad->getQuadWeights();
222 const Teuchos::Array< Teuchos::Array<double> >& values =
223 setup.gs_quad->getBasisAtQuadPoints();
224 Teuchos::SerialDenseMatrix<int,double> mat(
setup.gs_sz,
setup.gs_sz);
225 for (
int i=0; i<
setup.gs_sz; i++) {
227 for (
int k=0; k<weights.size(); k++)
228 mat(i,
j) += weights[k]*values[k][i]*values[k][
j];
229 mat(i,
j) /= std::sqrt(norms[i]*norms[
j]);
234 success = mat.normInf() < tol;
236 out <<
"\n Error, mat.normInf() < tol = " << mat.normInf()
237 <<
" < " << tol <<
": failed!\n";
238 out <<
"mat = " <<
printMat(mat) << std::endl;
247 success = comparePCEs(
setup.w,
"w", w2,
"w2",
setup.rtol,
setup.atol, out);
252 success = Teuchos::testRelErr(
"w.mean()",
setup.w.mean(),
253 "w_gs.mean()",
setup.w_gs.mean(),
256 Teuchos::Ptr<std::ostream>(out.getOStream().get()));
261 success = Teuchos::testRelErr(
"w.standard_deviation()",
262 setup.w.standard_deviation(),
263 "w_gs.standard_devaition()",
264 setup.w_gs.standard_deviation(),
267 Teuchos::Ptr<std::ostream>(out.getOStream().get()));
273 Teuchos::GlobalMPISession mpiSession(&argc, &
argv);
274 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...
Transforms a non-orthogonal multivariate basis to an orthogonal one using the Gram-Schmit procedure.
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.
value_type mean() const
Compute mean of expansion.
pointer coeff()
Return coefficient array.
reference term(ordinal_type dimension, ordinal_type order)
Get coefficient term for given dimension and order.
value_type evaluate(const Teuchos::Array< value_type > &point) const
Evaluate polynomial approximation at a point.
Abstract base class for multivariate orthogonal polynomials.
Orthogonal polynomial expansions based on numerical quadrature.
void times(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a, const OrthogPolyApprox< ordinal_type, value_type, node_type > &b)
Generates three-term recurrence using the Discretized Stieltjes procedure applied to a polynomial cha...
Defines quadrature for a tensor product basis by tensor products of 1-D quadrature rules.
TEUCHOS_UNIT_TEST(Stokhos_GramSchmidtBasis, Map)
GramSchmidt_PCE_Setup< int, double > setup
Stokhos::OrthogPolyApprox< OrdinalType, ValueType > v_gs
Stokhos::OrthogPolyApprox< OrdinalType, ValueType > u_gs
Teuchos::RCP< const Stokhos::CompletePolynomialBasis< OrdinalType, ValueType > > basis
Stokhos::OrthogPolyApprox< OrdinalType, ValueType > w_gs
Stokhos::OrthogPolyApprox< OrdinalType, ValueType > u
Teuchos::RCP< const Stokhos::Quadrature< OrdinalType, ValueType > > gs_quad
Stokhos::OrthogPolyApprox< OrdinalType, ValueType > w
Teuchos::RCP< Stokhos::QuadOrthogPolyExpansion< OrdinalType, ValueType > > exp
Teuchos::RCP< const Stokhos::GramSchmidtBasis< OrdinalType, ValueType > > gs_basis
Stokhos::OrthogPolyApprox< OrdinalType, ValueType > v
gram_schmidt_pce_binary_quad_func(const Stokhos::OrthogPolyApprox< int, double > &pce_, const Stokhos::OrthogPolyBasis< int, double > &basis_)
const Stokhos::OrthogPolyApprox< int, double > & pce
double operator()(const double &a, const double &b) const
const Stokhos::OrthogPolyBasis< int, double > & basis
Teuchos::Array< double > vec
double operator()(const double &a) const
gram_schmidt_pce_unary_quad_func(const Stokhos::OrthogPolyApprox< int, double > &pce_, const Stokhos::OrthogPolyBasis< int, double > &basis_)
Teuchos::Array< double > vec
const Stokhos::OrthogPolyApprox< int, double > & pce
const Stokhos::OrthogPolyBasis< int, double > & basis
void printMat(const char *name, Epetra_IntSerialDenseMatrix &matrix)