Intrepid
test_16.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,
64 std::vector<int> order,
65 std::vector<EIntrepidBurkardt> rule) {
66
67 CubatureTensorSorted<long double> lineCub(dimension,order,rule,false);
68 order[0]--; order[1]--;
69 CubatureTensorSorted<long double> lineCub1(dimension,order,rule,false);
70 lineCub.update(1.0,lineCub1,1.0);
71
72 int size = lineCub.getNumPoints();
73 FieldContainer<long double> cubPoints(size,dimension);
74 FieldContainer<long double> cubWeights(size);
75 lineCub.getCubature(cubPoints,cubWeights);
76
77 int mid = size/2;
78 long double Q = 0.0;
79 if (size%2) {
80 Q = cubWeights(mid);
81 for (int i=0; i<dimension; i++) {
82 Q *= powl(cubPoints(mid,i),(long double)power[i]);
83 }
84 }
85
86 for (int i=0; i<mid; i++) {
87 long double value1 = cubWeights(i);
88 long double value2 = cubWeights(size-i-1);
89 for (int j=0; j<dimension; j++) {
90 value1 *= powl(cubPoints(i,j),(long double)power[j]);
91 value2 *= powl(cubPoints(size-i-1,j),(long double)power[j]);
92 }
93 Q += value1+value2;
94 }
95 return Q;
96}
97
98long double evalInt(int dimension, std::vector<int> power,
99 std::vector<EIntrepidBurkardt> rule) {
100 long double I = 1.0;
101
102 for (int i=0; i<dimension; i++) {
103 if (rule[i]==BURK_CLENSHAWCURTIS||rule[i]==BURK_FEJER2||
104 rule[i]==BURK_LEGENDRE||rule[i]==BURK_PATTERSON ||
105 rule[i]==BURK_TRAPEZOIDAL) {
106 if (power[i]%2)
107 I *= 0.0;
108 else
109 I *= 2.0/((long double)power[i]+1.0);
110 }
111 else if (rule[i]==BURK_LAGUERRE) {
112 I *= tgammal((long double)(power[i]+1));
113 }
114 else if (rule[i]==BURK_CHEBYSHEV1) {
115 long double bot, top;
116 if (!(power[i]%2)) {
117 top = 1; bot = 1;
118 for (int j=2;j<=power[i];j+=2) {
119 top *= (long double)(j-1);
120 bot *= (long double)j;
121 }
122 I *= M_PI*top/bot;
123 }
124 else {
125 I *= 0.0;
126 }
127 }
128 else if (rule[i]==BURK_CHEBYSHEV2) {
129 long double bot, top;
130 if (!(power[i]%2)) {
131 top = 1; bot = 1;
132 for (int j=2;j<=power[i];j+=2) {
133 top *= (long double)(j-1);
134 bot *= (long double)j;
135 }
136 bot *= (long double)(power[i]+2);
137 I *= M_PI*top/bot;
138 }
139 else {
140 I *= 0.0;
141 }
142 }
143 else if (rule[i]==BURK_HERMITE||rule[i]==BURK_GENZKEISTER) {
144 if (power[i]%2) {
145 I *= 0.0;
146 }
147 else {
148 long double value = 1.0;
149 if ((power[i]-1)>=1) {
150 int n_copy = power[i]-1;
151 while (1<n_copy) {
152 value *= (long double)n_copy;
153 n_copy -= 2;
154 }
155 }
156 I *= value*sqrt(M_PI)/powl(2.0,(long double)power[i]/2.0);
157 }
158 }
159 }
160 return I;
161}
162
163int main(int argc, char *argv[]) {
164
165 Teuchos::GlobalMPISession mpiSession(&argc, &argv);
166
167 // This little trick lets us print to std::cout only if
168 // a (dummy) command-line argument is provided.
169 int iprint = argc - 1;
170 Teuchos::RCP<std::ostream> outStream;
171 Teuchos::oblackholestream bhs; // outputs nothing
172 if (iprint > 0)
173 outStream = Teuchos::rcp(&std::cout, false);
174 else
175 outStream = Teuchos::rcp(&bhs, false);
176
177 // Save the format state of the original std::cout.
178 Teuchos::oblackholestream oldFormatState;
179 oldFormatState.copyfmt(std::cout);
180
181 *outStream \
182 << "===============================================================================\n" \
183 << "| |\n" \
184 << "| Unit Test (CubatureTensorSorted) |\n" \
185 << "| |\n" \
186 << "| 1) Computing integrals of monomials in 2D |\n" \
187 << "| |\n" \
188 << "| Questions? Contact Drew Kouri (dpkouri@sandia.gov) or |\n" \
189 << "| Denis Ridzal (dridzal@sandia.gov). |\n" \
190 << "| |\n" \
191 << "| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
192 << "| Trilinos website: http://trilinos.sandia.gov |\n" \
193 << "| |\n" \
194 << "===============================================================================\n"\
195 << "| TEST 12: integrals of monomials in 2D - Anisotropic but no growth rules |\n"\
196 << "===============================================================================\n";
197
198
199 // internal variables:
200 int dimension = 2;
201 int errorFlag = 0;
202 long double reltol = 1.0e+02*INTREPID_TOL;
203 int maxDeg = 0;
204 long double analyticInt = 0;
205 long double testInt = 0;
206 int maxOrder = 20;
207 std::vector<int> power(2,0);
208 std::vector<EIntrepidBurkardt> rule1(2,BURK_CLENSHAWCURTIS);
209 std::vector<int> order(2,0);
210
211 *outStream << "\nIntegrals of monomials on a reference line (edge):\n";
212 // compute and compare integrals
213 try {
214 for (int i=0; i<=maxOrder; i++) {
215 maxDeg = i-2;
216
217 order[0] = i; order[1] = i;
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 analyticInt = 2.0*evalInt(dimension, power, rule1);
223 testInt = evalQuad(power,dimension,order,rule1);
224
225 long double abstol = (analyticInt == 0.0 ? reltol : std::fabs(reltol*analyticInt) );
226 long double absdiff = std::fabs(analyticInt - testInt);
227 *outStream << "Cubature order " << std::setw(2) << std::left << i
228 << " integrating "
229 << "x^" << std::setw(2) << std::left << j << "y^"
230 << std::setw(2) << std::left
231 << k << ":" << " " << std::scientific
232 << std::setprecision(16) << testInt
233 << " " << analyticInt << " " << std::setprecision(4)
234 << absdiff << " "
235 << "<?" << " " << abstol << "\n";
236 if (absdiff > abstol) {
237 errorFlag++;
238 *outStream << std::right << std::setw(104) << "^^^^---FAILURE!\n";
239 }
240 } // end for k
241 *outStream << "\n";
242 } // end for j
243 *outStream << "\n";
244 } // end for i
245 }
246 catch (const std::logic_error & err) {
247 *outStream << err.what() << "\n";
248 errorFlag = -1;
249 };
250
251
252 if (errorFlag != 0)
253 std::cout << "End Result: TEST FAILED\n";
254 else
255 std::cout << "End Result: TEST PASSED\n";
256
257 // reset format state of std::cout
258 std::cout.copyfmt(oldFormatState);
259
260 return errorFlag;
261}
Header file for the Intrepid::CubatureTensorSorted class.
Intrepid utilities.