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 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-...
int size() const
Returns size of the FieldContainer defined as the product of its dimensions.
int dim1_
2nd dimension of the array
int dim2_
3rd dimension of the array
Teuchos::ArrayRCP< Scalar > data_
Array to store the multi-indexed quantity.
void dimensions(Vector &dimensions) const
Returns array with the dimensions of the container.
int getOperatorRank(const EFunctionSpace spaceType, const EOperator operatorType, const int spaceDim)
Returns rank of an operator.
EOperator
Enumeration of primitive operators available in Intrepid. Primitive operators act on reconstructed fu...
int getFieldRank(const EFunctionSpace spaceType)
Returns the rank of fields in a function space of the specified type.
Implementation of a templated lexicographical container for a multi-indexed scalar quantity...
int getDkCardinality(const EOperator operatorType, const int spaceDim)
Returns cardinality of Dk, i.e., the number of all derivatives of order k.
void resize(const int dim0)
Resizes FieldContainer to a rank-1 container with the specified dimension, initialized by 0...
FieldContainer()
Default constructor.
int isValidOperator(const EOperator operatorType)
Verifies validity of an operator enum.
int dim0_
1st dimension of the array
int dim4_
5th dimension of the array