Intrepid
test_22.cpp
Go to the documentation of this file.
1// @HEADER
2// ************************************************************************
3//
4// Intrepid Package
5// Copyright (2007) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Pavel Bochev (pbboche@sandia.gov)
38// Denis Ridzal (dridzal@sandia.gov), or
39// Kara Peterson (kjpeter@sandia.gov)
40//
41// ************************************************************************
42// @HEADER
43
44
53#include "Intrepid_Utils.hpp"
54#include "Teuchos_oblackholestream.hpp"
55#include "Teuchos_RCP.hpp"
56#include "Teuchos_RefCountPtr.hpp"
57#include "Teuchos_GlobalMPISession.hpp"
58
59using namespace Intrepid;
60
61/*
62 Computes integrals of monomials over a given reference cell.
63*/
64long double evalQuad(std::vector<int> power, int dimension, int order,
65 std::vector<EIntrepidBurkardt> rule,
66 std::vector<EIntrepidGrowth> growth) {
67
68 CubatureTensorSorted<long double> lineCub(0,dimension);
70 lineCub,dimension,
71 order,rule,
72 growth,false);
73
74 int size = lineCub.getNumPoints();
75 FieldContainer<long double> cubPoints(size,dimension);
76 FieldContainer<long double> cubWeights(size);
77 lineCub.getCubature(cubPoints,cubWeights);
78
79 int mid = size/2;
80 long double Q = 0.0;
81 if (size%2) {
82 Q = cubWeights(mid);
83 for (int i=0; i<dimension; i++) {
84 Q *= powl(cubPoints(mid,i),(long double)power[i]);
85 }
86 }
87
88 for (int i=0; i<mid; i++) {
89 long double value1 = cubWeights(i);
90 long double value2 = cubWeights(size-i-1);
91 for (int j=0; j<dimension; j++) {
92 value1 *= powl(cubPoints(i,j),(long double)power[j]);
93 value2 *= powl(cubPoints(size-i-1,j),(long double)power[j]);
94 }
95 Q += value1+value2;
96 }
97 return Q;
98}
99
100long double evalInt(int dimension, std::vector<int> power,
101 std::vector<EIntrepidBurkardt> rule) {
102 long double I = 1.0;
103
104 for (int i=0; i<dimension; i++) {
105 if (rule[i]==BURK_CLENSHAWCURTIS||rule[i]==BURK_FEJER2||
106 rule[i]==BURK_LEGENDRE||rule[i]==BURK_PATTERSON ||
107 rule[i]==BURK_TRAPEZOIDAL) {
108 if (power[i]%2)
109 I *= 0.0;
110 else
111 I *= 2.0/((long double)power[i]+1.0);
112 }
113 else if (rule[i]==BURK_LAGUERRE) {
114 I *= tgammal((long double)(power[i]+1));
115 }
116 else if (rule[i]==BURK_CHEBYSHEV1) {
117 long double bot, top;
118 if (!(power[i]%2)) {
119 top = 1; bot = 1;
120 for (int j=2;j<=power[i];j+=2) {
121 top *= (long double)(j-1);
122 bot *= (long double)j;
123 }
124 I *= M_PI*top/bot;
125 }
126 else {
127 I *= 0.0;
128 }
129 }
130 else if (rule[i]==BURK_CHEBYSHEV2) {
131 long double bot, top;
132 if (!(power[i]%2)) {
133 top = 1; bot = 1;
134 for (int j=2;j<=power[i];j+=2) {
135 top *= (long double)(j-1);
136 bot *= (long double)j;
137 }
138 bot *= (long double)(power[i]+2);
139 I *= M_PI*top/bot;
140 }
141 else {
142 I *= 0.0;
143 }
144 }
145 else if (rule[i]==BURK_HERMITE||rule[i]==BURK_GENZKEISTER) {
146 if (power[i]%2) {
147 I *= 0.0;
148 }
149 else {
150 long double value = 1.0;
151 if ((power[i]-1)>=1) {
152 int n_copy = power[i]-1;
153 while (1<n_copy) {
154 value *= (long double)n_copy;
155 n_copy -= 2;
156 }
157 }
158 I *= value*sqrt(M_PI)/powl(2.0,(long double)power[i]/2.0);
159 }
160 }
161 }
162 return I;
163}
164
165int main(int argc, char *argv[]) {
166
167 Teuchos::GlobalMPISession mpiSession(&argc, &argv);
168
169 // This little trick lets us print to std::cout only if
170 // a (dummy) command-line argument is provided.
171 int iprint = argc - 1;
172 Teuchos::RCP<std::ostream> outStream;
173 Teuchos::oblackholestream bhs; // outputs nothing
174 if (iprint > 0)
175 outStream = Teuchos::rcp(&std::cout, false);
176 else
177 outStream = Teuchos::rcp(&bhs, false);
178
179 // Save the format state of the original std::cout.
180 Teuchos::oblackholestream oldFormatState;
181 oldFormatState.copyfmt(std::cout);
182
183 *outStream \
184 << "===============================================================================\n" \
185 << "| |\n" \
186 << "| Unit Test (CubatureTensorSorted) |\n" \
187 << "| |\n" \
188 << "| 1) Computing integrals of monomials in 2D |\n" \
189 << "| |\n" \
190 << "| Questions? Contact Drew Kouri (dpkouri@sandia.gov) or |\n" \
191 << "| Denis Ridzal (dridzal@sandia.gov). |\n" \
192 << "| |\n" \
193 << "| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
194 << "| Trilinos website: http://trilinos.sandia.gov |\n" \
195 << "| |\n" \
196 << "===============================================================================\n"\
197 << "| TEST 22: integrals of monomials in 2D - Anisotropic but no growth rules |\n"\
198 << "===============================================================================\n";
199
200
201 // internal variables:
202 int dimension = 2;
203 int errorFlag = 0;
204 long double reltol = 1.0e+02*INTREPID_TOL;
205 int maxDeg = 0;
206 long double analyticInt = 0;
207 long double testInt = 0;
208 int maxOrder = 6;
209 std::vector<int> power(2,0);
210 std::vector<EIntrepidBurkardt> rule1(2,BURK_CLENSHAWCURTIS);
211 std::vector<EIntrepidGrowth> growth1(2,GROWTH_FULLEXP);
212
213 *outStream << "\nIntegrals of monomials on a reference line (edge):\n";
214 // compute and compare integrals
215 try {
216 for (int i=0; i<=maxOrder; i++) {
217 maxDeg = i-1;
218 for (int j=0; j <= maxDeg; j++) {
219 power[0] = j;
220 for (int k=0; k <= maxDeg; k++) {
221 power[1] = k;
222 if (j+k < maxDeg) {
223 analyticInt = evalInt(dimension, power, rule1);
224 testInt = evalQuad(power,dimension,maxOrder,rule1,growth1);
225
226 long double abstol = (analyticInt == 0.0 ? reltol : std::fabs(reltol*analyticInt) );
227 long double absdiff = std::fabs(analyticInt - testInt);
228 *outStream << "Cubature order " << std::setw(2) << std::left << i
229 << " integrating " << "x^" << std::setw(2)
230 << std::left << j << "y^" << std::setw(2) << std::left
231 << k << ":" << " " << std::scientific
232 << std::setprecision(16) << testInt
233 << " " << analyticInt << " "
234 << std::setprecision(4) << absdiff << " "
235 << "<?" << " " << abstol << "\n";
236 if (absdiff > abstol) {
237 errorFlag++;
238 *outStream << std::right << std::setw(104) << "^^^^---FAILURE!\n";
239 }
240 }
241 } // end for k
242 *outStream << "\n";
243 } // end for j
244 *outStream << "\n";
245 } // end for i
246 }
247 catch (const std::logic_error & err) {
248 *outStream << err.what() << "\n";
249 errorFlag = -1;
250 };
251
252
253 if (errorFlag != 0)
254 std::cout << "End Result: TEST FAILED\n";
255 else
256 std::cout << "End Result: TEST PASSED\n";
257
258 // reset format state of std::cout
259 std::cout.copyfmt(oldFormatState);
260
261 return errorFlag;
262}
Header file for the Intrepid::AdaptiveSparseGrid class.
Header file for the Intrepid::CubatureTensorSorted class.
Intrepid utilities.
Builds general adaptive sparse grid rules (Gerstner and Griebel) using the 1D cubature rules in the I...
Utilizes 1D cubature (integration) rules contained in the library sandia_rules (John Burkardt,...
Implementation of a templated lexicographical container for a multi-indexed scalar quantity....