Intrepid
test_03.cpp
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 
44 
52 #include "Shards_Array.hpp"
53 #include "Teuchos_oblackholestream.hpp"
54 #include "Teuchos_RCP.hpp"
55 #include "Teuchos_GlobalMPISession.hpp"
56 
57 
58 using namespace Intrepid;
59 
60 SHARDS_ARRAY_DIM_TAG_SIMPLE_DECLARATION( Cell )
61 SHARDS_ARRAY_DIM_TAG_SIMPLE_IMPLEMENTATION( Cell )
62 
63 SHARDS_ARRAY_DIM_TAG_SIMPLE_DECLARATION( Field )
64 SHARDS_ARRAY_DIM_TAG_SIMPLE_IMPLEMENTATION( Field )
65 
66 SHARDS_ARRAY_DIM_TAG_SIMPLE_DECLARATION( Point )
67 SHARDS_ARRAY_DIM_TAG_SIMPLE_IMPLEMENTATION( Point )
68 
69 SHARDS_ARRAY_DIM_TAG_SIMPLE_DECLARATION( Dim )
70 SHARDS_ARRAY_DIM_TAG_SIMPLE_IMPLEMENTATION( Dim )
71 
72 #define INTREPID_TEST_COMMAND( S ) \
73 { \
74  try { \
75  S ; \
76  } \
77  catch (std::logic_error err) { \
78  *outStream << "Expected Error ----------------------------------------------------------------\n"; \
79  *outStream << err.what() << '\n'; \
80  *outStream << "-------------------------------------------------------------------------------" << "\n\n"; \
81  }; \
82 }
83 
84 
85 int main(int argc, char *argv[]) {
86 
87  Teuchos::GlobalMPISession mpiSession(&argc, &argv);
88 
89  // This little trick lets us print to cout only if a (dummy) command-line argument is provided.
90  int iprint = argc - 1;
91 
92  Teuchos::RCP<std::ostream> outStream;
93  Teuchos::oblackholestream bhs; // outputs nothing
94 
95  if (iprint > 0)
96  outStream = Teuchos::rcp(&std::cout, false);
97  else
98  outStream = Teuchos::rcp(&bhs, false);
99 
100  // Save the format state of the original cout .
101  Teuchos::oblackholestream oldFormatState;
102  oldFormatState.copyfmt(std::cout);
103 
104  *outStream \
105  << "===============================================================================\n" \
106  << "| |\n" \
107  << "| Unit Test FieldContainer |\n" \
108  << "| |\n" \
109  << "| 1) Testing usage of various constructors / wrappers |\n" \
110  << "| 2) Testing usage of resize |\n" \
111  << "| |\n" \
112  << "| Questions? Contact Pavel Bochev (pbboche@sandia.gov) or |\n" \
113  << "| Denis Ridzal (dridzal@sandia.gov). |\n" \
114  << "| |\n" \
115  << "| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
116  << "| Trilinos website: http://trilinos.sandia.gov |\n" \
117  << "| |\n" \
118  << "===============================================================================\n";
119 
120 
121  // Set error flag.
122  int errorFlag = 0;
123 
124  double zero = INTREPID_TOL;
125 
126  try {
127 
128  // Dimensions array.
129  Teuchos::Array<int> dimensions;
130 
131  *outStream << "\n" \
132  << "===============================================================================\n"\
133  << "| TEST 1: Constructors / Wrappers for a particular rank-4 container |\n"\
134  << "===============================================================================\n\n";
135 
136  { // start variable scope
137 
138  // Initialize dimensions for rank-4 multi-index value
139  dimensions.resize(4);
140  dimensions[0] = 5;
141  dimensions[1] = 3;
142  dimensions[2] = 2;
143  dimensions[3] = 7;
144 
145  // Allocate data through a Teuchos::Array, initialized to 0.
146  Teuchos::Array<double> data(dimensions[0]*dimensions[1]*dimensions[2]*dimensions[3]);
147 
148  // Create a FieldContainer using a deep copy via Teuchos::Array.
149  FieldContainer<double> fc_array(dimensions, data);
150 
151  // modify the (1,1,1,1) entry
152  fc_array(1,1,1,1) = 1.0;
153  // verify that the data array has NOT changed
154  if (std::abs(data[dimensions[1]*dimensions[2]*dimensions[3] + dimensions[2]*dimensions[3] + dimensions[3] + 1]) > zero) {
155  *outStream << "\n\nError in constructor using Array (ArrayView) and deep copy.\n\n";
156  errorFlag = -1000;
157  }
158 
159  // test getData access function
160  if (std::abs((fc_array.getData())[dimensions[1]*dimensions[2]*dimensions[3] +
161  dimensions[2]*dimensions[3] + dimensions[3] + 1] - fc_array(1,1,1,1)) > zero) {
162  *outStream << "\n\nError in getData() member of FieldContainer.\n\n";
163  errorFlag = -1000;
164  }
165 
166  // Create a FieldContainer using a deep copy via Teuchos::ArrayRCP.
167  Teuchos::RCP< Teuchos::Array<double> > rcp_to_data = rcpFromRef(data); // first create RCP
168  Teuchos::ArrayRCP<double> arrayrcp = Teuchos::arcp(rcp_to_data); // now create ArrayRCP
169  FieldContainer<double> fc_arrayrcp_deep(dimensions, arrayrcp()); // IMPORTANT!!!: use '()' after arrayrcp
170  // for direct conversion to ArrayView
171  // modify the (1,1,1,1) entry
172  fc_arrayrcp_deep(1,1,1,1) = 1.0;
173  // verify that the data array has NOT changed
174  if (std::abs(data[dimensions[1]*dimensions[2]*dimensions[3] + dimensions[2]*dimensions[3] + dimensions[3] + 1]) > zero) {
175  *outStream << "\n\nError in constructor using ArrayRCP (ArrayView) and deep copy.\n\n";
176  errorFlag = -1000;
177  }
178 
179  // Create a FieldContainer using a shallow copy via Teuchos::ArrayRCP.
180  FieldContainer<double> fc_arrayrcp_shallow(dimensions, arrayrcp);
181  // modify the (1,1,1,1) entry
182  fc_arrayrcp_shallow(1,1,1,1) = 1.0;
183  // verify that the data array HAS changed
184  if (std::abs(data[dimensions[1]*dimensions[2]*dimensions[3] + dimensions[2]*dimensions[3] + dimensions[3] + 1] - fc_arrayrcp_shallow(1,1,1,1)) > zero) {
185  *outStream << "\n\nError in constructor using ArrayRCP and shallow copy.\n\n";
186  errorFlag = -1000;
187  }
188 
189  // Create a FieldContainer using a deep copy via Scalar*.
190  FieldContainer<double> fc_scalarptr_deep(dimensions, data.getRawPtr(), true);
191  // modify the (1,1,1,1) entry
192  fc_scalarptr_deep(1,1,1,1) = 2.0;
193  // verify that the data array has NOT changed
194  if (std::abs(data[dimensions[1]*dimensions[2]*dimensions[3] + dimensions[2]*dimensions[3] + dimensions[3] + 1] - 1.0) > zero) {
195  *outStream << "\n\nError in constructor using Scalar* and deep copy.\n\n";
196  errorFlag = -1000;
197  }
198 
199  // Create a FieldContainer using a shallow copy via Scalar*.
200  FieldContainer<double> fc_scalarptr_shallow(dimensions, data.getRawPtr());
201  // modify the (1,1,1,1) entry
202  fc_scalarptr_shallow(1,1,1,1) = 2.0;
203  // verify that the data array HAS changed
204  if (std::abs(data[dimensions[1]*dimensions[2]*dimensions[3] + dimensions[2]*dimensions[3] + dimensions[3] + 1] - fc_scalarptr_shallow(1,1,1,1)) > zero) {
205  *outStream << "\n\nError in constructor using Scalar* and shallow copy.\n\n";
206  errorFlag = -1000;
207  }
208 
209  // Create a FieldContainer using a deep copy via shards::Array.
210  shards::Array<double,shards::NaturalOrder,Cell,Field,Point,Dim> shards_array(data.getRawPtr(),dimensions[0],dimensions[1],dimensions[2],dimensions[3]);
211  FieldContainer<double> fc_shards_deep(shards_array, true);
212  // modify the (1,1,1,1) entry
213  fc_shards_deep(1,1,1,1) = 3.0;
214  // verify that the data array has NOT changed
215  if (std::abs(data[dimensions[1]*dimensions[2]*dimensions[3] + dimensions[2]*dimensions[3] + dimensions[3] + 1] - 2.0) > zero) {
216  *outStream << "\n\nError in constructor using shards::Array and deep copy.\n\n";
217  errorFlag = -1000;
218  }
219 
220  // Create a FieldContainer using a shallow copy via shards::Array.
221  FieldContainer<double> fc_shards_shallow(shards_array);
222  // modify the (1,1,1,1) entry
223  fc_shards_shallow(1,1,1,1) = 3.0;
224  // verify that the data array HAS changed
225  if (std::abs(data[dimensions[1]*dimensions[2]*dimensions[3] + dimensions[2]*dimensions[3] + dimensions[3] + 1] - fc_shards_shallow(1,1,1,1)) > zero) {
226  *outStream << "\n\nError in constructor using shards::Array and shallow copy.\n\n";
227  errorFlag = -1000;
228  }
229 
230  } // end variable scope
231 
232 
233  *outStream << "\n" \
234  << "===============================================================================\n"\
235  << "| TEST 1 cont'd: Run through constructors / wrappers of various ranks |\n"\
236  << "===============================================================================\n\n";
237 
238  for (int rank=0; rank<10; rank++) {
239  dimensions.resize(rank);
240  int total_size = 1;
241  if (rank == 0) {
242  total_size = 0;
243  }
244  for (int dim=0; dim<rank; dim++) {
245  dimensions[dim] = 2;
246  total_size *= dimensions[dim];
247  }
248 
249  // Allocate data through a Teuchos::Array, initialized to 0.
250  Teuchos::Array<double> data(total_size);
251 
252  // Create a FieldContainer using a deep copy via Teuchos::Array.
253  FieldContainer<double> fc_array(dimensions, data);
254  // Create a FieldContainer using a deep copy via Teuchos::ArrayRCP.
255  Teuchos::RCP< Teuchos::Array<double> > rcp_to_data = rcpFromRef(data); // first create RCP
256  Teuchos::ArrayRCP<double> arrayrcp = Teuchos::arcp(rcp_to_data); // now create ArrayRCP
257  FieldContainer<double> fc_arrayrcp_deep(dimensions, arrayrcp()); // IMPORTANT!!!: use '()' after arrayrcp
258  // for direct conversion to ArrayView
259  // Create a FieldContainer using a shallow copy via Teuchos::ArrayRCP.
260  FieldContainer<double> fc_arrayrcp_shallow(dimensions, arrayrcp);
261  // Create a FieldContainer using a deep copy via Scalar*.
262  FieldContainer<double> fc_scalarptr_deep(dimensions, data.getRawPtr(), true);
263  // Create a FieldContainer using a shallow copy via Scalar*.
264  FieldContainer<double> fc_scalarptr_shallow(dimensions, data.getRawPtr());
265  }
266 
267  { // start variable scope
268  // Create FieldContainers using a deep copy via shards::Array.
269  Teuchos::Array<double> data(2*2*2*2*2*2);
270  shards::Array<double,shards::NaturalOrder,Cell> shards_array_c(&data[0],2);
271  shards::Array<double,shards::NaturalOrder,Cell,Field> shards_array_cf(&data[0],2,2);
272  shards::Array<double,shards::NaturalOrder,Cell,Field,Point> shards_array_cfp(&data[0],2,2,2);
273  shards::Array<double,shards::NaturalOrder,Cell,Field,Point,Dim> shards_array_cfpd(&data[0],2,2,2,2);
274  shards::Array<double,shards::NaturalOrder,Cell,Field,Point,Dim,Dim> shards_array_cfpdd(&data[0],2,2,2,2,2);
275  shards::Array<double,shards::NaturalOrder,Cell,Field,Point,Dim,Dim,Dim> shards_array_cfpddd(&data[0],2,2,2,2,2,2);
276  FieldContainer<double> fc_shards_c_deep(shards_array_c, true);
277  FieldContainer<double> fc_shards_cf_deep(shards_array_cf, true);
278  FieldContainer<double> fc_shards_cfp_deep(shards_array_cfp, true);
279  FieldContainer<double> fc_shards_cfpd_deep(shards_array_cfpd, true);
280  FieldContainer<double> fc_shards_cfpdd_deep(shards_array_cfpdd, true);
281  FieldContainer<double> fc_shards_cfpddd_deep(shards_array_cfpddd, true);
282  // Create FieldContainers using a shallow copy via shards::Array.
283  FieldContainer<double> fc_shards_c_shallow(shards_array_c);
284  FieldContainer<double> fc_shards_cf_shallow(shards_array_cf);
285  FieldContainer<double> fc_shards_cfp_shallow(shards_array_cfp);
286  FieldContainer<double> fc_shards_cfpd_shallow(shards_array_cfpd);
287  FieldContainer<double> fc_shards_cfpdd_shallow(shards_array_cfpdd);
288  FieldContainer<double> fc_shards_cfpddd_shallow(shards_array_cfpddd);
289  } // end variable scope
290 
291 
292  *outStream << "\n" \
293  << "===============================================================================\n"\
294  << "| TEST 2: Usage of resize |\n"\
295  << "===============================================================================\n\n";
296 
297  { // start variable scope
298  // Initialize dimensions for rank-4 multi-index value
299  dimensions.resize(5);
300  dimensions[0] = 5;
301  dimensions[1] = 3;
302  dimensions[2] = 2;
303  dimensions[3] = 7;
304  dimensions[4] = 2;
305 
306  // temporary ints
307  int d0, d1, d2, d3;
308 
309  // Allocate data through a Teuchos::Array, initialized to 0.
310  Teuchos::Array<double> data(dimensions[0]*dimensions[1]*dimensions[2]*dimensions[3]*dimensions[4]);
311 
312  // Create a FieldContainer using a deep copy via Teuchos::Array.
313  FieldContainer<double> fc_array(dimensions, data);
314  // modify the (1,1,1,1) entry
315  double mod_entry = 1.0;
316  fc_array(1,1,1,1,1) = mod_entry;
317  int enumeration = fc_array.getEnumeration(1,1,1,1,1);
318 
319  // Resize into a (dimensions[0]*dimensions[1]) x (dimensions[2]) x (dimensions[3]) x (dimensions[4]) rank-4 array.
320  // This is effectively a reshape, the data should not be touched.
321  fc_array.resize(dimensions[0]*dimensions[1], dimensions[2], dimensions[3], dimensions[4]);
322  // verify that the data array has NOT changed
323  fc_array.getMultiIndex(d0,d1,d2,d3, enumeration);
324  if (std::abs(fc_array(d0,d1,d2,d3) - mod_entry) > zero) {
325  *outStream << "\n\nError in resize.\n\n";
326  errorFlag = -1000;
327  }
328 
329  // Resize into a (dimensions[0]*dimensions[1]*dimensions[2]) x (dimensions[3]) x (dimensions[4]) rank-3 array.
330  // This is effectively a reshape, the data should not be touched.
331  fc_array.resize(dimensions[0]*dimensions[1]*dimensions[2], dimensions[3], dimensions[4]);
332  // verify that the data array has NOT changed
333  fc_array.getMultiIndex(d0,d1,d2, enumeration);
334  if (std::abs(fc_array(d0,d1,d2) - mod_entry) > zero) {
335  *outStream << "\n\nError in resize.\n\n";
336  errorFlag = -1000;
337  }
338 
339  // Resize into a (dimensions[0]*dimensions[1]*dimensions[2]*dimensions[3]) x (dimensions[4]) rank-2 array.
340  // This is effectively a reshape, the data should not be touched.
341  fc_array.resize(dimensions[0]*dimensions[1]*dimensions[2]*dimensions[3], dimensions[4]);
342  // verify that the data array has NOT changed
343  fc_array.getMultiIndex(d0,d1, enumeration);
344  if (std::abs(fc_array(d0,d1) - mod_entry) > zero) {
345  *outStream << "\n\nError in resize.\n\n";
346  errorFlag = -1000;
347  }
348 
349  // Resize into a (dimensions[0]*dimensions[1]*dimensions[2]*dimensions[3]*dimensions[4]) rank-1 array.
350  // This is effectively a reshape, the data should not be touched.
351  fc_array.resize(dimensions[0]*dimensions[1]*dimensions[2]*dimensions[3]*dimensions[4]);
352  // verify that the data array has NOT changed
353  fc_array.getMultiIndex(d0, enumeration);
354  if (std::abs(fc_array(d0) - mod_entry) > zero) {
355  *outStream << "\n\nError in resize.\n\n";
356  errorFlag = -1000;
357  }
358 
359  // More advanced example allocating new memory, with shards::Array.
360  data.assign(dimensions[0]*dimensions[1]*dimensions[2]*dimensions[3]*dimensions[4], 3.0);
361  shards::Array<double,shards::NaturalOrder,Cell,Field,Point,Dim,Dim>
362  shards_array(data.getRawPtr(),dimensions[0],dimensions[1],dimensions[2],dimensions[3],dimensions[4]);
363  // Create a FieldContainer using a shallow copy via shards::Array.
364  FieldContainer<double> fc_shards_shallow(shards_array);
365  fc_shards_shallow.resize(4,4,4,4,4); // keep old entries + allocate new memory and fill with zeros
366  fc_shards_shallow.resize(4*4*4*4*4); // reshape into rank-1 vector
367  if (std::abs(RealSpaceTools<double>::vectorNorm(data.getRawPtr(), fc_array.size(), NORM_ONE) -
368  RealSpaceTools<double>::vectorNorm(fc_shards_shallow, NORM_ONE)) > zero) {
369  *outStream << "\n\nError in resize.\n\n";
370  errorFlag = -1000;
371  }
372 
373 
374  } // end variable scope
375 
376 
377  } // outer try block
378  catch (std::logic_error err) {
379  *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
380  *outStream << err.what() << "\n";
381  *outStream << "-------------------------------------------------------------------------------" << "\n\n";
382  errorFlag = -1000;
383  }
384 
385  if (errorFlag != 0)
386  std::cout << "End Result: TEST FAILED\n";
387  else
388  std::cout << "End Result: TEST PASSED\n";
389 
390  // reset format state of std::cout
391  std::cout.copyfmt(oldFormatState);
392 
393  return errorFlag;
394 }
Implementation of basic linear algebra functionality in Euclidean space.
Header file for utility class to provide multidimensional containers.
Header file for classes providing basic linear algebra functionality in 1D, 2D and 3D...