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