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