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