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
49namespace Intrepid {
50
51
52 //--------------------------------------------------------------------------------------------//
53 // //
54 // Member function definitions of the class FieldContainer //
55 // //
56 //--------------------------------------------------------------------------------------------//
57
58
59template<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
80template<class Scalar, int ArrayTypeId>
81FieldContainer<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
98template<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
120template<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
145template<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();
171
172
173
174template<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_;
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
207template<class Scalar, int ArrayTypeId>
208FieldContainer<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;
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
271
272
273
274template<class Scalar, int ArrayTypeId>
275FieldContainer<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}
345
346
347template<class Scalar, int ArrayTypeId>
348FieldContainer<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];
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
420template<class Scalar, int ArrayTypeId>
421FieldContainer<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
495template<class Scalar, int ArrayTypeId>
496FieldContainer<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);
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
593template<class Scalar, int ArrayTypeId>
595 return dimensions_.size();
596}
597
598
599
600template<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;
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;
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
649template<class Scalar, int ArrayTypeId>
650template<class Vector>
651inline void FieldContainer<Scalar, ArrayTypeId>::dimensions(Vector& dimVec) const {
652 dimVec.assign(dimensions_.begin(),dimensions_.end());
653}
654
655
656
657template<class Scalar, int ArrayTypeId>
658inline 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
670template<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
683template<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
698
699template<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
718template<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}
738
739
740template<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
766template<class Scalar, int ArrayTypeId>
767int 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
862template<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
878template<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
897template<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
925template<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
959template<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
997template<class Scalar, int ArrayTypeId>
998template<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
1047template<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
1065template<class Scalar, int ArrayTypeId>
1066void 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
1136template<class Scalar, int ArrayTypeId>
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
1151template<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
1168template<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
1187template<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
1208template<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
1231template<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
1287template<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
1373template<class Scalar, int ArrayTypeId>
1375 for (int i=0; i < this->size(); i++) {
1376 data_[i] = value;
1377 }
1378}
1379
1380
1381
1382template<class Scalar, int ArrayTypeId>
1383inline Scalar FieldContainer<Scalar, ArrayTypeId>::getValue(const Teuchos::Array<int>& multiIndex) const {
1384 return data_[this -> getEnumeration(multiIndex)];
1385}
1386
1387
1388
1389template<class Scalar, int ArrayTypeId>
1390inline void FieldContainer<Scalar, ArrayTypeId>::setValue(const Scalar dataValue,
1391 const Teuchos::Array<int>& multiIndex) {
1392 data_[this -> getEnumeration(multiIndex)] = dataValue;
1393}
1394
1395
1396
1397template<class Scalar, int ArrayTypeId>
1398inline void FieldContainer<Scalar, ArrayTypeId>::setValue(const Scalar dataValue,
1399 const int order) {
1400 data_[order] = dataValue;
1401}
1402
1403
1404
1405template<class Scalar, int ArrayTypeId>
1406void 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
1418template<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
1433template<class Scalar, int ArrayTypeId>
1434inline 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
1446template<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
1460template<class Scalar, int ArrayTypeId>
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
1476template<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
1493template<class Scalar, int ArrayTypeId>
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
1511template<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
1531template<class Scalar, int ArrayTypeId>
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
1552template<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
1574template<class Scalar, int ArrayTypeId>
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
1597template<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
1622template<class Scalar, int ArrayTypeId>
1623const 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
1634template<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
1646template<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
1673template<class Scalar, int ArrayTypeId>
1674std::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 isValidOperator(const EOperator operatorType)
Verifies validity of an operator enum.
int getFieldRank(const EFunctionSpace spaceType)
Returns the rank of fields in a function space of the specified type.
int getOperatorRank(const EFunctionSpace spaceType, const EOperator operatorType, const int spaceDim)
Returns rank of an operator.
int getDkCardinality(const EOperator operatorType, const int spaceDim)
Returns cardinality of Dk, i.e., the number of all derivatives of order k.
Implementation of a templated lexicographical container for a multi-indexed scalar quantity....
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 dim3_
4th dimension of the array
const Scalar & operator()(const int i0) const
Overloaded () operators for rank-1 containers. Data cannot be modified.
void dimensions(Vector &dimensions) const
Returns array with the dimensions of the container.
void resize(const int dim0)
Resizes FieldContainer to a rank-1 container with the specified dimension, initialized by 0.
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.
Teuchos::ArrayRCP< Scalar > data_
Array to store the multi-indexed quantity.
void clear()
Clears FieldContainer to trivial container (one with rank = 0 and size = 0)
void setValue(const Scalar dataValue, const Teuchos::Array< int > &multiIndex)
Assign value by its multi-index.
int dim0_
1st dimension of the array
int getEnumeration(const int i0) const
Returns enumeration of a value (its order relative to the container), based on its multi-index,...
int rank() const
Return rank of the FieldContainer = number of indices used to tag the multi-indexed value.
int dim4_
5th dimension of the array
FieldContainer & operator=(const FieldContainer &right)
Assignment operator *this = right.
const Scalar & operator[](const int address) const
Overloaded [] operator. Returns value based on its enumeration. Data cannot be modified.
void setValues(const Teuchos::ArrayView< Scalar > &dataArray)
Fills an existing FieldContainer with Scalars stored in a Teuchos::Array without changing rank and di...
void initialize(const Scalar value=0)
Initializes a field container by assigning value to all its elements.
int dimension(const int whichDim) const
Returns the specified dimension.
FieldContainer()
Default constructor.
int dim1_
2nd dimension of the array
int dim2_
3rd 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 [].