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 
49 #include "Intrepid_ArrayTools.hpp"
50 #include "Intrepid_FieldContainer.hpp"
51 #include "Intrepid_RealSpaceTools.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 (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 (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 (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 }
Intrepid::ArrayTools
Utility class that provides methods for higher-order algebraic manipulation of user-defined arrays,...
Definition: Intrepid_ArrayTools.hpp:71
Intrepid::ArrayTools::cloneFields
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,...
Definition: Intrepid_ArrayToolsDefCloneScale.hpp:51
Intrepid_FieldContainer.hpp
Header file for utility class to provide multidimensional containers.
Intrepid_RealSpaceTools.hpp
Header file for classes providing basic linear algebra functionality in 1D, 2D and 3D.
Intrepid::RealSpaceTools
Implementation of basic linear algebra functionality in Euclidean space.
Definition: Intrepid_RealSpaceTools.hpp:68
Intrepid::ArrayTools::cloneScaleFields
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,...
Definition: Intrepid_ArrayToolsDefCloneScale.hpp:136
Intrepid::FieldContainer< double >
Intrepid_ArrayTools.hpp
Header file for utility class to provide array tools, such as tensor contractions,...
Intrepid::ArrayTools::scaleFields
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,...
Definition: Intrepid_ArrayToolsDefCloneScale.hpp:228
Generated by   doxygen 1.8.16