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 
58 using namespace Intrepid;
59 
60 /*
61  Computes integrals of monomials over a given reference cell.
62 */
63 long 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 
98 long 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 
163 int 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 (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 }
Intrepid utilities.
Implementation of a templated lexicographical container for a multi-indexed scalar quantity...
Header file for the Intrepid::CubatureTensorSorted class.
Utilizes 1D cubature (integration) rules contained in the library sandia_rules (John Burkardt...