44#include "sandia_sgmga.hpp"
45#include "sandia_rules.hpp"
46#include "Teuchos_TimeMonitor.hpp"
48template <
typename ordinal_type,
typename value_type>
49Stokhos::AnisoSparseGridQuadrature<ordinal_type, value_type>::
50AnisoSparseGridQuadrature(
51 const Teuchos::RCP<
const ProductBasis<ordinal_type,value_type> >& product_basis,
52 ordinal_type sparse_grid_level, value_type dim_weights[],
53 value_type duplicate_tol,
54 ordinal_type growth_rate) :
55 coordinate_bases(product_basis->getCoordinateBases())
57#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
58 TEUCHOS_FUNC_TIME_MONITOR(
"Stokhos::AnisoSparseGridQuadrature -- Quad Grid Generation");
67 level =
static_cast<ordinal_type>(std::ceil(0.5*(p+d-1)));
72 std::cout <<
"Sparse grid level = " << level << std::endl;
75 Teuchos::Array<typename OneDOrthogPolyBasis<ordinal_type,value_type>::LevelToOrderFnPtr> growth_rules(d);
77 Teuchos::Array< void (*) (
int order,
int dim,
double x[] ) > compute1DPoints(d);
78 Teuchos::Array< void (*) (
int order,
int dim,
double w[] ) > compute1DWeights(d);
79 for (ordinal_type i=0; i<d; i++) {
80 compute1DPoints[i] = &(getMyPoints);
81 compute1DWeights[i] = &(getMyWeights);
82 growth_rules[i] = coordinate_bases[i]->getSparseGridGrowthRule();
92 webbur::sandia_sgmga_size_total(d,&dim_weights[0], level, growth_rate,
96 webbur::sandia_sgmga_size(d,&dim_weights[0],level,
97 &compute1DPoints[0], duplicate_tol, growth_rate,
100 Teuchos::Array<int> sparse_order(ntot*d);
101 Teuchos::Array<int> sparse_index(ntot*d);
102 Teuchos::Array<int> sparse_unique_index(num_total_pts);
103 quad_points.resize(ntot);
104 quad_weights.resize(ntot);
105 quad_values.resize(ntot);
106 Teuchos::Array<value_type> gp(ntot*d);
108 webbur::sandia_sgmga_unique_index(d, &dim_weights[0], level,
110 duplicate_tol, ntot, num_total_pts,
111 growth_rate, &growth_rules[0],
112 &sparse_unique_index[0] );
115 webbur::sandia_sgmga_index(d, &dim_weights[0], level,
117 &sparse_unique_index[0],
118 growth_rate, &growth_rules[0],
119 &sparse_order[0], &sparse_index[0]);
121 webbur::sandia_sgmga_weight(d,&dim_weights[0],level,
122 &compute1DWeights[0],
123 ntot, num_total_pts, &sparse_unique_index[0],
124 growth_rate, &growth_rules[0],
127 webbur::sandia_sgmga_point(d, &dim_weights[0], level,
129 ntot, &sparse_order[0], &sparse_index[0],
130 growth_rate, &growth_rules[0],
133 for (ordinal_type i=0; i<ntot; i++) {
134 quad_values[i].resize(sz);
135 quad_points[i].resize(d);
136 for (ordinal_type
j=0;
j<d;
j++)
137 quad_points[i][
j] = gp[i*d+
j];
138 product_basis->evaluateBases(quad_points[i], quad_values[i]);
141 std::cout <<
"Number of quadrature points = " << ntot << std::endl;
144template <
typename ordinal_type,
typename value_type>
145const Teuchos::Array< Teuchos::Array<value_type> >&
146Stokhos::AnisoSparseGridQuadrature<ordinal_type, value_type>::
152template <
typename ordinal_type,
typename value_type>
153const Teuchos::Array<value_type>&
154Stokhos::AnisoSparseGridQuadrature<ordinal_type, value_type>::
155getQuadWeights()
const
160template <
typename ordinal_type,
typename value_type>
161const Teuchos::Array< Teuchos::Array<value_type> >&
162Stokhos::AnisoSparseGridQuadrature<ordinal_type, value_type>::
163getBasisAtQuadPoints()
const
168template <
typename ordinal_type,
typename value_type>
170Stokhos::AnisoSparseGridQuadrature<ordinal_type,value_type>::
171getMyPoints(
int order,
int dim,
double x[] )
173 Teuchos::Array<double> quad_points;
174 Teuchos::Array<double> quad_weights;
175 Teuchos::Array< Teuchos::Array<double> > quad_values;
176 ordinal_type p = sgq->coordinate_bases[dim]->quadDegreeOfExactness(order);
177 sgq->coordinate_bases[dim]->getQuadPoints(p, quad_points,
178 quad_weights, quad_values);
179 TEUCHOS_ASSERT(quad_points.size() == order);
180 for (
int i = 0; i<quad_points.size(); i++) {
181 x[i] = quad_points[i];
185template <
typename ordinal_type,
typename value_type>
187Stokhos::AnisoSparseGridQuadrature<ordinal_type,value_type>::
188getMyWeights(
int order,
int dim,
double w[] )
190 Teuchos::Array<double> quad_points;
191 Teuchos::Array<double> quad_weights;
192 Teuchos::Array< Teuchos::Array<double> > quad_values;
193 ordinal_type p = sgq->coordinate_bases[dim]->quadDegreeOfExactness(order);
194 sgq->coordinate_bases[dim]->getQuadPoints(p, quad_points,
195 quad_weights, quad_values);
196 TEUCHOS_ASSERT(quad_points.size() == order);
197 for (
int i = 0; i<quad_points.size(); i++) {
198 w[i] = quad_weights[i];
const IndexType const IndexType const IndexType const IndexType const ValueType const ValueType * x