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