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 // @HEADER
2 // *****************************************************************************
3 // Stokhos Package
4 //
5 // Copyright 2009 NTESS and the Stokhos contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 // @HEADER
9 
10 template <typename ordinal_type, typename value_type>
11 Stokhos::PecosOneDOrthogPolyBasis<ordinal_type, value_type>::
12 PecosOneDOrthogPolyBasis(
13  const Teuchos::RCP<Pecos::OrthogonalPolynomial>& pecosPoly_,
14  const std::string& name_, ordinal_type p_) :
15  pecosPoly(pecosPoly_),
16  name(name_),
17  p(p_),
18  sparse_grid_growth_rule(webbur::level_to_order_linear_wn),
19  norms(p+1, value_type(0.0))
20 {
21  for (ordinal_type i=0; i<=p; i++)
22  norms[i] = pecosPoly->norm_squared(i);
23 }
24 
25 template <typename ordinal_type, typename value_type>
26 Stokhos::PecosOneDOrthogPolyBasis<ordinal_type, value_type>::
27 PecosOneDOrthogPolyBasis(
28  ordinal_type p_, const PecosOneDOrthogPolyBasis& basis) :
29  pecosPoly(basis.pecosPoly),
30  name(basis.name),
31  p(p_),
32  sparse_grid_growth_rule(basis.sparse_grid_growth_rule),
33  norms(p+1, value_type(0.0))
34 {
35  for (ordinal_type i=0; i<=p; i++)
36  norms[i] = pecosPoly->norm_squared(i);
37 }
38 
39 template <typename ordinal_type, typename value_type>
40 Stokhos::PecosOneDOrthogPolyBasis<ordinal_type, value_type>::
41 ~PecosOneDOrthogPolyBasis()
42 {
43 }
44 
45 template <typename ordinal_type, typename value_type>
46 ordinal_type
47 Stokhos::PecosOneDOrthogPolyBasis<ordinal_type, value_type>::
48 order() const
49 {
50  return p;
51 }
52 
53 template <typename ordinal_type, typename value_type>
54 ordinal_type
55 Stokhos::PecosOneDOrthogPolyBasis<ordinal_type, value_type>::
56 size() const
57 {
58  return p+1;
59 }
60 
61 template <typename ordinal_type, typename value_type>
62 const Teuchos::Array<value_type>&
63 Stokhos::PecosOneDOrthogPolyBasis<ordinal_type, value_type>::
64 norm_squared() const
65 {
66  return norms;
67 }
68 
69 template <typename ordinal_type, typename value_type>
70 const value_type&
71 Stokhos::PecosOneDOrthogPolyBasis<ordinal_type, value_type>::
72 norm_squared(ordinal_type i) const
73 {
74  return norms[i];
75 }
76 
77 template <typename ordinal_type, typename value_type>
78 Teuchos::RCP< Stokhos::Dense3Tensor<ordinal_type, value_type> >
79 Stokhos::PecosOneDOrthogPolyBasis<ordinal_type, value_type>::
80 computeTripleProductTensor() const
81 {
82  // Compute Cijk = < \Psi_i \Psi_j \Psi_k >
83  ordinal_type sz = size();
84  Teuchos::RCP< Stokhos::Dense3Tensor<ordinal_type, value_type> > Cijk =
85  Teuchos::rcp(new Dense3Tensor<ordinal_type, value_type>(sz));
86  Teuchos::Array<value_type> points, weights;
87  Teuchos::Array< Teuchos::Array<value_type> > values;
88  getQuadPoints(3*p, points, weights, values);
89 
90  for (ordinal_type i=0; i<sz; i++) {
91  for (ordinal_type j=0; j<sz; j++) {
92  for (ordinal_type k=0; k<sz; k++) {
93  value_type triple_product = 0;
94  for (ordinal_type l=0; l<static_cast<ordinal_type>(points.size());
95  l++){
96  triple_product +=
97  weights[l]*(values[l][i])*(values[l][j])*(values[l][k]);
98  }
99  (*Cijk)(i,j,k) = triple_product;
100  }
101  }
102  }
103 
104  return Cijk;
105 }
106 
107 template <typename ordinal_type, typename value_type>
108 Teuchos::RCP< Stokhos::Sparse3Tensor<ordinal_type, value_type> >
109 Stokhos::PecosOneDOrthogPolyBasis<ordinal_type, value_type>::
110 computeSparseTripleProductTensor(ordinal_type order) const
111 {
112  // Compute Cijk = < \Psi_i \Psi_j \Psi_k >
113  value_type sparse_tol = 1.0e-15;
114  ordinal_type sz = size();
115  Teuchos::RCP< Stokhos::Sparse3Tensor<ordinal_type, value_type> > Cijk =
116  Teuchos::rcp(new Sparse3Tensor<ordinal_type, value_type>());
117  Teuchos::Array<value_type> points, weights;
118  Teuchos::Array< Teuchos::Array<value_type> > values;
119  getQuadPoints(3*p, points, weights, values);
120 
121  for (ordinal_type i=0; i<sz; i++) {
122  for (ordinal_type j=0; j<sz; j++) {
123  for (ordinal_type k=0; k<order; k++) {
124  value_type triple_product = 0;
125  for (ordinal_type l=0; l<static_cast<ordinal_type>(points.size());
126  l++){
127  triple_product +=
128  weights[l]*(values[l][i])*(values[l][j])*(values[l][k]);
129  }
130  if (std::abs(triple_product/norms[i]) > sparse_tol)
131  Cijk->add_term(i,j,k,triple_product);
132  }
133  }
134  }
135  Cijk->fillComplete();
136 
137  return Cijk;
138 }
139 
140 template <typename ordinal_type, typename value_type>
141 Teuchos::RCP< Teuchos::SerialDenseMatrix<ordinal_type, value_type> >
142 Stokhos::PecosOneDOrthogPolyBasis<ordinal_type, value_type>::
143 computeDerivDoubleProductTensor() const
144 {
145  // Compute Bij = < \Psi_i' \Psi_j >
146  Teuchos::Array<value_type> points, weights;
147  Teuchos::Array< Teuchos::Array<value_type> > values, derivs;
148  getQuadPoints(2*p, points, weights, values);
149  ordinal_type nqp = weights.size();
150  derivs.resize(nqp);
151  ordinal_type sz = size();
152  for (ordinal_type i=0; i<nqp; i++) {
153  derivs[i].resize(sz);
154  evaluateBasesAndDerivatives(points[i], values[i], derivs[i]);
155  }
156  Teuchos::RCP< Teuchos::SerialDenseMatrix<ordinal_type, value_type> > Bij =
157  Teuchos::rcp(new Teuchos::SerialDenseMatrix<ordinal_type, value_type>(sz,sz));
158  for (ordinal_type i=0; i<sz; i++) {
159  for (ordinal_type j=0; j<sz; j++) {
160  value_type b = value_type(0.0);
161  for (int qp=0; qp<nqp; qp++)
162  b += weights[qp]*derivs[qp][i]*values[qp][j];
163  (*Bij)(i,j) = b;
164  }
165  }
166 
167  return Bij;
168 }
169 
170 template <typename ordinal_type, typename value_type>
171 void
172 Stokhos::PecosOneDOrthogPolyBasis<ordinal_type,value_type>::
173 evaluateBases(const value_type& x, Teuchos::Array<value_type>& basis_pts) const
174 {
175  for (ordinal_type i=0; i<=p; i++)
176  basis_pts[i] = pecosPoly->type1_value(x, i);
177 }
178 
179 template <typename ordinal_type, typename value_type>
180 void
181 Stokhos::PecosOneDOrthogPolyBasis<ordinal_type,value_type>::
182 evaluateBasesAndDerivatives(const value_type& x,
183  Teuchos::Array<value_type>& vals,
184  Teuchos::Array<value_type>& derivs) const
185 {
186  for (ordinal_type i=0; i<=p; i++) {
187  vals[i] = pecosPoly->type1_value(x, i);
188  derivs[i] = pecosPoly->type1_gradient(x, i);
189  }
190 }
191 
192 template <typename ordinal_type, typename value_type>
193 value_type
194 Stokhos::PecosOneDOrthogPolyBasis<ordinal_type,value_type>::
195 evaluate(const value_type& x, ordinal_type k) const
196 {
197  return pecosPoly->type1_value(x, k);
198 }
199 
200 template <typename ordinal_type, typename value_type>
201 void
202 Stokhos::PecosOneDOrthogPolyBasis<ordinal_type, value_type>::
203 print(std::ostream& os) const
204 {
205  os << "Pecos " << name << " basis of order " << p << "." << std::endl;
206  os << "Basis polynomial norms (squared):\n\t";
207  for (ordinal_type i=0; i<=p; i++)
208  os << norms[i] << " ";
209  os << std::endl;
210 }
211 
212 template <typename ordinal_type, typename value_type>
213 const std::string&
214 Stokhos::PecosOneDOrthogPolyBasis<ordinal_type, value_type>::
215 getName() const
216 {
217  return name;
218 }
219 
220 template <typename ordinal_type, typename value_type>
221 void
222 Stokhos::PecosOneDOrthogPolyBasis<ordinal_type,value_type>::
223 getQuadPoints(ordinal_type quad_order,
224  Teuchos::Array<value_type>& quad_points,
225  Teuchos::Array<value_type>& quad_weights,
226  Teuchos::Array< Teuchos::Array<value_type> >& quad_values) const
227 {
228  // This assumes Gaussian quadrature
229  ordinal_type num_points =
230  static_cast<ordinal_type>(std::ceil((quad_order+1)/2.0));
231  const Pecos::RealArray& gp = pecosPoly->collocation_points(num_points);
232  const Pecos::RealArray& gw = pecosPoly->type1_collocation_weights(num_points);
233  quad_points.resize(num_points);
234  quad_weights.resize(num_points);
235  for (ordinal_type i=0; i<num_points; i++) {
236  quad_points[i] = gp[i];
237  quad_weights[i] = gw[i];
238  }
239 
240  // Evalute basis at gauss points
241  quad_values.resize(num_points);
242  for (ordinal_type i=0; i<num_points; i++) {
243  quad_values[i].resize(p+1);
244  evaluateBases(quad_points[i], quad_values[i]);
245  }
246 }
247 
248 template <typename ordinal_type, typename value_type>
249 ordinal_type
250 Stokhos::PecosOneDOrthogPolyBasis<ordinal_type,value_type>::
251 quadDegreeOfExactness(ordinal_type n) const
252 {
253  return ordinal_type(2)*n-ordinal_type(1);
254 }
255 
256 template <typename ordinal_type, typename value_type>
257 Teuchos::RCP<Stokhos::OneDOrthogPolyBasis<ordinal_type,value_type> >
258 Stokhos::PecosOneDOrthogPolyBasis<ordinal_type,value_type>::
259 cloneWithOrder(ordinal_type p) const
260 {
261  return
262  Teuchos::rcp(new Stokhos::PecosOneDOrthogPolyBasis<ordinal_type,value_type>(p,*this));
263 }
264 
265 template <typename ordinal_type, typename value_type>
266 ordinal_type
267 Stokhos::PecosOneDOrthogPolyBasis<ordinal_type,value_type>::
268 coefficientGrowth(ordinal_type n) const
269 {
270  return n;
271 }
272 
273 template <typename ordinal_type, typename value_type>
274 ordinal_type
275 Stokhos::PecosOneDOrthogPolyBasis<ordinal_type,value_type>::
276 pointGrowth(ordinal_type n) const
277 {
278  if (n % ordinal_type(2) == ordinal_type(1))
279  return n + ordinal_type(1);
280  return n;
281 }
j
j
Definition: Sacado_Fad_Exp_MP_Vector.hpp:495
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:1160
Sacado::UQ::abs
KOKKOS_INLINE_FUNCTION PCE< Storage > abs(const PCE< Storage > &a)
Definition: Sacado_UQ_PCE_Imp.hpp:1150
Teuchos::Array::resize
void resize(size_type new_size, const value_type &x=value_type())
Teuchos::Array::size
size_type size() const
Generated on Fri Nov 22 2024 09:18:13 for Stokhos Package Browser (Single Doxygen Collection) by   doxygen 1.8.5