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"
50 #include "Intrepid_FieldContainer.hpp"
51 #include "Intrepid_RealSpaceTools.hpp"
52 #include "Teuchos_oblackholestream.hpp"
53 #include "Teuchos_RCP.hpp"
54 #include "Teuchos_ScalarTraits.hpp"
55 
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 }
Intrepid::RealSpaceTools
Implementation of basic linear algebra functionality in Euclidean space.
Definition: Intrepid_RealSpaceTools.hpp:65
Intrepid::ArrayTools::matmatProductDataData
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...
Definition: Intrepid_ArrayToolsDefTensor.hpp:1976
Intrepid_FieldContainer.hpp
Header file for utility class to provide multidimensional containers.
Intrepid_ArrayTools.hpp
Header file for utility class to provide array tools, such as tensor contractions, etc.
Intrepid::ArrayTools::crossProductDataField
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...
Definition: Intrepid_ArrayToolsDefTensor.hpp:55
Intrepid_RealSpaceTools.hpp
Header file for classes providing basic linear algebra functionality in 1D, 2D and 3D...
Intrepid::ArrayTools
Utility class that provides methods for higher-order algebraic manipulation of user-defined arrays...
Definition: Intrepid_ArrayTools.hpp:67
Intrepid::ArrayTools::crossProductDataData
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...
Definition: Intrepid_ArrayToolsDefTensor.hpp:238
Intrepid::ArrayTools::matmatProductDataField
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 ...
Definition: Intrepid_ArrayToolsDefTensor.hpp:1497
Intrepid::ArrayTools::outerProductDataData
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...
Definition: Intrepid_ArrayToolsDefTensor.hpp:532
Intrepid::ArrayTools::outerProductDataField
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...
Definition: Intrepid_ArrayToolsDefTensor.hpp:417
Intrepid::ArrayTools::matvecProductDataData
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...
Definition: Intrepid_ArrayToolsDefTensor.hpp:1085
Intrepid::ArrayTools::matvecProductDataField
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 ...
Definition: Intrepid_ArrayToolsDefTensor.hpp:642
Generated by   doxygen 1.8.5