Intrepid
/usr/src/RPM/BUILD/trilinos10-10.6.4/packages/intrepid/test/Shared/FieldContainer/test_03.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 // 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 
00030 
00036 #include "Intrepid_FieldContainer.hpp"
00037 #include "Intrepid_RealSpaceTools.hpp"
00038 #include "Shards_Array.hpp"
00039 #include "Teuchos_oblackholestream.hpp"
00040 #include "Teuchos_RCP.hpp"
00041 #include "Teuchos_GlobalMPISession.hpp"
00042 
00043 
00044 using namespace Intrepid;
00045 
00046 SHARDS_ARRAY_DIM_TAG_SIMPLE_DECLARATION( Cell )
00047 SHARDS_ARRAY_DIM_TAG_SIMPLE_IMPLEMENTATION( Cell )
00048 
00049 SHARDS_ARRAY_DIM_TAG_SIMPLE_DECLARATION( Field )
00050 SHARDS_ARRAY_DIM_TAG_SIMPLE_IMPLEMENTATION( Field )
00051 
00052 SHARDS_ARRAY_DIM_TAG_SIMPLE_DECLARATION( Point )
00053 SHARDS_ARRAY_DIM_TAG_SIMPLE_IMPLEMENTATION( Point )
00054 
00055 SHARDS_ARRAY_DIM_TAG_SIMPLE_DECLARATION( Dim )
00056 SHARDS_ARRAY_DIM_TAG_SIMPLE_IMPLEMENTATION( Dim )
00057 
00058 #define INTREPID_TEST_COMMAND( S )                                                                                  \
00059 {                                                                                                                   \
00060   try {                                                                                                             \
00061     S ;                                                                                                             \
00062   }                                                                                                                 \
00063   catch (std::logic_error err) {                                                                                    \
00064       *outStream << "Expected Error ----------------------------------------------------------------\n";            \
00065       *outStream << err.what() << '\n';                                                                             \
00066       *outStream << "-------------------------------------------------------------------------------" << "\n\n";    \
00067   };                                                                                                                \
00068 }
00069 
00070 
00071 int main(int argc, char *argv[]) {
00072 
00073   Teuchos::GlobalMPISession mpiSession(&argc, &argv);
00074   
00075   // This little trick lets us print to cout only if a (dummy) command-line argument is provided.
00076   int iprint     = argc - 1;
00077   
00078   Teuchos::RCP<std::ostream> outStream;
00079   Teuchos::oblackholestream bhs; // outputs nothing
00080   
00081   if (iprint > 0)
00082     outStream = Teuchos::rcp(&std::cout, false);
00083   else
00084     outStream = Teuchos::rcp(&bhs, false);
00085   
00086   // Save the format state of the original cout .
00087   Teuchos::oblackholestream oldFormatState;
00088   oldFormatState.copyfmt(std::cout);
00089   
00090   *outStream  \
00091     << "===============================================================================\n" \
00092     << "|                                                                             |\n" \
00093     << "|                           Unit Test FieldContainer                          |\n" \
00094     << "|                                                                             |\n" \
00095     << "|     1) Testing usage of various constructors / wrappers                     |\n" \
00096     << "|     2) Testing usage of resize                                              |\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   // Set error flag.
00108   int errorFlag  = 0;
00109 
00110   double zero = INTREPID_TOL;
00111 
00112   try {
00113   
00114     // Dimensions array.
00115     Teuchos::Array<int> dimensions;
00116   
00117     *outStream << "\n" \
00118       << "===============================================================================\n"\
00119       << "| TEST 1: Constructors / Wrappers for a particular rank-4 container           |\n"\
00120       << "===============================================================================\n\n";
00121 
00122     { // start variable scope
00123 
00124       // Initialize dimensions for rank-4 multi-index value
00125       dimensions.resize(4);
00126       dimensions[0] = 5;
00127       dimensions[1] = 3;
00128       dimensions[2] = 2;
00129       dimensions[3] = 7;
00130 
00131       // Allocate data through a Teuchos::Array, initialized to 0.
00132       Teuchos::Array<double> data(dimensions[0]*dimensions[1]*dimensions[2]*dimensions[3]);
00133     
00134       // Create a FieldContainer using a deep copy via Teuchos::Array.
00135       FieldContainer<double> fc_array(dimensions, data);
00136 
00137       // modify the (1,1,1,1) entry
00138       fc_array(1,1,1,1) = 1.0;
00139       // verify that the data array has NOT changed
00140       if (std::abs(data[dimensions[1]*dimensions[2]*dimensions[3] + dimensions[2]*dimensions[3] + dimensions[3] + 1]) > zero) {
00141         *outStream << "\n\nError in constructor using Array (ArrayView) and deep copy.\n\n";
00142         errorFlag = -1000;
00143       }
00144 
00145       // test getData access function
00146       if (std::abs((fc_array.getData())[dimensions[1]*dimensions[2]*dimensions[3] +
00147                                         dimensions[2]*dimensions[3] + dimensions[3] + 1] - fc_array(1,1,1,1)) > zero) {
00148         *outStream << "\n\nError in getData() member of FieldContainer.\n\n";
00149         errorFlag = -1000;
00150       }
00151 
00152       // Create a FieldContainer using a deep copy via Teuchos::ArrayRCP.
00153       Teuchos::RCP< Teuchos::Array<double>  > rcp_to_data = rcpFromRef(data); // first create RCP
00154       Teuchos::ArrayRCP<double> arrayrcp = Teuchos::arcp(rcp_to_data);        // now create ArrayRCP
00155       FieldContainer<double> fc_arrayrcp_deep(dimensions, arrayrcp());        // IMPORTANT!!!: use '()' after arrayrcp
00156                                                                               //    for direct conversion to ArrayView
00157       // modify the (1,1,1,1) entry
00158       fc_arrayrcp_deep(1,1,1,1) = 1.0;
00159       // verify that the data array has NOT changed
00160       if (std::abs(data[dimensions[1]*dimensions[2]*dimensions[3] + dimensions[2]*dimensions[3] + dimensions[3] + 1]) > zero) {
00161         *outStream << "\n\nError in constructor using ArrayRCP (ArrayView) and deep copy.\n\n";
00162         errorFlag = -1000;
00163       }
00164 
00165       // Create a FieldContainer using a shallow copy via Teuchos::ArrayRCP.
00166       FieldContainer<double> fc_arrayrcp_shallow(dimensions, arrayrcp);
00167       // modify the (1,1,1,1) entry
00168       fc_arrayrcp_shallow(1,1,1,1) = 1.0;
00169       // verify that the data array HAS changed
00170       if (std::abs(data[dimensions[1]*dimensions[2]*dimensions[3] + dimensions[2]*dimensions[3] + dimensions[3] + 1] - fc_arrayrcp_shallow(1,1,1,1)) > zero) {
00171         *outStream << "\n\nError in constructor using ArrayRCP and shallow copy.\n\n";
00172         errorFlag = -1000;
00173       }
00174 
00175       // Create a FieldContainer using a deep copy via Scalar*.
00176       FieldContainer<double> fc_scalarptr_deep(dimensions, data.getRawPtr(), true);
00177       // modify the (1,1,1,1) entry
00178       fc_scalarptr_deep(1,1,1,1) = 2.0;
00179       // verify that the data array has NOT changed
00180       if (std::abs(data[dimensions[1]*dimensions[2]*dimensions[3] + dimensions[2]*dimensions[3] + dimensions[3] + 1] - 1.0) > zero) {
00181         *outStream << "\n\nError in constructor using Scalar* and deep copy.\n\n";
00182         errorFlag = -1000;
00183       }
00184 
00185       // Create a FieldContainer using a shallow copy via Scalar*.
00186       FieldContainer<double> fc_scalarptr_shallow(dimensions, data.getRawPtr());
00187       // modify the (1,1,1,1) entry
00188       fc_scalarptr_shallow(1,1,1,1) = 2.0;
00189       // verify that the data array HAS changed
00190       if (std::abs(data[dimensions[1]*dimensions[2]*dimensions[3] + dimensions[2]*dimensions[3] + dimensions[3] + 1] - fc_scalarptr_shallow(1,1,1,1)) > zero) {
00191         *outStream << "\n\nError in constructor using Scalar* and shallow copy.\n\n";
00192         errorFlag = -1000;
00193       }
00194 
00195       // Create a FieldContainer using a deep copy via shards::Array.
00196       shards::Array<double,shards::NaturalOrder,Cell,Field,Point,Dim> shards_array(data.getRawPtr(),dimensions[0],dimensions[1],dimensions[2],dimensions[3]);
00197       FieldContainer<double> fc_shards_deep(shards_array, true);
00198       // modify the (1,1,1,1) entry
00199       fc_shards_deep(1,1,1,1) = 3.0;
00200       // verify that the data array has NOT changed
00201       if (std::abs(data[dimensions[1]*dimensions[2]*dimensions[3] + dimensions[2]*dimensions[3] + dimensions[3] + 1] - 2.0) > zero) {
00202         *outStream << "\n\nError in constructor using shards::Array and deep copy.\n\n";
00203         errorFlag = -1000;
00204       }
00205 
00206       // Create a FieldContainer using a shallow copy via shards::Array.
00207       FieldContainer<double> fc_shards_shallow(shards_array);
00208       // modify the (1,1,1,1) entry
00209       fc_shards_shallow(1,1,1,1) = 3.0;
00210       // verify that the data array HAS changed
00211       if (std::abs(data[dimensions[1]*dimensions[2]*dimensions[3] + dimensions[2]*dimensions[3] + dimensions[3] + 1] - fc_shards_shallow(1,1,1,1)) > zero) {
00212         *outStream << "\n\nError in constructor using shards::Array and shallow copy.\n\n";
00213         errorFlag = -1000;
00214       }
00215 
00216     } // end variable scope
00217 
00218 
00219     *outStream << "\n" \
00220       << "===============================================================================\n"\
00221       << "| TEST 1 cont'd: Run through constructors / wrappers of various ranks         |\n"\
00222       << "===============================================================================\n\n";
00223 
00224     for (int rank=0; rank<10; rank++) {
00225       dimensions.resize(rank);
00226       int total_size = 1;
00227       if (rank == 0) {
00228         total_size = 0;
00229       }
00230       for (int dim=0; dim<rank; dim++) {
00231         dimensions[dim] = 2;
00232         total_size *= dimensions[dim];
00233       }
00234 
00235       // Allocate data through a Teuchos::Array, initialized to 0.
00236       Teuchos::Array<double> data(total_size);
00237 
00238       // Create a FieldContainer using a deep copy via Teuchos::Array.
00239       FieldContainer<double> fc_array(dimensions, data);
00240       // Create a FieldContainer using a deep copy via Teuchos::ArrayRCP.
00241       Teuchos::RCP< Teuchos::Array<double>  > rcp_to_data = rcpFromRef(data); // first create RCP
00242       Teuchos::ArrayRCP<double> arrayrcp = Teuchos::arcp(rcp_to_data);        // now create ArrayRCP
00243       FieldContainer<double> fc_arrayrcp_deep(dimensions, arrayrcp());        // IMPORTANT!!!: use '()' after arrayrcp
00244                                                                               //    for direct conversion to ArrayView
00245       // Create a FieldContainer using a shallow copy via Teuchos::ArrayRCP.
00246       FieldContainer<double> fc_arrayrcp_shallow(dimensions, arrayrcp);
00247       // Create a FieldContainer using a deep copy via Scalar*.
00248       FieldContainer<double> fc_scalarptr_deep(dimensions, data.getRawPtr(), true);
00249       // Create a FieldContainer using a shallow copy via Scalar*.
00250       FieldContainer<double> fc_scalarptr_shallow(dimensions, data.getRawPtr());
00251     }
00252 
00253     { // start variable scope
00254       // Create FieldContainers using a deep copy via shards::Array.
00255       Teuchos::Array<double> data(2*2*2*2*2*2);
00256       shards::Array<double,shards::NaturalOrder,Cell> shards_array_c(&data[0],2);
00257       shards::Array<double,shards::NaturalOrder,Cell,Field> shards_array_cf(&data[0],2,2);
00258       shards::Array<double,shards::NaturalOrder,Cell,Field,Point> shards_array_cfp(&data[0],2,2,2);
00259       shards::Array<double,shards::NaturalOrder,Cell,Field,Point,Dim> shards_array_cfpd(&data[0],2,2,2,2);
00260       shards::Array<double,shards::NaturalOrder,Cell,Field,Point,Dim,Dim> shards_array_cfpdd(&data[0],2,2,2,2,2);
00261       shards::Array<double,shards::NaturalOrder,Cell,Field,Point,Dim,Dim,Dim> shards_array_cfpddd(&data[0],2,2,2,2,2,2);
00262       FieldContainer<double> fc_shards_c_deep(shards_array_c, true);
00263       FieldContainer<double> fc_shards_cf_deep(shards_array_cf, true);
00264       FieldContainer<double> fc_shards_cfp_deep(shards_array_cfp, true);
00265       FieldContainer<double> fc_shards_cfpd_deep(shards_array_cfpd, true);
00266       FieldContainer<double> fc_shards_cfpdd_deep(shards_array_cfpdd, true);
00267       FieldContainer<double> fc_shards_cfpddd_deep(shards_array_cfpddd, true);
00268       // Create FieldContainers using a shallow copy via shards::Array.
00269       FieldContainer<double> fc_shards_c_shallow(shards_array_c);
00270       FieldContainer<double> fc_shards_cf_shallow(shards_array_cf);
00271       FieldContainer<double> fc_shards_cfp_shallow(shards_array_cfp);
00272       FieldContainer<double> fc_shards_cfpd_shallow(shards_array_cfpd);
00273       FieldContainer<double> fc_shards_cfpdd_shallow(shards_array_cfpdd);
00274       FieldContainer<double> fc_shards_cfpddd_shallow(shards_array_cfpddd);
00275     } // end variable scope
00276 
00277 
00278     *outStream << "\n" \
00279       << "===============================================================================\n"\
00280       << "| TEST 2: Usage of resize                                                     |\n"\
00281       << "===============================================================================\n\n";
00282 
00283     { // start variable scope
00284       // Initialize dimensions for rank-4 multi-index value
00285       dimensions.resize(5);
00286       dimensions[0] = 5;
00287       dimensions[1] = 3;
00288       dimensions[2] = 2;
00289       dimensions[3] = 7;
00290       dimensions[4] = 2;
00291 
00292       // temporary ints
00293       int d0, d1, d2, d3;
00294 
00295       // Allocate data through a Teuchos::Array, initialized to 0.
00296       Teuchos::Array<double> data(dimensions[0]*dimensions[1]*dimensions[2]*dimensions[3]*dimensions[4]);
00297 
00298       // Create a FieldContainer using a deep copy via Teuchos::Array.
00299       FieldContainer<double> fc_array(dimensions, data);
00300       // modify the (1,1,1,1) entry
00301       double mod_entry    = 1.0;
00302       fc_array(1,1,1,1,1) = mod_entry;
00303       int enumeration = fc_array.getEnumeration(1,1,1,1,1);
00304 
00305       // Resize into a (dimensions[0]*dimensions[1]) x (dimensions[2]) x (dimensions[3]) x (dimensions[4]) rank-4 array.
00306       // This is effectively a reshape, the data should not be touched.
00307       fc_array.resize(dimensions[0]*dimensions[1], dimensions[2], dimensions[3], dimensions[4]);
00308       // verify that the data array has NOT changed
00309       fc_array.getMultiIndex(d0,d1,d2,d3, enumeration);
00310       if (std::abs(fc_array(d0,d1,d2,d3) - mod_entry) > zero) {
00311         *outStream << "\n\nError in resize.\n\n";
00312         errorFlag = -1000;
00313       }
00314 
00315       // Resize into a (dimensions[0]*dimensions[1]*dimensions[2]) x (dimensions[3]) x (dimensions[4]) rank-3 array.
00316       // This is effectively a reshape, the data should not be touched.
00317       fc_array.resize(dimensions[0]*dimensions[1]*dimensions[2], dimensions[3], dimensions[4]);
00318       // verify that the data array has NOT changed
00319       fc_array.getMultiIndex(d0,d1,d2, enumeration);
00320       if (std::abs(fc_array(d0,d1,d2) - mod_entry) > zero) {
00321         *outStream << "\n\nError in resize.\n\n";
00322         errorFlag = -1000;
00323       }
00324 
00325       // Resize into a (dimensions[0]*dimensions[1]*dimensions[2]*dimensions[3]) x (dimensions[4]) rank-2 array.
00326       // This is effectively a reshape, the data should not be touched.
00327       fc_array.resize(dimensions[0]*dimensions[1]*dimensions[2]*dimensions[3], dimensions[4]);
00328       // verify that the data array has NOT changed
00329       fc_array.getMultiIndex(d0,d1, enumeration);
00330       if (std::abs(fc_array(d0,d1) - mod_entry) > zero) {
00331         *outStream << "\n\nError in resize.\n\n";
00332         errorFlag = -1000;
00333       }
00334 
00335       // Resize into a (dimensions[0]*dimensions[1]*dimensions[2]*dimensions[3]*dimensions[4]) rank-1 array.
00336       // This is effectively a reshape, the data should not be touched.
00337       fc_array.resize(dimensions[0]*dimensions[1]*dimensions[2]*dimensions[3]*dimensions[4]);
00338       // verify that the data array has NOT changed
00339       fc_array.getMultiIndex(d0, enumeration);
00340       if (std::abs(fc_array(d0) - mod_entry) > zero) {
00341         *outStream << "\n\nError in resize.\n\n";
00342         errorFlag = -1000;
00343       }
00344 
00345       // More advanced example allocating new memory, with shards::Array.
00346       data.assign(dimensions[0]*dimensions[1]*dimensions[2]*dimensions[3]*dimensions[4], 3.0);
00347       shards::Array<double,shards::NaturalOrder,Cell,Field,Point,Dim,Dim>
00348         shards_array(data.getRawPtr(),dimensions[0],dimensions[1],dimensions[2],dimensions[3],dimensions[4]);
00349       // Create a FieldContainer using a shallow copy via shards::Array.
00350       FieldContainer<double> fc_shards_shallow(shards_array);
00351       fc_shards_shallow.resize(4,4,4,4,4);  // keep old entries + allocate new memory and fill with zeros
00352       fc_shards_shallow.resize(4*4*4*4*4);  // reshape into rank-1 vector
00353       if (std::abs(RealSpaceTools<double>::vectorNorm(data.getRawPtr(), fc_array.size(), NORM_ONE) -
00354                    RealSpaceTools<double>::vectorNorm(fc_shards_shallow, NORM_ONE)) > zero) {
00355         *outStream << "\n\nError in resize.\n\n";
00356         errorFlag = -1000;
00357       }
00358 
00359 
00360     } // end variable scope
00361 
00362 
00363   } // outer try block
00364   catch (std::logic_error err) {
00365     *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
00366     *outStream  << err.what() << "\n";
00367     *outStream << "-------------------------------------------------------------------------------" << "\n\n";
00368     errorFlag = -1000;
00369   }
00370   
00371   if (errorFlag != 0)
00372     std::cout << "End Result: TEST FAILED\n";
00373   else
00374     std::cout << "End Result: TEST PASSED\n";
00375   
00376   // reset format state of std::cout
00377   std::cout.copyfmt(oldFormatState);
00378   
00379   return errorFlag;
00380 }