Intrepid: http://trilinos.sandia.gov/packages/docs/r10.0/packages/intrepid/src/Shared/Intrepid

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