Intrepid
example_01.cpp
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 
50 #include "Teuchos_Time.hpp"
51 #include "Teuchos_GlobalMPISession.hpp"
52 
53 using namespace std;
54 using namespace Intrepid;
55 
56 int main(int argc, char *argv[]) {
57 
58  Teuchos::GlobalMPISession mpiSession(&argc, &argv);
59 
60  std::cout \
61  << "===============================================================================\n" \
62  << "| |\n" \
63  << "| Example use of the FieldContainer class |\n" \
64  << "| |\n" \
65  << "| 1) Creating and filling FieldContainer objects |\n" \
66  << "| 2) Accessing elements in FieldContainer objects |\n" \
67  << "| |\n" \
68  << "| Questions? Contact Pavel Bochev (pbboche@sandia.gov) or |\n" \
69  << "| Denis Ridzal (dridzal@sandia.gov). |\n" \
70  << "| |\n" \
71  << "| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
72  << "| Trilinos website: http://trilinos.sandia.gov |\n" \
73  << "| |\n" \
74  << "===============================================================================\n\n";
75 
76  // Define variables to create and use FieldContainers
77  Teuchos::Array<int> dimension;
78  Teuchos::Array<int> multiIndex;
79 
80  std::cout \
81  << "===============================================================================\n"\
82  << "| EXAMPLE 1: rank 2 multi-index: {u(p,i) | 0 <= p < 5; 0 <= i < 3 } |\n"\
83  << "===============================================================================\n\n";
84  // This rank 2 multi-indexed value can be used to store values of 3D vector field
85  // evaluated at 5 points, or values of gradient of scalar field evaluated at the points
86 
87  // Resize dimension and multiIndex for rank-2 multi-indexed value
88  dimension.resize(2);
89  multiIndex.resize(2);
90 
91  // Load upper ranges for the two indices in the multi-indexed value
92  dimension[0] = 5;
93  dimension[1] = 3;
94 
95  // Create FieldContainer that can hold the rank-2 multi-indexed value
96  FieldContainer<double> myContainer(dimension);
97 
98  // Fill with some data: leftmost index changes last, rightmost index changes first!
99  for(int p = 0; p < dimension[0]; p++){
100  multiIndex[0] = p;
101 
102  for(int i = 0; i < dimension[1]; i++){
103  multiIndex[1] = i;
104 
105  // Load value with multi-index {p,i} to container
106  myContainer.setValue((double)(i+p), multiIndex);
107  }
108  }
109 
110  // Show container contents
111  std::cout << myContainer;
112 
113  // Access by overloaded (), multiindex and []:
114  multiIndex[0] = 3;
115  multiIndex[1] = 1;
116  int enumeration = myContainer.getEnumeration(multiIndex);
117 
118  std::cout << "Access by (): myContainer(" << 3 <<"," << 1 << ") = " << myContainer(3,1) << "\n";
119  std::cout << "Access by multi-index: myContainer{" << multiIndex[0] << multiIndex[1] << "} = "<< myContainer.getValue(multiIndex) <<"\n";
120  std::cout << "Access by enumeration: myContainer[" << enumeration << "] = " << myContainer[enumeration] <<"\n";
121 
122  std::cout << "Assigning value by (): \n old value at (3,1) = " << myContainer(3,1) <<"\n";
123  myContainer(3,1) = 999.99;
124  std::cout << " new value at (3,1) = " << myContainer(3,1) <<"\n";
125 
126 
127  std::cout << "\n" \
128  << "===============================================================================\n"\
129  << "| EXAMPLE 2: rank 3 multi-index: {u(p,i,j) | 0 <=p< 5; 0 <= i<2, 0<=j<3 |\n"\
130  << "===============================================================================\n\n";
131  // This rank-3 value can be used to store subset of second partial derivatives values
132  // of a scalar function at p points.
133 
134  // Resize dimension and multiIndex for rank-3 multi-indexed value
135  dimension.resize(3);
136  multiIndex.resize(3);
137 
138  // Define upper ranges for the three indices in the multi-indexed value
139  dimension[0] = 5;
140  dimension[1] = 2;
141  dimension[2] = 3;
142 
143  // Reset the existing container to accept rank-3 value with the specified index ranges
144  myContainer.resize(dimension);
145 
146  // Fill with some data
147  for(int p = 0; p < dimension[0]; p++){
148  multiIndex[0] = p;
149  for(int i = 0; i < dimension[1]; i++){
150  multiIndex[1] = i;
151  for(int j = 0; j < dimension[2]; j++){
152  multiIndex[2] = j;
153 
154  // Load value with multi-index {p,i} to container
155  myContainer.setValue((double)(p+i+j), multiIndex);
156  }
157  }
158  }
159 
160  // Display contents
161  std::cout << myContainer;
162 
163  // Access by overloaded (), multiindex and []:
164  multiIndex[0] = 3;
165  multiIndex[1] = 1;
166  multiIndex[2] = 2;
167  enumeration = myContainer.getEnumeration(multiIndex);
168 
169  std::cout << "Access by (): myContainer(" << 3 <<"," << 1 << "," << 2 << ") = " << myContainer(3,1,2) << "\n";
170  std::cout << "Access by multi-index: myContainer{" << multiIndex[0] << multiIndex[1] << multiIndex[2] << "} = "<< myContainer.getValue(multiIndex) <<"\n";
171  std::cout << "Access by enumeration: myContainer[" << enumeration << "] = " << myContainer[enumeration] <<"\n";
172 
173  std::cout << "Assigning value by (): \n old value at (3,1,2) = " << myContainer(3,1,2) <<"\n";
174  myContainer(3,1,2) = -999.999;
175  std::cout << " new value at (3,1,2) = " << myContainer(3,1,2) <<"\n";
176 
177 
178  std::cout << "\n" \
179  << "===============================================================================\n"\
180  << "| EXAMPLE 4: making rank-5 FieldContainer from data array and index range array |\n"\
181  << "===============================================================================\n\n";
182  // Initialize dimension for rank-5 multi-index value
183  dimension.resize(5);
184  dimension[0] = 5;
185  dimension[1] = 2;
186  dimension[2] = 3;
187  dimension[3] = 4;
188  dimension[4] = 6;
189 
190  // Define Teuchos::Array to store values with dimension equal to the number of multi-indexed values
191  Teuchos::Array<double> dataTeuchosArray(5*2*3*4*6);
192 
193  // Fill with data
194  int counter = 0;
195  for(int i=0; i < dimension[0]; i++){
196  for(int j=0; j < dimension[1]; j++){
197  for(int k=0; k < dimension[2]; k++){
198  for(int l = 0; l < dimension[3]; l++){
199  for(int m = 0; m < dimension[4]; m++){
200  dataTeuchosArray[counter] = (double)counter;
201  counter++;
202  }
203  }
204  }
205  }
206  }
207 
208  // Create FieldContainer from data array and index array and show it
209  FieldContainer<double> myNewContainer(dimension, dataTeuchosArray);
210  std::cout << myNewContainer;
211 
212  // Access by overloaded (), multiindex and []:
213  multiIndex.resize(myNewContainer.rank());
214  multiIndex[0] = 3;
215  multiIndex[1] = 1;
216  multiIndex[2] = 2;
217  multiIndex[3] = 2;
218  multiIndex[4] = 5;
219  enumeration = myNewContainer.getEnumeration(multiIndex);
220 
221  std::cout << "Access by (): myNewContainer(" << 3 <<"," << 1 << "," << 2 << "," << 2 << "," << 5 << ") = " << myNewContainer(3,1,2,2,5) << "\n";
222  std::cout << "Access by multi-index: myNewContainer{" << multiIndex[0] << multiIndex[1] << multiIndex[2] << multiIndex[3] << multiIndex[4] << "} = "<< myNewContainer.getValue(multiIndex) <<"\n";
223  std::cout << "Access by enumeration: myNewContainer[" << enumeration << "] = " << myNewContainer[enumeration] <<"\n";
224 
225  std::cout << "Assigning value by (): \n old value at (3,1,2,2,5) = " << myNewContainer(3,1,2,2,5) <<"\n";
226  myNewContainer(3,1,2,2,5) = -888.888;
227  std::cout << " new value at (3,1,2,2,5) = " << myNewContainer(3,1,2,2,5) <<"\n";
228 
229 
230  std::cout << "\n" \
231  << "===============================================================================\n"\
232  << "| EXAMPLE 5: making trivial FieldContainers and storing a single zero |\n"\
233  << "===============================================================================\n\n";
234 
235  // Make trivial container by resetting the index range to zero rank (no indices) and then
236  // using resize method
237  dimension.resize(0);
238  myContainer.resize(dimension);
239  std::cout << myContainer;
240 
241  // Make trivial container by using clear method:
242  myNewContainer.clear();
243  std::cout << myNewContainer;
244 
245  // Now use initialize() to reset the container to hold a single zero
246  myNewContainer.initialize();
247  std::cout << myNewContainer;
248 
249 
250  std::cout << "\n" \
251  << "===============================================================================\n"\
252  << "| EXAMPLE 6: Timing read and write operations using () and getValue |\n"\
253  << "===============================================================================\n\n";
254 
255  // Initialize dimensions for rank-5 multi-index value
256  int dim0 = 10; // number of cells
257  int dim1 = 50; // number of points
258  int dim2 = 27; // number of functions
259  int dim3 = 3; // 1st space dim
260  int dim4 = 3; // 2nd space dim
261 
262  FieldContainer<double> myTensorContainer(dim0, dim1, dim2, dim3, dim4);
263  multiIndex.resize(myTensorContainer.rank());
264 // double aValue;
265 
266  Teuchos::Time timerGetValue("Reading and writing from rank-5 container using getValue");
267  timerGetValue.start();
268  for(int i0 = 0; i0 < dim0; i0++){
269  multiIndex[0] = i0;
270  for(int i1 = 0; i1 < dim1; i1++){
271  multiIndex[1] = i1;
272  for(int i2 = 0; i2 < dim2; i2++) {
273  multiIndex[2] = i2;
274  for(int i3 = 0; i3 < dim3; i3++) {
275  multiIndex[3] = i3;
276  for(int i4 =0; i4 < dim4; i4++) {
277  multiIndex[4] = i4;
278 
279  //aValue = myTensorContainer.getValue(multiIndex);
280  myTensorContainer.getValue(multiIndex);
281  myTensorContainer.setValue(999.999,multiIndex);
282 
283  }
284  }
285  }
286  }
287  }
288  timerGetValue.stop();
289  std::cout << " Time to read and write from container using getValue: " << timerGetValue.totalElapsedTime() <<"\n";
290 
291  Teuchos::Time timerRound("Reading and writing from rank-5 container using ()");
292  timerRound.start();
293  for(int i0 = 0; i0 < dim0; i0++){
294  for(int i1 = 0; i1 < dim1; i1++) {
295  for(int i2 = 0; i2 < dim2; i2++) {
296  for(int i3 = 0; i3 < dim3; i3++) {
297  for(int i4 =0; i4 < dim4; i4++) {
298 
299  // aValue = myTensorContainer(i0,i1,i2,i3,i4);
300  myTensorContainer(i0,i1,i2,i3,i4) = 999.999;
301 
302  }
303  }
304  }
305  }
306  }
307  timerRound.stop();
308  std::cout << " Time to read and write from container using (): " << timerRound.totalElapsedTime() <<"\n";
309 
310  std::cout << "\n" \
311  << "===============================================================================\n"\
312  << "| EXAMPLE 6: Specialized methods of FieldContainer |\n"\
313  << "===============================================================================\n\n";
314 
315  return 0;
316 }
317 
318 
319 
320 
321 
Header file for utility class to provide multidimensional containers.