Intrepid
test_05.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
52#include "Teuchos_oblackholestream.hpp"
53#include "Teuchos_RCP.hpp"
54#include "Teuchos_ScalarTraits.hpp"
55#include "Teuchos_GlobalMPISession.hpp"
56
57using namespace std;
58using namespace Intrepid;
59
60#define INTREPID_TEST_COMMAND( S ) \
61{ \
62 try { \
63 S ; \
64 } \
65 catch (const std::logic_error & err) { \
66 *outStream << "Expected Error ----------------------------------------------------------------\n"; \
67 *outStream << err.what() << '\n'; \
68 *outStream << "-------------------------------------------------------------------------------" << "\n\n"; \
69 }; \
70}
71
72
73int main(int argc, char *argv[]) {
74
75 Teuchos::GlobalMPISession mpiSession(&argc, &argv);
76
77 // This little trick lets us print to std::cout only if
78 // a (dummy) command-line argument is provided.
79 int iprint = argc - 1;
80 Teuchos::RCP<std::ostream> outStream;
81 Teuchos::oblackholestream bhs; // outputs nothing
82 if (iprint > 0)
83 outStream = Teuchos::rcp(&std::cout, false);
84 else
85 outStream = Teuchos::rcp(&bhs, false);
86
87 // Save the format state of the original std::cout.
88 Teuchos::oblackholestream oldFormatState;
89 oldFormatState.copyfmt(std::cout);
90
91 *outStream \
92 << "===============================================================================\n" \
93 << "| |\n" \
94 << "| Unit Test (ArrayTools) |\n" \
95 << "| |\n" \
96 << "| 1) Array operations: clone / scale |\n" \
97 << "| |\n" \
98 << "| Questions? Contact Pavel Bochev (pbboche@sandia.gov) or |\n" \
99 << "| Denis Ridzal (dridzal@sandia.gov). |\n" \
100 << "| |\n" \
101 << "| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
102 << "| Trilinos website: http://trilinos.sandia.gov |\n" \
103 << "| |\n" \
104 << "===============================================================================\n";
105
106
107 int errorFlag = 0;
108#ifdef HAVE_INTREPID_DEBUG
109 int beginThrowNumber = Teuchos::TestForException_getThrowNumber();
110 int endThrowNumber = beginThrowNumber + 21;
111#endif
112
113 typedef ArrayTools art;
114 typedef RealSpaceTools<double> rst;
115#ifdef HAVE_INTREPID_DEBUG
116 ArrayTools atools;
117#endif
118
119 *outStream \
120 << "\n"
121 << "===============================================================================\n"\
122 << "| TEST 1: exceptions |\n"\
123 << "===============================================================================\n";
124
125 try{
126
127#ifdef HAVE_INTREPID_DEBUG
129 FieldContainer<double> a_9_2(9, 2);
130 FieldContainer<double> a_10_2(10, 2);
131 FieldContainer<double> a_10_3(10, 3);
132 FieldContainer<double> a_10_2_2(10, 2, 2);
133 FieldContainer<double> a_10_2_2_2(10, 2, 2, 2);
134 FieldContainer<double> a_10_3_2_2(10, 3, 2, 2);
135 FieldContainer<double> a_10_3_2(10, 3, 2);
136 FieldContainer<double> a_2_2(2, 2);
137 FieldContainer<double> a_2_3_2_2(2, 3, 2, 2);
138 FieldContainer<double> a_2_2_2_2(2, 2, 2, 2);
139 FieldContainer<double> a_10_2_2_2_2(10, 2, 2, 2, 2);
140 FieldContainer<double> a_10_3_2_2_2(10, 3, 2, 2, 2);
141 FieldContainer<double> a_10_2_3_2_2(10, 2, 3, 2, 2);
142 FieldContainer<double> a_10_2_2_3_2(10, 2, 2, 3, 2);
143 FieldContainer<double> a_10_2_2_2_3(10, 2, 2, 2, 3);
144
145 *outStream << "-> cloneFields:\n";
146 INTREPID_TEST_COMMAND( atools.cloneFields<double>(a_10_2_2_2, a_2) );
147 INTREPID_TEST_COMMAND( atools.cloneFields<double>(a_10_2_2_2, a_10_2) );
148 INTREPID_TEST_COMMAND( atools.cloneFields<double>(a_10_3_2, a_2_2) );
149 INTREPID_TEST_COMMAND( atools.cloneFields<double>(a_10_2_2_2_2, a_2_3_2_2) );
150 INTREPID_TEST_COMMAND( atools.cloneFields<double>(a_10_2_2_3_2, a_2_2_2_2) );
151 INTREPID_TEST_COMMAND( atools.cloneFields<double>(a_10_2_2_2_3, a_2_2_2_2) );
152 INTREPID_TEST_COMMAND( atools.cloneFields<double>(a_10_2_2, a_2_2) );
153 INTREPID_TEST_COMMAND( atools.cloneFields<double>(a_10_2_2_2_2, a_2_2_2_2) );
154
155 *outStream << "-> cloneScaleFields:\n";
156 INTREPID_TEST_COMMAND( atools.cloneScaleFields<double>(a_10_2_2_2, a_2, a_2) );
157 INTREPID_TEST_COMMAND( atools.cloneScaleFields<double>(a_10_2_2_2, a_10_2, a_2) );
158 INTREPID_TEST_COMMAND( atools.cloneScaleFields<double>(a_10_2_2_2, a_10_2, a_10_2) );
159 INTREPID_TEST_COMMAND( atools.cloneScaleFields<double>(a_10_2_2, a_9_2, a_10_2) );
160 INTREPID_TEST_COMMAND( atools.cloneScaleFields<double>(a_10_2_2, a_10_3, a_10_2) );
161 INTREPID_TEST_COMMAND( atools.cloneScaleFields<double>(a_10_3_2_2_2, a_10_3, a_2_2_2_2) );
162 INTREPID_TEST_COMMAND( atools.cloneScaleFields<double>(a_10_2_3_2_2, a_10_2, a_2_2_2_2) );
163 INTREPID_TEST_COMMAND( atools.cloneScaleFields<double>(a_10_2_2_3_2, a_10_2, a_2_2_2_2) );
164 INTREPID_TEST_COMMAND( atools.cloneScaleFields<double>(a_10_2_2_2_3, a_10_2, a_2_2_2_2) );
165 INTREPID_TEST_COMMAND( atools.cloneScaleFields<double>(a_10_2_2, a_10_2, a_2_2) );
166 INTREPID_TEST_COMMAND( atools.cloneScaleFields<double>(a_10_2_2_2_2, a_10_2, a_2_2_2_2) );
167
168 *outStream << "-> scaleFields:\n";
169 INTREPID_TEST_COMMAND( atools.scaleFields<double>(a_10_2_2_2, a_2) );
170 INTREPID_TEST_COMMAND( atools.scaleFields<double>(a_10_2, a_2_2) );
171 INTREPID_TEST_COMMAND( atools.scaleFields<double>(a_10_2_2, a_2_2) );
172 INTREPID_TEST_COMMAND( atools.scaleFields<double>(a_10_3_2, a_10_2) );
173 INTREPID_TEST_COMMAND( atools.scaleFields<double>(a_10_3_2_2, a_10_2) );
174 INTREPID_TEST_COMMAND( atools.scaleFields<double>(a_10_3_2_2_2, a_10_2) );
175 INTREPID_TEST_COMMAND( atools.scaleFields<double>(a_10_2_2_2_2, a_10_2) );
176#endif
177
178 }
179 catch (const std::logic_error & err) {
180 *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
181 *outStream << err.what() << '\n';
182 *outStream << "-------------------------------------------------------------------------------" << "\n\n";
183 errorFlag = -1000;
184 };
185
186#ifdef HAVE_INTREPID_DEBUG
187 if (Teuchos::TestForException_getThrowNumber() != endThrowNumber)
188 errorFlag++;
189#endif
190
191 *outStream \
192 << "\n"
193 << "===============================================================================\n"\
194 << "| TEST 2: correctness of math operations |\n"\
195 << "===============================================================================\n";
196
197 outStream->precision(20);
198
199 try {
200 { // start scope
201 *outStream << "\n************ Checking cloneFields ************\n";
202
203 int c=5, p=9, f=7, d1=7, d2=13;
204
205 FieldContainer<double> in_f_p(f, p);
206 FieldContainer<double> in_f_p_d(f, p, d1);
207 FieldContainer<double> in_f_p_d_d(f, p, d1, d2);
208 FieldContainer<double> in_c_f_p(c, f, p);
209 FieldContainer<double> in_c_f_p_d(c, f, p, d1);
210 FieldContainer<double> in_c_f_p_d_d(c, f, p, d1, d2);
211 FieldContainer<double> data_c_p_one(c, p);
212 FieldContainer<double> out_c_f_p(c, f, p);
213 FieldContainer<double> out_c_f_p_d(c, f, p, d1);
214 FieldContainer<double> out_c_f_p_d_d(c, f, p, d1, d2);
215 double zero = INTREPID_TOL*100.0;
216
217 // fill with random numbers
218 for (int i=0; i<in_f_p.size(); i++) {
219 in_f_p[i] = Teuchos::ScalarTraits<double>::random();
220 }
221 for (int i=0; i<in_f_p_d.size(); i++) {
222 in_f_p_d[i] = Teuchos::ScalarTraits<double>::random();
223 }
224 for (int i=0; i<in_f_p_d_d.size(); i++) {
225 in_f_p_d_d[i] = Teuchos::ScalarTraits<double>::random();
226 }
227 for (int i=0; i<data_c_p_one.size(); i++) {
228 data_c_p_one[i] = 1.0;
229 }
230
231 art::cloneFields<double>(out_c_f_p, in_f_p);
232 art::scalarMultiplyDataField<double>(in_c_f_p, data_c_p_one, in_f_p);
233 rst::subtract(&out_c_f_p[0], &in_c_f_p[0], out_c_f_p.size());
234 if (rst::vectorNorm(&out_c_f_p[0], out_c_f_p.size(), NORM_ONE) > zero) {
235 *outStream << "\n\nINCORRECT cloneFields (1): check multiplyScalarData vs. cloneFields\n\n";
236 errorFlag = -1000;
237 }
238 art::cloneFields<double>(out_c_f_p_d, in_f_p_d);
239 art::scalarMultiplyDataField<double>(in_c_f_p_d, data_c_p_one, in_f_p_d);
240 rst::subtract(&out_c_f_p_d[0], &in_c_f_p_d[0], out_c_f_p_d.size());
241 if (rst::vectorNorm(&out_c_f_p_d[0], out_c_f_p_d.size(), NORM_ONE) > zero) {
242 *outStream << "\n\nINCORRECT cloneFields (2): check multiplyScalarData vs. cloneFields\n\n";
243 errorFlag = -1000;
244 }
245 art::cloneFields<double>(out_c_f_p_d_d, in_f_p_d_d);
246 art::scalarMultiplyDataField<double>(in_c_f_p_d_d, data_c_p_one, in_f_p_d_d);
247 rst::subtract(&out_c_f_p_d_d[0], &in_c_f_p_d_d[0], out_c_f_p_d_d.size());
248 if (rst::vectorNorm(&out_c_f_p_d_d[0], out_c_f_p_d_d.size(), NORM_ONE) > zero) {
249 *outStream << "\n\nINCORRECT cloneFields (3): check multiplyScalarData vs. cloneFields\n\n";
250 errorFlag = -1000;
251 }
252 } // end scope
253
254 { // start scope
255 *outStream << "\n************ Checking cloneScaleFields ************\n";
256 int c=5, p=9, f=7, d1=7, d2=13;
257
258 FieldContainer<double> in_f_p(f, p);
259 FieldContainer<double> in_f_p_d(f, p, d1);
260 FieldContainer<double> in_f_p_d_d(f, p, d1, d2);
261 FieldContainer<double> data_c_f(c, f);
262 FieldContainer<double> datainv_c_f(c, f);
263 FieldContainer<double> c_f_p_one(c, f, p);
264 FieldContainer<double> c_f_p_d_one(c, f, p, d1);
265 FieldContainer<double> c_f_p_d_d_one(c, f, p, d1, d2);
266 FieldContainer<double> out_c_f_p(c, f, p);
267 FieldContainer<double> outi_c_f_p(c, f, p);
268 FieldContainer<double> out_c_f_p_d(c, f, p, d1);
269 FieldContainer<double> outi_c_f_p_d(c, f, p, d1);
270 FieldContainer<double> out_c_f_p_d_d(c, f, p, d1, d2);
271 FieldContainer<double> outi_c_f_p_d_d(c, f, p, d1, d2);
272 double zero = INTREPID_TOL*100.0;
273
274 // fill with 1's
275 for (int i=0; i<in_f_p.size(); i++) {
276 in_f_p[i] = 1.0;
277 }
278 for (int i=0; i<in_f_p_d.size(); i++) {
279 in_f_p_d[i] = 1.0;
280 }
281 for (int i=0; i<in_f_p_d_d.size(); i++) {
282 in_f_p_d_d[i] = 1.0;
283 }
284 for (int i=0; i<c_f_p_one.size(); i++) {
285 c_f_p_one[i] = 1.0;
286 }
287 for (int i=0; i<c_f_p_d_one.size(); i++) {
288 c_f_p_d_one[i] = 1.0;
289 }
290 for (int i=0; i<c_f_p_d_d_one.size(); i++) {
291 c_f_p_d_d_one[i] = 1.0;
292 }
293 // fill with random numbers
294 for (int i=0; i<data_c_f.size(); i++) {
295 data_c_f[i] = Teuchos::ScalarTraits<double>::random();
296 datainv_c_f[i] = 1.0 / data_c_f[i];
297 }
298
299 art::cloneScaleFields<double>(out_c_f_p, data_c_f, in_f_p);
300 art::cloneScaleFields<double>(outi_c_f_p, datainv_c_f, in_f_p);
301 for (int i=0; i<out_c_f_p.size(); i++) {
302 out_c_f_p[i] *= outi_c_f_p[i];
303 }
304 rst::subtract(&out_c_f_p[0], &c_f_p_one[0], out_c_f_p.size());
305 if (rst::vectorNorm(&out_c_f_p[0], out_c_f_p.size(), NORM_ONE) > zero) {
306 *outStream << "\n\nINCORRECT cloneScaleValue (1): check scalar inverse property\n\n";
307 errorFlag = -1000;
308 }
309
310 art::cloneScaleFields<double>(out_c_f_p_d, data_c_f, in_f_p_d);
311 art::cloneScaleFields<double>(outi_c_f_p_d, datainv_c_f, in_f_p_d);
312 for (int i=0; i<out_c_f_p_d.size(); i++) {
313 out_c_f_p_d[i] *= outi_c_f_p_d[i];
314 }
315 rst::subtract(&out_c_f_p_d[0], &c_f_p_d_one[0], out_c_f_p_d.size());
316 if (rst::vectorNorm(&out_c_f_p_d[0], out_c_f_p_d.size(), NORM_ONE) > zero) {
317 *outStream << "\n\nINCORRECT cloneScaleValue (2): check scalar inverse property\n\n";
318 errorFlag = -1000;
319 }
320
321 art::cloneScaleFields<double>(out_c_f_p_d_d, data_c_f, in_f_p_d_d);
322 art::cloneScaleFields<double>(outi_c_f_p_d_d, datainv_c_f, in_f_p_d_d);
323 for (int i=0; i<out_c_f_p_d_d.size(); i++) {
324 out_c_f_p_d_d[i] *= outi_c_f_p_d_d[i];
325 }
326 rst::subtract(&out_c_f_p_d_d[0], &c_f_p_d_d_one[0], out_c_f_p_d_d.size());
327 if (rst::vectorNorm(&out_c_f_p_d_d[0], out_c_f_p_d_d.size(), NORM_ONE) > zero) {
328 *outStream << "\n\nINCORRECT cloneScaleValue (3): check scalar inverse property\n\n";
329 errorFlag = -1000;
330 }
331 } // end scope
332
333 { // start scope
334 *outStream << "\n************ Checking scaleFields ************\n";
335 int c=5, p=9, f=7, d1=7, d2=13;
336
337 FieldContainer<double> data_c_f(c, f);
338 FieldContainer<double> datainv_c_f(c, f);
339 FieldContainer<double> out_c_f_p(c, f, p);
340 FieldContainer<double> outi_c_f_p(c, f, p);
341 FieldContainer<double> out_c_f_p_d(c, f, p, d1);
342 FieldContainer<double> outi_c_f_p_d(c, f, p, d1);
343 FieldContainer<double> out_c_f_p_d_d(c, f, p, d1, d2);
344 FieldContainer<double> outi_c_f_p_d_d(c, f, p, d1, d2);
345 double zero = INTREPID_TOL*100.0;
346
347 // fill with random numbers
348 for (int i=0; i<out_c_f_p.size(); i++) {
349 out_c_f_p[i] = Teuchos::ScalarTraits<double>::random();
350 outi_c_f_p[i] = out_c_f_p[i];
351 }
352 for (int i=0; i<out_c_f_p_d.size(); i++) {
353 out_c_f_p_d[i] = Teuchos::ScalarTraits<double>::random();
354 outi_c_f_p_d[i] = out_c_f_p_d[i];
355 }
356 for (int i=0; i<out_c_f_p_d_d.size(); i++) {
357 out_c_f_p_d_d[i] = Teuchos::ScalarTraits<double>::random();
358 outi_c_f_p_d_d[i] = out_c_f_p_d_d[i];
359 }
360 for (int i=0; i<data_c_f.size(); i++) {
361 data_c_f[i] = Teuchos::ScalarTraits<double>::random();
362 datainv_c_f[i] = 1.0 / data_c_f[i];
363 }
364
365 art::scaleFields<double>(out_c_f_p, data_c_f);
366 art::scaleFields<double>(out_c_f_p, datainv_c_f);
367 rst::subtract(&out_c_f_p[0], &outi_c_f_p[0], out_c_f_p.size());
368 if (rst::vectorNorm(&out_c_f_p[0], out_c_f_p.size(), NORM_ONE) > zero) {
369 *outStream << "\n\nINCORRECT scaleValue (1): check scalar inverse property\n\n";
370 errorFlag = -1000;
371 }
372
373 art::scaleFields<double>(out_c_f_p_d, data_c_f);
374 art::scaleFields<double>(out_c_f_p_d, datainv_c_f);
375 rst::subtract(&out_c_f_p_d[0], &outi_c_f_p_d[0], out_c_f_p_d.size());
376 if (rst::vectorNorm(&out_c_f_p_d[0], out_c_f_p_d.size(), NORM_ONE) > zero) {
377 *outStream << "\n\nINCORRECT scaleValue (2): check scalar inverse property\n\n";
378 errorFlag = -1000;
379 }
380
381 art::scaleFields<double>(out_c_f_p_d_d, data_c_f);
382 art::scaleFields<double>(out_c_f_p_d_d, datainv_c_f);
383 rst::subtract(&out_c_f_p_d_d[0], &outi_c_f_p_d_d[0], out_c_f_p_d_d.size());
384 if (rst::vectorNorm(&out_c_f_p_d_d[0], out_c_f_p_d_d.size(), NORM_ONE) > zero) {
385 *outStream << "\n\nINCORRECT cloneScaleValue (3): check scalar inverse property\n\n";
386 errorFlag = -1000;
387 }
388 } // end scope
389
390 /******************************************/
391 *outStream << "\n";
392 }
393 catch (const std::logic_error & err) {
394 *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
395 *outStream << err.what() << '\n';
396 *outStream << "-------------------------------------------------------------------------------" << "\n\n";
397 errorFlag = -1000;
398 };
399
400
401 if (errorFlag != 0)
402 std::cout << "End Result: TEST FAILED\n";
403 else
404 std::cout << "End Result: TEST PASSED\n";
405
406 // reset format state of std::cout
407 std::cout.copyfmt(oldFormatState);
408
409 return errorFlag;
410}
Header file for utility class to provide array tools, such as tensor contractions,...
Header file for utility class to provide multidimensional containers.
Header file for classes providing basic linear algebra functionality in 1D, 2D and 3D.
Utility class that provides methods for higher-order algebraic manipulation of user-defined arrays,...
static void scaleFields(ArrayInOutFields &inoutFields, const ArrayInFactors &inputFactors)
Multiplies, in place, a rank-2, 3, or 4 container with dimensions (C,F,P), (C,F,P,...
static void cloneFields(ArrayOutFields &outputFields, const ArrayInFields &inputFields)
Replicates a rank-2, 3, or 4 container with dimensions (F,P), (F,P,D1) or (F,P,D1,...
static void cloneScaleFields(ArrayOutFields &outputFields, const ArrayInFactors &inputFactors, const ArrayInFields &inputFields)
Multiplies a rank-2, 3, or 4 container with dimensions (F,P), (F,P,D1) or (F,P,D1,...
Implementation of a templated lexicographical container for a multi-indexed scalar quantity....
Implementation of basic linear algebra functionality in Euclidean space.