Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
Stokhos_TensorProductQuadratureImp.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#include "Teuchos_TimeMonitor.hpp"
44
45template <typename ordinal_type, typename value_type>
47TensorProductQuadrature(const Teuchos::RCP<const ProductBasis<ordinal_type,value_type> >& product_basis)
48{
49#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
50 TEUCHOS_FUNC_TIME_MONITOR("Stokhos::TensorProductQuadrature -- Quad Grid Generation");
51#endif
52 ordinal_type d = product_basis->dimension();
53 ordinal_type sz = product_basis->size();
54
55 Teuchos::Array< Teuchos::RCP<const OneDOrthogPolyBasis<ordinal_type,value_type> > > coordinate_bases = product_basis->getCoordinateBases();
56
57 // Compute quad points, weights, values
58 Teuchos::Array< Teuchos::Array<value_type> > gp(d);
59 Teuchos::Array< Teuchos::Array<value_type> > gw(d);
60 Teuchos::Array< Teuchos::Array< Teuchos::Array<value_type> > > gv(d);
61 Teuchos::Array<ordinal_type> n(d);
62 ordinal_type ntot = 1;
63 for (ordinal_type i=0; i<d; i++) {
64 coordinate_bases[i]->getQuadPoints(2*(coordinate_bases[i]->order()),
65 gp[i], gw[i], gv[i]);
66 n[i] = gp[i].size();
67 ntot *= n[i];
68 }
69 quad_points.resize(ntot);
70 quad_weights.resize(ntot);
71 quad_values.resize(ntot);
72 Teuchos::Array<ordinal_type> index(d);
73 for (ordinal_type i=0; i<d; i++)
74 index[i] = 0;
75 ordinal_type cnt = 0;
76 while (cnt < ntot) {
77 quad_points[cnt].resize(d);
78 quad_weights[cnt] = value_type(1.0);
79 quad_values[cnt].resize(sz);
80 for (ordinal_type j=0; j<d; j++) {
81 quad_points[cnt][j] = gp[j][index[j]];
82 quad_weights[cnt] *= gw[j][index[j]];
83 }
84 for (ordinal_type k=0; k<sz; k++) {
85 quad_values[cnt][k] = value_type(1.0);
86 const MultiIndex<ordinal_type>& term = product_basis->term(k);
87 for (ordinal_type j=0; j<d; j++)
88 quad_values[cnt][k] *= gv[j][index[j]][term[j]];
89 }
90 ++index[0];
91 ordinal_type i = 0;
92 while (i < d-1 && index[i] == n[i]) {
93 index[i] = 0;
94 ++i;
95 ++index[i];
96 }
97 ++cnt;
98 }
99
100 //std::cout << "Number of quadrature points = " << ntot << std::endl;
101}
102
103template <typename ordinal_type, typename value_type>
105TensorProductQuadrature(const Teuchos::RCP<const ProductBasis<ordinal_type,value_type> >& product_basis,const ordinal_type& quad_order)
106{
107#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
108 TEUCHOS_FUNC_TIME_MONITOR("Stokhos::TensorProductQuadrature -- Quad Grid Generation");
109#endif
110 ordinal_type d = product_basis->dimension();
111 ordinal_type sz = product_basis->size();
112
113 Teuchos::Array< Teuchos::RCP<const OneDOrthogPolyBasis<ordinal_type,value_type> > > coordinate_bases = product_basis->getCoordinateBases();
114
115 // Compute quad points, weights, values
116 Teuchos::Array< Teuchos::Array<value_type> > gp(d);
117 Teuchos::Array< Teuchos::Array<value_type> > gw(d);
118 Teuchos::Array< Teuchos::Array< Teuchos::Array<value_type> > > gv(d);
119 Teuchos::Array<ordinal_type> n(d);
120 ordinal_type ntot = 1;
121 for (ordinal_type i=0; i<d; i++) {
122 coordinate_bases[i]->getQuadPoints(quad_order,
123 gp[i], gw[i], gv[i]);
124 n[i] = gp[i].size();
125 ntot *= n[i];
126 }
127 quad_points.resize(ntot);
128 quad_weights.resize(ntot);
129 quad_values.resize(ntot);
130 Teuchos::Array<ordinal_type> index(d);
131 for (ordinal_type i=0; i<d; i++)
132 index[i] = 0;
133 ordinal_type cnt = 0;
134 while (cnt < ntot) {
135 quad_points[cnt].resize(d);
136 quad_weights[cnt] = value_type(1.0);
137 quad_values[cnt].resize(sz);
138 for (ordinal_type j=0; j<d; j++) {
139 quad_points[cnt][j] = gp[j][index[j]];
140 quad_weights[cnt] *= gw[j][index[j]];
141 }
142 for (ordinal_type k=0; k<sz; k++) {
143 quad_values[cnt][k] = value_type(1.0);
144 MultiIndex<ordinal_type> term = product_basis->term(k);
145 for (ordinal_type j=0; j<d; j++)
146 quad_values[cnt][k] *= gv[j][index[j]][term[j]];
147 }
148 ++index[0];
149 ordinal_type i = 0;
150 while (i < d-1 && index[i] == n[i]) {
151 index[i] = 0;
152 ++i;
153 ++index[i];
154 }
155 ++cnt;
156 }
157
158 //std::cout << "Number of quadrature points = " << ntot << std::endl;
159}
160
161template <typename ordinal_type, typename value_type>
162const Teuchos::Array< Teuchos::Array<value_type> >&
164getQuadPoints() const
165{
166 return quad_points;
167}
168
169template <typename ordinal_type, typename value_type>
170const Teuchos::Array<value_type>&
172getQuadWeights() const
173{
174 return quad_weights;
175}
176
177template <typename ordinal_type, typename value_type>
178const Teuchos::Array< Teuchos::Array<value_type> >&
181{
182 return quad_values;
183}
184
185template <typename ordinal_type, typename value_type>
186std::ostream&
188print(std::ostream& os) const
189{
190 ordinal_type nqp = quad_weights.size();
191 os << "Tensor Product Quadrature with " << nqp << " points:"
192 << std::endl << "Weight : Points" << std::endl;
193 for (ordinal_type i=0; i<nqp; i++) {
194 os << i << ": " << quad_weights[i] << " : ";
195 for (ordinal_type j=0; j<static_cast<ordinal_type>(quad_points[i].size());
196 j++)
197 os << quad_points[i][j] << " ";
198 os << std::endl;
199 }
200 os << "Basis values at quadrature points:" << std::endl;
201 for (ordinal_type i=0; i<nqp; i++) {
202 os << i << " " << ": ";
203 for (ordinal_type j=0; j<static_cast<ordinal_type>(quad_values[i].size());
204 j++)
205 os << quad_values[i][j] << " ";
206 os << std::endl;
207 }
208
209 return os;
210}
A multidimensional index.
Abstract base class for multivariate orthogonal polynomials generated from tensor products of univari...
virtual const Teuchos::Array< Teuchos::Array< value_type > > & getQuadPoints() const
Get quadrature points.
virtual const Teuchos::Array< value_type > & getQuadWeights() const
Get quadrature weights.
virtual std::ostream & print(std::ostream &os) const
Print quadrature data.
TensorProductQuadrature(const Teuchos::RCP< const ProductBasis< ordinal_type, value_type > > &product_basis)
Constructor.
virtual const Teuchos::Array< Teuchos::Array< value_type > > & getBasisAtQuadPoints() const
Get values of basis at quadrature points.