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
52#include "Teuchos_oblackholestream.hpp"
53#include "Teuchos_RCP.hpp"
54#include "Teuchos_ScalarTraits.hpp"
55
56using namespace std;
57using 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
72int 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}
Header file for utility class to provide array tools, such as tensor contractions,...
Header file for utility class to provide multidimensional containers.
Header file for classes providing basic linear algebra functionality in 1D, 2D and 3D.
Utility class that provides methods for higher-order algebraic manipulation of user-defined arrays,...
static void matvecProductDataData(ArrayOutData &outputData, const ArrayInDataLeft &inputDataLeft, const ArrayInDataRight &inputDataRight, const char transpose='N')
There are two use cases: (1) matrix-vector product of a rank-3 container inputDataRight with dimensio...
static void 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,...
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 ...
static void crossProductDataField(ArrayOutFields &outputFields, const ArrayInData &inputData, const ArrayInFields &inputFields)
There are two use cases: (1) cross product of a rank-4 container inputFields with dimensions (C,...
static void matmatProductDataData(ArrayOutData &outputData, const ArrayInDataLeft &inputDataLeft, const ArrayInDataRight &inputDataRight, const char transpose='N')
There are two use cases: (1) matrix-matrix product of a rank-4 container inputDataRight with dimensio...
static void 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 ...
static void outerProductDataData(ArrayOutData &outputData, const ArrayInDataLeft &inputDataLeft, const ArrayInDataRight &inputDataRight)
There are two use cases: (1) outer product of a rank-3 container inputDataRight with dimensions (C,...
static void 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,...
Implementation of a templated lexicographical container for a multi-indexed scalar quantity....
Implementation of basic linear algebra functionality in Euclidean space.
static void inverse(Scalar *inverseMat, const Scalar *inMat, const size_t dim)
Computes inverse of the square matrix inMat of size dim by dim.
static void transpose(Scalar *transposeMat, const Scalar *inMat, const size_t dim)
Computes transpose of the square matrix inMat of size dim by dim.