Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
Stokhos_AnisoSparseGridQuadratureImp.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 "sandia_sgmga.hpp"
45#include "sandia_rules.hpp"
46#include "Teuchos_TimeMonitor.hpp"
47
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())
56{
57#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
58 TEUCHOS_FUNC_TIME_MONITOR("Stokhos::AnisoSparseGridQuadrature -- Quad Grid Generation");
59#endif
60 ordinal_type d = product_basis->dimension();
61 ordinal_type p = product_basis->order();
62 ordinal_type sz = product_basis->size();
63 ordinal_type level = sparse_grid_level;
64
65 // Mike's heuristic formula for computing the level
66 if (level == 0) {
67 level = static_cast<ordinal_type>(std::ceil(0.5*(p+d-1)));
68 if (level < d)
69 level = p;
70 }
71
72 std::cout << "Sparse grid level = " << level << std::endl;
73
74 // Compute quad points, weights, values
75 Teuchos::Array<typename OneDOrthogPolyBasis<ordinal_type,value_type>::LevelToOrderFnPtr> growth_rules(d);
76
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();
83 }
84
85 // Set the static sparse grid quadrature pointer to this
86 // (this will cause a conflict if another sparse grid quadrature object
87 // is trying to access the VPISparseGrid library, but that's all we can
88 // do at this point).
89 sgq = this;
90
91 int num_total_pts =
92 webbur::sandia_sgmga_size_total(d,&dim_weights[0], level, growth_rate,
93 &growth_rules[0]);
94
95 ordinal_type ntot =
96 webbur::sandia_sgmga_size(d,&dim_weights[0],level,
97 &compute1DPoints[0], duplicate_tol, growth_rate,
98 &growth_rules[0]);
99
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);
107
108 webbur::sandia_sgmga_unique_index(d, &dim_weights[0], level,
109 &compute1DPoints[0],
110 duplicate_tol, ntot, num_total_pts,
111 growth_rate, &growth_rules[0],
112 &sparse_unique_index[0] );
113
114
115 webbur::sandia_sgmga_index(d, &dim_weights[0], level,
116 ntot, num_total_pts,
117 &sparse_unique_index[0],
118 growth_rate, &growth_rules[0],
119 &sparse_order[0], &sparse_index[0]);
120
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],
125 &quad_weights[0]);
126
127 webbur::sandia_sgmga_point(d, &dim_weights[0], level,
128 &compute1DPoints[0],
129 ntot, &sparse_order[0], &sparse_index[0],
130 growth_rate, &growth_rules[0],
131 &gp[0]);
132
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]);
139 }
140
141 std::cout << "Number of quadrature points = " << ntot << std::endl;
142}
143
144template <typename ordinal_type, typename value_type>
145const Teuchos::Array< Teuchos::Array<value_type> >&
146Stokhos::AnisoSparseGridQuadrature<ordinal_type, value_type>::
147getQuadPoints() const
148{
149 return quad_points;
150}
151
152template <typename ordinal_type, typename value_type>
153const Teuchos::Array<value_type>&
154Stokhos::AnisoSparseGridQuadrature<ordinal_type, value_type>::
155getQuadWeights() const
156{
157 return quad_weights;
158}
159
160template <typename ordinal_type, typename value_type>
161const Teuchos::Array< Teuchos::Array<value_type> >&
162Stokhos::AnisoSparseGridQuadrature<ordinal_type, value_type>::
163getBasisAtQuadPoints() const
164{
165 return quad_values;
166}
167
168template <typename ordinal_type, typename value_type>
169void
170Stokhos::AnisoSparseGridQuadrature<ordinal_type,value_type>::
171getMyPoints( int order, int dim, double x[] )
172{
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];
182 }
183}
184
185template <typename ordinal_type, typename value_type>
186void
187Stokhos::AnisoSparseGridQuadrature<ordinal_type,value_type>::
188getMyWeights( int order, int dim, double w[] )
189{
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];
199 }
200}
const IndexType const IndexType const IndexType const IndexType const ValueType const ValueType * x
Definition: csr_vector.h:265