Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
Stokhos_GaussPattersonLegendreBasisImp.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#ifdef HAVE_STOKHOS_DAKOTA
45#include "sandia_rules.hpp"
46#endif
47#include "Teuchos_TestForException.hpp"
48
49template <typename ordinal_type, typename value_type>
51GaussPattersonLegendreBasis(ordinal_type p, bool normalize, bool isotropic_) :
52 LegendreBasis<ordinal_type, value_type>(p, normalize),
53 isotropic(isotropic_)
54{
55#ifdef HAVE_STOKHOS_DAKOTA
56 this->setSparseGridGrowthRule(webbur::level_to_order_exp_gp);
57#endif
58}
59
60template <typename ordinal_type, typename value_type>
62GaussPattersonLegendreBasis(ordinal_type p,
63 const GaussPattersonLegendreBasis& basis) :
64 LegendreBasis<ordinal_type, value_type>(p, basis),
65 isotropic(basis.isotropic)
66{
67}
68
69template <typename ordinal_type, typename value_type>
74
75template <typename ordinal_type, typename value_type>
76void
78getQuadPoints(ordinal_type quad_order,
79 Teuchos::Array<value_type>& quad_points,
80 Teuchos::Array<value_type>& quad_weights,
81 Teuchos::Array< Teuchos::Array<value_type> >& quad_values) const
82{
83#ifdef HAVE_STOKHOS_DAKOTA
84 // Gauss-Patterson points have the following structure
85 // (cf. http://people.sc.fsu.edu/~jburkardt/f_src/patterson_rule/patterson_rule.html):
86 // Level l Num points n Precision p
87 // -----------------------------------
88 // 0 1 1
89 // 1 3 5
90 // 2 7 11
91 // 3 15 23
92 // 4 31 47
93 // 5 63 95
94 // 6 127 191
95 // 7 255 383
96 // Thus for l > 0, n = 2^{l+1}-1 and p = 3*2^l-1. So for a given quadrature
97 // order p, we find the smallest l s.t. 3*s^l-1 >= p and then compute the
98 // number of points n from the above. In this case, l = ceil(log2((p+1)/3))
99 ordinal_type num_points;
100 if (quad_order <= ordinal_type(1))
101 num_points = 1;
102 else {
103 ordinal_type l = std::ceil(std::log((quad_order+1.0)/3.0)/std::log(2.0));
104 num_points = (1 << (l+1)) - 1; // std::pow(2,l+1)-1;
105 }
106
107 quad_points.resize(num_points);
108 quad_weights.resize(num_points);
109 quad_values.resize(num_points);
110
111 webbur::patterson_lookup(num_points, &quad_points[0], &quad_weights[0]);
112
113 for (ordinal_type i=0; i<num_points; i++) {
114 quad_weights[i] *= 0.5; // scale to unit measure
115 quad_values[i].resize(this->p+1);
116 this->evaluateBases(quad_points[i], quad_values[i]);
117 }
118
119#else
120 TEUCHOS_TEST_FOR_EXCEPTION(
121 true, std::logic_error, "Clenshaw-Curtis requires TriKota to be enabled!");
122#endif
123}
124
125template <typename ordinal_type, typename value_type>
126ordinal_type
128quadDegreeOfExactness(ordinal_type n) const
129{
130 // Based on the above structure, we find the largest l s.t. 2^{l+1}-1 <= n,
131 // which is floor(log2(n+1)-1) and compute p = 3*2^l-1
132 if (n == ordinal_type(1))
133 return 1;
134 ordinal_type l = std::floor(std::log(n+1.0)/std::log(2.0)-1.0);
135 return (3 << l) - 1; // 3*std::pow(2,l)-1;
136}
137
138template <typename ordinal_type, typename value_type>
139Teuchos::RCP<Stokhos::OneDOrthogPolyBasis<ordinal_type,value_type> >
146
147template <typename ordinal_type, typename value_type>
148ordinal_type
150coefficientGrowth(ordinal_type n) const
151{
152 // Gauss-Patterson rules have precision 3*2^l-1, which is odd.
153 // Since discrete orthogonality requires integrating polynomials of
154 // order 2*p, setting p = 3*2^{l-1}-1 will yield the largest p such that
155 // 2*p <= 3*2^l-1
156 if (n == 0)
157 return 0;
158 return (3 << (n-1)) - 1; // 3*std::pow(2,n-1) - 1;
159}
160
161template <typename ordinal_type, typename value_type>
162ordinal_type
Legendre polynomial basis using Gauss-Patterson quadrature points.
virtual void getQuadPoints(ordinal_type quad_order, Teuchos::Array< value_type > &points, Teuchos::Array< value_type > &weights, Teuchos::Array< Teuchos::Array< value_type > > &values) const
Compute quadrature points, weights, and values of basis polynomials at given set of points points.
GaussPattersonLegendreBasis(ordinal_type p, bool normalize=false, bool isotropic=false)
Constructor.
virtual ordinal_type pointGrowth(ordinal_type n) const
Evaluate point growth rule for Smolyak-type bases.
virtual Teuchos::RCP< OneDOrthogPolyBasis< ordinal_type, value_type > > cloneWithOrder(ordinal_type p) const
Clone this object with the option of building a higher order basis.
virtual ordinal_type coefficientGrowth(ordinal_type n) const
Evaluate coefficient growth rule for Smolyak-type bases.
virtual ordinal_type quadDegreeOfExactness(ordinal_type n) const
Legendre polynomial basis.
virtual void setSparseGridGrowthRule(LevelToOrderFnPtr ptr)
Set sparse grid rule.