Intrepid
Intrepid_FieldContainerDef.hpp
Go to the documentation of this file.
1 // @HEADER
2 // ************************************************************************
3 //
4 // Intrepid Package
5 // Copyright (2007) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Pavel Bochev (pbboche@sandia.gov)
38 // Denis Ridzal (dridzal@sandia.gov), or
39 // Kara Peterson (kjpeter@sandia.gov)
40 //
41 // ************************************************************************
42 // @HEADER
43 
49 namespace Intrepid {
50 
51 
52  //--------------------------------------------------------------------------------------------//
53  // //
54  // Member function definitions of the class FieldContainer //
55  // //
56  //--------------------------------------------------------------------------------------------//
57 
58 
59 template<class Scalar, int ArrayTypeId>
61 
62  // Copy dimensions and data values from right
63  dimensions_.assign(right.dimensions_.begin(),right.dimensions_.end());
64  data_.assign(right.data_.begin(),right.data_.end());
65  dim0_ = right.dim0_;
66  dim1_ = right.dim1_;
67  dim2_ = right.dim2_;
68  dim3_ = right.dim3_;
69  dim4_ = right.dim4_;
70  data_ptr_ = data_.begin();
71 }
72 
73 //--------------------------------------------------------------------------------------------//
74 // //
75 // Constructors of FieldContainer class //
76 // //
77 //--------------------------------------------------------------------------------------------//
78 
79 
80 template<class Scalar, int ArrayTypeId>
81 FieldContainer<Scalar, ArrayTypeId>::FieldContainer(const int dim0) : dim0_(dim0), dim1_(0), dim2_(0), dim3_(0), dim4_(0)
82 {
83  using Teuchos::as;
84  using Teuchos::Ordinal;
85 #ifdef HAVE_INTREPID_DEBUG
86  TEUCHOS_TEST_FOR_EXCEPTION( (0 > dim0), std::invalid_argument,
87  ">>> ERROR (FieldContainer): FieldContainer cannot have a negative dimension.");
88 
89 #endif
90  dimensions_.resize(as<Ordinal>(1));
91  dimensions_[0] = dim0_;
92  data_.assign(as<Ordinal>(dim0_), as<Scalar>(0));
93  data_ptr_ = data_.begin();
94 }
95 
96 
97 
98 template<class Scalar, int ArrayTypeId>
100  const int dim1) : dim0_(dim0), dim1_(dim1), dim2_(0), dim3_(0), dim4_(0)
101 {
102  using Teuchos::as;
103  using Teuchos::Ordinal;
104 #ifdef HAVE_INTREPID_DEBUG
105  TEUCHOS_TEST_FOR_EXCEPTION( (0 > dim0), std::invalid_argument,
106  ">>> ERROR (FieldContainer): FieldContainer cannot have a negative 1st dimension.");
107  TEUCHOS_TEST_FOR_EXCEPTION( (0 > dim1), std::invalid_argument,
108  ">>> ERROR (FieldContainer): FieldContainer cannot have a negative 2nd dimension.");
109 
110 #endif
111  dimensions_.resize(2);
112  dimensions_[0] = dim0_;
113  dimensions_[1] = dim1_;
114  data_.assign(as<Ordinal>(dim0_*dim1_), as<Scalar>(0));
115  data_ptr_ = data_.begin();
116 }
117 
118 
119 
120 template<class Scalar, int ArrayTypeId>
122  const int dim1,
123  const int dim2) : dim0_(dim0), dim1_(dim1), dim2_(dim2), dim3_(0), dim4_(0)
124 {
125  using Teuchos::as;
126  using Teuchos::Ordinal;
127 #ifdef HAVE_INTREPID_DEBUG
128  TEUCHOS_TEST_FOR_EXCEPTION( (0 > dim0), std::invalid_argument,
129  ">>> ERROR (FieldContainer): FieldContainer cannot have a negative 1st dimension.");
130  TEUCHOS_TEST_FOR_EXCEPTION( (0 > dim1), std::invalid_argument,
131  ">>> ERROR (FieldContainer): FieldContainer cannot have a negative 2nd dimension.");
132  TEUCHOS_TEST_FOR_EXCEPTION( (0 > dim2), std::invalid_argument,
133  ">>> ERROR (FieldContainer): FieldContainer cannot have a negative 3rd dimension.");
134 #endif
135  dimensions_.resize(3);
136  dimensions_[0] = dim0_;
137  dimensions_[1] = dim1_;
138  dimensions_[2] = dim2_;
139  data_.assign(as<Ordinal>(dim0_*dim1_*dim2_), as<Scalar>(0));
140  data_ptr_ = data_.begin();
141 }
142 
143 
144 
145 template<class Scalar, int ArrayTypeId>
147  const int dim1,
148  const int dim2,
149  const int dim3) : dim0_(dim0), dim1_(dim1), dim2_(dim2), dim3_(dim3), dim4_(0)
150 {
151  using Teuchos::as;
152  using Teuchos::Ordinal;
153 #ifdef HAVE_INTREPID_DEBUG
154  TEUCHOS_TEST_FOR_EXCEPTION( (0 > dim0), std::invalid_argument,
155  ">>> ERROR (FieldContainer): FieldContainer cannot have a negative 1st dimension.");
156  TEUCHOS_TEST_FOR_EXCEPTION( (0 > dim1), std::invalid_argument,
157  ">>> ERROR (FieldContainer): FieldContainer cannot have a negative 2nd dimension.");
158  TEUCHOS_TEST_FOR_EXCEPTION( (0 > dim2), std::invalid_argument,
159  ">>> ERROR (FieldContainer): FieldContainer cannot have a negative 3rd dimension.");
160  TEUCHOS_TEST_FOR_EXCEPTION( (0 > dim3), std::invalid_argument,
161  ">>> ERROR (FieldContainer): FieldContainer cannot have a negative 4th dimension.");
162 #endif
163  dimensions_.resize(4);
164  dimensions_[0] = dim0_;
165  dimensions_[1] = dim1_;
166  dimensions_[2] = dim2_;
167  dimensions_[3] = dim3_;
168  data_.assign(as<Ordinal>(dim0_*dim1_*dim2_*dim3_), as<Scalar>(0));
169  data_ptr_ = data_.begin();
170 }
171 
172 
173 
174 template<class Scalar, int ArrayTypeId>
176  const int dim1,
177  const int dim2,
178  const int dim3,
179  const int dim4) : dim0_(dim0), dim1_(dim1), dim2_(dim2), dim3_(dim3), dim4_(dim4)
180 {
181  using Teuchos::as;
182  using Teuchos::Ordinal;
183 #ifdef HAVE_INTREPID_DEBUG
184  TEUCHOS_TEST_FOR_EXCEPTION( (0 > dim0), std::invalid_argument,
185  ">>> ERROR (FieldContainer): FieldContainer cannot have a negative 1st dimension.");
186  TEUCHOS_TEST_FOR_EXCEPTION( (0 > dim1), std::invalid_argument,
187  ">>> ERROR (FieldContainer): FieldContainer cannot have a negative 2nd dimension.");
188  TEUCHOS_TEST_FOR_EXCEPTION( (0 > dim2), std::invalid_argument,
189  ">>> ERROR (FieldContainer): FieldContainer cannot have a negative 3rd dimension.");
190  TEUCHOS_TEST_FOR_EXCEPTION( (0 > dim3), std::invalid_argument,
191  ">>> ERROR (FieldContainer): FieldContainer cannot have a negative 4th dimension.");
192  TEUCHOS_TEST_FOR_EXCEPTION( (0 > dim4), std::invalid_argument,
193  ">>> ERROR (FieldContainer): FieldContainer cannot have a negative 5th dimension.");
194 #endif
195  dimensions_.resize(5);
196  dimensions_[0] = dim0_;
197  dimensions_[1] = dim1_;
198  dimensions_[2] = dim2_;
199  dimensions_[3] = dim3_;
200  dimensions_[4] = dim4_;
201  data_.assign(as<Ordinal>(dim0_*dim1_*dim2_*dim3_*dim4_), as<Scalar>(0));
202  data_ptr_ = data_.begin();
203 }
204 
205 
206 
207 template<class Scalar, int ArrayTypeId>
208 FieldContainer<Scalar, ArrayTypeId>::FieldContainer(const Teuchos::Array<int>& dimensionsArray) {
209 
210 #ifdef HAVE_INTREPID_DEBUG
211 // srkenno@sandia.gov 6/12/10: changed unsigned int to int - this was causing a warning on compilers that
212 // signed & unsigned int's were being comparied.
213  for( int dim = 0; dim < dimensionsArray.size(); dim++) {
214  TEUCHOS_TEST_FOR_EXCEPTION( (0 > dimensionsArray[dim] ), std::invalid_argument,
215  ">>> ERROR (FieldContainer): One or more negative dimensions");
216  }
217 #endif
218 
219  // Copy dimensions and resize container storage to match them
220  dimensions_.assign(dimensionsArray.begin(),dimensionsArray.end());
221 
222  // Copy first 5 dimensions to optimize class for low rank containers
223  unsigned int theRank = dimensions_.size();
224  switch(theRank) {
225  case 1:
226  dim0_ = dimensions_[0];
227  dim1_ = 0;
228  dim2_ = 0;
229  dim3_ = 0;
230  dim4_ = 0;
231  break;
232 
233  case 2:
234  dim0_ = dimensions_[0];
235  dim1_ = dimensions_[1];
236  dim2_ = 0;
237  dim3_ = 0;
238  dim4_ = 0;
239  break;
240 
241  case 3:
242  dim0_ = dimensions_[0];
243  dim1_ = dimensions_[1];
244  dim2_ = dimensions_[2];
245  dim3_ = 0;
246  dim4_ = 0;
247  break;
248 
249  case 4:
250  dim0_ = dimensions_[0];
251  dim1_ = dimensions_[1];
252  dim2_ = dimensions_[2];
253  dim3_ = dimensions_[3];
254  dim4_ = 0;
255  break;
256 
257  case 5:
258  default:
259  dim0_ = dimensions_[0];
260  dim1_ = dimensions_[1];
261  dim2_ = dimensions_[2];
262  dim3_ = dimensions_[3];
263  dim4_ = dimensions_[4];
264  }
265 
266  // resize data array according to specified dimensions
267  data_.assign( this -> size(), (Scalar)0);
268  data_ptr_ = data_.begin();
269 
270 }
271 
272 
273 
274 template<class Scalar, int ArrayTypeId>
275 FieldContainer<Scalar, ArrayTypeId>::FieldContainer(const Teuchos::Array<int>& dimensionsArray,
276  const Teuchos::ArrayView<Scalar>& data) {
277 
278  // Copy all dimensions
279  dimensions_.assign(dimensionsArray.begin(),dimensionsArray.end());
280 
281  // Copy first 5 dimensions to optimize class for low rank containers
282  unsigned int theRank = dimensions_.size();
283  switch (theRank) {
284  case 0:
285  dim0_ = 0;
286  dim1_ = 0;
287  dim2_ = 0;
288  dim3_ = 0;
289  dim4_ = 0;
290  break;
291 
292  case 1:
293  dim0_ = dimensions_[0];
294  dim1_ = 0;
295  dim2_ = 0;
296  dim3_ = 0;
297  dim4_ = 0;
298  break;
299 
300  case 2:
301  dim0_ = dimensions_[0];
302  dim1_ = dimensions_[1];
303  dim2_ = 0;
304  dim3_ = 0;
305  dim4_ = 0;
306  break;
307 
308  case 3:
309  dim0_ = dimensions_[0];
310  dim1_ = dimensions_[1];
311  dim2_ = dimensions_[2];
312  dim3_ = 0;
313  dim4_ = 0;
314  break;
315 
316  case 4:
317  dim0_ = dimensions_[0];
318  dim1_ = dimensions_[1];
319  dim2_ = dimensions_[2];
320  dim3_ = dimensions_[3];
321  dim4_ = 0;
322  break;
323 
324  case 5:
325  default:
326  dim0_ = dimensions_[0];
327  dim1_ = dimensions_[1];
328  dim2_ = dimensions_[2];
329  dim3_ = dimensions_[3];
330  dim4_ = dimensions_[4];
331  }
332 
333  // Validate input: size of data array must match container size specified by its dimensions
334 #ifdef HAVE_INTREPID_DEBUG
335  TEUCHOS_TEST_FOR_EXCEPTION( ( (int)data.size() != this -> size() ),
336  std::invalid_argument,
337  ">>> ERROR (FieldContainer): Size of input data does not match size of this container.");
338 #endif
339 
340  // Deep-copy ArrayView data.
341  data_.assign(data.begin(),data.end());
342  data_ptr_ = data_.begin();
343 }
344 
345 
346 
347 template<class Scalar, int ArrayTypeId>
348 FieldContainer<Scalar, ArrayTypeId>::FieldContainer(const Teuchos::Array<int>& dimensionsArray,
349  const Teuchos::ArrayRCP<Scalar>& data) {
350 
351  // Copy all dimensions
352  dimensions_.assign(dimensionsArray.begin(),dimensionsArray.end());
353 
354  // Copy first 5 dimensions to optimize class for low rank containers
355  unsigned int theRank = dimensions_.size();
356  switch(theRank) {
357  case 0:
358  dim0_ = 0;
359  dim1_ = 0;
360  dim2_ = 0;
361  dim3_ = 0;
362  dim4_ = 0;
363  break;
364 
365  case 1:
366  dim0_ = dimensions_[0];
367  dim1_ = 0;
368  dim2_ = 0;
369  dim3_ = 0;
370  dim4_ = 0;
371  break;
372 
373  case 2:
374  dim0_ = dimensions_[0];
375  dim1_ = dimensions_[1];
376  dim2_ = 0;
377  dim3_ = 0;
378  dim4_ = 0;
379  break;
380 
381  case 3:
382  dim0_ = dimensions_[0];
383  dim1_ = dimensions_[1];
384  dim2_ = dimensions_[2];
385  dim3_ = 0;
386  dim4_ = 0;
387  break;
388 
389  case 4:
390  dim0_ = dimensions_[0];
391  dim1_ = dimensions_[1];
392  dim2_ = dimensions_[2];
393  dim3_ = dimensions_[3];
394  dim4_ = 0;
395  break;
396 
397  case 5:
398  default:
399  dim0_ = dimensions_[0];
400  dim1_ = dimensions_[1];
401  dim2_ = dimensions_[2];
402  dim3_ = dimensions_[3];
403  dim4_ = dimensions_[4];
404  }
405 
406  // Validate input: size of data array must match container size specified by its dimensions
407 #ifdef HAVE_INTREPID_DEBUG
408  TEUCHOS_TEST_FOR_EXCEPTION( ( (int)data.size() != this -> size() ),
409  std::invalid_argument,
410  ">>> ERROR (FieldContainer): Size of input data does not match size of this container.");
411 #endif
412 
413  // Shallow-copy ArrayRCP data.
414  data_ = data;
415  data_ptr_ = data_.begin();
416 }
417 
418 
419 
420 template<class Scalar, int ArrayTypeId>
421 FieldContainer<Scalar, ArrayTypeId>::FieldContainer(const Teuchos::Array<int>& dimensionsArray,
422  Scalar* data,
423  const bool deep_copy,
424  const bool owns_mem) {
425 
426  // Copy all dimensions
427  dimensions_.assign(dimensionsArray.begin(),dimensionsArray.end());
428 
429  // Copy first 5 dimensions to optimize class for low rank containers
430  unsigned int theRank = dimensions_.size();
431  switch (theRank) {
432  case 0:
433  dim0_ = 0;
434  dim1_ = 0;
435  dim2_ = 0;
436  dim3_ = 0;
437  dim4_ = 0;
438  break;
439 
440  case 1:
441  dim0_ = dimensions_[0];
442  dim1_ = 0;
443  dim2_ = 0;
444  dim3_ = 0;
445  dim4_ = 0;
446  break;
447 
448  case 2:
449  dim0_ = dimensions_[0];
450  dim1_ = dimensions_[1];
451  dim2_ = 0;
452  dim3_ = 0;
453  dim4_ = 0;
454  break;
455 
456  case 3:
457  dim0_ = dimensions_[0];
458  dim1_ = dimensions_[1];
459  dim2_ = dimensions_[2];
460  dim3_ = 0;
461  dim4_ = 0;
462  break;
463 
464  case 4:
465  dim0_ = dimensions_[0];
466  dim1_ = dimensions_[1];
467  dim2_ = dimensions_[2];
468  dim3_ = dimensions_[3];
469  dim4_ = 0;
470  break;
471 
472  case 5:
473  default:
474  dim0_ = dimensions_[0];
475  dim1_ = dimensions_[1];
476  dim2_ = dimensions_[2];
477  dim3_ = dimensions_[3];
478  dim4_ = dimensions_[4];
479  }
480 
481 
482  if (deep_copy) {
483  Teuchos::ArrayRCP<Scalar> arrayrcp = Teuchos::arcp<Scalar>(data, 0, this -> size(), false);
484  data_.deepCopy(arrayrcp());
485  data_ptr_ = data_.begin();
486  }
487  else {
488  data_ = Teuchos::arcp<Scalar>(data, 0, this -> size(), owns_mem);
489  data_ptr_ = data_.begin();
490  }
491 }
492 
493 
494 
495 template<class Scalar, int ArrayTypeId>
496 FieldContainer<Scalar, ArrayTypeId>::FieldContainer(const shards::Array<Scalar,shards::NaturalOrder>& data,
497  const bool deep_copy,
498  const bool owns_mem) {
499 
500  // Copy all dimensions
501  dimensions_.resize(data.rank());
502 
503  // Copy first 5 dimensions to optimize class for low rank containers
504  unsigned int theRank = dimensions_.size();
505  switch(theRank) {
506  case 1:
507  dimensions_[0] = data.dimension(0);
508  dim0_ = dimensions_[0];
509  dim1_ = 0;
510  dim2_ = 0;
511  dim3_ = 0;
512  dim4_ = 0;
513  break;
514 
515  case 2:
516  dimensions_[0] = data.dimension(0);
517  dimensions_[1] = data.dimension(1);
518  dim0_ = dimensions_[0];
519  dim1_ = dimensions_[1];
520  dim2_ = 0;
521  dim3_ = 0;
522  dim4_ = 0;
523  break;
524 
525  case 3:
526  dimensions_[0] = data.dimension(0);
527  dimensions_[1] = data.dimension(1);
528  dimensions_[2] = data.dimension(2);
529  dim0_ = dimensions_[0];
530  dim1_ = dimensions_[1];
531  dim2_ = dimensions_[2];
532  dim3_ = 0;
533  dim4_ = 0;
534  break;
535 
536  case 4:
537  dimensions_[0] = data.dimension(0);
538  dimensions_[1] = data.dimension(1);
539  dimensions_[2] = data.dimension(2);
540  dimensions_[3] = data.dimension(3);
541  dim0_ = dimensions_[0];
542  dim1_ = dimensions_[1];
543  dim2_ = dimensions_[2];
544  dim3_ = dimensions_[3];
545  dim4_ = 0;
546  break;
547 
548  case 5:
549  dimensions_[0] = data.dimension(0);
550  dimensions_[1] = data.dimension(1);
551  dimensions_[2] = data.dimension(2);
552  dimensions_[3] = data.dimension(3);
553  dimensions_[4] = data.dimension(4);
554  dim0_ = dimensions_[0];
555  dim1_ = dimensions_[1];
556  dim2_ = dimensions_[2];
557  dim3_ = dimensions_[3];
558  dim4_ = dimensions_[4];
559  break;
560 
561  default:
562  for (int i=0; i<data.rank(); i++) {
563  dimensions_[i] = data.dimension(i);
564  }
565  dim0_ = dimensions_[0];
566  dim1_ = dimensions_[1];
567  dim2_ = dimensions_[2];
568  dim3_ = dimensions_[3];
569  dim4_ = dimensions_[4];
570  }
571 
572 
573  if (deep_copy) {
574  Teuchos::ArrayRCP<Scalar> arrayrcp = Teuchos::arcp<Scalar>(data.contiguous_data(), 0, this -> size(), false);
575  data_.deepCopy(arrayrcp());
576  data_ptr_ = data_.begin();
577  }
578  else {
579  data_ = Teuchos::arcp<Scalar>(data.contiguous_data(), 0, this -> size(), owns_mem);
580  data_ptr_ = data_.begin();
581  }
582 }
583 
584 
585 
586 //--------------------------------------------------------------------------------------------//
587 // //
588 // Access methods of FieldContainer class //
589 // //
590 //--------------------------------------------------------------------------------------------//
591 
592 
593 template<class Scalar, int ArrayTypeId>
595  return dimensions_.size();
596 }
597 
598 
599 
600 template<class Scalar, int ArrayTypeId>
602  // Important! This method is used by constructors to find out what is the needed size of data_
603  // based on the specified dimensions. Therefore, it cannot be implmented by returning data_.size
604  // and must be able to compute the size of the container based only on its specified dimensions
605 
606  // Size equals product of all dimensions stored in dimensions_
607  const int theRank = dimensions_.size ();
608 
609  // If container has no dimensions its size is zero
610  if (theRank == 0) {
611  return 0;
612  }
613  else {
614  int theSize = dim0_;
615 
616  // Compute size directly to optimize method for low rank (<=5) containers
617  switch(theRank) {
618  case 5:
619  theSize *= dim1_*dim2_*dim3_*dim4_;
620  break;
621 
622  case 4:
623  theSize *= dim1_*dim2_*dim3_;
624  break;
625 
626  case 3:
627  theSize *= dim1_*dim2_;
628  break;
629 
630  case 2:
631  theSize *= dim1_;
632  break;
633 
634  case 1:
635  break;
636 
637  // Compute size for containers with ranks hihger than 5
638  default:
639  for (int r = 1; r < theRank; ++r) {
640  theSize *= dimensions_[r];
641  }
642  }
643  return theSize;
644  }
645 }
646 
647 
648 
649 template<class Scalar, int ArrayTypeId>
650 template<class Vector>
651 inline void FieldContainer<Scalar, ArrayTypeId>::dimensions(Vector& dimVec) const {
652  dimVec.assign(dimensions_.begin(),dimensions_.end());
653 }
654 
655 
656 
657 template<class Scalar, int ArrayTypeId>
658 inline int FieldContainer<Scalar, ArrayTypeId>::dimension(const int whichDim) const {
659 #ifdef HAVE_INTREPID_DEBUG
660  TEUCHOS_TEST_FOR_EXCEPTION( (0 > whichDim), std::invalid_argument,
661  ">>> ERROR (FieldContainer): dimension order cannot be negative");
662  TEUCHOS_TEST_FOR_EXCEPTION( (whichDim >= this -> rank() ), std::invalid_argument,
663  ">>> ERROR (FieldContainer): dimension order cannot exceed rank of the container");
664 #endif
665  return dimensions_[whichDim];
666 }
667 
668 
669 
670 template<class Scalar, int ArrayTypeId>
672 #ifdef HAVE_INTREPID_DEBUG
673  TEUCHOS_TEST_FOR_EXCEPTION( ( this -> rank() != 1), std::invalid_argument,
674  ">>> ERROR (FieldContainer): Number of indices does not match rank of the container.");
675  TEUCHOS_TEST_FOR_EXCEPTION( ( (i0 < 0) || ( i0 >= dim0_) ), std::invalid_argument,
676  ">>> ERROR (FieldContainer): index is out of range.");
677 #endif
678  return i0;
679 }
680 
681 
682 
683 template<class Scalar, int ArrayTypeId>
685  const int i1) const {
686 #ifdef HAVE_INTREPID_DEBUG
687  TEUCHOS_TEST_FOR_EXCEPTION( ( this -> rank() != 2), std::invalid_argument,
688  ">>> ERROR (FieldContainer): Number of indices does not match rank of the container.");
689  TEUCHOS_TEST_FOR_EXCEPTION( ( (i0 < 0) || ( i0 >= dim0_) ), std::invalid_argument,
690  ">>> ERROR (FieldContainer): 1st index is out of range.");
691  TEUCHOS_TEST_FOR_EXCEPTION( ( (i1 < 0) || (i1 >= dim1_) ), std::invalid_argument,
692  ">>> ERROR (FieldContainer): 2nd index is out of range.");
693 #endif
694  return i0*dim1_ + i1;
695 }
696 
697 
698 
699 template<class Scalar, int ArrayTypeId>
701  const int i1,
702  const int i2) const {
703 #ifdef HAVE_INTREPID_DEBUG
704  TEUCHOS_TEST_FOR_EXCEPTION( ( this -> rank() != 3), std::invalid_argument,
705  ">>> ERROR (FieldContainer): Number of indices does not match rank of the container.");
706  TEUCHOS_TEST_FOR_EXCEPTION( ( (i0 < 0) || ( i0 >= dim0_) ), std::invalid_argument,
707  ">>> ERROR (FieldContainer): 1st index is out of range.");
708  TEUCHOS_TEST_FOR_EXCEPTION( ( (i1 < 0) || (i1 >= dim1_) ), std::invalid_argument,
709  ">>> ERROR (FieldContainer): 2nd index is out of range.");
710  TEUCHOS_TEST_FOR_EXCEPTION( ( (i2 < 0) || (i2 >= dim2_) ), std::invalid_argument,
711  ">>> ERROR (FieldContainer): 3rd index is out of range.");
712 #endif
713  return (i0*dim1_ + i1)*dim2_ + i2;
714 }
715 
716 
717 
718 template<class Scalar, int ArrayTypeId>
720  const int i1,
721  const int i2,
722  const int i3) const {
723 #ifdef HAVE_INTREPID_DEBUG
724  TEUCHOS_TEST_FOR_EXCEPTION( ( this -> rank() != 4), std::invalid_argument,
725  ">>> ERROR (FieldContainer): Number of indices does not match rank of the container.");
726  TEUCHOS_TEST_FOR_EXCEPTION( ( (i0 < 0) || ( i0 >= dim0_) ), std::invalid_argument,
727  ">>> ERROR (FieldContainer): 1st index is out of range.");
728  TEUCHOS_TEST_FOR_EXCEPTION( ( (i1 < 0) || (i1 >= dim1_) ), std::invalid_argument,
729  ">>> ERROR (FieldContainer): 2nd index is out of range.");
730  TEUCHOS_TEST_FOR_EXCEPTION( ( (i2 < 0) || (i2 >= dim2_) ), std::invalid_argument,
731  ">>> ERROR (FieldContainer): 3rd index is out of range.");
732  TEUCHOS_TEST_FOR_EXCEPTION( ( (i3 < 0) || (i3 >= dim3_) ), std::invalid_argument,
733  ">>> ERROR (FieldContainer): 4th index is out of range.");
734 #endif
735  return ( (i0*dim1_ + i1 )*dim2_ + i2 )*dim3_ + i3;
736 }
737 
738 
739 
740 template<class Scalar, int ArrayTypeId>
742  const int i1,
743  const int i2,
744  const int i3,
745  const int i4) const {
746 #ifdef HAVE_INTREPID_DEBUG
747  TEUCHOS_TEST_FOR_EXCEPTION( ( this -> rank() != 5), std::invalid_argument,
748  ">>> ERROR (FieldContainer): Number of indices does not match rank of the container.");
749  TEUCHOS_TEST_FOR_EXCEPTION( ( (i0 < 0) || ( i0 >= dim0_) ), std::invalid_argument,
750  ">>> ERROR (FieldContainer): 1st index is out of range.");
751  TEUCHOS_TEST_FOR_EXCEPTION( ( (i1 < 0) || (i1 >= dim1_) ), std::invalid_argument,
752  ">>> ERROR (FieldContainer): 2nd index is out of range.");
753  TEUCHOS_TEST_FOR_EXCEPTION( ( (i2 < 0) || (i2 >= dim2_) ), std::invalid_argument,
754  ">>> ERROR (FieldContainer): 3rd index is out of range.");
755  TEUCHOS_TEST_FOR_EXCEPTION( ( (i3 < 0) || (i3 >= dim3_) ), std::invalid_argument,
756  ">>> ERROR (FieldContainer): 4th index is out of range.");
757  TEUCHOS_TEST_FOR_EXCEPTION( ( (i4 < 0) || (i4 >= dim4_) ), std::invalid_argument,
758  ">>> ERROR (FieldContainer): 5th index is out of range.");
759 #endif
760  return ( ( (i0*dim1_ + i1 )*dim2_ + i2 )*dim3_ + i3 )*dim4_ + i4;
761 }
762 
763 
764 
765 
766 template<class Scalar, int ArrayTypeId>
767 int FieldContainer<Scalar, ArrayTypeId>::getEnumeration(const Teuchos::Array<int>& multiIndex) const {
768 
769 #ifdef HAVE_INTREPID_DEBUG
770  // Check if number of multi-indices matches rank of the FieldContainer object
771  TEUCHOS_TEST_FOR_EXCEPTION( ( multiIndex.size() != dimensions_.size() ),
772  std::invalid_argument,
773  ">>> ERROR (FieldContainer): Number of multi-indices does not match rank of container.");
774  TEUCHOS_TEST_FOR_EXCEPTION( ( ( multiIndex[0] < 0) || ( multiIndex[0] >= dim0_) ),
775  std::invalid_argument,
776  ">>> ERROR (FieldContainer): 1st index is out of range.");
777 #endif
778 
779  int theRank = dimensions_.size();
780  int address = 0;
781  switch (theRank) {
782 
783  // Optimize enumeration computation for low rank (<= 5) containers
784  case 5:
785 #ifdef HAVE_INTREPID_DEBUG
786  TEUCHOS_TEST_FOR_EXCEPTION( ( (multiIndex[4] < 0) || (multiIndex[4] >= dim4_) ),
787  std::invalid_argument,
788  ">>> ERROR (FieldContainer): 5th index is out of range.");
789  TEUCHOS_TEST_FOR_EXCEPTION( ( (multiIndex[3] < 0) || (multiIndex[3] >= dim3_) ),
790  std::invalid_argument,
791  ">>> ERROR (FieldContainer): 4th index is out of range.");
792  TEUCHOS_TEST_FOR_EXCEPTION( ( (multiIndex[2] < 0) || (multiIndex[2] >= dim2_) ),
793  std::invalid_argument,
794  ">>> ERROR (FieldContainer): 3rd index is out of range.");
795  TEUCHOS_TEST_FOR_EXCEPTION( ( (multiIndex[1] < 0) || (multiIndex[1] >= dim1_) ),
796  std::invalid_argument,
797  ">>> ERROR (FieldContainer): 2nd index is out of range.");
798 #endif
799  address = (((multiIndex[0]*dim1_ + multiIndex[1])*dim2_ + multiIndex[2])*dim3_ + multiIndex[3])*dim4_ + multiIndex[4];
800  break;
801 
802  case 4:
803 #ifdef HAVE_INTREPID_DEBUG
804  TEUCHOS_TEST_FOR_EXCEPTION( ( (multiIndex[3] < 0) || (multiIndex[3] >= dim3_) ),
805  std::invalid_argument,
806  ">>> ERROR (FieldContainer): 4th index is out of range.");
807  TEUCHOS_TEST_FOR_EXCEPTION( ( (multiIndex[2] < 0) || (multiIndex[2] >= dim2_) ),
808  std::invalid_argument,
809  ">>> ERROR (FieldContainer): 3rd index is out of range.");
810  TEUCHOS_TEST_FOR_EXCEPTION( ( (multiIndex[1] < 0) || (multiIndex[1] >= dim1_) ),
811  std::invalid_argument,
812  ">>> ERROR (FieldContainer): 2nd index is out of range.");
813 #endif
814  address = ((multiIndex[0]*dim1_ + multiIndex[1])*dim2_ + multiIndex[2])*dim3_ + multiIndex[3];
815  break;
816 
817  case 3:
818 #ifdef HAVE_INTREPID_DEBUG
819  TEUCHOS_TEST_FOR_EXCEPTION( ( (multiIndex[2] < 0) || (multiIndex[2] >= dim2_) ),
820  std::invalid_argument,
821  ">>> ERROR (FieldContainer): 3rd index is out of range.");
822  TEUCHOS_TEST_FOR_EXCEPTION( ( (multiIndex[1] < 0) || (multiIndex[1] >= dim1_) ),
823  std::invalid_argument,
824  ">>> ERROR (FieldContainer): 2nd index is out of range.");
825 #endif
826  address = (multiIndex[0]*dim1_ + multiIndex[1])*dim2_ + multiIndex[2];
827  break;
828 
829  case 2:
830 #ifdef HAVE_INTREPID_DEBUG
831  TEUCHOS_TEST_FOR_EXCEPTION( ( (multiIndex[1] < 0) || (multiIndex[1] >= dim1_) ),
832  std::invalid_argument,
833  ">>> ERROR (FieldContainer): 2nd index is out of range.");
834 #endif
835  address = multiIndex[0]*dim1_ + multiIndex[1];
836  break;
837 
838  case 1:
839  address = multiIndex[0];
840  break;
841 
842  default:
843 
844  // Arbitrary rank: compute address using Horner's nested scheme: intialize address to 0th index value
845  address = multiIndex[0];
846  for (int r = 0; r < theRank - 1; r++){
847 #ifdef HAVE_INTREPID_DEBUG
848  TEUCHOS_TEST_FOR_EXCEPTION( ( (multiIndex[r+1] < 0) || (multiIndex[r+1] >= dimensions_[r+1]) ),
849  std::invalid_argument,
850  ">>> ERROR (FieldContainer): Multi-index component out of range.");
851 #endif
852  // Add increment
853  address = address*dimensions_[r+1] + multiIndex[r+1];
854  }
855  } // end switch(theRank)
856 
857  return address;
858 }
859 
860 
861 
862 template<class Scalar, int ArrayTypeId>
864  const int valueEnum) const
865 {
866 #ifdef HAVE_INTREPID_DEBUG
867  TEUCHOS_TEST_FOR_EXCEPTION( ( this -> rank() != 1), std::invalid_argument,
868  ">>> ERROR (FieldContainer): Number of indices does not match rank of the container.");
869  TEUCHOS_TEST_FOR_EXCEPTION( ( (valueEnum < 0) || (valueEnum >= (int)data_.size()) ),
870  std::invalid_argument,
871  ">>> ERROR (FieldContainer): Value enumeration is out of range.");
872 #endif
873  i0 = valueEnum;
874 }
875 
876 
877 
878 template<class Scalar, int ArrayTypeId>
880  int & i1,
881  const int valueEnum) const
882 {
883 #ifdef HAVE_INTREPID_DEBUG
884  TEUCHOS_TEST_FOR_EXCEPTION( ( this -> rank() != 2), std::invalid_argument,
885  ">>> ERROR (FieldContainer): Number of indices does not match rank of the container.");
886  TEUCHOS_TEST_FOR_EXCEPTION( ( (valueEnum < 0) || (valueEnum >= (int)data_.size()) ),
887  std::invalid_argument,
888  ">>> ERROR (FieldContainer): Value enumeration is out of range.");
889 #endif
890 
891  i0 = valueEnum/dim1_;
892  i1 = valueEnum - i0*dim1_;
893 }
894 
895 
896 
897 template<class Scalar, int ArrayTypeId>
899  int & i1,
900  int & i2,
901  const int valueEnum) const
902 {
903 #ifdef HAVE_INTREPID_DEBUG
904  TEUCHOS_TEST_FOR_EXCEPTION( ( this -> rank() != 3), std::invalid_argument,
905  ">>> ERROR (FieldContainer): Number of indices does not match rank of the container.");
906  TEUCHOS_TEST_FOR_EXCEPTION( ( (valueEnum < 0) || (valueEnum >= (int)data_.size()) ),
907  std::invalid_argument,
908  ">>> ERROR (FieldContainer): Value enumeration is out of range.");
909 #endif
910  int tempDim = dim1_*dim2_;
911  int tempEnu = valueEnum;
912  i0 = tempEnu/tempDim;
913 
914  tempEnu -= i0*tempDim;
915  tempDim /= dim1_;
916  i1 = tempEnu/tempDim;
917 
918  tempEnu -= i1*tempDim;
919  tempDim /= dim2_;
920  i2 = tempEnu/tempDim;
921 }
922 
923 
924 
925 template<class Scalar, int ArrayTypeId>
927  int & i1,
928  int & i2,
929  int & i3,
930  const int valueEnum) const
931 {
932 #ifdef HAVE_INTREPID_DEBUG
933  TEUCHOS_TEST_FOR_EXCEPTION( ( this -> rank() != 4), std::invalid_argument,
934  ">>> ERROR (FieldContainer): Number of indices does not match rank of the container.");
935  TEUCHOS_TEST_FOR_EXCEPTION( ( (valueEnum < 0) || (valueEnum >= (int)data_.size()) ),
936  std::invalid_argument,
937  ">>> ERROR (FieldContainer): Value enumeration is out of range.");
938 #endif
939  int tempDim = dim1_*dim2_*dim3_;
940  int tempEnu = valueEnum;
941  i0 = tempEnu/tempDim;
942 
943  tempEnu -= i0*tempDim;
944  tempDim /= dim1_;
945  i1 = tempEnu/tempDim;
946 
947  tempEnu -= i1*tempDim;
948  tempDim /= dim2_;
949  i2 = tempEnu/tempDim;
950 
951  tempEnu -= i2*tempDim;
952  tempDim /= dim3_;
953  i3 = tempEnu/tempDim;
954 }
955 
956 
957 
958 
959 template<class Scalar, int ArrayTypeId>
961  int & i1,
962  int & i2,
963  int & i3,
964  int & i4,
965  const int valueEnum) const
966 {
967 #ifdef HAVE_INTREPID_DEBUG
968  TEUCHOS_TEST_FOR_EXCEPTION( ( this -> rank() != 5), std::invalid_argument,
969  ">>> ERROR (FieldContainer): Number of indices does not match rank of the container.");
970  TEUCHOS_TEST_FOR_EXCEPTION( ( (valueEnum < 0) || (valueEnum >= (int)data_.size()) ),
971  std::invalid_argument,
972  ">>> ERROR (FieldContainer): Value enumeration is out of range.");
973 #endif
974  int tempDim = dim1_*dim2_*dim3_*dim4_;
975  int tempEnu = valueEnum;
976  i0 = tempEnu/tempDim;
977 
978  tempEnu -= i0*tempDim;
979  tempDim /= dim1_;
980  i1 = tempEnu/tempDim;
981 
982  tempEnu -= i1*tempDim;
983  tempDim /= dim2_;
984  i2 = tempEnu/tempDim;
985 
986  tempEnu -= i2*tempDim;
987  tempDim /= dim3_;
988  i3 = tempEnu/tempDim;
989 
990  tempEnu -= i3*tempDim;
991  tempDim /= dim4_;
992  i4 = tempEnu/tempDim;
993 }
994 
995 
996 
997 template<class Scalar, int ArrayTypeId>
998 template<class Vector>
1000  const int valueEnum) const
1001 {
1002 
1003  // Verify address is in the admissible range for this FieldContainer
1004 #ifdef HAVE_INTREPID_DEBUG
1005  TEUCHOS_TEST_FOR_EXCEPTION( ( (valueEnum < 0) || (valueEnum >= (int)data_.size()) ),
1006  std::invalid_argument,
1007  ">>> ERROR (FieldContainer): Value enumeration is out of range.");
1008 #endif
1009 
1010  // make sure multiIndex has the right size to hold all multi-indices
1011  const int theRank = dimensions_.size ();
1012  multiIndex.resize (theRank);
1013 
1014  // Initializations
1015  int temp_enum = valueEnum;
1016  int temp_range = 1;
1017 
1018  // Compute product of all but the first upper bound
1019  for (int r = 1; r < theRank; ++r) {
1020  temp_range *= dimensions_[r];
1021  }
1022 
1023  // Index 0 is computed first using integer division
1024  multiIndex[0] = temp_enum/temp_range;
1025 
1026  // Indices 1 to (theRank - 2) are computed next; will be skipped if
1027  // theRank <= 2
1028  for (int r = 1; r < theRank - 1; ++r) {
1029  temp_enum -= multiIndex[r-1]*temp_range;
1030  temp_range /= dimensions_[r];
1031  multiIndex[r] = temp_enum/temp_range;
1032  }
1033 
1034  // Index (theRank - 1) is computed last, skip if theRank = 1 and
1035  // keep if theRank = 2
1036  if (theRank > 1) {
1037  multiIndex[theRank - 1] = temp_enum - multiIndex[theRank - 2] * temp_range;
1038  }
1039 }
1040 
1041 //--------------------------------------------------------------------------------------------//
1042 // //
1043 // Methods to shape (resize) a field container //
1044 // //
1045 //--------------------------------------------------------------------------------------------//
1046 
1047 template<class Scalar, int ArrayTypeId>
1049  dimensions_.resize(0);
1050 
1051  // Reset first five dimensions:
1052  dim0_ = 0;
1053  dim1_ = 0;
1054  dim2_ = 0;
1055  dim3_ = 0;
1056  dim4_ = 0;
1057 
1058  // Clears data array and sets to zero length
1059  data_.clear();
1060  data_ptr_ = data_.begin();
1061 }
1062 
1063 
1064 
1065 template<class Scalar, int ArrayTypeId>
1066 void FieldContainer<Scalar, ArrayTypeId>::resize(const Teuchos::Array<int>& newDimensions) {
1067 
1068  // First handle the trivial case of zero dimensions
1069  if( newDimensions.size() == 0) {
1070  dimensions_.resize(0);
1071  dim0_ = 0;
1072  dim1_ = 0;
1073  dim2_ = 0;
1074  dim3_ = 0;
1075  dim4_ = 0;
1076  data_.resize(0);
1077  data_ptr_ = data_.begin();
1078  }
1079  else {
1080 
1081  // Copy upper index bounds and resize container storage to match new upper bounds.
1082  dimensions_.assign(newDimensions.begin(),newDimensions.end());
1083 
1084  // Copy first 5 dimensions for faster access
1085  unsigned int theRank = dimensions_.size();
1086  switch (theRank) {
1087  case 1:
1088  dim0_ = dimensions_[0];
1089  dim1_ = 0;
1090  dim2_ = 0;
1091  dim3_ = 0;
1092  dim4_ = 0;
1093  break;
1094 
1095  case 2:
1096  dim0_ = dimensions_[0];
1097  dim1_ = dimensions_[1];
1098  dim2_ = 0;
1099  dim3_ = 0;
1100  dim4_ = 0;
1101  break;
1102 
1103  case 3:
1104  dim0_ = dimensions_[0];
1105  dim1_ = dimensions_[1];
1106  dim2_ = dimensions_[2];
1107  dim3_ = 0;
1108  dim4_ = 0;
1109  break;
1110 
1111  case 4:
1112  dim0_ = dimensions_[0];
1113  dim1_ = dimensions_[1];
1114  dim2_ = dimensions_[2];
1115  dim3_ = dimensions_[3];
1116  dim4_ = 0;
1117  break;
1118 
1119  case 5:
1120  default:
1121  dim0_ = dimensions_[0];
1122  dim1_ = dimensions_[1];
1123  dim2_ = dimensions_[2];
1124  dim3_ = dimensions_[3];
1125  dim4_ = dimensions_[4];
1126  }
1127 
1128  // Resize data array
1129  data_.resize(this -> size());
1130  data_ptr_ = data_.begin();
1131  }
1132 }
1133 
1134 
1135 
1136 template<class Scalar, int ArrayTypeId>
1137 inline void FieldContainer<Scalar, ArrayTypeId>::resize(const int dim0) {
1138  dim0_ = dim0;
1139  dim1_ = 0;
1140  dim2_ = 0;
1141  dim3_ = 0;
1142  dim4_ = 0;
1143  dimensions_.resize(1);
1144  dimensions_[0] = dim0_;
1145  data_.resize(dim0_);
1146  data_ptr_ = data_.begin();
1147 }
1148 
1149 
1150 
1151 template<class Scalar, int ArrayTypeId>
1153  const int dim1) {
1154  dim0_ = dim0;
1155  dim1_ = dim1;
1156  dim2_ = 0;
1157  dim3_ = 0;
1158  dim4_ = 0;
1159  dimensions_.resize(2);
1160  dimensions_[0] = dim0_;
1161  dimensions_[1] = dim1_;
1162  data_.resize(dim0_*dim1_);
1163  data_ptr_ = data_.begin();
1164 }
1165 
1166 
1167 
1168 template<class Scalar, int ArrayTypeId>
1170  const int dim1,
1171  const int dim2) {
1172  dim0_ = dim0;
1173  dim1_ = dim1;
1174  dim2_ = dim2;
1175  dim3_ = 0;
1176  dim4_ = 0;
1177  dimensions_.resize(3);
1178  dimensions_[0] = dim0_;
1179  dimensions_[1] = dim1_;
1180  dimensions_[2] = dim2_;
1181  data_.resize(dim0_*dim1_*dim2_);
1182  data_ptr_ = data_.begin();
1183 }
1184 
1185 
1186 
1187 template<class Scalar, int ArrayTypeId>
1189  const int dim1,
1190  const int dim2,
1191  const int dim3) {
1192  dim0_ = dim0;
1193  dim1_ = dim1;
1194  dim2_ = dim2;
1195  dim3_ = dim3;
1196  dim4_ = 0;
1197  dimensions_.resize(4);
1198  dimensions_[0] = dim0_;
1199  dimensions_[1] = dim1_;
1200  dimensions_[2] = dim2_;
1201  dimensions_[3] = dim3_;
1202  data_.resize(dim0_*dim1_*dim2_*dim3_);
1203  data_ptr_ = data_.begin();
1204 }
1205 
1206 
1207 
1208 template<class Scalar, int ArrayTypeId>
1210  const int dim1,
1211  const int dim2,
1212  const int dim3,
1213  const int dim4) {
1214  dim0_ = dim0;
1215  dim1_ = dim1;
1216  dim2_ = dim2;
1217  dim3_ = dim3;
1218  dim4_ = dim4;
1219  dimensions_.resize(5);
1220  dimensions_[0] = dim0_;
1221  dimensions_[1] = dim1_;
1222  dimensions_[2] = dim2_;
1223  dimensions_[3] = dim3_;
1224  dimensions_[4] = dim4_;
1225  data_.resize(dim0_*dim1_*dim2_*dim3_*dim4_);
1226  data_ptr_ = data_.begin();
1227 }
1228 
1229 
1230 
1231 template<class Scalar, int ArrayTypeId>
1233 
1234  // Copy dimensions from the specified container
1235  anotherContainer.dimensions(dimensions_);
1236  int newRank = dimensions_.size();
1237 
1238  // Copy first 5 dimensions for faster access
1239  switch(newRank) {
1240  case 1:
1241  dim0_ = dimensions_[0];
1242  dim1_ = 0;
1243  dim2_ = 0;
1244  dim3_ = 0;
1245  dim4_ = 0;
1246  break;
1247 
1248  case 2:
1249  dim0_ = dimensions_[0];
1250  dim1_ = dimensions_[1];
1251  dim2_ = 0;
1252  dim3_ = 0;
1253  dim4_ = 0;
1254  break;
1255 
1256  case 3:
1257  dim0_ = dimensions_[0];
1258  dim1_ = dimensions_[1];
1259  dim2_ = dimensions_[2];
1260  dim3_ = 0;
1261  dim4_ = 0;
1262  break;
1263 
1264  case 4:
1265  dim0_ = dimensions_[0];
1266  dim1_ = dimensions_[1];
1267  dim2_ = dimensions_[2];
1268  dim3_ = dimensions_[3];
1269  dim4_ = 0;
1270  break;
1271 
1272  case 5:
1273  default:
1274  dim0_ = dimensions_[0];
1275  dim1_ = dimensions_[1];
1276  dim2_ = dimensions_[2];
1277  dim3_ = dimensions_[3];
1278  dim4_ = dimensions_[4];
1279  }
1280 
1281  // Resize data array
1282  data_.resize(this->size());
1283  data_ptr_ = data_.begin();
1284 }
1285 
1286 
1287 template<class Scalar, int ArrayTypeId>
1289  const int numFields,
1290  const EFunctionSpace spaceType,
1291  const EOperator operatorType,
1292  const int spaceDim) {
1293  // Validate input
1294 #ifdef HAVE_INTREPID_DEBUG
1295  TEUCHOS_TEST_FOR_EXCEPTION( ( numPoints < 0),
1296  std::invalid_argument,
1297  ">>> ERROR (FieldContainer): Number of points cannot be negative!");
1298  TEUCHOS_TEST_FOR_EXCEPTION( ( numFields < 0),
1299  std::invalid_argument,
1300  ">>> ERROR (FieldContainer): Number of fields cannot be negative!");
1301  TEUCHOS_TEST_FOR_EXCEPTION( !( (1 <= spaceDim ) && ( spaceDim <= 3 ) ),
1302  std::invalid_argument,
1303  ">>> ERROR (FieldContainer): Invalid space dimension.");
1304 #endif
1305 
1306  // Find out field and operator ranks
1307  const int fieldRank = getFieldRank(spaceType);
1308  const int operatorRank = getOperatorRank(spaceType,operatorType,spaceDim);
1309 
1310  // Compute rank of the container = 1(numPoints) + 1(numFields) + fieldRank + operatorRank
1311  const int theRank = 1 + 1 + fieldRank + operatorRank;
1312 
1313  // Define temp array for the dimensions
1314  Teuchos::Array<int> newDimensions (theRank);
1315 
1316  // Dimensions 0 and 1 are number of points and number of fields, resp.
1317  newDimensions[0] = numPoints;
1318  newDimensions[1] = numFields;
1319 
1320  // The rest of the dimensions depend on whether we had VALUE, GRAD (D1), CURL, DIV or Dk, k>1
1321  switch (operatorType) {
1322 
1323  case OPERATOR_VALUE:
1324  case OPERATOR_GRAD:
1325  case OPERATOR_D1:
1326  case OPERATOR_CURL:
1327  case OPERATOR_DIV:
1328 
1329  // For these operators all dimensions from 2 to 2 + fieldRank + OperatorRank are bounded by spaceDim
1330  for (int i = 0; i < fieldRank + operatorRank; ++i) {
1331  newDimensions[2 + i] = spaceDim;
1332  }
1333  break;
1334 
1335  case OPERATOR_D2:
1336  case OPERATOR_D3:
1337  case OPERATOR_D4:
1338  case OPERATOR_D5:
1339  case OPERATOR_D6:
1340  case OPERATOR_D7:
1341  case OPERATOR_D8:
1342  case OPERATOR_D9:
1343  case OPERATOR_D10:
1344 
1345  // All dimensions from 2 to 2 + fieldRank, if any, are bounded by spaceDim
1346  for(int i = 0; i < fieldRank; i++){
1347  newDimensions[2 + i] = spaceDim;
1348  }
1349 
1350  // We know that for Dk operatorRank = 1 and so there's just one more dimension left
1351  // given by the cardinality of the set of all derivatives of order k
1352  newDimensions[2 + fieldRank] = getDkCardinality(operatorType,spaceDim);
1353  break;
1354 
1355  default:
1356  TEUCHOS_TEST_FOR_EXCEPTION(
1357  !(Intrepid::isValidOperator(operatorType) ), std::invalid_argument,
1358  ">>> ERROR (FieldContainer): Invalid operator type");
1359  }
1360 
1361  // Resize FieldContainer using the newDimensions in the array
1362  this->resize (newDimensions);
1363 }
1364 
1365 //--------------------------------------------------------------------------------------------//
1366 // //
1367 // Methods to read and write values to FieldContainer //
1368 // //
1369 //--------------------------------------------------------------------------------------------//
1370 
1371 
1372 
1373 template<class Scalar, int ArrayTypeId>
1374 inline void FieldContainer<Scalar, ArrayTypeId>::initialize(const Scalar value) {
1375  for (int i=0; i < this->size(); i++) {
1376  data_[i] = value;
1377  }
1378 }
1379 
1380 
1381 
1382 template<class Scalar, int ArrayTypeId>
1383 inline Scalar FieldContainer<Scalar, ArrayTypeId>::getValue(const Teuchos::Array<int>& multiIndex) const {
1384  return data_[this -> getEnumeration(multiIndex)];
1385 }
1386 
1387 
1388 
1389 template<class Scalar, int ArrayTypeId>
1390 inline void FieldContainer<Scalar, ArrayTypeId>::setValue(const Scalar dataValue,
1391  const Teuchos::Array<int>& multiIndex) {
1392  data_[this -> getEnumeration(multiIndex)] = dataValue;
1393 }
1394 
1395 
1396 
1397 template<class Scalar, int ArrayTypeId>
1398 inline void FieldContainer<Scalar, ArrayTypeId>::setValue(const Scalar dataValue,
1399  const int order) {
1400  data_[order] = dataValue;
1401 }
1402 
1403 
1404 
1405 template<class Scalar, int ArrayTypeId>
1406 void FieldContainer<Scalar, ArrayTypeId>::setValues(const Teuchos::ArrayView<Scalar>& dataArray) {
1407 #ifdef HAVE_INTREPID_DEBUG
1408  TEUCHOS_TEST_FOR_EXCEPTION( (dataArray.size() != (data_.size()) ),
1409  std::invalid_argument,
1410  ">>> ERROR (FieldContainer): Size of argument does not match the size of container.");
1411 #endif
1412  data_.assign(dataArray.begin(),dataArray.end());
1413  data_ptr_ = data_.begin();
1414 }
1415 
1416 
1417 
1418 template<class Scalar, int ArrayTypeId>
1420  const int numData)
1421 {
1422 #ifdef HAVE_INTREPID_DEBUG
1423  TEUCHOS_TEST_FOR_EXCEPTION( (numData != this -> size() ), std::invalid_argument,
1424  ">>> ERROR (FieldContainer): Number of data does not match the size of container.");
1425 
1426 #endif
1427  data_.assign(dataPtr, dataPtr + numData);
1428  data_ptr_ = data_.begin();
1429 }
1430 
1431 
1432 
1433 template<class Scalar, int ArrayTypeId>
1434 inline const Scalar& FieldContainer<Scalar, ArrayTypeId>::operator () (const int i0) const
1435 {
1436 #ifdef HAVE_INTREPID_DEBUG
1437  TEUCHOS_TEST_FOR_EXCEPTION( ( this -> rank() != 1), std::invalid_argument,
1438  ">>> ERROR (FieldContainer): Number of indices does not match rank of the container.");
1439  TEUCHOS_TEST_FOR_EXCEPTION( ( (i0 < 0) || (i0 >= dim0_) ), std::invalid_argument,
1440  ">>> ERROR (FieldContainer): index is out of range.");
1441 #endif
1442  return data_ptr_[i0];
1443 }
1444 
1445 
1446 template<class Scalar, int ArrayTypeId>
1448 {
1449 #ifdef HAVE_INTREPID_DEBUG
1450  TEUCHOS_TEST_FOR_EXCEPTION( ( this -> rank() != 1), std::invalid_argument,
1451  ">>> ERROR (FieldContainer): Number of indices does not match rank of the container.");
1452  TEUCHOS_TEST_FOR_EXCEPTION( ( (i0 < 0) || (i0 >= dim0_) ), std::invalid_argument,
1453  ">>> ERROR (FieldContainer): index is out of range.");
1454 #endif
1455  return data_ptr_[i0];
1456 }
1457 
1458 
1459 
1460 template<class Scalar, int ArrayTypeId>
1461 inline const Scalar& FieldContainer<Scalar, ArrayTypeId>::operator () (const int i0,
1462  const int i1) const
1463 {
1464 #ifdef HAVE_INTREPID_DEBUG
1465  TEUCHOS_TEST_FOR_EXCEPTION( ( this -> rank() != 2), std::invalid_argument,
1466  ">>> ERROR (FieldContainer): Number of indices does not match rank of the container.");
1467  TEUCHOS_TEST_FOR_EXCEPTION( ( (i0 < 0) || (i0 >= dim0_) ), std::invalid_argument,
1468  ">>> ERROR (FieldContainer): 1st index is out of range.");
1469  TEUCHOS_TEST_FOR_EXCEPTION( ( (i1 < 0) || (i1 >= dim1_) ), std::invalid_argument,
1470  ">>> ERROR (FieldContainer): 2nd index is out of range.");
1471 #endif
1472  return data_ptr_[i0*dim1_ + i1];
1473 }
1474 
1475 
1476 template<class Scalar, int ArrayTypeId>
1478  const int i1)
1479 {
1480 #ifdef HAVE_INTREPID_DEBUG
1481  TEUCHOS_TEST_FOR_EXCEPTION( ( this -> rank() != 2), std::invalid_argument,
1482  ">>> ERROR (FieldContainer): Number of indices does not match rank of the container.");
1483  TEUCHOS_TEST_FOR_EXCEPTION( ( (i0 < 0) || (i0 >= dim0_) ), std::invalid_argument,
1484  ">>> ERROR (FieldContainer): 1st index is out of range.");
1485  TEUCHOS_TEST_FOR_EXCEPTION( ( (i1 < 0) || (i1 >= dim1_) ), std::invalid_argument,
1486  ">>> ERROR (FieldContainer): 2nd index is out of range.");
1487 #endif
1488  return data_ptr_[i0*dim1_ + i1];
1489 }
1490 
1491 
1492 
1493 template<class Scalar, int ArrayTypeId>
1494 inline const Scalar& FieldContainer<Scalar, ArrayTypeId>::operator () (const int i0,
1495  const int i1,
1496  const int i2) const
1497 {
1498 #ifdef HAVE_INTREPID_DEBUG
1499  TEUCHOS_TEST_FOR_EXCEPTION( ( this -> rank() != 3), std::invalid_argument,
1500  ">>> ERROR (FieldContainer): Number of indices does not match rank of the container.");
1501  TEUCHOS_TEST_FOR_EXCEPTION( ( (i0 < 0) || (i0 >= dim0_) ), std::invalid_argument,
1502  ">>> ERROR (FieldContainer): 1st index is out of range.");
1503  TEUCHOS_TEST_FOR_EXCEPTION( ( (i1 < 0) || (i1 >= dim1_) ), std::invalid_argument,
1504  ">>> ERROR (FieldContainer): 2nd index is out of range.");
1505  TEUCHOS_TEST_FOR_EXCEPTION( ( (i2 < 0) || (i2 >= dim2_) ), std::invalid_argument,
1506  ">>> ERROR (FieldContainer): 3rd index is out of range.");
1507 #endif
1508  return data_ptr_[(i0*dim1_ + i1)*dim2_ + i2];
1509 }
1510 
1511 template<class Scalar, int ArrayTypeId>
1513  const int i1,
1514  const int i2)
1515 {
1516 #ifdef HAVE_INTREPID_DEBUG
1517  TEUCHOS_TEST_FOR_EXCEPTION( ( this -> rank() != 3), std::invalid_argument,
1518  ">>> ERROR (FieldContainer): Number of indices does not match rank of the container.");
1519  TEUCHOS_TEST_FOR_EXCEPTION( ( (i0 < 0) || (i0 >= dim0_) ), std::invalid_argument,
1520  ">>> ERROR (FieldContainer): 1st index is out of range.");
1521  TEUCHOS_TEST_FOR_EXCEPTION( ( (i1 < 0) || (i1 >= dim1_) ), std::invalid_argument,
1522  ">>> ERROR (FieldContainer): 2nd index is out of range.");
1523  TEUCHOS_TEST_FOR_EXCEPTION( ( (i2 < 0) || (i2 >= dim2_) ), std::invalid_argument,
1524  ">>> ERROR (FieldContainer): 3rd index is out of range.");
1525 #endif
1526  return data_ptr_[(i0*dim1_ + i1)*dim2_ + i2];
1527 }
1528 
1529 
1530 
1531 template<class Scalar, int ArrayTypeId>
1532 inline const Scalar& FieldContainer<Scalar, ArrayTypeId>::operator () (const int i0,
1533  const int i1,
1534  const int i2,
1535  const int i3) const {
1536 #ifdef HAVE_INTREPID_DEBUG
1537  TEUCHOS_TEST_FOR_EXCEPTION( ( this -> rank() != 4), std::invalid_argument,
1538  ">>> ERROR (FieldContainer): Number of indices does not match rank of the container.");
1539  TEUCHOS_TEST_FOR_EXCEPTION( ( (i0 < 0) || (i0 >= dim0_) ), std::invalid_argument,
1540  ">>> ERROR (FieldContainer): 1st index is out of range.");
1541  TEUCHOS_TEST_FOR_EXCEPTION( ( (i1 < 0) || (i1 >= dim1_) ), std::invalid_argument,
1542  ">>> ERROR (FieldContainer): 2nd index is out of range.");
1543  TEUCHOS_TEST_FOR_EXCEPTION( ( (i2 < 0) || (i2 >= dim2_) ), std::invalid_argument,
1544  ">>> ERROR (FieldContainer): 3rd index is out of range.");
1545  TEUCHOS_TEST_FOR_EXCEPTION( ( (i3 < 0) || (i3 >= dim3_) ), std::invalid_argument,
1546  ">>> ERROR (FieldContainer): 4th index is out of range.");
1547 #endif
1548  return data_ptr_[( (i0*dim1_ + i1 )*dim2_ + i2 )*dim3_ + i3];
1549 }
1550 
1551 
1552 template<class Scalar, int ArrayTypeId>
1554  const int i1,
1555  const int i2,
1556  const int i3) {
1557 #ifdef HAVE_INTREPID_DEBUG
1558  TEUCHOS_TEST_FOR_EXCEPTION( ( this -> rank() != 4), std::invalid_argument,
1559  ">>> ERROR (FieldContainer): Number of indices does not match rank of the container.");
1560  TEUCHOS_TEST_FOR_EXCEPTION( ( (i0 < 0) || (i0 >= dim0_) ), std::invalid_argument,
1561  ">>> ERROR (FieldContainer): 1st index is out of range.");
1562  TEUCHOS_TEST_FOR_EXCEPTION( ( (i1 < 0) || (i1 >= dim1_) ), std::invalid_argument,
1563  ">>> ERROR (FieldContainer): 2nd index is out of range.");
1564  TEUCHOS_TEST_FOR_EXCEPTION( ( (i2 < 0) || (i2 >= dim2_) ), std::invalid_argument,
1565  ">>> ERROR (FieldContainer): 3rd index is out of range.");
1566  TEUCHOS_TEST_FOR_EXCEPTION( ( (i3 < 0) || (i3 >= dim3_) ), std::invalid_argument,
1567  ">>> ERROR (FieldContainer): 4th index is out of range.");
1568 #endif
1569  return data_ptr_[( (i0*dim1_ + i1 )*dim2_ + i2 )*dim3_ + i3];
1570 }
1571 
1572 
1573 
1574 template<class Scalar, int ArrayTypeId>
1575 inline const Scalar& FieldContainer<Scalar, ArrayTypeId>::operator () (const int i0,
1576  const int i1,
1577  const int i2,
1578  const int i3,
1579  const int i4) const {
1580 #ifdef HAVE_INTREPID_DEBUG
1581  TEUCHOS_TEST_FOR_EXCEPTION( ( this -> rank() != 5), std::invalid_argument,
1582  ">>> ERROR (FieldContainer): Number of indices does not match rank of the container.");
1583  TEUCHOS_TEST_FOR_EXCEPTION( ( (i0 < 0) || (i0 >= dim0_) ), std::invalid_argument,
1584  ">>> ERROR (FieldContainer): 1st index is out of range.");
1585  TEUCHOS_TEST_FOR_EXCEPTION( ( (i1 < 0) || (i1 >= dim1_) ), std::invalid_argument,
1586  ">>> ERROR (FieldContainer): 2nd index is out of range.");
1587  TEUCHOS_TEST_FOR_EXCEPTION( ( (i2 < 0) || (i2 >= dim2_) ), std::invalid_argument,
1588  ">>> ERROR (FieldContainer): 3rd index is out of range.");
1589  TEUCHOS_TEST_FOR_EXCEPTION( ( (i3 < 0) || (i3 >= dim3_) ), std::invalid_argument,
1590  ">>> ERROR (FieldContainer): 4th index is out of range.");
1591  TEUCHOS_TEST_FOR_EXCEPTION( ( (i4 < 0) || (i4 >= dim4_) ), std::invalid_argument,
1592  ">>> ERROR (FieldContainer): 5th index is out of range.");
1593 #endif
1594  return data_ptr_[( ( (i0*dim1_ + i1 )*dim2_ + i2 )*dim3_ + i3 )*dim4_ + i4];
1595 }
1596 
1597 template<class Scalar, int ArrayTypeId>
1599  const int i1,
1600  const int i2,
1601  const int i3,
1602  const int i4) {
1603 #ifdef HAVE_INTREPID_DEBUG
1604  TEUCHOS_TEST_FOR_EXCEPTION( ( this -> rank() != 5), std::invalid_argument,
1605  ">>> ERROR (FieldContainer): Number of indices does not match rank of the container.");
1606  TEUCHOS_TEST_FOR_EXCEPTION( ( (i0 < 0) || (i0 >= dim0_) ), std::invalid_argument,
1607  ">>> ERROR (FieldContainer): 1st index is out of range.");
1608  TEUCHOS_TEST_FOR_EXCEPTION( ( (i1 < 0) || (i1 >= dim1_) ), std::invalid_argument,
1609  ">>> ERROR (FieldContainer): 2nd index is out of range.");
1610  TEUCHOS_TEST_FOR_EXCEPTION( ( (i2 < 0) || (i2 >= dim2_) ), std::invalid_argument,
1611  ">>> ERROR (FieldContainer): 3rd index is out of range.");
1612  TEUCHOS_TEST_FOR_EXCEPTION( ( (i3 < 0) || (i3 >= dim3_) ), std::invalid_argument,
1613  ">>> ERROR (FieldContainer): 4th index is out of range.");
1614  TEUCHOS_TEST_FOR_EXCEPTION( ( (i4 < 0) || (i4 >= dim4_) ), std::invalid_argument,
1615  ">>> ERROR (FieldContainer): 5th index is out of range.");
1616 #endif
1617  return data_ptr_[( ( (i0*dim1_ + i1 )*dim2_ + i2 )*dim3_ + i3 )*dim4_ + i4];
1618 }
1619 
1620 
1621 
1622 template<class Scalar, int ArrayTypeId>
1623 const Scalar& FieldContainer<Scalar, ArrayTypeId>::operator [] (const int address) const {
1624 #ifdef HAVE_INTREPID_DEBUG
1625  TEUCHOS_TEST_FOR_EXCEPTION( ( (address < 0) || (address >= (int)data_.size() ) ),
1626  std::invalid_argument,
1627  ">>> ERROR (FieldContainer): Specified address is out of range.");
1628 #endif
1629  return data_ptr_[address];
1630 }
1631 
1632 
1633 
1634 template<class Scalar, int ArrayTypeId>
1636 #ifdef HAVE_INTREPID_DEBUG
1637  TEUCHOS_TEST_FOR_EXCEPTION( ( (address < 0) || (address >= (int)data_.size() ) ),
1638  std::invalid_argument,
1639  ">>> ERROR (FieldContainer): Specified address is out of range.");
1640 #endif
1641  return data_ptr_[address];
1642 }
1643 
1644 
1645 
1646 template<class Scalar, int ArrayTypeId>
1648 {
1649 #ifdef HAVE_INTREPID_DEBUG
1650  TEUCHOS_TEST_FOR_EXCEPTION( ( this == &right ),
1651  std::invalid_argument,
1652  ">>> ERROR (FieldContainer): Invalid right-hand side to '='. Self-assignment prohibited.");
1653 #endif
1654  dim0_ = right.dim0_;
1655  dim1_ = right.dim1_;
1656  dim2_ = right.dim2_;
1657  dim3_ = right.dim3_;
1658  dim4_ = right.dim4_;
1659  data_.deepCopy(right.data_());
1660  data_ptr_ = data_.begin();
1661  dimensions_ = right.dimensions_;
1662  return *this;
1663 }
1664 
1665 
1666 //===========================================================================//
1667 // //
1668 // END of member definitions; START friends and related //
1669 // //
1670 //===========================================================================//
1671 
1672 
1673 template<class Scalar, int ArrayTypeId>
1674 std::ostream& operator << (std::ostream& os, const FieldContainer<Scalar, ArrayTypeId>& container) {
1675 
1676  // Save the format state of the original ostream os.
1677  Teuchos::oblackholestream oldFormatState;
1678  oldFormatState.copyfmt(os);
1679 
1680  os.setf(std::ios_base::scientific, std::ios_base::floatfield);
1681  os.setf(std::ios_base::right);
1682  int myprec = os.precision();
1683 
1684  int size = container.size();
1685  int rank = container.rank();
1686  Teuchos::Array<int> multiIndex(rank);
1687  Teuchos::Array<int> dimensions;
1688  container.dimensions(dimensions);
1689 
1690  os<< "===============================================================================\n"\
1691  << "\t Container size = " << size << "\n"
1692  << "\t Container rank = " << rank << "\n" ;
1693 
1694  if( (rank == 0 ) && (size == 0) ) {
1695  os<< "====================================================================================\n"\
1696  << "| *** This is an empty container **** |\n";
1697  }
1698  else {
1699  os<< "\t Dimensions = ";
1700 
1701  for(int r = 0; r < rank; r++){
1702  os << " (" << dimensions[r] <<") ";
1703  }
1704  os << "\n";
1705 
1706  os<< "====================================================================================\n"\
1707  << "| Multi-index Enumeration Value |\n"\
1708  << "====================================================================================\n";
1709  }
1710 
1711  for(int address = 0; address < size; address++){
1712  container.getMultiIndex(multiIndex,address);
1713  std::ostringstream mistring;
1714  for(int r = 0; r < rank; r++){
1715  mistring << multiIndex[r] << std::dec << " ";
1716  }
1717  os.setf(std::ios::right, std::ios::adjustfield);
1718  os << std::setw(27) << mistring.str();
1719  os << std::setw(20) << address;
1720  os << " ";
1721  os.setf(std::ios::left, std::ios::adjustfield);
1722  os << std::setw(myprec+8) << container[address] << "\n";
1723  }
1724 
1725  os<< "====================================================================================\n\n";
1726 
1727  // reset format state of os
1728  os.copyfmt(oldFormatState);
1729 
1730  return os;
1731 }
1732 // End member, friend, and related function definitions of class FieldContainer.
1733 
1734 } // end namespace Intrepid
int size() const
Returns size of the FieldContainer defined as the product of its dimensions.
const Scalar & operator[](const int address) const
Overloaded [] operator. Returns value based on its enumeration. Data cannot be modified.
FieldContainer & operator=(const FieldContainer &right)
Assignment operator *this = right.
int dim3_
4th dimension of the array
Teuchos::Array< int > dimensions_
Array to store dimensions (dimensions) for the multi-indices. Admissible range (dimension) for the k-...
void clear()
Clears FieldContainer to trivial container (one with rank = 0 and size = 0)
int dim1_
2nd dimension of the array
int dimension(const int whichDim) const
Returns the specified dimension.
int dim2_
3rd dimension of the array
Teuchos::ArrayRCP< Scalar > data_
Array to store the multi-indexed quantity.
void getMultiIndex(int &i0, const int valueEnum) const
Returns the multi-index of a value, based on its enumeration, as a list, for rank-1 containers...
int rank() const
Return rank of the FieldContainer = number of indices used to tag the multi-indexed value...
int getEnumeration(const int i0) const
Returns enumeration of a value (its order relative to the container), based on its multi-index...
const Scalar & operator()(const int i0) const
Overloaded () operators for rank-1 containers. Data cannot be modified.
void setValue(const Scalar dataValue, const Teuchos::Array< int > &multiIndex)
Assign value by its multi-index.
void setValues(const Teuchos::ArrayView< Scalar > &dataArray)
Fills an existing FieldContainer with Scalars stored in a Teuchos::Array without changing rank and di...
Implementation of a templated lexicographical container for a multi-indexed scalar quantity...
void resize(const int dim0)
Resizes FieldContainer to a rank-1 container with the specified dimension, initialized by 0...
FieldContainer()
Default constructor.
void initialize(const Scalar value=0)
Initializes a field container by assigning value to all its elements.
int dim0_
1st dimension of the array
int dim4_
5th dimension of the array
Scalar getValue(const Teuchos::Array< int > &multiIndex) const
Retrieve value by its multi-index. To retrieve it by enumeration use the overloaded []...
void dimensions(Vector &dimensions) const
Returns array with the dimensions of the container.