|
Intrepid
|
00001 // @HEADER 00002 // ************************************************************************ 00003 // 00004 // Intrepid Package 00005 // Copyright (2007) Sandia Corporation 00006 // 00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive 00008 // license for use of this work by or on behalf of the U.S. Government. 00009 // 00010 // Redistribution and use in source and binary forms, with or without 00011 // modification, are permitted provided that the following conditions are 00012 // met: 00013 // 00014 // 1. Redistributions of source code must retain the above copyright 00015 // notice, this list of conditions and the following disclaimer. 00016 // 00017 // 2. Redistributions in binary form must reproduce the above copyright 00018 // notice, this list of conditions and the following disclaimer in the 00019 // documentation and/or other materials provided with the distribution. 00020 // 00021 // 3. Neither the name of the Corporation nor the names of the 00022 // contributors may be used to endorse or promote products derived from 00023 // this software without specific prior written permission. 00024 // 00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY 00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR 00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE 00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, 00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, 00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR 00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF 00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING 00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS 00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 00036 // 00037 // Questions? Contact Pavel Bochev (pbboche@sandia.gov) 00038 // Denis Ridzal (dridzal@sandia.gov), or 00039 // Kara Peterson (kjpeter@sandia.gov) 00040 // 00041 // ************************************************************************ 00042 // @HEADER 00043 00044 00051 //#include "Intrepid_CubatureLineSorted.hpp" 00052 #include "Intrepid_CubatureLineSorted.hpp" 00053 #include "Intrepid_Utils.hpp" 00054 #include "Teuchos_oblackholestream.hpp" 00055 #include "Teuchos_RCP.hpp" 00056 #include "Teuchos_GlobalMPISession.hpp" 00057 00058 using namespace Intrepid; 00059 00060 /* 00061 Computes integrals of monomials over a given reference cell. 00062 */ 00063 long double evalQuad(int order, int power, EIntrepidBurkardt rule) { 00064 00065 CubatureLineSorted<long double> lineCub(order,rule,false); 00066 int size = lineCub.getNumPoints(); 00067 FieldContainer<long double> cubPoints(size); 00068 FieldContainer<long double> cubWeights(size); 00069 lineCub.getCubature(cubPoints,cubWeights); 00070 00071 int mid = size/2; 00072 long double Q = 0.0; 00073 if (size%2) 00074 Q = cubWeights(mid)*powl(cubPoints(mid),power); 00075 00076 for (int i=0; i<mid; i++) { 00077 Q += cubWeights(i)*powl(cubPoints(i),power)+ 00078 cubWeights(size-i-1)*powl(cubPoints(size-i-1),power); 00079 } 00080 return Q; 00081 } 00082 00083 long double evalInt(int power, EIntrepidBurkardt rule) { 00084 double I = 0; 00085 00086 if (rule==BURK_CLENSHAWCURTIS||rule==BURK_FEJER2|| 00087 rule==BURK_LEGENDRE||rule==BURK_PATTERSON || 00088 rule==BURK_TRAPEZOIDAL) { 00089 if (power%2) 00090 I = 0.0; 00091 else 00092 I = 2.0/((long double)power+1.0); 00093 } 00094 else if (rule==BURK_LAGUERRE) { 00095 I = tgammal((long double)(power+1)); 00096 } 00097 else if (rule==BURK_CHEBYSHEV1) { 00098 long double bot, top; 00099 if (!(power%2)) { 00100 top = 1; bot = 1; 00101 for (int i=2;i<=power;i+=2) { 00102 top *= (long double)(i-1); 00103 bot *= (long double)i; 00104 } 00105 I = M_PI*top/bot; 00106 } 00107 else { 00108 I = 0.0; 00109 } 00110 } 00111 else if (rule==BURK_CHEBYSHEV2) { 00112 long double bot, top; 00113 if (!(power%2)) { 00114 top = 1; bot = 1; 00115 for (int i=2;i<=power;i+=2) { 00116 top *= (long double)(i-1); 00117 bot *= (long double)i; 00118 } 00119 bot *= (long double)(power+2); 00120 I = M_PI*top/bot; 00121 } 00122 else { 00123 I = 0.0; 00124 } 00125 } 00126 else if (rule==BURK_HERMITE||rule==BURK_GENZKEISTER) { 00127 if (power%2) { 00128 I = 0.0; 00129 } 00130 else { 00131 long double value = 1.0; 00132 if ((power-1)>=1) { 00133 int n_copy = power-1; 00134 while (1<n_copy) { 00135 value *= (long double)n_copy; 00136 n_copy -= 2; 00137 } 00138 } 00139 I = value*sqrt(M_PI)/powl(2.0,(long double)power/2.0); 00140 } 00141 } 00142 return I; 00143 } 00144 00145 int main(int argc, char *argv[]) { 00146 00147 Teuchos::GlobalMPISession mpiSession(&argc, &argv); 00148 00149 // This little trick lets us print to std::cout only if 00150 // a (dummy) command-line argument is provided. 00151 int iprint = argc - 1; 00152 Teuchos::RCP<std::ostream> outStream; 00153 Teuchos::oblackholestream bhs; // outputs nothing 00154 if (iprint > 0) 00155 outStream = Teuchos::rcp(&std::cout, false); 00156 else 00157 outStream = Teuchos::rcp(&bhs, false); 00158 00159 // Save the format state of the original std::cout. 00160 Teuchos::oblackholestream oldFormatState; 00161 oldFormatState.copyfmt(std::cout); 00162 00163 *outStream \ 00164 << "===============================================================================\n" \ 00165 << "| |\n" \ 00166 << "| Unit Test (CubatureLineSorted) |\n" \ 00167 << "| |\n" \ 00168 << "| 1) Computing integrals of monomials in 1D |\n" \ 00169 << "| |\n" \ 00170 << "| Questions? Contact Drew Kouri (dpkouri@sandia.gov) or |\n" \ 00171 << "| Denis Ridzal (dridzal@sandia.gov). |\n" \ 00172 << "| |\n" \ 00173 << "| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \ 00174 << "| Trilinos website: http://trilinos.sandia.gov |\n" \ 00175 << "| |\n" \ 00176 << "===============================================================================\n"\ 00177 << "| TEST 11: integrals of monomials in 1D |\n"\ 00178 << "===============================================================================\n"; 00179 00180 // internal variables: 00181 int errorFlag = 0; 00182 long double reltol = 1.0e+05*INTREPID_TOL; 00183 int maxDeg = 0; 00184 long double analyticInt = 0; 00185 long double testInt = 0; 00186 int maxOrder = 30; 00187 00188 *outStream << "\nIntegrals of monomials on a reference line (edge):\n"; 00189 // compute and compare integrals 00190 try { 00191 for (EIntrepidBurkardt rule=BURK_CHEBYSHEV1; rule <= BURK_LAGUERRE; rule++) { 00192 *outStream << "Testing " << EIntrepidBurkardtToString(rule) << "\n"; 00193 // compute integrals 00194 if (rule==BURK_HERMITE) 00195 maxOrder = 10; 00196 else if (rule==BURK_TRAPEZOIDAL) 00197 maxOrder = 2; 00198 else 00199 maxOrder = 30; 00200 00201 if (rule!=BURK_PATTERSON&&rule!=BURK_GENZKEISTER) { 00202 for (int i=1; i <= maxOrder; i++) { 00203 if (rule==BURK_CHEBYSHEV1 || 00204 rule==BURK_CHEBYSHEV2 || 00205 rule==BURK_LEGENDRE || 00206 rule==BURK_LAGUERRE || 00207 rule==BURK_HERMITE) 00208 maxDeg = 2*i-1; 00209 else if (rule==BURK_CLENSHAWCURTIS || 00210 rule==BURK_FEJER2 || 00211 rule==BURK_TRAPEZOIDAL) 00212 maxDeg = i-1; 00213 00214 for (int j=0; j <= maxDeg; j++) { 00215 analyticInt = evalInt(j,rule); 00216 testInt = evalQuad(maxDeg,j,rule); 00217 00218 long double abstol = (analyticInt == 0.0 ? reltol : std::fabs(reltol*analyticInt) ); 00219 long double absdiff = std::fabs(analyticInt - testInt); 00220 *outStream << "Cubature order " << std::setw(2) << std::left 00221 << i << " integrating " 00222 << "x^" << std::setw(2) << std::left << j << ":" 00223 << " " 00224 << std::scientific << std::setprecision(16) << testInt 00225 << " " << analyticInt 00226 << " " << std::setprecision(4) << absdiff << " " 00227 << "<?" << " " << abstol 00228 << "\n"; 00229 if (absdiff > abstol) { 00230 errorFlag++; 00231 *outStream << std::right << std::setw(104) 00232 << "^^^^---FAILURE!\n"; 00233 } 00234 } // end for j 00235 *outStream << "\n"; 00236 } // end for i 00237 } 00238 else if (rule==BURK_PATTERSON) { 00239 for (int i=0; i < 8; i++) { 00240 int l = (int)std::pow(2.0,(double)i+1.0)-1; 00241 if (i==0) 00242 maxDeg = 1; 00243 else 00244 maxDeg = (int)(1.5*(double)l+0.5); 00245 for (int j=0; j <= maxDeg; j++) { 00246 analyticInt = evalInt(j,rule); 00247 testInt = evalQuad(maxDeg,j,rule); 00248 00249 long double abstol = (analyticInt == 0.0 ? reltol : std::fabs(reltol*analyticInt) ); 00250 long double absdiff = std::fabs(analyticInt - testInt); 00251 *outStream << "Cubature order " << std::setw(2) << std::left 00252 << l << " integrating " 00253 << "x^" << std::setw(2) << std::left << j << ":" 00254 << " " 00255 << std::scientific << std::setprecision(16) << testInt 00256 << " " << analyticInt 00257 << " " << std::setprecision(4) << absdiff << " " 00258 << "<?" << " " << abstol 00259 << "\n"; 00260 if (absdiff > abstol) { 00261 errorFlag++; 00262 *outStream << std::right << std::setw(104) 00263 << "^^^^---FAILURE!\n"; 00264 } 00265 } // end for j 00266 *outStream << "\n"; 00267 } // end for i 00268 } 00269 else if (rule==BURK_GENZKEISTER) { 00270 reltol *= 1.0e+02; 00271 int o_ghk[4] = {1,3, 9,19}; 00272 int p_ghk[4] = {1,5,15,29}; 00273 for (int i=0; i < 4; i++) { 00274 int l = o_ghk[i]; 00275 maxDeg = p_ghk[i]; 00276 for (int j=0; j <= maxDeg; j++) { 00277 analyticInt = evalInt(j,rule); 00278 testInt = evalQuad(maxDeg,j,rule); 00279 00280 long double abstol = (analyticInt == 0.0 ? reltol : std::fabs(reltol*analyticInt) ); 00281 long double absdiff = std::fabs(analyticInt - testInt); 00282 *outStream << "Cubature order " << std::setw(2) << std::left 00283 << l << " integrating " 00284 << "x^" << std::setw(2) << std::left << j << ":" 00285 << " " 00286 << std::scientific << std::setprecision(16) << testInt 00287 << " " << analyticInt 00288 << " " << std::setprecision(4) << absdiff << " " 00289 << "<?" << " " << abstol 00290 << "\n"; 00291 if (absdiff > abstol) { 00292 errorFlag++; 00293 *outStream << std::right << std::setw(104) 00294 << "^^^^---FAILURE!\n"; 00295 } 00296 } // end for j 00297 *outStream << "\n"; 00298 } // end for i 00299 } 00300 } // end for rule 00301 } 00302 catch (std::logic_error err) { 00303 *outStream << err.what() << "\n"; 00304 errorFlag = -1; 00305 }; 00306 00307 00308 if (errorFlag != 0) 00309 std::cout << "End Result: TEST FAILED\n"; 00310 else 00311 std::cout << "End Result: TEST PASSED\n"; 00312 00313 // reset format state of std::cout 00314 std::cout.copyfmt(oldFormatState); 00315 00316 return errorFlag; 00317 }
1.7.6.1