44#include "TriKota_ConfigDefs.hpp"
45#include "sandia_sgmg.hpp"
46#include "sandia_rules.hpp"
47#include "Teuchos_TimeMonitor.hpp"
49template <
typename ordinal_type,
typename value_type>
50Stokhos::SparseGridQuadrature<ordinal_type, value_type>::
52 const Teuchos::RCP<
const ProductBasis<ordinal_type,value_type> >& product_basis,
53 ordinal_type sparse_grid_level,
54 value_type duplicate_tol,
55 ordinal_type growth_rate) :
56 coordinate_bases(product_basis->getCoordinateBases())
58#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
59 TEUCHOS_FUNC_TIME_MONITOR(
"Stokhos: Sparse Grid Generation");
76 Teuchos::Array<typename OneDOrthogPolyBasis<ordinal_type,value_type>::LevelToOrderFnPtr> growth_rules(d);
78 Teuchos::Array< void (*) (
int order,
int dim,
double x[] ) > compute1DPoints(d);
79 Teuchos::Array< void (*) (
int order,
int dim,
double w[] ) > compute1DWeights(d);
80 for (ordinal_type i=0; i<d; i++) {
81 compute1DPoints[i] = &(getMyPoints);
82 compute1DWeights[i] = &(getMyWeights);
83 growth_rules[i] = coordinate_bases[i]->getSparseGridGrowthRule();
93 webbur::sgmg_size_total(d, level, growth_rate, &growth_rules[0]);
95 webbur::sgmg_size(d, level, &compute1DPoints[0], duplicate_tol,
96 growth_rate, &growth_rules[0]);
97 Teuchos::Array<int> sparse_order(ntot*d);
98 Teuchos::Array<int> sparse_index(ntot*d);
99 Teuchos::Array<int> sparse_unique_index(num_total_pts);
100 quad_points.resize(ntot);
101 quad_weights.resize(ntot);
102 quad_values.resize(ntot);
103 Teuchos::Array<value_type> gp(ntot*d);
105 webbur::sgmg_unique_index(d, level, &compute1DPoints[0],
106 duplicate_tol, ntot, num_total_pts,
107 growth_rate, &growth_rules[0],
108 &sparse_unique_index[0]);
109 webbur::sgmg_index(d, level, ntot, num_total_pts,
110 &sparse_unique_index[0], growth_rate,
112 &sparse_order[0], &sparse_index[0]);
113 webbur::sgmg_weight(d, level, &compute1DWeights[0],
115 &sparse_unique_index[0], growth_rate,
118 webbur::sgmg_point(d, level, &compute1DPoints[0],
119 ntot, &sparse_order[0],
120 &sparse_index[0], growth_rate,
124 for (ordinal_type i=0; i<ntot; i++) {
125 quad_values[i].resize(sz);
126 quad_points[i].resize(d);
127 for (ordinal_type
j=0;
j<d;
j++)
128 quad_points[i][
j] = gp[i*d+
j];
129 product_basis->evaluateBases(quad_points[i], quad_values[i]);
135template <
typename ordinal_type,
typename value_type>
136const Teuchos::Array< Teuchos::Array<value_type> >&
137Stokhos::SparseGridQuadrature<ordinal_type, value_type>::
143template <
typename ordinal_type,
typename value_type>
144const Teuchos::Array<value_type>&
145Stokhos::SparseGridQuadrature<ordinal_type, value_type>::
146getQuadWeights()
const
151template <
typename ordinal_type,
typename value_type>
152const Teuchos::Array< Teuchos::Array<value_type> >&
153Stokhos::SparseGridQuadrature<ordinal_type, value_type>::
154getBasisAtQuadPoints()
const
159template <
typename ordinal_type,
typename value_type>
161Stokhos::SparseGridQuadrature<ordinal_type,value_type>::
162getMyPoints(
int order,
int dim,
double x[] )
164 Teuchos::Array<double> quad_points;
165 Teuchos::Array<double> quad_weights;
166 Teuchos::Array< Teuchos::Array<double> > quad_values;
167 ordinal_type p = sgq->coordinate_bases[dim]->quadDegreeOfExactness(order);
168 sgq->coordinate_bases[dim]->getQuadPoints(p, quad_points,
169 quad_weights, quad_values);
170 TEUCHOS_ASSERT(quad_points.size() == order);
171 for (
int i = 0; i<quad_points.size(); i++) {
172 x[i] = quad_points[i];
176template <
typename ordinal_type,
typename value_type>
178Stokhos::SparseGridQuadrature<ordinal_type,value_type>::
179getMyWeights(
int order,
int dim,
double w[] )
181 Teuchos::Array<double> quad_points;
182 Teuchos::Array<double> quad_weights;
183 Teuchos::Array< Teuchos::Array<double> > quad_values;
184 ordinal_type p = sgq->coordinate_bases[dim]->quadDegreeOfExactness(order);
185 sgq->coordinate_bases[dim]->getQuadPoints(p, quad_points,
186 quad_weights, quad_values);
187 TEUCHOS_ASSERT(quad_points.size() == order);
188 for (
int i = 0; i<quad_points.size(); i++) {
189 w[i] = quad_weights[i];
193template <
typename ordinal_type,
typename value_type>
195Stokhos::SparseGridQuadrature<ordinal_type,value_type>::
196print(std::ostream& os)
const
199 os <<
"Sparse Grid Quadrature with " << nqp <<
" points:"
200 << std::endl <<
"Weight : Points" << std::endl;
201 for (ordinal_type i=0; i<nqp; i++) {
202 os << i <<
": " << quad_weights[i] <<
" : ";
203 for (ordinal_type
j=0; j<static_cast<ordinal_type>(quad_points[i].size());
205 os << quad_points[i][
j] <<
" ";
208 os <<
"Basis values at quadrature points:" << std::endl;
209 for (ordinal_type i=0; i<nqp; i++) {
210 os << i <<
" " <<
": ";
211 for (ordinal_type
j=0; j<static_cast<ordinal_type>(quad_values[i].size());
213 os << quad_values[i][
j] <<
" ";
const IndexType const IndexType const IndexType const IndexType const ValueType const ValueType * x