Intrepid
test_10.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_Utils.hpp"
53 #include "Teuchos_oblackholestream.hpp"
54 #include "Teuchos_RCP.hpp"
55 #include "Teuchos_GlobalMPISession.hpp"
56 
57 #define INTREPID_CUBATURE_LINE_MAX 61
58 
59 using namespace Intrepid;
60 
61 
62 /*
63  Monomial evaluation.
64  in 1D, for point p(x) : x^xDeg
65  in 2D, for point p(x,y) : x^xDeg * y^yDeg
66  in 3D, for point p(x,y,z): x^xDeg * y^yDeg * z^zDeg
67 */
68 double computeMonomial(FieldContainer<double> & p, int xDeg, int yDeg=0, int zDeg=0) {
69  double val = 1.0;
70  int polydeg[3];
71  polydeg[0] = xDeg; polydeg[1] = yDeg; polydeg[2] = zDeg;
72  for (int i=0; i<p.dimension(0); i++) {
73  val *= std::pow(p(i),polydeg[i]);
74  }
75  return val;
76 }
77 
78 
79 /*
80  Computes integrals of monomials over a given reference cell.
81 */
82 double computeIntegral(int cubDegree, int polyDegree, EIntrepidPLPoly poly_type) {
83 
84  CubaturePolylib<double> lineCub(cubDegree, poly_type);
85  double val = 0.0;
86 
87  int cubDim = lineCub.getDimension();
88 
89  int numCubPoints = lineCub.getNumPoints();
90 
91  FieldContainer<double> point(cubDim);
92  FieldContainer<double> cubPoints(numCubPoints, cubDim);
93  FieldContainer<double> cubWeights(numCubPoints);
94 
95  lineCub.getCubature(cubPoints, cubWeights);
96 
97  for (int i=0; i<numCubPoints; i++) {
98  for (int j=0; j<cubDim; j++) {
99  point(j) = cubPoints(i,j);
100  }
101  val += computeMonomial(point, polyDegree)*cubWeights(i);
102  }
103 
104  return val;
105 }
106 
107 
108 int main(int argc, char *argv[]) {
109 
110  Teuchos::GlobalMPISession mpiSession(&argc, &argv);
111 
112  // This little trick lets us print to std::cout only if
113  // a (dummy) command-line argument is provided.
114  int iprint = argc - 1;
115  Teuchos::RCP<std::ostream> outStream;
116  Teuchos::oblackholestream bhs; // outputs nothing
117  if (iprint > 0)
118  outStream = Teuchos::rcp(&std::cout, false);
119  else
120  outStream = Teuchos::rcp(&bhs, false);
121 
122  // Save the format state of the original std::cout.
123  Teuchos::oblackholestream oldFormatState;
124  oldFormatState.copyfmt(std::cout);
125 
126  *outStream \
127  << "===============================================================================\n" \
128  << "| |\n" \
129  << "| Unit Test (CubaturePolylib) |\n" \
130  << "| |\n" \
131  << "| 1) Computing integrals of monomials on reference cells in 1D |\n" \
132  << "| |\n" \
133  << "| Questions? Contact Pavel Bochev (pbboche@sandia.gov) or |\n" \
134  << "| Denis Ridzal (dridzal@sandia.gov). |\n" \
135  << "| |\n" \
136  << "| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
137  << "| Trilinos website: http://trilinos.sandia.gov |\n" \
138  << "| |\n" \
139  << "===============================================================================\n"\
140  << "| TEST 1: integrals of monomials in 1D |\n"\
141  << "===============================================================================\n";
142 
143  // internal variables:
144  int errorFlag = 0;
145  Teuchos::Array< Teuchos::Array<double> > testInt;
146  Teuchos::Array< Teuchos::Array<double> > analyticInt;
147  Teuchos::Array<double> tmparray(1);
148  double reltol = 1.0e+03 * INTREPID_TOL;
149  testInt.assign(INTREPID_CUBATURE_LINE_MAX+1, tmparray);
150  analyticInt.assign(INTREPID_CUBATURE_LINE_MAX+1, tmparray);
151 
152  // open file with analytic values
153  std::string basedir = "./data";
154  std::stringstream namestream;
155  std::string filename;
156  namestream << basedir << "/EDGE_integrals" << ".dat";
157  namestream >> filename;
158  std::ifstream filecompare(&filename[0]);
159 
160  *outStream << "\nIntegrals of monomials on a reference line (edge):\n";
161 
162  // compute and compare integrals
163  try {
164  for (EIntrepidPLPoly poly_type=PL_GAUSS; poly_type <= PL_GAUSS_LOBATTO; poly_type++) {
165  // compute integrals
166  for (int cubDeg=0; cubDeg <= INTREPID_CUBATURE_LINE_MAX; cubDeg++) {
167  testInt[cubDeg].resize(cubDeg+1);
168  for (int polyDeg=0; polyDeg <= cubDeg; polyDeg++) {
169  testInt[cubDeg][polyDeg] = computeIntegral(cubDeg, polyDeg, poly_type);
170  }
171  }
172  // get analytic values
173  if (filecompare.is_open()) {
174  getAnalytic(analyticInt, filecompare);
175  // close file
176  filecompare.close();
177  }
178  // perform comparison
179  for (int cubDeg=0; cubDeg <= INTREPID_CUBATURE_LINE_MAX; cubDeg++) {
180  for (int polyDeg=0; polyDeg <= cubDeg; polyDeg++) {
181  double abstol = ( analyticInt[polyDeg][0] == 0.0 ? reltol : std::fabs(reltol*analyticInt[polyDeg][0]) );
182  double absdiff = std::fabs(analyticInt[polyDeg][0] - testInt[cubDeg][polyDeg]);
183  *outStream << "Cubature order " << std::setw(2) << std::left << cubDeg << " integrating "
184  << "x^" << std::setw(2) << std::left << polyDeg << ":" << " "
185  << std::scientific << std::setprecision(16) << testInt[cubDeg][polyDeg] << " " << analyticInt[polyDeg][0] << " "
186  << std::setprecision(4) << absdiff << " " << "<?" << " " << abstol << "\n";
187  if (absdiff > abstol) {
188  errorFlag++;
189  *outStream << std::right << std::setw(104) << "^^^^---FAILURE!\n";
190  }
191  }
192  *outStream << "\n";
193  } // end for cubDeg
194  } // end for poly_type
195  }
196  catch (std::logic_error err) {
197  *outStream << err.what() << "\n";
198  errorFlag = -1;
199  };
200 
201 
202  if (errorFlag != 0)
203  std::cout << "End Result: TEST FAILED\n";
204  else
205  std::cout << "End Result: TEST PASSED\n";
206 
207  // reset format state of std::cout
208  std::cout.copyfmt(oldFormatState);
209 
210  return errorFlag;
211 }
int dimension(const int whichDim) const
Returns the specified dimension.
Intrepid utilities.
Utilizes cubature (integration) rules contained in the library Polylib (Spencer Sherwin, Aeronautics, Imperial College London) within Intrepid.
Header file for the Intrepid::CubaturePolylib class.