Intrepid
test_03.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
44
52#include "Shards_Array.hpp"
53#include "Teuchos_oblackholestream.hpp"
54#include "Teuchos_RCP.hpp"
55#include "Teuchos_GlobalMPISession.hpp"
56
57
58using namespace Intrepid;
59
60SHARDS_ARRAY_DIM_TAG_SIMPLE_DECLARATION( Cell )
61SHARDS_ARRAY_DIM_TAG_SIMPLE_IMPLEMENTATION( Cell )
62
63SHARDS_ARRAY_DIM_TAG_SIMPLE_DECLARATION( Field )
64SHARDS_ARRAY_DIM_TAG_SIMPLE_IMPLEMENTATION( Field )
65
66SHARDS_ARRAY_DIM_TAG_SIMPLE_DECLARATION( Point )
67SHARDS_ARRAY_DIM_TAG_SIMPLE_IMPLEMENTATION( Point )
68
69SHARDS_ARRAY_DIM_TAG_SIMPLE_DECLARATION( Dim )
70SHARDS_ARRAY_DIM_TAG_SIMPLE_IMPLEMENTATION( Dim )
71
72#define INTREPID_TEST_COMMAND( S ) \
73{ \
74 try { \
75 S ; \
76 } \
77 catch (const std::logic_error & err) { \
78 *outStream << "Expected Error ----------------------------------------------------------------\n"; \
79 *outStream << err.what() << '\n'; \
80 *outStream << "-------------------------------------------------------------------------------" << "\n\n"; \
81 }; \
82}
83
84
85int 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 (const 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}
Header file for utility class to provide multidimensional containers.
Header file for classes providing basic linear algebra functionality in 1D, 2D and 3D.
Implementation of a templated lexicographical container for a multi-indexed scalar quantity....
Implementation of basic linear algebra functionality in Euclidean space.