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 
49 namespace 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
139 ArrayWrapper<Scalar,ArrayOutFields, Rank<ArrayOutFields >::value, false>outputFieldsWrap(outputFields);
140 ArrayWrapper<Scalar,ArrayInData, Rank<ArrayInData >::value, true>inputDataWrap(inputData);
141 ArrayWrapper<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 
327 ArrayWrapper<Scalar,ArrayOutData, Rank<ArrayOutData >::value, false>outputDataWrap(outputData);
328 ArrayWrapper<Scalar,ArrayInDataLeft, Rank<ArrayInDataLeft >::value, true>inputDataLeftWrap(inputDataLeft);
329 ArrayWrapper<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
484 ArrayWrapper<Scalar,ArrayOutFields, Rank<ArrayOutFields >::value, false>outputFieldsWrap(outputFields);
485 ArrayWrapper<Scalar,ArrayInData, Rank<ArrayInData >::value, true>inputDataWrap(inputData);
486 ArrayWrapper<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 
599 ArrayWrapper<Scalar,ArrayOutData, Rank<ArrayOutData >::value, false>outputDataWrap(outputData);
600 ArrayWrapper<Scalar,ArrayInDataLeft, Rank<ArrayInDataLeft >::value, true>inputDataLeftWrap(inputDataLeft);
601 ArrayWrapper<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 
2444 
2445 #if defined(Intrepid_SHOW_DEPRECATED_WARNINGS)
2446 #ifdef __GNUC__
2447 #warning "The Intrepid package is deprecated"
2448 #endif
2449 #endif
2450 
Intrepid::ArrayTools::matmatProductDataData
static void matmatProductDataData(ArrayOutData &outputData, const ArrayInDataLeft &inputDataLeft, const ArrayInDataRight &inputDataRight, const char transpose= 'N')
There are two use cases: (1) matrix-matrix product of a rank-4 container inputDataRight with dimensio...
Definition: Intrepid_ArrayToolsDefTensor.hpp:1976
Intrepid::ArrayTools::crossProductDataField
static void crossProductDataField(ArrayOutFields &outputFields, const ArrayInData &inputData, const ArrayInFields &inputFields)
There are two use cases: (1) cross product of a rank-4 container inputFields with dimensions (C...
Definition: Intrepid_ArrayToolsDefTensor.hpp:55
Intrepid::ArrayTools::crossProductDataData
static void crossProductDataData(ArrayOutData &outputData, const ArrayInDataLeft &inputDataLeft, const ArrayInDataRight &inputDataRight)
There are two use cases: (1) cross product of a rank-3 container inputDataRight with dimensions (C...
Definition: Intrepid_ArrayToolsDefTensor.hpp:238
Intrepid::ArrayTools::matmatProductDataField
static void matmatProductDataField(ArrayOutFields &outputFields, const ArrayInData &inputData, const ArrayInFields &inputFields, const char transpose= 'N')
There are two use cases: (1) matrix-matrix product of a rank-5 container inputFields with dimensions ...
Definition: Intrepid_ArrayToolsDefTensor.hpp:1497
Intrepid::ArrayTools::outerProductDataData
static void outerProductDataData(ArrayOutData &outputData, const ArrayInDataLeft &inputDataLeft, const ArrayInDataRight &inputDataRight)
There are two use cases: (1) outer product of a rank-3 container inputDataRight with dimensions (C...
Definition: Intrepid_ArrayToolsDefTensor.hpp:532
Intrepid::ArrayTools::outerProductDataField
static void outerProductDataField(ArrayOutFields &outputFields, const ArrayInData &inputData, const ArrayInFields &inputFields)
There are two use cases: (1) outer product of a rank-4 container inputFields with dimensions (C...
Definition: Intrepid_ArrayToolsDefTensor.hpp:417
Intrepid::ArrayTools::matvecProductDataData
static void matvecProductDataData(ArrayOutData &outputData, const ArrayInDataLeft &inputDataLeft, const ArrayInDataRight &inputDataRight, const char transpose= 'N')
There are two use cases: (1) matrix-vector product of a rank-3 container inputDataRight with dimensio...
Definition: Intrepid_ArrayToolsDefTensor.hpp:1085
Intrepid::ArrayTools::matvecProductDataField
static void matvecProductDataField(ArrayOutFields &outputFields, const ArrayInData &inputData, const ArrayInFields &inputFields, const char transpose= 'N')
There are two use cases: (1) matrix-vector product of a rank-4 container inputFields with dimensions ...
Definition: Intrepid_ArrayToolsDefTensor.hpp:642
Generated by   doxygen 1.8.5