|
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_CubatureTensorSorted.hpp" 00052 #include "Intrepid_CubatureTensorSorted.hpp" 00053 //#include "Intrepid_CubatureLineSorted.hpp" 00054 #include "Intrepid_Utils.hpp" 00055 #include "Teuchos_oblackholestream.hpp" 00056 #include "Teuchos_RCP.hpp" 00057 #include "Teuchos_GlobalMPISession.hpp" 00058 00059 using namespace Intrepid; 00060 00061 /* 00062 Computes integrals of monomials over a given reference cell. 00063 */ 00064 long double evalQuad(std::vector<int> power, 00065 int dimension, std::vector<int> order, 00066 std::vector<EIntrepidBurkardt> rule) { 00067 00068 CubatureTensorSorted<long double> lineCub(dimension,order,rule,false); 00069 int size = lineCub.getNumPoints(); 00070 FieldContainer<long double> cubPoints(size,dimension); 00071 FieldContainer<long double> cubWeights(size); 00072 lineCub.getCubature(cubPoints,cubWeights); 00073 00074 int mid = size/2; 00075 long double Q = 0.0; 00076 if (size%2) { 00077 Q = cubWeights(mid); 00078 for (int i=0; i<dimension; i++) { 00079 Q *= powl(cubPoints(mid,i),power[i]); 00080 } 00081 } 00082 00083 for (int i=0; i<mid; i++) { 00084 long double value1 = cubWeights(i); 00085 long double value2 = cubWeights(size-i-1); 00086 for (int j=0; j<dimension; j++) { 00087 value1 *= powl(cubPoints(i,j),power[j]); 00088 value2 *= powl(cubPoints(size-i-1,j),power[j]); 00089 } 00090 Q += value1+value2; 00091 } 00092 return Q; 00093 } 00094 00095 long double evalInt(int dimension, std::vector<int> power, 00096 std::vector<EIntrepidBurkardt> rule) { 00097 long double I = 1.0; 00098 00099 for (int i=0; i<dimension; i++) { 00100 if (rule[i]==BURK_CLENSHAWCURTIS||rule[i]==BURK_FEJER2|| 00101 rule[i]==BURK_LEGENDRE||rule[i]==BURK_PATTERSON || 00102 rule[i]==BURK_TRAPEZOIDAL) { 00103 if (power[i]%2) 00104 I *= 0.0; 00105 else 00106 I *= 2.0/((long double)power[i]+1.0); 00107 } 00108 else if (rule[i]==BURK_LAGUERRE) { 00109 I *= tgammal((long double)(power[i]+1)); 00110 } 00111 else if (rule[i]==BURK_CHEBYSHEV1) { 00112 long double bot, top; 00113 if (!(power[i]%2)) { 00114 top = 1; bot = 1; 00115 for (int j=2;j<=power[i];j+=2) { 00116 top *= (long double)(j-1); 00117 bot *= (long double)j; 00118 } 00119 I *= M_PI*top/bot; 00120 } 00121 else { 00122 I *= 0.0; 00123 } 00124 } 00125 else if (rule[i]==BURK_CHEBYSHEV2) { 00126 long double bot, top; 00127 if (!(power[i]%2)) { 00128 top = 1; bot = 1; 00129 for (int j=2;j<=power[i];j+=2) { 00130 top *= (long double)(j-1); 00131 bot *= (long double)j; 00132 } 00133 bot *= (long double)(power[i]+2); 00134 I *= M_PI*top/bot; 00135 } 00136 else { 00137 I *= 0.0; 00138 } 00139 } 00140 else if (rule[i]==BURK_HERMITE||rule[i]==BURK_GENZKEISTER) { 00141 if (power[i]%2) { 00142 I *= 0.0; 00143 } 00144 else { 00145 long double value = 1.0; 00146 if ((power[i]-1)>=1) { 00147 int n_copy = power[i]-1; 00148 while (1<n_copy) { 00149 value *= (long double)n_copy; 00150 n_copy -= 2; 00151 } 00152 } 00153 I *= value*sqrt(M_PI)/powl(2.0,(long double)power[i]/2.0); 00154 } 00155 } 00156 } 00157 return I; 00158 } 00159 00160 int main(int argc, char *argv[]) { 00161 00162 Teuchos::GlobalMPISession mpiSession(&argc, &argv); 00163 00164 // This little trick lets us print to std::cout only if 00165 // a (dummy) command-line argument is provided. 00166 int iprint = argc - 1; 00167 Teuchos::RCP<std::ostream> outStream; 00168 Teuchos::oblackholestream bhs; // outputs nothing 00169 if (iprint > 0) 00170 outStream = Teuchos::rcp(&std::cout, false); 00171 else 00172 outStream = Teuchos::rcp(&bhs, false); 00173 00174 // Save the format state of the original std::cout. 00175 Teuchos::oblackholestream oldFormatState; 00176 oldFormatState.copyfmt(std::cout); 00177 00178 *outStream \ 00179 << "===============================================================================\n" \ 00180 << "| |\n" \ 00181 << "| Unit Test (CubatureTensorSorted) |\n" \ 00182 << "| |\n" \ 00183 << "| 1) Computing integrals of monomials in 2D |\n" \ 00184 << "| |\n" \ 00185 << "| Questions? Contact Drew Kouri (dpkouri@sandia.gov) or |\n" \ 00186 << "| Denis Ridzal (dridzal@sandia.gov). |\n" \ 00187 << "| |\n" \ 00188 << "| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \ 00189 << "| Trilinos website: http://trilinos.sandia.gov |\n" \ 00190 << "| |\n" \ 00191 << "===============================================================================\n"\ 00192 << "| TEST 12: integrals of monomials in 2D - Anisotropic but no growth rules |\n"\ 00193 << "===============================================================================\n"; 00194 00195 00196 // internal variables: 00197 int dimension = 2; 00198 int errorFlag = 0; 00199 long double reltol = 1.0e+05*INTREPID_TOL; 00200 int maxDeg = 0; 00201 long double analyticInt = 0; 00202 long double testInt = 0; 00203 int maxOrder = 14; 00204 std::vector<int> power(2,0); 00205 std::vector<EIntrepidBurkardt> rule1(2,BURK_CLENSHAWCURTIS); 00206 std::vector<int> order(2,0); 00207 00208 *outStream << "\nIntegrals of monomials on a reference line (edge):\n"; 00209 // compute and compare integrals 00210 try { 00211 for (EIntrepidBurkardt rule=BURK_CHEBYSHEV1;rule<=BURK_LAGUERRE;rule++) { 00212 *outStream << "Testing " << EIntrepidBurkardtToString(rule) << "\n"; 00213 // compute integrals 00214 if (rule==BURK_HERMITE) 00215 maxOrder = 8; 00216 else if (rule==BURK_TRAPEZOIDAL) 00217 maxOrder = 2; 00218 else 00219 maxOrder = 9; 00220 00221 rule1[0] = rule; rule1[1] = rule; 00222 if (rule!=BURK_PATTERSON&&rule!=BURK_GENZKEISTER) { 00223 for (int i=1; i <= maxOrder; i++) { 00224 if ( rule==BURK_CHEBYSHEV1 || 00225 rule==BURK_CHEBYSHEV2 || 00226 rule==BURK_LEGENDRE || 00227 rule==BURK_LAGUERRE || 00228 rule==BURK_HERMITE ) 00229 maxDeg = 2*i-1; 00230 else if ( rule==BURK_CLENSHAWCURTIS || 00231 rule==BURK_FEJER2 || 00232 rule==BURK_TRAPEZOIDAL ) 00233 maxDeg = i-1; 00234 00235 order[0] = i; order[1] = i; 00236 for (int j=0; j <= maxDeg; j++) { 00237 power[0] = j; 00238 for (int k=0; k <= maxDeg; k++) { 00239 power[1] = k; 00240 analyticInt = evalInt(dimension, power, rule1); 00241 testInt = evalQuad(power,dimension,order,rule1); 00242 00243 long double abstol = (analyticInt == 0.0 ? reltol : std::fabs(reltol*analyticInt) ); 00244 long double absdiff = std::fabs(analyticInt - testInt); 00245 *outStream << "Cubature order " << std::setw(2) 00246 << std::left << i << " integrating " 00247 << "x^" << std::setw(2) << std::left << j 00248 << "y^" << std::setw(2) << std::left 00249 << k << ":" << " " << std::scientific 00250 << std::setprecision(16) << testInt 00251 << " " << analyticInt << " " 00252 << std::setprecision(4) << absdiff << " " 00253 << "<?" << " " << abstol << "\n"; 00254 if (absdiff > abstol) { 00255 errorFlag++; 00256 *outStream << std::right << std::setw(104) 00257 << "^^^^---FAILURE!\n"; 00258 } 00259 } // end for k 00260 *outStream << "\n"; 00261 } // end for j 00262 *outStream << "\n"; 00263 } // end for i 00264 } 00265 else if (rule==BURK_PATTERSON) { 00266 for (int i=0; i < 3; i++) { 00267 int l = (int)std::pow(2.0,(double)i+1.0)-1; 00268 if (i==0) 00269 maxDeg = 1; 00270 else 00271 maxDeg = (int)(1.5*(double)l+0.5); 00272 00273 order[0] = l; order[1] = l; 00274 for (int j=0; j <= maxDeg; j++) { 00275 power[0] = j; 00276 for (int k=0; k <= maxDeg; k++) { 00277 power[1] = k; 00278 analyticInt = evalInt(dimension, power, rule1); 00279 testInt = evalQuad(power,dimension,order,rule1); 00280 00281 long double abstol = (analyticInt == 0.0 ? reltol : std::fabs(reltol*analyticInt) ); 00282 long double absdiff = std::fabs(analyticInt - testInt); 00283 *outStream << "Cubature order " << std::setw(2) 00284 << std::left << l << " integrating " 00285 << "x^" << std::setw(2) << std::left << j 00286 << "y^" << std::setw(2) << std::left 00287 << k << ":" << " " << std::scientific 00288 << std::setprecision(16) << testInt 00289 << " " << analyticInt << " " 00290 << std::setprecision(4) << absdiff << " " 00291 << "<?" << " " << abstol << "\n"; 00292 if (absdiff > abstol) { 00293 errorFlag++; 00294 *outStream << std::right << std::setw(104) 00295 << "^^^^---FAILURE!\n"; 00296 } 00297 } // end for k 00298 *outStream << "\n"; 00299 } // end for j 00300 *outStream << "\n"; 00301 } // end for i 00302 } 00303 else if (rule==BURK_GENZKEISTER) { 00304 int o_ghk[8] = {1,3,9,19,35,37,41,43}; 00305 for (int i=0; i < 3; i++) { 00306 int l = o_ghk[i]; 00307 if (i==0) 00308 maxDeg = 1; 00309 else 00310 maxDeg = (int)(1.5*(double)l+0.5); 00311 00312 order[0] = l; order[1] = l; 00313 for (int j=0; j <= maxDeg; j++) { 00314 power[0] = j; 00315 for (int k=0; k <= maxDeg; k++) { 00316 power[1] = k; 00317 analyticInt = evalInt(dimension, power, rule1); 00318 testInt = evalQuad(power,dimension,order,rule1); 00319 00320 long double abstol = (analyticInt == 0.0 ? reltol : std::fabs(reltol*analyticInt) ); 00321 long double absdiff = std::fabs(analyticInt - testInt); 00322 *outStream << "Cubature order " << std::setw(2) 00323 << std::left << l << " integrating " 00324 << "x^" << std::setw(2) << std::left << j 00325 << "y^" << std::setw(2) << std::left 00326 << k << ":" << " " << std::scientific 00327 << std::setprecision(16) << testInt 00328 << " " << analyticInt << " " 00329 << std::setprecision(4) << absdiff << " " 00330 << "<?" << " " << abstol << "\n"; 00331 if (absdiff > abstol) { 00332 errorFlag++; 00333 *outStream << std::right << std::setw(104) 00334 << "^^^^---FAILURE!\n"; 00335 } 00336 } // end for k 00337 *outStream << "\n"; 00338 } // end for j 00339 *outStream << "\n"; 00340 } // end for i 00341 } 00342 } // end for rule 00343 } 00344 catch (std::logic_error err) { 00345 *outStream << err.what() << "\n"; 00346 errorFlag = -1; 00347 }; 00348 00349 00350 if (errorFlag != 0) 00351 std::cout << "End Result: TEST FAILED\n"; 00352 else 00353 std::cout << "End Result: TEST PASSED\n"; 00354 00355 // reset format state of std::cout 00356 std::cout.copyfmt(oldFormatState); 00357 00358 return errorFlag; 00359 }
1.7.6.1