Intrepid
test_05.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 
49 #include "Intrepid_ArrayTools.hpp"
52 #include "Teuchos_oblackholestream.hpp"
53 #include "Teuchos_RCP.hpp"
54 #include "Teuchos_ScalarTraits.hpp"
55 #include "Teuchos_GlobalMPISession.hpp"
56 
57 using namespace std;
58 using 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 
73 int 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
128  FieldContainer<double> a_2(2);
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 }
Implementation of basic linear algebra functionality in Euclidean space.
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,D2), representing the values of a scalar, vector or a tensor field, F-componentwise with a scalar container indexed by (C,F), and stores the result in an output value container of size (C,F,P), (C,F,P,D1) or (C,F,P,D1,D2).
Header file for utility class to provide multidimensional containers.
Header file for utility class to provide array tools, such as tensor contractions, etc.
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 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,D2), representing the values of a scalar, vector or a tensor field, into an output value container of size (C,F,P), (C,F,P,D1) or (C,F,P,D1,D2).
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,D1) or (C,F,P,D1,D2), representing the values of a scalar, vector or a tensor field, F-componentwise with a scalar container indexed by (C,F).