Intrepid
Intrepid_ArrayToolsDefTensor.hpp
Go to the documentation of this file.
1// @HEADER
2// ************************************************************************
3//
4// Intrepid Package
5// Copyright (2007) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Pavel Bochev (pbboche@sandia.gov)
38// Denis Ridzal (dridzal@sandia.gov), or
39// Kara Peterson (kjpeter@sandia.gov)
40//
41// ************************************************************************
42// @HEADER
43
49namespace Intrepid {
50
51 template<class Scalar,
52 class ArrayOutFields,
53 class ArrayInData,
54 class ArrayInFields>
55 void ArrayTools::crossProductDataField(ArrayOutFields & outputFields,
56 const ArrayInData & inputData,
57 const ArrayInFields & inputFields){
58#ifdef HAVE_INTREPID_DEBUG
59 std::string errmsg = ">>> ERROR (ArrayTools::crossProductDataField):";
60 /*
61 * Check array rank and spatial dimension range (if applicable)
62 * (1) inputData(C,P,D);
63 * (2) inputFields(C,F,P,D) or (F,P,D);
64 * (3) outputFields(C,F,P,D) in 3D, or (C,F,P) in 2D
65 */
66 // (1) inputData is (C, P, D) and 2 <= D <= 3 is required
67 TEUCHOS_TEST_FOR_EXCEPTION( !requireRankRange(errmsg, inputData, 3, 3),
68 std::invalid_argument, errmsg);
69 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, inputData, 2, 2,3),
70 std::invalid_argument, errmsg);
71 // (2) inputFields is (C, F, P, D) or (F, P, D) and 2 <= (D=dimension(rank - 1)) <= 3 is required.
72 TEUCHOS_TEST_FOR_EXCEPTION( !requireRankRange(errmsg, inputFields, 3,4), std::invalid_argument, errmsg);
73 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, inputFields, getrank(inputFields)-1, 2,3),
74 std::invalid_argument, errmsg);
75 // (3) outputFields is (C,F,P,D) in 3D and (C,F,P) in 2D => rank = inputData.dimension(2) + 1
76 TEUCHOS_TEST_FOR_EXCEPTION( !requireRankRange(errmsg, outputFields, inputData.dimension(2)+1, inputData.dimension(2)+1),
77 std::invalid_argument, errmsg);
78 /*
79 * Dimension cross-checks:
80 * (1) inputData vs. inputFields
81 * (2) outputFields vs. inputData
82 * (3) outputFields vs. inputFields
83 *
84 * Cross-check (1):
85 */
86 if( getrank(inputFields) == 4) {
87 // inputData(C,P,D) vs. inputFields(C,F,P,D): dimensions C, P, D must match
88 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, inputData, 0,1,2, inputFields, 0,2,3),
89 std::invalid_argument, errmsg);
90 }
91 else{
92 // inputData(C,P,D) vs. inputFields(F,P,D): dimensions P, D must match
93 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, inputData, 1,2, inputFields, 1,2),
94 std::invalid_argument, errmsg);
95 }
96 /*
97 * Cross-check (2):
98 */
99 if(inputData.dimension(2) == 2) {
100 // in 2D: outputFields(C,F,P) vs. inputData(C,P,D): dimensions C,P must match
101 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, outputFields, 0,2, inputData, 0,1),
102 std::invalid_argument, errmsg);
103 }
104 else{
105 // in 3D: outputFields(C,F,P,D) vs. inputData(C,P,D): dimensions C,P,D must match
106 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, outputFields, 0,2,3, inputData, 0,1,2),
107 std::invalid_argument, errmsg);
108 }
109 /*
110 * Cross-check (3):
111 */
112 if(inputData.dimension(2) == 2) {
113 // In 2D:
114 if(getrank(inputFields) == 4){
115 // and rank-4 inputFields: outputFields(C,F,P) vs. inputFields(C,F,P,D): dimensions C,F,P must match
116 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, outputFields, 0,1,2, inputFields, 0,1,2),
117 std::invalid_argument, errmsg);
118 }
119 else{
120 // and rank-3 inputFields: outputFields(C,F,P) vs. inputFields(F,P,D): dimensions F,P must match
121 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, outputFields, 1,2, inputFields, 0,1),
122 std::invalid_argument, errmsg);
123 }
124 }
125 else{
126 // In 3D:
127 if(getrank(inputFields) == 4){
128 // and rank-4 inputFields: outputFields(C,F,P,D) vs. inputFields(C,F,P,D): all dimensions C,F,P,D must match
129 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, outputFields, inputFields),
130 std::invalid_argument, errmsg);
131 }
132 else{
133 // and rank-3 inputFields: outputFields(C,F,P,D) vs. inputFields(F,P,D): dimensions F,P,D must match
134 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, outputFields, 1,2,3, inputFields, 0,1,2),
135 std::invalid_argument, errmsg);
136 }
137 }
138#endif
139ArrayWrapper<Scalar,ArrayOutFields, Rank<ArrayOutFields >::value, false>outputFieldsWrap(outputFields);
140ArrayWrapper<Scalar,ArrayInData, Rank<ArrayInData >::value, true>inputDataWrap(inputData);
141ArrayWrapper<Scalar,ArrayInFields, Rank<ArrayInFields >::value, true>inputFieldsWrap(inputFields);
142
143 // 3D cross product
144 if(inputData.dimension(2) == 3) {
145
146 // inputFields is (C,F,P,D)
147 if(getrank(inputFields) == 4){
148
149 for(size_t cell = 0; cell < static_cast<size_t>(outputFields.dimension(0)); cell++){
150 for(size_t field = 0; field < static_cast<size_t>(outputFields.dimension(1)); field++){
151 for(size_t point = 0; point < static_cast<size_t>(outputFields.dimension(2)); point++){
152 //
153 outputFieldsWrap(cell, field, point, 0) = \
154 inputDataWrap(cell, point, 1)*inputFieldsWrap(cell, field, point, 2) -
155 inputDataWrap(cell, point, 2)*inputFieldsWrap(cell, field, point, 1);
156 //
157 outputFieldsWrap(cell, field, point, 1) = \
158 inputDataWrap(cell, point, 2)*inputFieldsWrap(cell, field, point, 0) -
159 inputDataWrap(cell, point, 0)*inputFieldsWrap(cell, field, point, 2);
160 //
161 outputFieldsWrap(cell, field, point, 2) = \
162 inputDataWrap(cell, point, 0)*inputFieldsWrap(cell, field, point, 1) -
163 inputDataWrap(cell, point, 1)*inputFieldsWrap(cell, field, point, 0);
164 }// point
165 }// field
166 } // cell
167 }// rank = 4
168 // inputFields is (F,P,D)
169 else if(getrank(inputFields) == 3){
170
171 for(size_t cell = 0; cell < static_cast<size_t>(outputFields.dimension(0)); cell++){
172 for(size_t field = 0; field < static_cast<size_t>(outputFields.dimension(1)); field++){
173 for(size_t point = 0; point < static_cast<size_t>(outputFields.dimension(2)); point++){
174 //
175 outputFieldsWrap(cell, field, point, 0) = \
176 inputDataWrap(cell, point, 1)*inputFieldsWrap(field, point, 2) -
177 inputDataWrap(cell, point, 2)*inputFieldsWrap(field, point, 1);
178 //
179 outputFieldsWrap(cell, field, point, 1) = \
180 inputDataWrap(cell, point, 2)*inputFieldsWrap(field, point, 0) -
181 inputDataWrap(cell, point, 0)*inputFieldsWrap(field, point, 2);
182 //
183 outputFieldsWrap(cell, field, point, 2) = \
184 inputDataWrap(cell, point, 0)*inputFieldsWrap(field, point, 1) -
185 inputDataWrap(cell, point, 1)*inputFieldsWrap(field, point, 0);
186 }// point
187 }// field
188 } // cell
189 }// rank = 3
190 else{
191 TEUCHOS_TEST_FOR_EXCEPTION( true, std::invalid_argument,
192 ">>> ERROR (ArrayTools::crossProductDataField): inputFields rank 3 or 4 required.")
193 }
194 }
195 // 2D cross product
196 else if(inputData.dimension(2) == 2){
197
198 // inputFields is (C,F,P,D)
199 if(getrank(inputFields) == 4){
200
201 for(size_t cell = 0; cell < static_cast<size_t>(outputFields.dimension(0)); cell++){
202 for(size_t field = 0; field < static_cast<size_t>(outputFields.dimension(1)); field++){
203 for(size_t point = 0; point < static_cast<size_t>(outputFields.dimension(2)); point++){
204 outputFieldsWrap(cell, field, point) = \
205 inputDataWrap(cell, point, 0)*inputFieldsWrap(cell, field, point, 1) -
206 inputDataWrap(cell, point, 1)*inputFieldsWrap(cell, field, point, 0);
207 }// point
208 }// field
209 } // cell
210 }// rank = 4
211 // inputFields is (F,P,D)
212 else if(getrank(inputFields) == 3) {
213
214 for(size_t cell = 0; cell < static_cast<size_t>(outputFields.dimension(0)); cell++){
215 for(size_t field = 0; field < static_cast<size_t>(outputFields.dimension(1)); field++){
216 for(size_t point = 0; point < static_cast<size_t>(outputFields.dimension(2)); point++){
217 outputFieldsWrap(cell, field, point) = \
218 inputDataWrap(cell, point, 0)*inputFieldsWrap(field, point, 1) -
219 inputDataWrap(cell, point, 1)*inputFieldsWrap(field, point, 0);
220 }// point
221 }// field
222 } // cell
223 }// rank = 3
224 }
225 // Error: wrong dimension
226 else {
227 TEUCHOS_TEST_FOR_EXCEPTION( true, std::invalid_argument,
228 ">>> ERROR (ArrayTools::crossProductDataField): spatial dimension 2 or 3 required.")
229 }
230 }
231
232
233
234 template<class Scalar,
235 class ArrayOutData,
236 class ArrayInDataLeft,
237 class ArrayInDataRight>
238 void ArrayTools::crossProductDataData(ArrayOutData & outputData,
239 const ArrayInDataLeft & inputDataLeft,
240 const ArrayInDataRight & inputDataRight){
241
242
243#ifdef HAVE_INTREPID_DEBUG
244 std::string errmsg = ">>> ERROR (ArrayTools::crossProductDataData):";
245 /*
246 * Check array rank and spatial dimension range (if applicable)
247 * (1) inputDataLeft(C,P,D);
248 * (2) inputDataRight(C,P,D) or (P,D);
249 * (3) outputData(C,P,D) in 3D, or (C,P) in 2D
250 */
251 // (1) inputDataLeft is (C, P, D) and 2 <= D <= 3 is required
252 TEUCHOS_TEST_FOR_EXCEPTION( !requireRankRange(errmsg, inputDataLeft, 3,3),
253 std::invalid_argument, errmsg);
254 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, inputDataLeft, 2, 2,3),
255 std::invalid_argument, errmsg);
256 // (2) inputDataRight is (C, P, D) or (P, D) and 2 <= (D=dimension(rank - 1)) <= 3 is required.
257 TEUCHOS_TEST_FOR_EXCEPTION( !requireRankRange(errmsg, inputDataRight, 2,3),
258 std::invalid_argument, errmsg);
259 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, inputDataRight, getrank(inputDataRight)-1, 2,3),
260 std::invalid_argument, errmsg);
261 // (3) outputData is (C,P,D) in 3D and (C,P) in 2D => rank = inputDataLeft.dimension(2)
262 TEUCHOS_TEST_FOR_EXCEPTION( !requireRankRange(errmsg, outputData, inputDataLeft.dimension(2), inputDataLeft.dimension(2)),
263 std::invalid_argument, errmsg);
264 /*
265 * Dimension cross-checks:
266 * (1) inputDataLeft vs. inputDataRight
267 * (2) outputData vs. inputDataLeft
268 * (3) outputData vs. inputDataRight
269 *
270 * Cross-check (1):
271 */
272 if( getrank(inputDataRight) == 3) {
273 // inputDataLeft(C,P,D) vs. inputDataRight(C,P,D): all dimensions C, P, D must match
274 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, inputDataLeft, inputDataRight),
275 std::invalid_argument, errmsg);
276 }
277 // inputDataLeft(C, P,D) vs. inputDataRight(P,D): dimensions P, D must match
278 else{
279 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, inputDataLeft, 1,2, inputDataRight, 0,1),
280 std::invalid_argument, errmsg);
281 }
282 /*
283 * Cross-check (2):
284 */
285 if(inputDataLeft.dimension(2) == 2){
286 // in 2D: outputData(C,P) vs. inputDataLeft(C,P,D): dimensions C, P must match
287 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, inputDataLeft, 0,1, outputData, 0,1),
288 std::invalid_argument, errmsg);
289 }
290 else{
291 // in 3D: outputData(C,P,D) vs. inputDataLeft(C,P,D): all dimensions C, P, D must match
292 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, inputDataLeft, outputData),
293 std::invalid_argument, errmsg);
294 }
295 /*
296 * Cross-check (3):
297 */
298 if(inputDataLeft.dimension(2) == 2) {
299 // In 2D:
300 if(getrank(inputDataRight) == 3){
301 // and rank-3 inputDataRight: outputData(C,P) vs. inputDataRight(C,P,D): dimensions C,P must match
302 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, outputData, 0,1, inputDataRight, 0,1),
303 std::invalid_argument, errmsg);
304 }
305 else{
306 // and rank-2 inputDataRight: outputData(C,P) vs. inputDataRight(P,D): dimension P must match
307 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, outputData, 1, inputDataRight, 0),
308 std::invalid_argument, errmsg);
309 }
310 }
311 else{
312 // In 3D:
313 if(getrank(inputDataRight) == 3){
314 // and rank-3 inputDataRight: outputData(C,P,D) vs. inputDataRight(C,P,D): all dimensions C,P,D must match
315 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, outputData, inputDataRight),
316 std::invalid_argument, errmsg);
317 }
318 else{
319 // and rank-2 inputDataRight: outputData(C,P,D) vs. inputDataRight(P,D): dimensions P, D must match
320 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, outputData, 1,2, inputDataRight, 0,1),
321 std::invalid_argument, errmsg);
322 }
323 }
324#endif
325
326
327ArrayWrapper<Scalar,ArrayOutData, Rank<ArrayOutData >::value, false>outputDataWrap(outputData);
328ArrayWrapper<Scalar,ArrayInDataLeft, Rank<ArrayInDataLeft >::value, true>inputDataLeftWrap(inputDataLeft);
329ArrayWrapper<Scalar,ArrayInDataRight, Rank<ArrayInDataRight >::value, true>inputDataRightWrap(inputDataRight);
330 // 3D cross product
331 if(inputDataLeft.dimension(2) == 3) {
332
333 // inputDataRight is (C,P,D)
334 if(getrank(inputDataRight) == 3){
335
336 for(size_t cell = 0; cell < (size_t)inputDataLeft.dimension(0); cell++){
337 for(size_t point = 0; point < (size_t)inputDataLeft.dimension(1); point++){
338 //
339 outputDataWrap(cell, point, 0) = \
340 inputDataLeftWrap(cell, point, 1)*inputDataRightWrap(cell, point, 2) -
341 inputDataLeftWrap(cell, point, 2)*inputDataRightWrap(cell, point, 1);
342 //
343 outputDataWrap(cell, point, 1) = \
344 inputDataLeftWrap(cell, point, 2)*inputDataRightWrap(cell, point, 0) -
345 inputDataLeftWrap(cell, point, 0)*inputDataRightWrap(cell, point, 2);
346 //
347 outputDataWrap(cell, point, 2) = \
348 inputDataLeftWrap(cell, point, 0)*inputDataRightWrap(cell, point, 1) -
349 inputDataLeftWrap(cell, point, 1)*inputDataRightWrap(cell, point, 0);
350 }// point
351 } // cell
352 }// rank = 3
353 // inputDataRight is (P,D)
354 else if(getrank(inputDataRight) == 2){
355
356 for(size_t cell = 0; cell < (size_t)inputDataLeft.dimension(0); cell++){
357 for(size_t point = 0; point < (size_t)inputDataLeft.dimension(1); point++){
358 //
359 outputDataWrap(cell, point, 0) = \
360 inputDataLeftWrap(cell, point, 1)*inputDataRightWrap(point, 2) -
361 inputDataLeftWrap(cell, point, 2)*inputDataRightWrap(point, 1);
362 //
363 outputDataWrap(cell, point, 1) = \
364 inputDataLeftWrap(cell, point, 2)*inputDataRightWrap(point, 0) -
365 inputDataLeftWrap(cell, point, 0)*inputDataRightWrap(point, 2);
366 //
367 outputDataWrap(cell, point, 2) = \
368 inputDataLeftWrap(cell, point, 0)*inputDataRightWrap(point, 1) -
369 inputDataLeftWrap(cell, point, 1)*inputDataRightWrap(point, 0);
370 }// point
371 } // cell
372 }// rank = 2
373 else{
374 TEUCHOS_TEST_FOR_EXCEPTION( true, std::invalid_argument,
375 ">>> ERROR (ArrayTools::crossProductDataData): inputDataRight rank 2 or 3 required.")
376 }
377 }
378 // 2D cross product
379 else if(inputDataLeft.dimension(2) == 2){
380
381 // inputDataRight is (C,P,D)
382 if(getrank(inputDataRight) == 3){
383
384 for(size_t cell = 0; cell < (size_t)inputDataLeft.dimension(0); cell++){
385 for(size_t point = 0; point < (size_t)inputDataLeft.dimension(1); point++){
386 outputDataWrap(cell, point) = \
387 inputDataLeftWrap(cell, point, 0)*inputDataRightWrap(cell, point, 1) -
388 inputDataLeftWrap(cell, point, 1)*inputDataRightWrap(cell, point, 0);
389 }// point
390 } // cell
391 }// rank = 3
392 // inputDataRight is (P,D)
393 else if(getrank(inputDataRight) == 2) {
394
395 for(size_t cell = 0; cell < (size_t)inputDataLeft.dimension(0); cell++){
396 for(size_t point = 0; point < (size_t)inputDataLeft.dimension(1); point++){
397 outputDataWrap(cell, point) = \
398 inputDataLeftWrap(cell, point, 0)*inputDataRightWrap(point, 1) -
399 inputDataLeftWrap(cell, point, 1)*inputDataRightWrap(point, 0);
400 }// point
401 } // cell
402 }// rank = 2
403 }
404 // Error: wrong dimension
405 else {
406 TEUCHOS_TEST_FOR_EXCEPTION( true, std::invalid_argument,
407 ">>> ERROR (ArrayTools::crossProductDataData): spatial dimension 2 or 3 required.")
408 }
409 }
410
411
412
413 template<class Scalar,
414 class ArrayOutFields,
415 class ArrayInData,
416 class ArrayInFields>
417 void ArrayTools::outerProductDataField(ArrayOutFields & outputFields,
418 const ArrayInData & inputData,
419 const ArrayInFields & inputFields){
420
421#ifdef HAVE_INTREPID_DEBUG
422 std::string errmsg = ">>> ERROR (ArrayTools::outerProductDataField):";
423 /*
424 * Check array rank and spatial dimension range (if applicable)
425 * (1) inputData(C,P,D);
426 * (2) inputFields(C,F,P,D) or (F,P,D);
427 * (3) outputFields(C,F,P,D,D)
428 */
429 // (1) inputData is (C, P, D) and 2 <= D <= 3 is required
430 TEUCHOS_TEST_FOR_EXCEPTION( !requireRankRange(errmsg, inputData, 3,3),
431 std::invalid_argument, errmsg);
432 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, inputData, 2, 2,3),
433 std::invalid_argument, errmsg);
434 // (2) inputFields is (C, F, P, D) or (F, P, D) and 2 <= (D=dimension(rank - 1)) <= 3 is required.
435 TEUCHOS_TEST_FOR_EXCEPTION( !requireRankRange(errmsg, inputFields, 3,4), std::invalid_argument, errmsg);
436 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, inputFields, getrank(inputFields)-1, 2,3),
437 std::invalid_argument, errmsg);
438 // (3) outputFields is (C,F,P,D,D)
439 TEUCHOS_TEST_FOR_EXCEPTION( !requireRankRange(errmsg, outputFields, 5,5), std::invalid_argument, errmsg);
440 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, outputFields, 3, 2,3),
441 std::invalid_argument, errmsg);
442 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, outputFields, 4, 2,3),
443 std::invalid_argument, errmsg);
444 /*
445 * Dimension cross-checks:
446 * (1) inputData vs. inputFields
447 * (2) outputFields vs. inputData
448 * (3) outputFields vs. inputFields
449 *
450 * Cross-check (2): outputFields(C,F,P,D,D) vs. inputData(C,P,D): dimensions C, P, D must match
451 */
452 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg,
453 outputFields, 0,2,3,4,
454 inputData, 0,1,2,2),
455 std::invalid_argument, errmsg);
456 /*
457 * Cross-checks (1,3):
458 */
459 if( getrank(inputFields) == 4) {
460 // Cross-check (1): inputData(C,P,D) vs. inputFields(C,F,P,D): dimensions C, P, D must match
461 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg,
462 inputData, 0,1,2,
463 inputFields, 0,2,3),
464 std::invalid_argument, errmsg);
465 // Cross-check (3): outputFields(C,F,P,D,D) vs. inputFields(C,F,P,D): dimensions C, F, P, D must match
466 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg,
467 outputFields, 0,1,2,3,4,
468 inputFields, 0,1,2,3,3),
469 std::invalid_argument, errmsg);
470 }
471 else{
472 // Cross-check (1): inputData(C,P,D) vs. inputFields(F,P,D): dimensions P, D must match
473 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg,
474 inputData, 1,2,
475 inputFields, 1,2),
476 std::invalid_argument, errmsg);
477 // Cross-check (3): outputFields(C,F,P,D,D) vs. inputFields(F,P,D): dimensions F, P, D must match
478 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg,
479 outputFields, 1,2,3,4,
480 inputFields, 0,1,2,2),
481 std::invalid_argument, errmsg);
482 }
483#endif
484ArrayWrapper<Scalar,ArrayOutFields, Rank<ArrayOutFields >::value, false>outputFieldsWrap(outputFields);
485ArrayWrapper<Scalar,ArrayInData, Rank<ArrayInData >::value, true>inputDataWrap(inputData);
486ArrayWrapper<Scalar,ArrayInFields, Rank<ArrayInFields >::value, true>inputFieldsWrap(inputFields);
487
488 // inputFields is (C,F,P,D)
489 if(getrank(inputFields) == 4){
490
491 for(size_t cell = 0; cell < (size_t)outputFields.dimension(0); cell++){
492 for(size_t field = 0; field < (size_t)outputFields.dimension(1); field++){
493 for(size_t point = 0; point < (size_t)outputFields.dimension(2); point++){
494 for(size_t row = 0; row < (size_t)outputFields.dimension(3); row++){
495 for(size_t col = 0; col < (size_t)outputFields.dimension(4); col++){
496 outputFieldsWrap(cell, field, point, row, col) = \
497 inputDataWrap(cell, point, row)*inputFieldsWrap(cell, field, point, col);
498 }// col
499 }// row
500 }// point
501 }// field
502 } // cell
503 }// rank = 4
504 // inputFields is (F,P,D)
505 else if(getrank(inputFields) == 3){
506
507 for(size_t cell = 0; cell < (size_t)outputFields.dimension(0); cell++){
508 for(size_t field = 0; field < (size_t)outputFields.dimension(1); field++){
509 for(size_t point = 0; point < (size_t)outputFields.dimension(2); point++){
510 for(size_t row = 0; row < (size_t)outputFields.dimension(3); row++){
511 for(size_t col = 0; col < (size_t)outputFields.dimension(4); col++){
512 outputFieldsWrap(cell, field, point, row, col) = \
513 inputDataWrap(cell, point, row)*inputFieldsWrap(field, point, col);
514 }// col
515 }// row
516 }// point
517 }// field
518 } // cell
519 }// rank = 3
520 else{
521 TEUCHOS_TEST_FOR_EXCEPTION( true, std::invalid_argument,
522 ">>> ERROR (ArrayTools::outerProductDataField): inputFields rank 3 or 4 required.")
523 }
524 }
525
526
527
528 template<class Scalar,
529 class ArrayOutData,
530 class ArrayInDataLeft,
531 class ArrayInDataRight>
532 void ArrayTools::outerProductDataData(ArrayOutData & outputData,
533 const ArrayInDataLeft & inputDataLeft,
534 const ArrayInDataRight & inputDataRight){
535
536#ifdef HAVE_INTREPID_DEBUG
537 std::string errmsg = ">>> ERROR (ArrayTools::outerProductDataData):";
538 /*
539 * Check array rank and spatial dimension range (if applicable)
540 * (1) inputDataLeft(C,P,D);
541 * (2) inputDataRight(C,P,D) or (P,D);
542 * (3) outputData(C,P,D,D)
543 */
544 // (1) inputDataLeft is (C, P, D) and 2 <= D <= 3 is required
545 TEUCHOS_TEST_FOR_EXCEPTION( !requireRankRange(errmsg, inputDataLeft, 3,3),
546 std::invalid_argument, errmsg);
547 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, inputDataLeft, 2, 2,3),
548 std::invalid_argument, errmsg);
549 // (2) inputDataRight is (C, P, D) or (P, D) and 2 <= (D=dimension(rank - 1)) <= 3 is required.
550 TEUCHOS_TEST_FOR_EXCEPTION( !requireRankRange(errmsg, inputDataRight, 2,3),
551 std::invalid_argument, errmsg);
552 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, inputDataRight, getrank(inputDataRight)-1, 2,3),
553 std::invalid_argument, errmsg);
554 // (3) outputData is (C,P,D,D)
555 TEUCHOS_TEST_FOR_EXCEPTION( !requireRankRange(errmsg, outputData, 4, 4), std::invalid_argument, errmsg);
556 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, outputData, 2, 2,3),
557 std::invalid_argument, errmsg);
558 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, outputData, 3, 2,3),
559 std::invalid_argument, errmsg);
560 /*
561 * Dimension cross-checks:
562 * (1) inputDataLeft vs. inputDataRight
563 * (2) outputData vs. inputDataLeft
564 * (3) outputData vs. inputDataRight
565 *
566 * Cross-check (2): outputData(C,P,D,D) vs. inputDataLeft(C,P,D): dimensions C, P, D must match
567 */
568 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg,
569 outputData, 0,1,2,3,
570 inputDataLeft, 0,1,2,2),
571 std::invalid_argument, errmsg);
572 /*
573 * Cross-checks (1,3):
574 */
575 if( getrank(inputDataRight) == 3) {
576 // Cross-check (1): inputDataLeft(C,P,D) vs. inputDataRight(C,P,D): all dimensions C, P, D must match
577 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, inputDataLeft, inputDataRight),
578 std::invalid_argument, errmsg);
579 // Cross-check (3): outputData(C,P,D,D) vs. inputDataRight(C,P,D): dimensions C, P, D must match
580 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg,
581 outputData, 0,1,2,3,
582 inputDataRight, 0,1,2,2),
583 std::invalid_argument, errmsg);
584 }
585 else{
586 // Cross-check (1): inputDataLeft(C,P,D) vs. inputDataRight(P,D): dimensions P, D must match
587 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg,
588 inputDataLeft, 1,2,
589 inputDataRight, 0,1),
590 std::invalid_argument, errmsg);
591 // Cross-check (3): outputData(C,P,D,D) vs. inputDataRight(P,D): dimensions P, D must match
592 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg,
593 outputData, 1,2,3,
594 inputDataRight, 0,1,1),
595 std::invalid_argument, errmsg);
596 }
597#endif
598
599ArrayWrapper<Scalar,ArrayOutData, Rank<ArrayOutData >::value, false>outputDataWrap(outputData);
600ArrayWrapper<Scalar,ArrayInDataLeft, Rank<ArrayInDataLeft >::value, true>inputDataLeftWrap(inputDataLeft);
601ArrayWrapper<Scalar,ArrayInDataRight, Rank<ArrayInDataRight >::value, true>inputDataRightWrap(inputDataRight);
602
603 // inputDataRight is (C,P,D)
604 if(getrank(inputDataRight) == 3){
605
606 for(size_t cell = 0; cell <(size_t) inputDataLeft.dimension(0); cell++){
607 for(size_t point = 0; point <(size_t) inputDataLeft.dimension(1); point++){
608 for(size_t row = 0; row <(size_t) inputDataLeft.dimension(2); row++){
609 for(size_t col = 0; col <(size_t) inputDataLeft.dimension(2); col++){
610
611 outputDataWrap(cell, point, row, col) = \
612 inputDataLeftWrap(cell, point, row)*inputDataRightWrap(cell, point, col);
613 }// col
614 }// row
615 }// point
616 } // cell
617 }// rank = 3
618 // inputDataRight is (P,D)
619 else if(getrank(inputDataRight) == 2){
620
621 for(size_t cell = 0; cell <(size_t) inputDataLeft.dimension(0); cell++){
622 for(size_t point = 0; point <(size_t) inputDataLeft.dimension(1); point++){
623 for(size_t row = 0; row <(size_t) inputDataLeft.dimension(2); row++){
624 for(size_t col = 0; col <(size_t) inputDataLeft.dimension(2); col++){
625 //
626 outputDataWrap(cell, point, row, col) = \
627 inputDataLeftWrap(cell, point, row)*inputDataRightWrap(point, col);
628 } // col
629 } // row
630 } // point
631 } // cell
632 }// rank = 2
633 else{
634 TEUCHOS_TEST_FOR_EXCEPTION( true, std::invalid_argument,
635 ">>> ERROR (ArrayTools::crossProductDataData): inputDataRight rank 2 or 3 required.")
636 }
637 }
638 template<class Scalar,
639 class ArrayOutFields,
640 class ArrayInData,
641 class ArrayInFields>
642 void ArrayTools::matvecProductDataField(ArrayOutFields & outputFields,
643 const ArrayInData & inputData,
644 const ArrayInFields & inputFields,
645 const char transpose){
646
647#ifdef HAVE_INTREPID_DEBUG
648 std::string errmsg = ">>> ERROR (ArrayTools::matvecProductDataField):";
649 /*
650 * Check array rank and spatial dimension range (if applicable)
651 * (1) inputData(C,P), (C,P,D) or (C,P,D,D); P=1 is admissible to allow multiply by const. data
652 * (2) inputFields(C,F,P,D) or (F,P,D);
653 * (3) outputFields(C,F,P,D)
654 */
655 // (1) inputData is (C,P), (C, P, D) or (C, P, D, D) and 1 <= D <= 3 is required
656 TEUCHOS_TEST_FOR_EXCEPTION( !requireRankRange(errmsg, inputData, 2,4),
657 std::invalid_argument, errmsg);
658 if(getrank(inputData) > 2) {
659 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, inputData, 2, 1,3),
660 std::invalid_argument, errmsg);
661 }
662 if(getrank(inputData) == 4) {
663 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, inputData, 3, 1,3),
664 std::invalid_argument, errmsg);
665 }
666 // (2) inputFields is (C, F, P, D) or (F, P, D) and 1 <= (D=dimension(rank - 1)) <= 3 is required.
667 TEUCHOS_TEST_FOR_EXCEPTION( !requireRankRange(errmsg, inputFields, 3,4), std::invalid_argument, errmsg);
668 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, inputFields, getrank(inputFields)-1, 1,3),
669 std::invalid_argument, errmsg);
670 // (3) outputFields is (C,F,P,D)
671 TEUCHOS_TEST_FOR_EXCEPTION( !requireRankRange(errmsg, outputFields, 4,4), std::invalid_argument, errmsg);
672 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, outputFields, 3, 1,3),
673 std::invalid_argument, errmsg);
674 /*
675 * Dimension cross-checks:
676 * (1) inputData vs. inputFields
677 * (2) outputFields vs. inputData
678 * (3) outputFields vs. inputFields
679 *
680 * Cross-check (2): outputFields(C,F,P,D) vs. inputData(C,P), (C,P,D) or (C,P,D,D):
681 * dimensions C, and D must match in all cases, dimension P must match only when non-constant
682 * data is specified (P>1). Do not check P dimensions with constant data, i.e., when P=1 in
683 * inputData(C,1,...)
684 */
685 if(inputData.dimension(1) > 1){ // check P dimension if P>1 in inputData
686 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg,
687 outputFields, 2,
688 inputData, 1),
689 std::invalid_argument, errmsg);
690 }
691 if(getrank(inputData) == 2) { // inputData(C,P) -> C match
692 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg,
693 outputFields, 0,
694 inputData, 0),
695 std::invalid_argument, errmsg);
696 }
697 if(getrank(inputData) == 3){ // inputData(C,P,D) -> C, D match
698 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg,
699 outputFields, 0,3,
700 inputData, 0,2),
701 std::invalid_argument, errmsg);
702
703 }
704 if(getrank(inputData) == 4){ // inputData(C,P,D,D) -> C, D, D match
705 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg,
706 outputFields, 0,3,3,
707 inputData, 0,2,3),
708 std::invalid_argument, errmsg);
709 }
710 /*
711 * Cross-checks (1,3):
712 */
713 if(getrank(inputFields) == 4) {
714 // 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
715 if(inputData.dimension(1) > 1){ // check P dimension if P>1 in inputData
716 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg,
717 inputFields, 2,
718 inputData, 1),
719 std::invalid_argument, errmsg);
720 }
721 if(getrank(inputData) == 2){ // inputData(C,P) -> C match
722 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg,
723 inputData, 0,
724 inputFields, 0),
725 std::invalid_argument, errmsg);
726 }
727 if(getrank(inputData) == 3){ // inputData(C,P,D) -> C, D match
728 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg,
729 inputData, 0,2,
730 inputFields, 0,3),
731 std::invalid_argument, errmsg);
732 }
733 if(getrank(inputData) == 4){ // inputData(C,P,D,D) -> C, D, D match
734 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg,
735 inputData, 0,2,3,
736 inputFields, 0,3,3),
737 std::invalid_argument, errmsg);
738 }
739 // Cross-check (3): outputFields(C,F,P,D) vs. inputFields(C,F,P,D): all dimensions C, F, P, D must match
740 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, outputFields, inputFields),
741 std::invalid_argument, errmsg);
742 }
743 else{
744 // Cross-check (1): inputData(C,P), (C,P,D) or (C,P,D,D) vs. inputFields(F,P,D): dimensions P, D must match
745 if(inputData.dimension(1) > 1){ // check P if P>1 in inputData
746 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg,
747 inputData, 1,
748 inputFields, 1),
749 std::invalid_argument, errmsg);
750 }
751 if(getrank(inputData) == 3){
752 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg,
753 inputData, 2,
754 inputFields, 2),
755 std::invalid_argument, errmsg);
756 }
757 if(getrank(inputData) == 4){
758 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg,
759 inputData, 2,3,
760 inputFields, 2,2),
761 std::invalid_argument, errmsg);
762 }
763 // Cross-check (3): outputFields(C,F,P,D) vs. inputFields(F,P,D): dimensions F, P, D must match
764 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg,
765 outputFields, 1,2,3,
766 inputFields, 0,1,2),
767 std::invalid_argument, errmsg);
768 }
769#endif
770
771 ArrayWrapper<Scalar,ArrayOutFields, Rank<ArrayOutFields >::value, false>outputFieldsWrap(outputFields);
772 ArrayWrapper<Scalar,ArrayInData, Rank<ArrayInData >::value, true>inputDataWrap(inputData);
773 ArrayWrapper<Scalar,ArrayInFields, Rank<ArrayInFields >::value,true>inputFieldsWrap(inputFields);
774 size_t dataRank = getrank(inputData);
775 size_t numDataPts = inputData.dimension(1);
776 size_t inRank = getrank(inputFields);
777 size_t numCells = outputFields.dimension(0);
778 size_t numFields = outputFields.dimension(1);
779 size_t numPoints = outputFields.dimension(2);
780 size_t matDim = outputFields.dimension(3);
781 /*********************************************************************************************
782 * inputFields is (C,F,P,D) *
783 *********************************************************************************************/
784 if(inRank == 4){
785 if(numDataPts != 1){ // non-constant data
786
787 switch(dataRank){
788 case 2:
789 for(size_t cell = 0; cell < numCells; cell++) {
790 for(size_t field = 0; field < numFields; field++) {
791 for(size_t point = 0; point < numPoints; point++) {
792 for( size_t row = 0; row < matDim; row++) {
793 outputFieldsWrap(cell, field, point, row) = \
794 inputDataWrap(cell, point)*inputFieldsWrap(cell, field, point, row);
795 } // Row-loop
796 } // P-loop
797 } // F-loop
798 }// C-loop
799 break;
800
801 case 3:
802 for(size_t cell = 0; cell < numCells; cell++) {
803 for(size_t field = 0; field < numFields; field++) {
804 for(size_t point = 0; point < numPoints; point++) {
805 for(size_t row = 0; row < matDim; row++) {
806 outputFieldsWrap(cell, field, point, row) = \
807 inputDataWrap(cell, point, row)*inputFieldsWrap(cell, field, point, row);
808 } // Row-loop
809 } // P-loop
810 } // F-loop
811 }// C-loop
812 break;
813
814 case 4:
815 if ((transpose == 'n') || (transpose == 'N')) {
816 for(size_t cell = 0; cell < numCells; cell++){
817 for(size_t field = 0; field < numFields; field++){
818 for(size_t point = 0; point < numPoints; point++){
819 for(size_t row = 0; row < matDim; row++){
820 outputFieldsWrap(cell, field, point, row) = 0.0;
821 for(size_t col = 0; col < matDim; col++){
822 outputFieldsWrap(cell, field, point, row) += \
823 inputDataWrap(cell, point, row, col)*inputFieldsWrap(cell, field, point, col);
824 }// col
825 } //row
826 }// point
827 }// field
828 }// cell
829 } // no transpose
830 else if ((transpose == 't') || (transpose == 'T')) {
831 for(size_t cell = 0; cell < numCells; cell++){
832 for(size_t field = 0; field < numFields; field++){
833 for(size_t point = 0; point < numPoints; point++){
834 for(size_t row = 0; row < matDim; row++){
835 outputFieldsWrap(cell, field, point, row) = 0.0;
836 for(size_t col = 0; col < matDim; col++){
837 outputFieldsWrap(cell, field, point, row) += \
838 inputDataWrap(cell, point, col, row)*inputFieldsWrap(cell, field, point, col);
839 }// col
840 } //row
841 }// point
842 }// field
843 }// cell
844 } //transpose
845 else {
846 TEUCHOS_TEST_FOR_EXCEPTION( !( (transpose == 'n') || (transpose == 'N') || (transpose == 't') || (transpose == 'T') ), std::invalid_argument,
847 ">>> ERROR (ArrayTools::matvecProductDataField): The transpose flag must be 'n', 'N', 't' or 'T'.");
848 }
849 break;
850
851 default:
852 TEUCHOS_TEST_FOR_EXCEPTION( !( (dataRank == 2) || (dataRank == 3) || (dataRank == 4) ), std::invalid_argument,
853 ">>> ERROR (ArrayTools::matvecProductDataField): inputData rank 2, 3 or 4 required.")
854 } // switch inputData rank
855 }
856 else{ // constant data case
857 switch(dataRank){
858 case 2:
859 for(size_t cell = 0; cell < numCells; cell++) {
860 for(size_t field = 0; field < numFields; field++) {
861 for(size_t point = 0; point < numPoints; point++) {
862 for(size_t row = 0; row < matDim; row++) {
863 outputFieldsWrap(cell, field, point, row) = \
864 inputDataWrap(cell, 0)*inputFieldsWrap(cell, field, point, row);
865 } // Row-loop
866 } // P-loop
867 } // F-loop
868 }// C-loop
869 break;
870
871 case 3:
872 for(size_t cell = 0; cell < numCells; cell++) {
873 for(size_t field = 0; field < numFields; field++) {
874 for(size_t point = 0; point < numPoints; point++) {
875 for(size_t row = 0; row < matDim; row++) {
876 outputFieldsWrap(cell, field, point, row) = \
877 inputDataWrap(cell, 0, row)*inputFieldsWrap(cell, field, point, row);
878 } // Row-loop
879 } // P-loop
880 } // F-loop
881 }// C-loop
882 break;
883
884 case 4:
885 if ((transpose == 'n') || (transpose == 'N')) {
886 for(size_t cell = 0; cell < numCells; cell++){
887 for(size_t field = 0; field < numFields; field++){
888 for(size_t point = 0; point < numPoints; point++){
889 for(size_t row = 0; row < matDim; row++){
890 outputFieldsWrap(cell, field, point, row) = 0.0;
891 for(size_t col = 0; col < matDim; col++){
892 outputFieldsWrap(cell, field, point, row) += \
893 inputDataWrap(cell, 0, row, col)*inputFieldsWrap(cell, field, point, col);
894 }// col
895 } //row
896 }// point
897 }// field
898 }// cell
899 } // no transpose
900 else if ((transpose == 't') || (transpose == 'T')) {
901 for(size_t cell = 0; cell < numCells; cell++){
902 for(size_t field = 0; field < numFields; field++){
903 for(size_t point = 0; point < numPoints; point++){
904 for(size_t row = 0; row < matDim; row++){
905 outputFieldsWrap(cell, field, point, row) = 0.0;
906 for(size_t col = 0; col < matDim; col++){
907 outputFieldsWrap(cell, field, point, row) += \
908 inputDataWrap(cell, 0, col, row)*inputFieldsWrap(cell, field, point, col);
909 }// col
910 } //row
911 }// point
912 }// field
913 }// cell
914 } //transpose
915 else {
916 TEUCHOS_TEST_FOR_EXCEPTION( !( (transpose == 'n') || (transpose == 'N') || (transpose == 't') || (transpose == 'T') ), std::invalid_argument,
917 ">>> ERROR (ArrayTools::matvecProductDataField): The transpose flag must be 'n', 'N', 't' or 'T'.");
918 }
919 break;
920
921 default:
922 TEUCHOS_TEST_FOR_EXCEPTION( !( (dataRank == 2) || (dataRank == 3) || (dataRank == 4) ), std::invalid_argument,
923 ">>> ERROR (ArrayTools::matvecProductDataField): inputData rank 2, 3 or 4 required.")
924 } // switch inputData rank
925 } // end constant data case
926 } // inputFields rank 4
927 /*********************************************************************************************
928 * inputFields is (F,P,D) *
929 *********************************************************************************************/
930 else if(inRank == 3) {
931 if(numDataPts != 1){ // non-constant data
932
933 switch(dataRank){
934 case 2:
935 for(size_t cell = 0; cell < numCells; cell++) {
936 for(size_t field = 0; field < numFields; field++) {
937 for(size_t point = 0; point < numPoints; point++) {
938 for(size_t row = 0; row < matDim; row++) {
939 outputFieldsWrap(cell, field, point, row) = \
940 inputDataWrap(cell, point)*inputFieldsWrap(field, point, row);
941 } // Row-loop
942 } // P-loop
943 } // F-loop
944 }// C-loop
945 break;
946
947 case 3:
948 for(size_t cell = 0; cell < numCells; cell++) {
949 for(size_t field = 0; field < numFields; field++) {
950 for(size_t point = 0; point < numPoints; point++) {
951 for(size_t row = 0; row < matDim; row++) {
952 outputFieldsWrap(cell, field, point, row) = \
953 inputDataWrap(cell, point, row)*inputFieldsWrap(field, point, row);
954 } // Row-loop
955 } // P-loop
956 } // F-loop
957 }// C-loop
958 break;
959
960 case 4:
961 if ((transpose == 'n') || (transpose == 'N')) {
962 for(size_t cell = 0; cell < numCells; cell++){
963 for(size_t field = 0; field < numFields; field++){
964 for(size_t point = 0; point < numPoints; point++){
965 for(size_t row = 0; row < matDim; row++){
966 outputFieldsWrap(cell, field, point, row) = 0.0;
967 for(size_t col = 0; col < matDim; col++){
968 outputFieldsWrap(cell, field, point, row) += \
969 inputDataWrap(cell, point, row, col)*inputFieldsWrap(field, point, col);
970 }// col
971 } //row
972 }// point
973 }// field
974 }// cell
975 } // no transpose
976 else if ((transpose == 't') || (transpose == 'T')) {
977 for(size_t cell = 0; cell < numCells; cell++){
978 for(size_t field = 0; field < numFields; field++){
979 for(size_t point = 0; point < numPoints; point++){
980 for(size_t row = 0; row < matDim; row++){
981 outputFieldsWrap(cell, field, point, row) = 0.0;
982 for(size_t col = 0; col < matDim; col++){
983 outputFieldsWrap(cell, field, point, row) += \
984 inputDataWrap(cell, point, col, row)*inputFieldsWrap(field, point, col);
985 }// col
986 } //row
987 }// point
988 }// field
989 }// cell
990 } //transpose
991 else {
992 TEUCHOS_TEST_FOR_EXCEPTION( !( (transpose == 'n') || (transpose == 'N') || (transpose == 't') || (transpose == 'T') ), std::invalid_argument,
993 ">>> ERROR (ArrayTools::matvecProductDataField): The transpose flag must be 'n', 'N', 't' or 'T'.");
994 }
995 break;
996
997 default:
998 TEUCHOS_TEST_FOR_EXCEPTION( !( (dataRank == 2) || (dataRank == 3) || (dataRank == 4) ), std::invalid_argument,
999 ">>> ERROR (ArrayTools::matvecProductDataField): inputData rank 2, 3 or 4 required.")
1000 } // switch inputData rank
1001 }
1002 else{ // constant data case
1003 switch(dataRank){
1004 case 2:
1005 for(size_t cell = 0; cell < numCells; cell++) {
1006 for(size_t field = 0; field < numFields; field++) {
1007 for(size_t point = 0; point < numPoints; point++) {
1008 for(size_t row = 0; row < matDim; row++) {
1009 outputFieldsWrap(cell, field, point, row) = \
1010 inputDataWrap(cell, 0)*inputFieldsWrap(field, point, row);
1011 } // Row-loop
1012 } // P-loop
1013 } // F-loop
1014 }// C-loop
1015 break;
1016
1017 case 3:
1018 for(size_t cell = 0; cell < numCells; cell++) {
1019 for(size_t field = 0; field < numFields; field++) {
1020 for(size_t point = 0; point < numPoints; point++) {
1021 for(size_t row = 0; row < matDim; row++) {
1022 outputFieldsWrap(cell, field, point, row) = \
1023 inputDataWrap(cell, 0, row)*inputFieldsWrap(field, point, row);
1024 } // Row-loop
1025 } // P-loop
1026 } // F-loop
1027 }// C-loop
1028 break;
1029
1030 case 4:
1031 if ((transpose == 'n') || (transpose == 'N')) {
1032 for(size_t cell = 0; cell < numCells; cell++){
1033 for(size_t field = 0; field < numFields; field++){
1034 for(size_t point = 0; point < numPoints; point++){
1035 for(size_t row = 0; row < matDim; row++){
1036 outputFieldsWrap(cell, field, point, row) = 0.0;
1037 for(size_t col = 0; col < matDim; col++){
1038 outputFieldsWrap(cell, field, point, row) += \
1039 inputDataWrap(cell, 0, row, col)*inputFieldsWrap(field, point, col);
1040 }// col
1041 } //row
1042 }// point
1043 }// field
1044 }// cell
1045 } // no transpose
1046 else if ((transpose == 't') || (transpose == 'T')) {
1047 for(size_t cell = 0; cell < numCells; cell++){
1048 for(size_t field = 0; field < numFields; field++){
1049 for(size_t point = 0; point < numPoints; point++){
1050 for(size_t row = 0; row < matDim; row++){
1051 outputFieldsWrap(cell, field, point, row) = 0.0;
1052 for(size_t col = 0; col < matDim; col++){
1053 outputFieldsWrap(cell, field, point, row) += \
1054 inputDataWrap(cell, 0, col, row)*inputFieldsWrap(field, point, col);
1055 }// col
1056 } //row
1057 }// point
1058 }// field
1059 }// cell
1060 } //transpose
1061 else {
1062 TEUCHOS_TEST_FOR_EXCEPTION( !( (transpose == 'n') || (transpose == 'N') || (transpose == 't') || (transpose == 'T') ), std::invalid_argument,
1063 ">>> ERROR (ArrayTools::matvecProductDataField): The transpose flag must be 'n', 'N', 't' or 'T'.");
1064 }
1065 break;
1066
1067 default:
1068 TEUCHOS_TEST_FOR_EXCEPTION( !( (dataRank == 2) || (dataRank == 3) || (dataRank == 4) ), std::invalid_argument,
1069 ">>> ERROR (ArrayTools::matvecProductDataField): inputData rank 2, 3 or 4 required.")
1070 } // switch inputData rank
1071 } // end constant data case
1072 } // inputFields rank 3
1073 else {
1074 TEUCHOS_TEST_FOR_EXCEPTION( true, std::invalid_argument,
1075 ">>> ERROR (ArrayTools::matvecProductDataField): inputFields rank 3 or 4 required.")
1076 }// rank error
1077 }
1078
1079
1080
1081 template<class Scalar,
1082 class ArrayOutData,
1083 class ArrayInDataLeft,
1084 class ArrayInDataRight>
1085 void ArrayTools::matvecProductDataData(ArrayOutData & outputData,
1086 const ArrayInDataLeft & inputDataLeft,
1087 const ArrayInDataRight & inputDataRight,
1088 const char transpose){
1089
1090#ifdef HAVE_INTREPID_DEBUG
1091 std::string errmsg = ">>> ERROR (ArrayTools::matvecProductDataData):";
1092 /*
1093 * Check array rank and spatial dimension range (if applicable)
1094 * (1) inputDataLeft(C,P), (C,P,D) or (C,P,D,D); P=1 is admissible to allow multiply by const. left data
1095 * (2) inputDataRight(C,P,D) or (P,D);
1096 * (3) outputData(C,P,D)
1097 */
1098 // (1) inputDataLeft is (C,P), (C,P,D) or (C,P,D,D) and 1 <= D <= 3 is required
1099 TEUCHOS_TEST_FOR_EXCEPTION( !requireRankRange(errmsg, inputDataLeft, 2,4),
1100 std::invalid_argument, errmsg);
1101 if(getrank(inputDataLeft) > 2) {
1102 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, inputDataLeft, 2, 1,3),
1103 std::invalid_argument, errmsg);
1104 }
1105 if(getrank(inputDataLeft) == 4) {
1106 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, inputDataLeft, 3, 1,3),
1107 std::invalid_argument, errmsg);
1108 }
1109 // (2) inputDataRight is (C, P, D) or (P, D) and 1 <= (D=dimension(rank - 1)) <= 3 is required.
1110 TEUCHOS_TEST_FOR_EXCEPTION( !requireRankRange(errmsg, inputDataRight, 2,3),
1111 std::invalid_argument, errmsg);
1112 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, inputDataRight, getrank(inputDataRight)-1, 1,3),
1113 std::invalid_argument, errmsg);
1114 // (3) outputData is (C,P,D)
1115 TEUCHOS_TEST_FOR_EXCEPTION( !requireRankRange(errmsg, outputData, 3,3), std::invalid_argument, errmsg);
1116 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, outputData, 2, 1,3),
1117 std::invalid_argument, errmsg);
1118 /*
1119 * Dimension cross-checks:
1120 * (1) inputDataLeft vs. inputDataRight
1121 * (2) outputData vs. inputDataLeft
1122 * (3) outputData vs. inputDataRight
1123 *
1124 * Cross-check (2): outputData(C,P,D) vs. inputDataLeft(C,P), (C,P,D) or (C,P,D,D):
1125 * dimensions C, and D must match in all cases, dimension P must match only when non-constant
1126 * data is specified (P>1). Do not check P dimensions with constant left data, i.e., when P=1 in
1127 * inputDataLeft(C,1,...)
1128 */
1129 if(inputDataLeft.dimension(1) > 1){ // check P dimension if P>1 in inputDataLeft
1130 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg,
1131 outputData, 1,
1132 inputDataLeft, 1),
1133 std::invalid_argument, errmsg);
1134 }
1135 if(getrank(inputDataLeft) == 2){ // inputDataLeft(C,P): check C
1136 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg,
1137 outputData, 0,
1138 inputDataLeft, 0),
1139 std::invalid_argument, errmsg);
1140 }
1141 if(getrank(inputDataLeft) == 3){ // inputDataLeft(C,P,D): check C and D
1142 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg,
1143 outputData, 0,2,
1144 inputDataLeft, 0,2),
1145 std::invalid_argument, errmsg);
1146 }
1147 if(getrank(inputDataLeft) == 4){ // inputDataLeft(C,P,D,D): check C and D
1148 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg,
1149 outputData, 0,2,2,
1150 inputDataLeft, 0,2,3),
1151 std::invalid_argument, errmsg);
1152 }
1153 /*
1154 * Cross-checks (1,3):
1155 */
1156 if( getrank(inputDataRight) == 3) {
1157 // 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
1158 if(inputDataLeft.dimension(1) > 1){ // check P dimension if P>1 in inputDataLeft
1159 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg,
1160 inputDataLeft, 1,
1161 inputDataRight, 1),
1162 std::invalid_argument, errmsg);
1163 }
1164 if(getrank(inputDataLeft) == 2){ // inputDataLeft(C,P): check C
1165 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg,
1166 inputDataLeft, 0,
1167 inputDataRight, 0),
1168 std::invalid_argument, errmsg);
1169 }
1170 if(getrank(inputDataLeft) == 3){ // inputDataLeft(C,P,D): check C and D
1171 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg,
1172 inputDataLeft, 0,2,
1173 inputDataRight, 0,2),
1174 std::invalid_argument, errmsg);
1175 }
1176 if(getrank(inputDataLeft) == 4){ // inputDataLeft(C,P,D,D): check C and D
1177 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg,
1178 inputDataLeft, 0,2,3,
1179 inputDataRight, 0,2,2),
1180 std::invalid_argument, errmsg);
1181 }
1182
1183 // Cross-check (3): outputData(C,P,D) vs. inputDataRight(C,P,D): all dimensions C, P, D must match
1184 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, outputData, inputDataRight),
1185 std::invalid_argument, errmsg);
1186 }
1187 else{
1188 // Cross-check (1): inputDataLeft(C,P), (C,P,D), or (C,P,D,D) vs. inputDataRight(P,D): dimensions P, D must match
1189 if(inputDataLeft.dimension(1) > 1){ // check P if P>1 in inputData
1190 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg,
1191 inputDataLeft, 1,
1192 inputDataRight, 0),
1193 std::invalid_argument, errmsg);
1194 }
1195 if(getrank(inputDataLeft) == 3){ // inputDataLeft(C,P,D): check D
1196 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg,
1197 inputDataLeft, 2,
1198 inputDataRight, 1),
1199 std::invalid_argument, errmsg);
1200 }
1201 if(getrank(inputDataLeft) == 4){ // inputDataLeft(C,P,D,D): check D
1202 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg,
1203 inputDataLeft, 2,3,
1204 inputDataRight, 1,1),
1205 std::invalid_argument, errmsg);
1206 }
1207 // Cross-check (3): outputData(C,P,D) vs. inputDataRight(P,D): dimensions P, D must match
1208 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg,
1209 outputData, 1,2,
1210 inputDataRight, 0,1),
1211 std::invalid_argument, errmsg);
1212 }
1213#endif
1214
1215 ArrayWrapper<Scalar,ArrayOutData, Rank<ArrayOutData >::value, false>outputDataWrap(outputData);
1216 ArrayWrapper<Scalar,ArrayInDataLeft, Rank<ArrayInDataLeft >::value, true>inputDataLeftWrap(inputDataLeft);
1217 ArrayWrapper<Scalar,ArrayInDataRight, Rank<ArrayInDataRight >::value, true>inputDataRightWrap(inputDataRight);
1218 size_t dataLeftRank = getrank(inputDataLeft);
1219 size_t numDataLeftPts = inputDataLeft.dimension(1);
1220 size_t dataRightRank = getrank(inputDataRight);
1221 size_t numCells = outputData.dimension(0);
1222 size_t numPoints = outputData.dimension(1);
1223 size_t matDim = outputData.dimension(2);
1224
1225 /*********************************************************************************************
1226 * inputDataRight is (C,P,D) *
1227 *********************************************************************************************/
1228 if(dataRightRank == 3){
1229 if(numDataLeftPts != 1){ // non-constant left data
1230
1231 switch(dataLeftRank){
1232 case 2:
1233 for(size_t cell = 0; cell < numCells; cell++) {
1234 for(size_t point = 0; point < numPoints; point++) {
1235 for(size_t row = 0; row < matDim; row++) {
1236 outputDataWrap(cell, point, row) = \
1237 inputDataLeftWrap(cell, point)*inputDataRightWrap(cell, point, row);
1238 } // Row-loop
1239 } // P-loop
1240 }// C-loop
1241 break;
1242
1243 case 3:
1244 for(size_t cell = 0; cell < numCells; cell++) {
1245 for(size_t point = 0; point < numPoints; point++) {
1246 for(size_t row = 0; row < matDim; row++) {
1247 outputDataWrap(cell, point, row) = \
1248 inputDataLeftWrap(cell, point, row)*inputDataRightWrap(cell, point, row);
1249 } // Row-loop
1250 } // P-loop
1251 }// C-loop
1252 break;
1253
1254 case 4:
1255 if ((transpose == 'n') || (transpose == 'N')) {
1256 for(size_t cell = 0; cell < numCells; cell++){
1257 for(size_t point = 0; point < numPoints; point++){
1258 for(size_t row = 0; row < matDim; row++){
1259 outputDataWrap(cell, point, row) = 0.0;
1260 for(size_t col = 0; col < matDim; col++){
1261 outputDataWrap(cell, point, row) += \
1262 inputDataLeftWrap(cell, point, row, col)*inputDataRightWrap(cell, point, col);
1263 }// col
1264 } //row
1265 }// point
1266 }// cell
1267 } // no transpose
1268 else if ((transpose == 't') || (transpose == 'T')) {
1269 for(size_t cell = 0; cell < numCells; cell++){
1270 for(size_t point = 0; point < numPoints; point++){
1271 for(size_t row = 0; row < matDim; row++){
1272 outputDataWrap(cell, point, row) = 0.0;
1273 for(size_t col = 0; col < matDim; col++){
1274 outputDataWrap(cell, point, row) += \
1275 inputDataLeftWrap(cell, point, col, row)*inputDataRightWrap(cell, point, col);
1276 }// col
1277 } //row
1278 }// point
1279 }// cell
1280 } //transpose
1281 else {
1282 TEUCHOS_TEST_FOR_EXCEPTION( !( (transpose == 'n') || (transpose == 'N') || (transpose == 't') || (transpose == 'T') ), std::invalid_argument,
1283 ">>> ERROR (ArrayTools::matvecProductDataData): The transpose flag must be 'n', 'N', 't' or 'T'.");
1284 }
1285 break;
1286
1287 default:
1288 TEUCHOS_TEST_FOR_EXCEPTION( !( (dataLeftRank == 2) || (dataLeftRank == 3) || (dataLeftRank == 4) ), std::invalid_argument,
1289 ">>> ERROR (ArrayTools::matvecProductDataData): inputDataLeft rank 2, 3 or 4 required.")
1290 } // switch inputDataLeft rank
1291 }
1292 else{ // constant data case
1293 switch(dataLeftRank){
1294 case 2:
1295 for(size_t cell = 0; cell < numCells; cell++) {
1296 for(size_t point = 0; point < numPoints; point++) {
1297 for(size_t row = 0; row < matDim; row++) {
1298 outputDataWrap(cell, point, row) = \
1299 inputDataLeftWrap(cell, 0)*inputDataRightWrap(cell, point, row);
1300 } // Row-loop
1301 } // F-loop
1302 }// C-loop
1303 break;
1304
1305 case 3:
1306 for(size_t cell = 0; cell < numCells; cell++) {
1307 for(size_t point = 0; point < numPoints; point++) {
1308 for(size_t row = 0; row < matDim; row++) {
1309 outputDataWrap(cell, point, row) = \
1310 inputDataLeftWrap(cell, 0, row)*inputDataRightWrap(cell, point, row);
1311 } // Row-loop
1312 } // P-loop
1313 }// C-loop
1314 break;
1315
1316 case 4:
1317 if ((transpose == 'n') || (transpose == 'N')) {
1318 for(size_t cell = 0; cell < numCells; cell++){
1319 for(size_t point = 0; point < numPoints; point++){
1320 for(size_t row = 0; row < matDim; row++){
1321 outputDataWrap(cell, point, row) = 0.0;
1322 for(size_t col = 0; col < matDim; col++){
1323 outputDataWrap(cell, point, row) += \
1324 inputDataLeftWrap(cell, 0, row, col)*inputDataRightWrap(cell, point, col);
1325 }// col
1326 } //row
1327 }// point
1328 }// cell
1329 } // no transpose
1330 else if ((transpose == 't') || (transpose == 'T')) {
1331 for(size_t cell = 0; cell < numCells; cell++){
1332 for(size_t point = 0; point < numPoints; point++){
1333 for(size_t row = 0; row < matDim; row++){
1334 outputDataWrap(cell, point, row) = 0.0;
1335 for(size_t col = 0; col < matDim; col++){
1336 outputDataWrap(cell, point, row) += \
1337 inputDataLeftWrap(cell, 0, col, row)*inputDataRightWrap(cell, point, col);
1338 }// col
1339 } //row
1340 }// point
1341 }// cell
1342 } //transpose
1343 else {
1344 TEUCHOS_TEST_FOR_EXCEPTION( !( (transpose == 'n') || (transpose == 'N') || (transpose == 't') || (transpose == 'T') ), std::invalid_argument,
1345 ">>> ERROR (ArrayTools::matvecProductDataData): The transpose flag must be 'n', 'N', 't' or 'T'.");
1346 }
1347 break;
1348
1349 default:
1350 TEUCHOS_TEST_FOR_EXCEPTION( !( (dataLeftRank == 2) || (dataLeftRank == 3) || (dataLeftRank == 4) ), std::invalid_argument,
1351 ">>> ERROR (ArrayTools::matvecProductDataData): inputDataLeft rank 2, 3 or 4 required.")
1352 } // switch inputDataLeft rank
1353 } // end constant data case
1354 } // inputDataRight rank 4
1355 /*********************************************************************************************
1356 * inputDataRight is (P,D) *
1357 *********************************************************************************************/
1358 else if(dataRightRank == 2) {
1359 if(numDataLeftPts != 1){ // non-constant data
1360
1361 switch(dataLeftRank){
1362 case 2:
1363 for(size_t cell = 0; cell < numCells; cell++) {
1364 for(size_t point = 0; point < numPoints; point++) {
1365 for(size_t row = 0; row < matDim; row++) {
1366 outputDataWrap(cell, point, row) = \
1367 inputDataLeftWrap(cell, point)*inputDataRightWrap(point, row);
1368 } // Row-loop
1369 } // P-loop
1370 }// C-loop
1371 break;
1372
1373 case 3:
1374 for(size_t cell = 0; cell < numCells; cell++) {
1375 for(size_t point = 0; point < numPoints; point++) {
1376 for(size_t row = 0; row < matDim; row++) {
1377 outputDataWrap(cell, point, row) = \
1378 inputDataLeftWrap(cell, point, row)*inputDataRightWrap(point, row);
1379 } // Row-loop
1380 } // P-loop
1381 }// C-loop
1382 break;
1383
1384 case 4:
1385 if ((transpose == 'n') || (transpose == 'N')) {
1386 for(size_t cell = 0; cell < numCells; cell++){
1387 for(size_t point = 0; point < numPoints; point++){
1388 for(size_t row = 0; row < matDim; row++){
1389 outputDataWrap(cell, point, row) = 0.0;
1390 for(size_t col = 0; col < matDim; col++){
1391 outputDataWrap(cell, point, row) += \
1392 inputDataLeftWrap(cell, point, row, col)*inputDataRightWrap(point, col);
1393 }// col
1394 } //row
1395 }// point
1396 }// cell
1397 } // no transpose
1398 else if ((transpose == 't') || (transpose == 'T')) {
1399 for(size_t cell = 0; cell < numCells; cell++){
1400 for(size_t point = 0; point < numPoints; point++){
1401 for(size_t row = 0; row < matDim; row++){
1402 outputDataWrap(cell, point, row) = 0.0;
1403 for(size_t col = 0; col < matDim; col++){
1404 outputDataWrap(cell, point, row) += \
1405 inputDataLeftWrap(cell, point, col, row)*inputDataRightWrap(point, col);
1406 }// col
1407 } //row
1408 }// point
1409 }// cell
1410 } //transpose
1411 else {
1412 TEUCHOS_TEST_FOR_EXCEPTION( !( (transpose == 'n') || (transpose == 'N') || (transpose == 't') || (transpose == 'T') ), std::invalid_argument,
1413 ">>> ERROR (ArrayTools::matvecProductDataData): The transpose flag must be 'n', 'N', 't' or 'T'.");
1414 }
1415 break;
1416
1417 default:
1418 TEUCHOS_TEST_FOR_EXCEPTION( !( (dataLeftRank == 2) || (dataLeftRank == 3) || (dataLeftRank == 4) ), std::invalid_argument,
1419 ">>> ERROR (ArrayTools::matvecProductDataData): inputDataLeft rank 2, 3 or 4 required.")
1420 } // switch inputDataLeft rank
1421 }
1422 else{ // constant data case
1423 switch(dataLeftRank){
1424 case 2:
1425 for(size_t cell = 0; cell < numCells; cell++) {
1426 for(size_t point = 0; point < numPoints; point++) {
1427 for(size_t row = 0; row < matDim; row++) {
1428 outputDataWrap(cell, point, row) = \
1429 inputDataLeftWrap(cell, 0)*inputDataRightWrap(point, row);
1430 } // Row-loop
1431 } // P-loop
1432 }// C-loop
1433 break;
1434
1435 case 3:
1436 for(size_t cell = 0; cell < numCells; cell++) {
1437 for(size_t point = 0; point < numPoints; point++) {
1438 for(size_t row = 0; row < matDim; row++) {
1439 outputDataWrap(cell, point, row) = \
1440 inputDataLeftWrap(cell, 0, row)*inputDataRightWrap(point, row);
1441 } // Row-loop
1442 } // P-loop
1443 }// C-loop
1444 break;
1445
1446 case 4:
1447 if ((transpose == 'n') || (transpose == 'N')) {
1448 for(size_t cell = 0; cell < numCells; cell++){
1449 for(size_t point = 0; point < numPoints; point++){
1450 for(size_t row = 0; row < matDim; row++){
1451 outputDataWrap(cell, point, row) = 0.0;
1452 for(size_t col = 0; col < matDim; col++){
1453 outputDataWrap(cell, point, row) += \
1454 inputDataLeftWrap(cell, 0, row, col)*inputDataRightWrap(point, col);
1455 }// col
1456 } //row
1457 }// point
1458 }// cell
1459 } // no transpose
1460 else if ((transpose == 't') || (transpose == 'T')) {
1461 for(size_t cell = 0; cell < numCells; cell++){
1462 for(size_t point = 0; point < numPoints; point++){
1463 for(size_t row = 0; row < matDim; row++){
1464 outputDataWrap(cell, point, row) = 0.0;
1465 for(size_t col = 0; col < matDim; col++){
1466 outputDataWrap(cell, point, row) += \
1467 inputDataLeftWrap(cell, 0, col, row)*inputDataRightWrap(point, col);
1468 }// col
1469 } //row
1470 }// point
1471 }// cell
1472 } //transpose
1473 else {
1474 TEUCHOS_TEST_FOR_EXCEPTION( !( (transpose == 'n') || (transpose == 'N') || (transpose == 't') || (transpose == 'T') ), std::invalid_argument,
1475 ">>> ERROR (ArrayTools::matvecProductDataData): The transpose flag must be 'n', 'N', 't' or 'T'.");
1476 }
1477 break;
1478
1479 default:
1480 TEUCHOS_TEST_FOR_EXCEPTION( !( (dataLeftRank == 2) || (dataLeftRank == 3) || (dataLeftRank == 4) ), std::invalid_argument,
1481 ">>> ERROR (ArrayTools::matvecProductDataData): inputDataLeft rank 2, 3 or 4 required.")
1482 } // switch inputDataLeft rank
1483 } // end constant inputDataLeft case
1484 } // inputDataRight rank 2
1485 else {
1486 TEUCHOS_TEST_FOR_EXCEPTION( true, std::invalid_argument,
1487 ">>> ERROR (ArrayTools::matvecProductDataData): inputDataRight rank 2 or 3 required.")
1488 }// rank error
1489 }
1490
1491
1492
1493 template<class Scalar,
1494 class ArrayOutFields,
1495 class ArrayInData,
1496 class ArrayInFields>
1497 void ArrayTools::matmatProductDataField(ArrayOutFields & outputFields,
1498 const ArrayInData & inputData,
1499 const ArrayInFields & inputFields,
1500 const char transpose){
1501#ifdef HAVE_INTREPID_DEBUG
1502 std::string errmsg = ">>> ERROR (ArrayTools::matmatProductDataField):";
1503 /*
1504 * Check array rank and spatial dimension range (if applicable)
1505 * (1) inputData(C,P), (C,P,D) or (C,P,D,D); P=1 is admissible to allow multiply by const. data
1506 * (2) inputFields(C,F,P,D,D) or (F,P,D,D);
1507 * (3) outputFields(C,F,P,D,D)
1508 */
1509 // (1) inputData is (C,P), (C, P, D) or (C, P, D, D) and 1 <= D <= 3 is required
1510 TEUCHOS_TEST_FOR_EXCEPTION( !requireRankRange(errmsg, inputData, 2,4),
1511 std::invalid_argument, errmsg);
1512 if(getrank(inputData) > 2) {
1513 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, inputData, 2, 1,3),
1514 std::invalid_argument, errmsg);
1515 }
1516 if(getrank(inputData) == 4) {
1517 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, inputData, 3, 1,3),
1518 std::invalid_argument, errmsg);
1519 }
1520 // (2) inputFields is (C,F,P,D,D) or (F,P,D,D) and 1 <= (dimension(rank-1), (rank-2)) <= 3 is required.
1521 TEUCHOS_TEST_FOR_EXCEPTION( !requireRankRange(errmsg, inputFields, 4,5), std::invalid_argument, errmsg);
1522 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, inputFields, getrank(inputFields)-1, 1,3),
1523 std::invalid_argument, errmsg);
1524 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, inputFields, getrank(inputFields)-2, 1,3),
1525 std::invalid_argument, errmsg);
1526 // (3) outputFields is (C,F,P,D,D)
1527 TEUCHOS_TEST_FOR_EXCEPTION( !requireRankRange(errmsg, outputFields, 5,5), std::invalid_argument, errmsg);
1528 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, outputFields, 3, 1,3),
1529 std::invalid_argument, errmsg);
1530 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, outputFields, 4, 1,3),
1531 std::invalid_argument, errmsg);
1532 /*
1533 * Dimension cross-checks:
1534 * (1) inputData vs. inputFields
1535 * (2) outputFields vs. inputData
1536 * (3) outputFields vs. inputFields
1537 *
1538 * Cross-check (2): outputFields(C,F,P,D,D) vs. inputData(C,P), (C,P,D) or (C,P,D,D):
1539 * dimensions C, and D must match in all cases, dimension P must match only when non-constant
1540 * data is specified (P>1). Do not check P dimensions with constant data, i.e., when P=1 in
1541 * inputData(C,1,...)
1542 */
1543 if(inputData.dimension(1) > 1){ // check P dimension if P>1 in inputData
1544 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg,
1545 outputFields, 2,
1546 inputData, 1),
1547 std::invalid_argument, errmsg);
1548 }
1549 if(getrank(inputData) == 2) { // inputData(C,P) -> C match
1550 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg,
1551 outputFields, 0,
1552 inputData, 0),
1553 std::invalid_argument, errmsg);
1554 }
1555 if(getrank(inputData) == 3){ // inputData(C,P,D) -> C, D match
1556 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg,
1557 outputFields, 0,3,4,
1558 inputData, 0,2,2),
1559 std::invalid_argument, errmsg);
1560
1561 }
1562 if(getrank(inputData) == 4){ // inputData(C,P,D,D) -> C, D, D match
1563 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg,
1564 outputFields, 0,3,4,
1565 inputData, 0,2,3),
1566 std::invalid_argument, errmsg);
1567 }
1568 /*
1569 * Cross-checks (1,3):
1570 */
1571 if( getrank(inputFields) == 5) {
1572 // 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
1573 if(inputData.dimension(1) > 1){ // check P dimension if P>1 in inputData
1574 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg,
1575 inputData, 1,
1576 inputFields, 2),
1577 std::invalid_argument, errmsg);
1578 }
1579 if(getrank(inputData) == 2){ // inputData(C,P) -> C match
1580 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg,
1581 inputData, 0,
1582 inputFields, 0),
1583 std::invalid_argument, errmsg);
1584 }
1585 if(getrank(inputData) == 3){ // inputData(C,P,D) -> C, D match
1586 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg,
1587 inputData, 0,2,2,
1588 inputFields, 0,3,4),
1589 std::invalid_argument, errmsg);
1590 }
1591 if(getrank(inputData) == 4){ // inputData(C,P,D,D) -> C, D, D match
1592 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg,
1593 inputData, 0,2,3,
1594 inputFields, 0,3,4),
1595 std::invalid_argument, errmsg);
1596 }
1597 // Cross-check (3): outputFields(C,F,P,D,D) vs. inputFields(C,F,P,D,D): all dimensions C, F, P, D must match
1598 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, outputFields, inputFields),
1599 std::invalid_argument, errmsg);
1600 }
1601 else{
1602 // 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
1603 if(inputData.dimension(1) > 1){ // check P if P>1 in inputData
1604 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg,
1605 inputData, 1,
1606 inputFields, 1),
1607 std::invalid_argument, errmsg);
1608 }
1609 if(getrank(inputData) == 3){
1610 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg,
1611 inputData, 2,2,
1612 inputFields, 2,3),
1613 std::invalid_argument, errmsg);
1614 }
1615 if(getrank(inputData) == 4){
1616 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg,
1617 inputData, 2,3,
1618 inputFields, 2,3),
1619 std::invalid_argument, errmsg);
1620 }
1621 // Cross-check (3): outputFields(C,F,P,D,D) vs. inputFields(F,P,D,D): dimensions F, P, D must match
1622 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg,
1623 outputFields, 1,2,3,4,
1624 inputFields, 0,1,2,3),
1625 std::invalid_argument, errmsg);
1626 }
1627#endif
1628 ArrayWrapper<Scalar,ArrayOutFields, Rank<ArrayOutFields >::value, false>outputFieldsWrap(outputFields);
1629 ArrayWrapper<Scalar,ArrayInData, Rank<ArrayInData >::value, true>inputDataWrap(inputData);
1630 ArrayWrapper<Scalar,ArrayInFields, Rank<ArrayInFields >::value, true>inputFieldsWrap(inputFields);
1631
1632
1633 size_t dataRank = getrank(inputData);
1634 size_t numDataPts = inputData.dimension(1);
1635 size_t inRank = getrank(inputFields);
1636 size_t numCells = outputFields.dimension(0);
1637 size_t numFields = outputFields.dimension(1);
1638 size_t numPoints = outputFields.dimension(2);
1639 size_t matDim = outputFields.dimension(3);
1640
1641 /*********************************************************************************************
1642 * inputFields is (C,F,P,D,D) *
1643 *********************************************************************************************/
1644 if(inRank == 5){
1645 if(numDataPts != 1){ // non-constant data
1646
1647 switch(dataRank){
1648 case 2:
1649 for(size_t cell = 0; cell < numCells; cell++) {
1650 for(size_t field = 0; field < numFields; field++) {
1651 for(size_t point = 0; point < numPoints; point++) {
1652 for(size_t row = 0; row < matDim; row++) {
1653 for(size_t col = 0; col < matDim; col++) {
1654 outputFieldsWrap(cell, field, point, row, col) = \
1655 inputDataWrap(cell, point)*inputFieldsWrap(cell, field, point, row, col);
1656 }// Col-loop
1657 } // Row-loop
1658 } // P-loop
1659 } // F-loop
1660 }// C-loop
1661 break;
1662
1663 case 3:
1664 for(size_t cell = 0; cell < numCells; cell++) {
1665 for(size_t field = 0; field < numFields; field++) {
1666 for(size_t point = 0; point < numPoints; point++) {
1667 for(size_t row = 0; row < matDim; row++) {
1668 for(size_t col = 0; col < matDim; col++) {
1669 outputFieldsWrap(cell, field, point, row, col) = \
1670 inputDataWrap(cell, point, row)*inputFieldsWrap(cell, field, point, row, col);
1671 }// Col-loop
1672 } // Row-loop
1673 } // P-loop
1674 } // F-loop
1675 }// C-loop
1676 break;
1677
1678 case 4:
1679 if ((transpose == 'n') || (transpose == 'N')) {
1680 for(size_t cell = 0; cell < numCells; cell++){
1681 for(size_t field = 0; field < numFields; field++){
1682 for(size_t point = 0; point < numPoints; point++){
1683 for(size_t row = 0; row < matDim; row++){
1684 for(size_t col = 0; col < matDim; col++){
1685 outputFieldsWrap(cell, field, point, row, col) = 0.0;
1686 for(size_t i = 0; i < matDim; i++){
1687 outputFieldsWrap(cell, field, point, row, col) += \
1688 inputDataWrap(cell, point, row, i)*inputFieldsWrap(cell, field, point, i, col);
1689 }// i
1690 } // col
1691 } //row
1692 }// point
1693 }// field
1694 }// cell
1695 } // no transpose
1696 else if ((transpose == 't') || (transpose == 'T')) {
1697 for(size_t cell = 0; cell < numCells; cell++){
1698 for(size_t field = 0; field < numFields; field++){
1699 for(size_t point = 0; point < numPoints; point++){
1700 for(size_t row = 0; row < matDim; row++){
1701 for(size_t col = 0; col < matDim; col++){
1702 outputFieldsWrap(cell, field, point, row, col) = 0.0;
1703 for(size_t i = 0; i < matDim; i++){
1704 outputFieldsWrap(cell, field, point, row, col) += \
1705 inputDataWrap(cell, point, i, row)*inputFieldsWrap(cell, field, point, i, col);
1706 }// i
1707 } // col
1708 } //row
1709 }// point
1710 }// field
1711 }// cell
1712 } //transpose
1713 else {
1714 TEUCHOS_TEST_FOR_EXCEPTION( !( (transpose == 'n') || (transpose == 'N') || (transpose == 't') || (transpose == 'T') ), std::invalid_argument,
1715 ">>> ERROR (ArrayTools::matmatProductDataField): The transpose flag must be 'n', 'N', 't' or 'T'.");
1716 }
1717 break;
1718
1719 default:
1720 TEUCHOS_TEST_FOR_EXCEPTION( !( (dataRank == 2) || (dataRank == 3) || (dataRank == 4) ), std::invalid_argument,
1721 ">>> ERROR (ArrayTools::matmatProductDataField): inputData rank 2, 3 or 4 required.")
1722 } // switch inputData rank
1723 }
1724 else{ // constant data case
1725 switch(dataRank){
1726 case 2:
1727 for(size_t cell = 0; cell < numCells; cell++) {
1728 for(size_t field = 0; field < numFields; field++) {
1729 for(size_t point = 0; point < numPoints; point++) {
1730 for(size_t row = 0; row < matDim; row++) {
1731 for(size_t col = 0; col < matDim; col++) {
1732 outputFieldsWrap(cell, field, point, row, col) = \
1733 inputDataWrap(cell, 0)*inputFieldsWrap(cell, field, point, row, col);
1734 }// Col-loop
1735 } // Row-loop
1736 } // P-loop
1737 } // F-loop
1738 }// C-loop
1739 break;
1740
1741 case 3:
1742 for(size_t cell = 0; cell < numCells; cell++) {
1743 for(size_t field = 0; field < numFields; field++) {
1744 for(size_t point = 0; point < numPoints; point++) {
1745 for(size_t row = 0; row < matDim; row++) {
1746 for(size_t col = 0; col < matDim; col++) {
1747 outputFieldsWrap(cell, field, point, row, col) = \
1748 inputDataWrap(cell, 0, row)*inputFieldsWrap(cell, field, point, row, col);
1749 }// Col-loop
1750 } // Row-loop
1751 } // P-loop
1752 } // F-loop
1753 }// C-loop
1754 break;
1755
1756 case 4:
1757 if ((transpose == 'n') || (transpose == 'N')) {
1758 for(size_t cell = 0; cell < numCells; cell++){
1759 for(size_t field = 0; field < numFields; field++){
1760 for(size_t point = 0; point < numPoints; point++){
1761 for(size_t row = 0; row < matDim; row++){
1762 for(size_t col = 0; col < matDim; col++){
1763 outputFieldsWrap(cell, field, point, row, col) = 0.0;
1764 for(size_t i = 0; i < matDim; i++){
1765 outputFieldsWrap(cell, field, point, row, col) += \
1766 inputDataWrap(cell, 0, row, i)*inputFieldsWrap(cell, field, point, i, col);
1767 }// i
1768 } // col
1769 } //row
1770 }// point
1771 }// field
1772 }// cell
1773 } // no transpose
1774 else if ((transpose == 't') || (transpose == 'T')) {
1775 for(size_t cell = 0; cell < numCells; cell++){
1776 for(size_t field = 0; field < numFields; field++){
1777 for(size_t point = 0; point < numPoints; point++){
1778 for(size_t row = 0; row < matDim; row++){
1779 for(size_t col = 0; col < matDim; col++){
1780 outputFieldsWrap(cell, field, point, row, col) = 0.0;
1781 for(size_t i = 0; i < matDim; i++){
1782 outputFieldsWrap(cell, field, point, row, col) += \
1783 inputDataWrap(cell, 0, i, row)*inputFieldsWrap(cell, field, point, i, col);
1784 }// i
1785 } // col
1786 } //row
1787 }// point
1788 }// field
1789 }// cell
1790 } //transpose
1791 else {
1792 TEUCHOS_TEST_FOR_EXCEPTION( !( (transpose == 'n') || (transpose == 'N') || (transpose == 't') || (transpose == 'T') ), std::invalid_argument,
1793 ">>> ERROR (ArrayTools::matmatProductDataField): The transpose flag must be 'n', 'N', 't' or 'T'.");
1794 }
1795 break;
1796
1797 default:
1798 TEUCHOS_TEST_FOR_EXCEPTION( !( (dataRank == 2) || (dataRank == 3) || (dataRank == 4) ), std::invalid_argument,
1799 ">>> ERROR (ArrayTools::matmatProductDataField): inputData rank 2, 3 or 4 required.")
1800 } // switch inputData rank
1801 } // end constant data case
1802 } // inputFields rank 5
1803 /**********************************************************************************************
1804 * inputFields is (F,P,D,D) *
1805 *********************************************************************************************/
1806 else if(inRank == 4) {
1807 if(numDataPts != 1){ // non-constant data
1808
1809 switch(dataRank){
1810 case 2:
1811 for(size_t cell = 0; cell < numCells; cell++) {
1812 for(size_t field = 0; field < numFields; field++) {
1813 for(size_t point = 0; point < numPoints; point++) {
1814 for(size_t row = 0; row < matDim; row++) {
1815 for(size_t col = 0; col < matDim; col++) {
1816 outputFieldsWrap(cell, field, point, row, col) = \
1817 inputDataWrap(cell, point)*inputFieldsWrap(field, point, row, col);
1818 }// Col-loop
1819 } // Row-loop
1820 } // P-loop
1821 } // F-loop
1822 }// C-loop
1823 break;
1824
1825 case 3:
1826 for(size_t cell = 0; cell < numCells; cell++) {
1827 for(size_t field = 0; field < numFields; field++) {
1828 for(size_t point = 0; point < numPoints; point++) {
1829 for(size_t row = 0; row < matDim; row++) {
1830 for(size_t col = 0; col < matDim; col++) {
1831 outputFieldsWrap(cell, field, point, row, col) = \
1832 inputDataWrap(cell, point, row)*inputFieldsWrap(field, point, row, col);
1833 }// Col-loop
1834 } // Row-loop
1835 } // P-loop
1836 } // F-loop
1837 }// C-loop
1838 break;
1839
1840 case 4:
1841 if ((transpose == 'n') || (transpose == 'N')) {
1842 for(size_t cell = 0; cell < numCells; cell++){
1843 for(size_t field = 0; field < numFields; field++){
1844 for(size_t point = 0; point < numPoints; point++){
1845 for(size_t row = 0; row < matDim; row++){
1846 for(size_t col = 0; col < matDim; col++){
1847 outputFieldsWrap(cell, field, point, row, col) = 0.0;
1848 for(size_t i = 0; i < matDim; i++){
1849 outputFieldsWrap(cell, field, point, row, col) += \
1850 inputDataWrap(cell, point, row, i)*inputFieldsWrap(field, point, i, col);
1851 }// i
1852 } // col
1853 } //row
1854 }// point
1855 }// field
1856 }// cell
1857 } // no transpose
1858 else if ((transpose == 't') || (transpose == 'T')) {
1859 for(size_t cell = 0; cell < numCells; cell++){
1860 for(size_t field = 0; field < numFields; field++){
1861 for(size_t point = 0; point < numPoints; point++){
1862 for(size_t row = 0; row < matDim; row++){
1863 for(size_t col = 0; col < matDim; col++){
1864 outputFieldsWrap(cell, field, point, row, col) = 0.0;
1865 for(size_t i = 0; i < matDim; i++){
1866 outputFieldsWrap(cell, field, point, row, col) += \
1867 inputDataWrap(cell, point, i, row)*inputFieldsWrap(field, point, i, col);
1868 }// i
1869 } // col
1870 } //row
1871 }// point
1872 }// field
1873 }// cell
1874 } //transpose
1875 else {
1876 TEUCHOS_TEST_FOR_EXCEPTION( !( (transpose == 'n') || (transpose == 'N') || (transpose == 't') || (transpose == 'T') ), std::invalid_argument,
1877 ">>> ERROR (ArrayTools::matmatProductDataField): The transpose flag must be 'n', 'N', 't' or 'T'.");
1878 }
1879 break;
1880
1881 default:
1882 TEUCHOS_TEST_FOR_EXCEPTION( !( (dataRank == 2) || (dataRank == 3) || (dataRank == 4) ), std::invalid_argument,
1883 ">>> ERROR (ArrayTools::matmatProductDataField): inputData rank 2, 3 or 4 required.")
1884 } // switch inputData rank
1885 }
1886 else{ // constant data case
1887 switch(dataRank){
1888 case 2:
1889 for(size_t cell = 0; cell < numCells; cell++) {
1890 for(size_t field = 0; field < numFields; field++) {
1891 for(size_t point = 0; point < numPoints; point++) {
1892 for(size_t row = 0; row < matDim; row++) {
1893 for(size_t col = 0; col < matDim; col++) {
1894 outputFieldsWrap(cell, field, point, row, col) = \
1895 inputDataWrap(cell, 0)*inputFieldsWrap(field, point, row, col);
1896 }// Col-loop
1897 } // Row-loop
1898 } // P-loop
1899 } // F-loop
1900 }// C-loop
1901 break;
1902
1903 case 3:
1904 for(size_t cell = 0; cell < numCells; cell++) {
1905 for(size_t field = 0; field < numFields; field++) {
1906 for(size_t point = 0; point < numPoints; point++) {
1907 for(size_t row = 0; row < matDim; row++) {
1908 for(size_t col = 0; col < matDim; col++) {
1909 outputFieldsWrap(cell, field, point, row, col) = \
1910 inputDataWrap(cell, 0, row)*inputFieldsWrap(field, point, row, col);
1911 }// Col-loop
1912 } // Row-loop
1913 } // P-loop
1914 } // F-loop
1915 }// C-loop
1916 break;
1917
1918 case 4:
1919 if ((transpose == 'n') || (transpose == 'N')) {
1920 for(size_t cell = 0; cell < numCells; cell++){
1921 for(size_t field = 0; field < numFields; field++){
1922 for(size_t point = 0; point < numPoints; point++){
1923 for(size_t row = 0; row < matDim; row++){
1924 for(size_t col = 0; col < matDim; col++){
1925 outputFieldsWrap(cell, field, point, row, col) = 0.0;
1926 for(size_t i = 0; i < matDim; i++){
1927 outputFieldsWrap(cell, field, point, row, col) += \
1928 inputDataWrap(cell, 0, row, i)*inputFieldsWrap(field, point, i, col);
1929 }// i
1930 } // col
1931 } //row
1932 }// point
1933 }// field
1934 }// cell
1935 } // no transpose
1936 else if ((transpose == 't') || (transpose == 'T')) {
1937 for(size_t cell = 0; cell < numCells; cell++){
1938 for(size_t field = 0; field < numFields; field++){
1939 for(size_t point = 0; point < numPoints; point++){
1940 for(size_t row = 0; row < matDim; row++){
1941 for(size_t col = 0; col < matDim; col++){
1942 outputFieldsWrap(cell, field, point, row, col) = 0.0;
1943 for(size_t i = 0; i < matDim; i++){
1944 outputFieldsWrap(cell, field, point, row, col) += \
1945 inputDataWrap(cell, 0, i, row)*inputFieldsWrap(field, point, i, col);
1946 }// i
1947 } // col
1948 } //row
1949 }// point
1950 }// field
1951 }// cell
1952 } //transpose
1953 else {
1954 TEUCHOS_TEST_FOR_EXCEPTION( !( (transpose == 'n') || (transpose == 'N') || (transpose == 't') || (transpose == 'T') ), std::invalid_argument,
1955 ">>> ERROR (ArrayTools::matmatProductDataField): The transpose flag must be 'n', 'N', 't' or 'T'.");
1956 }
1957 break;
1958
1959 default:
1960 TEUCHOS_TEST_FOR_EXCEPTION( !( (dataRank == 2) || (dataRank == 3) || (dataRank == 4) ), std::invalid_argument,
1961 ">>> ERROR (ArrayTools::matmatProductDataField): inputData rank 2, 3 or 4 required.")
1962 } // switch inputData rank
1963 } // end constant data case
1964 } // inputFields rank 4
1965 else {
1966 TEUCHOS_TEST_FOR_EXCEPTION( true, std::invalid_argument,
1967 ">>> ERROR (ArrayTools::matmatProductDataField): inputFields rank 4 or 5 required.")
1968 }// rank error
1969 }
1970
1971
1972 template<class Scalar,
1973 class ArrayOutData,
1974 class ArrayInDataLeft,
1975 class ArrayInDataRight>
1976 void ArrayTools::matmatProductDataData(ArrayOutData & outputData,
1977 const ArrayInDataLeft & inputDataLeft,
1978 const ArrayInDataRight & inputDataRight,
1979 const char transpose){
1980
1981#ifdef HAVE_INTREPID_DEBUG
1982 std::string errmsg = ">>> ERROR (ArrayTools::matmatProductDataData):";
1983 /*
1984 * Check array rank and spatial dimension range (if applicable)
1985 * (1) inputDataLeft(C,P), (C,P,D) or (C,P,D,D); P=1 is admissible to allow multiply by const. left data
1986 * (2) inputDataRight(C,P,D,D) or (P,D,D);
1987 * (3) outputData(C,P,D,D)
1988 */
1989 // (1) inputDataLeft is (C,P), (C,P,D) or (C,P,D,D) and 1 <= D <= 3 is required
1990 TEUCHOS_TEST_FOR_EXCEPTION( !requireRankRange(errmsg, inputDataLeft, 2,4),
1991 std::invalid_argument, errmsg);
1992 if(getrank(inputDataLeft) > 2) {
1993 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, inputDataLeft, 2, 1,3),
1994 std::invalid_argument, errmsg);
1995 }
1996 if(getrank(inputDataLeft) == 4) {
1997 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, inputDataLeft, 3, 1,3),
1998 std::invalid_argument, errmsg);
1999 }
2000 // (2) inputDataRight is (C,P,D,D) or (P,D,D) and 1 <= (D=dimension(rank-1),(rank-2)) <= 3 is required.
2001 TEUCHOS_TEST_FOR_EXCEPTION( !requireRankRange(errmsg, inputDataRight, 3,4),
2002 std::invalid_argument, errmsg);
2003 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, inputDataRight, getrank(inputDataRight)-1, 1,3),
2004 std::invalid_argument, errmsg);
2005 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, inputDataRight, getrank(inputDataRight)-2, 1,3),
2006 std::invalid_argument, errmsg);
2007 // (3) outputData is (C,P,D,D)
2008 TEUCHOS_TEST_FOR_EXCEPTION( !requireRankRange(errmsg, outputData, 4, 4), std::invalid_argument, errmsg);
2009 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, outputData, 2, 1,3),
2010 std::invalid_argument, errmsg);
2011 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, outputData, 3, 1,3),
2012 std::invalid_argument, errmsg);
2013 /*
2014 * Dimension cross-checks:
2015 * (1) inputDataLeft vs. inputDataRight
2016 * (2) outputData vs. inputDataLeft
2017 * (3) outputData vs. inputDataRight
2018 *
2019 * Cross-check (2): outputData(C,P,D,D) vs. inputDataLeft(C,P), (C,P,D) or (C,P,D,D):
2020 * dimensions C, and D must match in all cases, dimension P must match only when non-constant
2021 * data is specified (P>1). Do not check P dimensions with constant left data, i.e., when P=1 in
2022 * inputDataLeft(C,1,...)
2023 */
2024 if(inputDataLeft.dimension(1) > 1){ // check P dimension if P>1 in inputDataLeft
2025 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg,
2026 outputData, 1,
2027 inputDataLeft, 1),
2028 std::invalid_argument, errmsg);
2029 }
2030 if(getrank(inputDataLeft) == 2){ // inputDataLeft(C,P): check C
2031 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg,
2032 outputData, 0,
2033 inputDataLeft, 0),
2034 std::invalid_argument, errmsg);
2035 }
2036 if(getrank(inputDataLeft) == 3){ // inputDataLeft(C,P,D): check C and D
2037 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg,
2038 outputData, 0,2,3,
2039 inputDataLeft, 0,2,2),
2040 std::invalid_argument, errmsg);
2041 }
2042 if(getrank(inputDataLeft) == 4){ // inputDataLeft(C,P,D,D): check C and D
2043 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg,
2044 outputData, 0,2,3,
2045 inputDataLeft, 0,2,3),
2046 std::invalid_argument, errmsg);
2047 }
2048 /*
2049 * Cross-checks (1,3):
2050 */
2051 if( getrank(inputDataRight) == 4) {
2052 // 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
2053 if(inputDataLeft.dimension(1) > 1){ // check P dimension if P>1 in inputDataLeft
2054 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg,
2055 inputDataLeft, 1,
2056 inputDataRight, 1),
2057 std::invalid_argument, errmsg);
2058 }
2059 if(getrank(inputDataLeft) == 2){ // inputDataLeft(C,P): check C
2060 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg,
2061 inputDataLeft, 0,
2062 inputDataRight, 0),
2063 std::invalid_argument, errmsg);
2064 }
2065 if(getrank(inputDataLeft) == 3){ // inputDataLeft(C,P,D): check C and D
2066 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg,
2067 inputDataLeft, 0,2,2,
2068 inputDataRight, 0,2,3),
2069 std::invalid_argument, errmsg);
2070 }
2071 if(getrank(inputDataLeft) == 4){ // inputDataLeft(C,P,D,D): check C and D
2072 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg,
2073 inputDataLeft, 0,2,3,
2074 inputDataRight, 0,2,3),
2075 std::invalid_argument, errmsg);
2076 }
2077 // Cross-check (3): outputData(C,P,D,D) vs. inputDataRight(C,P,D,D): all dimensions C, P, D must match
2078 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, outputData, inputDataRight),
2079 std::invalid_argument, errmsg);
2080 }
2081 else{
2082 // Cross-check (1): inputDataLeft(C,P), (C,P,D), or (C,P,D,D) vs. inputDataRight(P,D,D): dimensions P, D must match
2083 if(inputDataLeft.dimension(1) > 1){ // check P if P>1 in inputData
2084 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg,
2085 inputDataLeft, 1,
2086 inputDataRight, 0),
2087 std::invalid_argument, errmsg);
2088 }
2089 if(getrank(inputDataLeft) == 3){ // inputDataLeft(C,P,D): check D
2090 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg,
2091 inputDataLeft, 2,2,
2092 inputDataRight, 1,2),
2093 std::invalid_argument, errmsg);
2094 }
2095 if(getrank(inputDataLeft) == 4){ // inputDataLeft(C,P,D,D): check D
2096 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg,
2097 inputDataLeft, 2,3,
2098 inputDataRight, 1,2),
2099 std::invalid_argument, errmsg);
2100 }
2101 // Cross-check (3): outputData(C,P,D,D) vs. inputDataRight(P,D,D): dimensions P, D must match
2102 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg,
2103 outputData, 1,2,3,
2104 inputDataRight, 0,1,2),
2105 std::invalid_argument, errmsg);
2106 }
2107#endif
2108
2109
2110 ArrayWrapper<Scalar,ArrayOutData, Rank<ArrayOutData >::value, false>outputDataWrap(outputData);
2111 ArrayWrapper<Scalar,ArrayInDataLeft, Rank<ArrayInDataLeft >::value, true>inputDataLeftWrap(inputDataLeft);
2112 ArrayWrapper<Scalar,ArrayInDataRight, Rank<ArrayInDataRight >::value, true>inputDataRightWrap(inputDataRight);
2113
2114 size_t dataLeftRank = getrank(inputDataLeft);
2115 size_t numDataLeftPts = inputDataLeft.dimension(1);
2116 size_t dataRightRank = getrank(inputDataRight);
2117 size_t numCells = outputData.dimension(0);
2118 size_t numPoints = outputData.dimension(1);
2119 size_t matDim = outputData.dimension(2);
2120
2121 /*********************************************************************************************
2122 * inputDataRight is (C,P,D,D) *
2123 *********************************************************************************************/
2124 if(dataRightRank == 4){
2125 if(numDataLeftPts != 1){ // non-constant data
2126
2127 switch(dataLeftRank){
2128 case 2:
2129 for(size_t cell = 0; cell < numCells; cell++) {
2130 for(size_t point = 0; point < numPoints; point++) {
2131 for(size_t row = 0; row < matDim; row++) {
2132 for(size_t col = 0; col < matDim; col++) {
2133 outputDataWrap(cell, point, row, col) = \
2134 inputDataLeftWrap(cell, point)*inputDataRightWrap(cell, point, row, col);
2135 }// Col-loop
2136 } // Row-loop
2137 } // P-loop
2138 }// C-loop
2139 break;
2140
2141 case 3:
2142 for(size_t cell = 0; cell < numCells; cell++) {
2143 for(size_t point = 0; point < numPoints; point++) {
2144 for(size_t row = 0; row < matDim; row++) {
2145 for(size_t col = 0; col < matDim; col++) {
2146 outputDataWrap(cell, point, row, col) = \
2147 inputDataLeftWrap(cell, point, row)*inputDataRightWrap(cell, point, row, col);
2148 }// Col-loop
2149 } // Row-loop
2150 } // P-loop
2151 }// C-loop
2152 break;
2153
2154 case 4:
2155 if ((transpose == 'n') || (transpose == 'N')) {
2156 for(size_t cell = 0; cell < numCells; cell++){
2157 for(size_t point = 0; point < numPoints; point++){
2158 for(size_t row = 0; row < matDim; row++){
2159 for(size_t col = 0; col < matDim; col++){
2160 outputDataWrap(cell, point, row, col) = 0.0;
2161 for(size_t i = 0; i < matDim; i++){
2162 outputDataWrap(cell, point, row, col) += \
2163 inputDataLeftWrap(cell, point, row, i)*inputDataRightWrap(cell, point, i, col);
2164 }// i
2165 } // col
2166 } //row
2167 }// point
2168 }// cell
2169 } // no transpose
2170 else if ((transpose == 't') || (transpose == 'T')) {
2171 for(size_t cell = 0; cell < numCells; cell++){
2172 for(size_t point = 0; point < numPoints; point++){
2173 for(size_t row = 0; row < matDim; row++){
2174 for(size_t col = 0; col < matDim; col++){
2175 outputDataWrap(cell, point, row, col) = 0.0;
2176 for(size_t i = 0; i < matDim; i++){
2177 outputDataWrap(cell, point, row, col) += \
2178 inputDataLeftWrap(cell, point, i, row)*inputDataRightWrap(cell, point, i, col);
2179 }// i
2180 } // col
2181 } //row
2182 }// point
2183 }// cell
2184 } //transpose
2185 else {
2186 TEUCHOS_TEST_FOR_EXCEPTION( !( (transpose == 'n') || (transpose == 'N') || (transpose == 't') || (transpose == 'T') ), std::invalid_argument,
2187 ">>> ERROR (ArrayTools::matmatProductDataData): The transpose flag must be 'n', 'N', 't' or 'T'.");
2188 }
2189 break;
2190
2191 default:
2192 TEUCHOS_TEST_FOR_EXCEPTION( !( (dataLeftRank == 2) || (dataLeftRank == 3) || (dataLeftRank == 4) ), std::invalid_argument,
2193 ">>> ERROR (ArrayTools::matmatProductDataData): inputDataLeft rank 2, 3 or 4 required.")
2194 } // switch inputData rank
2195 }
2196 else{ // constant data case
2197 switch(dataLeftRank){
2198 case 2:
2199 for(size_t cell = 0; cell < numCells; cell++) {
2200 for(size_t point = 0; point < numPoints; point++) {
2201 for(size_t row = 0; row < matDim; row++) {
2202 for(size_t col = 0; col < matDim; col++) {
2203 outputDataWrap(cell, point, row, col) = \
2204 inputDataLeftWrap(cell, 0)*inputDataRightWrap(cell, point, row, col);
2205 }// Col-loop
2206 } // Row-loop
2207 } // P-loop
2208 }// C-loop
2209 break;
2210
2211 case 3:
2212 for(size_t cell = 0; cell < numCells; cell++) {
2213 for(size_t point = 0; point < numPoints; point++) {
2214 for(size_t row = 0; row < matDim; row++) {
2215 for(size_t col = 0; col < matDim; col++) {
2216 outputDataWrap(cell, point, row, col) = \
2217 inputDataLeftWrap(cell, 0, row)*inputDataRightWrap(cell, point, row, col);
2218 }// Col-loop
2219 } // Row-loop
2220 } // P-loop
2221 }// C-loop
2222 break;
2223
2224 case 4:
2225 if ((transpose == 'n') || (transpose == 'N')) {
2226 for(size_t cell = 0; cell < numCells; cell++){
2227 for(size_t point = 0; point < numPoints; point++){
2228 for(size_t row = 0; row < matDim; row++){
2229 for(size_t col = 0; col < matDim; col++){
2230 outputDataWrap(cell, point, row, col) = 0.0;
2231 for(size_t i = 0; i < matDim; i++){
2232 outputDataWrap(cell, point, row, col) += \
2233 inputDataLeftWrap(cell, 0, row, i)*inputDataRightWrap(cell, point, i, col);
2234 }// i
2235 } // col
2236 } //row
2237 }// point
2238 }// cell
2239 } // no transpose
2240 else if ((transpose == 't') || (transpose == 'T')) {
2241 for(size_t cell = 0; cell < numCells; cell++){
2242 for(size_t point = 0; point < numPoints; point++){
2243 for(size_t row = 0; row < matDim; row++){
2244 for(size_t col = 0; col < matDim; col++){
2245 outputDataWrap(cell, point, row, col) = 0.0;
2246 for(size_t i = 0; i < matDim; i++){
2247 outputDataWrap(cell, point, row, col) += \
2248 inputDataLeftWrap(cell, 0, i, row)*inputDataRightWrap(cell, point, i, col);
2249 }// i
2250 } // col
2251 } //row
2252 }// point
2253 }// cell
2254 } //transpose
2255 else {
2256 TEUCHOS_TEST_FOR_EXCEPTION( !( (transpose == 'n') || (transpose == 'N') || (transpose == 't') || (transpose == 'T') ), std::invalid_argument,
2257 ">>> ERROR (ArrayTools::matmatProductDataData): The transpose flag must be 'n', 'N', 't' or 'T'.");
2258 }
2259 break;
2260
2261 default:
2262 TEUCHOS_TEST_FOR_EXCEPTION( !( (dataLeftRank == 2) || (dataLeftRank == 3) || (dataLeftRank == 4) ), std::invalid_argument,
2263 ">>> ERROR (ArrayTools::matmatProductDataData): inputDataLeft rank 2, 3 or 4 required.")
2264 } // switch inputDataLeft rank
2265 } // end constant data case
2266 } // inputDataRight rank 4
2267 /**********************************************************************************************
2268 * inputDataRight is (P,D,D) *
2269 *********************************************************************************************/
2270 else if(dataRightRank == 3) {
2271 if(numDataLeftPts != 1){ // non-constant data
2272
2273 switch(dataLeftRank){
2274 case 2:
2275 for(size_t cell = 0; cell < numCells; cell++) {
2276 for(size_t point = 0; point < numPoints; point++) {
2277 for(size_t row = 0; row < matDim; row++) {
2278 for(size_t col = 0; col < matDim; col++) {
2279 outputDataWrap(cell, point, row, col) = \
2280 inputDataLeftWrap(cell, point)*inputDataRightWrap(point, row, col);
2281 }// Col-loop
2282 } // Row-loop
2283 } // P-loop
2284 }// C-loop
2285 break;
2286
2287 case 3:
2288 for(size_t cell = 0; cell < numCells; cell++) {
2289 for(size_t point = 0; point < numPoints; point++) {
2290 for(size_t row = 0; row < matDim; row++) {
2291 for(size_t col = 0; col < matDim; col++) {
2292 outputDataWrap(cell, point, row, col) = \
2293 inputDataLeftWrap(cell, point, row)*inputDataRightWrap(point, row, col);
2294 }// Col-loop
2295 } // Row-loop
2296 } // P-loop
2297 }// C-loop
2298 break;
2299
2300 case 4:
2301 if ((transpose == 'n') || (transpose == 'N')) {
2302 for(size_t cell = 0; cell < numCells; cell++){
2303 for(size_t point = 0; point < numPoints; point++){
2304 for(size_t row = 0; row < matDim; row++){
2305 for(size_t col = 0; col < matDim; col++){
2306 outputDataWrap(cell, point, row, col) = 0.0;
2307 for(size_t i = 0; i < matDim; i++){
2308 outputDataWrap(cell, point, row, col) += \
2309 inputDataLeftWrap(cell, point, row, i)*inputDataRightWrap(point, i, col);
2310 }// i
2311 } // col
2312 } //row
2313 }// point
2314 }// cell
2315 } // no transpose
2316 else if ((transpose == 't') || (transpose == 'T')) {
2317 for(size_t cell = 0; cell < numCells; cell++){
2318 for(size_t point = 0; point < numPoints; point++){
2319 for(size_t row = 0; row < matDim; row++){
2320 for(size_t col = 0; col < matDim; col++){
2321 outputDataWrap(cell, point, row, col) = 0.0;
2322 for(size_t i = 0; i < matDim; i++){
2323 outputDataWrap(cell, point, row, col) += \
2324 inputDataLeftWrap(cell, point, i, row)*inputDataRightWrap(point, i, col);
2325 }// i
2326 } // col
2327 } //row
2328 }// point
2329 }// cell
2330 } //transpose
2331 else {
2332 TEUCHOS_TEST_FOR_EXCEPTION( !( (transpose == 'n') || (transpose == 'N') || (transpose == 't') || (transpose == 'T') ), std::invalid_argument,
2333 ">>> ERROR (ArrayTools::matmatProductDataData): The transpose flag must be 'n', 'N', 't' or 'T'.");
2334 }
2335 break;
2336
2337 default:
2338 TEUCHOS_TEST_FOR_EXCEPTION( !( (dataLeftRank == 2) || (dataLeftRank == 3) || (dataLeftRank == 4) ), std::invalid_argument,
2339 ">>> ERROR (ArrayTools::matmatProductDataData): inputDataLeft rank 2, 3 or 4 required.")
2340 } // switch inputDataLeft rank
2341 }
2342 else{ // constant data case
2343 switch(dataLeftRank){
2344 case 2:
2345 for(size_t cell = 0; cell < numCells; cell++) {
2346 for(size_t point = 0; point < numPoints; point++) {
2347 for(size_t row = 0; row < matDim; row++) {
2348 for(size_t col = 0; col < matDim; col++) {
2349 outputDataWrap(cell, point, row, col) = \
2350 inputDataLeftWrap(cell, 0)*inputDataRightWrap(point, row, col);
2351 }// Col-loop
2352 } // Row-loop
2353 } // P-loop
2354 }// C-loop
2355 break;
2356
2357 case 3:
2358 for(size_t cell = 0; cell < numCells; cell++) {
2359 for(size_t point = 0; point < numPoints; point++) {
2360 for(size_t row = 0; row < matDim; row++) {
2361 for(size_t col = 0; col < matDim; col++) {
2362 outputDataWrap(cell, point, row, col) = \
2363 inputDataLeftWrap(cell, 0, row)*inputDataRightWrap(point, row, col);
2364 }// Col-loop
2365 } // Row-loop
2366 } // P-loop
2367 }// C-loop
2368 break;
2369
2370 case 4:
2371 if ((transpose == 'n') || (transpose == 'N')) {
2372 for(size_t cell = 0; cell < numCells; cell++){
2373 for(size_t point = 0; point < numPoints; point++){
2374 for(size_t row = 0; row < matDim; row++){
2375 for(size_t col = 0; col < matDim; col++){
2376 outputDataWrap(cell, point, row, col) = 0.0;
2377 for(size_t i = 0; i < matDim; i++){
2378 outputDataWrap(cell, point, row, col) += \
2379 inputDataLeftWrap(cell, 0, row, i)*inputDataRightWrap(point, i, col);
2380 }// i
2381 } // col
2382 } //row
2383 }// point
2384 }// cell
2385 } // no transpose
2386 else if ((transpose == 't') || (transpose == 'T')) {
2387 for(size_t cell = 0; cell < numCells; cell++){
2388 for(size_t point = 0; point < numPoints; point++){
2389 for(size_t row = 0; row < matDim; row++){
2390 for(size_t col = 0; col < matDim; col++){
2391 outputDataWrap(cell, point, row, col) = 0.0;
2392 for(size_t i = 0; i < matDim; i++){
2393 outputDataWrap(cell, point, row, col) += \
2394 inputDataLeftWrap(cell, 0, i, row)*inputDataRightWrap(point, i, col);
2395 }// i
2396 } // col
2397 } //row
2398 }// point
2399 }// cell
2400 } //transpose
2401 else {
2402 TEUCHOS_TEST_FOR_EXCEPTION( !( (transpose == 'n') || (transpose == 'N') || (transpose == 't') || (transpose == 'T') ), std::invalid_argument,
2403 ">>> ERROR (ArrayTools::matmatProductDataData): The transpose flag must be 'n', 'N', 't' or 'T'.");
2404 }
2405 break;
2406
2407 default:
2408 TEUCHOS_TEST_FOR_EXCEPTION( !( (dataLeftRank == 2) || (dataLeftRank == 3) || (dataLeftRank == 4) ), std::invalid_argument,
2409 ">>> ERROR (ArrayTools::matmatProductDataData): inputDataLeft rank 2, 3 or 4 required.")
2410 } // switch inputDataLeft rank
2411 } // end constant data case
2412 } // inputDataRight rank 3
2413 else {
2414 TEUCHOS_TEST_FOR_EXCEPTION( true, std::invalid_argument,
2415 ">>> ERROR (ArrayTools::matmatProductDataData): inputDataRight rank 3 or 4 required.")
2416 }// rank error
2417 }
2418
2419
2420
2421
2422} // end namespace Intrepid
2423
2424
2425
2426
2427
2428
2429
2430
2431
2432
2433
2434
2435
2436
2437
2438
2439
2440
2441
2442
2443
bool requireRankRange(std::string &errmsg, const Array &array, const int lowerBound, const int upperBound)
Checks if the rank of the array argument is in the required range.
bool requireDimensionRange(std::string &errmsg, const Array &array, const int dim, const int lowerBound, const int upperBound)
Checks if the specified array dimension is in the required range.
bool requireDimensionMatch(std::string &errmsg, const Array1 &array1, const int a1_dim0, const Array2 &array2, const int a2_dim0)
Checks arrays for a single matching dimension.
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,...