Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
Stokhos_SparseGridQuadratureImp.hpp
Go to the documentation of this file.
1// $Id$
2// $Source$
3// @HEADER
4// ***********************************************************************
5//
6// Stokhos Package
7// Copyright (2009) Sandia Corporation
8//
9// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
10// license for use of this work by or on behalf of the U.S. Government.
11//
12// Redistribution and use in source and binary forms, with or without
13// modification, are permitted provided that the following conditions are
14// met:
15//
16// 1. Redistributions of source code must retain the above copyright
17// notice, this list of conditions and the following disclaimer.
18//
19// 2. Redistributions in binary form must reproduce the above copyright
20// notice, this list of conditions and the following disclaimer in the
21// documentation and/or other materials provided with the distribution.
22//
23// 3. Neither the name of the Corporation nor the names of the
24// contributors may be used to endorse or promote products derived from
25// this software without specific prior written permission.
26//
27// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
28// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
29// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
30// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
31// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
32// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
33// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
34// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
35// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
36// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
37// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
38//
39// Questions? Contact Eric T. Phipps (etphipp@sandia.gov).
40//
41// ***********************************************************************
42// @HEADER
43
44#include "TriKota_ConfigDefs.hpp"
45#include "sandia_sgmg.hpp"
46#include "sandia_rules.hpp"
47#include "Teuchos_TimeMonitor.hpp"
48
49template <typename ordinal_type, typename value_type>
50Stokhos::SparseGridQuadrature<ordinal_type, value_type>::
51SparseGridQuadrature(
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())
57{
58#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
59 TEUCHOS_FUNC_TIME_MONITOR("Stokhos: Sparse Grid Generation");
60#endif
61
62 ordinal_type d = product_basis->dimension();
63 ordinal_type p = product_basis->order();
64 ordinal_type sz = product_basis->size();
65 ordinal_type level = sparse_grid_level;
66
67 // Make level = order the default, which is correct for Gaussian abscissas
68 // slow, linear growth, and total-order basis
69 if (level == 0) {
70 level = p;
71 }
72
73 //std::cout << "Sparse grid level = " << level << std::endl;
74
75 // Compute quad points, weights, values
76 Teuchos::Array<typename OneDOrthogPolyBasis<ordinal_type,value_type>::LevelToOrderFnPtr> growth_rules(d);
77
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();
84 }
85
86 // Set the static sparse grid quadrature pointer to this
87 // (this will cause a conflict if another sparse grid quadrature object
88 // is trying to access the VPISparseGrid library, but that's all we can
89 // do at this point).
90 sgq = this;
91
92 int num_total_pts =
93 webbur::sgmg_size_total(d, level, growth_rate, &growth_rules[0]);
94 ordinal_type ntot =
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);
104
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,
111 &growth_rules[0],
112 &sparse_order[0], &sparse_index[0]);
113 webbur::sgmg_weight(d, level, &compute1DWeights[0],
114 ntot, num_total_pts,
115 &sparse_unique_index[0], growth_rate,
116 &growth_rules[0],
117 &quad_weights[0]);
118 webbur::sgmg_point(d, level, &compute1DPoints[0],
119 ntot, &sparse_order[0],
120 &sparse_index[0], growth_rate,
121 &growth_rules[0],
122 &gp[0]);
123
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]);
130 }
131
132 //std::cout << "Number of quadrature points = " << ntot << std::endl;
133}
134
135template <typename ordinal_type, typename value_type>
136const Teuchos::Array< Teuchos::Array<value_type> >&
137Stokhos::SparseGridQuadrature<ordinal_type, value_type>::
138getQuadPoints() const
139{
140 return quad_points;
141}
142
143template <typename ordinal_type, typename value_type>
144const Teuchos::Array<value_type>&
145Stokhos::SparseGridQuadrature<ordinal_type, value_type>::
146getQuadWeights() const
147{
148 return quad_weights;
149}
150
151template <typename ordinal_type, typename value_type>
152const Teuchos::Array< Teuchos::Array<value_type> >&
153Stokhos::SparseGridQuadrature<ordinal_type, value_type>::
154getBasisAtQuadPoints() const
155{
156 return quad_values;
157}
158
159template <typename ordinal_type, typename value_type>
160void
161Stokhos::SparseGridQuadrature<ordinal_type,value_type>::
162getMyPoints( int order, int dim, double x[] )
163{
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];
173 }
174}
175
176template <typename ordinal_type, typename value_type>
177void
178Stokhos::SparseGridQuadrature<ordinal_type,value_type>::
179getMyWeights( int order, int dim, double w[] )
180{
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];
190 }
191}
192
193template <typename ordinal_type, typename value_type>
194std::ostream&
195Stokhos::SparseGridQuadrature<ordinal_type,value_type>::
196print(std::ostream& os) const
197{
198 ordinal_type nqp = quad_weights.size();
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());
204 j++)
205 os << quad_points[i][j] << " ";
206 os << std::endl;
207 }
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());
212 j++)
213 os << quad_values[i][j] << " ";
214 os << std::endl;
215 }
216
217 return os;
218}
const IndexType const IndexType const IndexType const IndexType const ValueType const ValueType * x
Definition: csr_vector.h:265