Intrepid
http://trilinos.sandia.gov/packages/docs/dev/packages/intrepid/src/Shared/Intrepid_FieldContainerDef.hpp
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 namespace Intrepid {
00050 
00051   
00052   //--------------------------------------------------------------------------------------------//
00053   //                                                                                            //
00054   //                 Member function definitions of the class FieldContainer                    //
00055   //                                                                                            //
00056   //--------------------------------------------------------------------------------------------//
00057   
00058 
00059 template<class Scalar, int ArrayTypeId>
00060 FieldContainer<Scalar, ArrayTypeId>::FieldContainer(const FieldContainer<Scalar, ArrayTypeId>& right) {
00061   
00062   // Copy dimensions and data values from right
00063   dimensions_.assign(right.dimensions_.begin(),right.dimensions_.end());  
00064   data_.assign(right.data_.begin(),right.data_.end());
00065   dim0_ = right.dim0_;
00066   dim1_ = right.dim1_;
00067   dim2_ = right.dim2_;
00068   dim3_ = right.dim3_;
00069   dim4_ = right.dim4_;
00070   data_ptr_ = data_.begin();
00071 }
00072 
00073 //--------------------------------------------------------------------------------------------//
00074 //                                                                                            //
00075 //                              Constructors of FieldContainer class                          //
00076 //                                                                                            //
00077 //--------------------------------------------------------------------------------------------//
00078 
00079 
00080 template<class Scalar, int ArrayTypeId>
00081 FieldContainer<Scalar, ArrayTypeId>::FieldContainer(const int dim0) : dim0_(dim0), dim1_(0), dim2_(0), dim3_(0), dim4_(0) 
00082 {
00083   using Teuchos::as;
00084   using Teuchos::Ordinal;
00085 #ifdef HAVE_INTREPID_DEBUG
00086   TEUCHOS_TEST_FOR_EXCEPTION( (0 > dim0), std::invalid_argument, 
00087                       ">>> ERROR (FieldContainer): FieldContainer cannot have a negative dimension.");
00088 
00089 #endif
00090   dimensions_.resize(as<Ordinal>(1)); 
00091   dimensions_[0] = dim0_;
00092   data_.assign(as<Ordinal>(dim0_), as<Scalar>(0));
00093   data_ptr_ = data_.begin();
00094 }
00095 
00096 
00097 
00098 template<class Scalar, int ArrayTypeId>
00099 FieldContainer<Scalar, ArrayTypeId>::FieldContainer(const int dim0,
00100                                        const int dim1) : dim0_(dim0), dim1_(dim1), dim2_(0), dim3_(0), dim4_(0)
00101 {
00102   using Teuchos::as;
00103   using Teuchos::Ordinal;
00104 #ifdef HAVE_INTREPID_DEBUG
00105   TEUCHOS_TEST_FOR_EXCEPTION( (0 > dim0), std::invalid_argument, 
00106                       ">>> ERROR (FieldContainer): FieldContainer cannot have a negative 1st dimension.");
00107   TEUCHOS_TEST_FOR_EXCEPTION( (0 > dim1), std::invalid_argument, 
00108                       ">>> ERROR (FieldContainer): FieldContainer cannot have a negative 2nd dimension.");
00109   
00110 #endif
00111   dimensions_.resize(2); 
00112   dimensions_[0] = dim0_;   
00113   dimensions_[1] = dim1_;  
00114   data_.assign(as<Ordinal>(dim0_*dim1_), as<Scalar>(0));
00115   data_ptr_ = data_.begin();
00116 }
00117 
00118 
00119 
00120 template<class Scalar, int ArrayTypeId>
00121 FieldContainer<Scalar, ArrayTypeId>::FieldContainer(const int dim0,
00122                                        const int dim1,
00123                                        const int dim2) : dim0_(dim0), dim1_(dim1), dim2_(dim2), dim3_(0), dim4_(0)
00124 {
00125   using Teuchos::as;
00126   using Teuchos::Ordinal;
00127 #ifdef HAVE_INTREPID_DEBUG
00128   TEUCHOS_TEST_FOR_EXCEPTION( (0 > dim0), std::invalid_argument, 
00129                       ">>> ERROR (FieldContainer): FieldContainer cannot have a negative 1st dimension.");
00130   TEUCHOS_TEST_FOR_EXCEPTION( (0 > dim1), std::invalid_argument, 
00131                       ">>> ERROR (FieldContainer): FieldContainer cannot have a negative 2nd dimension.");
00132   TEUCHOS_TEST_FOR_EXCEPTION( (0 > dim2), std::invalid_argument, 
00133                       ">>> ERROR (FieldContainer): FieldContainer cannot have a negative 3rd dimension.");
00134 #endif
00135   dimensions_.resize(3); 
00136   dimensions_[0] = dim0_;  
00137   dimensions_[1] = dim1_; 
00138   dimensions_[2] = dim2_;  
00139   data_.assign(as<Ordinal>(dim0_*dim1_*dim2_), as<Scalar>(0));
00140   data_ptr_ = data_.begin();
00141 }
00142 
00143 
00144 
00145 template<class Scalar, int ArrayTypeId>
00146 FieldContainer<Scalar, ArrayTypeId>::FieldContainer(const int dim0,
00147                                        const int dim1,
00148                                        const int dim2,
00149                                        const int dim3) : dim0_(dim0), dim1_(dim1), dim2_(dim2), dim3_(dim3), dim4_(0)
00150 {
00151   using Teuchos::as;
00152   using Teuchos::Ordinal;
00153 #ifdef HAVE_INTREPID_DEBUG
00154   TEUCHOS_TEST_FOR_EXCEPTION( (0 > dim0), std::invalid_argument, 
00155                       ">>> ERROR (FieldContainer): FieldContainer cannot have a negative 1st dimension.");
00156   TEUCHOS_TEST_FOR_EXCEPTION( (0 > dim1), std::invalid_argument, 
00157                       ">>> ERROR (FieldContainer): FieldContainer cannot have a negative 2nd dimension.");
00158   TEUCHOS_TEST_FOR_EXCEPTION( (0 > dim2), std::invalid_argument, 
00159                       ">>> ERROR (FieldContainer): FieldContainer cannot have a negative 3rd dimension.");
00160   TEUCHOS_TEST_FOR_EXCEPTION( (0 > dim3), std::invalid_argument, 
00161                       ">>> ERROR (FieldContainer): FieldContainer cannot have a negative 4th dimension.");  
00162 #endif
00163   dimensions_.resize(4); 
00164   dimensions_[0] = dim0_;  
00165   dimensions_[1] = dim1_; 
00166   dimensions_[2] = dim2_;  
00167   dimensions_[3] = dim3_; 
00168   data_.assign(as<Ordinal>(dim0_*dim1_*dim2_*dim3_), as<Scalar>(0));
00169   data_ptr_ = data_.begin();
00170 }
00171 
00172 
00173 
00174 template<class Scalar, int ArrayTypeId>
00175 FieldContainer<Scalar, ArrayTypeId>::FieldContainer(const int dim0,
00176                                        const int dim1,
00177                                        const int dim2,
00178                                        const int dim3,
00179                                        const int dim4) : dim0_(dim0), dim1_(dim1), dim2_(dim2), dim3_(dim3), dim4_(dim4)
00180 {
00181   using Teuchos::as;
00182   using Teuchos::Ordinal;
00183 #ifdef HAVE_INTREPID_DEBUG
00184   TEUCHOS_TEST_FOR_EXCEPTION( (0 > dim0), std::invalid_argument, 
00185                       ">>> ERROR (FieldContainer): FieldContainer cannot have a negative 1st dimension.");
00186   TEUCHOS_TEST_FOR_EXCEPTION( (0 > dim1), std::invalid_argument, 
00187                       ">>> ERROR (FieldContainer): FieldContainer cannot have a negative 2nd dimension.");
00188   TEUCHOS_TEST_FOR_EXCEPTION( (0 > dim2), std::invalid_argument, 
00189                       ">>> ERROR (FieldContainer): FieldContainer cannot have a negative 3rd dimension.");
00190   TEUCHOS_TEST_FOR_EXCEPTION( (0 > dim3), std::invalid_argument, 
00191                       ">>> ERROR (FieldContainer): FieldContainer cannot have a negative 4th dimension.");  
00192   TEUCHOS_TEST_FOR_EXCEPTION( (0 > dim4), std::invalid_argument, 
00193                       ">>> ERROR (FieldContainer): FieldContainer cannot have a negative 5th dimension.");  
00194 #endif
00195   dimensions_.resize(5); 
00196   dimensions_[0] = dim0_;
00197   dimensions_[1] = dim1_;
00198   dimensions_[2] = dim2_;  
00199   dimensions_[3] = dim3_;  
00200   dimensions_[4] = dim4_;  
00201   data_.assign(as<Ordinal>(dim0_*dim1_*dim2_*dim3_*dim4_), as<Scalar>(0));
00202   data_ptr_ = data_.begin();
00203 }
00204 
00205 
00206 
00207 template<class Scalar, int ArrayTypeId>
00208 FieldContainer<Scalar, ArrayTypeId>::FieldContainer(const Teuchos::Array<int>& dimensions) {  
00209   
00210 #ifdef HAVE_INTREPID_DEBUG
00211 // srkenno@sandia.gov 6/12/10: changed unsigned int to int - this was causing a warning on compilers that 
00212 //   signed & unsigned int's were being comparied.
00213   for( int dim = 0; dim < dimensions.size(); dim++) {
00214     TEUCHOS_TEST_FOR_EXCEPTION( (0 > dimensions[dim] ), std::invalid_argument,  
00215                         ">>> ERROR (FieldContainer): One or more negative dimensions");  
00216   }
00217 #endif
00218   
00219   // Copy dimensions and resize container storage to match them
00220   dimensions_.assign(dimensions.begin(),dimensions.end());  
00221   
00222   // Copy first 5 dimensions to optimize class for low rank containers
00223   unsigned int rank = dimensions_.size();
00224   switch(rank) {
00225     case 1:
00226       dim0_ = dimensions_[0]; 
00227       dim1_ = 0;
00228       dim2_ = 0;
00229       dim3_ = 0;
00230       dim4_ = 0;
00231       break;
00232       
00233     case 2:
00234       dim0_ = dimensions_[0]; 
00235       dim1_ = dimensions_[1]; 
00236       dim2_ = 0;
00237       dim3_ = 0;
00238       dim4_ = 0;
00239       break;
00240       
00241     case 3:
00242       dim0_ = dimensions_[0]; 
00243       dim1_ = dimensions_[1]; 
00244       dim2_ = dimensions_[2]; 
00245       dim3_ = 0;
00246       dim4_ = 0;
00247       break;
00248       
00249     case 4:
00250       dim0_ = dimensions_[0]; 
00251       dim1_ = dimensions_[1]; 
00252       dim2_ = dimensions_[2]; 
00253       dim3_ = dimensions_[3]; 
00254       dim4_ = 0;
00255       break;
00256       
00257     case 5:
00258     default:
00259       dim0_ = dimensions_[0]; 
00260       dim1_ = dimensions_[1]; 
00261       dim2_ = dimensions_[2]; 
00262       dim3_ = dimensions_[3]; 
00263       dim4_ = dimensions_[4]; 
00264   }
00265 
00266   // resize data array according to specified dimensions
00267   data_.assign( this -> size(), (Scalar)0);
00268   data_ptr_ = data_.begin();
00269   
00270 }
00271 
00272 
00273 
00274 template<class Scalar, int ArrayTypeId>
00275 FieldContainer<Scalar, ArrayTypeId>::FieldContainer(const Teuchos::Array<int>&         dimensions,
00276                                        const Teuchos::ArrayView<Scalar>&  data) {
00277  
00278   // Copy all dimensions
00279   dimensions_.assign(dimensions.begin(),dimensions.end());
00280   
00281   // Copy first 5 dimensions to optimize class for low rank containers
00282   unsigned int rank = dimensions_.size();
00283   switch(rank) {
00284     case 0:
00285       dim0_ = 0; 
00286       dim1_ = 0;
00287       dim2_ = 0;
00288       dim3_ = 0;
00289       dim4_ = 0;
00290       break;
00291 
00292     case 1:
00293       dim0_ = dimensions_[0]; 
00294       dim1_ = 0;
00295       dim2_ = 0;
00296       dim3_ = 0;
00297       dim4_ = 0;
00298       break;
00299       
00300     case 2:
00301       dim0_ = dimensions_[0]; 
00302       dim1_ = dimensions_[1]; 
00303       dim2_ = 0;
00304       dim3_ = 0;
00305       dim4_ = 0;
00306       break;
00307       
00308     case 3:
00309       dim0_ = dimensions_[0]; 
00310       dim1_ = dimensions_[1]; 
00311       dim2_ = dimensions_[2]; 
00312       dim3_ = 0;
00313       dim4_ = 0;
00314       break;
00315       
00316     case 4:
00317       dim0_ = dimensions_[0]; 
00318       dim1_ = dimensions_[1]; 
00319       dim2_ = dimensions_[2]; 
00320       dim3_ = dimensions_[3]; 
00321       dim4_ = 0;
00322       break;
00323       
00324     case 5:
00325     default:
00326       dim0_ = dimensions_[0]; 
00327       dim1_ = dimensions_[1]; 
00328       dim2_ = dimensions_[2]; 
00329       dim3_ = dimensions_[3]; 
00330       dim4_ = dimensions_[4]; 
00331   }
00332   
00333     // Validate input: size of data array must match container size specified by its dimensions
00334 #ifdef HAVE_INTREPID_DEBUG
00335   TEUCHOS_TEST_FOR_EXCEPTION( ( (int)data.size() != this -> size() ),
00336                       std::invalid_argument,
00337                       ">>> ERROR (FieldContainer): Size of input data does not match size of this container.");
00338 #endif
00339   
00340   // Deep-copy ArrayView data.
00341   data_.assign(data.begin(),data.end());
00342   data_ptr_ = data_.begin();
00343 }
00344 
00345 
00346 
00347 template<class Scalar, int ArrayTypeId>
00348 FieldContainer<Scalar, ArrayTypeId>::FieldContainer(const Teuchos::Array<int>&        dimensions,
00349                                        const Teuchos::ArrayRCP<Scalar>&  data) {
00350  
00351   // Copy all dimensions
00352   dimensions_.assign(dimensions.begin(),dimensions.end());
00353   
00354   // Copy first 5 dimensions to optimize class for low rank containers
00355   unsigned int rank = dimensions_.size();
00356   switch(rank) {
00357     case 0:
00358       dim0_ = 0; 
00359       dim1_ = 0;
00360       dim2_ = 0;
00361       dim3_ = 0;
00362       dim4_ = 0;
00363       break;
00364 
00365     case 1:
00366       dim0_ = dimensions_[0]; 
00367       dim1_ = 0;
00368       dim2_ = 0;
00369       dim3_ = 0;
00370       dim4_ = 0;
00371       break;
00372       
00373     case 2:
00374       dim0_ = dimensions_[0]; 
00375       dim1_ = dimensions_[1]; 
00376       dim2_ = 0;
00377       dim3_ = 0;
00378       dim4_ = 0;
00379       break;
00380       
00381     case 3:
00382       dim0_ = dimensions_[0]; 
00383       dim1_ = dimensions_[1]; 
00384       dim2_ = dimensions_[2]; 
00385       dim3_ = 0;
00386       dim4_ = 0;
00387       break;
00388       
00389     case 4:
00390       dim0_ = dimensions_[0]; 
00391       dim1_ = dimensions_[1]; 
00392       dim2_ = dimensions_[2]; 
00393       dim3_ = dimensions_[3]; 
00394       dim4_ = 0;
00395       break;
00396       
00397     case 5:
00398     default:
00399       dim0_ = dimensions_[0]; 
00400       dim1_ = dimensions_[1]; 
00401       dim2_ = dimensions_[2]; 
00402       dim3_ = dimensions_[3]; 
00403       dim4_ = dimensions_[4]; 
00404   }
00405   
00406     // Validate input: size of data array must match container size specified by its dimensions
00407 #ifdef HAVE_INTREPID_DEBUG
00408   TEUCHOS_TEST_FOR_EXCEPTION( ( (int)data.size() != this -> size() ),
00409                       std::invalid_argument,
00410                       ">>> ERROR (FieldContainer): Size of input data does not match size of this container.");
00411 #endif
00412   
00413   // Shallow-copy ArrayRCP data.
00414   data_ = data;
00415   data_ptr_ = data_.begin();
00416 }
00417 
00418 
00419 
00420 template<class Scalar, int ArrayTypeId>
00421 FieldContainer<Scalar, ArrayTypeId>::FieldContainer(const Teuchos::Array<int>&    dimensions,
00422                                        Scalar*                       data,
00423                                        const bool                    deep_copy,
00424                                        const bool                    owns_mem) {
00425  
00426   // Copy all dimensions
00427   dimensions_.assign(dimensions.begin(),dimensions.end());
00428   
00429   // Copy first 5 dimensions to optimize class for low rank containers
00430   unsigned int rank = dimensions_.size();
00431   switch(rank) {
00432     case 0:
00433       dim0_ = 0; 
00434       dim1_ = 0;
00435       dim2_ = 0;
00436       dim3_ = 0;
00437       dim4_ = 0;
00438       break;
00439 
00440     case 1:
00441       dim0_ = dimensions_[0]; 
00442       dim1_ = 0;
00443       dim2_ = 0;
00444       dim3_ = 0;
00445       dim4_ = 0;
00446       break;
00447       
00448     case 2:
00449       dim0_ = dimensions_[0]; 
00450       dim1_ = dimensions_[1]; 
00451       dim2_ = 0;
00452       dim3_ = 0;
00453       dim4_ = 0;
00454       break;
00455       
00456     case 3:
00457       dim0_ = dimensions_[0]; 
00458       dim1_ = dimensions_[1]; 
00459       dim2_ = dimensions_[2]; 
00460       dim3_ = 0;
00461       dim4_ = 0;
00462       break;
00463       
00464     case 4:
00465       dim0_ = dimensions_[0]; 
00466       dim1_ = dimensions_[1]; 
00467       dim2_ = dimensions_[2]; 
00468       dim3_ = dimensions_[3]; 
00469       dim4_ = 0;
00470       break;
00471       
00472     case 5:
00473     default:
00474       dim0_ = dimensions_[0]; 
00475       dim1_ = dimensions_[1]; 
00476       dim2_ = dimensions_[2]; 
00477       dim3_ = dimensions_[3]; 
00478       dim4_ = dimensions_[4]; 
00479   }
00480   
00481 
00482   if (deep_copy) {
00483     Teuchos::ArrayRCP<Scalar> arrayrcp = Teuchos::arcp<Scalar>(data, 0, this -> size(), false);
00484     data_.deepCopy(arrayrcp());
00485     data_ptr_ = data_.begin();
00486   }
00487   else {
00488     data_ = Teuchos::arcp<Scalar>(data, 0, this -> size(), owns_mem);
00489     data_ptr_ = data_.begin();
00490   }
00491 }
00492 
00493 
00494 
00495 template<class Scalar, int ArrayTypeId>
00496 FieldContainer<Scalar, ArrayTypeId>::FieldContainer(const shards::Array<Scalar,shards::NaturalOrder>&  data,
00497                                        const bool                                         deep_copy,
00498                                        const bool                                         owns_mem) {
00499  
00500   // Copy all dimensions
00501   dimensions_.resize(data.rank());
00502   
00503   // Copy first 5 dimensions to optimize class for low rank containers
00504   unsigned int rank = dimensions_.size();
00505   switch(rank) {
00506     case 1:
00507       dimensions_[0] = data.dimension(0);
00508       dim0_ = dimensions_[0]; 
00509       dim1_ = 0;
00510       dim2_ = 0;
00511       dim3_ = 0;
00512       dim4_ = 0;
00513       break;
00514       
00515     case 2:
00516       dimensions_[0] = data.dimension(0);
00517       dimensions_[1] = data.dimension(1);
00518       dim0_ = dimensions_[0]; 
00519       dim1_ = dimensions_[1]; 
00520       dim2_ = 0;
00521       dim3_ = 0;
00522       dim4_ = 0;
00523       break;
00524       
00525     case 3:
00526       dimensions_[0] = data.dimension(0);
00527       dimensions_[1] = data.dimension(1);
00528       dimensions_[2] = data.dimension(2);
00529       dim0_ = dimensions_[0]; 
00530       dim1_ = dimensions_[1]; 
00531       dim2_ = dimensions_[2]; 
00532       dim3_ = 0;
00533       dim4_ = 0;
00534       break;
00535       
00536     case 4:
00537       dimensions_[0] = data.dimension(0);
00538       dimensions_[1] = data.dimension(1);
00539       dimensions_[2] = data.dimension(2);
00540       dimensions_[3] = data.dimension(3);
00541       dim0_ = dimensions_[0]; 
00542       dim1_ = dimensions_[1]; 
00543       dim2_ = dimensions_[2]; 
00544       dim3_ = dimensions_[3]; 
00545       dim4_ = 0;
00546       break;
00547       
00548     case 5:
00549       dimensions_[0] = data.dimension(0);
00550       dimensions_[1] = data.dimension(1);
00551       dimensions_[2] = data.dimension(2);
00552       dimensions_[3] = data.dimension(3);
00553       dimensions_[4] = data.dimension(4);
00554       dim0_ = dimensions_[0]; 
00555       dim1_ = dimensions_[1]; 
00556       dim2_ = dimensions_[2]; 
00557       dim3_ = dimensions_[3]; 
00558       dim4_ = dimensions_[4]; 
00559       break;
00560 
00561     default:
00562       for (int i=0; i<data.rank(); i++) {
00563         dimensions_[i] = data.dimension(i);
00564       }
00565       dim0_ = dimensions_[0]; 
00566       dim1_ = dimensions_[1]; 
00567       dim2_ = dimensions_[2]; 
00568       dim3_ = dimensions_[3]; 
00569       dim4_ = dimensions_[4]; 
00570   }
00571   
00572 
00573   if (deep_copy) {
00574     Teuchos::ArrayRCP<Scalar> arrayrcp = Teuchos::arcp<Scalar>(data.contiguous_data(), 0, this -> size(), false);
00575     data_.deepCopy(arrayrcp());
00576     data_ptr_ = data_.begin();
00577   }
00578   else {
00579     data_ = Teuchos::arcp<Scalar>(data.contiguous_data(), 0, this -> size(), owns_mem);
00580     data_ptr_ = data_.begin();
00581   }
00582 }
00583 
00584 
00585 
00586 //--------------------------------------------------------------------------------------------//
00587 //                                                                                            //
00588 //                            Access methods of FieldContainer class                          //
00589 //                                                                                            //
00590 //--------------------------------------------------------------------------------------------//
00591 
00592 
00593 template<class Scalar, int ArrayTypeId>
00594 inline int FieldContainer<Scalar, ArrayTypeId>::rank() const {
00595   return dimensions_.size();
00596 }
00597   
00598 
00599 
00600 template<class Scalar, int ArrayTypeId>
00601 int FieldContainer<Scalar, ArrayTypeId>::size() const {
00602   // Important! This method is used by constructors to find out what is the needed size of data_
00603   // based on the specified dimensions. Therefore, it cannot be implmented by returning data_.size
00604   // and must be able to compute the size of the container based only on its specified dimensions
00605   
00606   // Size equals product of all dimensions stored in dimensions_
00607   int rank = dimensions_.size();
00608   
00609   // If container has no dimensions its size is zero
00610   if(rank == 0) {
00611     return 0;
00612   }
00613   else {
00614     int size = dim0_;
00615     
00616     // Compute size directly to optimize method for low rank (<=5) containers
00617     switch(rank) {
00618       case 5:
00619         size *= dim1_*dim2_*dim3_*dim4_;
00620         break;
00621         
00622       case 4:
00623         size *= dim1_*dim2_*dim3_;
00624         break;
00625         
00626       case 3:
00627         size *= dim1_*dim2_;
00628         break;
00629         
00630       case 2:
00631         size *= dim1_;
00632         break;
00633         
00634       case 1:
00635         break;
00636         
00637         // Compute size for containers with ranks hihger than 5
00638       default:
00639         for(int r = 1; r < rank ; r++){
00640           size *= dimensions_[r];
00641         }
00642     }
00643     return size;
00644   }
00645 }
00646 
00647 
00648 
00649 template<class Scalar, int ArrayTypeId>
00650 template<class Vector>
00651 inline void FieldContainer<Scalar, ArrayTypeId>::dimensions(Vector& dimensions) const {
00652   dimensions.assign(dimensions_.begin(),dimensions_.end());
00653 }
00654 
00655 
00656 
00657 template<class Scalar, int ArrayTypeId>
00658 inline int FieldContainer<Scalar, ArrayTypeId>::dimension(const int whichDim) const {
00659 #ifdef HAVE_INTREPID_DEBUG
00660   TEUCHOS_TEST_FOR_EXCEPTION( (0 > whichDim), std::invalid_argument,
00661                       ">>> ERROR (FieldContainer): dimension order cannot be negative");
00662   TEUCHOS_TEST_FOR_EXCEPTION( (whichDim >= this -> rank() ), std::invalid_argument,
00663                       ">>> ERROR (FieldContainer): dimension order cannot exceed rank of the container");
00664 #endif
00665   return dimensions_[whichDim];
00666 }
00667 
00668 
00669 
00670 template<class Scalar, int ArrayTypeId>
00671 inline int FieldContainer<Scalar, ArrayTypeId>::getEnumeration(const int i0) const {
00672 #ifdef HAVE_INTREPID_DEBUG
00673   TEUCHOS_TEST_FOR_EXCEPTION( ( this -> rank() != 1), std::invalid_argument, 
00674                       ">>> ERROR (FieldContainer): Number of indices does not match rank of the container.");  
00675   TEUCHOS_TEST_FOR_EXCEPTION( ( (i0 < 0) || ( i0 >= dim0_) ), std::invalid_argument,
00676                       ">>> ERROR (FieldContainer): index is out of range.");
00677 #endif
00678   return i0;
00679 }
00680 
00681 
00682 
00683 template<class Scalar, int ArrayTypeId>
00684 inline int FieldContainer<Scalar, ArrayTypeId>::getEnumeration(const int i0,
00685                                                   const int i1) const {
00686 #ifdef HAVE_INTREPID_DEBUG
00687   TEUCHOS_TEST_FOR_EXCEPTION( ( this -> rank() != 2), std::invalid_argument, 
00688                       ">>> ERROR (FieldContainer): Number of indices does not match rank of the container.");  
00689   TEUCHOS_TEST_FOR_EXCEPTION( ( (i0 < 0) || ( i0 >= dim0_) ), std::invalid_argument,
00690                       ">>> ERROR (FieldContainer): 1st index is out of range.");    
00691   TEUCHOS_TEST_FOR_EXCEPTION( ( (i1 < 0) || (i1 >= dim1_) ), std::invalid_argument,
00692                       ">>> ERROR (FieldContainer): 2nd index is out of range.");    
00693 #endif
00694   return i0*dim1_ + i1;
00695 }
00696 
00697 
00698 
00699 template<class Scalar, int ArrayTypeId>
00700 inline int FieldContainer<Scalar, ArrayTypeId>::getEnumeration(const int i0,
00701                                                   const int i1,
00702                                                   const int i2) const {
00703 #ifdef HAVE_INTREPID_DEBUG
00704   TEUCHOS_TEST_FOR_EXCEPTION( ( this -> rank() != 3), std::invalid_argument, 
00705                       ">>> ERROR (FieldContainer): Number of indices does not match rank of the container.");  
00706   TEUCHOS_TEST_FOR_EXCEPTION( ( (i0 < 0) || ( i0 >= dim0_) ), std::invalid_argument,
00707                       ">>> ERROR (FieldContainer): 1st index is out of range.");    
00708   TEUCHOS_TEST_FOR_EXCEPTION( ( (i1 < 0) || (i1 >= dim1_) ), std::invalid_argument,
00709                       ">>> ERROR (FieldContainer): 2nd index is out of range.");    
00710   TEUCHOS_TEST_FOR_EXCEPTION( ( (i2 < 0) || (i2 >= dim2_) ), std::invalid_argument,
00711                       ">>> ERROR (FieldContainer): 3rd index is out of range.");    
00712 #endif
00713   return (i0*dim1_ + i1)*dim2_ + i2;
00714 }
00715 
00716 
00717 
00718 template<class Scalar, int ArrayTypeId>
00719 inline int FieldContainer<Scalar, ArrayTypeId>::getEnumeration(const int i0,
00720                                                   const int i1,
00721                                                   const int i2,
00722                                                   const int i3) const {
00723 #ifdef HAVE_INTREPID_DEBUG
00724   TEUCHOS_TEST_FOR_EXCEPTION( ( this -> rank() != 4), std::invalid_argument, 
00725                       ">>> ERROR (FieldContainer): Number of indices does not match rank of the container.");  
00726   TEUCHOS_TEST_FOR_EXCEPTION( ( (i0 < 0) || ( i0 >= dim0_) ), std::invalid_argument,
00727                       ">>> ERROR (FieldContainer): 1st index is out of range.");    
00728   TEUCHOS_TEST_FOR_EXCEPTION( ( (i1 < 0) || (i1 >= dim1_) ), std::invalid_argument,
00729                       ">>> ERROR (FieldContainer): 2nd index is out of range.");    
00730   TEUCHOS_TEST_FOR_EXCEPTION( ( (i2 < 0) || (i2 >= dim2_) ), std::invalid_argument,
00731                       ">>> ERROR (FieldContainer): 3rd index is out of range.");    
00732   TEUCHOS_TEST_FOR_EXCEPTION( ( (i3 < 0) || (i3 >= dim3_) ), std::invalid_argument,
00733                       ">>> ERROR (FieldContainer): 4th index is out of range.");    
00734 #endif
00735   return ( (i0*dim1_ + i1 )*dim2_ + i2 )*dim3_ + i3;
00736 }
00737 
00738 
00739 
00740 template<class Scalar, int ArrayTypeId>
00741 inline int FieldContainer<Scalar, ArrayTypeId>::getEnumeration(const int i0,
00742                                                   const int i1,
00743                                                   const int i2,
00744                                                   const int i3,
00745                                                   const int i4) const {
00746 #ifdef HAVE_INTREPID_DEBUG
00747   TEUCHOS_TEST_FOR_EXCEPTION( ( this -> rank() != 5), std::invalid_argument, 
00748                       ">>> ERROR (FieldContainer): Number of indices does not match rank of the container.");  
00749   TEUCHOS_TEST_FOR_EXCEPTION( ( (i0 < 0) || ( i0 >= dim0_) ), std::invalid_argument,
00750                       ">>> ERROR (FieldContainer): 1st index is out of range.");    
00751   TEUCHOS_TEST_FOR_EXCEPTION( ( (i1 < 0) || (i1 >= dim1_) ), std::invalid_argument,
00752                       ">>> ERROR (FieldContainer): 2nd index is out of range.");    
00753   TEUCHOS_TEST_FOR_EXCEPTION( ( (i2 < 0) || (i2 >= dim2_) ), std::invalid_argument,
00754                       ">>> ERROR (FieldContainer): 3rd index is out of range.");    
00755   TEUCHOS_TEST_FOR_EXCEPTION( ( (i3 < 0) || (i3 >= dim3_) ), std::invalid_argument,
00756                       ">>> ERROR (FieldContainer): 4th index is out of range.");    
00757   TEUCHOS_TEST_FOR_EXCEPTION( ( (i4 < 0) || (i4 >= dim4_) ), std::invalid_argument,
00758                       ">>> ERROR (FieldContainer): 5th index is out of range.");    
00759 #endif
00760   return ( ( (i0*dim1_ + i1 )*dim2_ + i2 )*dim3_ + i3 )*dim4_ + i4;
00761 }
00762 
00763 
00764 
00765 
00766 template<class Scalar, int ArrayTypeId>
00767 int FieldContainer<Scalar, ArrayTypeId>::getEnumeration(const Teuchos::Array<int>& multiIndex) const {
00768 
00769 #ifdef HAVE_INTREPID_DEBUG
00770   // Check if number of multi-indices matches rank of the FieldContainer object
00771   TEUCHOS_TEST_FOR_EXCEPTION( ( multiIndex.size() != dimensions_.size() ),
00772                       std::invalid_argument,
00773                       ">>> ERROR (FieldContainer): Number of multi-indices does not match rank of container.");
00774   TEUCHOS_TEST_FOR_EXCEPTION( ( ( multiIndex[0] < 0) || ( multiIndex[0] >= dim0_) ),
00775                       std::invalid_argument,
00776                       ">>> ERROR (FieldContainer): 1st index is out of range.");    
00777 #endif  
00778 
00779   int rank = dimensions_.size();
00780   int address = 0;
00781   switch (rank) {
00782     
00783     // Optimize enumeration computation for low rank (<= 5) containers
00784     case 5:
00785 #ifdef HAVE_INTREPID_DEBUG
00786       TEUCHOS_TEST_FOR_EXCEPTION( ( (multiIndex[4] < 0) || (multiIndex[4] >= dim4_) ),
00787                           std::invalid_argument,
00788                           ">>> ERROR (FieldContainer): 5th index is out of range.");    
00789       TEUCHOS_TEST_FOR_EXCEPTION( ( (multiIndex[3] < 0) || (multiIndex[3] >= dim3_) ),
00790                           std::invalid_argument,
00791                           ">>> ERROR (FieldContainer): 4th index is out of range.");    
00792       TEUCHOS_TEST_FOR_EXCEPTION( ( (multiIndex[2] < 0) || (multiIndex[2] >= dim2_) ),
00793                           std::invalid_argument,
00794                           ">>> ERROR (FieldContainer): 3rd index is out of range.");    
00795       TEUCHOS_TEST_FOR_EXCEPTION( ( (multiIndex[1] < 0) || (multiIndex[1] >= dim1_) ),
00796                           std::invalid_argument,
00797                           ">>> ERROR (FieldContainer): 2nd index is out of range.");    
00798 #endif
00799       address = (((multiIndex[0]*dim1_ + multiIndex[1])*dim2_ + multiIndex[2])*dim3_ + multiIndex[3])*dim4_ + multiIndex[4];
00800       break;
00801       
00802     case 4:
00803 #ifdef HAVE_INTREPID_DEBUG
00804       TEUCHOS_TEST_FOR_EXCEPTION( ( (multiIndex[3] < 0) || (multiIndex[3] >= dim3_) ),
00805                           std::invalid_argument,
00806                           ">>> ERROR (FieldContainer): 4th index is out of range.");    
00807       TEUCHOS_TEST_FOR_EXCEPTION( ( (multiIndex[2] < 0) || (multiIndex[2] >= dim2_) ),
00808                           std::invalid_argument,
00809                           ">>> ERROR (FieldContainer): 3rd index is out of range.");    
00810       TEUCHOS_TEST_FOR_EXCEPTION( ( (multiIndex[1] < 0) || (multiIndex[1] >= dim1_) ),
00811                           std::invalid_argument,
00812                           ">>> ERROR (FieldContainer): 2nd index is out of range.");    
00813 #endif
00814       address = ((multiIndex[0]*dim1_ + multiIndex[1])*dim2_ + multiIndex[2])*dim3_ + multiIndex[3];
00815       break;
00816 
00817     case 3:
00818 #ifdef HAVE_INTREPID_DEBUG
00819       TEUCHOS_TEST_FOR_EXCEPTION( ( (multiIndex[2] < 0) || (multiIndex[2] >= dim2_) ),
00820                           std::invalid_argument,
00821                           ">>> ERROR (FieldContainer): 3rd index is out of range.");    
00822       TEUCHOS_TEST_FOR_EXCEPTION( ( (multiIndex[1] < 0) || (multiIndex[1] >= dim1_) ),
00823                           std::invalid_argument,
00824                           ">>> ERROR (FieldContainer): 2nd index is out of range.");    
00825 #endif
00826       address = (multiIndex[0]*dim1_ + multiIndex[1])*dim2_ + multiIndex[2];
00827       break;
00828 
00829     case 2:
00830 #ifdef HAVE_INTREPID_DEBUG
00831       TEUCHOS_TEST_FOR_EXCEPTION( ( (multiIndex[1] < 0) || (multiIndex[1] >= dim1_) ),
00832                           std::invalid_argument,
00833                           ">>> ERROR (FieldContainer): 2nd index is out of range.");    
00834 #endif
00835       address = multiIndex[0]*dim1_ + multiIndex[1];
00836       break;
00837       
00838     case 1:
00839       address = multiIndex[0];
00840       break;
00841 
00842     default:
00843       
00844       // Arbitrary rank: compute address using Horner's nested scheme: intialize address to 0th index value
00845       address = multiIndex[0];
00846       for (int r = 0; r < rank - 1; r++){
00847 #ifdef HAVE_INTREPID_DEBUG
00848         TEUCHOS_TEST_FOR_EXCEPTION( ( (multiIndex[r+1] < 0) || (multiIndex[r+1] >= dimensions_[r+1]) ),
00849                             std::invalid_argument,
00850                             ">>> ERROR (FieldContainer): Multi-index component out of range.");    
00851 #endif
00852         // Add increment
00853         address = address*dimensions_[r+1] + multiIndex[r+1];
00854       }
00855   } // end switch(rank)
00856 
00857   return address;
00858 }
00859 
00860 
00861 
00862 template<class Scalar, int ArrayTypeId>
00863 void FieldContainer<Scalar, ArrayTypeId>::getMultiIndex(int & i0,
00864                                            const int valueEnum) const 
00865 {
00866 #ifdef HAVE_INTREPID_DEBUG
00867   TEUCHOS_TEST_FOR_EXCEPTION( ( this -> rank() != 1), std::invalid_argument, 
00868                       ">>> ERROR (FieldContainer): Number of indices does not match rank of the container.");  
00869   TEUCHOS_TEST_FOR_EXCEPTION( ( (valueEnum < 0) || (valueEnum >= (int)data_.size()) ),
00870                       std::invalid_argument,
00871                       ">>> ERROR (FieldContainer): Value enumeration is out of range.");    
00872 #endif
00873   i0 = valueEnum;
00874 }
00875 
00876 
00877 
00878 template<class Scalar, int ArrayTypeId>
00879 void FieldContainer<Scalar, ArrayTypeId>::getMultiIndex(int & i0,
00880                                            int & i1,
00881                                            const int valueEnum) const 
00882 {
00883 #ifdef HAVE_INTREPID_DEBUG
00884   TEUCHOS_TEST_FOR_EXCEPTION( ( this -> rank() != 2), std::invalid_argument, 
00885                       ">>> ERROR (FieldContainer): Number of indices does not match rank of the container.");  
00886   TEUCHOS_TEST_FOR_EXCEPTION( ( (valueEnum < 0) || (valueEnum >= (int)data_.size()) ),
00887                       std::invalid_argument,
00888                       ">>> ERROR (FieldContainer): Value enumeration is out of range.");    
00889 #endif
00890   
00891   i0 = valueEnum/dim1_;
00892   i1 = valueEnum - i0*dim1_;
00893 }
00894 
00895 
00896 
00897 template<class Scalar, int ArrayTypeId>
00898 void FieldContainer<Scalar, ArrayTypeId>::getMultiIndex(int & i0,
00899                                            int & i1,
00900                                            int & i2,
00901                                            const int valueEnum) const 
00902 {
00903 #ifdef HAVE_INTREPID_DEBUG
00904   TEUCHOS_TEST_FOR_EXCEPTION( ( this -> rank() != 3), std::invalid_argument, 
00905                       ">>> ERROR (FieldContainer): Number of indices does not match rank of the container.");  
00906   TEUCHOS_TEST_FOR_EXCEPTION( ( (valueEnum < 0) || (valueEnum >= (int)data_.size()) ),
00907                       std::invalid_argument,
00908                       ">>> ERROR (FieldContainer): Value enumeration is out of range.");    
00909 #endif
00910   int tempDim = dim1_*dim2_;
00911   int tempEnu = valueEnum;
00912   i0 = tempEnu/tempDim;
00913   
00914   tempEnu -= i0*tempDim;
00915   tempDim /= dim1_;
00916   i1 = tempEnu/tempDim;
00917   
00918   tempEnu -= i1*tempDim;
00919   tempDim /= dim2_;
00920   i2 = tempEnu/tempDim;
00921 }
00922 
00923 
00924 
00925 template<class Scalar, int ArrayTypeId>
00926 void FieldContainer<Scalar, ArrayTypeId>::getMultiIndex(int & i0,
00927                                            int & i1,
00928                                            int & i2,
00929                                            int & i3,
00930                                            const int valueEnum) const 
00931 {
00932 #ifdef HAVE_INTREPID_DEBUG
00933   TEUCHOS_TEST_FOR_EXCEPTION( ( this -> rank() != 4), std::invalid_argument, 
00934                       ">>> ERROR (FieldContainer): Number of indices does not match rank of the container.");  
00935   TEUCHOS_TEST_FOR_EXCEPTION( ( (valueEnum < 0) || (valueEnum >= (int)data_.size()) ),
00936                       std::invalid_argument,
00937                       ">>> ERROR (FieldContainer): Value enumeration is out of range.");    
00938 #endif
00939   int tempDim = dim1_*dim2_*dim3_;
00940   int tempEnu = valueEnum;
00941   i0 = tempEnu/tempDim;
00942   
00943   tempEnu -= i0*tempDim;
00944   tempDim /= dim1_;
00945   i1 = tempEnu/tempDim;
00946   
00947   tempEnu -= i1*tempDim;
00948   tempDim /= dim2_;
00949   i2 = tempEnu/tempDim;
00950   
00951   tempEnu -= i2*tempDim;
00952   tempDim /= dim3_;
00953   i3 = tempEnu/tempDim;
00954 }
00955 
00956 
00957 
00958 
00959 template<class Scalar, int ArrayTypeId>
00960 void FieldContainer<Scalar, ArrayTypeId>::getMultiIndex(int & i0,
00961                                            int & i1,
00962                                            int & i2,
00963                                            int & i3,
00964                                            int & i4,
00965                                            const int valueEnum) const 
00966 {
00967 #ifdef HAVE_INTREPID_DEBUG
00968   TEUCHOS_TEST_FOR_EXCEPTION( ( this -> rank() != 5), std::invalid_argument, 
00969                       ">>> ERROR (FieldContainer): Number of indices does not match rank of the container.");  
00970   TEUCHOS_TEST_FOR_EXCEPTION( ( (valueEnum < 0) || (valueEnum >= (int)data_.size()) ),
00971                       std::invalid_argument,
00972                       ">>> ERROR (FieldContainer): Value enumeration is out of range.");    
00973 #endif
00974   int tempDim = dim1_*dim2_*dim3_*dim4_;
00975   int tempEnu = valueEnum;
00976   i0 = tempEnu/tempDim;
00977   
00978   tempEnu -= i0*tempDim;
00979   tempDim /= dim1_;
00980   i1 = tempEnu/tempDim;
00981   
00982   tempEnu -= i1*tempDim;
00983   tempDim /= dim2_;
00984   i2 = tempEnu/tempDim;
00985   
00986   tempEnu -= i2*tempDim;
00987   tempDim /= dim3_;
00988   i3 = tempEnu/tempDim;
00989 
00990   tempEnu -= i3*tempDim;
00991   tempDim /= dim4_;
00992   i4 = tempEnu/tempDim;
00993 }
00994 
00995 
00996 
00997 template<class Scalar, int ArrayTypeId>
00998 template<class Vector>
00999 void FieldContainer<Scalar, ArrayTypeId>::getMultiIndex(Vector &             multiIndex,
01000                                            const int            valueEnum) const 
01001 {
01002   
01003   // Verify address is in the admissible range for this FieldContainer
01004 #ifdef HAVE_INTREPID_DEBUG
01005   TEUCHOS_TEST_FOR_EXCEPTION( ( (valueEnum < 0) || (valueEnum >= (int)data_.size()) ),
01006                       std::invalid_argument,
01007                       ">>> ERROR (FieldContainer): Value enumeration is out of range.");    
01008 #endif
01009   
01010   // make sure multiIndex has the right size to hold all multi-indices
01011   int rank = dimensions_.size();
01012   multiIndex.resize(rank);
01013   
01014   // Initializations
01015   int temp_enum = valueEnum;
01016   int temp_range = 1;
01017   
01018   // Compute product of all but the first upper bound
01019   for(int r = 1; r < rank ; r++){
01020     temp_range *=dimensions_[r];
01021   }
01022   
01023   // Index 0 is computed first using integer division
01024   multiIndex[0] = temp_enum/temp_range;
01025   
01026   // Indices 1 to (rank - 2) are computed next; will be skipped if rank <=2
01027   for(int r = 1; r < rank - 1; r++){
01028     temp_enum  -= multiIndex[r-1]*temp_range;
01029     temp_range /= dimensions_[r];
01030     multiIndex[r] = temp_enum/temp_range;
01031   }
01032   
01033   // Index (rank - 1) is computed last, skip if rank = 1 and keep if rank = 2
01034   if(rank > 1) {
01035     multiIndex[rank - 1] = temp_enum - multiIndex[rank - 2]*temp_range;
01036   }
01037 }
01038 
01039 //--------------------------------------------------------------------------------------------//
01040 //                                                                                            //
01041 //                          Methods to shape (resize) a field container                       //
01042 //                                                                                            //
01043 //--------------------------------------------------------------------------------------------//
01044 
01045 template<class Scalar, int ArrayTypeId>
01046 inline void FieldContainer<Scalar, ArrayTypeId>::clear() {
01047   dimensions_.resize(0);
01048   
01049   // Reset first five dimensions:
01050   dim0_ = 0;
01051   dim1_ = 0;
01052   dim2_ = 0;
01053   dim3_ = 0;
01054   dim4_ = 0;
01055   
01056   // Clears data array and sets to zero length
01057   data_.clear();
01058   data_ptr_ = data_.begin();
01059 }
01060 
01061 
01062 
01063 template<class Scalar, int ArrayTypeId>
01064 void FieldContainer<Scalar, ArrayTypeId>::resize(const Teuchos::Array<int>& newDimensions) {
01065   
01066   // First handle the trivial case of zero dimensions
01067   if( newDimensions.size() == 0) {
01068     dimensions_.resize(0);
01069     dim0_ = 0;
01070     dim1_ = 0;
01071     dim2_ = 0;
01072     dim3_ = 0;
01073     dim4_ = 0;
01074     data_.resize(0);
01075     data_ptr_ = data_.begin();
01076   }
01077   else {
01078     
01079     // Copy upper index bounds and resize container storage to match new upper bounds.
01080     dimensions_.assign(newDimensions.begin(),newDimensions.end());  
01081     
01082     // Copy first 5 dimensions for faster access
01083     unsigned int rank = dimensions_.size();
01084     switch(rank) {
01085       case 1:
01086         dim0_ = dimensions_[0]; 
01087         dim1_ = 0;
01088         dim2_ = 0;
01089         dim3_ = 0;
01090         dim4_ = 0;
01091         break;
01092         
01093       case 2:
01094         dim0_ = dimensions_[0]; 
01095         dim1_ = dimensions_[1]; 
01096         dim2_ = 0;
01097         dim3_ = 0;
01098         dim4_ = 0;
01099         break;
01100         
01101       case 3:
01102         dim0_ = dimensions_[0]; 
01103         dim1_ = dimensions_[1]; 
01104         dim2_ = dimensions_[2]; 
01105         dim3_ = 0;
01106         dim4_ = 0;
01107         break;
01108         
01109       case 4:
01110         dim0_ = dimensions_[0]; 
01111         dim1_ = dimensions_[1]; 
01112         dim2_ = dimensions_[2]; 
01113         dim3_ = dimensions_[3]; 
01114         dim4_ = 0;
01115         break;
01116         
01117       case 5:
01118       default:
01119         dim0_ = dimensions_[0]; 
01120         dim1_ = dimensions_[1]; 
01121         dim2_ = dimensions_[2]; 
01122         dim3_ = dimensions_[3]; 
01123         dim4_ = dimensions_[4]; 
01124     }
01125     
01126     // Resize data array
01127     data_.resize(this -> size());
01128     data_ptr_ = data_.begin();
01129   }
01130 }
01131 
01132 
01133 
01134 template<class Scalar, int ArrayTypeId>
01135 inline void FieldContainer<Scalar, ArrayTypeId>::resize(const int dim0) {
01136   dim0_ = dim0;  
01137   dim1_ = 0;  
01138   dim2_ = 0;  
01139   dim3_ = 0;  
01140   dim4_ = 0;
01141   dimensions_.resize(1);  
01142   dimensions_[0] = dim0_; 
01143   data_.resize(dim0_); 
01144   data_ptr_ = data_.begin();
01145 }
01146 
01147 
01148 
01149 template<class Scalar, int ArrayTypeId>
01150 inline void FieldContainer<Scalar, ArrayTypeId>::resize(const int dim0,
01151                                            const int dim1) {
01152   dim0_ = dim0;  
01153   dim1_ = dim1;  
01154   dim2_ = 0;  
01155   dim3_ = 0;  
01156   dim4_ = 0;
01157   dimensions_.resize(2);  
01158   dimensions_[0] = dim0_;  
01159   dimensions_[1] = dim1_;  
01160   data_.resize(dim0_*dim1_); 
01161   data_ptr_ = data_.begin();
01162 }
01163 
01164 
01165 
01166 template<class Scalar, int ArrayTypeId>
01167 inline void FieldContainer<Scalar, ArrayTypeId>::resize(const int dim0,
01168                                            const int dim1,
01169                                            const int dim2) {
01170   dim0_ = dim0;
01171   dim1_ = dim1;
01172   dim2_ = dim2;
01173   dim3_ = 0;  
01174   dim4_ = 0;
01175   dimensions_.resize(3);  
01176   dimensions_[0] = dim0_; 
01177   dimensions_[1] = dim1_;  
01178   dimensions_[2] = dim2_;  
01179   data_.resize(dim0_*dim1_*dim2_); 
01180   data_ptr_ = data_.begin();
01181 }
01182 
01183 
01184 
01185 template<class Scalar, int ArrayTypeId>
01186 inline void FieldContainer<Scalar, ArrayTypeId>::resize(const int dim0,
01187                                            const int dim1,
01188                                            const int dim2,
01189                                            const int dim3) {
01190   dim0_ = dim0;
01191   dim1_ = dim1;
01192   dim2_ = dim2;
01193   dim3_ = dim3;
01194   dim4_ = 0;
01195   dimensions_.resize(4);  
01196   dimensions_[0] = dim0_;  
01197   dimensions_[1] = dim1_;  
01198   dimensions_[2] = dim2_;  
01199   dimensions_[3] = dim3_;  
01200   data_.resize(dim0_*dim1_*dim2_*dim3_); 
01201   data_ptr_ = data_.begin();
01202 }
01203 
01204 
01205 
01206 template<class Scalar, int ArrayTypeId>
01207 inline void FieldContainer<Scalar, ArrayTypeId>::resize(const int dim0,
01208                                            const int dim1,
01209                                            const int dim2,
01210                                            const int dim3,
01211                                            const int dim4) {
01212   dim0_ = dim0;
01213   dim1_ = dim1;
01214   dim2_ = dim2;
01215   dim3_ = dim3;
01216   dim4_ = dim4;
01217   dimensions_.resize(5);  
01218   dimensions_[0] = dim0_;  
01219   dimensions_[1] = dim1_;  
01220   dimensions_[2] = dim2_;  
01221   dimensions_[3] = dim3_;  
01222   dimensions_[4] = dim4_;  
01223   data_.resize(dim0_*dim1_*dim2_*dim3_*dim4_); 
01224   data_ptr_ = data_.begin();
01225 }
01226 
01227 
01228 
01229 template<class Scalar, int ArrayTypeId>
01230 inline void FieldContainer<Scalar, ArrayTypeId>::resize(const FieldContainer<Scalar, ArrayTypeId>& anotherContainer) {
01231   
01232   // Copy dimensions from the specified container
01233   anotherContainer.dimensions(dimensions_);
01234   int newRank = dimensions_.size();
01235   
01236   // Copy first 5 dimensions for faster access
01237   switch(newRank) {
01238     case 1:
01239       dim0_ = dimensions_[0]; 
01240       dim1_ = 0;
01241       dim2_ = 0;
01242       dim3_ = 0;
01243       dim4_ = 0;
01244       break;
01245       
01246     case 2:
01247       dim0_ = dimensions_[0]; 
01248       dim1_ = dimensions_[1]; 
01249       dim2_ = 0;
01250       dim3_ = 0;
01251       dim4_ = 0;
01252       break;
01253       
01254     case 3:
01255       dim0_ = dimensions_[0]; 
01256       dim1_ = dimensions_[1]; 
01257       dim2_ = dimensions_[2]; 
01258       dim3_ = 0;
01259       dim4_ = 0;
01260       break;
01261       
01262     case 4:
01263       dim0_ = dimensions_[0]; 
01264       dim1_ = dimensions_[1]; 
01265       dim2_ = dimensions_[2]; 
01266       dim3_ = dimensions_[3]; 
01267       dim4_ = 0;
01268       break;
01269       
01270     case 5:
01271     default:
01272       dim0_ = dimensions_[0]; 
01273       dim1_ = dimensions_[1]; 
01274       dim2_ = dimensions_[2]; 
01275       dim3_ = dimensions_[3]; 
01276       dim4_ = dimensions_[4]; 
01277   }
01278   
01279   // Resize data array
01280   data_.resize(this->size());  
01281   data_ptr_ = data_.begin();
01282 }
01283 
01284 
01285 template<class Scalar, int ArrayTypeId>
01286 void FieldContainer<Scalar, ArrayTypeId>::resize(const int             numPoints,
01287                                     const int             numFields,
01288                                     const EFunctionSpace  spaceType,
01289                                     const EOperator       operatorType,
01290                                     const int             spaceDim) {  
01291   // Validate input
01292 #ifdef HAVE_INTREPID_DEBUG
01293   TEUCHOS_TEST_FOR_EXCEPTION( ( numPoints < 0),
01294                       std::invalid_argument,
01295                       ">>> ERROR (FieldContainer): Number of points cannot be negative!");  
01296   TEUCHOS_TEST_FOR_EXCEPTION( ( numFields < 0),
01297                       std::invalid_argument,
01298                       ">>> ERROR (FieldContainer): Number of fields cannot be negative!");  
01299   TEUCHOS_TEST_FOR_EXCEPTION( !( (1 <=  spaceDim ) && ( spaceDim <= 3  ) ),
01300                       std::invalid_argument,
01301                       ">>> ERROR (FieldContainer): Invalid space dimension.");  
01302 #endif  
01303   
01304   // Find out field and operator ranks
01305   int fieldRank    = getFieldRank(spaceType);
01306   int operatorRank = getOperatorRank(spaceType,operatorType,spaceDim);
01307   
01308   // Compute rank of the container = 1(numPoints) + 1(numFields) + fieldRank + operatorRank
01309   int rank = 1 + 1 + fieldRank + operatorRank;
01310   
01311   // Define temp array for the dimensions
01312   Teuchos::Array<int> newDimensions(rank);
01313   
01314   // Dimensions 0 and 1 are number of points and number of fields, resp.
01315   newDimensions[0] = numPoints;
01316   newDimensions[1] = numFields;
01317   
01318   // The rest of the dimensions depend on whether we had VALUE, GRAD (D1), CURL, DIV or Dk, k>1
01319   switch(operatorType) {
01320     
01321     case OPERATOR_VALUE:
01322     case OPERATOR_GRAD:
01323     case OPERATOR_D1:
01324     case OPERATOR_CURL:
01325     case OPERATOR_DIV:
01326       
01327       // For these operators all dimensions from 2 to 2 + fieldRank + OperatorRank are bounded by spaceDim
01328       for(int i = 0; i < fieldRank + operatorRank; i++){
01329         newDimensions[2 + i] = spaceDim; 
01330       }
01331       break;
01332       
01333     case OPERATOR_D2:
01334     case OPERATOR_D3:
01335     case OPERATOR_D4:
01336     case OPERATOR_D5:
01337     case OPERATOR_D6:
01338     case OPERATOR_D7:
01339     case OPERATOR_D8:
01340     case OPERATOR_D9:
01341     case OPERATOR_D10:
01342       
01343       // All dimensions from 2 to 2 + fieldRank, if any, are bounded by spaceDim
01344       for(int i = 0; i < fieldRank; i++){
01345         newDimensions[2 + i] = spaceDim; 
01346       }
01347       
01348       // We know that for Dk operatorRank = 1 and so there's just one more dimension left
01349       // given by the cardinality of the set of all derivatives of order k
01350       newDimensions[2 + fieldRank] = getDkCardinality(operatorType,spaceDim);
01351       break;
01352       
01353     default:
01354       TEUCHOS_TEST_FOR_EXCEPTION( !(Intrepid::isValidOperator(operatorType) ), std::invalid_argument,
01355                           ">>> ERROR (FieldContainer): Invalid operator type");    
01356   }
01357   
01358   // Resize FieldContainer using the newDimensions in the array
01359   this -> resize(newDimensions);
01360 }
01361 
01362 //--------------------------------------------------------------------------------------------//
01363 //                                                                                            //
01364 //                     Methods to read and write values to FieldContainer                     //
01365 //                                                                                            //
01366 //--------------------------------------------------------------------------------------------//
01367 
01368 
01369 
01370 template<class Scalar, int ArrayTypeId>
01371 inline void FieldContainer<Scalar, ArrayTypeId>::initialize(const Scalar value) {
01372   for (int i=0; i < this->size(); i++) {
01373     data_[i] = value;
01374   } 
01375 }
01376 
01377 
01378 
01379 template<class Scalar, int ArrayTypeId>
01380 inline Scalar FieldContainer<Scalar, ArrayTypeId>::getValue(const Teuchos::Array<int>& multiIndex) const {
01381   return data_[this -> getEnumeration(multiIndex)];
01382 }
01383 
01384 
01385 
01386 template<class Scalar, int ArrayTypeId>
01387 inline void FieldContainer<Scalar, ArrayTypeId>::setValue(const Scalar dataValue, 
01388                                              const Teuchos::Array<int>& multiIndex) {
01389   data_[this -> getEnumeration(multiIndex)] = dataValue; 
01390 }
01391 
01392 
01393 
01394 template<class Scalar, int ArrayTypeId>
01395 inline void FieldContainer<Scalar, ArrayTypeId>::setValue(const Scalar dataValue, 
01396                                              const int    order) {
01397   data_[order] = dataValue; 
01398 }
01399 
01400 
01401 
01402 template<class Scalar, int ArrayTypeId>
01403 void FieldContainer<Scalar, ArrayTypeId>::setValues(const Teuchos::ArrayView<Scalar>& dataArray) {
01404 #ifdef HAVE_INTREPID_DEBUG
01405   TEUCHOS_TEST_FOR_EXCEPTION( (dataArray.size() != (data_.size()) ),
01406                       std::invalid_argument,
01407                       ">>> ERROR (FieldContainer): Size of argument does not match the size of container.");  
01408 #endif  
01409   data_.assign(dataArray.begin(),dataArray.end());
01410   data_ptr_ = data_.begin();
01411 }
01412 
01413 
01414 
01415 template<class Scalar, int ArrayTypeId>
01416 void FieldContainer<Scalar, ArrayTypeId>::setValues(const Scalar* dataPtr, 
01417                                        const int numData) 
01418 {
01419 #ifdef HAVE_INTREPID_DEBUG
01420   TEUCHOS_TEST_FOR_EXCEPTION( (numData != this -> size() ), std::invalid_argument,
01421                       ">>> ERROR (FieldContainer): Number of data does not match the size of container.");  
01422 
01423 #endif
01424   data_.assign(dataPtr, dataPtr + numData);  
01425   data_ptr_ = data_.begin();
01426 }
01427 
01428 
01429 
01430 template<class Scalar, int ArrayTypeId>
01431 inline const Scalar& FieldContainer<Scalar, ArrayTypeId>::operator () (const int i0) const 
01432 {
01433 #ifdef HAVE_INTREPID_DEBUG
01434   TEUCHOS_TEST_FOR_EXCEPTION( ( this -> rank() != 1), std::invalid_argument, 
01435                       ">>> ERROR (FieldContainer): Number of indices does not match rank of the container.");  
01436   TEUCHOS_TEST_FOR_EXCEPTION( ( (i0 < 0) || (i0 >= dim0_) ), std::invalid_argument,
01437                       ">>> ERROR (FieldContainer): index is out of range.");    
01438 #endif
01439   return data_ptr_[i0]; 
01440 }
01441 
01442 
01443 template<class Scalar, int ArrayTypeId>
01444 inline Scalar& FieldContainer<Scalar, ArrayTypeId>::operator () (const int i0)  
01445 {
01446 #ifdef HAVE_INTREPID_DEBUG
01447   TEUCHOS_TEST_FOR_EXCEPTION( ( this -> rank() != 1), std::invalid_argument, 
01448                       ">>> ERROR (FieldContainer): Number of indices does not match rank of the container.");  
01449   TEUCHOS_TEST_FOR_EXCEPTION( ( (i0 < 0) || (i0 >= dim0_) ), std::invalid_argument,
01450                       ">>> ERROR (FieldContainer): index is out of range.");    
01451 #endif
01452   return data_ptr_[i0]; 
01453 }
01454 
01455 
01456 
01457 template<class Scalar, int ArrayTypeId>
01458 inline const Scalar& FieldContainer<Scalar, ArrayTypeId>::operator () (const int i0,
01459                                                           const int i1) const 
01460 {
01461 #ifdef HAVE_INTREPID_DEBUG
01462   TEUCHOS_TEST_FOR_EXCEPTION( ( this -> rank() != 2), std::invalid_argument, 
01463                       ">>> ERROR (FieldContainer): Number of indices does not match rank of the container.");  
01464   TEUCHOS_TEST_FOR_EXCEPTION( ( (i0 < 0) || (i0 >= dim0_) ), std::invalid_argument,
01465                       ">>> ERROR (FieldContainer): 1st index is out of range.");    
01466   TEUCHOS_TEST_FOR_EXCEPTION( ( (i1 < 0) || (i1 >= dim1_) ), std::invalid_argument,
01467                       ">>> ERROR (FieldContainer): 2nd index is out of range.");    
01468 #endif
01469   return data_ptr_[i0*dim1_ + i1]; 
01470 }
01471 
01472 
01473 template<class Scalar, int ArrayTypeId>
01474 inline Scalar& FieldContainer<Scalar, ArrayTypeId>::operator () (const int i0,
01475                                                     const int i1)  
01476 {
01477 #ifdef HAVE_INTREPID_DEBUG
01478   TEUCHOS_TEST_FOR_EXCEPTION( ( this -> rank() != 2), std::invalid_argument, 
01479                       ">>> ERROR (FieldContainer): Number of indices does not match rank of the container.");  
01480   TEUCHOS_TEST_FOR_EXCEPTION( ( (i0 < 0) || (i0 >= dim0_) ), std::invalid_argument,
01481                       ">>> ERROR (FieldContainer): 1st index is out of range.");    
01482   TEUCHOS_TEST_FOR_EXCEPTION( ( (i1 < 0) || (i1 >= dim1_) ), std::invalid_argument,
01483                       ">>> ERROR (FieldContainer): 2nd index is out of range.");    
01484 #endif
01485   return data_ptr_[i0*dim1_ + i1]; 
01486 }
01487 
01488 
01489 
01490 template<class Scalar, int ArrayTypeId>
01491 inline const Scalar& FieldContainer<Scalar, ArrayTypeId>::operator () (const int i0,
01492                                                           const int i1,
01493                                                           const int i2) const 
01494 {
01495 #ifdef HAVE_INTREPID_DEBUG
01496   TEUCHOS_TEST_FOR_EXCEPTION( ( this -> rank() != 3), std::invalid_argument, 
01497                       ">>> ERROR (FieldContainer): Number of indices does not match rank of the container.");  
01498   TEUCHOS_TEST_FOR_EXCEPTION( ( (i0 < 0) || (i0 >= dim0_) ), std::invalid_argument,
01499                       ">>> ERROR (FieldContainer): 1st index is out of range.");    
01500   TEUCHOS_TEST_FOR_EXCEPTION( ( (i1 < 0) || (i1 >= dim1_) ), std::invalid_argument,
01501                       ">>> ERROR (FieldContainer): 2nd index is out of range.");    
01502   TEUCHOS_TEST_FOR_EXCEPTION( ( (i2 < 0) || (i2 >= dim2_) ), std::invalid_argument,
01503                       ">>> ERROR (FieldContainer): 3rd index is out of range.");    
01504 #endif
01505   return data_ptr_[(i0*dim1_ + i1)*dim2_ + i2]; 
01506 }
01507 
01508 template<class Scalar, int ArrayTypeId>
01509 inline Scalar& FieldContainer<Scalar, ArrayTypeId>::operator () (const int i0,
01510                                                     const int i1,
01511                                                     const int i2) 
01512 {
01513 #ifdef HAVE_INTREPID_DEBUG
01514   TEUCHOS_TEST_FOR_EXCEPTION( ( this -> rank() != 3), std::invalid_argument, 
01515                       ">>> ERROR (FieldContainer): Number of indices does not match rank of the container.");  
01516   TEUCHOS_TEST_FOR_EXCEPTION( ( (i0 < 0) || (i0 >= dim0_) ), std::invalid_argument,
01517                       ">>> ERROR (FieldContainer): 1st index is out of range.");    
01518   TEUCHOS_TEST_FOR_EXCEPTION( ( (i1 < 0) || (i1 >= dim1_) ), std::invalid_argument,
01519                       ">>> ERROR (FieldContainer): 2nd index is out of range.");    
01520   TEUCHOS_TEST_FOR_EXCEPTION( ( (i2 < 0) || (i2 >= dim2_) ), std::invalid_argument,
01521                       ">>> ERROR (FieldContainer): 3rd index is out of range.");    
01522 #endif
01523   return data_ptr_[(i0*dim1_ + i1)*dim2_ + i2]; 
01524 }
01525 
01526 
01527 
01528 template<class Scalar, int ArrayTypeId>
01529 inline const Scalar& FieldContainer<Scalar, ArrayTypeId>::operator ()  (const int i0,
01530                                                            const int i1,
01531                                                            const int i2,
01532                                                            const int i3) const {
01533 #ifdef HAVE_INTREPID_DEBUG
01534   TEUCHOS_TEST_FOR_EXCEPTION( ( this -> rank() != 4), std::invalid_argument, 
01535                       ">>> ERROR (FieldContainer): Number of indices does not match rank of the container.");  
01536   TEUCHOS_TEST_FOR_EXCEPTION( ( (i0 < 0) || (i0 >= dim0_) ), std::invalid_argument,
01537                       ">>> ERROR (FieldContainer): 1st index is out of range.");    
01538   TEUCHOS_TEST_FOR_EXCEPTION( ( (i1 < 0) || (i1 >= dim1_) ), std::invalid_argument,
01539                       ">>> ERROR (FieldContainer): 2nd index is out of range.");    
01540   TEUCHOS_TEST_FOR_EXCEPTION( ( (i2 < 0) || (i2 >= dim2_) ), std::invalid_argument,
01541                       ">>> ERROR (FieldContainer): 3rd index is out of range.");    
01542   TEUCHOS_TEST_FOR_EXCEPTION( ( (i3 < 0) || (i3 >= dim3_) ), std::invalid_argument,
01543                       ">>> ERROR (FieldContainer): 4th index is out of range.");    
01544 #endif
01545   return data_ptr_[( (i0*dim1_ + i1 )*dim2_ + i2 )*dim3_ + i3];
01546 }
01547 
01548 
01549 template<class Scalar, int ArrayTypeId>
01550 inline Scalar& FieldContainer<Scalar, ArrayTypeId>::operator ()  (const int i0,
01551                                                      const int i1,
01552                                                      const int i2,
01553                                                      const int i3) {
01554 #ifdef HAVE_INTREPID_DEBUG
01555   TEUCHOS_TEST_FOR_EXCEPTION( ( this -> rank() != 4), std::invalid_argument, 
01556                       ">>> ERROR (FieldContainer): Number of indices does not match rank of the container.");  
01557   TEUCHOS_TEST_FOR_EXCEPTION( ( (i0 < 0) || (i0 >= dim0_) ), std::invalid_argument,
01558                       ">>> ERROR (FieldContainer): 1st index is out of range.");    
01559   TEUCHOS_TEST_FOR_EXCEPTION( ( (i1 < 0) || (i1 >= dim1_) ), std::invalid_argument,
01560                       ">>> ERROR (FieldContainer): 2nd index is out of range.");    
01561   TEUCHOS_TEST_FOR_EXCEPTION( ( (i2 < 0) || (i2 >= dim2_) ), std::invalid_argument,
01562                       ">>> ERROR (FieldContainer): 3rd index is out of range.");    
01563   TEUCHOS_TEST_FOR_EXCEPTION( ( (i3 < 0) || (i3 >= dim3_) ), std::invalid_argument,
01564                       ">>> ERROR (FieldContainer): 4th index is out of range.");    
01565 #endif
01566   return data_ptr_[( (i0*dim1_ + i1 )*dim2_ + i2 )*dim3_ + i3];
01567 }
01568 
01569 
01570 
01571 template<class Scalar, int ArrayTypeId>
01572 inline const Scalar& FieldContainer<Scalar, ArrayTypeId>::operator ()  (const int i0,
01573                                                            const int i1,
01574                                                            const int i2,
01575                                                            const int i3,
01576                                                            const int i4) const {
01577 #ifdef HAVE_INTREPID_DEBUG
01578   TEUCHOS_TEST_FOR_EXCEPTION( ( this -> rank() != 5), std::invalid_argument, 
01579                       ">>> ERROR (FieldContainer): Number of indices does not match rank of the container.");  
01580   TEUCHOS_TEST_FOR_EXCEPTION( ( (i0 < 0) || (i0 >= dim0_) ), std::invalid_argument,
01581                       ">>> ERROR (FieldContainer): 1st index is out of range.");    
01582   TEUCHOS_TEST_FOR_EXCEPTION( ( (i1 < 0) || (i1 >= dim1_) ), std::invalid_argument,
01583                       ">>> ERROR (FieldContainer): 2nd index is out of range.");    
01584   TEUCHOS_TEST_FOR_EXCEPTION( ( (i2 < 0) || (i2 >= dim2_) ), std::invalid_argument,
01585                       ">>> ERROR (FieldContainer): 3rd index is out of range.");    
01586   TEUCHOS_TEST_FOR_EXCEPTION( ( (i3 < 0) || (i3 >= dim3_) ), std::invalid_argument,
01587                       ">>> ERROR (FieldContainer): 4th index is out of range.");    
01588   TEUCHOS_TEST_FOR_EXCEPTION( ( (i4 < 0) || (i4 >= dim4_) ), std::invalid_argument,
01589                       ">>> ERROR (FieldContainer): 5th index is out of range.");    
01590 #endif
01591   return data_ptr_[( ( (i0*dim1_ + i1 )*dim2_ + i2 )*dim3_ + i3 )*dim4_ + i4];
01592 }
01593 
01594 template<class Scalar, int ArrayTypeId>
01595 inline Scalar& FieldContainer<Scalar, ArrayTypeId>::operator ()  (const int i0,
01596                                                      const int i1,
01597                                                      const int i2,
01598                                                      const int i3,
01599                                                      const int i4) {
01600 #ifdef HAVE_INTREPID_DEBUG
01601   TEUCHOS_TEST_FOR_EXCEPTION( ( this -> rank() != 5), std::invalid_argument, 
01602                       ">>> ERROR (FieldContainer): Number of indices does not match rank of the container.");  
01603   TEUCHOS_TEST_FOR_EXCEPTION( ( (i0 < 0) || (i0 >= dim0_) ), std::invalid_argument,
01604                       ">>> ERROR (FieldContainer): 1st index is out of range.");    
01605   TEUCHOS_TEST_FOR_EXCEPTION( ( (i1 < 0) || (i1 >= dim1_) ), std::invalid_argument,
01606                       ">>> ERROR (FieldContainer): 2nd index is out of range.");    
01607   TEUCHOS_TEST_FOR_EXCEPTION( ( (i2 < 0) || (i2 >= dim2_) ), std::invalid_argument,
01608                       ">>> ERROR (FieldContainer): 3rd index is out of range.");    
01609   TEUCHOS_TEST_FOR_EXCEPTION( ( (i3 < 0) || (i3 >= dim3_) ), std::invalid_argument,
01610                       ">>> ERROR (FieldContainer): 4th index is out of range.");    
01611   TEUCHOS_TEST_FOR_EXCEPTION( ( (i4 < 0) || (i4 >= dim4_) ), std::invalid_argument,
01612                       ">>> ERROR (FieldContainer): 5th index is out of range.");    
01613 #endif
01614   return data_ptr_[( ( (i0*dim1_ + i1 )*dim2_ + i2 )*dim3_ + i3 )*dim4_ + i4];
01615 }
01616 
01617 
01618 
01619 template<class Scalar, int ArrayTypeId>
01620 const Scalar& FieldContainer<Scalar, ArrayTypeId>::operator [] (const int address) const {
01621 #ifdef HAVE_INTREPID_DEBUG
01622   TEUCHOS_TEST_FOR_EXCEPTION( ( (address < 0) || (address >= (int)data_.size() ) ),
01623                       std::invalid_argument,
01624                       ">>> ERROR (FieldContainer): Specified address is out of range.");
01625 #endif
01626   return data_ptr_[address];
01627 }
01628 
01629 
01630 
01631 template<class Scalar, int ArrayTypeId>
01632 Scalar& FieldContainer<Scalar, ArrayTypeId>::operator [] (const int address) {
01633 #ifdef HAVE_INTREPID_DEBUG
01634   TEUCHOS_TEST_FOR_EXCEPTION( ( (address < 0) || (address >= (int)data_.size() ) ),
01635                       std::invalid_argument,
01636                       ">>> ERROR (FieldContainer): Specified address is out of range.");
01637 #endif
01638   return data_ptr_[address];
01639 }
01640 
01641 
01642 
01643 template<class Scalar, int ArrayTypeId>
01644 inline FieldContainer<Scalar, ArrayTypeId>& FieldContainer<Scalar, ArrayTypeId>::operator = (const FieldContainer<Scalar, ArrayTypeId>& right)
01645 {
01646 #ifdef HAVE_INTREPID_DEBUG
01647   TEUCHOS_TEST_FOR_EXCEPTION( ( this == &right ),
01648                       std::invalid_argument,
01649                       ">>> ERROR (FieldContainer): Invalid right-hand side to '='. Self-assignment prohibited.");
01650 #endif
01651   dim0_ = right.dim0_;
01652   dim1_ = right.dim1_;
01653   dim2_ = right.dim2_;
01654   dim3_ = right.dim3_;
01655   dim4_ = right.dim4_;
01656   data_.deepCopy(right.data_());
01657   data_ptr_ = data_.begin();
01658   dimensions_ = right.dimensions_; 
01659   return *this;
01660 }
01661 
01662 
01663 //===========================================================================//
01664 //                                                                           //
01665 //           END of member definitions; START friends and related            //
01666 //                                                                           //
01667 //===========================================================================//
01668 
01669 
01670 template<class Scalar, int ArrayTypeId>
01671 std::ostream& operator << (std::ostream& os, const FieldContainer<Scalar, ArrayTypeId>& container) {
01672   
01673   // Save the format state of the original ostream os.
01674   Teuchos::oblackholestream oldFormatState;
01675   oldFormatState.copyfmt(os);  
01676 
01677   os.setf(std::ios_base::scientific, std::ios_base::floatfield);
01678   os.setf(std::ios_base::right);
01679   int myprec = os.precision();
01680   
01681   int size = container.size();
01682   int rank = container.rank();
01683   Teuchos::Array<int> multiIndex(rank);
01684   Teuchos::Array<int> dimensions;
01685   container.dimensions(dimensions);
01686   
01687   os<< "===============================================================================\n"\
01688     << "\t Container size = " << size << "\n"
01689     << "\t Container rank = " << rank << "\n" ;
01690   
01691   if( (rank == 0 ) && (size == 0) ) {
01692     os<< "====================================================================================\n"\
01693       << "|                        *** This is an empty container ****                       |\n";
01694   }
01695   else {
01696     os<< "\t Dimensions     = ";
01697     
01698     for(int r = 0; r < rank; r++){
01699       os << " (" << dimensions[r] <<") ";
01700     }
01701     os << "\n";
01702     
01703     os<< "====================================================================================\n"\
01704       << "|              Multi-index          Enumeration             Value                  |\n"\
01705       << "====================================================================================\n";
01706   }
01707   
01708   for(int address = 0; address < size; address++){
01709     container.getMultiIndex(multiIndex,address);
01710     std::ostringstream mistring;
01711     for(int r = 0; r < rank; r++){
01712       mistring <<  multiIndex[r] << std::dec << " "; 
01713     }
01714     os.setf(std::ios::right, std::ios::adjustfield);
01715     os << std::setw(27) << mistring.str(); 
01716     os << std::setw(20) << address;
01717     os << "             ";
01718     os.setf(std::ios::left, std::ios::adjustfield);
01719     os << std::setw(myprec+8) << container[address] << "\n";
01720   }
01721   
01722   os<< "====================================================================================\n\n";
01723 
01724   // reset format state of os
01725   os.copyfmt(oldFormatState);
01726 
01727   return os;
01728 }
01729 // End member, friend, and related function definitions of class FieldContainer.
01730 
01731 } // end namespace Intrepid