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 basic linear algebra functionality in Euclidean space.