Intrepid
http://trilinos.sandia.gov/packages/docs/r10.10/packages/intrepid/example/FieldContainer/example_01.cpp
Go to the documentation of this file.
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_FieldContainer.hpp"
00050 #include "Teuchos_Time.hpp"
00051 #include "Teuchos_GlobalMPISession.hpp"
00052 
00053 using namespace std;
00054 using namespace Intrepid;
00055 
00056 int main(int argc, char *argv[]) {
00057 
00058   Teuchos::GlobalMPISession mpiSession(&argc, &argv);
00059 
00060   std::cout \
00061   << "===============================================================================\n" \
00062   << "|                                                                             |\n" \
00063   << "|                  Example use of the FieldContainer class                    |\n" \
00064   << "|                                                                             |\n" \
00065   << "|    1) Creating and filling FieldContainer objects                           |\n" \
00066   << "|    2) Accessing elements in FieldContainer objects                          |\n" \
00067   << "|                                                                             |\n" \
00068   << "|  Questions? Contact  Pavel Bochev (pbboche@sandia.gov) or                   |\n" \
00069   << "|                      Denis Ridzal (dridzal@sandia.gov).                     |\n" \
00070   << "|                                                                             |\n" \
00071   << "|  Intrepid's website: http://trilinos.sandia.gov/packages/intrepid           |\n" \
00072   << "|  Trilinos website:   http://trilinos.sandia.gov                             |\n" \
00073   << "|                                                                             |\n" \
00074   << "===============================================================================\n\n";
00075   
00076   // Define variables to create and use FieldContainers
00077   Teuchos::Array<int> dimension;
00078   Teuchos::Array<int> multiIndex;
00079   
00080   std::cout \
00081     << "===============================================================================\n"\
00082     << "| EXAMPLE 1: rank 2 multi-index: {u(p,i) | 0 <= p < 5; 0 <= i < 3 }           |\n"\
00083     << "===============================================================================\n\n";
00084   // This rank 2 multi-indexed value can be used to store values of 3D vector field
00085   // evaluated at 5 points, or values of gradient of scalar field evaluated at the points
00086   
00087   // Resize dimension and multiIndex for rank-2 multi-indexed value 
00088   dimension.resize(2);
00089   multiIndex.resize(2);
00090   
00091   // Load upper ranges for the two indices in the multi-indexed value
00092   dimension[0] = 5;
00093   dimension[1] = 3;
00094     
00095   // Create FieldContainer that can hold the rank-2 multi-indexed value
00096   FieldContainer<double> myContainer(dimension);
00097   
00098   // Fill with some data: leftmost index changes last, rightmost index changes first!  
00099   for(int p = 0; p < dimension[0]; p++){
00100     multiIndex[0] = p;
00101     
00102     for(int i = 0; i < dimension[1]; i++){
00103       multiIndex[1] = i;
00104       
00105       // Load value with multi-index {p,i} to container
00106       myContainer.setValue((double)(i+p), multiIndex);
00107     }
00108   }
00109   
00110   // Show container contents
00111   std::cout << myContainer;
00112   
00113   // Access by overloaded (), multiindex and []:
00114   multiIndex[0] = 3; 
00115   multiIndex[1] = 1;
00116   int enumeration = myContainer.getEnumeration(multiIndex);
00117 
00118   std::cout << "Access by ():          myContainer(" << 3 <<"," << 1 << ") = " << myContainer(3,1) << "\n";  
00119   std::cout << "Access by multi-index: myContainer{" << multiIndex[0] << multiIndex[1] << "} = "<< myContainer.getValue(multiIndex) <<"\n";
00120   std::cout << "Access by enumeration: myContainer[" << enumeration << "] = " << myContainer[enumeration] <<"\n";
00121   
00122   std::cout << "Assigning value by (): \n old value at (3,1) = " << myContainer(3,1) <<"\n";
00123   myContainer(3,1) = 999.99;
00124   std::cout << " new value at (3,1) = " << myContainer(3,1) <<"\n";
00125       
00126   
00127   std::cout << "\n" \
00128   << "===============================================================================\n"\
00129   << "| EXAMPLE 2: rank 3 multi-index: {u(p,i,j) | 0 <=p< 5; 0 <= i<2, 0<=j<3       |\n"\
00130   << "===============================================================================\n\n";
00131   // This rank-3 value can be used to store subset of second partial derivatives  values
00132   // of a scalar function at p points.
00133   
00134   // Resize dimension and multiIndex for rank-3 multi-indexed value 
00135   dimension.resize(3);
00136   multiIndex.resize(3);
00137   
00138   // Define upper ranges for the three indices in the multi-indexed value
00139   dimension[0] = 5;
00140   dimension[1] = 2;
00141   dimension[2] = 3;
00142   
00143   // Reset the existing container to accept rank-3 value with the specified index ranges
00144   myContainer.resize(dimension);
00145   
00146   // Fill with some data
00147   for(int p = 0; p < dimension[0]; p++){
00148     multiIndex[0] = p;
00149     for(int i = 0; i < dimension[1]; i++){
00150       multiIndex[1] = i;
00151       for(int j = 0; j < dimension[2]; j++){
00152         multiIndex[2] = j;
00153         
00154         // Load value with multi-index {p,i} to container
00155         myContainer.setValue((double)(p+i+j), multiIndex);
00156       }
00157     }
00158   }
00159   
00160   // Display contents
00161   std::cout << myContainer;
00162   
00163   // Access by overloaded (), multiindex and []:
00164   multiIndex[0] = 3; 
00165   multiIndex[1] = 1;
00166   multiIndex[2] = 2;
00167   enumeration = myContainer.getEnumeration(multiIndex);
00168     
00169   std::cout << "Access by ():          myContainer(" << 3 <<"," << 1 << "," << 2 << ") = " << myContainer(3,1,2) << "\n";  
00170   std::cout << "Access by multi-index: myContainer{" << multiIndex[0] << multiIndex[1] << multiIndex[2] << "} = "<< myContainer.getValue(multiIndex) <<"\n";
00171   std::cout << "Access by enumeration: myContainer[" << enumeration << "] = " << myContainer[enumeration] <<"\n";
00172 
00173   std::cout << "Assigning value by (): \n old value at (3,1,2) = " << myContainer(3,1,2) <<"\n";
00174   myContainer(3,1,2) = -999.999;
00175   std::cout << " new value at (3,1,2) = " << myContainer(3,1,2) <<"\n";
00176   
00177   
00178   std::cout << "\n" \
00179   << "===============================================================================\n"\
00180   << "| EXAMPLE 4: making rank-5 FieldContainer from data array and index range array |\n"\
00181   << "===============================================================================\n\n";
00182   // Initialize dimension for rank-5 multi-index value
00183   dimension.resize(5);
00184   dimension[0] = 5;
00185   dimension[1] = 2;
00186   dimension[2] = 3;
00187   dimension[3] = 4;
00188   dimension[4] = 6;
00189   
00190   // Define Teuchos::Array to store values with dimension equal to the number of multi-indexed values
00191   Teuchos::Array<double> dataTeuchosArray(5*2*3*4*6);
00192   
00193   // Fill with data
00194   int counter = 0;
00195   for(int i=0; i < dimension[0]; i++){
00196     for(int j=0; j < dimension[1]; j++){
00197       for(int k=0; k < dimension[2]; k++){
00198         for(int l = 0; l < dimension[3]; l++){
00199           for(int m = 0; m < dimension[4]; m++){
00200             dataTeuchosArray[counter] = (double)counter;
00201             counter++;
00202           }
00203         }
00204       }
00205     }
00206   }
00207   
00208   // Create FieldContainer from data array and index array and show it
00209   FieldContainer<double> myNewContainer(dimension, dataTeuchosArray);
00210   std::cout << myNewContainer;
00211 
00212   // Access by overloaded (), multiindex and []:
00213   multiIndex.resize(myNewContainer.rank());
00214   multiIndex[0] = 3; 
00215   multiIndex[1] = 1;
00216   multiIndex[2] = 2;
00217   multiIndex[3] = 2;
00218   multiIndex[4] = 5;
00219   enumeration = myNewContainer.getEnumeration(multiIndex);
00220     
00221   std::cout << "Access by ():          myNewContainer(" << 3 <<"," << 1 << "," << 2 << "," << 2 << "," << 5 << ") = " << myNewContainer(3,1,2,2,5) << "\n";  
00222   std::cout << "Access by multi-index: myNewContainer{" << multiIndex[0] << multiIndex[1] << multiIndex[2] << multiIndex[3] << multiIndex[4] << "} = "<< myNewContainer.getValue(multiIndex) <<"\n";
00223   std::cout << "Access by enumeration: myNewContainer[" << enumeration << "] = " << myNewContainer[enumeration] <<"\n";
00224 
00225   std::cout << "Assigning value by (): \n old value at (3,1,2,2,5) = " << myNewContainer(3,1,2,2,5) <<"\n";
00226   myNewContainer(3,1,2,2,5) = -888.888;
00227   std::cout << " new value at (3,1,2,2,5) = " << myNewContainer(3,1,2,2,5) <<"\n";
00228   
00229   
00230   std::cout << "\n" \
00231     << "===============================================================================\n"\
00232     << "| EXAMPLE 5: making trivial FieldContainers and storing a single zero         |\n"\
00233     << "===============================================================================\n\n";
00234   
00235   // Make trivial container by resetting the index range to zero rank (no indices) and then
00236   // using resize method
00237   dimension.resize(0);
00238   myContainer.resize(dimension);
00239   std::cout << myContainer;
00240   
00241   // Make trivial container by using clear method:
00242   myNewContainer.clear();
00243   std::cout << myNewContainer;
00244   
00245   // Now use initialize() to reset the container to hold a single zero
00246   myNewContainer.initialize();
00247   std::cout << myNewContainer;
00248   
00249   
00250   std::cout << "\n" \
00251     << "===============================================================================\n"\
00252     << "| EXAMPLE 6: Timing read and write operations using () and getValue           |\n"\
00253     << "===============================================================================\n\n";
00254   
00255   // Initialize dimensions for rank-5 multi-index value
00256   int dim0 = 10;     // number of cells
00257   int dim1 = 50;     // number of points
00258   int dim2 = 27;     // number of functions
00259   int dim3 = 3;      // 1st space dim
00260   int dim4 = 3;      // 2nd space dim
00261   
00262   FieldContainer<double> myTensorContainer(dim0, dim1, dim2, dim3, dim4);
00263   multiIndex.resize(myTensorContainer.rank());
00264   double aValue;
00265   
00266   Teuchos::Time timerGetValue("Reading and writing from rank-5 container using getValue");
00267   timerGetValue.start();
00268   for(int i0 = 0; i0 < dim0; i0++){
00269     multiIndex[0] = i0;
00270     for(int i1 = 0; i1 < dim1; i1++){ 
00271       multiIndex[1] = i1;
00272       for(int i2 = 0; i2 < dim2; i2++) {
00273         multiIndex[2] = i2;
00274         for(int i3 = 0; i3 < dim3; i3++) {
00275           multiIndex[3] = i3;
00276           for(int i4 =0; i4 < dim4; i4++) { 
00277             multiIndex[4] = i4;
00278             
00279             aValue = myTensorContainer.getValue(multiIndex);
00280             myTensorContainer.setValue(999.999,multiIndex);
00281             
00282           }
00283         }
00284       }
00285     }
00286   }
00287   timerGetValue.stop();
00288   std::cout << " Time to read and write from container using getValue: " << timerGetValue.totalElapsedTime() <<"\n";
00289   
00290   Teuchos::Time timerRound("Reading and writing from rank-5 container using ()");
00291   timerRound.start();
00292   for(int i0 = 0; i0 < dim0; i0++){
00293     for(int i1 = 0; i1 < dim1; i1++) { 
00294       for(int i2 = 0; i2 < dim2; i2++) {
00295         for(int i3 = 0; i3 < dim3; i3++) {
00296           for(int i4 =0; i4 < dim4; i4++) { 
00297             
00298             aValue = myTensorContainer(i0,i1,i2,i3,i4);
00299             myTensorContainer(i0,i1,i2,i3,i4) = 999.999;
00300             
00301           }
00302         }
00303       }
00304     }
00305   }
00306   timerRound.stop();
00307   std::cout << " Time to read and write from container using (): " << timerRound.totalElapsedTime() <<"\n";
00308     
00309   std::cout << "\n" \
00310     << "===============================================================================\n"\
00311     << "| EXAMPLE 6: Specialized methods of FieldContainer                            |\n"\
00312     << "===============================================================================\n\n";
00313   
00314    return 0;
00315 }
00316 
00317 
00318 
00319 
00320