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);
69 AdaptiveSparseGrid<long double,std::vector<long double> >::buildSparseGrid(
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.