Intrepid
test_04.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 
56 using namespace std;
57 using namespace Intrepid;
58 
59 #define INTREPID_TEST_COMMAND( S ) \
60 { \
61  try { \
62  S ; \
63  } \
64  catch (const std::logic_error & err) { \
65  *outStream << "Expected Error ----------------------------------------------------------------\n"; \
66  *outStream << err.what() << '\n'; \
67  *outStream << "-------------------------------------------------------------------------------" << "\n\n"; \
68  }; \
69 }
70 
71 
72 int main(int argc, char *argv[]) {
73 
74  // This little trick lets us print to std::cout only if
75  // a (dummy) command-line argument is provided.
76  int iprint = argc - 1;
77  Teuchos::RCP<std::ostream> outStream;
78  Teuchos::oblackholestream bhs; // outputs nothing
79  if (iprint > 0)
80  outStream = Teuchos::rcp(&std::cout, false);
81  else
82  outStream = Teuchos::rcp(&bhs, false);
83 
84  // Save the format state of the original std::cout.
85  Teuchos::oblackholestream oldFormatState;
86  oldFormatState.copyfmt(std::cout);
87 
88  *outStream \
89  << "===============================================================================\n" \
90  << "| |\n" \
91  << "| Unit Test (ArrayTools) |\n" \
92  << "| |\n" \
93  << "| 1) Array operations: multiplication, contractions |\n" \
94  << "| |\n" \
95  << "| Questions? Contact Pavel Bochev (pbboche@sandia.gov) or |\n" \
96  << "| Denis Ridzal (dridzal@sandia.gov). |\n" \
97  << "| |\n" \
98  << "| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
99  << "| Trilinos website: http://trilinos.sandia.gov |\n" \
100  << "| |\n" \
101  << "===============================================================================\n";
102 
103 
104  int errorFlag = 0;
105 #ifdef HAVE_INTREPID_DEBUG
106  int beginThrowNumber = Teuchos::TestForException_getThrowNumber();
107  int endThrowNumber = beginThrowNumber + 39 + 18 + 28 + 26 + 34 + 37 + 46 + 45;
108 #endif
109 
110  typedef ArrayTools art;
111  typedef RealSpaceTools<double> rst;
112 
113 
114 #ifdef HAVE_INTREPID_DEBUG
115  ArrayTools atools;
116  /************************************************************************************************
117  * *
118  * Exception tests: should only run when compiled in DEBUG mode *
119  * *
120  ***********************************************************************************************/
121  int C = 12, C1 = 10;
122  int P = 8, P1 = 16;
123  int F = 4, F1 = 8;
124  int D1 = 1, D2 = 2, D3 = 3, D4 = 4;
125 
126  FieldContainer<double> fc_P (P);
127  FieldContainer<double> fc_P_D1 (P, D1);
128  FieldContainer<double> fc_P_D2 (P, D2);
129  FieldContainer<double> fc_P_D3 (P, D3);
130  FieldContainer<double> fc_P_D2_D2 (P, D2, D2);
131  FieldContainer<double> fc_P1_D3 (P1, D3);
132  FieldContainer<double> fc_P1_D2_D2 (P1, D2, D2);
133  FieldContainer<double> fc_P1_D3_D3 (P1, D3, D3);
134 
135  FieldContainer<double> fc_C (C);
136  FieldContainer<double> fc_C1_P (C1,P);
137  FieldContainer<double> fc_C_P1 (C, P1);
138  FieldContainer<double> fc_C_P (C, P);
139  FieldContainer<double> fc_C_P_D1 (C, P, D1);
140  FieldContainer<double> fc_C_P_D2 (C, P, D2);
141  FieldContainer<double> fc_C_P_D3 (C, P, D3);
142  FieldContainer<double> fc_C_P_D4 (C, P, D4);
143  FieldContainer<double> fc_C1_P_D2 (C1,P, D2);
144  FieldContainer<double> fc_C1_P_D3 (C1,P, D3);
145  FieldContainer<double> fc_C_P1_D1 (C, P1,D1);
146  FieldContainer<double> fc_C_P1_D2 (C, P1,D2);
147  FieldContainer<double> fc_C_P1_D3 (C, P1,D3);
148  FieldContainer<double> fc_C_P_D1_D1 (C, P, D1, D1);
149  FieldContainer<double> fc_C_P_D1_D2 (C, P, D1, D2);
150  FieldContainer<double> fc_C_P_D1_D3 (C, P, D1, D3);
151  FieldContainer<double> fc_C_P_D2_D2 (C, P, D2, D2);
152  FieldContainer<double> fc_C_P_D3_D1 (C, P, D3, D1);
153  FieldContainer<double> fc_C_P_D3_D2 (C, P, D3, D2);
154  FieldContainer<double> fc_C_P_D2_D3 (C, P, D2, D3);
155  FieldContainer<double> fc_C_P_D3_D3 (C, P, D3, D3);
156  FieldContainer<double> fc_C1_P_D3_D3 (C1,P, D3, D3);
157  FieldContainer<double> fc_C1_P_D2_D2 (C1,P, D2, D2);
158  FieldContainer<double> fc_C_P1_D2_D2 (C, P1,D2, D2);
159  FieldContainer<double> fc_C_P1_D3_D3 (C, P1,D3, D3);
160  FieldContainer<double> fc_C_P_D3_D3_D3 (C, P, D3, D3, D3);
161 
162  FieldContainer<double> fc_F_P (F, P);
163  FieldContainer<double> fc_F_P_D1 (F, P, D1);
164  FieldContainer<double> fc_F_P_D2 (F, P, D2);
165  FieldContainer<double> fc_F_P1_D2 (F, P1, D2);
166  FieldContainer<double> fc_F_P_D3 (F, P, D3);
167  FieldContainer<double> fc_F_P_D3_D3 (F, P, D3, D3);
168  FieldContainer<double> fc_F1_P_D2 (F1,P, D2);
169  FieldContainer<double> fc_F1_P_D3 (F1,P, D3);
170  FieldContainer<double> fc_F1_P_D3_D3 (F1,P, D3, D3);
171  FieldContainer<double> fc_F_P1_D3 (F, P1,D3);
172  FieldContainer<double> fc_F_P1_D3_D3 (F, P1,D3, D3);
173  FieldContainer<double> fc_F_P_D2_D2 (F, P, D2, D2);
174  FieldContainer<double> fc_F_P1_D2_D2 (F, P1,D2, D2);
175  FieldContainer<double> fc_C_F_P (C, F, P);
176  FieldContainer<double> fc_C1_F_P (C1, F, P);
177  FieldContainer<double> fc_C_F1_P (C, F1,P);
178  FieldContainer<double> fc_C_F_P1 (C, F, P1);
179  FieldContainer<double> fc_C_F_P_D1 (C, F, P, D1);
180  FieldContainer<double> fc_C_F_P_D2 (C, F, P, D2);
181  FieldContainer<double> fc_C_F_P_D3 (C, F, P, D3);
182  FieldContainer<double> fc_C1_F_P_D2 (C1, F, P,D2);
183  FieldContainer<double> fc_C1_F_P_D3 (C1, F, P,D3);
184  FieldContainer<double> fc_C_F1_P_D2 (C, F1,P, D2);
185  FieldContainer<double> fc_C_F1_P_D3 (C, F1,P, D3);
186  FieldContainer<double> fc_C_F1_P_D3_D3 (C, F1,P, D3, D3);
187  FieldContainer<double> fc_C_F_P1_D2 (C, F, P1,D2);
188  FieldContainer<double> fc_C_F_P1_D3 (C, F, P1,D3);
189  FieldContainer<double> fc_C_F_P_D1_D1 (C, F, P, D1, D1);
190  FieldContainer<double> fc_C_F_P_D2_D2 (C, F, P, D2, D2);
191  FieldContainer<double> fc_C_F_P_D1_D3 (C, F, P, D1, D3);
192  FieldContainer<double> fc_C_F_P_D2_D3 (C, F, P, D2, D3);
193  FieldContainer<double> fc_C_F_P_D3_D1 (C, F, P, D3, D1);
194  FieldContainer<double> fc_C_F_P_D3_D2 (C, F, P, D3, D2);
195  FieldContainer<double> fc_C_F_P_D3_D3 (C, F, P, D3, D3);
196  FieldContainer<double> fc_C_F_P1_D2_D2 (C, F, P1,D2, D2);
197  FieldContainer<double> fc_C_F_P1_D3_D3 (C, F, P1,D3, D3);
198  FieldContainer<double> fc_C1_F_P_D2_D2 (C1,F, P, D2, D2);
199  FieldContainer<double> fc_C1_F_P1_D2_D2 (C1,F, P1,D2, D2);
200  FieldContainer<double> fc_C1_F_P_D3_D3 (C1,F, P, D3, D3);
201 
202  Teuchos::Array<int> dimensions(6);
203  dimensions[0] = C;
204  dimensions[1] = F;
205  dimensions[2] = P;
206  dimensions[3] = D3;
207  dimensions[4] = D3;
208  dimensions[5] = D3;
209  FieldContainer<double> fc_C_F_P_D3_D3_D3 (dimensions);
210 
211 
212  *outStream \
213  << "\n"
214  << "===============================================================================\n"\
215  << "| TEST 1: crossProductDataField exceptions |\n"\
216  << "===============================================================================\n";
217 
218  try{
219  // 39 exceptions
220  // Test rank and D dimension of inputData
221  INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C_F_P_D3, fc_C_P, fc_C_F_P_D3) );
222  INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C_F_P_D3, fc_C_P_D2_D2, fc_C_F_P_D3) );
223  INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C_F_P_D3, fc_C_P_D1, fc_C_F_P_D3) );
224 
225  // Test rank and D dimension of inputFields
226  INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C_F_P_D3, fc_C_P_D3, fc_F_P) );
227  INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C_F_P_D3, fc_C_P_D3, fc_C_F_P_D3_D3) );
228  INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C_F_P_D3, fc_C_P_D3, fc_C_F_P_D1) );
229  INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C_F_P_D3, fc_C_P_D3, fc_F_P_D1) );
230 
231  // Test rank of outputFields
232  INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C_F_P_D3, fc_C_P_D2, fc_C_F_P_D2) );
233  INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C_F_P, fc_C_P_D3, fc_C_F_P_D3) );
234 
235  // Dimension cross-check: (1) inputData vs. inputFields
236  INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C_F_P_D3, fc_C_P_D3, fc_C1_F_P_D3) );
237  INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C_F_P_D3, fc_C_P_D3, fc_C_F_P1_D3) );
238  INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C_F_P_D3, fc_C_P_D3, fc_C_F_P_D2) );
239  INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C_F_P, fc_C_P_D2, fc_C1_F_P_D2) );
240  INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C_F_P, fc_C_P_D2, fc_C_F_P1_D2) );
241  INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C_F_P, fc_C_P_D2, fc_C_F_P_D3) );
242  INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C_F_P_D3, fc_C_P_D3, fc_F_P1_D3) );
243  INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C_F_P_D3, fc_C_P_D3, fc_F_P_D2) );
244  INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C_F_P, fc_C_P_D2, fc_F_P1_D2) );
245  INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C_F_P, fc_C_P_D2, fc_F_P_D3) );
246 
247  // Dimension cross-check: (2) outputFields vs. inputData
248  INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C1_F_P, fc_C_P_D2, fc_C_F_P_D2) );
249  INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C_F_P1, fc_C_P_D2, fc_C_F_P_D2) );
250  INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C1_F_P, fc_C_P_D2, fc_F_P_D2) );
251  INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C_F_P1, fc_C_P_D2, fc_F_P_D2) );
252  INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C1_F_P_D3, fc_C_P_D3, fc_C_F_P_D3) );
253  INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C_F_P1_D3, fc_C_P_D3, fc_C_F_P_D3) );
254  INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C_F_P_D2, fc_C_P_D3, fc_C_F_P_D3) );
255  INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C1_F_P_D3, fc_C_P_D3, fc_F_P_D3) );
256  INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C_F_P1_D3, fc_C_P_D3, fc_F_P_D3) );
257  INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C_F_P_D2, fc_C_P_D3, fc_F_P_D3) );
258 
259  // Dimension cross-check: (3) outputFields vs. inputFields
260  INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C_F_P, fc_C1_P_D2, fc_C1_F_P_D2) );
261  INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C_F_P, fc_C_P1_D2, fc_C_F_P1_D2) );
262  INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C_F_P, fc_C_P_D2, fc_C_F1_P_D2) );
263  INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C_F_P, fc_C_P1_D2, fc_F_P1_D2) );
264  INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C_F_P, fc_C_P_D2, fc_F1_P_D2) );
265  INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C_F_P_D3, fc_C1_P_D3, fc_C1_F_P_D3) );
266  INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C_F_P_D3, fc_C_P1_D3, fc_C_F_P1_D3) );
267  INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C_F_P_D3, fc_C_P_D3, fc_C_F1_P_D3) );
268  INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C_F_P_D3, fc_C_P1_D3, fc_F_P1_D3) );
269  INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C_F_P_D3, fc_C_P_D3, fc_F1_P_D3) );
270 
271  *outStream \
272  << "\n"
273  << "===============================================================================\n"\
274  << "| TEST 2: crossProductDataData exceptions |\n"\
275  << "===============================================================================\n";
276 
277  // 18 exceptions
278  // inputDataL is (C, P, D) and 2 <= D <= 3 is required
279  INTREPID_TEST_COMMAND(atools.crossProductDataData<double>(fc_C_P_D3, fc_C_P, fc_C_P_D3) );
280  INTREPID_TEST_COMMAND(atools.crossProductDataData<double>(fc_C_P_D3, fc_C_P_D2_D2, fc_C_P_D3) );
281  INTREPID_TEST_COMMAND(atools.crossProductDataData<double>(fc_C_P_D3, fc_C_P_D1, fc_C_P_D3) );
282 
283  // inputDataRight is (C, P, D) or (P, D) and 2 <= (D=dimension(rank - 1)) <= 3 is required.
284  INTREPID_TEST_COMMAND(atools.crossProductDataData<double>(fc_C_P_D3, fc_C_P_D3, fc_C) );
285  INTREPID_TEST_COMMAND(atools.crossProductDataData<double>(fc_C_P_D3, fc_C_P_D3, fc_C_P_D2_D2) );
286  INTREPID_TEST_COMMAND(atools.crossProductDataData<double>(fc_C_P_D3, fc_C_P_D3, fc_C_P_D1) );
287 
288  // outputData is (C,P,D) in 3D and (C,P) in 2D => rank = inputDataLeft.dimension(2)
289  INTREPID_TEST_COMMAND(atools.crossProductDataData<double>(fc_C_P_D3, fc_C_P_D2, fc_C_P_D2) );
290  INTREPID_TEST_COMMAND(atools.crossProductDataData<double>(fc_C_P, fc_C_P_D3, fc_C_P_D3) );
291 
292  // Dimension cross-check (1):
293  // inputDataLeft(C,P,D) vs. inputDataRight(C,P,D): C, P, D must match
294  INTREPID_TEST_COMMAND(atools.crossProductDataData<double>(fc_C_P_D3, fc_C1_P_D3, fc_C_P_D3) );
295  INTREPID_TEST_COMMAND(atools.crossProductDataData<double>(fc_C_P_D3, fc_C_P1_D3, fc_C_P_D3) );
296  INTREPID_TEST_COMMAND(atools.crossProductDataData<double>(fc_C_P_D3, fc_C_P_D3, fc_C_P_D2) );
297  // inputDataLeft(C, P,D) vs. inputDataRight(P,D): P, D must match
298  INTREPID_TEST_COMMAND(atools.crossProductDataData<double>(fc_C_P_D3, fc_C_P1_D3, fc_P_D3) );
299  INTREPID_TEST_COMMAND(atools.crossProductDataData<double>(fc_C_P_D3, fc_C_P_D3, fc_P_D2) );
300 
301  // Dimension cross-check (2):
302  // in 2D: outputData(C,P) vs. inputDataLeft(C,P,D): dimensions C, P must match
303  INTREPID_TEST_COMMAND(atools.crossProductDataData<double>(fc_C1_P, fc_C_P_D2, fc_P_D2) );
304  INTREPID_TEST_COMMAND(atools.crossProductDataData<double>(fc_C_P1, fc_C_P_D2, fc_P_D2) );
305  // in 3D: outputData(C,P,D) vs. inputDataLeft(C,P,D): all dimensions C, P, D must match
306  INTREPID_TEST_COMMAND(atools.crossProductDataData<double>(fc_C1_P_D3, fc_C_P_D3, fc_P_D3) );
307  INTREPID_TEST_COMMAND(atools.crossProductDataData<double>(fc_C_P1_D3, fc_C_P_D3, fc_P_D3) );
308  INTREPID_TEST_COMMAND(atools.crossProductDataData<double>(fc_C_P_D2, fc_C_P_D3, fc_P_D3) );
309 
310  *outStream \
311  << "\n"
312  << "===============================================================================\n"\
313  << "| TEST 3: outerProductDataField exceptions |\n"\
314  << "===============================================================================\n";
315  // 28 exceptions
316  // Test rank and D dimension: inputData(C, P, D) and 2 <= D <= 3 is required
317  INTREPID_TEST_COMMAND(atools.outerProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P, fc_C_F_P_D3) );
318  INTREPID_TEST_COMMAND(atools.outerProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D2_D2, fc_C_F_P_D3) );
319  INTREPID_TEST_COMMAND(atools.outerProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D1, fc_C_F_P_D3) );
320 
321  // Test rank and D dimension: inputFields(C,F,P,D)/(F,P,D) and 2 <= (D=dimension(rank - 1)) <= 3 is required.
322  INTREPID_TEST_COMMAND(atools.outerProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3, fc_F_P) );
323  INTREPID_TEST_COMMAND(atools.outerProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3, fc_C_F_P_D3_D3) );
324  INTREPID_TEST_COMMAND(atools.outerProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3, fc_C_F_P_D1) );
325 
326  // Test rank and D dimension: outputFields(C,F,P,D,D)
327  INTREPID_TEST_COMMAND(atools.outerProductDataField<double>(fc_C_F_P, fc_C_P_D3, fc_C_F_P_D3) );
328  INTREPID_TEST_COMMAND(atools.outerProductDataField<double>(fc_C_F_P_D3, fc_C_P_D3, fc_C_F_P_D3) );
329  INTREPID_TEST_COMMAND(atools.outerProductDataField<double>(fc_C_F_P_D3_D3_D3,fc_C_P_D3, fc_C_F_P_D3) );
330  INTREPID_TEST_COMMAND(atools.outerProductDataField<double>(fc_C_F_P_D1_D3, fc_C_P_D3, fc_C_F_P_D3) );
331  INTREPID_TEST_COMMAND(atools.outerProductDataField<double>(fc_C_F_P_D3_D1, fc_C_P_D3, fc_C_F_P_D3) );
332  INTREPID_TEST_COMMAND(atools.outerProductDataField<double>(fc_C_F_P_D1_D1, fc_C_P_D3, fc_C_F_P_D3) );
333 
334  // Cross-check (2): outputFields(C,F,P,D,D) vs. inputData(C,P,D): dimensions C, P, D must match
335  INTREPID_TEST_COMMAND(atools.outerProductDataField<double>(fc_C_F_P_D3_D3, fc_C1_P_D3, fc_C_F_P_D3) );
336  INTREPID_TEST_COMMAND(atools.outerProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P1_D3, fc_C_F_P_D3) );
337  INTREPID_TEST_COMMAND(atools.outerProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D2, fc_C_F_P_D2) );
338 
339  // Cross-check (1): inputData(C,P,D) vs. inputFields(C,F,P,D): dimensions C, P, D must match
340  INTREPID_TEST_COMMAND(atools.outerProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3, fc_C1_F_P_D3) );
341  INTREPID_TEST_COMMAND(atools.outerProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3, fc_C_F_P1_D3) );
342  INTREPID_TEST_COMMAND(atools.outerProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3, fc_C_F_P_D2) );
343  // Cross-check (1): inputData(C,P,D) vs. inputFields(F,P,D): dimensions P, D must match
344  INTREPID_TEST_COMMAND(atools.outerProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3, fc_F_P1_D3) );
345  INTREPID_TEST_COMMAND(atools.outerProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3, fc_F_P_D2) );
346 
347  // Cross-check (3): outputFields(C,F,P,D,D) vs. inputFields(C,F,P,D): dimensions C, F, P, D must match
348  INTREPID_TEST_COMMAND(atools.outerProductDataField<double>(fc_C_F_P_D3_D3, fc_C1_P_D3, fc_C1_F_P_D3) );
349  INTREPID_TEST_COMMAND(atools.outerProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3, fc_C_F1_P_D3) );
350  INTREPID_TEST_COMMAND(atools.outerProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P1_D3, fc_C_F_P1_D3) );
351  INTREPID_TEST_COMMAND(atools.outerProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D2, fc_C_F_P_D2) );
352  INTREPID_TEST_COMMAND(atools.outerProductDataField<double>(fc_C_F_P_D2_D3, fc_C_P_D2, fc_C_F_P_D2) );
353  // Cross-check (3): outputFields(C,F,P,D,D) vs. inputFields(F,P,D): dimensions F, P, D must match
354  INTREPID_TEST_COMMAND(atools.outerProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3, fc_F1_P_D3) );
355  INTREPID_TEST_COMMAND(atools.outerProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P1_D3, fc_F_P1_D3) );
356  INTREPID_TEST_COMMAND(atools.outerProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D2, fc_F_P_D2) );
357  *outStream \
358  << "\n"
359  << "===============================================================================\n"\
360  << "| TEST 4: outerProductDataData exceptions |\n"\
361  << "===============================================================================\n";
362  // 26 exceptions
363  // (1) inputDataLeft is (C, P, D) and 2 <= D <= 3 is required
364  INTREPID_TEST_COMMAND(atools.outerProductDataData<double>(fc_C_P_D3_D3, fc_C_P, fc_C_P_D3) );
365  INTREPID_TEST_COMMAND(atools.outerProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D3_D3, fc_C_P_D3) );
366  INTREPID_TEST_COMMAND(atools.outerProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D1, fc_C_P_D3) );
367 
368  // (2) inputDataRight is (C, P, D) or (P, D) and 2 <= (D=dimension(rank - 1)) <= 3 is required.
369  INTREPID_TEST_COMMAND(atools.outerProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D3, fc_C) );
370  INTREPID_TEST_COMMAND(atools.outerProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D3, fc_C_P_D3_D3) );
371  INTREPID_TEST_COMMAND(atools.outerProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D3, fc_C_P_D1) );
372  INTREPID_TEST_COMMAND(atools.outerProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D3, fc_P_D1) );
373 
374  // (3) outputData is (C,P,D,D)
375  INTREPID_TEST_COMMAND(atools.outerProductDataData<double>(fc_C_P_D3, fc_C_P_D3, fc_C_P_D3) );
376  INTREPID_TEST_COMMAND(atools.outerProductDataData<double>(fc_C_P_D3_D3_D3,fc_C_P_D3, fc_C_P_D3) );
377  INTREPID_TEST_COMMAND(atools.outerProductDataData<double>(fc_C_P_D3_D1, fc_C_P_D3, fc_C_P_D3) );
378  INTREPID_TEST_COMMAND(atools.outerProductDataData<double>(fc_C_P_D1_D2, fc_C_P_D3, fc_C_P_D3) );
379 
380  // Cross-check (2): outputData(C,P,D,D) vs. inputDataLeft(C,P,D): dimensions C, P, D must match
381  INTREPID_TEST_COMMAND(atools.outerProductDataData<double>(fc_C_P_D3_D3, fc_C1_P_D3, fc_P_D3) );
382  INTREPID_TEST_COMMAND(atools.outerProductDataData<double>(fc_C_P_D3_D3, fc_C_P1_D3, fc_P_D3) );
383  INTREPID_TEST_COMMAND(atools.outerProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D2, fc_P_D3) );
384  INTREPID_TEST_COMMAND(atools.outerProductDataData<double>(fc_C_P_D2_D3, fc_C_P_D2, fc_P_D3) );
385  INTREPID_TEST_COMMAND(atools.outerProductDataData<double>(fc_C_P_D3_D2, fc_C_P_D2, fc_P_D3) );
386 
387  // Cross-check (1): inputDataLeft(C,P,D) vs. inputDataRight(C,P,D): all dimensions C, P, D must match
388  INTREPID_TEST_COMMAND(atools.outerProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D3, fc_C1_P_D3) );
389  INTREPID_TEST_COMMAND(atools.outerProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D3, fc_C_P1_D3) );
390  INTREPID_TEST_COMMAND(atools.outerProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D3, fc_C_P_D2) );
391  // Cross-check (1): inputDataLeft(C,P,D) vs. inputDataRight(P,D): dimensions P, D must match
392  INTREPID_TEST_COMMAND(atools.outerProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D3, fc_P1_D3) );
393  INTREPID_TEST_COMMAND(atools.outerProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D3, fc_P_D2) );
394 
395  // Cross-check (3): outputData(C,P,D,D) vs. inputDataRight(C,P,D): dimensions C, P, D must match
396  INTREPID_TEST_COMMAND(atools.outerProductDataData<double>(fc_C_P_D3_D3, fc_C1_P_D3, fc_C1_P_D3) );
397  INTREPID_TEST_COMMAND(atools.outerProductDataData<double>(fc_C_P_D3_D3, fc_C_P1_D3, fc_C_P1_D3) );
398  INTREPID_TEST_COMMAND(atools.outerProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D2, fc_C_P_D2) );
399  // Cross-check (3): outputData(C,P,D,D) vs. inputDataRight(P,D): dimensions P, D must match
400  INTREPID_TEST_COMMAND(atools.outerProductDataData<double>(fc_C_P_D3_D3, fc_C_P1_D3, fc_P1_D3) );
401  INTREPID_TEST_COMMAND(atools.outerProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D2, fc_P_D2) );
402 
403  *outStream \
404  << "\n"
405  << "===============================================================================\n"\
406  << "| TEST 5: matvecProductDataField exceptions |\n"\
407  << "===============================================================================\n";
408  // 34 exceptions
409  // (1) inputData is (C,P), (C,P,D) or (C, P, D, D) and 1 <= D <= 3 is required
410  INTREPID_TEST_COMMAND(atools.matvecProductDataField<double>(fc_C_F_P_D3, fc_C, fc_C_F_P_D3) );
411  INTREPID_TEST_COMMAND(atools.matvecProductDataField<double>(fc_C_F_P_D3, fc_C_P_D4, fc_C_F_P_D3) );
412  INTREPID_TEST_COMMAND(atools.matvecProductDataField<double>(fc_C_F_P_D3, fc_C_P_D3_D3_D3, fc_C_F_P_D3) );
413  INTREPID_TEST_COMMAND(atools.matvecProductDataField<double>(fc_C_F_P_D3, fc_C_P_D3_D1, fc_C_F_P_D3) );
414  INTREPID_TEST_COMMAND(atools.matvecProductDataField<double>(fc_C_F_P_D3, fc_C_P_D1_D3, fc_C_F_P_D3) );
415 
416  // (2) inputFields is (C, F, P, D) or (F, P, D) and 1 <= (D=dimension(rank - 1)) <= 3 is required.
417  INTREPID_TEST_COMMAND(atools.matvecProductDataField<double>(fc_C_F_P_D3, fc_C_P_D3_D3, fc_P_D3) );
418  INTREPID_TEST_COMMAND(atools.matvecProductDataField<double>(fc_C_F_P_D3, fc_C_P_D3_D3, fc_C_F_P_D3_D3) );
419  INTREPID_TEST_COMMAND(atools.matvecProductDataField<double>(fc_C_F_P_D3, fc_C_P_D3_D3, fc_C_F_P_D1) );
420  INTREPID_TEST_COMMAND(atools.matvecProductDataField<double>(fc_C_F_P_D3, fc_C_P_D3_D3, fc_F_P_D1) );
421  // (3) outputFields is (C,F,P,D)
422  INTREPID_TEST_COMMAND(atools.matvecProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3_D3, fc_C_F_P_D3) );
423  INTREPID_TEST_COMMAND(atools.matvecProductDataField<double>(fc_C_F_P, fc_C_P_D3_D3, fc_C_F_P_D3) );
424  INTREPID_TEST_COMMAND(atools.matvecProductDataField<double>(fc_C_F_P_D1, fc_C_P_D3_D3, fc_C_F_P_D3) );
425 
426  // Cross-check (2): outputFields(C,F,P,D) vs. inputData(C,P,D) and (C,P,D,D): dimensions C, P, D must match
427  INTREPID_TEST_COMMAND(atools.matvecProductDataField<double>(fc_C_F_P_D3, fc_C1_P_D3_D3, fc_C_F_P_D3) );
428  INTREPID_TEST_COMMAND(atools.matvecProductDataField<double>(fc_C_F_P_D3, fc_C_P1_D3_D3, fc_C_F_P_D3) );
429  INTREPID_TEST_COMMAND(atools.matvecProductDataField<double>(fc_C_F_P_D3, fc_C_P_D2_D2, fc_C_F_P_D3) );
430  INTREPID_TEST_COMMAND(atools.matvecProductDataField<double>(fc_C_F_P_D3, fc_C1_P, fc_C_F_P_D3) );
431  INTREPID_TEST_COMMAND(atools.matvecProductDataField<double>(fc_C_F_P_D3, fc_C_P1, fc_C_F_P_D3) );
432  INTREPID_TEST_COMMAND(atools.matvecProductDataField<double>(fc_C_F_P_D3, fc_C1_P_D3, fc_C_F_P_D3) );
433  INTREPID_TEST_COMMAND(atools.matvecProductDataField<double>(fc_C_F_P_D3, fc_C_P1_D3, fc_C_F_P_D3) );
434  INTREPID_TEST_COMMAND(atools.matvecProductDataField<double>(fc_C_F_P_D3, fc_C_P_D2, fc_C_F_P_D3) );
435 
436  // Cross-check (1): inputData(C,P), (C,P,D) or (C,P,D,D) vs. inputFields(C,F,P,D): dimensions C, P, D must match
437  INTREPID_TEST_COMMAND(atools.matvecProductDataField<double>(fc_C_F_P_D3, fc_C_P_D3_D3, fc_C1_F_P_D3) );
438  INTREPID_TEST_COMMAND(atools.matvecProductDataField<double>(fc_C_F_P_D3, fc_C_P_D3_D3, fc_C_F_P1_D3) );
439  INTREPID_TEST_COMMAND(atools.matvecProductDataField<double>(fc_C_F_P_D3, fc_C_P_D3_D3, fc_C_F_P_D2) );
440  INTREPID_TEST_COMMAND(atools.matvecProductDataField<double>(fc_C_F_P_D3, fc_C_P_D3, fc_C_F_P_D2) );
441 
442  // Cross-check (1): inputData(C,P), (C,P,D) or (C,P,D,D) vs. inputFields(F,P,D): dimensions P, D must match
443  INTREPID_TEST_COMMAND(atools.matvecProductDataField<double>(fc_C_F_P_D3, fc_C_P_D3_D3, fc_F_P1_D3) );
444  INTREPID_TEST_COMMAND(atools.matvecProductDataField<double>(fc_C_F_P_D3, fc_C_P_D3_D3, fc_F_P_D2) );
445  INTREPID_TEST_COMMAND(atools.matvecProductDataField<double>(fc_C_F_P_D3, fc_C_P_D3, fc_F_P_D2) );
446 
447  // Cross-check (3): outputFields(C,F,P,D) vs. inputFields(C,F,P,D): all dimensions C, F, P, D must match
448  INTREPID_TEST_COMMAND(atools.matvecProductDataField<double>(fc_C_F_P_D3, fc_C_P_D3_D3, fc_C1_F_P_D3) );
449  INTREPID_TEST_COMMAND(atools.matvecProductDataField<double>(fc_C_F_P_D3, fc_C_P_D3_D3, fc_C_F1_P_D3) );
450  INTREPID_TEST_COMMAND(atools.matvecProductDataField<double>(fc_C_F_P_D3, fc_C_P1_D3_D3, fc_C_F_P1_D3) );
451  INTREPID_TEST_COMMAND(atools.matvecProductDataField<double>(fc_C_F_P_D2, fc_C_P_D2_D2, fc_C_F_P_D3) );
452 
453  // Cross-check (3): outputFields(C,F,P,D) vs. inputFields(F,P,D): dimensions F, P, D must match
454  INTREPID_TEST_COMMAND(atools.matvecProductDataField<double>(fc_C_F_P_D3, fc_C_P_D3_D3, fc_C_F1_P_D3) );
455  INTREPID_TEST_COMMAND(atools.matvecProductDataField<double>(fc_C_F_P_D3, fc_C_P1_D3_D3, fc_C_F_P1_D3) );
456  INTREPID_TEST_COMMAND(atools.matvecProductDataField<double>(fc_C_F_P_D3, fc_C_P_D3_D3, fc_C_F_P_D2) );
457 
458  *outStream \
459  << "\n"
460  << "===============================================================================\n"\
461  << "| TEST 6: matvecProductDataData exceptions |\n"\
462  << "===============================================================================\n";
463  // 37 exceptions
464  // (1) inputDataLeft is (C,P), (C,P,D) or (C,P,D,D) and 1 <= D <= 3 is required
465  INTREPID_TEST_COMMAND(atools.matvecProductDataData<double>(fc_C_P_D3, fc_C, fc_C_P_D3) );
466  INTREPID_TEST_COMMAND(atools.matvecProductDataData<double>(fc_C_P_D3, fc_C_P_D3_D3_D3, fc_C_P_D3) );
467  INTREPID_TEST_COMMAND(atools.matvecProductDataData<double>(fc_C_P_D3, fc_C_P_D3_D1, fc_C_P_D3) );
468  INTREPID_TEST_COMMAND(atools.matvecProductDataData<double>(fc_C_P_D3, fc_C_P_D1_D3, fc_C_P_D3) );
469 
470  // (2) inputDataRight is (C, P, D) or (P, D) and 1 <= (D=dimension(rank - 1)) <= 3 is required.
471  INTREPID_TEST_COMMAND(atools.matvecProductDataData<double>(fc_C_P_D3, fc_C_P_D3_D3, fc_C_P_D3_D3) );
472  INTREPID_TEST_COMMAND(atools.matvecProductDataData<double>(fc_C_P_D3, fc_C_P_D3_D3, fc_P) );
473  INTREPID_TEST_COMMAND(atools.matvecProductDataData<double>(fc_C_P_D3, fc_C_P_D3_D3, fc_P_D1) );
474  INTREPID_TEST_COMMAND(atools.matvecProductDataData<double>(fc_C_P_D3, fc_C_P_D3_D3, fc_C_P_D1) );
475 
476  // (3) outputData is (C,P,D)
477  INTREPID_TEST_COMMAND(atools.matvecProductDataData<double>(fc_C_P, fc_C_P_D3_D3, fc_C_P_D3) );
478  INTREPID_TEST_COMMAND(atools.matvecProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D3_D3, fc_C_P_D3) );
479  INTREPID_TEST_COMMAND(atools.matvecProductDataData<double>(fc_C_P_D1, fc_C_P_D3_D3, fc_C_P_D3) );
480 
481  // Cross-check (2): outputData(C,P,D) vs. inputDataLeft(C,P), (C,P,D), (C,P,D,D): dimensions C, P, D must match
482  INTREPID_TEST_COMMAND(atools.matvecProductDataData<double>(fc_C1_P_D3, fc_C_P_D3_D3, fc_C_P_D3) );
483  INTREPID_TEST_COMMAND(atools.matvecProductDataData<double>(fc_C_P1_D3, fc_C_P_D3_D3, fc_C_P_D3) );
484  INTREPID_TEST_COMMAND(atools.matvecProductDataData<double>(fc_C_P_D2, fc_C_P_D3_D3, fc_C_P_D3) );
485  INTREPID_TEST_COMMAND(atools.matvecProductDataData<double>(fc_C1_P_D3, fc_C_P, fc_C_P_D3) );
486  INTREPID_TEST_COMMAND(atools.matvecProductDataData<double>(fc_C_P1_D3, fc_C_P, fc_C_P_D3) );
487  INTREPID_TEST_COMMAND(atools.matvecProductDataData<double>(fc_C1_P_D3, fc_C_P_D3, fc_C_P_D3) );
488  INTREPID_TEST_COMMAND(atools.matvecProductDataData<double>(fc_C_P1_D3, fc_C_P_D3, fc_C_P_D3) );
489  INTREPID_TEST_COMMAND(atools.matvecProductDataData<double>(fc_C_P_D2, fc_C_P_D3, fc_C_P_D3) );
490 
491  // Cross-check (1): inputDataLeft(C,P), (C,P,D), or (C,P,D,D) vs. inputDataRight(C,P,D): dimensions C, P, D must match
492  INTREPID_TEST_COMMAND(atools.matvecProductDataData<double>(fc_C_P_D3, fc_C_P_D3_D3, fc_C1_P_D3) );
493  INTREPID_TEST_COMMAND(atools.matvecProductDataData<double>(fc_C_P_D3, fc_C_P_D3_D3, fc_C_P1_D3) );
494  INTREPID_TEST_COMMAND(atools.matvecProductDataData<double>(fc_C_P_D3, fc_C_P_D3_D3, fc_C_P_D2) );
495  INTREPID_TEST_COMMAND(atools.matvecProductDataData<double>(fc_C_P_D3, fc_C_P, fc_C1_P_D3) );
496  INTREPID_TEST_COMMAND(atools.matvecProductDataData<double>(fc_C_P_D3, fc_C_P, fc_C_P1_D3) );
497  INTREPID_TEST_COMMAND(atools.matvecProductDataData<double>(fc_C_P_D3, fc_C_P_D3, fc_C1_P_D3) );
498  INTREPID_TEST_COMMAND(atools.matvecProductDataData<double>(fc_C_P_D3, fc_C_P_D3, fc_C_P1_D3) );
499  INTREPID_TEST_COMMAND(atools.matvecProductDataData<double>(fc_C_P_D3, fc_C_P_D3, fc_C_P_D2) );
500 
501  // Cross-check (1): inputDataLeft(C,P), (C,P,D), or (C,P,D,D) vs. inputDataRight(P,D): dimensions P, D must match
502  INTREPID_TEST_COMMAND(atools.matvecProductDataData<double>(fc_C_P_D3, fc_C_P_D3_D3, fc_P1_D3) );
503  INTREPID_TEST_COMMAND(atools.matvecProductDataData<double>(fc_C_P_D3, fc_C_P_D3_D3, fc_P_D2) );
504  INTREPID_TEST_COMMAND(atools.matvecProductDataData<double>(fc_C_P_D3, fc_C_P, fc_P1_D3) );
505  INTREPID_TEST_COMMAND(atools.matvecProductDataData<double>(fc_C_P_D3, fc_C_P_D3, fc_P1_D3) );
506  INTREPID_TEST_COMMAND(atools.matvecProductDataData<double>(fc_C_P_D3, fc_C_P_D3, fc_P_D2) );
507 
508  // Cross-check (3): outputData(C,P,D) vs. inputDataRight(C,P,D): all dimensions C, P, D must match
509  INTREPID_TEST_COMMAND(atools.matvecProductDataData<double>(fc_C1_P_D3, fc_C1_P_D3_D3, fc_C_P_D3) );
510  INTREPID_TEST_COMMAND(atools.matvecProductDataData<double>(fc_C_P1_D3, fc_C_P1_D3_D3, fc_C_P_D3) );
511  INTREPID_TEST_COMMAND(atools.matvecProductDataData<double>(fc_C_P_D2, fc_C_P_D3_D3, fc_C_P_D3) );
512 
513  // Cross-check (3): outputData(C,P,D) vs. inputDataRight(P,D): dimensions P, D must match
514  INTREPID_TEST_COMMAND(atools.matvecProductDataData<double>(fc_C_P1_D3, fc_C_P1_D3_D3, fc_P_D3) );
515  INTREPID_TEST_COMMAND(atools.matvecProductDataData<double>(fc_C_P_D3, fc_C_P_D3_D3, fc_P_D2) );
516 
517  *outStream \
518  << "\n"
519  << "===============================================================================\n"\
520  << "| TEST 7: matmatProductDataField exceptions |\n"\
521  << "===============================================================================\n";
522  // 46 exceptions
523  // (1) inputData is (C,P), (C,P,D), or (C,P,D,D) and 1 <= D <= 3 is required
524  INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C, fc_C_F_P_D3_D3) );
525  INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3_D3_D3, fc_C_F_P_D3_D3) );
526  INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D1_D3, fc_C_F_P_D3_D3) );
527  INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3_D1, fc_C_F_P_D3_D3) );
528 
529  // (2) inputFields is (C,F,P,D,D) or (F,P,D,D) and 1 <= (dimension(rank-1), (rank-2)) <= 3 is required.
530  INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3_D3, fc_C_F_P) );
531  INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3_D3, fc_C_F_P_D3_D3_D3) );
532  INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3_D3, fc_C_F_P_D1_D3) );
533  INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3_D3, fc_C_F_P_D3_D1) );
534 
535  // (3) outputFields is (C,F,P,D,D)
536  INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3, fc_C_P_D3_D3, fc_C_F_P_D3_D3) );
537  INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3_D3, fc_C_P_D3_D3, fc_C_F_P_D3_D3) );
538  INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D1_D3, fc_C_P_D3_D3, fc_C_F_P_D3_D3) );
539  INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D1, fc_C_P_D3_D3, fc_C_F_P_D3_D3) );
540 
541  // Cross-check (2): outputFields(C,F,P,D,D) vs. inputData(C,P), (C,P,D), or (C,P,D,D): dimensions C, P, D must match
542  INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C1_P_D3_D3, fc_C_F_P_D3_D3) );
543  INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P1_D3_D3, fc_C_F_P_D3_D3) );
544  INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D2_D2, fc_C_F_P_D3_D3) );
545  INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P1_D2_D2, fc_C_F_P_D3_D3) );
546  INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D2_D2, fc_C_P_D3_D3, fc_C_F_P_D3_D3) );
547  INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C1_P, fc_C_F_P_D3_D3) );
548  INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P1, fc_C_F_P_D3_D3) );
549  INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C1_P_D3, fc_C_F_P_D3_D3) );
550  INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P1_D3, fc_C_F_P_D3_D3) );
551  INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D1, fc_C_F_P_D3_D3) );
552 
553  // Cross-check (1): inputData(C,P), (C,P,D) or (C,P,D,D) vs. inputFields(C,F,P,D,D): dimensions C, P, D must match
554  INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3_D3, fc_C1_F_P_D3_D3) );
555  INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3_D3, fc_C_F_P1_D3_D3) );
556  INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3_D3, fc_C_F_P_D2_D2) );
557  INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3_D3, fc_C_F_P1_D2_D2) );
558  INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3_D3, fc_C1_F_P_D2_D2) );
559  INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3_D3, fc_C1_F_P1_D2_D2) );
560  INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P, fc_C1_F_P_D3_D3) );
561  INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P, fc_C_F_P1_D3_D3) );
562  INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3, fc_C1_F_P_D3_D3) );
563  INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3, fc_C_F_P1_D3_D3) );
564  INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3, fc_C_F_P_D2_D2) );
565 
566  // Cross-check (1): inputData(C,P), (C,P,D), or (C,P,D,D) vs. inputFields(F,P,D,D): dimensions P, D must match
567  INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3_D3, fc_F_P1_D3_D3) );
568  INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3_D3, fc_F_P_D2_D2) );
569  INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3_D3, fc_F_P1_D2_D2) );
570  INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P, fc_F_P1_D3_D3) );
571  INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3, fc_F_P1_D3_D3) );
572  INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3, fc_F_P_D2_D2) );
573 
574  // Cross-check (3): outputFields(C,F,P,D,D) vs. inputFields(C,F,P,D,D): all dimensions C, F, P, D must match
575  INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3_D3, fc_C1_F_P_D3_D3) );
576  INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3_D3, fc_C_F1_P_D3_D3) );
577  INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3_D3, fc_C_F_P1_D3_D3) );
578  INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3_D3, fc_C_F_P_D2_D2) );
579 
580  // Cross-check (3): outputFields(C,F,P,D,D) vs. inputFields(F,P,D,D): dimensions F, P, D must match
581  INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3_D3, fc_F1_P_D3_D3) );
582  INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3_D3, fc_F_P1_D3_D3) );
583  INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3_D3, fc_F_P_D2_D2) );
584  *outStream \
585  << "\n"
586  << "===============================================================================\n"\
587  << "| TEST 8: matmatProductDataData exceptions |\n"\
588  << "===============================================================================\n";
589  // 45 exceptions
590  // (1) inputDataLeft is (C,P), (C,P,D), or (C,P,D,D) and 2 <= D <= 3 is required
591  INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3, fc_C, fc_C_P_D3_D3) );
592  INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D3_D3_D3, fc_C_P_D3_D3) );
593  INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D1_D3, fc_C_P_D3_D3) );
594  INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D3_D1, fc_C_P_D3_D3) );
595 
596  // (2) inputDataRight is (C,P,D,D) or (P,D,D) and 1 <= (dimension(rank-1), (rank-2)) <= 3 is required.
597  INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3, fc_C_P, fc_C_P) );
598  INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D3, fc_C_P_D3_D3_D3) );
599  INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D3_D3, fc_C_P_D1_D3) );
600  INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D3_D3, fc_C_P_D3_D1) );
601 
602  // (3) outputData is (C,P,D,D)
603  INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3, fc_C_P, fc_C_P_D3_D3) );
604  INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3_D3, fc_C_P_D3, fc_C_P_D3_D3) );
605  INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D1_D3, fc_C_P_D3_D3, fc_C_P_D3_D3) );
606  INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D1, fc_C_P_D3_D3, fc_C_P_D3_D3) );
607 
608  // Cross-check (2): outputData(C,P,D,D) vs. inputDataLeft(C,P), (C,P,D), or (C,P,D,D): dimensions C, P, D must match
609  INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3, fc_C1_P_D3_D3, fc_C_P_D3_D3) );
610  INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3, fc_C_P1_D3_D3, fc_C_P_D3_D3) );
611  INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D2_D2, fc_C_P_D3_D3) );
612  INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3, fc_C_P1_D2_D2, fc_C_P_D3_D3) );
613  INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D2_D2, fc_C_P_D3_D3, fc_C_P_D3_D3) );
614  INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3, fc_C1_P, fc_C_P_D3_D3) );
615  INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3, fc_C_P1, fc_C_P_D3_D3) );
616  INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3, fc_C1_P_D3, fc_C_P_D3_D3) );
617  INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3, fc_C_P1_D3, fc_C_P_D3_D3) );
618  INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D1, fc_C_P_D3_D3) );
619 
620  // Cross-check (1): inputDataLeft(C,P), (C,P,D) or (C,P,D,D) vs. inputDataRight(C,P,D,D): dimensions C, P, D must match
621  INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D3_D3, fc_C1_P_D3_D3) );
622  INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D3_D3, fc_C_P1_D3_D3) );
623  INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D3_D3, fc_C_P_D2_D2) );
624  INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D3_D3, fc_C_P1_D2_D2) );
625  INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D3_D3, fc_C1_P_D2_D2) );
626  INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D3_D3, fc_C_P1_D2_D2) );
627  INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3, fc_C_P, fc_C1_P_D3_D3) );
628  INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3, fc_C_P, fc_C_P1_D3_D3) );
629  INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D3, fc_C1_P_D3_D3) );
630  INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D3, fc_C_P1_D3_D3) );
631  INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D3, fc_C_P_D2_D2) );
632 
633  // Cross-check (1): inputDataLeft(C,P), (C,P,D), or (C,P,D,D) vs. inputDataRight(P,D,D): dimensions P, D must match
634  INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D3_D3, fc_P1_D3_D3) );
635  INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D3_D3, fc_P_D2_D2) );
636  INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D3_D3, fc_P1_D2_D2) );
637  INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_P_D3_D3, fc_C_P, fc_P1_D3_D3) );
638  INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_P_D3_D3, fc_C_P_D3, fc_P1_D3_D3) );
639  INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_P_D3_D3, fc_C_P_D3, fc_P_D2_D2) );
640 
641  // Cross-check (3): outputData(C,P,D,D) vs. inputDataRight(C,P,D,D): all dimensions C, P, D must match
642  INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D3_D3, fc_C1_P_D3_D3) );
643  INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D3_D3, fc_C1_P_D3_D3) );
644  INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D3_D3, fc_C_P1_D3_D3) );
645  INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D3_D3, fc_C_P_D2_D2) );
646 
647  // Cross-check (3): outputData(C,P,D,D) vs. inputDataRight(P,D,D): dimensions P, D must match
648  INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D3_D3, fc_P_D2_D2) );
649  INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D3_D3, fc_P1_D3_D3) );
650  }
651 
652  catch (const std::logic_error & err) {
653  *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
654  *outStream << err.what() << '\n';
655  *outStream << "-------------------------------------------------------------------------------" << "\n\n";
656  errorFlag = -1000;
657  };
658 
659  if (Teuchos::TestForException_getThrowNumber() != endThrowNumber)
660  errorFlag++;
661 #endif
662 
663  /************************************************************************************************
664  * *
665  * Operation corectness tests *
666  * *
667  ***********************************************************************************************/
668 
669  try{
670  *outStream \
671  << "\n"
672  << "===============================================================================\n"\
673  << "| TEST 1.a: 3D crossProductDataField operations: (C,P,D) and (C,F,P,D) |\n"\
674  << "===============================================================================\n";
675  /*
676  * ijkData_1a ijkFields_1a Expected result in (C,F,P,D) array
677  * i,i i,i j,j, k,k 0, 0 k, k -j,-j
678  * j,j i,i j,j, k,k -k,-k 0, 0 i, i
679  * k,k i,i j,j, k,k j, j -i,-i 0, 0
680  */
681 
682  // input data is (C,P,D)
683  FieldContainer<double> ijkData_1a(3, 2, 3);
684  // C=0 contains i
685  ijkData_1a(0, 0, 0) = 1.0; ijkData_1a(0, 0, 1) = 0.0; ijkData_1a(0, 0, 2) = 0.0;
686  ijkData_1a(0, 1, 0) = 1.0; ijkData_1a(0, 1, 1) = 0.0; ijkData_1a(0, 1, 2) = 0.0;
687  // C=1 contains j
688  ijkData_1a(1, 0, 0) = 0.0; ijkData_1a(1, 0, 1) = 1.0; ijkData_1a(1, 0, 2) = 0.0;
689  ijkData_1a(1, 1, 0) = 0.0; ijkData_1a(1, 1, 1) = 1.0; ijkData_1a(1, 1, 2) = 0.0;
690  // C=2 contains k
691  ijkData_1a(2, 0, 0) = 0.0; ijkData_1a(2, 0, 1) = 0.0; ijkData_1a(2, 0, 2) = 1.0;
692  ijkData_1a(2, 1, 0) = 0.0; ijkData_1a(2, 1, 1) = 0.0; ijkData_1a(2, 1, 2) = 1.0;
693 
694 
695  FieldContainer<double> ijkFields_1a(3, 3, 2, 3);
696  // C=0, F=0 is i
697  ijkFields_1a(0, 0, 0, 0) = 1.0; ijkFields_1a(0, 0, 0, 1) = 0.0; ijkFields_1a(0, 0, 0, 2) = 0.0;
698  ijkFields_1a(0, 0, 1, 0) = 1.0; ijkFields_1a(0, 0, 1, 1) = 0.0; ijkFields_1a(0, 0, 1, 2) = 0.0;
699  // C=0, F=1 is j
700  ijkFields_1a(0, 1, 0, 0) = 0.0; ijkFields_1a(0, 1, 0, 1) = 1.0; ijkFields_1a(0, 1, 0, 2) = 0.0;
701  ijkFields_1a(0, 1, 1, 0) = 0.0; ijkFields_1a(0, 1, 1, 1) = 1.0; ijkFields_1a(0, 1, 1, 2) = 0.0;
702  // C=0, F=2 is k
703  ijkFields_1a(0, 2, 0, 0) = 0.0; ijkFields_1a(0, 2, 0, 1) = 0.0; ijkFields_1a(0, 2, 0, 2) = 1.0;
704  ijkFields_1a(0, 2, 1, 0) = 0.0; ijkFields_1a(0, 2, 1, 1) = 0.0; ijkFields_1a(0, 2, 1, 2) = 1.0;
705 
706  // C=1, F=0 is i
707  ijkFields_1a(1, 0, 0, 0) = 1.0; ijkFields_1a(1, 0, 0, 1) = 0.0; ijkFields_1a(1, 0, 0, 2) = 0.0;
708  ijkFields_1a(1, 0, 1, 0) = 1.0; ijkFields_1a(1, 0, 1, 1) = 0.0; ijkFields_1a(1, 0, 1, 2) = 0.0;
709  // C=1, F=1 is j
710  ijkFields_1a(1, 1, 0, 0) = 0.0; ijkFields_1a(1, 1, 0, 1) = 1.0; ijkFields_1a(1, 1, 0, 2) = 0.0;
711  ijkFields_1a(1, 1, 1, 0) = 0.0; ijkFields_1a(1, 1, 1, 1) = 1.0; ijkFields_1a(1, 1, 1, 2) = 0.0;
712  // C=1, F=2 is k
713  ijkFields_1a(1, 2, 0, 0) = 0.0; ijkFields_1a(1, 2, 0, 1) = 0.0; ijkFields_1a(1, 2, 0, 2) = 1.0;
714  ijkFields_1a(1, 2, 1, 0) = 0.0; ijkFields_1a(1, 2, 1, 1) = 0.0; ijkFields_1a(1, 2, 1, 2) = 1.0;
715 
716  // C=2, F=0 is i
717  ijkFields_1a(2, 0, 0, 0) = 1.0; ijkFields_1a(2, 0, 0, 1) = 0.0; ijkFields_1a(2, 0, 0, 2) = 0.0;
718  ijkFields_1a(2, 0, 1, 0) = 1.0; ijkFields_1a(2, 0, 1, 1) = 0.0; ijkFields_1a(2, 0, 1, 2) = 0.0;
719  // C=2, F=1 is j
720  ijkFields_1a(2, 1, 0, 0) = 0.0; ijkFields_1a(2, 1, 0, 1) = 1.0; ijkFields_1a(2, 1, 0, 2) = 0.0;
721  ijkFields_1a(2, 1, 1, 0) = 0.0; ijkFields_1a(2, 1, 1, 1) = 1.0; ijkFields_1a(2, 1, 1, 2) = 0.0;
722  // C=2, F=2 is k
723  ijkFields_1a(2, 2, 0, 0) = 0.0; ijkFields_1a(2, 2, 0, 1) = 0.0; ijkFields_1a(2, 2, 0, 2) = 1.0;
724  ijkFields_1a(2, 2, 1, 0) = 0.0; ijkFields_1a(2, 2, 1, 1) = 0.0; ijkFields_1a(2, 2, 1, 2) = 1.0;
725 
726 
727  FieldContainer<double> outFields(3, 3, 2, 3);
728  art::crossProductDataField<double>(outFields, ijkData_1a, ijkFields_1a);
729 
730  // checks for C = 0
731  if( !(outFields(0,0,0,0)==0.0 && outFields(0,0,0,1)==0.0 && outFields(0,0,0,2)==0.0 &&
732  outFields(0,0,1,0)==0.0 && outFields(0,0,1,1)==0.0 && outFields(0,0,1,2)==0.0 ) ) {
733  *outStream << "\n\nINCORRECT crossProductDataField (1): i x i != 0; ";
734  errorFlag++;
735  }
736  if( !(outFields(0,1,0,0)==0.0 && outFields(0,1,0,1)==0.0 && outFields(0,1,0,2)==1.0 &&
737  outFields(0,1,1,0)==0.0 && outFields(0,1,1,1)==0.0 && outFields(0,1,1,2)==1.0 ) ) {
738  *outStream << "\n\nINCORRECT crossProductDataField (2): i x j != k; ";
739  errorFlag++;
740  }
741  if( !(outFields(0,2,0,0)==0.0 && outFields(0,2,0,1)==-1.0 && outFields(0,2,0,2)==0.0 &&
742  outFields(0,2,1,0)==0.0 && outFields(0,2,1,1)==-1.0 && outFields(0,2,1,2)==0.0 ) ) {
743  *outStream << "\n\nINCORRECT crossProductDataField (3): i x k != -j; ";
744  errorFlag++;
745  }
746 
747  // checks for C = 1
748  if( !(outFields(1,0,0,0)==0.0 && outFields(1,0,0,1)==0.0 && outFields(1,0,0,2)==-1.0 &&
749  outFields(1,0,1,0)==0.0 && outFields(1,0,1,1)==0.0 && outFields(1,0,1,2)==-1.0 ) ) {
750  *outStream << "\n\nINCORRECT crossProductDataField (4): j x i != -k; ";
751  errorFlag++;
752  }
753  if( !(outFields(1,1,0,0)==0.0 && outFields(1,1,0,1)==0.0 && outFields(1,1,0,2)==0.0 &&
754  outFields(1,1,1,0)==0.0 && outFields(1,1,1,1)==0.0 && outFields(1,1,1,2)==0.0 ) ) {
755  *outStream << "\n\nINCORRECT crossProductDataField (5): j x j != 0; ";
756  errorFlag++;
757  }
758  if( !(outFields(1,2,0,0)==1.0 && outFields(1,2,0,1)==0.0 && outFields(1,2,0,2)==0.0 &&
759  outFields(1,2,1,0)==1.0 && outFields(1,2,1,1)==0.0 && outFields(1,2,1,2)==0.0 ) ) {
760  *outStream << "\n\nINCORRECT crossProductDataField (6): j x k != i; ";
761  errorFlag++;
762  }
763 
764  // checks for C = 2
765  if( !(outFields(2,0,0,0)==0.0 && outFields(2,0,0,1)==1.0 && outFields(2,0,0,2)==0.0 &&
766  outFields(2,0,1,0)==0.0 && outFields(2,0,1,1)==1.0 && outFields(2,0,1,2)==0.0 ) ) {
767  *outStream << "\n\nINCORRECT crossProductDataField (7): k x i != j; ";
768  errorFlag++;
769  }
770  if( !(outFields(2,1,0,0)==-1.0 && outFields(2,1,0,1)==0.0 && outFields(2,1,0,2)==0.0 &&
771  outFields(2,1,1,0)==-1.0 && outFields(2,1,1,1)==0.0 && outFields(2,1,1,2)==0.0 ) ) {
772  *outStream << "\n\nINCORRECT crossProductDataField (8): k x j != -i; ";
773  errorFlag++;
774  }
775  if( !(outFields(2,2,0,0)==0.0 && outFields(2,2,0,1)==0.0 && outFields(2,2,0,2)==0.0 &&
776  outFields(2,2,1,0)==0.0 && outFields(2,2,1,1)==0.0 && outFields(2,2,1,2)==0.0 ) ) {
777  *outStream << "\n\nINCORRECT crossProductDataField (9): k x k != 0; ";
778  errorFlag++;
779  }
780 
781  *outStream \
782  << "\n"
783  << "===============================================================================\n"\
784  << "| TEST 1.b: 3D crossProductDataField operations: (C,P,D) and (F,P,D) |\n"\
785  << "===============================================================================\n";
786  /*
787  * ijkData_1b ijkFields_1b expected result in (C,F,P,D) array
788  * i,i,i i,j,k 0, k,-j
789  * j,j,j -k, 0, i
790  * k,k,k j,-i, 0
791  */
792 
793  // input data is (C,P,D)
794  FieldContainer<double> ijkData_1b(3, 3, 3);
795  // C=0 contains i
796  ijkData_1b(0, 0, 0) = 1.0; ijkData_1b(0, 0, 1) = 0.0; ijkData_1b(0, 0, 2) = 0.0;
797  ijkData_1b(0, 1, 0) = 1.0; ijkData_1b(0, 1, 1) = 0.0; ijkData_1b(0, 1, 2) = 0.0;
798  ijkData_1b(0, 2, 0) = 1.0; ijkData_1b(0, 2, 1) = 0.0; ijkData_1b(0, 2, 2) = 0.0;
799  // C=1 contains j
800  ijkData_1b(1, 0, 0) = 0.0; ijkData_1b(1, 0, 1) = 1.0; ijkData_1b(1, 0, 2) = 0.0;
801  ijkData_1b(1, 1, 0) = 0.0; ijkData_1b(1, 1, 1) = 1.0; ijkData_1b(1, 1, 2) = 0.0;
802  ijkData_1b(1, 2, 0) = 0.0; ijkData_1b(1, 2, 1) = 1.0; ijkData_1b(1, 2, 2) = 0.0;
803  // C=2 contains k
804  ijkData_1b(2, 0, 0) = 0.0; ijkData_1b(2, 0, 1) = 0.0; ijkData_1b(2, 0, 2) = 1.0;
805  ijkData_1b(2, 1, 0) = 0.0; ijkData_1b(2, 1, 1) = 0.0; ijkData_1b(2, 1, 2) = 1.0;
806  ijkData_1b(2, 2, 0) = 0.0; ijkData_1b(2, 2, 1) = 0.0; ijkData_1b(2, 2, 2) = 1.0;
807 
808  // input fields are (F,P,D)
809  FieldContainer<double> ijkFields_1b(1, 3, 3);
810  // F=0 at 3 points is (i,j,k)
811  ijkFields_1b(0, 0, 0) = 1.0; ijkFields_1b(0, 0, 1) = 0.0; ijkFields_1b(0, 0, 2) = 0.0;
812  ijkFields_1b(0, 1, 0) = 0.0; ijkFields_1b(0, 1, 1) = 1.0; ijkFields_1b(0, 1, 2) = 0.0;
813  ijkFields_1b(0, 2, 0) = 0.0; ijkFields_1b(0, 2, 1) = 0.0; ijkFields_1b(0, 2, 2) = 1.0;
814 
815  // Output array is (C,F,P,D)
816  outFields.resize(3, 1, 3, 3);
817  art::crossProductDataField<double>(outFields, ijkData_1b, ijkFields_1b);
818 
819  // checks for C = 0
820  if( !(outFields(0,0,0,0)==0.0 && outFields(0,0,0,1)==0.0 && outFields(0,0,0,2)==0.0) ) {
821  *outStream << "\n\nINCORRECT crossProductDataField (10): i x i != 0; ";
822  errorFlag++;
823  }
824  if( !(outFields(0,0,1,0)==0.0 && outFields(0,0,1,1)==0.0 && outFields(0,0,1,2)==1.0) ) {
825  *outStream << "\n\nINCORRECT crossProductDataField (11): i x j != k; ";
826  errorFlag++;
827  }
828  if( !(outFields(0,0,2,0)==0.0 && outFields(0,0,2,1)==-1.0 && outFields(0,0,2,2)==0.0) ) {
829  *outStream << "\n\nINCORRECT crossProductDataField (12): i x k != -j; ";
830  errorFlag++;
831  }
832 
833  // checks for C = 1
834  if( !(outFields(1,0,0,0)==0.0 && outFields(1,0,0,1)==0.0 && outFields(1,0,0,2)==-1.0) ) {
835  *outStream << "\n\nINCORRECT crossProductDataField (13): j x i != -k; ";
836  errorFlag++;
837  }
838  if( !(outFields(1,0,1,0)==0.0 && outFields(1,0,1,1)==0.0 && outFields(1,0,1,2)==0.0) ) {
839  *outStream << "\n\nINCORRECT crossProductDataField (14): j x j != 0; ";
840  errorFlag++;
841  }
842  if( !(outFields(1,0,2,0)==1.0 && outFields(1,0,2,1)==0.0 && outFields(1,0,2,2)==0.0) ) {
843  *outStream << "\n\nINCORRECT crossProductDataField (15): j x k != i; ";
844  errorFlag++;
845  }
846 
847  // checks for C = 2
848  if( !(outFields(2,0,0,0)==0.0 && outFields(2,0,0,1)==1.0 && outFields(2,0,0,2)==0.0) ) {
849  *outStream << "\n\nINCORRECT crossProductDataField (16): k x i != j; ";
850  errorFlag++;
851  }
852  if( !(outFields(2,0,1,0)==-1.0 && outFields(2,0,1,1)==0.0 && outFields(2,0,1,2)==0.0) ) {
853  *outStream << "\n\nINCORRECT crossProductDataField (17): k x j != -i; ";
854  errorFlag++;
855  }
856  if( !(outFields(2,0,2,0)==0.0 && outFields(2,0,2,1)==0.0 && outFields(2,0,2,2)==0.0) ) {
857  *outStream << "\n\nINCORRECT crossProductDataField (18): k x k != 0; ";
858  errorFlag++;
859  }
860 
861  *outStream \
862  << "\n"
863  << "===============================================================================\n"\
864  << "| TEST 1.c: 2D crossProductDataField operations: (C,P,D) and (C,F,P,D) |\n"\
865  << "===============================================================================\n";
866 
867  /*
868  * ijData_1c ijFields_1c expected result in (C,F,P) array
869  * i,i i,i j,j 0, 0 1, 1
870  * j,j i,i j,j -1,-1 0, 0
871  */
872  // input data is (C,P,D)
873  FieldContainer<double> ijData_1c(2, 2, 2);
874  // C=0 contains i
875  ijData_1c(0, 0, 0) = 1.0; ijData_1c(0, 0, 1) = 0.0;
876  ijData_1c(0, 1, 0) = 1.0; ijData_1c(0, 1, 1) = 0.0;
877  // C=1 contains j
878  ijData_1c(1, 0, 0) = 0.0; ijData_1c(1, 0, 1) = 1.0;
879  ijData_1c(1, 1, 0) = 0.0; ijData_1c(1, 1, 1) = 1.0;
880 
881 
882  FieldContainer<double> ijFields_1c(2, 2, 2, 2);
883  // C=0, F=0 is i
884  ijFields_1c(0, 0, 0, 0) = 1.0; ijFields_1c(0, 0, 0, 1) = 0.0;
885  ijFields_1c(0, 0, 1, 0) = 1.0; ijFields_1c(0, 0, 1, 1) = 0.0;
886  // C=0, F=1 is j
887  ijFields_1c(0, 1, 0, 0) = 0.0; ijFields_1c(0, 1, 0, 1) = 1.0;
888  ijFields_1c(0, 1, 1, 0) = 0.0; ijFields_1c(0, 1, 1, 1) = 1.0;
889 
890  // C=1, F=0 is i
891  ijFields_1c(1, 0, 0, 0) = 1.0; ijFields_1c(1, 0, 0, 1) = 0.0;
892  ijFields_1c(1, 0, 1, 0) = 1.0; ijFields_1c(1, 0, 1, 1) = 0.0;
893  // C=1, F=1 is j
894  ijFields_1c(1, 1, 0, 0) = 0.0; ijFields_1c(1, 1, 0, 1) = 1.0;
895  ijFields_1c(1, 1, 1, 0) = 0.0; ijFields_1c(1, 1, 1, 1) = 1.0;
896 
897  // Output array is (C,F,P)
898  outFields.resize(2, 2, 2);
899  art::crossProductDataField<double>(outFields, ijData_1c, ijFields_1c);
900 
901  if( !(outFields(0,0,0)==0.0 && outFields(0,0,1)==0.0 ) ) {
902  *outStream << "\n\nINCORRECT crossProductDataField (19): i x i != 0; ";
903  errorFlag++;
904  }
905  if( !(outFields(0,1,0)==1.0 && outFields(0,1,1)==1.0 ) ) {
906  *outStream << "\n\nINCORRECT crossProductDataField (20): i x j != 1; ";
907  errorFlag++;
908  }
909 
910  if( !(outFields(1,0,0)==-1.0 && outFields(1,0,1)==-1.0 ) ) {
911  *outStream << "\n\nINCORRECT crossProductDataField (21): j x i != -1; ";
912  errorFlag++;
913  }
914  if( !(outFields(1,1,0)==0.0 && outFields(1,1,1)==0.0 ) ) {
915  *outStream << "\n\nINCORRECT crossProductDataField (22): j x j != 0; ";
916  errorFlag++;
917  }
918 
919  *outStream \
920  << "\n"
921  << "===============================================================================\n"\
922  << "| TEST 1.d: 2D crossProductDataField operations: (C,P,D) and (F,P,D) |\n"\
923  << "===============================================================================\n";
924  /*
925  * ijData_1c ijFields_1d expected result in (C,F,P) array
926  * i,i i,j 0, 1
927  * j,j -1, 0
928  */
929  // inputFields is (F,P,D)
930  FieldContainer<double> ijFields_1d(1, 2, 2);
931  // F=0 at 2 points is i,j
932  ijFields_1d(0, 0, 0) = 1.0; ijFields_1d(0, 0, 1) = 0.0;
933  ijFields_1d(0, 1, 0) = 0.0; ijFields_1d(0, 1, 1) = 1.0;
934 
935  // Output array is (C,F,P)
936  outFields.resize(2, 1, 2);
937  art::crossProductDataField<double>(outFields, ijData_1c, ijFields_1d);
938 
939  if( !(outFields(0,0,0)==0.0 ) ) {
940  *outStream << "\n\nINCORRECT crossProductDataField (23): i x i != 0; ";
941  errorFlag++;
942  }
943  if( !(outFields(0,0,1)==1.0 ) ) {
944  *outStream << "\n\nINCORRECT crossProductDataField (24): i x j != 1; ";
945  errorFlag++;
946  }
947  if( !(outFields(1,0,0)==-1.0 ) ) {
948  *outStream << "\n\nINCORRECT crossProductDataField (25): j x i != -1; ";
949  errorFlag++;
950  }
951  if( !(outFields(1,0,1)==0.0 ) ) {
952  *outStream << "\n\nINCORRECT crossProductDataField (26): j x j != 0; ";
953  errorFlag++;
954  }
955 
956 
957  *outStream \
958  << "\n"
959  << "===============================================================================\n"\
960  << "| TEST 2.a: 3D crossProductDataData operations: (C,P,D) and (C,P,D) |\n"\
961  << "===============================================================================\n";
962  /*
963  * ijkData_1a jkiData_2a kijData_2a
964  * i,i j,j k,k
965  * j,j k,k i,i
966  * k,k i,i j,j
967  */
968  FieldContainer<double> jkiData_2a(3, 2, 3);
969  // C=0 contains j
970  jkiData_2a(0, 0, 0) = 0.0; jkiData_2a(0, 0, 1) = 1.0; jkiData_2a(0, 0, 2) = 0.0;
971  jkiData_2a(0, 1, 0) = 0.0; jkiData_2a(0, 1, 1) = 1.0; jkiData_2a(0, 1, 2) = 0.0;
972  // C=1 contains k
973  jkiData_2a(1, 0, 0) = 0.0; jkiData_2a(1, 0, 1) = 0.0; jkiData_2a(1, 0, 2) = 1.0;
974  jkiData_2a(1, 1, 0) = 0.0; jkiData_2a(1, 1, 1) = 0.0; jkiData_2a(1, 1, 2) = 1.0;
975  // C=2 contains i
976  jkiData_2a(2, 0, 0) = 1.0; jkiData_2a(2, 0, 1) = 0.0; jkiData_2a(2, 0, 2) = 0.0;
977  jkiData_2a(2, 1, 0) = 1.0; jkiData_2a(2, 1, 1) = 0.0; jkiData_2a(2, 1, 2) = 0.0;
978 
979  FieldContainer<double> kijData_2a(3, 2, 3);
980  // C=0 contains k
981  kijData_2a(0, 0, 0) = 0.0; kijData_2a(0, 0, 1) = 0.0; kijData_2a(0, 0, 2) = 1.0;
982  kijData_2a(0, 1, 0) = 0.0; kijData_2a(0, 1, 1) = 0.0; kijData_2a(0, 1, 2) = 1.0;
983  // C=1 contains i
984  kijData_2a(1, 0, 0) = 1.0; kijData_2a(1, 0, 1) = 0.0; kijData_2a(1, 0, 2) = 0.0;
985  kijData_2a(1, 1, 0) = 1.0; kijData_2a(1, 1, 1) = 0.0; kijData_2a(1, 1, 2) = 0.0;
986  // C=2 contains j
987  kijData_2a(2, 0, 0) = 0.0; kijData_2a(2, 0, 1) = 1.0; kijData_2a(2, 0, 2) = 0.0;
988  kijData_2a(2, 1, 0) = 0.0; kijData_2a(2, 1, 1) = 1.0; kijData_2a(2, 1, 2) = 0.0;
989 
990 
991  // ijkData_1a x ijkData_1a: outData should contain ixi=0, jxj=0, kxk=0
992  FieldContainer<double> outData(3,2,3);
993  art::crossProductDataData<double>(outData, ijkData_1a, ijkData_1a);
994 
995  for(int i = 0; i < outData.size(); i++){
996  if(outData[i] != 0) {
997  *outStream << "\n\nINCORRECT crossProductDataData (1): i x i, j x j, or k x k != 0; ";
998  errorFlag++;
999  }
1000  }
1001 
1002 
1003  // ijkData_1a x jkiData_2a
1004  art::crossProductDataData<double>(outData, ijkData_1a, jkiData_2a);
1005 
1006  // cell 0 should contain i x j = k
1007  if( !( outData(0,0,0)==0.0 && outData(0,0,1)==0.0 && outData(0,0,2)==1.0 &&
1008  outData(0,1,0)==0.0 && outData(0,1,1)==0.0 && outData(0,1,2)==1.0) ) {
1009  *outStream << "\n\nINCORRECT crossProductDataData (2): i x j != k; ";
1010  errorFlag++;
1011  }
1012 
1013  // cell 1 should contain j x k = i
1014  if( !( outData(1,0,0)==1.0 && outData(1,0,1)==0.0 && outData(1,0,2)==0.0 &&
1015  outData(1,1,0)==1.0 && outData(1,1,1)==0.0 && outData(1,1,2)==0.0) ) {
1016  *outStream << "\n\nINCORRECT crossProductDataData (3): j x k != i; ";
1017  errorFlag++;
1018  }
1019 
1020  // cell 2 should contain k x i = j
1021  if( !( outData(2,0,0)==0.0 && outData(2,0,1)==1.0 && outData(2,0,2)==0.0 &&
1022  outData(2,1,0)==0.0 && outData(2,1,1)==1.0 && outData(2,1,2)==0.0) ) {
1023  *outStream << "\n\nINCORRECT crossProductDataData (4): k x i != j; ";
1024  errorFlag++;
1025  }
1026 
1027 
1028  // ijkData_1a x kijData_2a
1029  art::crossProductDataData<double>(outData, ijkData_1a, kijData_2a);
1030 
1031  // cell 0 should contain i x k = -j
1032  if( !( outData(0,0,0)==0.0 && outData(0,0,1)==-1.0 && outData(0,0,2)==0.0 &&
1033  outData(0,1,0)==0.0 && outData(0,1,1)==-1.0 && outData(0,1,2)==0.0) ) {
1034  *outStream << "\n\nINCORRECT crossProductDataData (5): i x k != -j; ";
1035  errorFlag++;
1036  }
1037 
1038  // cell 1 should contain j x i = -k
1039  if( !( outData(1,0,0)==0.0 && outData(1,0,1)==0.0 && outData(1,0,2)==-1.0 &&
1040  outData(1,1,0)==0.0 && outData(1,1,1)==0.0 && outData(1,1,2)==-1.0) ) {
1041  *outStream << "\n\nINCORRECT crossProductDataData (6): j x i != -k; ";
1042  errorFlag++;
1043  }
1044 
1045  // cell 2 should contain k x j = -i
1046  if( !( outData(2,0,0)==-1.0 && outData(2,0,1)==0.0 && outData(2,0,2)==0.0 &&
1047  outData(2,1,0)==-1.0 && outData(2,1,1)==0.0 && outData(2,1,2)==0.0) ) {
1048  *outStream << "\n\nINCORRECT crossProductDataData (7): k x j != -i; ";
1049  errorFlag++;
1050  }
1051 
1052 
1053  *outStream \
1054  << "\n"
1055  << "===============================================================================\n"\
1056  << "| TEST 2.b: 3D crossProductDataData operations: (C,P,D) and (P,D) |\n"\
1057  << "===============================================================================\n";
1058  /*
1059  * ijkData_1b ijkData_2b expected result in (C,P,D) array
1060  * i,i,i i,j,k 0, k,-j
1061  * j,j,j -k, 0, i
1062  * k,k,k j,-i, 0
1063  */
1064  // input data is (P,D)
1065  FieldContainer<double> ijkData_2b(3, 3);
1066  // F=0 at 3 points is (i,j,k)
1067  ijkData_2b(0, 0) = 1.0; ijkData_2b(0, 1) = 0.0; ijkData_2b(0, 2) = 0.0;
1068  ijkData_2b(1, 0) = 0.0; ijkData_2b(1, 1) = 1.0; ijkData_2b(1, 2) = 0.0;
1069  ijkData_2b(2, 0) = 0.0; ijkData_2b(2, 1) = 0.0; ijkData_2b(2, 2) = 1.0;
1070 
1071  // Output array is (C,P,D)
1072  outData.resize(3, 3, 3);
1073  art::crossProductDataData<double>(outData, ijkData_1b, ijkData_2b);
1074 
1075  // checks for C = 0
1076  if( !(outData(0,0,0)==0.0 && outData(0,0,1)==0.0 && outData(0,0,2)==0.0) ) {
1077  *outStream << "\n\nINCORRECT crossProductDataData (5): i x i != 0; ";
1078  errorFlag++;
1079  }
1080  if( !(outData(0,1,0)==0.0 && outData(0,1,1)==0.0 && outData(0,1,2)==1.0) ) {
1081  *outStream << "\n\nINCORRECT crossProductDataData (6): i x j != k; ";
1082  errorFlag++;
1083  }
1084  if( !(outData(0,2,0)==0.0 && outData(0,2,1)==-1.0 && outData(0,2,2)==0.0) ) {
1085  *outStream << "\n\nINCORRECT crossProductDataData (7): i x k != -j; ";
1086  errorFlag++;
1087  }
1088 
1089  // checks for C = 1
1090  if( !(outData(1,0,0)==0.0 && outData(1,0,1)==0.0 && outData(1,0,2)==-1.0) ) {
1091  *outStream << "\n\nINCORRECT crossProductDataData (8): j x i != -k; ";
1092  errorFlag++;
1093  }
1094  if( !(outData(1,1,0)==0.0 && outData(1,1,1)==0.0 && outData(1,1,2)==0.0) ) {
1095  *outStream << "\n\nINCORRECT crossProductDataData (9): j x j != 0; ";
1096  errorFlag++;
1097  }
1098  if( !(outData(1,2,0)==1.0 && outData(1,2,1)==0.0 && outData(1,2,2)==0.0) ) {
1099  *outStream << "\n\nINCORRECT crossProductDataData (10): j x k != i; ";
1100  errorFlag++;
1101  }
1102 
1103  // checks for C = 2
1104  if( !(outData(2,0,0)==0.0 && outData(2,0,1)==1.0 && outData(2,0,2)==0.0) ) {
1105  *outStream << "\n\nINCORRECT crossProductDataData (11): k x i != j; ";
1106  errorFlag++;
1107  }
1108  if( !(outData(2,1,0)==-1.0 && outData(2,1,1)==0.0 && outData(2,1,2)==0.0) ) {
1109  *outStream << "\n\nINCORRECT crossProductDataData (12): k x j != -i; ";
1110  errorFlag++;
1111  }
1112  if( !(outData(2,2,0)==0.0 && outData(2,2,1)==0.0 && outData(2,2,2)==0.0) ) {
1113  *outStream << "\n\nINCORRECT crossProductDataData (13): k x k != 0; ";
1114  errorFlag++;
1115  }
1116 
1117 
1118  *outStream \
1119  << "\n"
1120  << "===============================================================================\n"\
1121  << "| TEST 2.c: 2D crossProductDataData operations: (C,P,D) and (C,P,D) |\n"\
1122  << "===============================================================================\n";
1123  /*
1124  * ijData_1c jiData_2c
1125  * i,i j,j
1126  * j,j i,i
1127  */
1128  FieldContainer<double> jiData_2c(2, 2, 2);
1129  // C=0 contains j
1130  jiData_2c(0, 0, 0) = 0.0; jiData_2c(0, 0, 1) = 1.0;
1131  jiData_2c(0, 1, 0) = 0.0; jiData_2c(0, 1, 1) = 1.0;
1132  // C=1 contains i
1133  jiData_2c(1, 0, 0) = 1.0; jiData_2c(1, 0, 1) = 0.0;
1134  jiData_2c(1, 1, 0) = 1.0; jiData_2c(1, 1, 1) = 0.0;
1135 
1136 
1137  // ijData_1c x ijData_1c: outData should contain ixi=0, jxj=0
1138  outData.resize(2,2);
1139  art::crossProductDataData<double>(outData, ijData_1c, ijData_1c);
1140 
1141  for(int i = 0; i < outData.size(); i++){
1142  if(outData[i] != 0) {
1143  *outStream << "\n\nINCORRECT crossProductDataData (14): i x i or j x j != 0; ";
1144  errorFlag++;
1145  }
1146  }
1147 
1148  // ijData_1c x jiData_1c: outData should contain ixi=0, jxj=0
1149  art::crossProductDataData<double>(outData, ijData_1c, jiData_2c);
1150 
1151  if( !(outData(0,0)==1.0 && outData(0,1)==1.0 ) ) {
1152  *outStream << "\n\nINCORRECT crossProductDataData (15): i x j != 1; ";
1153  errorFlag++;
1154  }
1155  if( !(outData(1,0)==-1.0 && outData(1,1)==-1.0 ) ) {
1156  *outStream << "\n\nINCORRECT crossProductDataData (16): j x i != -1; ";
1157  errorFlag++;
1158  }
1159 
1160  *outStream \
1161  << "\n"
1162  << "===============================================================================\n"\
1163  << "| TEST 2.d: 2D crossProductDataData operations: (C,P,D) and (P,D) |\n"\
1164  << "===============================================================================\n";
1165  /*
1166  * ijData_1c ijData_2d expected result in (C,P) array
1167  * i,i i,j 0, 1
1168  * j,j -1, 0
1169  */
1170  FieldContainer<double> ijData_2d(2, 2);
1171  ijData_2d(0, 0) = 1.0; ijData_2d(0, 1) = 0.0;
1172  ijData_2d(1, 0) = 0.0; ijData_2d(1, 1) = 1.0;
1173 
1174  outData.resize(2,2);
1175  art::crossProductDataData<double>(outData, ijData_1c, ijData_2d);
1176 
1177  if( !(outData(0,0)==0.0 ) ) {
1178  *outStream << "\n\nINCORRECT crossProductDataData (17): i x i != 0; ";
1179  errorFlag++;
1180  }
1181  if( !(outData(0,1)==1.0 ) ) {
1182  *outStream << "\n\nINCORRECT crossProductDataData (18): i x j != 1; ";
1183  errorFlag++;
1184  }
1185  if( !(outData(1,0)==-1.0 ) ) {
1186  *outStream << "\n\nINCORRECT crossProductDataData (19): j x i != -1; ";
1187  errorFlag++;
1188  }
1189  if( !(outData(1,1)==0.0 ) ) {
1190  *outStream << "\n\nINCORRECT crossProductDataData (20): j x j != 0; ";
1191  errorFlag++;
1192  }
1193 
1194 
1195  *outStream \
1196  << "\n"
1197  << "===============================================================================\n"\
1198  << "| TEST 3.a: outerProductDataField operations: (C,P,D) and (C,F,P,D) |\n"\
1199  << "===============================================================================\n";
1200  /*
1201  * ijkData_1a ijkFields_1a Expected result in (C,F,P,D,D) array:
1202  * i,i i,i j,j, k,k (0,0) (0,1) (0,2)
1203  * j,j i,i j,j, k,k (1,0) (1,1) (1,2)
1204  * k,k i,i j,j, k,k (2,0) (2,1) (2,2)
1205  * Indicates the only non-zero element of (C,F,P,*,*), specifically,
1206  * element with row = cell and col = field should equal 1; all other should equal 0
1207  */
1208 
1209  outFields.resize(3, 3, 2, 3, 3);
1210  art::outerProductDataField<double>(outFields, ijkData_1a, ijkFields_1a);
1211 
1212  for(int cell = 0; cell < ijkData_1a.dimension(0); cell++){
1213  for(int field = 0; field < ijkFields_1a.dimension(1); field++){
1214  for(int point = 0; point < ijkData_1a.dimension(1); point++){
1215  for(int row = 0; row < ijkData_1a.dimension(2); row++){
1216  for(int col = 0; col < ijkData_1a.dimension(2); col++){
1217 
1218  // element with row = cell and col = field should equal 1; all other should equal 0
1219  if( (row == cell && col == field) ){
1220  if(outFields(cell, field, point, row, col) != 1.0) {
1221  *outStream << "\n\nINCORRECT outerProductDataField (1): computed value is "
1222  << outFields(cell, field, point, row, col) << " whereas correct value is 1.0";
1223  errorFlag++;
1224  }
1225  }
1226  else {
1227  if(outFields(cell, field, point, row, col) != 0.0) {
1228  *outStream << "\n\nINCORRECT outerProductDataField (2): computed value is "
1229  << outFields(cell, field, point, row, col) << " whereas correct value is 0.0";
1230  errorFlag++;
1231  }
1232  } // if
1233  }// col
1234  }// row
1235  }// point
1236  }// field
1237  }// cell
1238 
1239  *outStream \
1240  << "\n"
1241  << "===============================================================================\n"\
1242  << "| TEST 3.b: outerProductDataField operations: (C,P,D) and (F,P,D) |\n"\
1243  << "===============================================================================\n";
1244  /*
1245  * ijkData_1b ijkFields_1b expected result in (C,F,P,D,D) array
1246  * i,i,i i,j,k (0,0) (0,1) (0,2)
1247  * j,j,j (1,0) (1,1) (1,2)
1248  * k,k,k (2,0) (2,1) (2,2)
1249  * Indicates the only non-zero element of (C,F,P,*,*), specifically,
1250  * element with row = cell and col = point should equal 1; all other should equal 0
1251  */
1252 
1253  outFields.resize(3, 1, 3, 3, 3);
1254  art::outerProductDataField<double>(outFields, ijkData_1b, ijkFields_1b);
1255 
1256  for(int cell = 0; cell < ijkData_1b.dimension(0); cell++){
1257  for(int field = 0; field < ijkFields_1b.dimension(0); field++){
1258  for(int point = 0; point < ijkData_1b.dimension(1); point++){
1259  for(int row = 0; row < ijkData_1b.dimension(2); row++){
1260  for(int col = 0; col < ijkData_1b.dimension(2); col++){
1261 
1262  // element with row = cell and col = point should equal 1; all other should equal 0
1263  if( (row == cell && col == point) ){
1264  if(outFields(cell, field, point, row, col) != 1.0) {
1265  *outStream << "\n\nINCORRECT outerProductDataField (3): computed value is "
1266  << outFields(cell, field, point, row, col) << " whereas correct value is 1.0";
1267  errorFlag++;
1268 
1269  }
1270  }
1271  else {
1272  if(outFields(cell, field, point, row, col) != 0.0) {
1273  *outStream << "\n\nINCORRECT outerProductDataField (4): computed value is "
1274  << outFields(cell, field, point, row, col) << " whereas correct value is 0.0";
1275  errorFlag++;
1276  }
1277  } // if
1278  }// col
1279  }// row
1280  }// point
1281  }// field
1282  }// cell
1283 
1284  *outStream \
1285  << "\n"
1286  << "===============================================================================\n"\
1287  << "| TEST 4.a: outerProductDataData operations: (C,P,D) and (C,P,D) |\n"\
1288  << "===============================================================================\n";
1289  /*
1290  * ijkData_1a jkiData_2a kijData_2a
1291  * i,i j,j k,k
1292  * j,j k,k i,i
1293  * k,k i,i j,j
1294  *
1295  * Expected results are stated with each test case.
1296  */
1297  outData.resize(3, 2, 3, 3);
1298 
1299  art::outerProductDataData<double>(outData, ijkData_1a, ijkData_1a);
1300  for(int cell = 0; cell < ijkData_1a.dimension(0); cell++){
1301  for(int point = 0; point < ijkData_1a.dimension(1); point++){
1302  for(int row = 0; row < ijkData_1a.dimension(2); row++){
1303  for(int col = 0; col < ijkData_1a.dimension(2); col++){
1304 
1305  // element with row = cell and col = cell should equal 1; all other should equal 0
1306  if( (row == cell && col == cell) ){
1307  if(outData(cell, point, row, col) != 1.0) {
1308  *outStream << "\n\nINCORRECT outerProductDataData (1): computed value is "
1309  << outData(cell, point, row, col) << " whereas correct value is 1.0";
1310  errorFlag++;
1311  }
1312  }
1313  else {
1314  if(outData(cell, point, row, col) != 0.0) {
1315  *outStream << "\n\nINCORRECT outerProductDataData (2): computed value is "
1316  << outData(cell, point, row, col) << " whereas correct value is 0.0";
1317  errorFlag++;
1318  }
1319  } // if
1320  }// col
1321  }// row
1322  }// point
1323  }// cell
1324 
1325  outData.initialize();
1326  art::outerProductDataData<double>(outData, ijkData_1a, jkiData_2a);
1327  for(int cell = 0; cell < ijkData_1a.dimension(0); cell++){
1328  for(int point = 0; point < ijkData_1a.dimension(1); point++){
1329  for(int row = 0; row < ijkData_1a.dimension(2); row++){
1330  for(int col = 0; col < ijkData_1a.dimension(2); col++){
1331 
1332  // element with row = cell and col = cell + 1 (mod 3) should equal 1; all other should equal 0
1333  if( (row == cell && col == (cell + 1) % 3) ){
1334  if(outData(cell, point, row, col) != 1.0) {
1335  *outStream << "\n\nINCORRECT outerProductDataData (3): computed value is "
1336  << outData(cell, point, row, col) << " whereas correct value is 1.0";
1337  errorFlag++;
1338  }
1339  }
1340  else {
1341  if(outData(cell, point, row, col) != 0.0) {
1342  *outStream << "\n\nINCORRECT outerProductDataData (4): computed value is "
1343  << outData(cell, point, row, col) << " whereas correct value is 0.0";
1344  errorFlag++;
1345  }
1346  } // if
1347  }// col
1348  }// row
1349  }// point
1350  }// cell
1351 
1352 
1353  outData.initialize();
1354  art::outerProductDataData<double>(outData, ijkData_1a, kijData_2a);
1355  for(int cell = 0; cell < ijkData_1a.dimension(0); cell++){
1356  for(int point = 0; point < ijkData_1a.dimension(1); point++){
1357  for(int row = 0; row < ijkData_1a.dimension(2); row++){
1358  for(int col = 0; col < ijkData_1a.dimension(2); col++){
1359 
1360  // element with row = cell and col = cell + 2 (mod 3) should equal 1; all other should equal 0
1361  if( (row == cell && col == (cell + 2) % 3) ){
1362  if(outData(cell, point, row, col) != 1.0) {
1363  *outStream << "\n\nINCORRECT outerProductDataData (5): computed value is "
1364  << outData(cell, point, row, col) << " whereas correct value is 1.0";
1365  errorFlag++;
1366  }
1367  }
1368  else {
1369  if(outData(cell, point, row, col) != 0.0) {
1370  *outStream << "\n\nINCORRECT outerProductDataData (6): computed value is "
1371  << outData(cell, point, row, col) << " whereas correct value is 0.0";
1372  errorFlag++;
1373  }
1374  } // if
1375  }// col
1376  }// row
1377  }// point
1378  }// cell
1379 
1380 
1381  *outStream \
1382  << "\n"
1383  << "===============================================================================\n"\
1384  << "| TEST 4.b: outerProductDataData operations: (C,P,D) and (P,D) |\n"\
1385  << "===============================================================================\n";
1386  /*
1387  * ijkData_1b ijkData_2b expected result in (C,P,D,D) array
1388  * i,i,i i,j,k (0,0) (0,1) (0,2)
1389  * j,j,j (1,0) (1,1) (1,2)
1390  * k,k,k (2,0) (2,1) (2,2)
1391  * Indicates the only non-zero element of (C,P,*,*), specifically,
1392  * element with row = cell and col = point should equal 1; all other should equal 0
1393  */
1394  outData.resize(3,3,3,3);
1395  art::outerProductDataData<double>(outData, ijkData_1b, ijkData_2b);
1396  for(int cell = 0; cell < ijkData_1b.dimension(0); cell++){
1397  for(int point = 0; point < ijkData_1b.dimension(1); point++){
1398  for(int row = 0; row < ijkData_1b.dimension(2); row++){
1399  for(int col = 0; col < ijkData_1b.dimension(2); col++){
1400 
1401  // element with row = cell and col = cell + 2 (mod 3) should equal 1; all other should equal 0
1402  if( (row == cell && col == point) ){
1403  if(outData(cell, point, row, col) != 1.0) {
1404  *outStream << "\n\nINCORRECT outerProductDataData (7): computed value is "
1405  << outData(cell, point, row, col) << " whereas correct value is 1.0";
1406  errorFlag++;
1407  }
1408  }
1409  else {
1410  if(outData(cell, point, row, col) != 0.0) {
1411  *outStream << "\n\nINCORRECT outerProductDataData (8): computed value is "
1412  << outData(cell, point, row, col) << " whereas correct value is 0.0";
1413  errorFlag++;
1414  }
1415  } // if
1416  }// col
1417  }// row
1418  }// point
1419  }// cell
1420 
1421  *outStream \
1422  << "\n"
1423  << "===============================================================================\n"\
1424  << "| TEST 5.a: matvecProductDataField operations: (C,P,D,D) and (C,F,P,D) |\n"\
1425  << "===============================================================================\n";
1426  /*
1427  * inputMat inputVecFields outFields
1428  * 1 1 1 0 0 0 0 1 -1 -1 0 3 0 0
1429  * -1 2 -1 -1 -2 -3 0 1 -1 1 0 0 6 2
1430  * 1 2 3 -2 6 -4 0 1 -1 -1 0 6 0 12
1431  */
1432 
1433  // (C,P,D,D)
1434  FieldContainer<double> inputMat(2,1,3,3);
1435  // cell 0
1436  inputMat(0,0,0,0) = 1.0; inputMat(0,0,0,1) = 1.0; inputMat(0,0,0,2) = 1.0;
1437  inputMat(0,0,1,0) =-1.0; inputMat(0,0,1,1) = 2.0; inputMat(0,0,1,2) =-1.0;
1438  inputMat(0,0,2,0) = 1.0; inputMat(0,0,2,1) = 2.0; inputMat(0,0,2,2) = 3.0;
1439  // cell 1
1440  inputMat(1,0,0,0) = 0.0; inputMat(1,0,0,1) = 0.0; inputMat(1,0,0,2) = 0.0;
1441  inputMat(1,0,1,0) =-1.0; inputMat(1,0,1,1) =-2.0; inputMat(1,0,1,2) =-3.0;
1442  inputMat(1,0,2,0) =-2.0; inputMat(1,0,2,1) = 6.0; inputMat(1,0,2,2) =-4.0;
1443 
1444  // (C,F,P,D)
1445  FieldContainer<double> inputVecFields(2,2,1,3);
1446  // cell 0; fields 0,1
1447  inputVecFields(0,0,0,0) = 0.0; inputVecFields(0,0,0,1) = 0.0; inputVecFields(0,0,0,2) = 0.0;
1448  inputVecFields(0,1,0,0) = 1.0; inputVecFields(0,1,0,1) = 1.0; inputVecFields(0,1,0,2) = 1.0;
1449  // cell 1; fields 0,1
1450  inputVecFields(1,0,0,0) =-1.0; inputVecFields(1,0,0,1) =-1.0; inputVecFields(1,0,0,2) =-1.0;
1451  inputVecFields(1,1,0,0) =-1.0; inputVecFields(1,1,0,1) = 1.0; inputVecFields(1,1,0,2) =-1.0;
1452 
1453  // (C,F,P,D) - true
1454  FieldContainer<double> outFieldsCorrect(2,2,1,3);
1455  // cell 0; fields 0,1
1456  outFieldsCorrect(0,0,0,0) = 0.0; outFieldsCorrect(0,0,0,1) = 0.0; outFieldsCorrect(0,0,0,2) = 0.0;
1457  outFieldsCorrect(0,1,0,0) = 3.0; outFieldsCorrect(0,1,0,1) = 0.0; outFieldsCorrect(0,1,0,2) = 6.0;
1458  // cell 1; fields 0,1
1459  outFieldsCorrect(1,0,0,0) = 0.0; outFieldsCorrect(1,0,0,1) = 6.0; outFieldsCorrect(1,0,0,2) = 0.0;
1460  outFieldsCorrect(1,1,0,0) = 0.0; outFieldsCorrect(1,1,0,1) = 2.0; outFieldsCorrect(1,1,0,2) = 12.0;
1461 
1462  // (C,F,P,D)
1463  outFields.resize(2,2,1,3);
1464  art::matvecProductDataField<double>(outFields, inputMat, inputVecFields);
1465 
1466  // test loop
1467  for(int cell = 0; cell < outFields.dimension(0); cell++){
1468  for(int field = 0; field < outFields.dimension(1); field++){
1469  for(int point = 0; point < outFields.dimension(2); point++){
1470  for(int row = 0; row < outFields.dimension(3); row++){
1471  if(outFields(cell, field, point, row) != outFieldsCorrect(cell, field, point, row)) {
1472  *outStream << "\n\nINCORRECT matvecProductDataField (1): \n value at multi-index ("
1473  << cell << "," << field << "," << point << "," << row << ") = "
1474  << outFields(cell, field, point, row) << " but correct value is "
1475  << outFieldsCorrect(cell, field, point, row) <<"\n";
1476  errorFlag++;
1477  }
1478  }//row
1479  }// point
1480  }// field
1481  }// cell
1482 
1483 
1484  *outStream \
1485  << "\n"
1486  << "===============================================================================\n"\
1487  << "| TEST 5.b: matvecProductDataField operations: (C,P,D,D) and (F,P,D) |\n"\
1488  << "===============================================================================\n";
1489  /*
1490  * inputMat inputVecFields outFields
1491  * 1 1 1 0 0 0 0 1 -1 -1 0 3 -3 -1 0 0 0 0
1492  * -1 2 -1 -1 -2 -3 0 1 -1 1 0 0 0 4 0 -6 6 2
1493  * 1 2 3 -2 6 -4 0 1 -1 -1 0 6 -6 -2 0 0 0 12
1494  * Use the same 4 vector fields as above but formatted as (F,P,D) array
1495  */
1496  // (C,F,P,D)
1497  inputVecFields.resize(4,1,3);
1498  // fields 0,1,2,3
1499  inputVecFields(0,0,0) = 0.0; inputVecFields(0,0,1) = 0.0; inputVecFields(0,0,2) = 0.0;
1500  inputVecFields(1,0,0) = 1.0; inputVecFields(1,0,1) = 1.0; inputVecFields(1,0,2) = 1.0;
1501  inputVecFields(2,0,0) =-1.0; inputVecFields(2,0,1) =-1.0; inputVecFields(2,0,2) =-1.0;
1502  inputVecFields(3,0,0) =-1.0; inputVecFields(3,0,1) = 1.0; inputVecFields(3,0,2) =-1.0;
1503 
1504  // (C,F,P,D) - true
1505  outFieldsCorrect.resize(2,4,1,3);
1506  // cell 0; fields 0,1,2,3
1507  outFieldsCorrect(0,0,0,0) = 0.0; outFieldsCorrect(0,0,0,1) = 0.0; outFieldsCorrect(0,0,0,2) = 0.0;
1508  outFieldsCorrect(0,1,0,0) = 3.0; outFieldsCorrect(0,1,0,1) = 0.0; outFieldsCorrect(0,1,0,2) = 6.0;
1509  outFieldsCorrect(0,2,0,0) =-3.0; outFieldsCorrect(0,2,0,1) = 0.0; outFieldsCorrect(0,2,0,2) =-6.0;
1510  outFieldsCorrect(0,3,0,0) =-1.0; outFieldsCorrect(0,3,0,1) = 4.0; outFieldsCorrect(0,3,0,2) =-2.0;
1511  // cell 1; fields 0,1,2,3
1512  outFieldsCorrect(1,0,0,0) = 0.0; outFieldsCorrect(1,0,0,1) = 0.0; outFieldsCorrect(1,0,0,2) = 0.0;
1513  outFieldsCorrect(1,1,0,0) = 0.0; outFieldsCorrect(1,1,0,1) =-6.0; outFieldsCorrect(1,1,0,2) = 0.0;
1514  outFieldsCorrect(1,2,0,0) = 0.0; outFieldsCorrect(1,2,0,1) = 6.0; outFieldsCorrect(1,2,0,2) = 0.0;
1515  outFieldsCorrect(1,3,0,0) = 0.0; outFieldsCorrect(1,3,0,1) = 2.0; outFieldsCorrect(1,3,0,2) =12.0;
1516 
1517  // (C,F,P,D)
1518  outFields.resize(2,4,1,3);
1519  art::matvecProductDataField<double>(outFields, inputMat, inputVecFields);
1520 
1521  // test loop
1522  for(int cell = 0; cell < outFields.dimension(0); cell++){
1523  for(int field = 0; field < outFields.dimension(1); field++){
1524  for(int point = 0; point < outFields.dimension(2); point++){
1525  for(int row = 0; row < outFields.dimension(3); row++){
1526  if(outFields(cell, field, point, row) != outFieldsCorrect(cell, field, point, row)) {
1527  *outStream << "\n\nINCORRECT matvecProductDataField (2): \n value at multi-index ("
1528  << cell << "," << field << "," << point << "," << row << ") = "
1529  << outFields(cell, field, point, row) << " but correct value is "
1530  << outFieldsCorrect(cell, field, point, row) <<"\n";
1531  errorFlag++;
1532  }
1533  }//row
1534  }// point
1535  }// field
1536  }// cell
1537 
1538 
1539  *outStream \
1540  << "\n"
1541  << "===============================================================================\n"\
1542  << "| TEST 5.c: matvecProductDataField random tests: branch inputFields(C,F,P,D) |\n"\
1543  << "===============================================================================\n";
1544  /*
1545  * d1 is spatial dimension and should be 1, 2 or 3. If d1>3, the RealSpaceTools function 'inverse' will fail
1546  */
1547  {// test 5.c scope
1548  int c=5, p=9, f=7, d1=3;
1549  double zero = INTREPID_TOL*10000.0;
1550 
1551  FieldContainer<double> in_c_f_p_d(c, f, p, d1);
1552  FieldContainer<double> out_c_f_p_d(c, f, p, d1);
1553  FieldContainer<double> outi_c_f_p_d(c, f, p, d1);
1554 
1555  FieldContainer<double> data_c_p(c, p);
1556  FieldContainer<double> datainv_c_p(c, p);
1557  FieldContainer<double> data_c_1(c, 1);
1558  FieldContainer<double> datainv_c_1(c, 1);
1559  FieldContainer<double> data_c_p_d(c, p, d1);
1560  FieldContainer<double> datainv_c_p_d(c, p, d1);
1561  FieldContainer<double> data_c_1_d(c, 1, d1);
1562  FieldContainer<double> datainv_c_1_d(c, 1, d1);
1563  FieldContainer<double> data_c_p_d_d(c, p, d1, d1);
1564  FieldContainer<double> datainv_c_p_d_d(c, p, d1, d1);
1565  FieldContainer<double> datainvtrn_c_p_d_d(c, p, d1, d1);
1566  FieldContainer<double> data_c_1_d_d(c, 1, d1, d1);
1567  FieldContainer<double> datainv_c_1_d_d(c, 1, d1, d1);
1568  FieldContainer<double> datainvtrn_c_1_d_d(c, 1, d1, d1);
1569  /***********************************************************************************************
1570  * Constant diagonal tensor: inputData(C,P) *
1571  **********************************************************************************************/
1572  for (int i=0; i<in_c_f_p_d.size(); i++) {
1573  in_c_f_p_d[i] = Teuchos::ScalarTraits<double>::random();
1574  }
1575  for (int i=0; i<data_c_p.size(); i++) {
1576  data_c_p[i] = Teuchos::ScalarTraits<double>::random();
1577  datainv_c_p[i] = 1.0 / data_c_p[i];
1578  }
1579  for (int i=0; i<data_c_1.size(); i++) {
1580  data_c_1[i] = Teuchos::ScalarTraits<double>::random();
1581  datainv_c_1[i] = 1.0 / data_c_1[i];
1582  }
1583  // Tensor values vary by point:
1584  art::matvecProductDataField<double>(out_c_f_p_d, data_c_p, in_c_f_p_d);
1585  art::matvecProductDataField<double>(out_c_f_p_d, datainv_c_p, out_c_f_p_d);
1586  rst::subtract(&out_c_f_p_d[0], &in_c_f_p_d[0], out_c_f_p_d.size());
1587  if (rst::vectorNorm(&out_c_f_p_d[0], out_c_f_p_d.size(), NORM_ONE) > zero) {
1588  *outStream << "\n\nINCORRECT matvecProductDataField (3): check scalar inverse property\n\n";
1589  errorFlag = -1000;
1590  }
1591  // Tensor values do not vary by point:
1592  art::matvecProductDataField<double>(out_c_f_p_d, data_c_1, in_c_f_p_d);
1593  art::matvecProductDataField<double>(out_c_f_p_d, datainv_c_1, out_c_f_p_d);
1594  rst::subtract(&out_c_f_p_d[0], &in_c_f_p_d[0], out_c_f_p_d.size());
1595  if (rst::vectorNorm(&out_c_f_p_d[0], out_c_f_p_d.size(), NORM_ONE) > zero) {
1596  *outStream << "\n\nINCORRECT matvecProductDataField (4): check scalar inverse property\n\n";
1597  errorFlag = -1000;
1598  }
1599  /***********************************************************************************************
1600  * Non-onstant diagonal tensor: inputData(C,P,D) *
1601  **********************************************************************************************/
1602  for (int i=0; i<in_c_f_p_d.size(); i++) {
1603  in_c_f_p_d[i] = Teuchos::ScalarTraits<double>::random();
1604  }
1605  for (int i=0; i<data_c_p_d.size(); i++) {
1606  data_c_p_d[i] = Teuchos::ScalarTraits<double>::random();
1607  datainv_c_p_d[i] = 1.0 / data_c_p_d[i];
1608  }
1609  for (int i=0; i<data_c_1_d.size(); i++) {
1610  data_c_1_d[i] = Teuchos::ScalarTraits<double>::random();
1611  datainv_c_1_d[i] = 1.0 / data_c_1_d[i];
1612  }
1613  // Tensor values vary by point:
1614  art::matvecProductDataField<double>(out_c_f_p_d, data_c_p_d, in_c_f_p_d);
1615  art::matvecProductDataField<double>(out_c_f_p_d, datainv_c_p_d, out_c_f_p_d);
1616  rst::subtract(&out_c_f_p_d[0], &in_c_f_p_d[0], out_c_f_p_d.size());
1617  if (rst::vectorNorm(&out_c_f_p_d[0], out_c_f_p_d.size(), NORM_ONE) > zero) {
1618  *outStream << "\n\nINCORRECT matvecProductDataField (5): check scalar inverse property\n\n";
1619  errorFlag = -1000;
1620  }
1621  // Tensor values do not vary by point
1622  art::matvecProductDataField<double>(out_c_f_p_d, data_c_1_d, in_c_f_p_d);
1623  art::matvecProductDataField<double>(out_c_f_p_d, datainv_c_1_d, out_c_f_p_d);
1624  rst::subtract(&out_c_f_p_d[0], &in_c_f_p_d[0], out_c_f_p_d.size());
1625  if (rst::vectorNorm(&out_c_f_p_d[0], out_c_f_p_d.size(), NORM_ONE) > zero) {
1626  *outStream << "\n\nINCORRECT matvecProductDataField (6): check scalar inverse property\n\n";
1627  errorFlag = -1000;
1628  }
1629  /***********************************************************************************************
1630  * Full tensor: inputData(C,P,D,D) *
1631  **********************************************************************************************/
1632  for (int i=0; i<in_c_f_p_d.size(); i++) {
1633  in_c_f_p_d[i] = Teuchos::ScalarTraits<double>::random();
1634  }
1635  for (int ic=0; ic < c; ic++) {
1636  for (int ip=0; ip < p; ip++) {
1637  for (int i=0; i<d1*d1; i++) {
1638  (&data_c_p_d_d(ic, ip, 0, 0))[i] = Teuchos::ScalarTraits<double>::random();
1639  }
1640  RealSpaceTools<double>::inverse(&datainv_c_p_d_d(ic, ip, 0, 0), &data_c_p_d_d(ic, ip, 0, 0), d1);
1641  }
1642  }
1643  for (int ic=0; ic < c; ic++) {
1644  for (int ip=0; ip < 1; ip++) {
1645  for (int i=0; i<d1*d1; i++) {
1646  (&data_c_1_d_d(ic, 0, 0, 0))[i] = Teuchos::ScalarTraits<double>::random();
1647  }
1648  RealSpaceTools<double>::inverse(&datainv_c_1_d_d(ic, 0, 0, 0), &data_c_1_d_d(ic, 0, 0, 0), d1);
1649  }
1650  }
1651  // Tensor values vary by point: test "N" and "T" options (no transpose/transpose)
1652  art::matvecProductDataField<double>(out_c_f_p_d, data_c_p_d_d, in_c_f_p_d);
1653  art::matvecProductDataField<double>(outi_c_f_p_d, datainv_c_p_d_d, out_c_f_p_d);
1654  rst::subtract(&outi_c_f_p_d[0], &in_c_f_p_d[0], outi_c_f_p_d.size());
1655  if (rst::vectorNorm(&outi_c_f_p_d[0], outi_c_f_p_d.size(), NORM_ONE) > zero) {
1656  *outStream << "\n\nINCORRECT matvecProductDataField (7): check matrix inverse property\n\n";
1657  errorFlag = -1000;
1658  }
1659  art::matvecProductDataField<double>(out_c_f_p_d, data_c_p_d_d, in_c_f_p_d, 't');
1660  art::matvecProductDataField<double>(outi_c_f_p_d, datainv_c_p_d_d, out_c_f_p_d, 't');
1661  rst::subtract(&outi_c_f_p_d[0], &in_c_f_p_d[0], outi_c_f_p_d.size());
1662  if (rst::vectorNorm(&outi_c_f_p_d[0], outi_c_f_p_d.size(), NORM_ONE) > zero) {
1663  *outStream << "\n\nINCORRECT matvecProductDataField (8): check matrix inverse property, w/ double transpose\n\n";
1664  errorFlag = -1000;
1665  }
1666  // Tensor values do not vary by point: test "N" and "T" options (no transpose/transpose)
1667  art::matvecProductDataField<double>(out_c_f_p_d, data_c_1_d_d, in_c_f_p_d);
1668  art::matvecProductDataField<double>(outi_c_f_p_d, datainv_c_1_d_d, out_c_f_p_d);
1669  rst::subtract(&outi_c_f_p_d[0], &in_c_f_p_d[0], outi_c_f_p_d.size());
1670  if (rst::vectorNorm(&outi_c_f_p_d[0], outi_c_f_p_d.size(), NORM_ONE) > zero) {
1671  *outStream << "\n\nINCORRECT matvecProductDataField (9): check matrix inverse property\n\n";
1672  errorFlag = -1000;
1673  }
1674  art::matvecProductDataField<double>(out_c_f_p_d, data_c_1_d_d, in_c_f_p_d, 't');
1675  art::matvecProductDataField<double>(outi_c_f_p_d, datainv_c_1_d_d, out_c_f_p_d, 't');
1676  rst::subtract(&outi_c_f_p_d[0], &in_c_f_p_d[0], outi_c_f_p_d.size());
1677  if (rst::vectorNorm(&outi_c_f_p_d[0], outi_c_f_p_d.size(), NORM_ONE) > zero) {
1678  *outStream << "\n\nINCORRECT matvecProductDataField (10): check matrix inverse property, w/ double transpose\n\n";
1679  errorFlag = -1000;
1680  }
1681  /***********************************************************************************************
1682  * Full tensor: inputData(C,P,D,D) test inverse transpose *
1683  **********************************************************************************************/
1684  for (int i=0; i<in_c_f_p_d.size(); i++) {
1685  in_c_f_p_d[i] = Teuchos::ScalarTraits<double>::random();
1686  }
1687  for (int ic=0; ic < c; ic++) {
1688  for (int ip=0; ip < p; ip++) {
1689  for (int i=0; i<d1*d1; i++) {
1690  (&data_c_p_d_d(ic, ip, 0, 0))[i] = Teuchos::ScalarTraits<double>::random();
1691  }
1692  RealSpaceTools<double>::inverse(&datainv_c_p_d_d(ic, ip, 0, 0), &data_c_p_d_d(ic, ip, 0, 0), d1);
1693  RealSpaceTools<double>::transpose(&datainvtrn_c_p_d_d(ic, ip, 0, 0), &datainv_c_p_d_d(ic, ip, 0, 0), d1);
1694  }
1695  }
1696  for (int ic=0; ic < c; ic++) {
1697  for (int ip=0; ip < 1; ip++) {
1698  for (int i=0; i<d1*d1; i++) {
1699  (&data_c_1_d_d(ic, 0, 0, 0))[i] = Teuchos::ScalarTraits<double>::random();
1700  }
1701  RealSpaceTools<double>::inverse(&datainv_c_1_d_d(ic, 0, 0, 0), &data_c_1_d_d(ic, 0, 0, 0), d1);
1702  RealSpaceTools<double>::transpose(&datainvtrn_c_1_d_d(ic, 0, 0, 0), &datainv_c_1_d_d(ic, 0, 0, 0), d1);
1703  }
1704  }
1705  // Tensor values vary by point:
1706  art::matvecProductDataField<double>(out_c_f_p_d, data_c_p_d_d, in_c_f_p_d);
1707  art::matvecProductDataField<double>(outi_c_f_p_d, datainvtrn_c_p_d_d, out_c_f_p_d, 't');
1708  rst::subtract(&outi_c_f_p_d[0], &in_c_f_p_d[0], outi_c_f_p_d.size());
1709  if (rst::vectorNorm(&outi_c_f_p_d[0], outi_c_f_p_d.size(), NORM_ONE) > zero) {
1710  *outStream << "\n\nINCORRECT matvecProductDataField (11): check matrix inverse transpose property\n\n";
1711  errorFlag = -1000;
1712  }
1713  // Tensor values do not vary by point
1714  art::matvecProductDataField<double>(out_c_f_p_d, data_c_1_d_d, in_c_f_p_d);
1715  art::matvecProductDataField<double>(outi_c_f_p_d, datainvtrn_c_1_d_d, out_c_f_p_d, 't');
1716  rst::subtract(&outi_c_f_p_d[0], &in_c_f_p_d[0], outi_c_f_p_d.size());
1717  if (rst::vectorNorm(&outi_c_f_p_d[0], outi_c_f_p_d.size(), NORM_ONE) > zero) {
1718  *outStream << "\n\nINCORRECT matvecProductDataField (12): check matrix inverse transpose property\n\n";
1719  errorFlag = -1000;
1720  }
1721  }// test 5.c scope
1722 
1723 
1724  *outStream \
1725  << "\n"
1726  << "===============================================================================\n"\
1727  << "| TEST 5.d: matvecProductDataField random tests: branch inputFields(F,P,D) |\n"\
1728  << "===============================================================================\n";
1729  /*
1730  * d1 is the spatial dimension and should be 1, 2 or 3. If d1>3, the RealSpaceTools function 'inverse' will fail
1731  */
1732  {// test 5.d scope
1733  int c=5, p=9, f=7, d1=3;
1734  double zero = INTREPID_TOL*10000.0;
1735 
1736  FieldContainer<double> in_f_p_d(f, p, d1);
1737  FieldContainer<double> in_c_f_p_d(c, f, p, d1);
1738  FieldContainer<double> data_c_p(c, p);
1739  FieldContainer<double> datainv_c_p(c, p);
1740  FieldContainer<double> data_c_1(c, 1);
1741  FieldContainer<double> datainv_c_1(c, 1);
1742  FieldContainer<double> data_c_p_d(c, p, d1);
1743  FieldContainer<double> datainv_c_p_d(c, p, d1);
1744  FieldContainer<double> data_c_1_d(c, 1, d1);
1745  FieldContainer<double> datainv_c_1_d(c, 1, d1);
1746  FieldContainer<double> data_c_p_d_d(c, p, d1, d1);
1747  FieldContainer<double> datainv_c_p_d_d(c, p, d1, d1);
1748  FieldContainer<double> datainvtrn_c_p_d_d(c, p, d1, d1);
1749  FieldContainer<double> data_c_1_d_d(c, 1, d1, d1);
1750  FieldContainer<double> datainv_c_1_d_d(c, 1, d1, d1);
1751  FieldContainer<double> datainvtrn_c_1_d_d(c, 1, d1, d1);
1752  FieldContainer<double> data_c_p_one(c, p);
1753  FieldContainer<double> data_c_1_one(c, 1);
1754  FieldContainer<double> out_c_f_p_d(c, f, p, d1);
1755  FieldContainer<double> outi_c_f_p_d(c, f, p, d1);
1756  /***********************************************************************************************
1757  * Constant diagonal tensor: inputData(C,P) *
1758  **********************************************************************************************/
1759  for (int i=0; i<in_f_p_d.size(); i++) {
1760  in_f_p_d[i] = Teuchos::ScalarTraits<double>::random();
1761  }
1762  for (int i=0; i<data_c_p.size(); i++) {
1763  data_c_p[i] = Teuchos::ScalarTraits<double>::random();
1764  datainv_c_p[i] = 1.0 / data_c_p[i];
1765  data_c_p_one[i] = 1.0;
1766  }
1767  for (int i=0; i<data_c_1.size(); i++) {
1768  data_c_1[i] = Teuchos::ScalarTraits<double>::random();
1769  datainv_c_1[i] = 1.0 / data_c_1[i];
1770  }
1771  // Tensor values vary by point
1772  art::matvecProductDataField<double>(in_c_f_p_d, data_c_p_one, in_f_p_d);
1773  art::matvecProductDataField<double>(out_c_f_p_d, data_c_p, in_f_p_d);
1774  art::matvecProductDataField<double>(outi_c_f_p_d, datainv_c_p, out_c_f_p_d);
1775  rst::subtract(&outi_c_f_p_d[0], &in_c_f_p_d[0], outi_c_f_p_d.size());
1776  if (rst::vectorNorm(&outi_c_f_p_d[0], outi_c_f_p_d.size(), NORM_ONE) > zero) {
1777  *outStream << "\n\nINCORRECT matvecProductDataField (13): check scalar inverse property\n\n";
1778  errorFlag = -1000;
1779  }
1780  // Tensor values do not vary by point
1781  art::matvecProductDataField<double>(in_c_f_p_d, data_c_p_one, in_f_p_d);
1782  art::matvecProductDataField<double>(out_c_f_p_d, data_c_1, in_f_p_d);
1783  art::matvecProductDataField<double>(outi_c_f_p_d, datainv_c_1, out_c_f_p_d);
1784  rst::subtract(&outi_c_f_p_d[0], &in_c_f_p_d[0], outi_c_f_p_d.size());
1785  if (rst::vectorNorm(&outi_c_f_p_d[0], outi_c_f_p_d.size(), NORM_ONE) > zero) {
1786  *outStream << "\n\nINCORRECT matvecProductDataField (14): check scalar inverse property\n\n";
1787  errorFlag = -1000;
1788  }
1789  /***********************************************************************************************
1790  * Non-constant diagonal tensor: inputData(C,P,D) *
1791  **********************************************************************************************/
1792  for (int i=0; i<in_f_p_d.size(); i++) {
1793  in_f_p_d[i] = Teuchos::ScalarTraits<double>::random();
1794  }
1795  for (int i=0; i<data_c_p_d.size(); i++) {
1796  data_c_p_d[i] = Teuchos::ScalarTraits<double>::random();
1797  datainv_c_p_d[i] = 1.0 / data_c_p_d[i];
1798  }
1799  for (int i=0; i<data_c_1_d.size(); i++) {
1800  data_c_1_d[i] = Teuchos::ScalarTraits<double>::random();
1801  datainv_c_1_d[i] = 1.0 / data_c_1_d[i];
1802  }
1803  // Tensor values vary by point:
1804  art::matvecProductDataField<double>(in_c_f_p_d, data_c_p_one, in_f_p_d);
1805  art::matvecProductDataField<double>(out_c_f_p_d, data_c_p_d, in_f_p_d);
1806  art::matvecProductDataField<double>(outi_c_f_p_d, datainv_c_p_d, out_c_f_p_d);
1807  rst::subtract(&outi_c_f_p_d[0], &in_c_f_p_d[0], outi_c_f_p_d.size());
1808  if (rst::vectorNorm(&outi_c_f_p_d[0], outi_c_f_p_d.size(), NORM_ONE) > zero) {
1809  *outStream << "\n\nINCORRECT matvecProductDataField (15): check scalar inverse property\n\n";
1810  errorFlag = -1000;
1811  }
1812  // Tensor values do not vary by point:
1813  art::matvecProductDataField<double>(in_c_f_p_d, data_c_p_one, in_f_p_d);
1814  art::matvecProductDataField<double>(out_c_f_p_d, data_c_1_d, in_f_p_d);
1815  art::matvecProductDataField<double>(outi_c_f_p_d, datainv_c_1_d, out_c_f_p_d);
1816  rst::subtract(&outi_c_f_p_d[0], &in_c_f_p_d[0], outi_c_f_p_d.size());
1817  if (rst::vectorNorm(&outi_c_f_p_d[0], outi_c_f_p_d.size(), NORM_ONE) > zero) {
1818  *outStream << "\n\nINCORRECT matvecProductDataField (16): check scalar inverse property\n\n";
1819  errorFlag = -1000;
1820  }
1821  /***********************************************************************************************
1822  * Full tensor: inputData(C,P,D,D) *
1823  **********************************************************************************************/
1824  for (int i=0; i<in_f_p_d.size(); i++) {
1825  in_f_p_d[i] = Teuchos::ScalarTraits<double>::random();
1826  }
1827  for (int ic=0; ic < c; ic++) {
1828  for (int ip=0; ip < p; ip++) {
1829  for (int i=0; i<d1*d1; i++) {
1830  (&data_c_p_d_d(ic, ip, 0, 0))[i] = Teuchos::ScalarTraits<double>::random();
1831  }
1832  RealSpaceTools<double>::inverse(&datainv_c_p_d_d(ic, ip, 0, 0), &data_c_p_d_d(ic, ip, 0, 0), d1);
1833  }
1834  }
1835  for (int ic=0; ic < c; ic++) {
1836  for (int ip=0; ip < 1; ip++) {
1837  for (int i=0; i<d1*d1; i++) {
1838  (&data_c_1_d_d(ic, 0, 0, 0))[i] = Teuchos::ScalarTraits<double>::random();
1839  }
1840  RealSpaceTools<double>::inverse(&datainv_c_1_d_d(ic, 0, 0, 0), &data_c_1_d_d(ic, 0, 0, 0), d1);
1841  }
1842  }
1843  // Tensor values vary by point: test "N" and "T" (no-transpose/transpose) options
1844  art::matvecProductDataField<double>(in_c_f_p_d, data_c_p_one, in_f_p_d);
1845  art::matvecProductDataField<double>(out_c_f_p_d, data_c_p_d_d, in_f_p_d);
1846  art::matvecProductDataField<double>(outi_c_f_p_d, datainv_c_p_d_d, out_c_f_p_d);
1847  rst::subtract(&outi_c_f_p_d[0], &in_c_f_p_d[0], outi_c_f_p_d.size());
1848  if (rst::vectorNorm(&outi_c_f_p_d[0], outi_c_f_p_d.size(), NORM_ONE) > zero) {
1849  *outStream << "\n\nINCORRECT matvecProductDataField (17): check matrix inverse property\n\n";
1850  errorFlag = -1000;
1851  }
1852  art::matvecProductDataField<double>(in_c_f_p_d, data_c_p_one, in_f_p_d);
1853  art::matvecProductDataField<double>(out_c_f_p_d, data_c_p_d_d, in_f_p_d, 't');
1854  art::matvecProductDataField<double>(outi_c_f_p_d, datainv_c_p_d_d, out_c_f_p_d, 't');
1855  rst::subtract(&outi_c_f_p_d[0], &in_c_f_p_d[0], outi_c_f_p_d.size());
1856  if (rst::vectorNorm(&outi_c_f_p_d[0], outi_c_f_p_d.size(), NORM_ONE) > zero) {
1857  *outStream << "\n\nINCORRECT matvecProductDataField (18): check matrix inverse property, w/ double transpose\n\n";
1858  errorFlag = -1000;
1859  }
1860  // Tensor values do not vary by point: test "N" and "T" (no-transpose/transpose) options
1861  art::matvecProductDataField<double>(in_c_f_p_d, data_c_p_one, in_f_p_d);
1862  art::matvecProductDataField<double>(out_c_f_p_d, data_c_1_d_d, in_f_p_d);
1863  art::matvecProductDataField<double>(outi_c_f_p_d, datainv_c_1_d_d, out_c_f_p_d);
1864  rst::subtract(&outi_c_f_p_d[0], &in_c_f_p_d[0], outi_c_f_p_d.size());
1865  if (rst::vectorNorm(&outi_c_f_p_d[0], outi_c_f_p_d.size(), NORM_ONE) > zero) {
1866  *outStream << "\n\nINCORRECT matvecProductDataField (19): check matrix inverse property\n\n";
1867  errorFlag = -1000;
1868  }
1869  art::matvecProductDataField<double>(in_c_f_p_d, data_c_p_one, in_f_p_d);
1870  art::matvecProductDataField<double>(out_c_f_p_d, data_c_1_d_d, in_f_p_d, 't');
1871  art::matvecProductDataField<double>(outi_c_f_p_d, datainv_c_1_d_d, out_c_f_p_d, 't');
1872  rst::subtract(&outi_c_f_p_d[0], &in_c_f_p_d[0], outi_c_f_p_d.size());
1873  if (rst::vectorNorm(&outi_c_f_p_d[0], outi_c_f_p_d.size(), NORM_ONE) > zero) {
1874  *outStream << "\n\nINCORRECT matvecProductDataField (20): check matrix inverse property, w/ double transpose\n\n";
1875  errorFlag = -1000;
1876  }
1877  /***********************************************************************************************
1878  * Full tensor: inputData(C,P,D,D) test inverse transpose *
1879  **********************************************************************************************/
1880  for (int i=0; i<in_f_p_d.size(); i++) {
1881  in_f_p_d[i] = Teuchos::ScalarTraits<double>::random();
1882  }
1883  for (int ic=0; ic < c; ic++) {
1884  for (int ip=0; ip < p; ip++) {
1885  for (int i=0; i<d1*d1; i++) {
1886  (&data_c_p_d_d(ic, ip, 0, 0))[i] = Teuchos::ScalarTraits<double>::random();
1887  }
1888  RealSpaceTools<double>::inverse(&datainv_c_p_d_d(ic, ip, 0, 0), &data_c_p_d_d(ic, ip, 0, 0), d1);
1889  RealSpaceTools<double>::transpose(&datainvtrn_c_p_d_d(ic, ip, 0, 0), &datainv_c_p_d_d(ic, ip, 0, 0), d1);
1890  }
1891  }
1892  for (int ic=0; ic < c; ic++) {
1893  for (int ip=0; ip < 1; ip++) {
1894  for (int i=0; i<d1*d1; i++) {
1895  (&data_c_1_d_d(ic, 0, 0, 0))[i] = Teuchos::ScalarTraits<double>::random();
1896  }
1897  RealSpaceTools<double>::inverse(&datainv_c_1_d_d(ic, 0, 0, 0), &data_c_1_d_d(ic, 0, 0, 0), d1);
1898  RealSpaceTools<double>::transpose(&datainvtrn_c_1_d_d(ic, 0, 0, 0), &datainv_c_1_d_d(ic, 0, 0, 0), d1);
1899  }
1900  }
1901  // Tensor values vary by point:
1902  art::matvecProductDataField<double>(in_c_f_p_d, data_c_p_one, in_f_p_d);
1903  art::matvecProductDataField<double>(out_c_f_p_d, data_c_p_d_d, in_f_p_d);
1904  art::matvecProductDataField<double>(outi_c_f_p_d, datainvtrn_c_p_d_d, out_c_f_p_d, 't');
1905  rst::subtract(&outi_c_f_p_d[0], &in_c_f_p_d[0], outi_c_f_p_d.size());
1906  if (rst::vectorNorm(&outi_c_f_p_d[0], outi_c_f_p_d.size(), NORM_ONE) > zero) {
1907  *outStream << "\n\nINCORRECT matvecProductDataField (21): check matrix inverse transpose property\n\n";
1908  errorFlag = -1000;
1909  }
1910  // Tensor values do not vary by point:
1911  art::matvecProductDataField<double>(in_c_f_p_d, data_c_p_one, in_f_p_d);
1912  art::matvecProductDataField<double>(out_c_f_p_d, data_c_1_d_d, in_f_p_d);
1913  art::matvecProductDataField<double>(outi_c_f_p_d, datainvtrn_c_1_d_d, out_c_f_p_d, 't');
1914  rst::subtract(&outi_c_f_p_d[0], &in_c_f_p_d[0], outi_c_f_p_d.size());
1915  if (rst::vectorNorm(&outi_c_f_p_d[0], outi_c_f_p_d.size(), NORM_ONE) > zero) {
1916  *outStream << "\n\nINCORRECT matvecProductDataField (22): check matrix inverse transpose property\n\n";
1917  errorFlag = -1000;
1918  }
1919  }// test 5.d scope
1920 
1921 
1922 
1923  *outStream \
1924  << "\n"
1925  << "===============================================================================\n"\
1926  << "| TEST 6.a: matvecProductDataData operations: (C,P,D,D) and (C,P,D) |\n"\
1927  << "===============================================================================\n";
1928  /*
1929  * inputMat inputVecFields outFields
1930  * 1 1 1 0 0 0 1 1 1 0 0 0 0 1 -1 -1 0 0 -3 0
1931  * -1 2 -1 -1 -2 -3 -1 2 -1 -1 -2 -3 0 1 -1 1 0 -6 0 2
1932  * 1 2 3 -2 6 -4 1 2 3 -2 6 -4 0 1 -1 -1 0 0 -6 12
1933  */
1934 
1935  // (C,P,D,D)
1936  inputMat.resize(4,1,3,3);
1937  // cell 0
1938  inputMat(0,0,0,0) = 1.0; inputMat(0,0,0,1) = 1.0; inputMat(0,0,0,2) = 1.0;
1939  inputMat(0,0,1,0) =-1.0; inputMat(0,0,1,1) = 2.0; inputMat(0,0,1,2) =-1.0;
1940  inputMat(0,0,2,0) = 1.0; inputMat(0,0,2,1) = 2.0; inputMat(0,0,2,2) = 3.0;
1941  // cell 1
1942  inputMat(1,0,0,0) = 0.0; inputMat(1,0,0,1) = 0.0; inputMat(1,0,0,2) = 0.0;
1943  inputMat(1,0,1,0) =-1.0; inputMat(1,0,1,1) =-2.0; inputMat(1,0,1,2) =-3.0;
1944  inputMat(1,0,2,0) =-2.0; inputMat(1,0,2,1) = 6.0; inputMat(1,0,2,2) =-4.0;
1945  // cell 2
1946  inputMat(2,0,0,0) = 1.0; inputMat(2,0,0,1) = 1.0; inputMat(2,0,0,2) = 1.0;
1947  inputMat(2,0,1,0) =-1.0; inputMat(2,0,1,1) = 2.0; inputMat(2,0,1,2) =-1.0;
1948  inputMat(2,0,2,0) = 1.0; inputMat(2,0,2,1) = 2.0; inputMat(2,0,2,2) = 3.0;
1949  // cell 3
1950  inputMat(3,0,0,0) = 0.0; inputMat(3,0,0,1) = 0.0; inputMat(3,0,0,2) = 0.0;
1951  inputMat(3,0,1,0) =-1.0; inputMat(3,0,1,1) =-2.0; inputMat(3,0,1,2) =-3.0;
1952  inputMat(3,0,2,0) =-2.0; inputMat(3,0,2,1) = 6.0; inputMat(3,0,2,2) =-4.0;
1953 
1954  // (C,P,D)
1955  inputVecFields.resize(4,1,3);
1956  inputVecFields(0,0,0) = 0.0; inputVecFields(0,0,1) = 0.0; inputVecFields(0,0,2) = 0.0;
1957  inputVecFields(1,0,0) = 1.0; inputVecFields(1,0,1) = 1.0; inputVecFields(1,0,2) = 1.0;
1958  inputVecFields(2,0,0) =-1.0; inputVecFields(2,0,1) =-1.0; inputVecFields(2,0,2) =-1.0;
1959  inputVecFields(3,0,0) =-1.0; inputVecFields(3,0,1) = 1.0; inputVecFields(3,0,2) =-1.0;
1960 
1961  // (C,P,D) - true
1962  outFieldsCorrect.resize(4,1,3);
1963  outFieldsCorrect(0,0,0) = 0.0; outFieldsCorrect(0,0,1) = 0.0; outFieldsCorrect(0,0,2) = 0.0;
1964  outFieldsCorrect(1,0,0) = 0.0; outFieldsCorrect(1,0,1) =-6.0; outFieldsCorrect(1,0,2) = 0.0;
1965  outFieldsCorrect(2,0,0) =-3.0; outFieldsCorrect(2,0,1) = 0.0; outFieldsCorrect(2,0,2) =-6.0;
1966  outFieldsCorrect(3,0,0) = 0.0; outFieldsCorrect(3,0,1) = 2.0; outFieldsCorrect(3,0,2) = 12.0;
1967 
1968  // (C,P,D)
1969  outFields.resize(4,1,3);
1970  art::matvecProductDataData<double>(outFields, inputMat, inputVecFields);
1971 
1972  // test loop
1973  for(int cell = 0; cell < outFields.dimension(0); cell++){
1974  for(int point = 0; point < outFields.dimension(1); point++){
1975  for(int row = 0; row < outFields.dimension(2); row++){
1976  if(outFields(cell, point, row) != outFieldsCorrect(cell, point, row)) {
1977  *outStream << "\n\nINCORRECT matvecProductDataData (1): \n value at multi-index ("
1978  << cell << "," << point << "," << row << ") = "
1979  << outFields(cell, point, row) << " but correct value is "
1980  << outFieldsCorrect(cell, point, row) <<"\n";
1981  errorFlag++;
1982  }
1983  }//row
1984  }// point
1985  }// cell
1986 
1987 
1988  *outStream \
1989  << "\n"
1990  << "===============================================================================\n"\
1991  << "| TEST 6.b: matvecProductDataData operations: (C,P,D,D) and (P,D) |\n"\
1992  << "===============================================================================\n";
1993  /*
1994  * inputMat inputVecFields outFields
1995  * 1 1 1 0 0 0 1 1 1 0 0 0 0 1 -1 -1 0 0 -3 0
1996  * -1 2 -1 -1 -2 -3 -1 2 -1 -1 -2 -3 0 1 -1 1 0 -6 0 2
1997  * 1 2 3 -2 6 -4 1 2 3 -2 6 -4 0 1 -1 -1 0 0 -6 12
1998  */
1999  // (C,P,D,D)
2000  inputMat.resize(1,4,3,3);
2001  // point 0
2002  inputMat(0,0,0,0) = 1.0; inputMat(0,0,0,1) = 1.0; inputMat(0,0,0,2) = 1.0;
2003  inputMat(0,0,1,0) =-1.0; inputMat(0,0,1,1) = 2.0; inputMat(0,0,1,2) =-1.0;
2004  inputMat(0,0,2,0) = 1.0; inputMat(0,0,2,1) = 2.0; inputMat(0,0,2,2) = 3.0;
2005  // point 1
2006  inputMat(0,1,0,0) = 0.0; inputMat(0,1,0,1) = 0.0; inputMat(0,1,0,2) = 0.0;
2007  inputMat(0,1,1,0) =-1.0; inputMat(0,1,1,1) =-2.0; inputMat(0,1,1,2) =-3.0;
2008  inputMat(0,1,2,0) =-2.0; inputMat(0,1,2,1) = 6.0; inputMat(0,1,2,2) =-4.0;
2009  // point 2
2010  inputMat(0,2,0,0) = 1.0; inputMat(0,2,0,1) = 1.0; inputMat(0,2,0,2) = 1.0;
2011  inputMat(0,2,1,0) =-1.0; inputMat(0,2,1,1) = 2.0; inputMat(0,2,1,2) =-1.0;
2012  inputMat(0,2,2,0) = 1.0; inputMat(0,2,2,1) = 2.0; inputMat(0,2,2,2) = 3.0;
2013  // point 3
2014  inputMat(0,3,0,0) = 0.0; inputMat(0,3,0,1) = 0.0; inputMat(0,3,0,2) = 0.0;
2015  inputMat(0,3,1,0) =-1.0; inputMat(0,3,1,1) =-2.0; inputMat(0,3,1,2) =-3.0;
2016  inputMat(0,3,2,0) =-2.0; inputMat(0,3,2,1) = 6.0; inputMat(0,3,2,2) =-4.0;
2017 
2018  // (P,D)
2019  inputVecFields.resize(4,3);
2020  //
2021  inputVecFields(0,0) = 0.0; inputVecFields(0,1) = 0.0; inputVecFields(0,2) = 0.0;
2022  inputVecFields(1,0) = 1.0; inputVecFields(1,1) = 1.0; inputVecFields(1,2) = 1.0;
2023  inputVecFields(2,0) =-1.0; inputVecFields(2,1) =-1.0; inputVecFields(2,2) =-1.0;
2024  inputVecFields(3,0) =-1.0; inputVecFields(3,1) = 1.0; inputVecFields(3,2) =-1.0;
2025 
2026  // (C,P,D) - true
2027  outFieldsCorrect.resize(1,4,3);
2028  outFieldsCorrect(0,0,0) = 0.0; outFieldsCorrect(0,0,1) = 0.0; outFieldsCorrect(0,0,2) = 0.0;
2029  outFieldsCorrect(0,1,0) = 0.0; outFieldsCorrect(0,1,1) =-6.0; outFieldsCorrect(0,1,2) = 0.0;
2030  outFieldsCorrect(0,2,0) =-3.0; outFieldsCorrect(0,2,1) = 0.0; outFieldsCorrect(0,2,2) =-6.0;
2031  outFieldsCorrect(0,3,0) = 0.0; outFieldsCorrect(0,3,1) = 2.0; outFieldsCorrect(0,3,2) = 12.0;
2032 
2033  // (C,P,D)
2034  outFields.resize(1,4,3);
2035  art::matvecProductDataData<double>(outFields, inputMat, inputVecFields);
2036 
2037  // test loop
2038  for(int cell = 0; cell < outFields.dimension(0); cell++){
2039  for(int point = 0; point < outFields.dimension(1); point++){
2040  for(int row = 0; row < outFields.dimension(2); row++){
2041  if(outFields(cell, point, row) != outFieldsCorrect(cell, point, row)) {
2042  *outStream << "\n\nINCORRECT matvecProductDataData (2): \n value at multi-index ("
2043  << cell << "," << point << "," << row << ") = "
2044  << outFields(cell, point, row) << " but correct value is "
2045  << outFieldsCorrect(cell, point, row) <<"\n";
2046  errorFlag++;
2047  }
2048  }//row
2049  }// point
2050  }// cell
2051 
2052 
2053 
2054  *outStream \
2055  << "\n"
2056  << "===============================================================================\n"\
2057  << "| TEST 6.c: matvecProductDataData random tests: branch inputDataRight(C,P,D) |\n"\
2058  << "===============================================================================\n";
2059  /*
2060  * Test derived from Test 5.c
2061  */
2062  {// test 6.c scope
2063  int c=5, p=9, d1=3;
2064  double zero = INTREPID_TOL*10000.0;
2065 
2066  FieldContainer<double> in_c_p_d(c, p, d1);
2067  FieldContainer<double> out_c_p_d(c, p, d1);
2068  FieldContainer<double> outi_c_p_d(c, p, d1);
2069 
2070  FieldContainer<double> data_c_p(c, p);
2071  FieldContainer<double> datainv_c_p(c, p);
2072  FieldContainer<double> data_c_1(c, 1);
2073  FieldContainer<double> datainv_c_1(c, 1);
2074  FieldContainer<double> data_c_p_d(c, p, d1);
2075  FieldContainer<double> datainv_c_p_d(c, p, d1);
2076  FieldContainer<double> data_c_1_d(c, 1, d1);
2077  FieldContainer<double> datainv_c_1_d(c, 1, d1);
2078  FieldContainer<double> data_c_p_d_d(c, p, d1, d1);
2079  FieldContainer<double> datainv_c_p_d_d(c, p, d1, d1);
2080  FieldContainer<double> datainvtrn_c_p_d_d(c, p, d1, d1);
2081  FieldContainer<double> data_c_1_d_d(c, 1, d1, d1);
2082  FieldContainer<double> datainv_c_1_d_d(c, 1, d1, d1);
2083  FieldContainer<double> datainvtrn_c_1_d_d(c, 1, d1, d1);
2084  /***********************************************************************************************
2085  * Constant diagonal tensor: inputDataLeft(C,P) *
2086  **********************************************************************************************/
2087  for (int i=0; i<in_c_p_d.size(); i++) {
2088  in_c_p_d[i] = Teuchos::ScalarTraits<double>::random();
2089  }
2090  for (int i=0; i<data_c_p.size(); i++) {
2091  data_c_p[i] = Teuchos::ScalarTraits<double>::random();
2092  datainv_c_p[i] = 1.0 / data_c_p[i];
2093  }
2094  for (int i=0; i<data_c_1.size(); i++) {
2095  data_c_1[i] = Teuchos::ScalarTraits<double>::random();
2096  datainv_c_1[i] = 1.0 / data_c_1[i];
2097  }
2098  // Tensor values vary by point:
2099  art::matvecProductDataData<double>(out_c_p_d, data_c_p, in_c_p_d);
2100  art::matvecProductDataData<double>(out_c_p_d, datainv_c_p, out_c_p_d);
2101  rst::subtract(&out_c_p_d[0], &in_c_p_d[0], out_c_p_d.size());
2102  if (rst::vectorNorm(&out_c_p_d[0], out_c_p_d.size(), NORM_ONE) > zero) {
2103  *outStream << "\n\nINCORRECT matvecProductDataData (3): check scalar inverse property\n\n";
2104  errorFlag = -1000;
2105  }
2106  // Tensor values do not vary by point:
2107  art::matvecProductDataData<double>(out_c_p_d, data_c_1, in_c_p_d);
2108  art::matvecProductDataData<double>(out_c_p_d, datainv_c_1, out_c_p_d);
2109  rst::subtract(&out_c_p_d[0], &in_c_p_d[0], out_c_p_d.size());
2110  if (rst::vectorNorm(&out_c_p_d[0], out_c_p_d.size(), NORM_ONE) > zero) {
2111  *outStream << "\n\nINCORRECT matvecProductDataData (4): check scalar inverse property\n\n";
2112  errorFlag = -1000;
2113  }
2114  /***********************************************************************************************
2115  * Non-onstant diagonal tensor: inputDataLeft(C,P,D) *
2116  **********************************************************************************************/
2117  for (int i=0; i<in_c_p_d.size(); i++) {
2118  in_c_p_d[i] = Teuchos::ScalarTraits<double>::random();
2119  }
2120  for (int i=0; i<data_c_p_d.size(); i++) {
2121  data_c_p_d[i] = Teuchos::ScalarTraits<double>::random();
2122  datainv_c_p_d[i] = 1.0 / data_c_p_d[i];
2123  }
2124  for (int i=0; i<data_c_1_d.size(); i++) {
2125  data_c_1_d[i] = Teuchos::ScalarTraits<double>::random();
2126  datainv_c_1_d[i] = 1.0 / data_c_1_d[i];
2127  }
2128  // Tensor values vary by point:
2129  art::matvecProductDataData<double>(out_c_p_d, data_c_p_d, in_c_p_d);
2130  art::matvecProductDataData<double>(out_c_p_d, datainv_c_p_d, out_c_p_d);
2131  rst::subtract(&out_c_p_d[0], &in_c_p_d[0], out_c_p_d.size());
2132  if (rst::vectorNorm(&out_c_p_d[0], out_c_p_d.size(), NORM_ONE) > zero) {
2133  *outStream << "\n\nINCORRECT matvecProductDataData (5): check scalar inverse property\n\n";
2134  errorFlag = -1000;
2135  }
2136  // Tensor values do not vary by point
2137  art::matvecProductDataData<double>(out_c_p_d, data_c_1_d, in_c_p_d);
2138  art::matvecProductDataData<double>(out_c_p_d, datainv_c_1_d, out_c_p_d);
2139  rst::subtract(&out_c_p_d[0], &in_c_p_d[0], out_c_p_d.size());
2140  if (rst::vectorNorm(&out_c_p_d[0], out_c_p_d.size(), NORM_ONE) > zero) {
2141  *outStream << "\n\nINCORRECT matvecProductDataData (6): check scalar inverse property\n\n";
2142  errorFlag = -1000;
2143  }
2144  /***********************************************************************************************
2145  * Full tensor: inputDataLeft(C,P,D,D) *
2146  **********************************************************************************************/
2147  for (int i=0; i<in_c_p_d.size(); i++) {
2148  in_c_p_d[i] = Teuchos::ScalarTraits<double>::random();
2149  }
2150  for (int ic=0; ic < c; ic++) {
2151  for (int ip=0; ip < p; ip++) {
2152  for (int i=0; i<d1*d1; i++) {
2153  (&data_c_p_d_d(ic, ip, 0, 0))[i] = Teuchos::ScalarTraits<double>::random();
2154  }
2155  RealSpaceTools<double>::inverse(&datainv_c_p_d_d(ic, ip, 0, 0), &data_c_p_d_d(ic, ip, 0, 0), d1);
2156  }
2157  }
2158  for (int ic=0; ic < c; ic++) {
2159  for (int ip=0; ip < 1; ip++) {
2160  for (int i=0; i<d1*d1; i++) {
2161  (&data_c_1_d_d(ic, 0, 0, 0))[i] = Teuchos::ScalarTraits<double>::random();
2162  }
2163  RealSpaceTools<double>::inverse(&datainv_c_1_d_d(ic, 0, 0, 0), &data_c_1_d_d(ic, 0, 0, 0), d1);
2164  }
2165  }
2166  // Tensor values vary by point: test "N" and "T" options (no transpose/transpose)
2167  art::matvecProductDataData<double>(out_c_p_d, data_c_p_d_d, in_c_p_d);
2168  art::matvecProductDataData<double>(outi_c_p_d, datainv_c_p_d_d, out_c_p_d);
2169  rst::subtract(&outi_c_p_d[0], &in_c_p_d[0], outi_c_p_d.size());
2170  if (rst::vectorNorm(&outi_c_p_d[0], outi_c_p_d.size(), NORM_ONE) > zero) {
2171  *outStream << "\n\nINCORRECT matvecProductDataData (7): check matrix inverse property\n\n";
2172  errorFlag = -1000;
2173  }
2174  art::matvecProductDataData<double>(out_c_p_d, data_c_p_d_d, in_c_p_d, 't');
2175  art::matvecProductDataData<double>(outi_c_p_d, datainv_c_p_d_d, out_c_p_d, 't');
2176  rst::subtract(&outi_c_p_d[0], &in_c_p_d[0], outi_c_p_d.size());
2177  if (rst::vectorNorm(&outi_c_p_d[0], outi_c_p_d.size(), NORM_ONE) > zero) {
2178  *outStream << "\n\nINCORRECT matvecProductDataData (8): check matrix inverse property, w/ double transpose\n\n";
2179  errorFlag = -1000;
2180  }
2181  // Tensor values do not vary by point: test "N" and "T" options (no transpose/transpose)
2182  art::matvecProductDataData<double>(out_c_p_d, data_c_1_d_d, in_c_p_d);
2183  art::matvecProductDataData<double>(outi_c_p_d, datainv_c_1_d_d, out_c_p_d);
2184  rst::subtract(&outi_c_p_d[0], &in_c_p_d[0], outi_c_p_d.size());
2185  if (rst::vectorNorm(&outi_c_p_d[0], outi_c_p_d.size(), NORM_ONE) > zero) {
2186  *outStream << "\n\nINCORRECT matvecProductDataData (9): check matrix inverse property\n\n";
2187  errorFlag = -1000;
2188  }
2189  art::matvecProductDataData<double>(out_c_p_d, data_c_1_d_d, in_c_p_d, 't');
2190  art::matvecProductDataData<double>(outi_c_p_d, datainv_c_1_d_d, out_c_p_d, 't');
2191  rst::subtract(&outi_c_p_d[0], &in_c_p_d[0], outi_c_p_d.size());
2192  if (rst::vectorNorm(&outi_c_p_d[0], outi_c_p_d.size(), NORM_ONE) > zero) {
2193  *outStream << "\n\nINCORRECT matvecProductDataData (10): check matrix inverse property, w/ double transpose\n\n";
2194  errorFlag = -1000;
2195  }
2196  /***********************************************************************************************
2197  * Full tensor: inputDataLeft(C,P,D,D) test inverse transpose *
2198  **********************************************************************************************/
2199  for (int i=0; i<in_c_p_d.size(); i++) {
2200  in_c_p_d[i] = Teuchos::ScalarTraits<double>::random();
2201  }
2202  for (int ic=0; ic < c; ic++) {
2203  for (int ip=0; ip < p; ip++) {
2204  for (int i=0; i<d1*d1; i++) {
2205  (&data_c_p_d_d(ic, ip, 0, 0))[i] = Teuchos::ScalarTraits<double>::random();
2206  }
2207  RealSpaceTools<double>::inverse(&datainv_c_p_d_d(ic, ip, 0, 0), &data_c_p_d_d(ic, ip, 0, 0), d1);
2208  RealSpaceTools<double>::transpose(&datainvtrn_c_p_d_d(ic, ip, 0, 0), &datainv_c_p_d_d(ic, ip, 0, 0), d1);
2209  }
2210  }
2211  for (int ic=0; ic < c; ic++) {
2212  for (int ip=0; ip < 1; ip++) {
2213  for (int i=0; i<d1*d1; i++) {
2214  (&data_c_1_d_d(ic, 0, 0, 0))[i] = Teuchos::ScalarTraits<double>::random();
2215  }
2216  RealSpaceTools<double>::inverse(&datainv_c_1_d_d(ic, 0, 0, 0), &data_c_1_d_d(ic, 0, 0, 0), d1);
2217  RealSpaceTools<double>::transpose(&datainvtrn_c_1_d_d(ic, 0, 0, 0), &datainv_c_1_d_d(ic, 0, 0, 0), d1);
2218  }
2219  }
2220  // Tensor values vary by point:
2221  art::matvecProductDataData<double>(out_c_p_d, data_c_p_d_d, in_c_p_d);
2222  art::matvecProductDataData<double>(outi_c_p_d, datainvtrn_c_p_d_d, out_c_p_d, 't');
2223  rst::subtract(&outi_c_p_d[0], &in_c_p_d[0], outi_c_p_d.size());
2224  if (rst::vectorNorm(&outi_c_p_d[0], outi_c_p_d.size(), NORM_ONE) > zero) {
2225  *outStream << "\n\nINCORRECT matvecProductDataData (11): check matrix inverse transpose property\n\n";
2226  errorFlag = -1000;
2227  }
2228  // Tensor values do not vary by point
2229  art::matvecProductDataData<double>(out_c_p_d, data_c_1_d_d, in_c_p_d);
2230  art::matvecProductDataData<double>(outi_c_p_d, datainvtrn_c_1_d_d, out_c_p_d, 't');
2231  rst::subtract(&outi_c_p_d[0], &in_c_p_d[0], outi_c_p_d.size());
2232  if (rst::vectorNorm(&outi_c_p_d[0], outi_c_p_d.size(), NORM_ONE) > zero) {
2233  *outStream << "\n\nINCORRECT matvecProductDataData (12): check matrix inverse transpose property\n\n";
2234  errorFlag = -1000;
2235  }
2236  }// test 6.c scope
2237 
2238 
2239  *outStream \
2240  << "\n"
2241  << "===============================================================================\n"\
2242  << "| TEST 6.d: matvecProductDataData random tests: branch inputDataRight(P,D) |\n"\
2243  << "===============================================================================\n";
2244  /*
2245  * Test derived from Test 5.d
2246  */
2247  {// test 6.d scope
2248  int c=5, p=9, d1=3;
2249  double zero = INTREPID_TOL*10000.0;
2250 
2251  FieldContainer<double> in_p_d(p, d1);
2252  FieldContainer<double> in_c_p_d(c, p, d1);
2253  FieldContainer<double> data_c_p(c, p);
2254  FieldContainer<double> datainv_c_p(c, p);
2255  FieldContainer<double> data_c_1(c, 1);
2256  FieldContainer<double> datainv_c_1(c, 1);
2257  FieldContainer<double> data_c_p_d(c, p, d1);
2258  FieldContainer<double> datainv_c_p_d(c, p, d1);
2259  FieldContainer<double> data_c_1_d(c, 1, d1);
2260  FieldContainer<double> datainv_c_1_d(c, 1, d1);
2261  FieldContainer<double> data_c_p_d_d(c, p, d1, d1);
2262  FieldContainer<double> datainv_c_p_d_d(c, p, d1, d1);
2263  FieldContainer<double> datainvtrn_c_p_d_d(c, p, d1, d1);
2264  FieldContainer<double> data_c_1_d_d(c, 1, d1, d1);
2265  FieldContainer<double> datainv_c_1_d_d(c, 1, d1, d1);
2266  FieldContainer<double> datainvtrn_c_1_d_d(c, 1, d1, d1);
2267  FieldContainer<double> data_c_p_one(c, p);
2268  FieldContainer<double> data_c_1_one(c, 1);
2269  FieldContainer<double> out_c_p_d(c, p, d1);
2270  FieldContainer<double> outi_c_p_d(c, p, d1);
2271  /***********************************************************************************************
2272  * Constant diagonal tensor: inputData(C,P) *
2273  **********************************************************************************************/
2274  for (int i=0; i<in_p_d.size(); i++) {
2275  in_p_d[i] = Teuchos::ScalarTraits<double>::random();
2276  }
2277  for (int i=0; i<data_c_p.size(); i++) {
2278  data_c_p[i] = Teuchos::ScalarTraits<double>::random();
2279  datainv_c_p[i] = 1.0 / data_c_p[i];
2280  data_c_p_one[i] = 1.0;
2281  }
2282  for (int i=0; i<data_c_1.size(); i++) {
2283  data_c_1[i] = Teuchos::ScalarTraits<double>::random();
2284  datainv_c_1[i] = 1.0 / data_c_1[i];
2285  }
2286  // Tensor values vary by point
2287  art::matvecProductDataData<double>(in_c_p_d, data_c_p_one, in_p_d);
2288  art::matvecProductDataData<double>(out_c_p_d, data_c_p, in_p_d);
2289  art::matvecProductDataData<double>(outi_c_p_d, datainv_c_p, out_c_p_d);
2290  rst::subtract(&outi_c_p_d[0], &in_c_p_d[0], outi_c_p_d.size());
2291  if (rst::vectorNorm(&outi_c_p_d[0], outi_c_p_d.size(), NORM_ONE) > zero) {
2292  *outStream << "\n\nINCORRECT matvecProductDataData (13): check scalar inverse property\n\n";
2293  errorFlag = -1000;
2294  }
2295  // Tensor values do not vary by point
2296  art::matvecProductDataData<double>(in_c_p_d, data_c_p_one, in_p_d);
2297  art::matvecProductDataData<double>(out_c_p_d, data_c_1, in_p_d);
2298  art::matvecProductDataData<double>(outi_c_p_d, datainv_c_1, out_c_p_d);
2299  rst::subtract(&outi_c_p_d[0], &in_c_p_d[0], outi_c_p_d.size());
2300  if (rst::vectorNorm(&outi_c_p_d[0], outi_c_p_d.size(), NORM_ONE) > zero) {
2301  *outStream << "\n\nINCORRECT matvecProductDataData (14): check scalar inverse property\n\n";
2302  errorFlag = -1000;
2303  }
2304  /***********************************************************************************************
2305  * Non-constant diagonal tensor: inputData(C,P,D) *
2306  **********************************************************************************************/
2307  for (int i=0; i<in_p_d.size(); i++) {
2308  in_p_d[i] = Teuchos::ScalarTraits<double>::random();
2309  }
2310  for (int i=0; i<data_c_p_d.size(); i++) {
2311  data_c_p_d[i] = Teuchos::ScalarTraits<double>::random();
2312  datainv_c_p_d[i] = 1.0 / data_c_p_d[i];
2313  }
2314  for (int i=0; i<data_c_1_d.size(); i++) {
2315  data_c_1_d[i] = Teuchos::ScalarTraits<double>::random();
2316  datainv_c_1_d[i] = 1.0 / data_c_1_d[i];
2317  }
2318  // Tensor values vary by point:
2319  art::matvecProductDataData<double>(in_c_p_d, data_c_p_one, in_p_d);
2320  art::matvecProductDataData<double>(out_c_p_d, data_c_p_d, in_p_d);
2321  art::matvecProductDataData<double>(outi_c_p_d, datainv_c_p_d, out_c_p_d);
2322  rst::subtract(&outi_c_p_d[0], &in_c_p_d[0], outi_c_p_d.size());
2323  if (rst::vectorNorm(&outi_c_p_d[0], outi_c_p_d.size(), NORM_ONE) > zero) {
2324  *outStream << "\n\nINCORRECT matvecProductDataData (15): check scalar inverse property\n\n";
2325  errorFlag = -1000;
2326  }
2327  // Tensor values do not vary by point:
2328  art::matvecProductDataData<double>(in_c_p_d, data_c_p_one, in_p_d);
2329  art::matvecProductDataData<double>(out_c_p_d, data_c_1_d, in_p_d);
2330  art::matvecProductDataData<double>(outi_c_p_d, datainv_c_1_d, out_c_p_d);
2331  rst::subtract(&outi_c_p_d[0], &in_c_p_d[0], outi_c_p_d.size());
2332  if (rst::vectorNorm(&outi_c_p_d[0], outi_c_p_d.size(), NORM_ONE) > zero) {
2333  *outStream << "\n\nINCORRECT matvecProductDataData (16): check scalar inverse property\n\n";
2334  errorFlag = -1000;
2335  }
2336  /***********************************************************************************************
2337  * Full tensor: inputData(C,P,D,D) *
2338  **********************************************************************************************/
2339  for (int i=0; i<in_p_d.size(); i++) {
2340  in_p_d[i] = Teuchos::ScalarTraits<double>::random();
2341  }
2342  for (int ic=0; ic < c; ic++) {
2343  for (int ip=0; ip < p; ip++) {
2344  for (int i=0; i<d1*d1; i++) {
2345  (&data_c_p_d_d(ic, ip, 0, 0))[i] = Teuchos::ScalarTraits<double>::random();
2346  }
2347  RealSpaceTools<double>::inverse(&datainv_c_p_d_d(ic, ip, 0, 0), &data_c_p_d_d(ic, ip, 0, 0), d1);
2348  }
2349  }
2350  for (int ic=0; ic < c; ic++) {
2351  for (int ip=0; ip < 1; ip++) {
2352  for (int i=0; i<d1*d1; i++) {
2353  (&data_c_1_d_d(ic, 0, 0, 0))[i] = Teuchos::ScalarTraits<double>::random();
2354  }
2355  RealSpaceTools<double>::inverse(&datainv_c_1_d_d(ic, 0, 0, 0), &data_c_1_d_d(ic, 0, 0, 0), d1);
2356  }
2357  }
2358  // Tensor values vary by point: test "N" and "T" (no-transpose/transpose) options
2359  art::matvecProductDataData<double>(in_c_p_d, data_c_p_one, in_p_d);
2360  art::matvecProductDataData<double>(out_c_p_d, data_c_p_d_d, in_p_d);
2361  art::matvecProductDataData<double>(outi_c_p_d, datainv_c_p_d_d, out_c_p_d);
2362  rst::subtract(&outi_c_p_d[0], &in_c_p_d[0], outi_c_p_d.size());
2363  if (rst::vectorNorm(&outi_c_p_d[0], outi_c_p_d.size(), NORM_ONE) > zero) {
2364  *outStream << "\n\nINCORRECT matvecProductDataData (17): check matrix inverse property\n\n";
2365  errorFlag = -1000;
2366  }
2367  art::matvecProductDataData<double>(in_c_p_d, data_c_p_one, in_p_d);
2368  art::matvecProductDataData<double>(out_c_p_d, data_c_p_d_d, in_p_d, 't');
2369  art::matvecProductDataData<double>(outi_c_p_d, datainv_c_p_d_d, out_c_p_d, 't');
2370  rst::subtract(&outi_c_p_d[0], &in_c_p_d[0], outi_c_p_d.size());
2371  if (rst::vectorNorm(&outi_c_p_d[0], outi_c_p_d.size(), NORM_ONE) > zero) {
2372  *outStream << "\n\nINCORRECT matvecProductDataData (18): check matrix inverse property, w/ double transpose\n\n";
2373  errorFlag = -1000;
2374  }
2375  // Tensor values do not vary by point: test "N" and "T" (no-transpose/transpose) options
2376  art::matvecProductDataData<double>(in_c_p_d, data_c_p_one, in_p_d);
2377  art::matvecProductDataData<double>(out_c_p_d, data_c_1_d_d, in_p_d);
2378  art::matvecProductDataData<double>(outi_c_p_d, datainv_c_1_d_d, out_c_p_d);
2379  rst::subtract(&outi_c_p_d[0], &in_c_p_d[0], outi_c_p_d.size());
2380  if (rst::vectorNorm(&outi_c_p_d[0], outi_c_p_d.size(), NORM_ONE) > zero) {
2381  *outStream << "\n\nINCORRECT matvecProductDataData (19): check matrix inverse property\n\n";
2382  errorFlag = -1000;
2383  }
2384  art::matvecProductDataData<double>(in_c_p_d, data_c_p_one, in_p_d);
2385  art::matvecProductDataData<double>(out_c_p_d, data_c_1_d_d, in_p_d, 't');
2386  art::matvecProductDataData<double>(outi_c_p_d, datainv_c_1_d_d, out_c_p_d, 't');
2387  rst::subtract(&outi_c_p_d[0], &in_c_p_d[0], outi_c_p_d.size());
2388  if (rst::vectorNorm(&outi_c_p_d[0], outi_c_p_d.size(), NORM_ONE) > zero) {
2389  *outStream << "\n\nINCORRECT matvecProductDataData (20): check matrix inverse property, w/ double transpose\n\n";
2390  errorFlag = -1000;
2391  }
2392  /***********************************************************************************************
2393  * Full tensor: inputData(C,P,D,D) test inverse transpose *
2394  **********************************************************************************************/
2395  for (int i=0; i<in_p_d.size(); i++) {
2396  in_p_d[i] = Teuchos::ScalarTraits<double>::random();
2397  }
2398  for (int ic=0; ic < c; ic++) {
2399  for (int ip=0; ip < p; ip++) {
2400  for (int i=0; i<d1*d1; i++) {
2401  (&data_c_p_d_d(ic, ip, 0, 0))[i] = Teuchos::ScalarTraits<double>::random();
2402  }
2403  RealSpaceTools<double>::inverse(&datainv_c_p_d_d(ic, ip, 0, 0), &data_c_p_d_d(ic, ip, 0, 0), d1);
2404  RealSpaceTools<double>::transpose(&datainvtrn_c_p_d_d(ic, ip, 0, 0), &datainv_c_p_d_d(ic, ip, 0, 0), d1);
2405  }
2406  }
2407  for (int ic=0; ic < c; ic++) {
2408  for (int ip=0; ip < 1; ip++) {
2409  for (int i=0; i<d1*d1; i++) {
2410  (&data_c_1_d_d(ic, 0, 0, 0))[i] = Teuchos::ScalarTraits<double>::random();
2411  }
2412  RealSpaceTools<double>::inverse(&datainv_c_1_d_d(ic, 0, 0, 0), &data_c_1_d_d(ic, 0, 0, 0), d1);
2413  RealSpaceTools<double>::transpose(&datainvtrn_c_1_d_d(ic, 0, 0, 0), &datainv_c_1_d_d(ic, 0, 0, 0), d1);
2414  }
2415  }
2416  // Tensor values vary by point:
2417  art::matvecProductDataData<double>(in_c_p_d, data_c_p_one, in_p_d);
2418  art::matvecProductDataData<double>(out_c_p_d, data_c_p_d_d, in_p_d);
2419  art::matvecProductDataData<double>(outi_c_p_d, datainvtrn_c_p_d_d, out_c_p_d, 't');
2420  rst::subtract(&outi_c_p_d[0], &in_c_p_d[0], outi_c_p_d.size());
2421  if (rst::vectorNorm(&outi_c_p_d[0], outi_c_p_d.size(), NORM_ONE) > zero) {
2422  *outStream << "\n\nINCORRECT matvecProductDataData (21): check matrix inverse transpose property\n\n";
2423  errorFlag = -1000;
2424  }
2425  // Tensor values do not vary by point:
2426  art::matvecProductDataData<double>(in_c_p_d, data_c_p_one, in_p_d);
2427  art::matvecProductDataData<double>(out_c_p_d, data_c_1_d_d, in_p_d);
2428  art::matvecProductDataData<double>(outi_c_p_d, datainvtrn_c_1_d_d, out_c_p_d, 't');
2429  rst::subtract(&outi_c_p_d[0], &in_c_p_d[0], outi_c_p_d.size());
2430  if (rst::vectorNorm(&outi_c_p_d[0], outi_c_p_d.size(), NORM_ONE) > zero) {
2431  *outStream << "\n\nINCORRECT matvecProductDataData (22): check matrix inverse transpose property\n\n";
2432  errorFlag = -1000;
2433  }
2434  }// test 6.d scope
2435 
2436 
2437 
2438  *outStream \
2439  << "\n"
2440  << "===============================================================================\n"\
2441  << "| TEST 7.a: matmatProductDataField random tests: branch inputFields(C,F,P,D,D)|\n"\
2442  << "===============================================================================\n";
2443  {// Test 7.a scope
2444  int c=5, p=9, f=7, d1=3;
2445  double zero = INTREPID_TOL*10000.0;
2446 
2447  FieldContainer<double> in_c_f_p_d_d(c, f, p, d1, d1);
2448  FieldContainer<double> data_c_p(c, p);
2449  FieldContainer<double> datainv_c_p(c, p);
2450  FieldContainer<double> data_c_1(c, 1);
2451  FieldContainer<double> datainv_c_1(c, 1);
2452  FieldContainer<double> data_c_p_d(c, p, d1);
2453  FieldContainer<double> datainv_c_p_d(c, p, d1);
2454  FieldContainer<double> data_c_1_d(c, 1, d1);
2455  FieldContainer<double> datainv_c_1_d(c, 1, d1);
2456  FieldContainer<double> data_c_p_d_d(c, p, d1, d1);
2457  FieldContainer<double> datainv_c_p_d_d(c, p, d1, d1);
2458  FieldContainer<double> datainvtrn_c_p_d_d(c, p, d1, d1);
2459  FieldContainer<double> data_c_1_d_d(c, 1, d1, d1);
2460  FieldContainer<double> datainv_c_1_d_d(c, 1, d1, d1);
2461  FieldContainer<double> datainvtrn_c_1_d_d(c, 1, d1, d1);
2462  FieldContainer<double> out_c_f_p_d_d(c, f, p, d1, d1);
2463  FieldContainer<double> outi_c_f_p_d_d(c, f, p, d1, d1);
2464  /***********************************************************************************************
2465  * Constant diagonal tensor: inputData(C,P) *
2466  **********************************************************************************************/
2467  for (int i=0; i<in_c_f_p_d_d.size(); i++) {
2468  in_c_f_p_d_d[i] = Teuchos::ScalarTraits<double>::random();
2469  }
2470  for (int i=0; i<data_c_p.size(); i++) {
2471  data_c_p[i] = Teuchos::ScalarTraits<double>::random();
2472  datainv_c_p[i] = 1.0 / data_c_p[i];
2473  }
2474  for (int i=0; i<data_c_1.size(); i++) {
2475  data_c_1[i] = Teuchos::ScalarTraits<double>::random();
2476  datainv_c_1[i] = 1.0 / data_c_1[i];
2477  }
2478  // Tensor values vary by point:
2479  art::matmatProductDataField<double>(out_c_f_p_d_d, data_c_p, in_c_f_p_d_d);
2480  art::matmatProductDataField<double>(out_c_f_p_d_d, datainv_c_p, out_c_f_p_d_d);
2481  rst::subtract(&out_c_f_p_d_d[0], &in_c_f_p_d_d[0], out_c_f_p_d_d.size());
2482  if (rst::vectorNorm(&out_c_f_p_d_d[0], out_c_f_p_d_d.size(), NORM_ONE) > zero) {
2483  *outStream << "\n\nINCORRECT matmatProductDataField (1): check scalar inverse property\n\n";
2484  errorFlag = -1000;
2485  }
2486  // Tensor values do not vary by point:
2487  art::matmatProductDataField<double>(out_c_f_p_d_d, data_c_1, in_c_f_p_d_d);
2488  art::matmatProductDataField<double>(out_c_f_p_d_d, datainv_c_1, out_c_f_p_d_d);
2489  rst::subtract(&out_c_f_p_d_d[0], &in_c_f_p_d_d[0], out_c_f_p_d_d.size());
2490  if (rst::vectorNorm(&out_c_f_p_d_d[0], out_c_f_p_d_d.size(), NORM_ONE) > zero) {
2491  *outStream << "\n\nINCORRECT matmatProductDataField (2): check scalar inverse property\n\n";
2492  errorFlag = -1000;
2493  }
2494  /***********************************************************************************************
2495  * Non-onstant diagonal tensor: inputData(C,P,D) *
2496  **********************************************************************************************/
2497  for (int i=0; i<in_c_f_p_d_d.size(); i++) {
2498  in_c_f_p_d_d[i] = Teuchos::ScalarTraits<double>::random();
2499  }
2500  for (int i=0; i<data_c_p_d.size(); i++) {
2501  data_c_p_d[i] = Teuchos::ScalarTraits<double>::random();
2502  datainv_c_p_d[i] = 1.0 / data_c_p_d[i];
2503  }
2504  for (int i=0; i<data_c_1_d.size(); i++) {
2505  data_c_1_d[i] = Teuchos::ScalarTraits<double>::random();
2506  datainv_c_1_d[i] = 1.0 / data_c_1_d[i];
2507  }
2508  // Tensor values vary by point:
2509  art::matmatProductDataField<double>(out_c_f_p_d_d, data_c_p_d, in_c_f_p_d_d);
2510  art::matmatProductDataField<double>(out_c_f_p_d_d, datainv_c_p_d, out_c_f_p_d_d);
2511  rst::subtract(&out_c_f_p_d_d[0], &in_c_f_p_d_d[0], out_c_f_p_d_d.size());
2512  if (rst::vectorNorm(&out_c_f_p_d_d[0], out_c_f_p_d_d.size(), NORM_ONE) > zero) {
2513  *outStream << "\n\nINCORRECT matmatProductDataField (3): check scalar inverse property\n\n";
2514  errorFlag = -1000;
2515  }
2516  // Tensor values do not vary by point
2517  art::matmatProductDataField<double>(out_c_f_p_d_d, data_c_1_d, in_c_f_p_d_d);
2518  art::matmatProductDataField<double>(out_c_f_p_d_d, datainv_c_1_d, out_c_f_p_d_d);
2519  rst::subtract(&out_c_f_p_d_d[0], &in_c_f_p_d_d[0], out_c_f_p_d_d.size());
2520  if (rst::vectorNorm(&out_c_f_p_d_d[0], out_c_f_p_d_d.size(), NORM_ONE) > zero) {
2521  *outStream << "\n\nINCORRECT matmatProductDataField (4): check scalar inverse property\n\n";
2522  errorFlag = -1000;
2523  }
2524  /***********************************************************************************************
2525  * Full tensor: inputData(C,P,D,D) *
2526  **********************************************************************************************/
2527  for (int i=0; i<in_c_f_p_d_d.size(); i++) {
2528  in_c_f_p_d_d[i] = Teuchos::ScalarTraits<double>::random();
2529  }
2530  for (int ic=0; ic < c; ic++) {
2531  for (int ip=0; ip < p; ip++) {
2532  for (int i=0; i<d1*d1; i++) {
2533  (&data_c_p_d_d(ic, ip, 0, 0))[i] = Teuchos::ScalarTraits<double>::random();
2534  }
2535  RealSpaceTools<double>::inverse(&datainv_c_p_d_d(ic, ip, 0, 0), &data_c_p_d_d(ic, ip, 0, 0), d1);
2536  }
2537  }
2538  for (int ic=0; ic < c; ic++) {
2539  for (int ip=0; ip < 1; ip++) {
2540  for (int i=0; i<d1*d1; i++) {
2541  (&data_c_1_d_d(ic, 0, 0, 0))[i] = Teuchos::ScalarTraits<double>::random();
2542  }
2543  RealSpaceTools<double>::inverse(&datainv_c_1_d_d(ic, 0, 0, 0), &data_c_1_d_d(ic, 0, 0, 0), d1);
2544  }
2545  }
2546  // Tensor values vary by point: test "N" and "T" options (no transpose/transpose)
2547  art::matmatProductDataField<double>(out_c_f_p_d_d, data_c_p_d_d, in_c_f_p_d_d);
2548  art::matmatProductDataField<double>(outi_c_f_p_d_d, datainv_c_p_d_d, out_c_f_p_d_d);
2549  rst::subtract(&outi_c_f_p_d_d[0], &in_c_f_p_d_d[0], outi_c_f_p_d_d.size());
2550  if (rst::vectorNorm(&outi_c_f_p_d_d[0], outi_c_f_p_d_d.size(), NORM_ONE) > zero) {
2551  *outStream << "\n\nINCORRECT matmatProductDataField (5): check matrix inverse property\n\n";
2552  errorFlag = -1000;
2553  }
2554  art::matmatProductDataField<double>(out_c_f_p_d_d, data_c_p_d_d, in_c_f_p_d_d, 't');
2555  art::matmatProductDataField<double>(outi_c_f_p_d_d, datainv_c_p_d_d, out_c_f_p_d_d, 't');
2556  rst::subtract(&outi_c_f_p_d_d[0], &in_c_f_p_d_d[0], outi_c_f_p_d_d.size());
2557  if (rst::vectorNorm(&outi_c_f_p_d_d[0], outi_c_f_p_d_d.size(), NORM_ONE) > zero) {
2558  *outStream << "\n\nINCORRECT matmatProductDataField (6): check matrix inverse property, w/ double transpose\n\n";
2559  errorFlag = -1000;
2560  }
2561  // Tensor values do not vary by point: test "N" and "T" options (no transpose/transpose)
2562  art::matmatProductDataField<double>(out_c_f_p_d_d, data_c_1_d_d, in_c_f_p_d_d);
2563  art::matmatProductDataField<double>(outi_c_f_p_d_d, datainv_c_1_d_d, out_c_f_p_d_d);
2564  rst::subtract(&outi_c_f_p_d_d[0], &in_c_f_p_d_d[0], outi_c_f_p_d_d.size());
2565  if (rst::vectorNorm(&outi_c_f_p_d_d[0], outi_c_f_p_d_d.size(), NORM_ONE) > zero) {
2566  *outStream << "\n\nINCORRECT matmatProductDataField (7): check matrix inverse property\n\n";
2567  errorFlag = -1000;
2568  }
2569  art::matmatProductDataField<double>(out_c_f_p_d_d, data_c_1_d_d, in_c_f_p_d_d,'t');
2570  art::matmatProductDataField<double>(outi_c_f_p_d_d, datainv_c_1_d_d, out_c_f_p_d_d, 't');
2571  rst::subtract(&outi_c_f_p_d_d[0], &in_c_f_p_d_d[0], outi_c_f_p_d_d.size());
2572  if (rst::vectorNorm(&outi_c_f_p_d_d[0], outi_c_f_p_d_d.size(), NORM_ONE) > zero) {
2573  *outStream << "\n\nINCORRECT matmatProductDataField (8): check matrix inverse property, w/ double transpose\n\n";
2574  errorFlag = -1000;
2575  }
2576  /***********************************************************************************************
2577  * Full tensor: inputData(C,P,D,D) test inverse transpose *
2578  **********************************************************************************************/
2579  for (int i=0; i<in_c_f_p_d_d.size(); i++) {
2580  in_c_f_p_d_d[i] = Teuchos::ScalarTraits<double>::random();
2581  }
2582  for (int ic=0; ic < c; ic++) {
2583  for (int ip=0; ip < p; ip++) {
2584  for (int i=0; i<d1*d1; i++) {
2585  (&data_c_p_d_d(ic, ip, 0, 0))[i] = Teuchos::ScalarTraits<double>::random();
2586  }
2587  RealSpaceTools<double>::inverse(&datainv_c_p_d_d(ic, ip, 0, 0), &data_c_p_d_d(ic, ip, 0, 0), d1);
2588  RealSpaceTools<double>::transpose(&datainvtrn_c_p_d_d(ic, ip, 0, 0), &datainv_c_p_d_d(ic, ip, 0, 0), d1);
2589  }
2590  }
2591  for (int ic=0; ic < c; ic++) {
2592  for (int ip=0; ip < 1; ip++) {
2593  for (int i=0; i<d1*d1; i++) {
2594  (&data_c_1_d_d(ic, 0, 0, 0))[i] = Teuchos::ScalarTraits<double>::random();
2595  }
2596  RealSpaceTools<double>::inverse(&datainv_c_1_d_d(ic, 0, 0, 0), &data_c_1_d_d(ic, 0, 0, 0), d1);
2597  RealSpaceTools<double>::transpose(&datainvtrn_c_1_d_d(ic, 0, 0, 0), &datainv_c_1_d_d(ic, 0, 0, 0), d1);
2598  }
2599  }
2600  // Tensor values vary by point:
2601  art::matmatProductDataField<double>(out_c_f_p_d_d, data_c_p_d_d, in_c_f_p_d_d);
2602  art::matmatProductDataField<double>(outi_c_f_p_d_d, datainvtrn_c_p_d_d, out_c_f_p_d_d, 't');
2603  rst::subtract(&outi_c_f_p_d_d[0], &in_c_f_p_d_d[0], outi_c_f_p_d_d.size());
2604  if (rst::vectorNorm(&outi_c_f_p_d_d[0], outi_c_f_p_d_d.size(), NORM_ONE) > zero) {
2605  *outStream << "\n\nINCORRECT matmatProductDataField (9): check matrix inverse transpose property\n\n";
2606  errorFlag = -1000;
2607  }
2608  // Tensor values do not vary by point
2609  art::matmatProductDataField<double>(out_c_f_p_d_d, data_c_1_d_d, in_c_f_p_d_d);
2610  art::matmatProductDataField<double>(outi_c_f_p_d_d, datainvtrn_c_1_d_d, out_c_f_p_d_d, 't');
2611  rst::subtract(&outi_c_f_p_d_d[0], &in_c_f_p_d_d[0], outi_c_f_p_d_d.size());
2612  if (rst::vectorNorm(&outi_c_f_p_d_d[0], outi_c_f_p_d_d.size(), NORM_ONE) > zero) {
2613  *outStream << "\n\nINCORRECT matmatProductDataField (10): check matrix inverse transpose property\n\n";
2614  errorFlag = -1000;
2615  }
2616  }// test 7.a scope
2617 
2618 
2619 
2620  *outStream \
2621  << "\n"
2622  << "===============================================================================\n"\
2623  << "| TEST 7.b: matmatProductDataField random tests: branch inputFields(F,P,D,D) |\n"\
2624  << "===============================================================================\n";
2625  {// Test 7.b scope
2626  int c=5, p=9, f=7, d1=3;
2627  double zero = INTREPID_TOL*10000.0;
2628 
2629  FieldContainer<double> in_f_p_d_d(f, p, d1, d1);
2630  FieldContainer<double> in_c_f_p_d_d(c, f, p, d1, d1);
2631  FieldContainer<double> data_c_p(c, p);
2632  FieldContainer<double> datainv_c_p(c, p);
2633  FieldContainer<double> data_c_1(c, 1);
2634  FieldContainer<double> datainv_c_1(c, 1);
2635  FieldContainer<double> data_c_p_d(c, p, d1);
2636  FieldContainer<double> datainv_c_p_d(c, p, d1);
2637  FieldContainer<double> data_c_1_d(c, 1, d1);
2638  FieldContainer<double> datainv_c_1_d(c, 1, d1);
2639  FieldContainer<double> data_c_p_d_d(c, p, d1, d1);
2640  FieldContainer<double> datainv_c_p_d_d(c, p, d1, d1);
2641  FieldContainer<double> datainvtrn_c_p_d_d(c, p, d1, d1);
2642  FieldContainer<double> data_c_1_d_d(c, 1, d1, d1);
2643  FieldContainer<double> datainv_c_1_d_d(c, 1, d1, d1);
2644  FieldContainer<double> datainvtrn_c_1_d_d(c, 1, d1, d1);
2645  FieldContainer<double> data_c_p_one(c, p);
2646  FieldContainer<double> data_c_1_one(c, 1);
2647  FieldContainer<double> out_c_f_p_d_d(c, f, p, d1, d1);
2648  FieldContainer<double> outi_c_f_p_d_d(c, f, p, d1, d1);
2649  /***********************************************************************************************
2650  * Constant diagonal tensor: inputData(C,P) *
2651  **********************************************************************************************/
2652  for (int i=0; i<in_f_p_d_d.size(); i++) {
2653  in_f_p_d_d[i] = Teuchos::ScalarTraits<double>::random();
2654  }
2655  for (int i=0; i<data_c_p.size(); i++) {
2656  data_c_p[i] = Teuchos::ScalarTraits<double>::random();
2657  datainv_c_p[i] = 1.0 / data_c_p[i];
2658  data_c_p_one[i] = 1.0;
2659  }
2660  for (int i=0; i<data_c_1.size(); i++) {
2661  data_c_1[i] = Teuchos::ScalarTraits<double>::random();
2662  datainv_c_1[i] = 1.0 / data_c_1[i];
2663  }
2664 
2665  // Tensor values vary by point
2666  art::matmatProductDataField<double>(in_c_f_p_d_d, data_c_p_one, in_f_p_d_d);
2667  art::matmatProductDataField<double>(out_c_f_p_d_d, data_c_p, in_f_p_d_d);
2668  art::matmatProductDataField<double>(outi_c_f_p_d_d, datainv_c_p, out_c_f_p_d_d);
2669  rst::subtract(&outi_c_f_p_d_d[0], &in_c_f_p_d_d[0], outi_c_f_p_d_d.size());
2670  if (rst::vectorNorm(&outi_c_f_p_d_d[0], outi_c_f_p_d_d.size(), NORM_ONE) > zero) {
2671  *outStream << "\n\nINCORRECT matmatProductDataField (11): check scalar inverse property\n\n";
2672  errorFlag = -1000;
2673  }
2674  // Tensor values do not vary by point
2675  art::matmatProductDataField<double>(in_c_f_p_d_d, data_c_p_one, in_f_p_d_d);
2676  art::matmatProductDataField<double>(out_c_f_p_d_d, data_c_1, in_f_p_d_d);
2677  art::matmatProductDataField<double>(outi_c_f_p_d_d, datainv_c_1, out_c_f_p_d_d);
2678  rst::subtract(&outi_c_f_p_d_d[0], &in_c_f_p_d_d[0], outi_c_f_p_d_d.size());
2679  if (rst::vectorNorm(&outi_c_f_p_d_d[0], outi_c_f_p_d_d.size(), NORM_ONE) > zero) {
2680  *outStream << "\n\nINCORRECT matmatProductDataField (12): check scalar inverse property\n\n";
2681  errorFlag = -1000;
2682  }
2683  /***********************************************************************************************
2684  * Non-constant diagonal tensor: inputData(C,P,D) *
2685  **********************************************************************************************/
2686  for (int i=0; i<in_f_p_d_d.size(); i++) {
2687  in_f_p_d_d[i] = Teuchos::ScalarTraits<double>::random();
2688  }
2689  for (int i=0; i<data_c_p_d.size(); i++) {
2690  data_c_p_d[i] = Teuchos::ScalarTraits<double>::random();
2691  datainv_c_p_d[i] = 1.0 / data_c_p_d[i];
2692  }
2693  for (int i=0; i<data_c_1_d.size(); i++) {
2694  data_c_1_d[i] = Teuchos::ScalarTraits<double>::random();
2695  datainv_c_1_d[i] = 1.0 / data_c_1_d[i];
2696  }
2697  // Tensor values vary by point:
2698  art::matmatProductDataField<double>(in_c_f_p_d_d, data_c_p_one, in_f_p_d_d);
2699  art::matmatProductDataField<double>(out_c_f_p_d_d, data_c_p_d, in_f_p_d_d);
2700  art::matmatProductDataField<double>(outi_c_f_p_d_d, datainv_c_p_d, out_c_f_p_d_d);
2701  rst::subtract(&outi_c_f_p_d_d[0], &in_c_f_p_d_d[0], outi_c_f_p_d_d.size());
2702  if (rst::vectorNorm(&outi_c_f_p_d_d[0], outi_c_f_p_d_d.size(), NORM_ONE) > zero) {
2703  *outStream << "\n\nINCORRECT matmatProductDataField (13): check scalar inverse property\n\n";
2704  errorFlag = -1000;
2705  }
2706  // Tensor values do not vary by point:
2707  art::matmatProductDataField<double>(in_c_f_p_d_d, data_c_p_one, in_f_p_d_d);
2708  art::matmatProductDataField<double>(out_c_f_p_d_d, data_c_1_d, in_f_p_d_d);
2709  art::matmatProductDataField<double>(outi_c_f_p_d_d, datainv_c_1_d, out_c_f_p_d_d);
2710  rst::subtract(&outi_c_f_p_d_d[0], &in_c_f_p_d_d[0], outi_c_f_p_d_d.size());
2711  if (rst::vectorNorm(&outi_c_f_p_d_d[0], outi_c_f_p_d_d.size(), NORM_ONE) > zero) {
2712  *outStream << "\n\nINCORRECT matmatProductDataField (14): check scalar inverse property\n\n";
2713  errorFlag = -1000;
2714  }
2715  /***********************************************************************************************
2716  * Full tensor: inputData(C,P,D,D) *
2717  **********************************************************************************************/
2718  for (int i=0; i<in_f_p_d_d.size(); i++) {
2719  in_f_p_d_d[i] = Teuchos::ScalarTraits<double>::random();
2720  }
2721  for (int ic=0; ic < c; ic++) {
2722  for (int ip=0; ip < p; ip++) {
2723  for (int i=0; i<d1*d1; i++) {
2724  (&data_c_p_d_d(ic, ip, 0, 0))[i] = Teuchos::ScalarTraits<double>::random();
2725  }
2726  RealSpaceTools<double>::inverse(&datainv_c_p_d_d(ic, ip, 0, 0), &data_c_p_d_d(ic, ip, 0, 0), d1);
2727  }
2728  }
2729  for (int ic=0; ic < c; ic++) {
2730  for (int ip=0; ip < 1; ip++) {
2731  for (int i=0; i<d1*d1; i++) {
2732  (&data_c_1_d_d(ic, 0, 0, 0))[i] = Teuchos::ScalarTraits<double>::random();
2733  }
2734  RealSpaceTools<double>::inverse(&datainv_c_1_d_d(ic, 0, 0, 0), &data_c_1_d_d(ic, 0, 0, 0), d1);
2735  }
2736  }
2737  // Tensor values vary by point: test "N" and "T" (no-transpose/transpose) options
2738  art::matmatProductDataField<double>(in_c_f_p_d_d, data_c_p_one, in_f_p_d_d);
2739  art::matmatProductDataField<double>(out_c_f_p_d_d, data_c_p_d_d, in_f_p_d_d);
2740  art::matmatProductDataField<double>(outi_c_f_p_d_d, datainv_c_p_d_d, out_c_f_p_d_d);
2741  rst::subtract(&outi_c_f_p_d_d[0], &in_c_f_p_d_d[0], outi_c_f_p_d_d.size());
2742  if (rst::vectorNorm(&outi_c_f_p_d_d[0], outi_c_f_p_d_d.size(), NORM_ONE) > zero) {
2743  *outStream << "\n\nINCORRECT matmatProductDataField (15): check matrix inverse property\n\n";
2744  errorFlag = -1000;
2745  }
2746  art::matmatProductDataField<double>(in_c_f_p_d_d, data_c_p_one, in_f_p_d_d);
2747  art::matmatProductDataField<double>(out_c_f_p_d_d, data_c_p_d_d, in_f_p_d_d, 't');
2748  art::matmatProductDataField<double>(outi_c_f_p_d_d, datainv_c_p_d_d, out_c_f_p_d_d, 't');
2749  rst::subtract(&outi_c_f_p_d_d[0], &in_c_f_p_d_d[0], outi_c_f_p_d_d.size());
2750  if (rst::vectorNorm(&outi_c_f_p_d_d[0], outi_c_f_p_d_d.size(), NORM_ONE) > zero) {
2751  *outStream << "\n\nINCORRECT matmatProductDataField (16): check matrix inverse property, w/ double transpose\n\n";
2752  errorFlag = -1000;
2753  }
2754  // Tensor values do not vary by point: test "N" and "T" (no-transpose/transpose) options
2755  art::matmatProductDataField<double>(in_c_f_p_d_d, data_c_p_one, in_f_p_d_d);
2756  art::matmatProductDataField<double>(out_c_f_p_d_d, data_c_1_d_d, in_f_p_d_d);
2757  art::matmatProductDataField<double>(outi_c_f_p_d_d, datainv_c_1_d_d, out_c_f_p_d_d);
2758  rst::subtract(&outi_c_f_p_d_d[0], &in_c_f_p_d_d[0], outi_c_f_p_d_d.size());
2759  if (rst::vectorNorm(&outi_c_f_p_d_d[0], outi_c_f_p_d_d.size(), NORM_ONE) > zero) {
2760  *outStream << "\n\nINCORRECT matmatProductDataField (17): check matrix inverse property\n\n";
2761  errorFlag = -1000;
2762  }
2763  art::matmatProductDataField<double>(in_c_f_p_d_d, data_c_p_one, in_f_p_d_d);
2764  art::matmatProductDataField<double>(out_c_f_p_d_d, data_c_1_d_d, in_f_p_d_d, 't');
2765  art::matmatProductDataField<double>(outi_c_f_p_d_d, datainv_c_1_d_d, out_c_f_p_d_d, 't');
2766  rst::subtract(&outi_c_f_p_d_d[0], &in_c_f_p_d_d[0], outi_c_f_p_d_d.size());
2767  if (rst::vectorNorm(&outi_c_f_p_d_d[0], outi_c_f_p_d_d.size(), NORM_ONE) > zero) {
2768  *outStream << "\n\nINCORRECT matmatProductDataField (18): check matrix inverse property, w/ double transpose\n\n";
2769  errorFlag = -1000;
2770  }
2771  /***********************************************************************************************
2772  * Full tensor: inputData(C,P,D,D) test inverse transpose *
2773  **********************************************************************************************/
2774  for (int i=0; i<in_f_p_d_d.size(); i++) {
2775  in_f_p_d_d[i] = Teuchos::ScalarTraits<double>::random();
2776  }
2777  for (int ic=0; ic < c; ic++) {
2778  for (int ip=0; ip < p; ip++) {
2779  for (int i=0; i<d1*d1; i++) {
2780  (&data_c_p_d_d(ic, ip, 0, 0))[i] = Teuchos::ScalarTraits<double>::random();
2781  }
2782  RealSpaceTools<double>::inverse(&datainv_c_p_d_d(ic, ip, 0, 0), &data_c_p_d_d(ic, ip, 0, 0), d1);
2783  RealSpaceTools<double>::transpose(&datainvtrn_c_p_d_d(ic, ip, 0, 0), &datainv_c_p_d_d(ic, ip, 0, 0), d1);
2784  }
2785  }
2786  for (int ic=0; ic < c; ic++) {
2787  for (int ip=0; ip < 1; ip++) {
2788  for (int i=0; i<d1*d1; i++) {
2789  (&data_c_1_d_d(ic, 0, 0, 0))[i] = Teuchos::ScalarTraits<double>::random();
2790  }
2791  RealSpaceTools<double>::inverse(&datainv_c_1_d_d(ic, 0, 0, 0), &data_c_1_d_d(ic, 0, 0, 0), d1);
2792  RealSpaceTools<double>::transpose(&datainvtrn_c_1_d_d(ic, 0, 0, 0), &datainv_c_1_d_d(ic, 0, 0, 0), d1);
2793  }
2794  }
2795  // Tensor values vary by point:
2796  art::matmatProductDataField<double>(in_c_f_p_d_d, data_c_p_one, in_f_p_d_d);
2797  art::matmatProductDataField<double>(out_c_f_p_d_d, data_c_p_d_d, in_f_p_d_d);
2798  art::matmatProductDataField<double>(outi_c_f_p_d_d, datainvtrn_c_p_d_d, out_c_f_p_d_d, 't');
2799  rst::subtract(&outi_c_f_p_d_d[0], &in_c_f_p_d_d[0], outi_c_f_p_d_d.size());
2800  if (rst::vectorNorm(&outi_c_f_p_d_d[0], outi_c_f_p_d_d.size(), NORM_ONE) > zero) {
2801  *outStream << "\n\nINCORRECT matmatProductDataField (19): check matrix inverse transpose property\n\n";
2802  errorFlag = -1000;
2803  }
2804  // Tensor values do not vary by point:
2805  art::matmatProductDataField<double>(in_c_f_p_d_d, data_c_p_one, in_f_p_d_d);
2806  art::matmatProductDataField<double>(out_c_f_p_d_d, data_c_1_d_d, in_f_p_d_d);
2807  art::matmatProductDataField<double>(outi_c_f_p_d_d, datainvtrn_c_1_d_d, out_c_f_p_d_d, 't');
2808  rst::subtract(&outi_c_f_p_d_d[0], &in_c_f_p_d_d[0], outi_c_f_p_d_d.size());
2809  if (rst::vectorNorm(&outi_c_f_p_d_d[0], outi_c_f_p_d_d.size(), NORM_ONE) > zero) {
2810  *outStream << "\n\nINCORRECT matmatProductDataField (20): check matrix inverse transpose property\n\n";
2811  errorFlag = -1000;
2812  }
2813  }// test 7.b scope
2814 
2815 
2816 
2817  *outStream \
2818  << "\n"
2819  << "===============================================================================\n"\
2820  << "| TEST 8.a: matmatProductDataData random tests: branch inputDataRight(C,P,D,D)|\n"\
2821  << "===============================================================================\n";
2822  /*
2823  * Test derived from Test 7.a
2824  */
2825  {// test 8.a scope
2826  int c=5, p=9, d1=3;
2827  double zero = INTREPID_TOL*10000.0;
2828 
2829  FieldContainer<double> in_c_p_d_d(c, p, d1, d1);
2830  FieldContainer<double> out_c_p_d_d(c, p, d1, d1);
2831  FieldContainer<double> outi_c_p_d_d(c, p, d1, d1);
2832 
2833  FieldContainer<double> data_c_p(c, p);
2834  FieldContainer<double> datainv_c_p(c, p);
2835  FieldContainer<double> data_c_1(c, 1);
2836  FieldContainer<double> datainv_c_1(c, 1);
2837  FieldContainer<double> data_c_p_d(c, p, d1);
2838  FieldContainer<double> datainv_c_p_d(c, p, d1);
2839  FieldContainer<double> data_c_1_d(c, 1, d1);
2840  FieldContainer<double> datainv_c_1_d(c, 1, d1);
2841  FieldContainer<double> data_c_p_d_d(c, p, d1, d1);
2842  FieldContainer<double> datainv_c_p_d_d(c, p, d1, d1);
2843  FieldContainer<double> datainvtrn_c_p_d_d(c, p, d1, d1);
2844  FieldContainer<double> data_c_1_d_d(c, 1, d1, d1);
2845  FieldContainer<double> datainv_c_1_d_d(c, 1, d1, d1);
2846  FieldContainer<double> datainvtrn_c_1_d_d(c, 1, d1, d1);
2847  /***********************************************************************************************
2848  * Constant diagonal tensor: inputDataLeft(C,P) *
2849  **********************************************************************************************/
2850  for (int i=0; i<in_c_p_d_d.size(); i++) {
2851  in_c_p_d_d[i] = Teuchos::ScalarTraits<double>::random();
2852  }
2853  for (int i=0; i<data_c_p.size(); i++) {
2854  data_c_p[i] = Teuchos::ScalarTraits<double>::random();
2855  datainv_c_p[i] = 1.0 / data_c_p[i];
2856  }
2857  for (int i=0; i<data_c_1.size(); i++) {
2858  data_c_1[i] = Teuchos::ScalarTraits<double>::random();
2859  datainv_c_1[i] = 1.0 / data_c_1[i];
2860  }
2861  // Tensor values vary by point:
2862  art::matmatProductDataData<double>(out_c_p_d_d, data_c_p, in_c_p_d_d);
2863  art::matmatProductDataData<double>(out_c_p_d_d, datainv_c_p, out_c_p_d_d);
2864  rst::subtract(&out_c_p_d_d[0], &in_c_p_d_d[0], out_c_p_d_d.size());
2865  if (rst::vectorNorm(&out_c_p_d_d[0], out_c_p_d_d.size(), NORM_ONE) > zero) {
2866  *outStream << "\n\nINCORRECT matmatProductDataData (1): check scalar inverse property\n\n";
2867  errorFlag = -1000;
2868  }
2869  // Tensor values do not vary by point:
2870  art::matmatProductDataData<double>(out_c_p_d_d, data_c_1, in_c_p_d_d);
2871  art::matmatProductDataData<double>(out_c_p_d_d, datainv_c_1, out_c_p_d_d);
2872  rst::subtract(&out_c_p_d_d[0], &in_c_p_d_d[0], out_c_p_d_d.size());
2873  if (rst::vectorNorm(&out_c_p_d_d[0], out_c_p_d_d.size(), NORM_ONE) > zero) {
2874  *outStream << "\n\nINCORRECT matmatProductDataData (2): check scalar inverse property\n\n";
2875  errorFlag = -1000;
2876  }
2877  /***********************************************************************************************
2878  * Non-onstant diagonal tensor: inputDataLeft(C,P,D) *
2879  **********************************************************************************************/
2880  for (int i=0; i<in_c_p_d_d.size(); i++) {
2881  in_c_p_d_d[i] = Teuchos::ScalarTraits<double>::random();
2882  }
2883  for (int i=0; i<data_c_p_d.size(); i++) {
2884  data_c_p_d[i] = Teuchos::ScalarTraits<double>::random();
2885  datainv_c_p_d[i] = 1.0 / data_c_p_d[i];
2886  }
2887  for (int i=0; i<data_c_1_d.size(); i++) {
2888  data_c_1_d[i] = Teuchos::ScalarTraits<double>::random();
2889  datainv_c_1_d[i] = 1.0 / data_c_1_d[i];
2890  }
2891  // Tensor values vary by point:
2892  art::matmatProductDataData<double>(out_c_p_d_d, data_c_p_d, in_c_p_d_d);
2893  art::matmatProductDataData<double>(out_c_p_d_d, datainv_c_p_d, out_c_p_d_d);
2894  rst::subtract(&out_c_p_d_d[0], &in_c_p_d_d[0], out_c_p_d_d.size());
2895  if (rst::vectorNorm(&out_c_p_d_d[0], out_c_p_d_d.size(), NORM_ONE) > zero) {
2896  *outStream << "\n\nINCORRECT matmatProductDataData (3): check scalar inverse property\n\n";
2897  errorFlag = -1000;
2898  }
2899  // Tensor values do not vary by point
2900  art::matmatProductDataData<double>(out_c_p_d_d, data_c_1_d, in_c_p_d_d);
2901  art::matmatProductDataData<double>(out_c_p_d_d, datainv_c_1_d, out_c_p_d_d);
2902  rst::subtract(&out_c_p_d_d[0], &in_c_p_d_d[0], out_c_p_d_d.size());
2903  if (rst::vectorNorm(&out_c_p_d_d[0], out_c_p_d_d.size(), NORM_ONE) > zero) {
2904  *outStream << "\n\nINCORRECT matmatProductDataData (4): check scalar inverse property\n\n";
2905  errorFlag = -1000;
2906  }
2907  /***********************************************************************************************
2908  * Full tensor: inputDataLeft(C,P,D,D) *
2909  **********************************************************************************************/
2910  for (int i=0; i<in_c_p_d_d.size(); i++) {
2911  in_c_p_d_d[i] = Teuchos::ScalarTraits<double>::random();
2912  }
2913  for (int ic=0; ic < c; ic++) {
2914  for (int ip=0; ip < p; ip++) {
2915  for (int i=0; i<d1*d1; i++) {
2916  (&data_c_p_d_d(ic, ip, 0, 0))[i] = Teuchos::ScalarTraits<double>::random();
2917  }
2918  RealSpaceTools<double>::inverse(&datainv_c_p_d_d(ic, ip, 0, 0), &data_c_p_d_d(ic, ip, 0, 0), d1);
2919  }
2920  }
2921  for (int ic=0; ic < c; ic++) {
2922  for (int ip=0; ip < 1; ip++) {
2923  for (int i=0; i<d1*d1; i++) {
2924  (&data_c_1_d_d(ic, 0, 0, 0))[i] = Teuchos::ScalarTraits<double>::random();
2925  }
2926  RealSpaceTools<double>::inverse(&datainv_c_1_d_d(ic, 0, 0, 0), &data_c_1_d_d(ic, 0, 0, 0), d1);
2927  }
2928  }
2929  // Tensor values vary by point: test "N" and "T" options (no transpose/transpose)
2930  art::matmatProductDataData<double>(out_c_p_d_d, data_c_p_d_d, in_c_p_d_d);
2931  art::matmatProductDataData<double>(outi_c_p_d_d, datainv_c_p_d_d, out_c_p_d_d);
2932  rst::subtract(&outi_c_p_d_d[0], &in_c_p_d_d[0], outi_c_p_d_d.size());
2933  if (rst::vectorNorm(&outi_c_p_d_d[0], outi_c_p_d_d.size(), NORM_ONE) > zero) {
2934  *outStream << "\n\nINCORRECT matmatProductDataData (5): check matrix inverse property\n\n";
2935  errorFlag = -1000;
2936  }
2937  art::matmatProductDataData<double>(out_c_p_d_d, data_c_p_d_d, in_c_p_d_d, 't');
2938  art::matmatProductDataData<double>(outi_c_p_d_d, datainv_c_p_d_d, out_c_p_d_d, 't');
2939  rst::subtract(&outi_c_p_d_d[0], &in_c_p_d_d[0], outi_c_p_d_d.size());
2940  if (rst::vectorNorm(&outi_c_p_d_d[0], outi_c_p_d_d.size(), NORM_ONE) > zero) {
2941  *outStream << "\n\nINCORRECT matmatProductDataData (6): check matrix inverse property, w/ double transpose\n\n";
2942  errorFlag = -1000;
2943  }
2944  // Tensor values do not vary by point: test "N" and "T" options (no transpose/transpose)
2945  art::matmatProductDataData<double>(out_c_p_d_d, data_c_1_d_d, in_c_p_d_d);
2946  art::matmatProductDataData<double>(outi_c_p_d_d, datainv_c_1_d_d, out_c_p_d_d);
2947  rst::subtract(&outi_c_p_d_d[0], &in_c_p_d_d[0], outi_c_p_d_d.size());
2948  if (rst::vectorNorm(&outi_c_p_d_d[0], outi_c_p_d_d.size(), NORM_ONE) > zero) {
2949  *outStream << "\n\nINCORRECT matmatProductDataData (7): check matrix inverse property\n\n";
2950  errorFlag = -1000;
2951  }
2952  art::matmatProductDataData<double>(out_c_p_d_d, data_c_1_d_d, in_c_p_d_d, 't');
2953  art::matmatProductDataData<double>(outi_c_p_d_d, datainv_c_1_d_d, out_c_p_d_d, 't');
2954  rst::subtract(&outi_c_p_d_d[0], &in_c_p_d_d[0], outi_c_p_d_d.size());
2955  if (rst::vectorNorm(&outi_c_p_d_d[0], outi_c_p_d_d.size(), NORM_ONE) > zero) {
2956  *outStream << "\n\nINCORRECT matmatProductDataData (8): check matrix inverse property, w/ double transpose\n\n";
2957  errorFlag = -1000;
2958  }
2959  /***********************************************************************************************
2960  * Full tensor: inputDataLeft(C,P,D,D) test inverse transpose *
2961  **********************************************************************************************/
2962  for (int i=0; i<in_c_p_d_d.size(); i++) {
2963  in_c_p_d_d[i] = Teuchos::ScalarTraits<double>::random();
2964  }
2965  for (int ic=0; ic < c; ic++) {
2966  for (int ip=0; ip < p; ip++) {
2967  for (int i=0; i<d1*d1; i++) {
2968  (&data_c_p_d_d(ic, ip, 0, 0))[i] = Teuchos::ScalarTraits<double>::random();
2969  }
2970  RealSpaceTools<double>::inverse(&datainv_c_p_d_d(ic, ip, 0, 0), &data_c_p_d_d(ic, ip, 0, 0), d1);
2971  RealSpaceTools<double>::transpose(&datainvtrn_c_p_d_d(ic, ip, 0, 0), &datainv_c_p_d_d(ic, ip, 0, 0), d1);
2972  }
2973  }
2974  for (int ic=0; ic < c; ic++) {
2975  for (int ip=0; ip < 1; ip++) {
2976  for (int i=0; i<d1*d1; i++) {
2977  (&data_c_1_d_d(ic, 0, 0, 0))[i] = Teuchos::ScalarTraits<double>::random();
2978  }
2979  RealSpaceTools<double>::inverse(&datainv_c_1_d_d(ic, 0, 0, 0), &data_c_1_d_d(ic, 0, 0, 0), d1);
2980  RealSpaceTools<double>::transpose(&datainvtrn_c_1_d_d(ic, 0, 0, 0), &datainv_c_1_d_d(ic, 0, 0, 0), d1);
2981  }
2982  }
2983  // Tensor values vary by point:
2984  art::matmatProductDataData<double>(out_c_p_d_d, data_c_p_d_d, in_c_p_d_d);
2985  art::matmatProductDataData<double>(outi_c_p_d_d, datainvtrn_c_p_d_d, out_c_p_d_d, 't');
2986  rst::subtract(&outi_c_p_d_d[0], &in_c_p_d_d[0], outi_c_p_d_d.size());
2987  if (rst::vectorNorm(&outi_c_p_d_d[0], outi_c_p_d_d.size(), NORM_ONE) > zero) {
2988  *outStream << "\n\nINCORRECT matmatProductDataData (9): check matrix inverse transpose property\n\n";
2989  errorFlag = -1000;
2990  }
2991  // Tensor values do not vary by point
2992  art::matmatProductDataData<double>(out_c_p_d_d, data_c_1_d_d, in_c_p_d_d);
2993  art::matmatProductDataData<double>(outi_c_p_d_d, datainvtrn_c_1_d_d, out_c_p_d_d, 't');
2994  rst::subtract(&outi_c_p_d_d[0], &in_c_p_d_d[0], outi_c_p_d_d.size());
2995  if (rst::vectorNorm(&outi_c_p_d_d[0], outi_c_p_d_d.size(), NORM_ONE) > zero) {
2996  *outStream << "\n\nINCORRECT matmatProductDataData (10): check matrix inverse transpose property\n\n";
2997  errorFlag = -1000;
2998  }
2999  }// test 8.a scope
3000  *outStream \
3001  << "\n"
3002  << "===============================================================================\n"\
3003  << "| TEST 8.b: matmatProductDataData random tests: branch inputDataRight(P,D,D) |\n"\
3004  << "===============================================================================\n";
3005  /*
3006  * Test derived from Test 7.b
3007  */
3008  {// test 8.b scope
3009  int c=5, p=9, d1=3;
3010  double zero = INTREPID_TOL*10000.0;
3011 
3012  FieldContainer<double> in_p_d_d(p, d1, d1);
3013  FieldContainer<double> in_c_p_d_d(c, p, d1, d1);
3014  FieldContainer<double> out_c_p_d_d(c, p, d1, d1);
3015  FieldContainer<double> outi_c_p_d_d(c, p, d1, d1);
3016 
3017  FieldContainer<double> data_c_p(c, p);
3018  FieldContainer<double> datainv_c_p(c, p);
3019  FieldContainer<double> data_c_1(c, 1);
3020  FieldContainer<double> datainv_c_1(c, 1);
3021  FieldContainer<double> data_c_p_d(c, p, d1);
3022  FieldContainer<double> datainv_c_p_d(c, p, d1);
3023  FieldContainer<double> data_c_1_d(c, 1, d1);
3024  FieldContainer<double> datainv_c_1_d(c, 1, d1);
3025  FieldContainer<double> data_c_p_d_d(c, p, d1, d1);
3026  FieldContainer<double> datainv_c_p_d_d(c, p, d1, d1);
3027  FieldContainer<double> datainvtrn_c_p_d_d(c, p, d1, d1);
3028  FieldContainer<double> data_c_1_d_d(c, 1, d1, d1);
3029  FieldContainer<double> datainv_c_1_d_d(c, 1, d1, d1);
3030  FieldContainer<double> datainvtrn_c_1_d_d(c, 1, d1, d1);
3031  FieldContainer<double> data_c_p_one(c, p);
3032  FieldContainer<double> data_c_1_one(c, 1);
3033  /***********************************************************************************************
3034  * Constant diagonal tensor: inputData(C,P) *
3035  **********************************************************************************************/
3036  for (int i=0; i<in_p_d_d.size(); i++) {
3037  in_p_d_d[i] = Teuchos::ScalarTraits<double>::random();
3038  }
3039  for (int i=0; i<data_c_p.size(); i++) {
3040  data_c_p[i] = Teuchos::ScalarTraits<double>::random();
3041  datainv_c_p[i] = 1.0 / data_c_p[i];
3042  data_c_p_one[i] = 1.0;
3043  }
3044  for (int i=0; i<data_c_1.size(); i++) {
3045  data_c_1[i] = Teuchos::ScalarTraits<double>::random();
3046  datainv_c_1[i] = 1.0 / data_c_1[i];
3047  }
3048  // Tensor values vary by point
3049  art::matmatProductDataData<double>(in_c_p_d_d, data_c_p_one, in_p_d_d);
3050  art::matmatProductDataData<double>(out_c_p_d_d, data_c_p, in_p_d_d);
3051  art::matmatProductDataData<double>(outi_c_p_d_d, datainv_c_p, out_c_p_d_d);
3052  rst::subtract(&outi_c_p_d_d[0], &in_c_p_d_d[0], outi_c_p_d_d.size());
3053  if (rst::vectorNorm(&outi_c_p_d_d[0], outi_c_p_d_d.size(), NORM_ONE) > zero) {
3054  *outStream << "\n\nINCORRECT matmatProductDataData (11): check scalar inverse property\n\n";
3055  errorFlag = -1000;
3056  }
3057  // Tensor values do not vary by point
3058  art::matmatProductDataData<double>(in_c_p_d_d, data_c_p_one, in_p_d_d);
3059  art::matmatProductDataData<double>(out_c_p_d_d, data_c_1, in_p_d_d);
3060  art::matmatProductDataData<double>(outi_c_p_d_d, datainv_c_1, out_c_p_d_d);
3061  rst::subtract(&outi_c_p_d_d[0], &in_c_p_d_d[0], outi_c_p_d_d.size());
3062  if (rst::vectorNorm(&outi_c_p_d_d[0], outi_c_p_d_d.size(), NORM_ONE) > zero) {
3063  *outStream << "\n\nINCORRECT matmatProductDataData (12): check scalar inverse property\n\n";
3064  errorFlag = -1000;
3065  }
3066  /***********************************************************************************************
3067  * Non-constant diagonal tensor: inputData(C,P,D) *
3068  **********************************************************************************************/
3069  for (int i=0; i<in_p_d_d.size(); i++) {
3070  in_p_d_d[i] = Teuchos::ScalarTraits<double>::random();
3071  }
3072  for (int i=0; i<data_c_p_d.size(); i++) {
3073  data_c_p_d[i] = Teuchos::ScalarTraits<double>::random();
3074  datainv_c_p_d[i] = 1.0 / data_c_p_d[i];
3075  }
3076  for (int i=0; i<data_c_1_d.size(); i++) {
3077  data_c_1_d[i] = Teuchos::ScalarTraits<double>::random();
3078  datainv_c_1_d[i] = 1.0 / data_c_1_d[i];
3079  }
3080  // Tensor values vary by point:
3081  art::matmatProductDataData<double>(in_c_p_d_d, data_c_p_one, in_p_d_d);
3082  art::matmatProductDataData<double>(out_c_p_d_d, data_c_p_d, in_p_d_d);
3083  art::matmatProductDataData<double>(outi_c_p_d_d, datainv_c_p_d, out_c_p_d_d);
3084  rst::subtract(&outi_c_p_d_d[0], &in_c_p_d_d[0], outi_c_p_d_d.size());
3085  if (rst::vectorNorm(&outi_c_p_d_d[0], outi_c_p_d_d.size(), NORM_ONE) > zero) {
3086  *outStream << "\n\nINCORRECT matmatProductDataData (13): check scalar inverse property\n\n";
3087  errorFlag = -1000;
3088  }
3089  // Tensor values do not vary by point:
3090  art::matmatProductDataData<double>(in_c_p_d_d, data_c_p_one, in_p_d_d);
3091  art::matmatProductDataData<double>(out_c_p_d_d, data_c_1_d, in_p_d_d);
3092  art::matmatProductDataData<double>(outi_c_p_d_d, datainv_c_1_d, out_c_p_d_d);
3093  rst::subtract(&outi_c_p_d_d[0], &in_c_p_d_d[0], outi_c_p_d_d.size());
3094  if (rst::vectorNorm(&outi_c_p_d_d[0], outi_c_p_d_d.size(), NORM_ONE) > zero) {
3095  *outStream << "\n\nINCORRECT matmatProductDataData (14): check scalar inverse property\n\n";
3096  errorFlag = -1000;
3097  }
3098  /***********************************************************************************************
3099  * Full tensor: inputData(C,P,D,D) *
3100  **********************************************************************************************/
3101  for (int i=0; i<in_p_d_d.size(); i++) {
3102  in_p_d_d[i] = Teuchos::ScalarTraits<double>::random();
3103  }
3104  for (int ic=0; ic < c; ic++) {
3105  for (int ip=0; ip < p; ip++) {
3106  for (int i=0; i<d1*d1; i++) {
3107  (&data_c_p_d_d(ic, ip, 0, 0))[i] = Teuchos::ScalarTraits<double>::random();
3108  }
3109  RealSpaceTools<double>::inverse(&datainv_c_p_d_d(ic, ip, 0, 0), &data_c_p_d_d(ic, ip, 0, 0), d1);
3110  }
3111  }
3112  for (int ic=0; ic < c; ic++) {
3113  for (int ip=0; ip < 1; ip++) {
3114  for (int i=0; i<d1*d1; i++) {
3115  (&data_c_1_d_d(ic, 0, 0, 0))[i] = Teuchos::ScalarTraits<double>::random();
3116  }
3117  RealSpaceTools<double>::inverse(&datainv_c_1_d_d(ic, 0, 0, 0), &data_c_1_d_d(ic, 0, 0, 0), d1);
3118  }
3119  }
3120  // Tensor values vary by point: test "N" and "T" (no-transpose/transpose) options
3121  art::matmatProductDataData<double>(in_c_p_d_d, data_c_p_one, in_p_d_d);
3122  art::matmatProductDataData<double>(out_c_p_d_d, data_c_p_d_d, in_p_d_d);
3123  art::matmatProductDataData<double>(outi_c_p_d_d, datainv_c_p_d_d, out_c_p_d_d);
3124  rst::subtract(&outi_c_p_d_d[0], &in_c_p_d_d[0], outi_c_p_d_d.size());
3125  if (rst::vectorNorm(&outi_c_p_d_d[0], outi_c_p_d_d.size(), NORM_ONE) > zero) {
3126  *outStream << "\n\nINCORRECT matmatProductDataData (15): check matrix inverse property\n\n";
3127  errorFlag = -1000;
3128  }
3129  art::matmatProductDataData<double>(in_c_p_d_d, data_c_p_one, in_p_d_d);
3130  art::matmatProductDataData<double>(out_c_p_d_d, data_c_p_d_d, in_p_d_d, 't');
3131  art::matmatProductDataData<double>(outi_c_p_d_d, datainv_c_p_d_d, out_c_p_d_d, 't');
3132  rst::subtract(&outi_c_p_d_d[0], &in_c_p_d_d[0], outi_c_p_d_d.size());
3133  if (rst::vectorNorm(&outi_c_p_d_d[0], outi_c_p_d_d.size(), NORM_ONE) > zero) {
3134  *outStream << "\n\nINCORRECT matmatProductDataData (16): check matrix inverse property, w/ double transpose\n\n";
3135  errorFlag = -1000;
3136  }
3137  // Tensor values do not vary by point: test "N" and "T" (no-transpose/transpose) options
3138  art::matmatProductDataData<double>(in_c_p_d_d, data_c_p_one, in_p_d_d);
3139  art::matmatProductDataData<double>(out_c_p_d_d, data_c_1_d_d, in_p_d_d);
3140  art::matmatProductDataData<double>(outi_c_p_d_d, datainv_c_1_d_d, out_c_p_d_d);
3141  rst::subtract(&outi_c_p_d_d[0], &in_c_p_d_d[0], outi_c_p_d_d.size());
3142  if (rst::vectorNorm(&outi_c_p_d_d[0], outi_c_p_d_d.size(), NORM_ONE) > zero) {
3143  *outStream << "\n\nINCORRECT matmatProductDataData (17): check matrix inverse property\n\n";
3144  errorFlag = -1000;
3145  }
3146  art::matmatProductDataData<double>(in_c_p_d_d, data_c_p_one, in_p_d_d);
3147  art::matmatProductDataData<double>(out_c_p_d_d, data_c_1_d_d, in_p_d_d, 't');
3148  art::matmatProductDataData<double>(outi_c_p_d_d, datainv_c_1_d_d, out_c_p_d_d, 't');
3149  rst::subtract(&outi_c_p_d_d[0], &in_c_p_d_d[0], outi_c_p_d_d.size());
3150  if (rst::vectorNorm(&outi_c_p_d_d[0], outi_c_p_d_d.size(), NORM_ONE) > zero) {
3151  *outStream << "\n\nINCORRECT matmatProductDataData (18): check matrix inverse property, w/ double transpose\n\n";
3152  errorFlag = -1000;
3153  }
3154  /***********************************************************************************************
3155  * Full tensor: inputData(C,P,D,D) test inverse transpose *
3156  **********************************************************************************************/
3157  for (int i=0; i<in_p_d_d.size(); i++) {
3158  in_p_d_d[i] = Teuchos::ScalarTraits<double>::random();
3159  }
3160  for (int ic=0; ic < c; ic++) {
3161  for (int ip=0; ip < p; ip++) {
3162  for (int i=0; i<d1*d1; i++) {
3163  (&data_c_p_d_d(ic, ip, 0, 0))[i] = Teuchos::ScalarTraits<double>::random();
3164  }
3165  RealSpaceTools<double>::inverse(&datainv_c_p_d_d(ic, ip, 0, 0), &data_c_p_d_d(ic, ip, 0, 0), d1);
3166  RealSpaceTools<double>::transpose(&datainvtrn_c_p_d_d(ic, ip, 0, 0), &datainv_c_p_d_d(ic, ip, 0, 0), d1);
3167  }
3168  }
3169  for (int ic=0; ic < c; ic++) {
3170  for (int ip=0; ip < 1; ip++) {
3171  for (int i=0; i<d1*d1; i++) {
3172  (&data_c_1_d_d(ic, 0, 0, 0))[i] = Teuchos::ScalarTraits<double>::random();
3173  }
3174  RealSpaceTools<double>::inverse(&datainv_c_1_d_d(ic, 0, 0, 0), &data_c_1_d_d(ic, 0, 0, 0), d1);
3175  RealSpaceTools<double>::transpose(&datainvtrn_c_1_d_d(ic, 0, 0, 0), &datainv_c_1_d_d(ic, 0, 0, 0), d1);
3176  }
3177  }
3178  // Tensor values vary by point:
3179  art::matmatProductDataData<double>(in_c_p_d_d, data_c_p_one, in_p_d_d);
3180  art::matmatProductDataData<double>(out_c_p_d_d, data_c_p_d_d, in_p_d_d);
3181  art::matmatProductDataData<double>(outi_c_p_d_d, datainvtrn_c_p_d_d, out_c_p_d_d, 't');
3182  rst::subtract(&outi_c_p_d_d[0], &in_c_p_d_d[0], outi_c_p_d_d.size());
3183  if (rst::vectorNorm(&outi_c_p_d_d[0], outi_c_p_d_d.size(), NORM_ONE) > zero) {
3184  *outStream << "\n\nINCORRECT matmatProductDataData (19): check matrix inverse transpose property\n\n";
3185  errorFlag = -1000;
3186  }
3187  // Tensor values do not vary by point:
3188  art::matmatProductDataData<double>(in_c_p_d_d, data_c_p_one, in_p_d_d);
3189  art::matmatProductDataData<double>(out_c_p_d_d, data_c_1_d_d, in_p_d_d);
3190  art::matmatProductDataData<double>(outi_c_p_d_d, datainvtrn_c_1_d_d, out_c_p_d_d, 't');
3191  rst::subtract(&outi_c_p_d_d[0], &in_c_p_d_d[0], outi_c_p_d_d.size());
3192  if (rst::vectorNorm(&outi_c_p_d_d[0], outi_c_p_d_d.size(), NORM_ONE) > zero) {
3193  *outStream << "\n\nINCORRECT matmatProductDataData (20): check matrix inverse transpose property\n\n";
3194  errorFlag = -1000;
3195  }
3196  }// test 8.b scope
3197  }//try
3198 
3199  /*************************************************************************************************
3200  * Finish test *
3201  ************************************************************************************************/
3202 
3203  catch (const std::logic_error & err) {
3204  *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
3205  *outStream << err.what() << '\n';
3206  *outStream << "-------------------------------------------------------------------------------" << "\n\n";
3207  errorFlag = -1000;
3208  };
3209 
3210 
3211  if (errorFlag != 0)
3212  std::cout << "End Result: TEST FAILED\n";
3213  else
3214  std::cout << "End Result: TEST PASSED\n";
3215 
3216  // reset format state of std::cout
3217  std::cout.copyfmt(oldFormatState);
3218 
3219  return errorFlag;
3220 }
Implementation of basic linear algebra functionality in Euclidean space.
static void matvecProductDataData(ArrayOutData &outputData, const ArrayInDataLeft &inputDataLeft, const ArrayInDataRight &inputDataRight, const char transpose='N')
There are two use cases: (1) matrix-vector product of a rank-3 container inputDataRight with dimensio...
static void matvecProductDataField(ArrayOutFields &outputFields, const ArrayInData &inputData, const ArrayInFields &inputFields, const char transpose='N')
There are two use cases: (1) matrix-vector product of a rank-4 container inputFields with dimensions ...
Header file for utility class to provide multidimensional containers.
Header file for utility class to provide array tools, such as tensor contractions, etc.
static void matmatProductDataData(ArrayOutData &outputData, const ArrayInDataLeft &inputDataLeft, const ArrayInDataRight &inputDataRight, const char transpose='N')
There are two use cases: (1) matrix-matrix product of a rank-4 container inputDataRight with dimensio...
static void crossProductDataField(ArrayOutFields &outputFields, const ArrayInData &inputData, const ArrayInFields &inputFields)
There are two use cases: (1) cross product of a rank-4 container inputFields with dimensions (C...
static void matmatProductDataField(ArrayOutFields &outputFields, const ArrayInData &inputData, const ArrayInFields &inputFields, const char transpose='N')
There are two use cases: (1) matrix-matrix product of a rank-5 container inputFields with dimensions ...
Header file for classes providing basic linear algebra functionality in 1D, 2D and 3D...
static void crossProductDataData(ArrayOutData &outputData, const ArrayInDataLeft &inputDataLeft, const ArrayInDataRight &inputDataRight)
There are two use cases: (1) cross product of a rank-3 container inputDataRight with dimensions (C...
static void outerProductDataData(ArrayOutData &outputData, const ArrayInDataLeft &inputDataLeft, const ArrayInDataRight &inputDataRight)
There are two use cases: (1) outer product of a rank-3 container inputDataRight with dimensions (C...
static void outerProductDataField(ArrayOutFields &outputFields, const ArrayInData &inputData, const ArrayInFields &inputFields)
There are two use cases: (1) outer product of a rank-4 container inputFields with dimensions (C...