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