Intrepid: http://trilinos.sandia.gov/packages/docs/r10.0/packages/intrepid/test/Shared/ArrayTools/test

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 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Pavel Bochev (pbboche@sandia.gov) or
00025 //                    Denis Ridzal (dridzal@sandia.gov).
00026 //
00027 // ************************************************************************
00028 // @HEADER
00029 
00035 #include "Intrepid_ArrayTools.hpp"
00036 #include "Intrepid_FieldContainer.hpp"
00037 #include "Intrepid_RealSpaceTools.hpp"
00038 #include "Teuchos_oblackholestream.hpp"
00039 #include "Teuchos_RCP.hpp"
00040 #include "Teuchos_ScalarTraits.hpp"
00041 #include "Teuchos_GlobalMPISession.hpp"
00042 
00043 using namespace std;
00044 using namespace Intrepid;
00045 
00046 #define INTREPID_TEST_COMMAND( S )                                                                                  \
00047 {                                                                                                                   \
00048   try {                                                                                                             \
00049     S ;                                                                                                             \
00050   }                                                                                                                 \
00051   catch (std::logic_error err) {                                                                                    \
00052       *outStream << "Expected Error ----------------------------------------------------------------\n";            \
00053       *outStream << err.what() << '\n';                                                                             \
00054       *outStream << "-------------------------------------------------------------------------------" << "\n\n";    \
00055   };                                                                                                                \
00056 }
00057 
00058 
00059 int main(int argc, char *argv[]) {
00060 
00061   Teuchos::GlobalMPISession mpiSession(&argc, &argv);
00062 
00063   // This little trick lets us print to std::cout only if
00064   // a (dummy) command-line argument is provided.
00065   int iprint     = argc - 1;
00066   Teuchos::RCP<std::ostream> outStream;
00067   Teuchos::oblackholestream bhs; // outputs nothing
00068   if (iprint > 0)
00069     outStream = Teuchos::rcp(&std::cout, false);
00070   else
00071     outStream = Teuchos::rcp(&bhs, false);
00072 
00073   // Save the format state of the original std::cout.
00074   Teuchos::oblackholestream oldFormatState;
00075   oldFormatState.copyfmt(std::cout);
00076 
00077   *outStream \
00078   << "===============================================================================\n" \
00079   << "|                                                                             |\n" \
00080   << "|                       Unit Test (ArrayTools)                                |\n" \
00081   << "|                                                                             |\n" \
00082   << "|     1) Array operations: clone / scale                                      |\n" \
00083   << "|                                                                             |\n" \
00084   << "|  Questions? Contact  Pavel Bochev (pbboche@sandia.gov) or                   |\n" \
00085   << "|                      Denis Ridzal (dridzal@sandia.gov).                     |\n" \
00086   << "|                                                                             |\n" \
00087   << "|  Intrepid's website: http://trilinos.sandia.gov/packages/intrepid           |\n" \
00088   << "|  Trilinos website:   http://trilinos.sandia.gov                             |\n" \
00089   << "|                                                                             |\n" \
00090   << "===============================================================================\n";
00091 
00092 
00093   int errorFlag = 0;
00094 #ifdef HAVE_INTREPID_DEBUG
00095   int beginThrowNumber = TestForException_getThrowNumber();
00096   int endThrowNumber = beginThrowNumber + 21;
00097 #endif
00098 
00099   typedef ArrayTools art; 
00100   typedef RealSpaceTools<double> rst; 
00101 #ifdef HAVE_INTREPID_DEBUG
00102   ArrayTools atools;
00103 #endif
00104 
00105   *outStream \
00106   << "\n"
00107   << "===============================================================================\n"\
00108   << "| TEST 1: exceptions                                                          |\n"\
00109   << "===============================================================================\n";
00110 
00111   try{
00112 
00113 #ifdef HAVE_INTREPID_DEBUG
00114     FieldContainer<double> a_2(2);
00115     FieldContainer<double> a_9_2(9, 2);
00116     FieldContainer<double> a_10_2(10, 2);
00117     FieldContainer<double> a_10_3(10, 3);
00118     FieldContainer<double> a_10_2_2(10, 2, 2);
00119     FieldContainer<double> a_10_2_2_2(10, 2, 2, 2);
00120     FieldContainer<double> a_10_3_2_2(10, 3, 2, 2);
00121     FieldContainer<double> a_10_3_2(10, 3, 2);
00122     FieldContainer<double> a_2_2(2, 2);
00123     FieldContainer<double> a_2_3_2_2(2, 3, 2, 2);
00124     FieldContainer<double> a_2_2_2_2(2, 2, 2, 2);
00125     FieldContainer<double> a_10_2_2_2_2(10, 2, 2, 2, 2);
00126     FieldContainer<double> a_10_3_2_2_2(10, 3, 2, 2, 2);
00127     FieldContainer<double> a_10_2_3_2_2(10, 2, 3, 2, 2);
00128     FieldContainer<double> a_10_2_2_3_2(10, 2, 2, 3, 2);
00129     FieldContainer<double> a_10_2_2_2_3(10, 2, 2, 2, 3);
00130 
00131     *outStream << "-> cloneFields:\n";
00132     INTREPID_TEST_COMMAND( atools.cloneFields<double>(a_10_2_2_2, a_2) );
00133     INTREPID_TEST_COMMAND( atools.cloneFields<double>(a_10_2_2_2, a_10_2) );
00134     INTREPID_TEST_COMMAND( atools.cloneFields<double>(a_10_3_2, a_2_2) );
00135     INTREPID_TEST_COMMAND( atools.cloneFields<double>(a_10_2_2_2_2, a_2_3_2_2) );
00136     INTREPID_TEST_COMMAND( atools.cloneFields<double>(a_10_2_2_3_2, a_2_2_2_2) );
00137     INTREPID_TEST_COMMAND( atools.cloneFields<double>(a_10_2_2_2_3, a_2_2_2_2) );
00138     INTREPID_TEST_COMMAND( atools.cloneFields<double>(a_10_2_2, a_2_2) );
00139     INTREPID_TEST_COMMAND( atools.cloneFields<double>(a_10_2_2_2_2, a_2_2_2_2) );
00140 
00141     *outStream << "-> cloneScaleFields:\n";
00142     INTREPID_TEST_COMMAND( atools.cloneScaleFields<double>(a_10_2_2_2, a_2, a_2) );
00143     INTREPID_TEST_COMMAND( atools.cloneScaleFields<double>(a_10_2_2_2, a_10_2, a_2) );
00144     INTREPID_TEST_COMMAND( atools.cloneScaleFields<double>(a_10_2_2_2, a_10_2, a_10_2) );
00145     INTREPID_TEST_COMMAND( atools.cloneScaleFields<double>(a_10_2_2, a_9_2, a_10_2) );
00146     INTREPID_TEST_COMMAND( atools.cloneScaleFields<double>(a_10_2_2, a_10_3, a_10_2) );
00147     INTREPID_TEST_COMMAND( atools.cloneScaleFields<double>(a_10_3_2_2_2, a_10_3, a_2_2_2_2) );
00148     INTREPID_TEST_COMMAND( atools.cloneScaleFields<double>(a_10_2_3_2_2, a_10_2, a_2_2_2_2) );
00149     INTREPID_TEST_COMMAND( atools.cloneScaleFields<double>(a_10_2_2_3_2, a_10_2, a_2_2_2_2) );
00150     INTREPID_TEST_COMMAND( atools.cloneScaleFields<double>(a_10_2_2_2_3, a_10_2, a_2_2_2_2) );
00151     INTREPID_TEST_COMMAND( atools.cloneScaleFields<double>(a_10_2_2, a_10_2, a_2_2) );
00152     INTREPID_TEST_COMMAND( atools.cloneScaleFields<double>(a_10_2_2_2_2, a_10_2, a_2_2_2_2) );
00153 
00154     *outStream << "-> scaleFields:\n";
00155     INTREPID_TEST_COMMAND( atools.scaleFields<double>(a_10_2_2_2, a_2) );
00156     INTREPID_TEST_COMMAND( atools.scaleFields<double>(a_10_2, a_2_2) );
00157     INTREPID_TEST_COMMAND( atools.scaleFields<double>(a_10_2_2, a_2_2) );
00158     INTREPID_TEST_COMMAND( atools.scaleFields<double>(a_10_3_2, a_10_2) );
00159     INTREPID_TEST_COMMAND( atools.scaleFields<double>(a_10_3_2_2, a_10_2) );
00160     INTREPID_TEST_COMMAND( atools.scaleFields<double>(a_10_3_2_2_2, a_10_2) );
00161     INTREPID_TEST_COMMAND( atools.scaleFields<double>(a_10_2_2_2_2, a_10_2) );
00162 #endif
00163 
00164   }
00165   catch (std::logic_error err) {
00166     *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
00167     *outStream << err.what() << '\n';
00168     *outStream << "-------------------------------------------------------------------------------" << "\n\n";
00169     errorFlag = -1000;
00170   };
00171 
00172 #ifdef HAVE_INTREPID_DEBUG
00173   if (TestForException_getThrowNumber() != endThrowNumber)
00174     errorFlag++;
00175 #endif
00176 
00177   *outStream \
00178   << "\n"
00179   << "===============================================================================\n"\
00180   << "| TEST 2: correctness of math operations                                      |\n"\
00181   << "===============================================================================\n";
00182 
00183   outStream->precision(20);
00184 
00185   try {
00186       { // start scope
00187       *outStream << "\n************ Checking cloneFields ************\n";
00188 
00189       int c=5, p=9, f=7, d1=7, d2=13;
00190 
00191       FieldContainer<double> in_f_p(f, p);
00192       FieldContainer<double> in_f_p_d(f, p, d1);
00193       FieldContainer<double> in_f_p_d_d(f, p, d1, d2);
00194       FieldContainer<double> in_c_f_p(c, f, p);
00195       FieldContainer<double> in_c_f_p_d(c, f, p, d1);
00196       FieldContainer<double> in_c_f_p_d_d(c, f, p, d1, d2);
00197       FieldContainer<double> data_c_p_one(c, p);
00198       FieldContainer<double> out_c_f_p(c, f, p);
00199       FieldContainer<double> out_c_f_p_d(c, f, p, d1);
00200       FieldContainer<double> out_c_f_p_d_d(c, f, p, d1, d2);
00201       double zero = INTREPID_TOL*100.0;
00202 
00203       // fill with random numbers
00204       for (int i=0; i<in_f_p.size(); i++) {
00205         in_f_p[i] = Teuchos::ScalarTraits<double>::random();
00206       }
00207       for (int i=0; i<in_f_p_d.size(); i++) {
00208         in_f_p_d[i] = Teuchos::ScalarTraits<double>::random();
00209       }
00210       for (int i=0; i<in_f_p_d_d.size(); i++) {
00211         in_f_p_d_d[i] = Teuchos::ScalarTraits<double>::random();
00212       }
00213       for (int i=0; i<data_c_p_one.size(); i++) {
00214         data_c_p_one[i] = 1.0;
00215       }
00216 
00217       art::cloneFields<double>(out_c_f_p, in_f_p);
00218       art::scalarMultiplyDataField<double>(in_c_f_p, data_c_p_one, in_f_p);
00219       rst::subtract(&out_c_f_p[0], &in_c_f_p[0], out_c_f_p.size());
00220       if (rst::vectorNorm(&out_c_f_p[0], out_c_f_p.size(), NORM_ONE) > zero) {
00221         *outStream << "\n\nINCORRECT cloneFields (1): check multiplyScalarData vs. cloneFields\n\n";
00222         errorFlag = -1000;
00223       }
00224       art::cloneFields<double>(out_c_f_p_d, in_f_p_d);
00225       art::scalarMultiplyDataField<double>(in_c_f_p_d, data_c_p_one, in_f_p_d);
00226       rst::subtract(&out_c_f_p_d[0], &in_c_f_p_d[0], out_c_f_p_d.size());
00227       if (rst::vectorNorm(&out_c_f_p_d[0], out_c_f_p_d.size(), NORM_ONE) > zero) {
00228         *outStream << "\n\nINCORRECT cloneFields (2): check multiplyScalarData vs. cloneFields\n\n";
00229         errorFlag = -1000;
00230       }
00231       art::cloneFields<double>(out_c_f_p_d_d, in_f_p_d_d);
00232       art::scalarMultiplyDataField<double>(in_c_f_p_d_d, data_c_p_one, in_f_p_d_d);
00233       rst::subtract(&out_c_f_p_d_d[0], &in_c_f_p_d_d[0], out_c_f_p_d_d.size());
00234       if (rst::vectorNorm(&out_c_f_p_d_d[0], out_c_f_p_d_d.size(), NORM_ONE) > zero) {
00235         *outStream << "\n\nINCORRECT cloneFields (3): check multiplyScalarData vs. cloneFields\n\n";
00236         errorFlag = -1000;
00237       }
00238       } // end scope
00239 
00240       { // start scope
00241       *outStream << "\n************ Checking cloneScaleFields ************\n";
00242       int c=5, p=9, f=7, d1=7, d2=13;
00243 
00244       FieldContainer<double> in_f_p(f, p);
00245       FieldContainer<double> in_f_p_d(f, p, d1);
00246       FieldContainer<double> in_f_p_d_d(f, p, d1, d2);
00247       FieldContainer<double> data_c_f(c, f);
00248       FieldContainer<double> datainv_c_f(c, f);
00249       FieldContainer<double> c_f_p_one(c, f, p);
00250       FieldContainer<double> c_f_p_d_one(c, f, p, d1);
00251       FieldContainer<double> c_f_p_d_d_one(c, f, p, d1, d2);
00252       FieldContainer<double> out_c_f_p(c, f, p);
00253       FieldContainer<double> outi_c_f_p(c, f, p);
00254       FieldContainer<double> out_c_f_p_d(c, f, p, d1);
00255       FieldContainer<double> outi_c_f_p_d(c, f, p, d1);
00256       FieldContainer<double> out_c_f_p_d_d(c, f, p, d1, d2);
00257       FieldContainer<double> outi_c_f_p_d_d(c, f, p, d1, d2);
00258       double zero = INTREPID_TOL*100.0;
00259 
00260       // fill with 1's
00261       for (int i=0; i<in_f_p.size(); i++) {
00262         in_f_p[i] = 1.0;
00263       }
00264       for (int i=0; i<in_f_p_d.size(); i++) {
00265         in_f_p_d[i] = 1.0;
00266       }
00267       for (int i=0; i<in_f_p_d_d.size(); i++) {
00268         in_f_p_d_d[i] = 1.0;
00269       }
00270       for (int i=0; i<c_f_p_one.size(); i++) {
00271         c_f_p_one[i] = 1.0;
00272       }
00273       for (int i=0; i<c_f_p_d_one.size(); i++) {
00274         c_f_p_d_one[i] = 1.0;
00275       }
00276       for (int i=0; i<c_f_p_d_d_one.size(); i++) {
00277         c_f_p_d_d_one[i] = 1.0;
00278       }
00279       // fill with random numbers
00280       for (int i=0; i<data_c_f.size(); i++) {
00281         data_c_f[i] = Teuchos::ScalarTraits<double>::random();
00282         datainv_c_f[i] = 1.0 / data_c_f[i];
00283       }
00284 
00285       art::cloneScaleFields<double>(out_c_f_p, data_c_f, in_f_p);
00286       art::cloneScaleFields<double>(outi_c_f_p, datainv_c_f, in_f_p);
00287       for (int i=0; i<out_c_f_p.size(); i++) {
00288         out_c_f_p[i] *= outi_c_f_p[i];
00289       }
00290       rst::subtract(&out_c_f_p[0], &c_f_p_one[0], out_c_f_p.size());
00291       if (rst::vectorNorm(&out_c_f_p[0], out_c_f_p.size(), NORM_ONE) > zero) {
00292         *outStream << "\n\nINCORRECT cloneScaleValue (1): check scalar inverse property\n\n";
00293         errorFlag = -1000;
00294       }
00295 
00296       art::cloneScaleFields<double>(out_c_f_p_d, data_c_f, in_f_p_d);
00297       art::cloneScaleFields<double>(outi_c_f_p_d, datainv_c_f, in_f_p_d);
00298       for (int i=0; i<out_c_f_p_d.size(); i++) {
00299         out_c_f_p_d[i] *= outi_c_f_p_d[i];
00300       }
00301       rst::subtract(&out_c_f_p_d[0], &c_f_p_d_one[0], out_c_f_p_d.size());
00302       if (rst::vectorNorm(&out_c_f_p_d[0], out_c_f_p_d.size(), NORM_ONE) > zero) {
00303         *outStream << "\n\nINCORRECT cloneScaleValue (2): check scalar inverse property\n\n";
00304         errorFlag = -1000;
00305       }
00306 
00307       art::cloneScaleFields<double>(out_c_f_p_d_d, data_c_f, in_f_p_d_d);
00308       art::cloneScaleFields<double>(outi_c_f_p_d_d, datainv_c_f, in_f_p_d_d);
00309       for (int i=0; i<out_c_f_p_d_d.size(); i++) {
00310         out_c_f_p_d_d[i] *= outi_c_f_p_d_d[i];
00311       }
00312       rst::subtract(&out_c_f_p_d_d[0], &c_f_p_d_d_one[0], out_c_f_p_d_d.size());
00313       if (rst::vectorNorm(&out_c_f_p_d_d[0], out_c_f_p_d_d.size(), NORM_ONE) > zero) {
00314         *outStream << "\n\nINCORRECT cloneScaleValue (3): check scalar inverse property\n\n";
00315         errorFlag = -1000;
00316       }
00317       } // end scope
00318 
00319       { // start scope
00320       *outStream << "\n************ Checking scaleFields ************\n";
00321       int c=5, p=9, f=7, d1=7, d2=13;
00322 
00323       FieldContainer<double> data_c_f(c, f);
00324       FieldContainer<double> datainv_c_f(c, f);
00325       FieldContainer<double> out_c_f_p(c, f, p);
00326       FieldContainer<double> outi_c_f_p(c, f, p);
00327       FieldContainer<double> out_c_f_p_d(c, f, p, d1);
00328       FieldContainer<double> outi_c_f_p_d(c, f, p, d1);
00329       FieldContainer<double> out_c_f_p_d_d(c, f, p, d1, d2);
00330       FieldContainer<double> outi_c_f_p_d_d(c, f, p, d1, d2);
00331       double zero = INTREPID_TOL*100.0;
00332 
00333       // fill with random numbers
00334       for (int i=0; i<out_c_f_p.size(); i++) {
00335         out_c_f_p[i] = Teuchos::ScalarTraits<double>::random();
00336         outi_c_f_p[i] = out_c_f_p[i];
00337       }
00338       for (int i=0; i<out_c_f_p_d.size(); i++) {
00339         out_c_f_p_d[i] = Teuchos::ScalarTraits<double>::random();
00340         outi_c_f_p_d[i] = out_c_f_p_d[i];
00341       }
00342       for (int i=0; i<out_c_f_p_d_d.size(); i++) {
00343         out_c_f_p_d_d[i] = Teuchos::ScalarTraits<double>::random();
00344         outi_c_f_p_d_d[i] = out_c_f_p_d_d[i];
00345       }
00346       for (int i=0; i<data_c_f.size(); i++) {
00347         data_c_f[i] = Teuchos::ScalarTraits<double>::random();
00348         datainv_c_f[i] = 1.0 / data_c_f[i];
00349       }
00350 
00351       art::scaleFields<double>(out_c_f_p, data_c_f);
00352       art::scaleFields<double>(out_c_f_p, datainv_c_f);
00353       rst::subtract(&out_c_f_p[0], &outi_c_f_p[0], out_c_f_p.size());
00354       if (rst::vectorNorm(&out_c_f_p[0], out_c_f_p.size(), NORM_ONE) > zero) {
00355         *outStream << "\n\nINCORRECT scaleValue (1): check scalar inverse property\n\n";
00356         errorFlag = -1000;
00357       }
00358 
00359       art::scaleFields<double>(out_c_f_p_d, data_c_f);
00360       art::scaleFields<double>(out_c_f_p_d, datainv_c_f);
00361       rst::subtract(&out_c_f_p_d[0], &outi_c_f_p_d[0], out_c_f_p_d.size());
00362       if (rst::vectorNorm(&out_c_f_p_d[0], out_c_f_p_d.size(), NORM_ONE) > zero) {
00363         *outStream << "\n\nINCORRECT scaleValue (2): check scalar inverse property\n\n";
00364         errorFlag = -1000;
00365       }
00366 
00367       art::scaleFields<double>(out_c_f_p_d_d, data_c_f);
00368       art::scaleFields<double>(out_c_f_p_d_d, datainv_c_f);
00369       rst::subtract(&out_c_f_p_d_d[0], &outi_c_f_p_d_d[0], out_c_f_p_d_d.size());
00370       if (rst::vectorNorm(&out_c_f_p_d_d[0], out_c_f_p_d_d.size(), NORM_ONE) > zero) {
00371         *outStream << "\n\nINCORRECT cloneScaleValue (3): check scalar inverse property\n\n";
00372         errorFlag = -1000;
00373       }
00374       } // end scope
00375 
00376       /******************************************/
00377       *outStream << "\n";
00378   }
00379   catch (std::logic_error err) {
00380     *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
00381     *outStream << err.what() << '\n';
00382     *outStream << "-------------------------------------------------------------------------------" << "\n\n";
00383     errorFlag = -1000;
00384   };
00385 
00386 
00387   if (errorFlag != 0)
00388     std::cout << "End Result: TEST FAILED\n";
00389   else
00390     std::cout << "End Result: TEST PASSED\n";
00391 
00392   // reset format state of std::cout
00393   std::cout.copyfmt(oldFormatState);
00394 
00395   return errorFlag;
00396 }
http://trilinos.sandia.gov/packages/docs/r10.0/packages/intrepid/test/Shared/ArrayTools/test_05.cpp
