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