Stokhos Package Browser (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 
44 template <typename ordinal_type, typename value_type>
45 Stokhos::PecosOneDOrthogPolyBasis<ordinal_type, value_type>::
46 PecosOneDOrthogPolyBasis(
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 
59 template <typename ordinal_type, typename value_type>
60 Stokhos::PecosOneDOrthogPolyBasis<ordinal_type, value_type>::
61 PecosOneDOrthogPolyBasis(
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 
73 template <typename ordinal_type, typename value_type>
74 Stokhos::PecosOneDOrthogPolyBasis<ordinal_type, value_type>::
75 ~PecosOneDOrthogPolyBasis()
76 {
77 }
78 
79 template <typename ordinal_type, typename value_type>
80 ordinal_type
81 Stokhos::PecosOneDOrthogPolyBasis<ordinal_type, value_type>::
82 order() const
83 {
84  return p;
85 }
86 
87 template <typename ordinal_type, typename value_type>
88 ordinal_type
89 Stokhos::PecosOneDOrthogPolyBasis<ordinal_type, value_type>::
90 size() const
91 {
92  return p+1;
93 }
94 
95 template <typename ordinal_type, typename value_type>
96 const Teuchos::Array<value_type>&
97 Stokhos::PecosOneDOrthogPolyBasis<ordinal_type, value_type>::
98 norm_squared() const
99 {
100  return norms;
101 }
102 
103 template <typename ordinal_type, typename value_type>
104 const value_type&
105 Stokhos::PecosOneDOrthogPolyBasis<ordinal_type, value_type>::
106 norm_squared(ordinal_type i) const
107 {
108  return norms[i];
109 }
110 
111 template <typename ordinal_type, typename value_type>
112 Teuchos::RCP< Stokhos::Dense3Tensor<ordinal_type, value_type> >
113 Stokhos::PecosOneDOrthogPolyBasis<ordinal_type, value_type>::
114 computeTripleProductTensor() 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 
141 template <typename ordinal_type, typename value_type>
142 Teuchos::RCP< Stokhos::Sparse3Tensor<ordinal_type, value_type> >
143 Stokhos::PecosOneDOrthogPolyBasis<ordinal_type, value_type>::
144 computeSparseTripleProductTensor(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 
174 template <typename ordinal_type, typename value_type>
175 Teuchos::RCP< Teuchos::SerialDenseMatrix<ordinal_type, value_type> >
176 Stokhos::PecosOneDOrthogPolyBasis<ordinal_type, value_type>::
177 computeDerivDoubleProductTensor() 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 
204 template <typename ordinal_type, typename value_type>
205 void
206 Stokhos::PecosOneDOrthogPolyBasis<ordinal_type,value_type>::
207 evaluateBases(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 
213 template <typename ordinal_type, typename value_type>
214 void
215 Stokhos::PecosOneDOrthogPolyBasis<ordinal_type,value_type>::
216 evaluateBasesAndDerivatives(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 
226 template <typename ordinal_type, typename value_type>
227 value_type
228 Stokhos::PecosOneDOrthogPolyBasis<ordinal_type,value_type>::
229 evaluate(const value_type& x, ordinal_type k) const
230 {
231  return pecosPoly->type1_value(x, k);
232 }
233 
234 template <typename ordinal_type, typename value_type>
235 void
236 Stokhos::PecosOneDOrthogPolyBasis<ordinal_type, value_type>::
237 print(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 
246 template <typename ordinal_type, typename value_type>
247 const std::string&
248 Stokhos::PecosOneDOrthogPolyBasis<ordinal_type, value_type>::
249 getName() const
250 {
251  return name;
252 }
253 
254 template <typename ordinal_type, typename value_type>
255 void
256 Stokhos::PecosOneDOrthogPolyBasis<ordinal_type,value_type>::
257 getQuadPoints(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 
282 template <typename ordinal_type, typename value_type>
283 ordinal_type
284 Stokhos::PecosOneDOrthogPolyBasis<ordinal_type,value_type>::
285 quadDegreeOfExactness(ordinal_type n) const
286 {
287  return ordinal_type(2)*n-ordinal_type(1);
288 }
289 
290 template <typename ordinal_type, typename value_type>
291 Teuchos::RCP<Stokhos::OneDOrthogPolyBasis<ordinal_type,value_type> >
292 Stokhos::PecosOneDOrthogPolyBasis<ordinal_type,value_type>::
293 cloneWithOrder(ordinal_type p) const
294 {
295  return
296  Teuchos::rcp(new Stokhos::PecosOneDOrthogPolyBasis<ordinal_type,value_type>(p,*this));
297 }
298 
299 template <typename ordinal_type, typename value_type>
300 ordinal_type
301 Stokhos::PecosOneDOrthogPolyBasis<ordinal_type,value_type>::
302 coefficientGrowth(ordinal_type n) const
303 {
304  return n;
305 }
306 
307 template <typename ordinal_type, typename value_type>
308 ordinal_type
309 Stokhos::PecosOneDOrthogPolyBasis<ordinal_type,value_type>::
310 pointGrowth(ordinal_type n) const
311 {
312  if (n % ordinal_type(2) == ordinal_type(1))
313  return n + ordinal_type(1);
314  return n;
315 }
j
j
Definition: Sacado_Fad_Exp_MP_Vector.hpp:527
Teuchos::rcp
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Sacado::UQ::ceil
KOKKOS_INLINE_FUNCTION PCE< Storage > ceil(const PCE< Storage > &a)
Definition: Sacado_UQ_PCE_Imp.hpp:1191
Sacado::UQ::abs
KOKKOS_INLINE_FUNCTION PCE< Storage > abs(const PCE< Storage > &a)
Definition: Sacado_UQ_PCE_Imp.hpp:1181
Teuchos::Array::resize
void resize(size_type new_size, const value_type &x=value_type())
Teuchos::Array::size
size_type size() const
Generated by   doxygen 1.8.5