42 #include "Teuchos_Assert.hpp" 43 #include "Teuchos_TimeMonitor.hpp" 45 template <
typename ordinal_type,
typename value_type,
typename func_type>
49 const Teuchos::RCP<const func_type>& func_,
51 bool use_pce_quad_points_,
56 pce_weights(quad->getQuadWeights()),
57 basis_values(quad->getBasisAtQuadPoints()),
60 use_pce_quad_points(use_pce_quad_points_)
66 template <
typename ordinal_type,
typename value_type,
typename func_type>
72 template <
typename ordinal_type,
typename value_type,
typename func_type>
76 Teuchos::Array<value_type>& quad_points,
77 Teuchos::Array<value_type>& quad_weights,
78 Teuchos::Array< Teuchos::Array<value_type> >& quad_values)
const 80 #ifdef STOKHOS_TEUCHOS_TIME_MONITOR 81 TEUCHOS_FUNC_TIME_MONITOR(
"Stokhos::StieltjesBasis -- compute Gauss points");
85 if (use_pce_quad_points) {
86 quad_points = func_vals;
87 quad_weights = pce_weights;
88 quad_values = phi_vals;
94 static_cast<ordinal_type>(std::ceil((quad_order+1)/2.0));
97 if (quad_order > 2*this->p)
98 quad_order = 2*this->p;
105 if (quad_weights.size() < num_points) {
107 quad_weights.resize(num_points);
108 quad_points.resize(num_points);
109 quad_values.resize(num_points);
112 quad_points[i] = quad_points[0];
113 quad_values[i].resize(this->p+1);
114 this->evaluateBases(quad_points[i], quad_values[i]);
119 template <
typename ordinal_type,
typename value_type,
typename func_type>
123 Teuchos::Array<value_type>& alpha,
124 Teuchos::Array<value_type>& beta,
125 Teuchos::Array<value_type>& delta,
126 Teuchos::Array<value_type>& gamma)
const 129 Teuchos::Array<value_type> nrm(n);
130 Teuchos::Array< Teuchos::Array<value_type> > vals(nqp);
133 stieltjes(0, n, pce_weights, func_vals, alpha, beta, nrm, vals);
146 template <
typename ordinal_type,
typename value_type,
typename func_type>
152 const Teuchos::Array< Teuchos::Array<value_type> >& quad_points =
153 quad->getQuadPoints();
155 func_vals.resize(nqp);
156 phi_vals.resize(nqp);
158 func_vals[i] = (*func)(quad_points[i]);
159 phi_vals[i].resize(this->p+1);
165 template <
typename ordinal_type,
typename value_type,
typename func_type>
170 const Teuchos::Array<value_type>& weights,
171 const Teuchos::Array<value_type>& points,
172 Teuchos::Array<value_type>& a,
173 Teuchos::Array<value_type>& b,
174 Teuchos::Array<value_type>& nrm,
175 Teuchos::Array< Teuchos::Array<value_type> >& phi_vals)
const 177 #ifdef STOKHOS_TEUCHOS_TIME_MONITOR 178 TEUCHOS_FUNC_TIME_MONITOR(
"Stokhos::StieltjesBasis -- Discretized Stieltjes Procedure");
184 integrateBasisSquared(0, a, b, weights, points, phi_vals, val1, val2);
191 integrateBasisSquared(i, a, b, weights, points, phi_vals, val1, val2);
194 TEUCHOS_TEST_FOR_EXCEPTION(val1 < 0.0, std::logic_error,
195 "Stokhos::StieltjesBasis::stieltjes(): " 196 <<
" Polynomial " << i <<
" out of " << nfinish
197 <<
" has norm " << val1
198 <<
"! Try increasing number of quadrature points");
201 b[i] = nrm[i]/nrm[i-1];
205 template <
typename ordinal_type,
typename value_type,
typename func_type>
209 const Teuchos::Array<value_type>& a,
210 const Teuchos::Array<value_type>& b,
211 const Teuchos::Array<value_type>& weights,
212 const Teuchos::Array<value_type>& points,
213 Teuchos::Array< Teuchos::Array<value_type> >& phi_vals,
216 evaluateRecurrence(k, a, b, points, phi_vals);
221 val1 += weights[i]*phi_vals[i][k]*phi_vals[i][k];
222 val2 += weights[i]*phi_vals[i][k]*phi_vals[i][k]*points[i];
226 template <
typename ordinal_type,
typename value_type,
typename func_type>
230 const Teuchos::Array<value_type>& a,
231 const Teuchos::Array<value_type>& b,
232 const Teuchos::Array<value_type>& points,
233 Teuchos::Array< Teuchos::Array<value_type> >& values)
const 241 values[i][k] = points[i] - a[k-1];
245 (points[i] - a[k-1])*values[i][k-1] - b[k-1]*values[i][k-2];
248 template <
typename ordinal_type,
typename value_type,
typename func_type>
249 Teuchos::RCP<Stokhos::OneDOrthogPolyBasis<ordinal_type,value_type> >
256 template <
typename ordinal_type,
typename value_type,
typename func_type>
262 pce_weights(quad->getQuadWeights()),
263 basis_values(quad->getBasisAtQuadPoints()),
266 use_pce_quad_points(basis.use_pce_quad_points)
StieltjesBasis(ordinal_type p, const Teuchos::RCP< const func_type > &func, const Teuchos::RCP< const Stokhos::Quadrature< ordinal_type, value_type > > &quad, bool use_pce_quad_points, bool normalize=false)
Constructor.
virtual void getQuadPoints(ordinal_type quad_order, Teuchos::Array< value_type > &points, Teuchos::Array< value_type > &weights, Teuchos::Array< Teuchos::Array< value_type > > &values) const
Get Gauss quadrature points, weights, and values of basis at points.
Implementation of OneDOrthogPolyBasis based on the general three-term recurrence relationship: for ...
virtual void getQuadPoints(ordinal_type quad_order, Teuchos::Array< value_type > &points, Teuchos::Array< value_type > &weights, Teuchos::Array< Teuchos::Array< value_type > > &values) const
Compute quadrature points, weights, and values of basis polynomials at given set of points points...
Generates three-term recurrence using the Discretized Stieltjes procedure applied to a functional map...
virtual void setup()
Setup basis after computing recurrence coefficients.
~StieltjesBasis()
Destructor.
Abstract base class for quadrature methods.
virtual Teuchos::RCP< OneDOrthogPolyBasis< ordinal_type, value_type > > cloneWithOrder(ordinal_type p) const
Clone this object with the option of building a higher order basis.
virtual void setup()
Setup basis after computing recurrence coefficients.
void integrateBasisSquared(ordinal_type k, const Teuchos::Array< value_type > &a, const Teuchos::Array< value_type > &b, const Teuchos::Array< value_type > &weights, const Teuchos::Array< value_type > &points, Teuchos::Array< Teuchos::Array< value_type > > &phi_vals, value_type &val1, value_type &val2) const
Compute and .
void stieltjes(ordinal_type nstart, ordinal_type nfinish, const Teuchos::Array< value_type > &weights, const Teuchos::Array< value_type > &points, Teuchos::Array< value_type > &a, Teuchos::Array< value_type > &b, Teuchos::Array< value_type > &nrm, Teuchos::Array< Teuchos::Array< value_type > > &phi_vals) const
Compute 3-term recurrence using Stieljtes procedure.
virtual bool computeRecurrenceCoefficients(ordinal_type n, Teuchos::Array< value_type > &alpha, Teuchos::Array< value_type > &beta, Teuchos::Array< value_type > &delta, Teuchos::Array< value_type > &gamma) const
Compute recurrence coefficients.
void evaluateRecurrence(ordinal_type k, const Teuchos::Array< value_type > &a, const Teuchos::Array< value_type > &b, const Teuchos::Array< value_type > &points, Teuchos::Array< Teuchos::Array< value_type > > &values) const
Evaluate polynomials via 3-term recurrence.