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