Intrepid
Intrepid_ArrayToolsDefContractions.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 template<class Scalar, class ArrayOutFields, class ArrayInFieldsLeft, class ArrayInFieldsRight>
51 void ArrayTools::contractFieldFieldScalar(ArrayOutFields & outputFields,
52  const ArrayInFieldsLeft & leftFields,
53  const ArrayInFieldsRight & rightFields,
54  const ECompEngine compEngine,
55  const bool sumInto) {
56 
57 
58 #ifdef HAVE_INTREPID_DEBUG
59  TEUCHOS_TEST_FOR_EXCEPTION( (getrank(leftFields) != 3 ), std::invalid_argument,
60  ">>> ERROR (ArrayTools::contractFieldFieldScalar): Rank of the left input argument must equal 3!");
61  TEUCHOS_TEST_FOR_EXCEPTION( (getrank(rightFields) != 3 ), std::invalid_argument,
62  ">>> ERROR (ArrayTools::contractFieldFieldScalar): Rank of right input argument must equal 3!");
63  TEUCHOS_TEST_FOR_EXCEPTION( (getrank(outputFields) != 3 ), std::invalid_argument,
64  ">>> ERROR (ArrayTools::contractFieldFieldScalar): Rank of output argument must equal 3!");
65  TEUCHOS_TEST_FOR_EXCEPTION( (leftFields.dimension(0) != rightFields.dimension(0) ), std::invalid_argument,
66  ">>> ERROR (ArrayTools::contractFieldFieldScalar): Zeroth dimensions (number of integration domains) of the left and right input containers must agree!");
67  TEUCHOS_TEST_FOR_EXCEPTION( (leftFields.dimension(2) != rightFields.dimension(2) ), std::invalid_argument,
68  ">>> ERROR (ArrayTools::contractFieldFieldScalar): Second dimensions (numbers of integration points) of the left and right input containers must agree!");
69  TEUCHOS_TEST_FOR_EXCEPTION( (outputFields.dimension(0) != rightFields.dimension(0) ), std::invalid_argument,
70  ">>> ERROR (ArrayTools::contractFieldFieldScalar): Zeroth dimensions (numbers of integration domains) of the input and output containers must agree!");
71  TEUCHOS_TEST_FOR_EXCEPTION( (outputFields.dimension(1) != leftFields.dimension(1) ), std::invalid_argument,
72  ">>> ERROR (ArrayTools::contractFieldFieldScalar): First dimension of output container and first dimension of left input container must agree!");
73  TEUCHOS_TEST_FOR_EXCEPTION( (outputFields.dimension(2) != rightFields.dimension(1) ), std::invalid_argument,
74  ">>> ERROR (ArrayTools::contractFieldFieldScalar): Second dimension of output container and first dimension of right input container must agree!");
75 
76  TEUCHOS_TEST_FOR_EXCEPTION( ( compEngine!=COMP_CPP && compEngine!=COMP_BLAS ), std::invalid_argument,
77  ">>> ERROR (ArrayTools::contractFieldFieldScalar): Computational engine not defined!");
78 #endif
79 
80  // get sizes
81  int numCells = leftFields.dimension(0);
82  int numLeftFields = leftFields.dimension(1);
83  int numRightFields = rightFields.dimension(1);
84  int numPoints = leftFields.dimension(2);
85 
86 
87  if (sumInto) {
88  for (int cl = 0; cl < numCells; cl++) {
89  for (int lbf = 0; lbf < numLeftFields; lbf++) {
90  for (int rbf = 0; rbf < numRightFields; rbf++) {
91  Scalar tmpVal(0);
92  for (int qp = 0; qp < numPoints; qp++) {
93  tmpVal += leftFields(cl, lbf, qp)*rightFields(cl, rbf, qp);
94  } // P-loop
95  outputFields(cl, lbf, rbf) += tmpVal;
96  } // R-loop
97  } // L-loop
98  } // C-loop
99  }
100  else {
101  for (int cl = 0; cl < numCells; cl++) {
102  for (int lbf = 0; lbf < numLeftFields; lbf++) {
103  for (int rbf = 0; rbf < numRightFields; rbf++) {
104  Scalar tmpVal(0);
105  for (int qp = 0; qp < numPoints; qp++) {
106  tmpVal += leftFields(cl, lbf, qp)*rightFields(cl, rbf, qp);
107  } // P-loop
108  outputFields(cl, lbf, rbf) = tmpVal;
109  } // R-loop
110  } // L-loop
111  } // C-loop
112  }
113 } // end contractFieldFieldScalar
114 
115 
116 template<class Scalar, class ArrayOutFields, class ArrayInFieldsLeft, class ArrayInFieldsRight>
117 void ArrayTools::contractFieldFieldVector(ArrayOutFields & outputFields,
118  const ArrayInFieldsLeft & leftFields,
119  const ArrayInFieldsRight & rightFields,
120  const ECompEngine compEngine,
121  const bool sumInto) {
122 
123 #ifdef HAVE_INTREPID_DEBUG
124  TEUCHOS_TEST_FOR_EXCEPTION( (getrank(leftFields) != 4 ), std::invalid_argument,
125  ">>> ERROR (ArrayTools::contractFieldFieldVector): Rank of the left input argument must equal 4!");
126  TEUCHOS_TEST_FOR_EXCEPTION( (getrank(rightFields) != 4 ), std::invalid_argument,
127  ">>> ERROR (ArrayTools::contractFieldFieldVector): Rank of right input argument must equal 4!");
128  TEUCHOS_TEST_FOR_EXCEPTION( (getrank(outputFields) != 3 ), std::invalid_argument,
129  ">>> ERROR (ArrayTools::contractFieldFieldVector): Rank of output argument must equal 3!");
130  TEUCHOS_TEST_FOR_EXCEPTION( (leftFields.dimension(0) != rightFields.dimension(0) ), std::invalid_argument,
131  ">>> ERROR (ArrayTools::contractFieldFieldVector): Zeroth dimensions (number of integration domains) of the left and right input containers must agree!");
132  TEUCHOS_TEST_FOR_EXCEPTION( (leftFields.dimension(2) != rightFields.dimension(2) ), std::invalid_argument,
133  ">>> ERROR (ArrayTools::contractFieldFieldVector): Second dimensions (numbers of integration points) of the left and right input containers must agree!");
134  TEUCHOS_TEST_FOR_EXCEPTION( (leftFields.dimension(3) != rightFields.dimension(3) ), std::invalid_argument,
135  ">>> ERROR (ArrayTools::contractFieldFieldVector): Third dimensions (numbers of vector components) of the left and right input containers must agree!");
136  TEUCHOS_TEST_FOR_EXCEPTION( (outputFields.dimension(0) != rightFields.dimension(0) ), std::invalid_argument,
137  ">>> ERROR (ArrayTools::contractFieldFieldVector): Zeroth dimensions (numbers of integration domains) of the input and output containers must agree!");
138  TEUCHOS_TEST_FOR_EXCEPTION( (outputFields.dimension(1) != leftFields.dimension(1) ), std::invalid_argument,
139  ">>> ERROR (ArrayTools::contractFieldFieldVector): First dimension of output container and first dimension of left input container must agree!");
140  TEUCHOS_TEST_FOR_EXCEPTION( (outputFields.dimension(2) != rightFields.dimension(1) ), std::invalid_argument,
141  ">>> ERROR (ArrayTools::contractFieldFieldVector): Second dimension of output container and first dimension of right input container must agree!");
142  TEUCHOS_TEST_FOR_EXCEPTION( ( compEngine!=COMP_CPP && compEngine!=COMP_BLAS ), std::invalid_argument,
143  ">>> ERROR (ArrayTools::contractFieldFieldScalar): Computational engine not defined!");
144 
145 #endif
146 
147  // get sizes
148  int numCells = leftFields.dimension(0);
149  int numLeftFields = leftFields.dimension(1);
150  int numRightFields = rightFields.dimension(1);
151  int numPoints = leftFields.dimension(2);
152  int dimVec = leftFields.dimension(3);
153 
154 
155  if (sumInto) {
156  for (int cl = 0; cl < numCells; cl++) {
157  for (int lbf = 0; lbf < numLeftFields; lbf++) {
158  for (int rbf = 0; rbf < numRightFields; rbf++) {
159  Scalar tmpVal(0);
160  for (int qp = 0; qp < numPoints; qp++) {
161  for (int iVec = 0; iVec < dimVec; iVec++) {
162  tmpVal += leftFields(cl, lbf, qp, iVec)*rightFields(cl, rbf, qp, iVec);
163  } //D-loop
164  } // P-loop
165  outputFields(cl, lbf, rbf) += tmpVal;
166  } // R-loop
167  } // L-loop
168  } // C-loop
169  }
170  else {
171  for (int cl = 0; cl < numCells; cl++) {
172  for (int lbf = 0; lbf < numLeftFields; lbf++) {
173  for (int rbf = 0; rbf < numRightFields; rbf++) {
174  Scalar tmpVal(0);
175  for (int qp = 0; qp < numPoints; qp++) {
176  for (int iVec = 0; iVec < dimVec; iVec++) {
177  tmpVal += leftFields(cl, lbf, qp, iVec)*rightFields(cl, rbf, qp, iVec);
178  } //D-loop
179  } // P-loop
180  outputFields(cl, lbf, rbf) = tmpVal;
181  } // R-loop
182  } // L-loop
183  } // C-loop
184  }
185 } // end contractFieldFieldVector
186 
187 
188 template<class Scalar, class ArrayOutFields, class ArrayInFieldsLeft, class ArrayInFieldsRight>
189 void ArrayTools::contractFieldFieldTensor(ArrayOutFields & outputFields,
190  const ArrayInFieldsLeft & leftFields,
191  const ArrayInFieldsRight & rightFields,
192  const ECompEngine compEngine,
193  const bool sumInto) {
194 
195 #ifdef HAVE_INTREPID_DEBUG
196  TEUCHOS_TEST_FOR_EXCEPTION( (getrank(leftFields) != 5 ), std::invalid_argument,
197  ">>> ERROR (ArrayTools::contractFieldFieldTensor): Rank of the left input argument must equal 5!");
198  TEUCHOS_TEST_FOR_EXCEPTION( (getrank(rightFields) != 5 ), std::invalid_argument,
199  ">>> ERROR (ArrayTools::contractFieldFieldTensor): Rank of right input argument must equal 5!");
200  TEUCHOS_TEST_FOR_EXCEPTION( (getrank(outputFields) != 3 ), std::invalid_argument,
201  ">>> ERROR (ArrayTools::contractFieldFieldTensor): Rank of output argument must equal 3!");
202  TEUCHOS_TEST_FOR_EXCEPTION( (leftFields.dimension(0) != rightFields.dimension(0) ), std::invalid_argument,
203  ">>> ERROR (ArrayTools::contractFieldFieldTensor): Zeroth dimensions (number of integration domains) of the left and right input containers must agree!");
204  TEUCHOS_TEST_FOR_EXCEPTION( (leftFields.dimension(2) != rightFields.dimension(2) ), std::invalid_argument,
205  ">>> ERROR (ArrayTools::contractFieldFieldTensor): Second dimensions (numbers of integration points) of the left and right input containers must agree!");
206  TEUCHOS_TEST_FOR_EXCEPTION( (leftFields.dimension(3) != rightFields.dimension(3) ), std::invalid_argument,
207  ">>> ERROR (ArrayTools::contractFieldFieldTensor): Third dimensions (first tensor dimensions) of the left and right input containers must agree!");
208  TEUCHOS_TEST_FOR_EXCEPTION( (leftFields.dimension(4) != rightFields.dimension(4) ), std::invalid_argument,
209  ">>> ERROR (ArrayTools::contractFieldFieldTensor): Fourth dimensions (second tensor dimensions) of the left and right input containers must agree!");
210  TEUCHOS_TEST_FOR_EXCEPTION( (outputFields.dimension(0) != rightFields.dimension(0) ), std::invalid_argument,
211  ">>> ERROR (ArrayTools::contractFieldFieldTensor): Zeroth dimensions (numbers of integration domains) of the input and output containers must agree!");
212  TEUCHOS_TEST_FOR_EXCEPTION( (outputFields.dimension(1) != leftFields.dimension(1) ), std::invalid_argument,
213  ">>> ERROR (ArrayTools::contractFieldFieldTensor): First dimension of output container and first dimension of left input container must agree!");
214  TEUCHOS_TEST_FOR_EXCEPTION( (outputFields.dimension(2) != rightFields.dimension(1) ), std::invalid_argument,
215  ">>> ERROR (ArrayTools::contractFieldFieldTensor): Second dimension of output container and first dimension of right input container must agree!");
216  TEUCHOS_TEST_FOR_EXCEPTION( ( compEngine!=COMP_CPP && compEngine!=COMP_BLAS ), std::invalid_argument,
217  ">>> ERROR (ArrayTools::contractFieldFieldScalar): Computational engine not defined!");
218 
219 #endif
220 
221  // get sizes
222  int numCells = leftFields.dimension(0);
223  int numLeftFields = leftFields.dimension(1);
224  int numRightFields = rightFields.dimension(1);
225  int numPoints = leftFields.dimension(2);
226  int dim1Tensor = leftFields.dimension(3);
227  int dim2Tensor = leftFields.dimension(4);
228 
229  if (sumInto) {
230  for (int cl = 0; cl < numCells; cl++) {
231  for (int lbf = 0; lbf < numLeftFields; lbf++) {
232  for (int rbf = 0; rbf < numRightFields; rbf++) {
233  Scalar tmpVal(0);
234  for (int qp = 0; qp < numPoints; qp++) {
235  for (int iTens1 = 0; iTens1 < dim1Tensor; iTens1++) {
236  for (int iTens2 = 0; iTens2 < dim2Tensor; iTens2++) {
237  tmpVal += leftFields(cl, lbf, qp, iTens1, iTens2)*rightFields(cl, rbf, qp, iTens1, iTens2);
238  } // D2-loop
239  } // D1-loop
240  } // P-loop
241  outputFields(cl, lbf, rbf) += tmpVal;
242  } // R-loop
243  } // L-loop
244  } // C-loop
245  }
246  else {
247  for (int cl = 0; cl < numCells; cl++) {
248  for (int lbf = 0; lbf < numLeftFields; lbf++) {
249  for (int rbf = 0; rbf < numRightFields; rbf++) {
250  Scalar tmpVal(0);
251  for (int qp = 0; qp < numPoints; qp++) {
252  for (int iTens1 = 0; iTens1 < dim1Tensor; iTens1++) {
253  for (int iTens2 = 0; iTens2 < dim2Tensor; iTens2++) {
254  tmpVal += leftFields(cl, lbf, qp, iTens1, iTens2)*rightFields(cl, rbf, qp, iTens1, iTens2);
255  } // D2-loop
256  } // D1-loop
257  } // P-loop
258  outputFields(cl, lbf, rbf) = tmpVal;
259  } // R-loop
260  } // L-loop
261  } // C-loop
262  }
263 } // end contractFieldFieldTensor
264 
265 
266 template<class Scalar, class ArrayOutFields, class ArrayInData, class ArrayInFields>
267 void ArrayTools::contractDataFieldScalar(ArrayOutFields & outputFields,
268  const ArrayInData & inputData,
269  const ArrayInFields & inputFields,
270  const ECompEngine compEngine,
271  const bool sumInto) {
272 
273 
274 #ifdef HAVE_INTREPID_DEBUG
275  TEUCHOS_TEST_FOR_EXCEPTION( (getrank(inputFields) != 3 ), std::invalid_argument,
276  ">>> ERROR (ArrayTools::contractDataFieldScalar): Rank of the fields input argument must equal 3!");
277  TEUCHOS_TEST_FOR_EXCEPTION( (getrank(inputData) != 2 ), std::invalid_argument,
278  ">>> ERROR (ArrayTools::contractDataFieldScalar): Rank of the data input argument must equal 2!");
279  TEUCHOS_TEST_FOR_EXCEPTION( (getrank(outputFields) != 2 ), std::invalid_argument,
280  ">>> ERROR (ArrayTools::contractDataFieldScalar): Rank of output argument must equal 2!");
281  TEUCHOS_TEST_FOR_EXCEPTION( (inputFields.dimension(0) != inputData.dimension(0) ), std::invalid_argument,
282  ">>> ERROR (ArrayTools::contractDataFieldScalar): Zeroth dimensions (number of integration domains) of the fields and data input containers must agree!");
283  TEUCHOS_TEST_FOR_EXCEPTION( ( (inputFields.dimension(2) != inputData.dimension(1)) && (inputData.dimension(1) != 1) ), std::invalid_argument,
284  ">>> ERROR (ArrayTools::contractDataFieldScalar): Second dimension of fields input container and first dimension of data input container (number of integration points) must agree or first data dimension must be 1!");
285  TEUCHOS_TEST_FOR_EXCEPTION( (outputFields.dimension(0) != inputFields.dimension(0) ), std::invalid_argument,
286  ">>> ERROR (ArrayTools::contractDataFieldScalar): Zeroth dimensions (numbers of integration domains) of the fields input and output containers must agree!");
287  TEUCHOS_TEST_FOR_EXCEPTION( (outputFields.dimension(1) != inputFields.dimension(1) ), std::invalid_argument,
288  ">>> ERROR (ArrayTools::contractDataFieldScalar): First dimensions (number of fields) of the fields input and output containers must agree!");
289  TEUCHOS_TEST_FOR_EXCEPTION( ( compEngine!=COMP_CPP && compEngine!=COMP_BLAS ), std::invalid_argument,
290  ">>> ERROR (ArrayTools::contractFieldFieldScalar): Computational engine not defined!");
291 
292 #endif
293 
294  // get sizes
295  int numCells = inputFields.dimension(0);
296  int numFields = inputFields.dimension(1);
297  int numPoints = inputFields.dimension(2);
298  int numDataPoints = inputData.dimension(1);
299 
300 
301  if (sumInto) {
302  if (numDataPoints != 1) { // nonconstant data
303  for (int cl = 0; cl < numCells; cl++) {
304  for (int lbf = 0; lbf < numFields; lbf++) {
305  Scalar tmpVal(0);
306  for (int qp = 0; qp < numPoints; qp++) {
307  tmpVal += inputFields(cl, lbf, qp)*inputData(cl, qp);
308  } // P-loop
309  outputFields(cl, lbf) += tmpVal;
310  } // F-loop
311  } // C-loop
312  }
313  else { // constant data
314  for (int cl = 0; cl < numCells; cl++) {
315  for (int lbf = 0; lbf < numFields; lbf++) {
316  Scalar tmpVal(0);
317  for (int qp = 0; qp < numPoints; qp++) {
318  tmpVal += inputFields(cl, lbf, qp)*inputData(cl, 0);
319  } // P-loop
320  outputFields(cl, lbf) += tmpVal;
321  } // F-loop
322  } // C-loop
323  } // numDataPoints
324  }
325  else {
326  if (numDataPoints != 1) { // nonconstant data
327  for (int cl = 0; cl < numCells; cl++) {
328  for (int lbf = 0; lbf < numFields; lbf++) {
329  Scalar tmpVal(0);
330  for (int qp = 0; qp < numPoints; qp++) {
331  tmpVal += inputFields(cl, lbf, qp)*inputData(cl, qp);
332  } // P-loop
333  outputFields(cl, lbf) = tmpVal;
334  } // F-loop
335  } // C-loop
336  }
337  else { // constant data
338  for (int cl = 0; cl < numCells; cl++) {
339  for (int lbf = 0; lbf < numFields; lbf++) {
340  Scalar tmpVal(0);
341  for (int qp = 0; qp < numPoints; qp++) {
342  tmpVal += inputFields(cl, lbf, qp)*inputData(cl, 0);
343  } // P-loop
344  outputFields(cl, lbf) = tmpVal;
345  } // F-loop
346  } // C-loop
347  } // numDataPoints
348  }
349 
350 } // end contractDataFieldScalar
351 
352 
353 template<class Scalar, class ArrayOutFields, class ArrayInData, class ArrayInFields>
354 void ArrayTools::contractDataFieldVector(ArrayOutFields & outputFields,
355  const ArrayInData & inputData,
356  const ArrayInFields & inputFields,
357  const ECompEngine compEngine,
358  const bool sumInto) {
359 
360 #ifdef HAVE_INTREPID_DEBUG
361  TEUCHOS_TEST_FOR_EXCEPTION( (getrank(inputFields) != 4 ), std::invalid_argument,
362  ">>> ERROR (ArrayTools::contractDataFieldVector): Rank of the fields input argument must equal 4!");
363  TEUCHOS_TEST_FOR_EXCEPTION( (getrank(inputData) != 3 ), std::invalid_argument,
364  ">>> ERROR (ArrayTools::contractDataFieldVector): Rank of the data input argument must equal 3!");
365  TEUCHOS_TEST_FOR_EXCEPTION( (getrank(outputFields) != 2 ), std::invalid_argument,
366  ">>> ERROR (ArrayTools::contractDataFieldVector): Rank of output argument must equal 2!");
367  TEUCHOS_TEST_FOR_EXCEPTION( (inputFields.dimension(0) != inputData.dimension(0) ), std::invalid_argument,
368  ">>> ERROR (ArrayTools::contractDataFieldVector): Zeroth dimensions (number of integration domains) of the fields and data input containers must agree!");
369  TEUCHOS_TEST_FOR_EXCEPTION( ( (inputFields.dimension(2) != inputData.dimension(1)) && (inputData.dimension(1) != 1) ), std::invalid_argument,
370  ">>> ERROR (ArrayTools::contractDataFieldVector): Second dimension of the fields input container and first dimension of data input container (number of integration points) must agree or first data dimension must be 1!");
371  TEUCHOS_TEST_FOR_EXCEPTION( (inputFields.dimension(3) != inputData.dimension(2) ), std::invalid_argument,
372  ">>> ERROR (ArrayTools::contractDataFieldVector): Third dimension of the fields input container and second dimension of data input container (vector index) must agree!");
373  TEUCHOS_TEST_FOR_EXCEPTION( (outputFields.dimension(0) != inputFields.dimension(0) ), std::invalid_argument,
374  ">>> ERROR (ArrayTools::contractDataFieldVector): Zeroth dimensions (numbers of integration domains) of the fields input and output containers must agree!");
375  TEUCHOS_TEST_FOR_EXCEPTION( (outputFields.dimension(1) != inputFields.dimension(1) ), std::invalid_argument,
376  ">>> ERROR (ArrayTools::contractDataFieldVector): First dimensions of output container and fields input container (number of fields) must agree!");
377  TEUCHOS_TEST_FOR_EXCEPTION( ( compEngine!=COMP_CPP && compEngine!=COMP_BLAS ), std::invalid_argument,
378  ">>> ERROR (ArrayTools::contractFieldFieldScalar): Computational engine not defined!");
379 
380 
381 #endif
382 
383  // get sizes
384  int numCells = inputFields.dimension(0);
385  int numFields = inputFields.dimension(1);
386  int numPoints = inputFields.dimension(2);
387  int dimVec = inputFields.dimension(3);
388  int numDataPoints = inputData.dimension(1);
389 
390  if (sumInto) {
391  if (numDataPoints != 1) { // nonconstant data
392  for (int cl = 0; cl < numCells; cl++) {
393  for (int lbf = 0; lbf < numFields; lbf++) {
394  Scalar tmpVal(0);
395  for (int qp = 0; qp < numPoints; qp++) {
396  for (int iVec = 0; iVec < dimVec; iVec++) {
397  tmpVal += inputFields(cl, lbf, qp, iVec)*inputData(cl, qp, iVec);
398  } // D-loop
399  } // P-loop
400  outputFields(cl, lbf) += tmpVal;
401  } // F-loop
402  } // C-loop
403  }
404  else { // constant data
405  for (int cl = 0; cl < numCells; cl++) {
406  for (int lbf = 0; lbf < numFields; lbf++) {
407  Scalar tmpVal(0);
408  for (int qp = 0; qp < numPoints; qp++) {
409  for (int iVec = 0; iVec < dimVec; iVec++) {
410  tmpVal += inputFields(cl, lbf, qp, iVec)*inputData(cl, 0, iVec);
411  } //D-loop
412  } // P-loop
413  outputFields(cl, lbf) += tmpVal;
414  } // F-loop
415  } // C-loop
416  } // numDataPoints
417  }
418  else {
419  if (numDataPoints != 1) { // nonconstant data
420  for (int cl = 0; cl < numCells; cl++) {
421  for (int lbf = 0; lbf < numFields; lbf++) {
422  Scalar tmpVal(0);
423  for (int qp = 0; qp < numPoints; qp++) {
424  for (int iVec = 0; iVec < dimVec; iVec++) {
425  tmpVal += inputFields(cl, lbf, qp, iVec)*inputData(cl, qp, iVec);
426  } // D-loop
427  } // P-loop
428  outputFields(cl, lbf) = tmpVal;
429  } // F-loop
430  } // C-loop
431  }
432  else { // constant data
433  for (int cl = 0; cl < numCells; cl++) {
434  for (int lbf = 0; lbf < numFields; lbf++) {
435  Scalar tmpVal(0);
436  for (int qp = 0; qp < numPoints; qp++) {
437  for (int iVec = 0; iVec < dimVec; iVec++) {
438  tmpVal += inputFields(cl, lbf, qp, iVec)*inputData(cl, 0, iVec);
439  } //D-loop
440  } // P-loop
441  outputFields(cl, lbf) = tmpVal;
442  } // F-loop
443  } // C-loop
444  } // numDataPoints
445  }
446  } // end contractDataFieldVector
447 
448 
449 template<class Scalar, class ArrayOutFields, class ArrayInData, class ArrayInFields>
450 void ArrayTools::contractDataFieldTensor(ArrayOutFields & outputFields,
451  const ArrayInData & inputData,
452  const ArrayInFields & inputFields,
453  const ECompEngine compEngine,
454  const bool sumInto) {
455 
456 #ifdef HAVE_INTREPID_DEBUG
457  TEUCHOS_TEST_FOR_EXCEPTION( (getrank(inputFields) != 5 ), std::invalid_argument,
458  ">>> ERROR (ArrayTools::contractDataFieldTensor): Rank of the fields input argument must equal 5!");
459  TEUCHOS_TEST_FOR_EXCEPTION( (getrank(inputData) != 4 ), std::invalid_argument,
460  ">>> ERROR (ArrayTools::contractDataFieldTensor): Rank of the data input argument must equal 4!");
461  TEUCHOS_TEST_FOR_EXCEPTION( (getrank(outputFields) != 2 ), std::invalid_argument,
462  ">>> ERROR (ArrayTools::contractDataFieldTensor): Rank of output argument must equal 2!");
463  TEUCHOS_TEST_FOR_EXCEPTION( (inputFields.dimension(0) != inputData.dimension(0) ), std::invalid_argument,
464  ">>> ERROR (ArrayTools::contractDataFieldTensor): Zeroth dimensions (number of integration domains) of the fields and data input containers must agree!");
465  TEUCHOS_TEST_FOR_EXCEPTION( ( (inputFields.dimension(2) != inputData.dimension(1)) && (inputData.dimension(1) != 1) ), std::invalid_argument,
466  ">>> ERROR (ArrayTools::contractDataFieldTensor): Second dimension of the fields input container and first dimension of data input container (number of integration points) must agree or first data dimension must be 1!");
467  TEUCHOS_TEST_FOR_EXCEPTION( (inputFields.dimension(3) != inputData.dimension(2) ), std::invalid_argument,
468  ">>> ERROR (ArrayTools::contractDataFieldTensor): Third dimension of the fields input container and second dimension of data input container (first tensor dimension) must agree!");
469  TEUCHOS_TEST_FOR_EXCEPTION( (inputFields.dimension(4) != inputData.dimension(3) ), std::invalid_argument,
470  ">>> ERROR (ArrayTools::contractDataFieldTensor): Fourth dimension of the fields input container and third dimension of data input container (second tensor dimension) must agree!");
471  TEUCHOS_TEST_FOR_EXCEPTION( (outputFields.dimension(0) != inputFields.dimension(0) ), std::invalid_argument,
472  ">>> ERROR (ArrayTools::contractDataFieldTensor): Zeroth dimensions (numbers of integration domains) of the fields input and output containers must agree!");
473  TEUCHOS_TEST_FOR_EXCEPTION( (outputFields.dimension(1) != inputFields.dimension(1) ), std::invalid_argument,
474  ">>> ERROR (ArrayTools::contractDataFieldTensor): First dimensions (number of fields) of output container and fields input container must agree!");
475  TEUCHOS_TEST_FOR_EXCEPTION( ( compEngine!=COMP_CPP && compEngine!=COMP_BLAS ), std::invalid_argument,
476  ">>> ERROR (ArrayTools::contractFieldFieldScalar): Computational engine not defined!");
477 
478 #endif
479 
480  // get sizes
481  int numCells = inputFields.dimension(0);
482  int numFields = inputFields.dimension(1);
483  int numPoints = inputFields.dimension(2);
484  int dim1Tens = inputFields.dimension(3);
485  int dim2Tens = inputFields.dimension(4);
486  int numDataPoints = inputData.dimension(1);
487 
488 
489  if (sumInto) {
490  if (numDataPoints != 1) { // nonconstant data
491  for (int cl = 0; cl < numCells; cl++) {
492  for (int lbf = 0; lbf < numFields; lbf++) {
493  Scalar tmpVal(0);
494  for (int qp = 0; qp < numPoints; qp++) {
495  for (int iTens1 = 0; iTens1 < dim1Tens; iTens1++) {
496  for (int iTens2 =0; iTens2 < dim2Tens; iTens2++) {
497  tmpVal += inputFields(cl, lbf, qp, iTens1, iTens2)*inputData(cl, qp, iTens1, iTens2);
498  } // D2-loop
499  } // D1-loop
500  } // P-loop
501  outputFields(cl, lbf) += tmpVal;
502  } // F-loop
503  } // C-loop
504  }
505  else { // constant data
506  for (int cl = 0; cl < numCells; cl++) {
507  for (int lbf = 0; lbf < numFields; lbf++) {
508  Scalar tmpVal(0);
509  for (int qp = 0; qp < numPoints; qp++) {
510  for (int iTens1 = 0; iTens1 < dim1Tens; iTens1++) {
511  for (int iTens2 = 0; iTens2 < dim2Tens; iTens2++) {
512  tmpVal += inputFields(cl, lbf, qp, iTens1, iTens2)*inputData(cl, 0, iTens1, iTens2);
513  } // D2-loop
514  } // D1-loop
515  } // P-loop
516  outputFields(cl, lbf) += tmpVal;
517  } // F-loop
518  } // C-loop
519  } // numDataPoints
520  }
521  else {
522  if (numDataPoints != 1) { // nonconstant data
523  for (int cl = 0; cl < numCells; cl++) {
524  for (int lbf = 0; lbf < numFields; lbf++) {
525  Scalar tmpVal(0);
526  for (int qp = 0; qp < numPoints; qp++) {
527  for (int iTens1 = 0; iTens1 < dim1Tens; iTens1++) {
528  for (int iTens2 =0; iTens2 < dim2Tens; iTens2++) {
529  tmpVal += inputFields(cl, lbf, qp, iTens1, iTens2)*inputData(cl, qp, iTens1, iTens2);
530  } // D2-loop
531  } // D1-loop
532  } // P-loop
533  outputFields(cl, lbf) = tmpVal;
534  } // F-loop
535  } // C-loop
536  }
537  else { // constant data
538  for (int cl = 0; cl < numCells; cl++) {
539  for (int lbf = 0; lbf < numFields; lbf++) {
540  Scalar tmpVal(0);
541  for (int qp = 0; qp < numPoints; qp++) {
542  for (int iTens1 = 0; iTens1 < dim1Tens; iTens1++) {
543  for (int iTens2 = 0; iTens2 < dim2Tens; iTens2++) {
544  tmpVal += inputFields(cl, lbf, qp, iTens1, iTens2)*inputData(cl, 0, iTens1, iTens2);
545  } // D2-loop
546  } // D1-loop
547  } // P-loop
548  outputFields(cl, lbf) = tmpVal;
549  } // F-loop
550  } // C-loop
551  } // numDataPoints
552  }
553 } // end contractDataFieldTensor
554 
555 
556 template<class Scalar, class ArrayOutData, class ArrayInDataLeft, class ArrayInDataRight>
557 void ArrayTools::contractDataDataScalar(ArrayOutData & outputData,
558  const ArrayInDataLeft & inputDataLeft,
559  const ArrayInDataRight & inputDataRight,
560  const ECompEngine compEngine,
561  const bool sumInto) {
562 #ifdef HAVE_INTREPID_DEBUG
563  TEUCHOS_TEST_FOR_EXCEPTION( (getrank(inputDataLeft) != 2 ), std::invalid_argument,
564  ">>> ERROR (ArrayTools::contractDataDataScalar): Rank of the left input argument must equal 2!");
565  TEUCHOS_TEST_FOR_EXCEPTION( (getrank(inputDataRight) != 2 ), std::invalid_argument,
566  ">>> ERROR (ArrayTools::contractDataDataScalar): Rank of right input argument must equal 2!");
567  TEUCHOS_TEST_FOR_EXCEPTION( (getrank(outputData) != 1 ), std::invalid_argument,
568  ">>> ERROR (ArrayTools::contractDataDataScalar): Rank of output argument must equal 1!");
569  TEUCHOS_TEST_FOR_EXCEPTION( (inputDataLeft.dimension(0) != inputDataRight.dimension(0) ), std::invalid_argument,
570  ">>> ERROR (ArrayTools::contractDataDataScalar): Zeroth dimensions (number of integration domains) of the left and right input containers must agree!");
571  TEUCHOS_TEST_FOR_EXCEPTION( (inputDataLeft.dimension(1) != inputDataRight.dimension(1) ), std::invalid_argument,
572  ">>> ERROR (ArrayTools::contractDataDataScalar): First dimensions (numbers of integration points) of the left and right input containers must agree!");
573  TEUCHOS_TEST_FOR_EXCEPTION( (outputData.dimension(0) != inputDataRight.dimension(0) ), std::invalid_argument,
574  ">>> ERROR (ArrayTools::contractDataDataScalar): Zeroth dimensions (numbers of integration domains) of the input and output containers must agree!");
575  TEUCHOS_TEST_FOR_EXCEPTION( ( compEngine!=COMP_CPP && compEngine!=COMP_BLAS ), std::invalid_argument,
576  ">>> ERROR (ArrayTools::contractFieldFieldScalar): Computational engine not defined!");
577 
578 #endif
579  // get sizes
580  int numCells = inputDataLeft.dimension(0);
581  int numPoints = inputDataLeft.dimension(1);
582 
583  if (sumInto) {
584  for (int cl = 0; cl < numCells; cl++) {
585  Scalar tmpVal(0);
586  for (int qp = 0; qp < numPoints; qp++) {
587  tmpVal += inputDataLeft(cl, qp)*inputDataRight(cl, qp);
588  } // P-loop
589  outputData(cl) += tmpVal;
590  } // C-loop
591  }
592  else {
593  for (int cl = 0; cl < numCells; cl++) {
594  Scalar tmpVal(0);
595  for (int qp = 0; qp < numPoints; qp++) {
596  tmpVal += inputDataLeft(cl, qp)*inputDataRight(cl, qp);
597  } // P-loop
598  outputData(cl) = tmpVal;
599  } // C-loop
600  }
601 
602 
603 } // end contractDataDataScalar
604 
605 
606 template<class Scalar, class ArrayOutData, class ArrayInDataLeft, class ArrayInDataRight>
607 void ArrayTools::contractDataDataVector(ArrayOutData & outputData,
608  const ArrayInDataLeft & inputDataLeft,
609  const ArrayInDataRight & inputDataRight,
610  const ECompEngine compEngine,
611  const bool sumInto) {
612 
613 #ifdef HAVE_INTREPID_DEBUG
614  TEUCHOS_TEST_FOR_EXCEPTION( (getrank(inputDataLeft) != 3 ), std::invalid_argument,
615  ">>> ERROR (ArrayTools::contractDataDataVector): Rank of the left input argument must equal 3!");
616  TEUCHOS_TEST_FOR_EXCEPTION( (getrank(inputDataRight) != 3 ), std::invalid_argument,
617  ">>> ERROR (ArrayTools::contractDataDataVector): Rank of right input argument must equal 3!");
618  TEUCHOS_TEST_FOR_EXCEPTION( (getrank(outputData) != 1 ), std::invalid_argument,
619  ">>> ERROR (ArrayTools::contractDataDataVector): Rank of output argument must equal 1!");
620  TEUCHOS_TEST_FOR_EXCEPTION( (inputDataLeft.dimension(0) != inputDataRight.dimension(0) ), std::invalid_argument,
621  ">>> ERROR (ArrayTools::contractDataDataVector): Zeroth dimensions (number of integration domains) of the left and right input containers must agree!");
622  TEUCHOS_TEST_FOR_EXCEPTION( (inputDataLeft.dimension(1) != inputDataRight.dimension(1) ), std::invalid_argument,
623  ">>> ERROR (ArrayTools::contractDataDataVector): First dimensions (numbers of integration points) of the left and right input containers must agree!");
624  TEUCHOS_TEST_FOR_EXCEPTION( (inputDataLeft.dimension(2) != inputDataRight.dimension(2) ), std::invalid_argument,
625  ">>> ERROR (ArrayTools::contractDataDataVector): Second dimensions (numbers of vector components) of the left and right input containers must agree!");
626  TEUCHOS_TEST_FOR_EXCEPTION( (outputData.dimension(0) != inputDataRight.dimension(0) ), std::invalid_argument,
627  ">>> ERROR (ArrayTools::contractDataDataVector): Zeroth dimensions (numbers of integration domains) of the input and output containers must agree!");
628  TEUCHOS_TEST_FOR_EXCEPTION( ( compEngine!=COMP_CPP && compEngine!=COMP_BLAS ), std::invalid_argument,
629  ">>> ERROR (ArrayTools::contractFieldFieldScalar): Computational engine not defined!");
630 
631 #endif
632 
633  // get sizes
634  int numCells = inputDataLeft.dimension(0);
635  int numPoints = inputDataLeft.dimension(1);
636  int dimVec = inputDataLeft.dimension(2);
637 
638 
639  if (sumInto) {
640  for (int cl = 0; cl < numCells; cl++) {
641  Scalar tmpVal(0);
642  for (int qp = 0; qp < numPoints; qp++) {
643  for (int iVec = 0; iVec < dimVec; iVec++) {
644  tmpVal += inputDataLeft(cl, qp, iVec)*inputDataRight(cl, qp, iVec);
645  } // D-loop
646  } // P-loop
647  outputData(cl) += tmpVal;
648  } // C-loop
649  }
650  else {
651  for (int cl = 0; cl < numCells; cl++) {
652  Scalar tmpVal(0);
653  for (int qp = 0; qp < numPoints; qp++) {
654  for (int iVec = 0; iVec < dimVec; iVec++) {
655  tmpVal += inputDataLeft(cl, qp, iVec)*inputDataRight(cl, qp, iVec);
656  } // D-loop
657  } // P-loop
658  outputData(cl) = tmpVal;
659  } // C-loop
660  }
661  } // end contractDataDataVector
662 
663 
664 template<class Scalar, class ArrayOutData, class ArrayInDataLeft, class ArrayInDataRight>
665 void ArrayTools::contractDataDataTensor(ArrayOutData & outputData,
666  const ArrayInDataLeft & inputDataLeft,
667  const ArrayInDataRight & inputDataRight,
668  const ECompEngine compEngine,
669  const bool sumInto) {
670 
671 #ifdef HAVE_INTREPID_DEBUG
672  TEUCHOS_TEST_FOR_EXCEPTION( (getrank(inputDataLeft) != 4 ), std::invalid_argument,
673  ">>> ERROR (ArrayTools::contractDataDataTensor): Rank of the left input argument must equal 4");
674  TEUCHOS_TEST_FOR_EXCEPTION( (getrank(inputDataRight) != 4 ), std::invalid_argument,
675  ">>> ERROR (ArrayTools::contractDataDataTensor): Rank of right input argument must equal 4!");
676  TEUCHOS_TEST_FOR_EXCEPTION( (getrank(outputData) != 1 ), std::invalid_argument,
677  ">>> ERROR (ArrayTools::contractDataDataTensor): Rank of output argument must equal 1!");
678  TEUCHOS_TEST_FOR_EXCEPTION( (inputDataLeft.dimension(0) != inputDataRight.dimension(0) ), std::invalid_argument,
679  ">>> ERROR (ArrayTools::contractDataDataTensor): Zeroth dimensions (number of integration domains) of the left and right input containers must agree!");
680  TEUCHOS_TEST_FOR_EXCEPTION( (inputDataLeft.dimension(1) != inputDataRight.dimension(1) ), std::invalid_argument,
681  ">>> ERROR (ArrayTools::contractDataDataTensor): First dimensions (numbers of integration points) of the left and right input containers must agree!");
682  TEUCHOS_TEST_FOR_EXCEPTION( (inputDataLeft.dimension(2) != inputDataRight.dimension(2) ), std::invalid_argument,
683  ">>> ERROR (ArrayTools::contractDataDataTensor): Second dimensions (first tensor dimensions) of the left and right input containers must agree!");
684  TEUCHOS_TEST_FOR_EXCEPTION( (inputDataLeft.dimension(3) != inputDataRight.dimension(3) ), std::invalid_argument,
685  ">>> ERROR (ArrayTools::contractDataDataTensor): Third dimensions (second tensor dimensions) of the left and right input containers must agree!");
686  TEUCHOS_TEST_FOR_EXCEPTION( (outputData.dimension(0) != inputDataRight.dimension(0) ), std::invalid_argument,
687  ">>> ERROR (ArrayTools::contractDataDataTensor): Zeroth dimensions (numbers of integration domains) of the input and output containers must agree!");
688  TEUCHOS_TEST_FOR_EXCEPTION( ( compEngine!=COMP_CPP && compEngine!=COMP_BLAS ), std::invalid_argument,
689  ">>> ERROR (ArrayTools::contractFieldFieldScalar): Computational engine not defined!");
690 
691 #endif
692 
693  // get sizes
694  int numCells = inputDataLeft.dimension(0);
695  int numPoints = inputDataLeft.dimension(1);
696  int dim1Tensor = inputDataLeft.dimension(2);
697  int dim2Tensor = inputDataLeft.dimension(3);
698 
699 
700  if (sumInto) {
701  for (int cl = 0; cl < numCells; cl++) {
702  Scalar tmpVal(0);
703  for (int qp = 0; qp < numPoints; qp++) {
704  for (int iTens1 = 0; iTens1 < dim1Tensor; iTens1++) {
705  for (int iTens2 = 0; iTens2 < dim2Tensor; iTens2++) {
706  tmpVal += inputDataLeft(cl, qp, iTens1, iTens2)*inputDataRight(cl, qp, iTens1, iTens2);
707  } // D2-loop
708  } // D1-loop
709  } // P-loop
710  outputData(cl) += tmpVal;
711  } // C-loop
712  }
713  else {
714  for (int cl = 0; cl < numCells; cl++) {
715  Scalar tmpVal(0);
716  for (int qp = 0; qp < numPoints; qp++) {
717  for (int iTens1 = 0; iTens1 < dim1Tensor; iTens1++) {
718  for (int iTens2 = 0; iTens2 < dim2Tensor; iTens2++) {
719  tmpVal += inputDataLeft(cl, qp, iTens1, iTens2)*inputDataRight(cl, qp, iTens1, iTens2);
720  } // D2-loop
721  } // D1-loop
722  } // P-loop
723  outputData(cl) = tmpVal;
724  } // C-loop
725  }
726 
727  } // end contractDataDataTensor
728 
729 
730 }
static void contractDataDataVector(ArrayOutData &outputData, const ArrayInDataLeft &inputDataLeft, const ArrayInDataRight &inputDataRight, const ECompEngine compEngine, const bool sumInto=false)
Contracts the &quot;point&quot; and &quot;space&quot; dimensions P and D of rank-3 containers with dimensions (C...
static void contractDataDataScalar(ArrayOutData &outputData, const ArrayInDataLeft &inputDataLeft, const ArrayInDataRight &inputDataRight, const ECompEngine compEngine, const bool sumInto=false)
Contracts the &quot;point&quot; dimensions P of rank-2 containers with dimensions (C,P), and returns the result...
static void contractDataDataTensor(ArrayOutData &outputData, const ArrayInDataLeft &inputDataLeft, const ArrayInDataRight &inputDataRight, const ECompEngine compEngine, const bool sumInto=false)
Contracts the &quot;point&quot; and &quot;space&quot; dimensions P, D1 and D2 of rank-4 containers with dimensions (C...
static void contractFieldFieldTensor(ArrayOutFields &outputFields, const ArrayInFieldsLeft &leftFields, const ArrayInFieldsRight &rightFields, const ECompEngine compEngine, const bool sumInto=false)
Contracts the &quot;point&quot; and &quot;space&quot; dimensions P, D1, and D2 of two rank-5 containers with dimensions (...
static void contractFieldFieldScalar(ArrayOutFields &outputFields, const ArrayInFieldsLeft &leftFields, const ArrayInFieldsRight &rightFields, const ECompEngine compEngine, const bool sumInto=false)
Contracts the &quot;point&quot; dimension P of two rank-3 containers with dimensions (C,L,P) and (C...
static void contractFieldFieldVector(ArrayOutFields &outputFields, const ArrayInFieldsLeft &leftFields, const ArrayInFieldsRight &rightFields, const ECompEngine compEngine, const bool sumInto=false)
Contracts the &quot;point&quot; and &quot;space&quot; dimensions P and D1 of two rank-4 containers with dimensions (C...
static void contractDataFieldVector(ArrayOutFields &outputFields, const ArrayInData &inputData, const ArrayInFields &inputFields, const ECompEngine compEngine, const bool sumInto=false)
Contracts the &quot;point&quot; and &quot;space&quot; dimensions P and D of a rank-4 container and a rank-3 container wit...
static void contractDataFieldScalar(ArrayOutFields &outputFields, const ArrayInData &inputData, const ArrayInFields &inputFields, const ECompEngine compEngine, const bool sumInto=false)
Contracts the &quot;point&quot; dimensions P of a rank-3 containers and a rank-2 container with dimensions (C...
static void contractDataFieldTensor(ArrayOutFields &outputFields, const ArrayInData &inputData, const ArrayInFields &inputFields, const ECompEngine compEngine, const bool sumInto=false)
Contracts the &quot;point&quot; and &quot;space&quot; dimensions P, D1 and D2 of a rank-5 container and a rank-4 containe...