Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
Stokhos_PecosOneDOrthogPolyBasisImp.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
44template <typename ordinal_type, typename value_type>
45Stokhos::PecosOneDOrthogPolyBasis<ordinal_type, value_type>::
46PecosOneDOrthogPolyBasis(
47 const Teuchos::RCP<Pecos::OrthogonalPolynomial>& pecosPoly_,
48 const std::string& name_, ordinal_type p_) :
49 pecosPoly(pecosPoly_),
50 name(name_),
51 p(p_),
52 sparse_grid_growth_rule(webbur::level_to_order_linear_wn),
53 norms(p+1, value_type(0.0))
54{
55 for (ordinal_type i=0; i<=p; i++)
56 norms[i] = pecosPoly->norm_squared(i);
57}
58
59template <typename ordinal_type, typename value_type>
60Stokhos::PecosOneDOrthogPolyBasis<ordinal_type, value_type>::
61PecosOneDOrthogPolyBasis(
62 ordinal_type p_, const PecosOneDOrthogPolyBasis& basis) :
63 pecosPoly(basis.pecosPoly),
64 name(basis.name),
65 p(p_),
66 sparse_grid_growth_rule(basis.sparse_grid_growth_rule),
67 norms(p+1, value_type(0.0))
68{
69 for (ordinal_type i=0; i<=p; i++)
70 norms[i] = pecosPoly->norm_squared(i);
71}
72
73template <typename ordinal_type, typename value_type>
74Stokhos::PecosOneDOrthogPolyBasis<ordinal_type, value_type>::
75~PecosOneDOrthogPolyBasis()
76{
77}
78
79template <typename ordinal_type, typename value_type>
81Stokhos::PecosOneDOrthogPolyBasis<ordinal_type, value_type>::
82order() const
83{
84 return p;
85}
86
87template <typename ordinal_type, typename value_type>
89Stokhos::PecosOneDOrthogPolyBasis<ordinal_type, value_type>::
90size() const
91{
92 return p+1;
93}
94
95template <typename ordinal_type, typename value_type>
96const Teuchos::Array<value_type>&
97Stokhos::PecosOneDOrthogPolyBasis<ordinal_type, value_type>::
98norm_squared() const
99{
100 return norms;
101}
102
103template <typename ordinal_type, typename value_type>
104const value_type&
105Stokhos::PecosOneDOrthogPolyBasis<ordinal_type, value_type>::
106norm_squared(ordinal_type i) const
107{
108 return norms[i];
109}
110
111template <typename ordinal_type, typename value_type>
112Teuchos::RCP< Stokhos::Dense3Tensor<ordinal_type, value_type> >
113Stokhos::PecosOneDOrthogPolyBasis<ordinal_type, value_type>::
114computeTripleProductTensor() const
115{
116 // Compute Cijk = < \Psi_i \Psi_j \Psi_k >
117 ordinal_type sz = size();
118 Teuchos::RCP< Stokhos::Dense3Tensor<ordinal_type, value_type> > Cijk =
119 Teuchos::rcp(new Dense3Tensor<ordinal_type, value_type>(sz));
120 Teuchos::Array<value_type> points, weights;
121 Teuchos::Array< Teuchos::Array<value_type> > values;
122 getQuadPoints(3*p, points, weights, values);
123
124 for (ordinal_type i=0; i<sz; i++) {
125 for (ordinal_type j=0; j<sz; j++) {
126 for (ordinal_type k=0; k<sz; k++) {
127 value_type triple_product = 0;
128 for (ordinal_type l=0; l<static_cast<ordinal_type>(points.size());
129 l++){
130 triple_product +=
131 weights[l]*(values[l][i])*(values[l][j])*(values[l][k]);
132 }
133 (*Cijk)(i,j,k) = triple_product;
134 }
135 }
136 }
137
138 return Cijk;
139}
140
141template <typename ordinal_type, typename value_type>
142Teuchos::RCP< Stokhos::Sparse3Tensor<ordinal_type, value_type> >
143Stokhos::PecosOneDOrthogPolyBasis<ordinal_type, value_type>::
144computeSparseTripleProductTensor(ordinal_type order) const
145{
146 // Compute Cijk = < \Psi_i \Psi_j \Psi_k >
147 value_type sparse_tol = 1.0e-15;
148 ordinal_type sz = size();
149 Teuchos::RCP< Stokhos::Sparse3Tensor<ordinal_type, value_type> > Cijk =
150 Teuchos::rcp(new Sparse3Tensor<ordinal_type, value_type>());
151 Teuchos::Array<value_type> points, weights;
152 Teuchos::Array< Teuchos::Array<value_type> > values;
153 getQuadPoints(3*p, points, weights, values);
154
155 for (ordinal_type i=0; i<sz; i++) {
156 for (ordinal_type j=0; j<sz; j++) {
157 for (ordinal_type k=0; k<order; k++) {
158 value_type triple_product = 0;
159 for (ordinal_type l=0; l<static_cast<ordinal_type>(points.size());
160 l++){
161 triple_product +=
162 weights[l]*(values[l][i])*(values[l][j])*(values[l][k]);
163 }
164 if (std::abs(triple_product/norms[i]) > sparse_tol)
165 Cijk->add_term(i,j,k,triple_product);
166 }
167 }
168 }
169 Cijk->fillComplete();
170
171 return Cijk;
172}
173
174template <typename ordinal_type, typename value_type>
175Teuchos::RCP< Teuchos::SerialDenseMatrix<ordinal_type, value_type> >
176Stokhos::PecosOneDOrthogPolyBasis<ordinal_type, value_type>::
177computeDerivDoubleProductTensor() const
178{
179 // Compute Bij = < \Psi_i' \Psi_j >
180 Teuchos::Array<value_type> points, weights;
181 Teuchos::Array< Teuchos::Array<value_type> > values, derivs;
182 getQuadPoints(2*p, points, weights, values);
183 ordinal_type nqp = weights.size();
184 derivs.resize(nqp);
185 ordinal_type sz = size();
186 for (ordinal_type i=0; i<nqp; i++) {
187 derivs[i].resize(sz);
188 evaluateBasesAndDerivatives(points[i], values[i], derivs[i]);
189 }
190 Teuchos::RCP< Teuchos::SerialDenseMatrix<ordinal_type, value_type> > Bij =
191 Teuchos::rcp(new Teuchos::SerialDenseMatrix<ordinal_type, value_type>(sz,sz));
192 for (ordinal_type i=0; i<sz; i++) {
193 for (ordinal_type j=0; j<sz; j++) {
194 value_type b = value_type(0.0);
195 for (int qp=0; qp<nqp; qp++)
196 b += weights[qp]*derivs[qp][i]*values[qp][j];
197 (*Bij)(i,j) = b;
198 }
199 }
200
201 return Bij;
202}
203
204template <typename ordinal_type, typename value_type>
205void
206Stokhos::PecosOneDOrthogPolyBasis<ordinal_type,value_type>::
207evaluateBases(const value_type& x, Teuchos::Array<value_type>& basis_pts) const
208{
209 for (ordinal_type i=0; i<=p; i++)
210 basis_pts[i] = pecosPoly->type1_value(x, i);
211}
212
213template <typename ordinal_type, typename value_type>
214void
215Stokhos::PecosOneDOrthogPolyBasis<ordinal_type,value_type>::
216evaluateBasesAndDerivatives(const value_type& x,
217 Teuchos::Array<value_type>& vals,
218 Teuchos::Array<value_type>& derivs) const
219{
220 for (ordinal_type i=0; i<=p; i++) {
221 vals[i] = pecosPoly->type1_value(x, i);
222 derivs[i] = pecosPoly->type1_gradient(x, i);
223 }
224}
225
226template <typename ordinal_type, typename value_type>
228Stokhos::PecosOneDOrthogPolyBasis<ordinal_type,value_type>::
229evaluate(const value_type& x, ordinal_type k) const
230{
231 return pecosPoly->type1_value(x, k);
232}
233
234template <typename ordinal_type, typename value_type>
235void
236Stokhos::PecosOneDOrthogPolyBasis<ordinal_type, value_type>::
237print(std::ostream& os) const
238{
239 os << "Pecos " << name << " basis of order " << p << "." << std::endl;
240 os << "Basis polynomial norms (squared):\n\t";
241 for (ordinal_type i=0; i<=p; i++)
242 os << norms[i] << " ";
243 os << std::endl;
244}
245
246template <typename ordinal_type, typename value_type>
247const std::string&
248Stokhos::PecosOneDOrthogPolyBasis<ordinal_type, value_type>::
249getName() const
250{
251 return name;
252}
253
254template <typename ordinal_type, typename value_type>
255void
256Stokhos::PecosOneDOrthogPolyBasis<ordinal_type,value_type>::
257getQuadPoints(ordinal_type quad_order,
258 Teuchos::Array<value_type>& quad_points,
259 Teuchos::Array<value_type>& quad_weights,
260 Teuchos::Array< Teuchos::Array<value_type> >& quad_values) const
261{
262 // This assumes Gaussian quadrature
263 ordinal_type num_points =
264 static_cast<ordinal_type>(std::ceil((quad_order+1)/2.0));
265 const Pecos::RealArray& gp = pecosPoly->collocation_points(num_points);
266 const Pecos::RealArray& gw = pecosPoly->type1_collocation_weights(num_points);
267 quad_points.resize(num_points);
268 quad_weights.resize(num_points);
269 for (ordinal_type i=0; i<num_points; i++) {
270 quad_points[i] = gp[i];
271 quad_weights[i] = gw[i];
272 }
273
274 // Evalute basis at gauss points
275 quad_values.resize(num_points);
276 for (ordinal_type i=0; i<num_points; i++) {
277 quad_values[i].resize(p+1);
278 evaluateBases(quad_points[i], quad_values[i]);
279 }
280}
281
282template <typename ordinal_type, typename value_type>
284Stokhos::PecosOneDOrthogPolyBasis<ordinal_type,value_type>::
285quadDegreeOfExactness(ordinal_type n) const
286{
287 return ordinal_type(2)*n-ordinal_type(1);
288}
289
290template <typename ordinal_type, typename value_type>
291Teuchos::RCP<Stokhos::OneDOrthogPolyBasis<ordinal_type,value_type> >
292Stokhos::PecosOneDOrthogPolyBasis<ordinal_type,value_type>::
293cloneWithOrder(ordinal_type p) const
294{
295 return
296 Teuchos::rcp(new Stokhos::PecosOneDOrthogPolyBasis<ordinal_type,value_type>(p,*this));
297}
298
299template <typename ordinal_type, typename value_type>
301Stokhos::PecosOneDOrthogPolyBasis<ordinal_type,value_type>::
302coefficientGrowth(ordinal_type n) const
303{
304 return n;
305}
306
307template <typename ordinal_type, typename value_type>
309Stokhos::PecosOneDOrthogPolyBasis<ordinal_type,value_type>::
310pointGrowth(ordinal_type n) const
311{
312 if (n % ordinal_type(2) == ordinal_type(1))
313 return n + ordinal_type(1);
314 return n;
315}