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
49namespace Intrepid {
50template<class Scalar, class ArrayOutFields, class ArrayInFieldsLeft, class ArrayInFieldsRight>
51void 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
116template<class Scalar, class ArrayOutFields, class ArrayInFieldsLeft, class ArrayInFieldsRight>
117void 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
188template<class Scalar, class ArrayOutFields, class ArrayInFieldsLeft, class ArrayInFieldsRight>
189void 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
266template<class Scalar, class ArrayOutFields, class ArrayInData, class ArrayInFields>
267void 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
353template<class Scalar, class ArrayOutFields, class ArrayInData, class ArrayInFields>
354void 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
449template<class Scalar, class ArrayOutFields, class ArrayInData, class ArrayInFields>
450void 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
556template<class Scalar, class ArrayOutData, class ArrayInDataLeft, class ArrayInDataRight>
557void 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
606template<class Scalar, class ArrayOutData, class ArrayInDataLeft, class ArrayInDataRight>
607void 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
664template<class Scalar, class ArrayOutData, class ArrayInDataLeft, class ArrayInDataRight>
665void 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 contractFieldFieldTensor(ArrayOutFields &outputFields, const ArrayInFieldsLeft &leftFields, const ArrayInFieldsRight &rightFields, const ECompEngine compEngine, const bool sumInto=false)
Contracts the "point" and "space" dimensions P, D1, and D2 of two rank-5 containers with dimensions (...
static void contractDataDataScalar(ArrayOutData &outputData, const ArrayInDataLeft &inputDataLeft, const ArrayInDataRight &inputDataRight, const ECompEngine compEngine, const bool sumInto=false)
Contracts the "point" dimensions P of rank-2 containers with dimensions (C,P), and returns the result...
static void contractDataFieldScalar(ArrayOutFields &outputFields, const ArrayInData &inputData, const ArrayInFields &inputFields, const ECompEngine compEngine, const bool sumInto=false)
Contracts the "point" dimensions P of a rank-3 containers and a rank-2 container with dimensions (C,...
static void contractDataDataVector(ArrayOutData &outputData, const ArrayInDataLeft &inputDataLeft, const ArrayInDataRight &inputDataRight, const ECompEngine compEngine, const bool sumInto=false)
Contracts the "point" and "space" dimensions P and D of rank-3 containers with dimensions (C,...
static void contractFieldFieldVector(ArrayOutFields &outputFields, const ArrayInFieldsLeft &leftFields, const ArrayInFieldsRight &rightFields, const ECompEngine compEngine, const bool sumInto=false)
Contracts the "point" and "space" 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 "point" and "space" dimensions P and D of a rank-4 container and a rank-3 container wit...
static void contractDataDataTensor(ArrayOutData &outputData, const ArrayInDataLeft &inputDataLeft, const ArrayInDataRight &inputDataRight, const ECompEngine compEngine, const bool sumInto=false)
Contracts the "point" and "space" dimensions P, D1 and D2 of rank-4 containers with dimensions (C,...
static void contractFieldFieldScalar(ArrayOutFields &outputFields, const ArrayInFieldsLeft &leftFields, const ArrayInFieldsRight &rightFields, const ECompEngine compEngine, const bool sumInto=false)
Contracts the "point" dimension P of two rank-3 containers with dimensions (C,L,P) and (C,...
static void contractDataFieldTensor(ArrayOutFields &outputFields, const ArrayInData &inputData, const ArrayInFields &inputFields, const ECompEngine compEngine, const bool sumInto=false)
Contracts the "point" and "space" dimensions P, D1 and D2 of a rank-5 container and a rank-4 containe...