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