Intrepid
test_24.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
51//#include "Intrepid_CubatureLineSorted.hpp"
53#include "Intrepid_Utils.hpp"
54#include "Teuchos_oblackholestream.hpp"
55#include "Teuchos_RCP.hpp"
56#include "Teuchos_GlobalMPISession.hpp"
57
58using namespace Intrepid;
59
60/*
61 Computes integrals of monomials over a given reference cell.
62*/
63long double evalQuad(int order, int power, EIntrepidBurkardt rule) {
64
65 CubatureLineSorted<long double> lineCub(order,rule,false);
66 int size = lineCub.getNumPoints();
67 FieldContainer<long double> cubPoints(size);
68 FieldContainer<long double> cubWeights(size);
69 lineCub.getCubature(cubPoints,cubWeights);
70
71 int mid = size/2;
72 long double Q = 0.0;
73 if (size%2)
74 Q = cubWeights(mid)*powl(cubPoints(mid),power);
75
76 for (int i=0; i<mid; i++) {
77 Q += cubWeights(i)*powl(cubPoints(i),power)+
78 cubWeights(size-i-1)*powl(cubPoints(size-i-1),power);
79 }
80 return Q;
81}
82
83long double evalInt(int power, EIntrepidBurkardt rule) {
84 double I = 0;
85
86 if (rule==BURK_CLENSHAWCURTIS||rule==BURK_FEJER2||
87 rule==BURK_LEGENDRE||rule==BURK_PATTERSON ||
88 rule==BURK_TRAPEZOIDAL) {
89 if (power%2)
90 I = 0.0;
91 else
92 I = 2.0/((long double)power+1.0);
93 }
94 else if (rule==BURK_LAGUERRE) {
95 I = tgammal((long double)(power+1));
96 }
97 else if (rule==BURK_CHEBYSHEV1) {
98 long double bot, top;
99 if (!(power%2)) {
100 top = 1; bot = 1;
101 for (int i=2;i<=power;i+=2) {
102 top *= (long double)(i-1);
103 bot *= (long double)i;
104 }
105 I = M_PI*top/bot;
106 }
107 else {
108 I = 0.0;
109 }
110 }
111 else if (rule==BURK_CHEBYSHEV2) {
112 long double bot, top;
113 if (!(power%2)) {
114 top = 1; bot = 1;
115 for (int i=2;i<=power;i+=2) {
116 top *= (long double)(i-1);
117 bot *= (long double)i;
118 }
119 bot *= (long double)(power+2);
120 I = M_PI*top/bot;
121 }
122 else {
123 I = 0.0;
124 }
125 }
126 else if (rule==BURK_HERMITE||rule==BURK_GENZKEISTER) {
127 if (power%2) {
128 I = 0.0;
129 }
130 else {
131 long double value = 1.0;
132 if ((power-1)>=1) {
133 int n_copy = power-1;
134 while (1<n_copy) {
135 value *= (long double)n_copy;
136 n_copy -= 2;
137 }
138 }
139 I = value*sqrt(M_PI)/powl(2.0,(long double)power/2.0);
140 }
141 }
142 return I;
143}
144
145int main(int argc, char *argv[]) {
146
147 Teuchos::GlobalMPISession mpiSession(&argc, &argv);
148
149 // This little trick lets us print to std::cout only if
150 // a (dummy) command-line argument is provided.
151 int iprint = argc - 1;
152 Teuchos::RCP<std::ostream> outStream;
153 Teuchos::oblackholestream bhs; // outputs nothing
154 if (iprint > 0)
155 outStream = Teuchos::rcp(&std::cout, false);
156 else
157 outStream = Teuchos::rcp(&bhs, false);
158
159 // Save the format state of the original std::cout.
160 Teuchos::oblackholestream oldFormatState;
161 oldFormatState.copyfmt(std::cout);
162
163 *outStream \
164 << "===============================================================================\n" \
165 << "| |\n" \
166 << "| Unit Test (CubatureLineSorted) |\n" \
167 << "| |\n" \
168 << "| 1) Computing integrals of monomials in 1D |\n" \
169 << "| |\n" \
170 << "| Questions? Contact Drew Kouri (dpkouri@sandia.gov) or |\n" \
171 << "| Denis Ridzal (dridzal@sandia.gov). |\n" \
172 << "| |\n" \
173 << "| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
174 << "| Trilinos website: http://trilinos.sandia.gov |\n" \
175 << "| |\n" \
176 << "===============================================================================\n"\
177 << "| TEST 11: integrals of monomials in 1D |\n"\
178 << "===============================================================================\n";
179
180 // internal variables:
181 int errorFlag = 0;
182 long double reltol = 1.0e+05*INTREPID_TOL;
183 int maxDeg = 0;
184 long double analyticInt = 0;
185 long double testInt = 0;
186 int maxOrder = 30;
187
188 *outStream << "\nIntegrals of monomials on a reference line (edge):\n";
189 // compute and compare integrals
190 try {
191 for (EIntrepidBurkardt rule=BURK_CHEBYSHEV1; rule <= BURK_LAGUERRE; rule++) {
192 *outStream << "Testing " << EIntrepidBurkardtToString(rule) << "\n";
193 // compute integrals
194 if (rule==BURK_HERMITE)
195 maxOrder = 10;
196 else if (rule==BURK_TRAPEZOIDAL)
197 maxOrder = 2;
198 else
199 maxOrder = 30;
200
201 if (rule!=BURK_PATTERSON&&rule!=BURK_GENZKEISTER) {
202 for (int i=1; i <= maxOrder; i++) {
203 if (rule==BURK_CHEBYSHEV1 ||
204 rule==BURK_CHEBYSHEV2 ||
205 rule==BURK_LEGENDRE ||
206 rule==BURK_LAGUERRE ||
207 rule==BURK_HERMITE)
208 maxDeg = 2*i-1;
209 else if (rule==BURK_CLENSHAWCURTIS ||
210 rule==BURK_FEJER2 ||
211 rule==BURK_TRAPEZOIDAL)
212 maxDeg = i-1;
213
214 for (int j=0; j <= maxDeg; j++) {
215 analyticInt = evalInt(j,rule);
216 testInt = evalQuad(maxDeg,j,rule);
217
218 long double abstol = (analyticInt == 0.0 ? reltol : std::fabs(reltol*analyticInt) );
219 long double absdiff = std::fabs(analyticInt - testInt);
220 *outStream << "Cubature order " << std::setw(2) << std::left
221 << i << " integrating "
222 << "x^" << std::setw(2) << std::left << j << ":"
223 << " "
224 << std::scientific << std::setprecision(16) << testInt
225 << " " << analyticInt
226 << " " << std::setprecision(4) << absdiff << " "
227 << "<?" << " " << abstol
228 << "\n";
229 if (absdiff > abstol) {
230 errorFlag++;
231 *outStream << std::right << std::setw(104)
232 << "^^^^---FAILURE!\n";
233 }
234 } // end for j
235 *outStream << "\n";
236 } // end for i
237 }
238 else if (rule==BURK_PATTERSON) {
239 for (int i=0; i < 8; i++) {
240 int l = (int)std::pow(2.0,(double)i+1.0)-1;
241 if (i==0)
242 maxDeg = 1;
243 else
244 maxDeg = (int)(1.5*(double)l+0.5);
245 for (int j=0; j <= maxDeg; j++) {
246 analyticInt = evalInt(j,rule);
247 testInt = evalQuad(maxDeg,j,rule);
248
249 long double abstol = (analyticInt == 0.0 ? reltol : std::fabs(reltol*analyticInt) );
250 long double absdiff = std::fabs(analyticInt - testInt);
251 *outStream << "Cubature order " << std::setw(2) << std::left
252 << l << " integrating "
253 << "x^" << std::setw(2) << std::left << j << ":"
254 << " "
255 << std::scientific << std::setprecision(16) << testInt
256 << " " << analyticInt
257 << " " << std::setprecision(4) << absdiff << " "
258 << "<?" << " " << abstol
259 << "\n";
260 if (absdiff > abstol) {
261 errorFlag++;
262 *outStream << std::right << std::setw(104)
263 << "^^^^---FAILURE!\n";
264 }
265 } // end for j
266 *outStream << "\n";
267 } // end for i
268 }
269 else if (rule==BURK_GENZKEISTER) {
270 reltol *= 1.0e+02;
271 int o_ghk[4] = {1,3, 9,19};
272 int p_ghk[4] = {1,5,15,29};
273 for (int i=0; i < 4; i++) {
274 int l = o_ghk[i];
275 maxDeg = p_ghk[i];
276 for (int j=0; j <= maxDeg; j++) {
277 analyticInt = evalInt(j,rule);
278 testInt = evalQuad(maxDeg,j,rule);
279
280 long double abstol = (analyticInt == 0.0 ? reltol : std::fabs(reltol*analyticInt) );
281 long double absdiff = std::fabs(analyticInt - testInt);
282 *outStream << "Cubature order " << std::setw(2) << std::left
283 << l << " integrating "
284 << "x^" << std::setw(2) << std::left << j << ":"
285 << " "
286 << std::scientific << std::setprecision(16) << testInt
287 << " " << analyticInt
288 << " " << std::setprecision(4) << absdiff << " "
289 << "<?" << " " << abstol
290 << "\n";
291 if (absdiff > abstol) {
292 errorFlag++;
293 *outStream << std::right << std::setw(104)
294 << "^^^^---FAILURE!\n";
295 }
296 } // end for j
297 *outStream << "\n";
298 } // end for i
299 }
300 } // end for rule
301 }
302 catch (const std::logic_error & err) {
303 *outStream << err.what() << "\n";
304 errorFlag = -1;
305 };
306
307
308 if (errorFlag != 0)
309 std::cout << "End Result: TEST FAILED\n";
310 else
311 std::cout << "End Result: TEST PASSED\n";
312
313 // reset format state of std::cout
314 std::cout.copyfmt(oldFormatState);
315
316 return errorFlag;
317}
Header file for the Intrepid::CubatureLineSorted class.
Intrepid utilities.