Intrepid
Intrepid_ArrayToolsDefScalar.hpp
Go to the documentation of this file.
1#ifndef INTREPID_ARRAYTOOLSDEFSCALAR_HPP
2#define INTREPID_ARRAYTOOLSDEFSCALAR_HPP
3// @HEADER
4// ************************************************************************
5//
6// Intrepid Package
7// Copyright (2007) Sandia Corporation
8//
9// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
10// license for use of this work by or on behalf of the U.S. Government.
11//
12// Redistribution and use in source and binary forms, with or without
13// modification, are permitted provided that the following conditions are
14// met:
15//
16// 1. Redistributions of source code must retain the above copyright
17// notice, this list of conditions and the following disclaimer.
18//
19// 2. Redistributions in binary form must reproduce the above copyright
20// notice, this list of conditions and the following disclaimer in the
21// documentation and/or other materials provided with the distribution.
22//
23// 3. Neither the name of the Corporation nor the names of the
24// contributors may be used to endorse or promote products derived from
25// this software without specific prior written permission.
26//
27// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
28// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
29// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
30// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
31// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
32// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
33// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
34// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
35// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
36// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
37// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
38//
39// Questions? Contact Pavel Bochev (pbboche@sandia.gov)
40// Denis Ridzal (dridzal@sandia.gov), or
41// Kara Peterson (kjpeter@sandia.gov)
42//
43// ************************************************************************
44// @HEADER
45
52namespace Intrepid {
53
54template<class Scalar, class ArrayOutFields, class ArrayInData, class ArrayInFields>
55void ArrayTools::scalarMultiplyDataField(ArrayOutFields & outputFields,
56 const ArrayInData & inputData,
57 const ArrayInFields & inputFields,
58 const bool reciprocal) {
59#ifdef HAVE_INTREPID_DEBUG
60 TEUCHOS_TEST_FOR_EXCEPTION( (getrank(inputData) != 2), std::invalid_argument,
61 ">>> ERROR (ArrayTools::scalarMultiplyDataField): Input data container must have rank 2.");
62 if (getrank(outputFields) <= getrank(inputFields)) {
63 TEUCHOS_TEST_FOR_EXCEPTION( ( (getrank(inputFields) < 3) || (getrank(inputFields) > 5) ), std::invalid_argument,
64 ">>> ERROR (ArrayTools::scalarMultiplyDataField): Input fields container must have rank 3, 4, or 5.");
65 TEUCHOS_TEST_FOR_EXCEPTION( (getrank(outputFields) != getrank(inputFields)), std::invalid_argument,
66 ">>> ERROR (ArrayTools::scalarMultiplyDataField): Input and output fields containers must have the same rank.");
67 TEUCHOS_TEST_FOR_EXCEPTION( (inputFields.dimension(0) != inputData.dimension(0) ), std::invalid_argument,
68 ">>> ERROR (ArrayTools::scalarMultiplyDataField): Zeroth dimensions (number of integration domains) of the fields and data input containers must agree!");
69 TEUCHOS_TEST_FOR_EXCEPTION( ( (inputFields.dimension(2) != inputData.dimension(1)) && (inputData.dimension(1) != 1) ), std::invalid_argument,
70 ">>> ERROR (ArrayTools::scalarMultiplyDataField): 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!");
71 for (size_t i=0; i<getrank(inputFields); i++) {
72 std::string errmsg = ">>> ERROR (ArrayTools::scalarMultiplyDataField): Dimension ";
73 errmsg += (char)(48+i);
74 errmsg += " of the input and output fields containers must agree!";
75 TEUCHOS_TEST_FOR_EXCEPTION( (inputFields.dimension(i) != outputFields.dimension(i)), std::invalid_argument, errmsg );
76 }
77 }
78 else {
79 TEUCHOS_TEST_FOR_EXCEPTION( ( (getrank(inputFields) < 2) || (getrank(inputFields) > 4) ), std::invalid_argument,
80 ">>> ERROR (ArrayTools::scalarMultiplyDataField): Input fields container must have rank 2, 3, or 4.");
81 TEUCHOS_TEST_FOR_EXCEPTION( (getrank(outputFields) != getrank(inputFields)+1), std::invalid_argument,
82 ">>> ERROR (ArrayTools::scalarMultiplyDataField): The rank of the input fields container must be one less than the rank of the output fields container.");
83 TEUCHOS_TEST_FOR_EXCEPTION( ( (inputFields.dimension(1) != inputData.dimension(1)) && (inputData.dimension(1) != 1) ), std::invalid_argument,
84 ">>> ERROR (ArrayTools::scalarMultiplyDataField): First dimensions of fields input container and data input container (number of integration points) must agree or first data dimension must be 1!");
85 TEUCHOS_TEST_FOR_EXCEPTION( ( inputData.dimension(0) != outputFields.dimension(0) ), std::invalid_argument,
86 ">>> ERROR (ArrayTools::scalarMultiplyDataField): Zeroth dimensions of fields output container and data input containers (number of integration domains) must agree!");
87 for (size_t i=0; i<getrank(inputFields); i++) {
88 std::string errmsg = ">>> ERROR (ArrayTools::scalarMultiplyDataField): Dimensions ";
89 errmsg += (char)(48+i);
90 errmsg += " and ";
91 errmsg += (char)(48+i+1);
92 errmsg += " of the input and output fields containers must agree!";
93 TEUCHOS_TEST_FOR_EXCEPTION( (inputFields.dimension(i) != outputFields.dimension(i+1)), std::invalid_argument, errmsg );
94 }
95 }
96#endif
97 ArrayWrapper<Scalar,ArrayOutFields, Rank<ArrayOutFields >::value, false>outputFieldsWrap(outputFields);
98 ArrayWrapper<Scalar,ArrayInData, Rank<ArrayInData >::value, true>inputDataWrap(inputData);
99 ArrayWrapper<Scalar,ArrayInFields, Rank<ArrayInFields >::value,true>inputFieldsWrap(inputFields);
100
101 // get sizes
102 size_t invalRank = getrank(inputFields);
103 size_t outvalRank = getrank(outputFields);
104 int numCells = outputFields.dimension(0);
105 int numFields = outputFields.dimension(1);
106 int numPoints = outputFields.dimension(2);
107 int numDataPoints = inputData.dimension(1);
108 int dim1Tens = 0;
109 int dim2Tens = 0;
110 if (outvalRank > 3) {
111 dim1Tens = outputFields.dimension(3);
112 if (outvalRank > 4) {
113 dim2Tens = outputFields.dimension(4);
114 }
115 }
116
117 if (outvalRank == invalRank) {
118
119 if (numDataPoints != 1) { // nonconstant data
120 switch(invalRank) {
121 case 3: {
122 if (reciprocal) {
123 for(int cl = 0; cl < numCells; cl++) {
124 for(int bf = 0; bf < numFields; bf++) {
125 for(int pt = 0; pt < numPoints; pt++) {
126 outputFieldsWrap(cl, bf, pt) = inputFieldsWrap(cl, bf, pt)/inputDataWrap(cl, pt);
127 } // P-loop
128 } // F-loop
129 } // C-loop
130 }
131 else {
132
133
134 for(int cl = 0; cl < numCells; cl++) {
135 for(int bf = 0; bf < numFields; bf++) {
136 for(int pt = 0; pt < numPoints; pt++) {
137 outputFieldsWrap(cl, bf, pt) = inputFieldsWrap(cl, bf, pt)*inputDataWrap(cl, pt);
138 } // P-loop
139 } // F-loop
140 } // C-loop
141
142 }
143 }// case 3
144 break;
145
146 case 4: {
147 if (reciprocal) {
148 for(int cl = 0; cl < numCells; cl++) {
149 for(int bf = 0; bf < numFields; bf++) {
150 for(int pt = 0; pt < numPoints; pt++) {
151 for( int iVec = 0; iVec < dim1Tens; iVec++) {
152 outputFieldsWrap(cl, bf, pt, iVec) = inputFieldsWrap(cl, bf, pt, iVec)/inputDataWrap(cl, pt);
153 } // D1-loop
154 } // P-loop
155 } // F-loop
156 } // C-loop
157 }
158 else {
159 for(int cl = 0; cl < numCells; cl++) {
160 for(int bf = 0; bf < numFields; bf++) {
161 for(int pt = 0; pt < numPoints; pt++) {
162 for( int iVec = 0; iVec < dim1Tens; iVec++) {
163 outputFieldsWrap(cl, bf, pt, iVec) = inputFieldsWrap(cl, bf, pt, iVec)*inputDataWrap(cl, pt);
164 } // D1-loop
165 } // P-loop
166 } // F-loop
167 } // C-loop
168 }
169 }// case 4
170 break;
171
172 case 5: {
173 if (reciprocal) {
174 for(int cl = 0; cl < numCells; cl++) {
175 for(int bf = 0; bf < numFields; bf++) {
176 for(int pt = 0; pt < numPoints; pt++) {
177 for( int iTens1 = 0; iTens1 < dim1Tens; iTens1++) {
178 for( int iTens2 = 0; iTens2 < dim2Tens; iTens2++) {
179 outputFieldsWrap(cl, bf, pt, iTens1, iTens2) = inputFieldsWrap(cl, bf, pt, iTens1, iTens2)/inputDataWrap(cl, pt);
180 } // D2-loop
181 } // D1-loop
182 } // F-loop
183 } // P-loop
184 } // C-loop
185 }
186 else {
187 for(int cl = 0; cl < numCells; cl++) {
188 for(int bf = 0; bf < numFields; bf++) {
189 for(int pt = 0; pt < numPoints; pt++) {
190 for( int iTens1 = 0; iTens1 < dim1Tens; iTens1++) {
191 for( int iTens2 = 0; iTens2 < dim2Tens; iTens2++) {
192 outputFieldsWrap(cl, bf, pt, iTens1, iTens2) = inputFieldsWrap(cl, bf, pt, iTens1, iTens2)*inputDataWrap(cl, pt);
193 } // D2-loop
194 } // D1-loop
195 } // F-loop
196 } // P-loop
197 } // C-loop
198 }
199 }// case 5
200 break;
201
202 default:
203 TEUCHOS_TEST_FOR_EXCEPTION( !( (invalRank == 3) || (invalRank == 4) || (invalRank == 5) ), std::invalid_argument,
204 ">>> ERROR (ArrayTools::scalarMultiplyDataField): This branch of the method is defined only for rank-3,4 or 5 containers.");
205 }// invalRank
206
207 }
208 else { //constant data
209
210 switch(invalRank) {
211 case 3: {
212 if (reciprocal) {
213 for(int cl = 0; cl < numCells; cl++) {
214 for(int bf = 0; bf < numFields; bf++) {
215 for(int pt = 0; pt < numPoints; pt++) {
216 outputFieldsWrap(cl, bf, pt) = inputFieldsWrap(cl, bf, pt)/inputDataWrap(cl, 0);
217 } // P-loop
218 } // F-loop
219 } // C-loop
220 }
221 else {
222 for(int cl = 0; cl < numCells; cl++) {
223 for(int bf = 0; bf < numFields; bf++) {
224 for(int pt = 0; pt < numPoints; pt++) {
225 outputFieldsWrap(cl, bf, pt) = inputFieldsWrap(cl, bf, pt)*inputDataWrap(cl, 0);
226 } // P-loop
227 } // F-loop
228 } // C-loop
229 }
230 }// case 3
231 break;
232
233 case 4: {
234 if (reciprocal) {
235 for(int cl = 0; cl < numCells; cl++) {
236 for(int bf = 0; bf < numFields; bf++) {
237 for(int pt = 0; pt < numPoints; pt++) {
238 for( int iVec = 0; iVec < dim1Tens; iVec++) {
239 outputFieldsWrap(cl, bf, pt, iVec) = inputFieldsWrap(cl, bf, pt, iVec)/inputDataWrap(cl, 0);
240 } // D1-loop
241 } // P-loop
242 } // F-loop
243 } // C-loop
244 }
245 else {
246 for(int cl = 0; cl < numCells; cl++) {
247 for(int bf = 0; bf < numFields; bf++) {
248 for(int pt = 0; pt < numPoints; pt++) {
249 for( int iVec = 0; iVec < dim1Tens; iVec++) {
250 outputFieldsWrap(cl, bf, pt, iVec) = inputFieldsWrap(cl, bf, pt, iVec)*inputDataWrap(cl, 0);
251 } // D1-loop
252 } // P-loop
253 } // F-loop
254 } // C-loop
255 }
256 }// case 4
257 break;
258
259 case 5: {
260 if (reciprocal) {
261 for(int cl = 0; cl < numCells; cl++) {
262 for(int bf = 0; bf < numFields; bf++) {
263 for(int pt = 0; pt < numPoints; pt++) {
264 for( int iTens1 = 0; iTens1 < dim1Tens; iTens1++) {
265 for( int iTens2 = 0; iTens2 < dim2Tens; iTens2++) {
266 outputFieldsWrap(cl, bf, pt, iTens1, iTens2) = inputFieldsWrap(cl, bf, pt, iTens1, iTens2)/inputDataWrap(cl, 0);
267 } // D2-loop
268 } // D1-loop
269 } // F-loop
270 } // P-loop
271 } // C-loop
272 }
273 else {
274 for(int cl = 0; cl < numCells; cl++) {
275 for(int bf = 0; bf < numFields; bf++) {
276 for(int pt = 0; pt < numPoints; pt++) {
277 for( int iTens1 = 0; iTens1 < dim1Tens; iTens1++) {
278 for( int iTens2 = 0; iTens2 < dim2Tens; iTens2++) {
279 outputFieldsWrap(cl, bf, pt, iTens1, iTens2) = inputFieldsWrap(cl, bf, pt, iTens1, iTens2)*inputDataWrap(cl, 0);
280 } // D2-loop
281 } // D1-loop
282 } // F-loop
283 } // P-loop
284 } // C-loop
285 }
286 }// case 5
287 break;
288
289 default:
290 TEUCHOS_TEST_FOR_EXCEPTION( !( (invalRank == 3) || (invalRank == 4) || (invalRank == 5) ), std::invalid_argument,
291 ">>> ERROR (ArrayTools::scalarMultiplyDataField): This branch of the method is defined only for rank-3, 4 or 5 input containers.");
292
293 } // invalRank
294 } // numDataPoints
295
296 }
297 else {
298
299 if (numDataPoints != 1) { // nonconstant data
300
301 switch(invalRank) {
302 case 2: {
303 if (reciprocal) {
304 for(int cl = 0; cl < numCells; cl++) {
305 for(int bf = 0; bf < numFields; bf++) {
306 for(int pt = 0; pt < numPoints; pt++) {
307 outputFieldsWrap(cl, bf, pt) = inputFieldsWrap(bf, pt)/inputDataWrap(cl, pt);
308 } // P-loop
309 } // F-loop
310 } // C-loop
311 }
312 else {
313 for(int cl = 0; cl < numCells; cl++) {
314 for(int bf = 0; bf < numFields; bf++) {
315 for(int pt = 0; pt < numPoints; pt++) {
316 outputFieldsWrap(cl, bf, pt) = inputFieldsWrap(bf, pt)*inputDataWrap(cl, pt);
317 } // P-loop
318 } // F-loop
319 } // C-loop
320 }
321 }// case 2
322 break;
323
324 case 3: {
325 if (reciprocal) {
326 for(int cl = 0; cl < numCells; cl++) {
327 for(int bf = 0; bf < numFields; bf++) {
328 for(int pt = 0; pt < numPoints; pt++) {
329 for( int iVec = 0; iVec < dim1Tens; iVec++) {
330 outputFieldsWrap(cl, bf, pt, iVec) = inputFieldsWrap(bf, pt, iVec)/inputDataWrap(cl, pt);
331 } // D1-loop
332 } // P-loop
333 } // F-loop
334 } // C-loop
335 }
336 else {
337 for(int cl = 0; cl < numCells; cl++) {
338 for(int bf = 0; bf < numFields; bf++) {
339 for(int pt = 0; pt < numPoints; pt++) {
340 for( int iVec = 0; iVec < dim1Tens; iVec++) {
341 outputFieldsWrap(cl, bf, pt, iVec) = inputFieldsWrap(bf, pt, iVec)*inputDataWrap(cl, pt);
342 } // D1-loop
343 } // P-loop
344 } // F-loop
345 } // C-loop
346 }
347 }// case 3
348 break;
349
350 case 4: {
351 if (reciprocal) {
352 for(int cl = 0; cl < numCells; cl++) {
353 for(int bf = 0; bf < numFields; bf++) {
354 for(int pt = 0; pt < numPoints; pt++) {
355 for( int iTens1 = 0; iTens1 < dim1Tens; iTens1++) {
356 for( int iTens2 = 0; iTens2 < dim2Tens; iTens2++) {
357 outputFieldsWrap(cl, bf, pt, iTens1, iTens2) = inputFieldsWrap(bf, pt, iTens1, iTens2)/inputDataWrap(cl, pt);
358 } // D2-loop
359 } // D1-loop
360 } // F-loop
361 } // P-loop
362 } // C-loop
363 }
364 else {
365 for(int cl = 0; cl < numCells; cl++) {
366 for(int bf = 0; bf < numFields; bf++) {
367 for(int pt = 0; pt < numPoints; pt++) {
368 for( int iTens1 = 0; iTens1 < dim1Tens; iTens1++) {
369 for( int iTens2 = 0; iTens2 < dim2Tens; iTens2++) {
370 outputFieldsWrap(cl, bf, pt, iTens1, iTens2) = inputFieldsWrap(bf, pt, iTens1, iTens2)*inputDataWrap(cl, pt);
371 } // D2-loop
372 } // D1-loop
373 } // F-loop
374 } // P-loop
375 } // C-loop
376 }
377 }// case 4
378 break;
379
380 default:
381 TEUCHOS_TEST_FOR_EXCEPTION( !( (invalRank == 2) || (invalRank == 3) || (invalRank == 4) ), std::invalid_argument,
382 ">>> ERROR (ArrayTools::scalarMultiplyDataField): This branch of the method is defined only for rank-2, 3 or 4 input containers.");
383 }// invalRank
384
385 }
386 else { //constant data
387
388 switch(invalRank) {
389 case 2: {
390 if (reciprocal) {
391 for(int cl = 0; cl < numCells; cl++) {
392 for(int bf = 0; bf < numFields; bf++) {
393 for(int pt = 0; pt < numPoints; pt++) {
394 outputFieldsWrap(cl, bf, pt) = inputFieldsWrap(bf, pt)/inputDataWrap(cl, 0);
395 } // P-loop
396 } // F-loop
397 } // C-loop
398 }
399 else {
400 for(int cl = 0; cl < numCells; cl++) {
401 for(int bf = 0; bf < numFields; bf++) {
402 for(int pt = 0; pt < numPoints; pt++) {
403 outputFieldsWrap(cl, bf, pt) = inputFieldsWrap(bf, pt)*inputDataWrap(cl, 0);
404 } // P-loop
405 } // F-loop
406 } // C-loop
407 }
408 }// case 2
409 break;
410
411 case 3: {
412 if (reciprocal) {
413 for(int cl = 0; cl < numCells; cl++) {
414 for(int bf = 0; bf < numFields; bf++) {
415 for(int pt = 0; pt < numPoints; pt++) {
416 for( int iVec = 0; iVec < dim1Tens; iVec++) {
417 outputFieldsWrap(cl, bf, pt, iVec) = inputFieldsWrap(bf, pt, iVec)/inputDataWrap(cl, 0);
418 } // D1-loop
419 } // P-loop
420 } // F-loop
421 } // C-loop
422 }
423 else {
424 for(int cl = 0; cl < numCells; cl++) {
425 for(int bf = 0; bf < numFields; bf++) {
426 for(int pt = 0; pt < numPoints; pt++) {
427 for( int iVec = 0; iVec < dim1Tens; iVec++) {
428 outputFieldsWrap(cl, bf, pt, iVec) = inputFieldsWrap(bf, pt, iVec)*inputDataWrap(cl, 0);
429 } // D1-loop
430 } // P-loop
431 } // F-loop
432 } // C-loop
433 }
434 }// case 3
435 break;
436
437 case 4: {
438 if (reciprocal) {
439 for(int cl = 0; cl < numCells; cl++) {
440 for(int bf = 0; bf < numFields; bf++) {
441 for(int pt = 0; pt < numPoints; pt++) {
442 for( int iTens1 = 0; iTens1 < dim1Tens; iTens1++) {
443 for( int iTens2 = 0; iTens2 < dim2Tens; iTens2++) {
444 outputFieldsWrap(cl, bf, pt, iTens1, iTens2) = inputFieldsWrap(bf, pt, iTens1, iTens2)/inputDataWrap(cl, 0);
445 } // D2-loop
446 } // D1-loop
447 } // F-loop
448 } // P-loop
449 } // C-loop
450 }
451 else {
452 for(int cl = 0; cl < numCells; cl++) {
453 for(int bf = 0; bf < numFields; bf++) {
454 for(int pt = 0; pt < numPoints; pt++) {
455 for( int iTens1 = 0; iTens1 < dim1Tens; iTens1++) {
456 for( int iTens2 = 0; iTens2 < dim2Tens; iTens2++) {
457 outputFieldsWrap(cl, bf, pt, iTens1, iTens2) = inputFieldsWrap(bf, pt, iTens1, iTens2)*inputDataWrap(cl, 0);
458 } // D2-loop
459 } // D1-loop
460 } // F-loop
461 } // P-loop
462 } // C-loop
463 }
464 }// case 4
465 break;
466
467 default:
468 TEUCHOS_TEST_FOR_EXCEPTION( !( (invalRank == 2) || (invalRank == 3) || (invalRank == 3) ), std::invalid_argument,
469 ">>> ERROR (ArrayTools::scalarMultiplyDataField): This branch of the method is defined only for rank-2, 3 or 4 input containers.");
470
471 } // invalRank
472 } // numDataPoints
473
474 } // end if (outvalRank = invalRank)
475
476} // scalarMultiplyDataField
477
478template<class Scalar, class ArrayOutData, class ArrayInDataLeft, class ArrayInDataRight>
479void ArrayTools::scalarMultiplyDataData(ArrayOutData & outputData,
480 const ArrayInDataLeft & inputDataLeft,
481 const ArrayInDataRight & inputDataRight,
482 const bool reciprocal) {
483
484#ifdef HAVE_INTREPID_DEBUG
485 TEUCHOS_TEST_FOR_EXCEPTION( (getrank(inputDataLeft) != 2), std::invalid_argument,
486 ">>> ERROR (ArrayTools::scalarMultiplyDataData): Left input data container must have rank 2.");
487 if (getrank(outputData) <= getrank(inputDataRight)) {
488 TEUCHOS_TEST_FOR_EXCEPTION( ( (getrank(inputDataRight) < 2) || (getrank(inputDataRight) > 4) ), std::invalid_argument,
489 ">>> ERROR (ArrayTools::scalarMultiplyDataData): Right input data container must have rank 2, 3, or 4.");
490 TEUCHOS_TEST_FOR_EXCEPTION( (getrank(outputData) != getrank(inputDataRight)), std::invalid_argument,
491 ">>> ERROR (ArrayTools::scalarMultiplyDataData): Right input and output data containers must have the same rank.");
492 TEUCHOS_TEST_FOR_EXCEPTION( (inputDataRight.dimension(0) != inputDataLeft.dimension(0) ), std::invalid_argument,
493 ">>> ERROR (ArrayTools::scalarMultiplyDataData): Zeroth dimensions (number of integration domains) of the left and right data input containers must agree!");
494 TEUCHOS_TEST_FOR_EXCEPTION( ( (inputDataRight.dimension(1) != inputDataLeft.dimension(1)) && (inputDataLeft.dimension(1) != 1) ), std::invalid_argument,
495 ">>> ERROR (ArrayTools::scalarMultiplyDataData): First dimensions of the left and right data input containers (number of integration points) must agree or first dimension of the left data input container must be 1!");
496 for (size_t i=0; i<getrank(inputDataRight); i++) {
497 std::string errmsg = ">>> ERROR (ArrayTools::scalarMultiplyDataData): Dimension ";
498 errmsg += (char)(48+i);
499 errmsg += " of the right input and output data containers must agree!";
500 TEUCHOS_TEST_FOR_EXCEPTION( (inputDataRight.dimension(i) != outputData.dimension(i)), std::invalid_argument, errmsg );
501 }
502 }
503 else {
504 TEUCHOS_TEST_FOR_EXCEPTION( ( (getrank(inputDataRight) < 1) || (getrank(inputDataRight) > 3) ), std::invalid_argument,
505 ">>> ERROR (ArrayTools::scalarMultiplyDataData): Right input data container must have rank 1, 2, or 3.");
506 TEUCHOS_TEST_FOR_EXCEPTION( (getrank(outputData) != getrank(inputDataRight)+1), std::invalid_argument,
507 ">>> ERROR (ArrayTools::scalarMultiplyDataData): The rank of the right input data container must be one less than the rank of the output data container.");
508 TEUCHOS_TEST_FOR_EXCEPTION( ( (inputDataRight.dimension(0) != inputDataLeft.dimension(1)) && (inputDataLeft.dimension(1) != 1) ), std::invalid_argument,
509 ">>> ERROR (ArrayTools::scalarMultiplyDataData): Zeroth dimension of the right input data container and first dimension of the left data input container (number of integration points) must agree or first dimension of the left data input container must be 1!");
510 TEUCHOS_TEST_FOR_EXCEPTION( ( inputDataLeft.dimension(0) != outputData.dimension(0) ), std::invalid_argument,
511 ">>> ERROR (ArrayTools::scalarMultiplyDataData): Zeroth dimensions of data output and left data input containers (number of integration domains) must agree!");
512 for (size_t i=0; i<getrank(inputDataRight); i++) {
513 std::string errmsg = ">>> ERROR (ArrayTools::scalarMultiplyDataData): Dimensions ";
514 errmsg += (char)(48+i);
515 errmsg += " and ";
516 errmsg += (char)(48+i+1);
517 errmsg += " of the right input and output data containers must agree!";
518 TEUCHOS_TEST_FOR_EXCEPTION( (inputDataRight.dimension(i) != outputData.dimension(i+1)), std::invalid_argument, errmsg );
519 }
520 }
521#endif
522
523
524 ArrayWrapper<Scalar,ArrayOutData, Rank<ArrayOutData >::value, false>outputDataWrap(outputData);
525 ArrayWrapper<Scalar,ArrayInDataLeft, Rank<ArrayInDataLeft >::value, true>inputDataLeftWrap(inputDataLeft);
526 ArrayWrapper<Scalar,ArrayInDataRight, Rank<ArrayInDataRight >::value,true>inputDataRightWrap(inputDataRight);
527
528
529 // get sizes
530 size_t invalRank = getrank(inputDataRight);
531 size_t outvalRank = getrank(outputData);
532 int numCells = outputData.dimension(0);
533 int numPoints = outputData.dimension(1);
534 int numDataPoints = inputDataLeft.dimension(1);
535 int dim1Tens = 0;
536 int dim2Tens = 0;
537 if (outvalRank > 2) {
538 dim1Tens = outputData.dimension(2);
539 if (outvalRank > 3) {
540 dim2Tens = outputData.dimension(3);
541 }
542 }
543
544 if (outvalRank == invalRank) {
545
546 if (numDataPoints != 1) { // nonconstant data
547
548 switch(invalRank) {
549 case 2: {
550 if (reciprocal) {
551 for(int cl = 0; cl < numCells; cl++) {
552 for(int pt = 0; pt < numPoints; pt++) {
553 outputDataWrap(cl, pt) = inputDataRightWrap(cl, pt)/inputDataLeftWrap(cl, pt);
554 } // P-loop
555 } // C-loop
556 }
557 else {
558 for(int cl = 0; cl < numCells; cl++) {
559 for(int pt = 0; pt < numPoints; pt++) {
560 outputDataWrap(cl, pt) = inputDataRightWrap(cl, pt)*inputDataLeftWrap(cl, pt);
561 } // P-loop
562 } // C-loop
563 }
564 }// case 2
565 break;
566
567 case 3: {
568 if (reciprocal) {
569 for(int cl = 0; cl < numCells; cl++) {
570 for(int pt = 0; pt < numPoints; pt++) {
571 for( int iVec = 0; iVec < dim1Tens; iVec++) {
572 outputDataWrap(cl, pt, iVec) = inputDataRightWrap(cl, pt, iVec)/inputDataLeftWrap(cl, pt);
573 } // D1-loop
574 } // P-loop
575 } // C-loop
576 }
577 else {
578 for(int cl = 0; cl < numCells; cl++) {
579 for(int pt = 0; pt < numPoints; pt++) {
580 for( int iVec = 0; iVec < dim1Tens; iVec++) {
581 outputDataWrap(cl, pt, iVec) = inputDataRightWrap(cl, pt, iVec)*inputDataLeftWrap(cl, pt);
582 } // D1-loop
583 } // P-loop
584 } // C-loop
585 }
586 }// case 3
587 break;
588
589 case 4: {
590 if (reciprocal) {
591 for(int cl = 0; cl < numCells; cl++) {
592 for(int pt = 0; pt < numPoints; pt++) {
593 for( int iTens1 = 0; iTens1 < dim1Tens; iTens1++) {
594 for( int iTens2 = 0; iTens2 < dim2Tens; iTens2++) {
595 outputDataWrap(cl, pt, iTens1, iTens2) = inputDataRightWrap(cl, pt, iTens1, iTens2)/inputDataLeftWrap(cl, pt);
596 } // D2-loop
597 } // D1-loop
598 } // P-loop
599 } // C-loop
600 }
601 else {
602 for(int cl = 0; cl < numCells; cl++) {
603 for(int pt = 0; pt < numPoints; pt++) {
604 for( int iTens1 = 0; iTens1 < dim1Tens; iTens1++) {
605 for( int iTens2 = 0; iTens2 < dim2Tens; iTens2++) {
606 outputDataWrap(cl, pt, iTens1, iTens2) = inputDataRightWrap(cl, pt, iTens1, iTens2)*inputDataLeftWrap(cl, pt);
607 } // D2-loop
608 } // D1-loop
609 } // P-loop
610 } // C-loop
611 }
612 }// case 4
613 break;
614
615 default:
616 TEUCHOS_TEST_FOR_EXCEPTION( !( (invalRank == 2) || (invalRank == 3) || (invalRank == 4) ), std::invalid_argument,
617 ">>> ERROR (ArrayTools::scalarMultiplyDataData): This branch of the method is defined only for rank-2, 3 or 4 containers.");
618 }// invalRank
619
620 }
621 else { // constant left data
622
623 switch(invalRank) {
624 case 2: {
625 if (reciprocal) {
626 for(int cl = 0; cl < numCells; cl++) {
627 for(int pt = 0; pt < numPoints; pt++) {
628 outputDataWrap(cl, pt) = inputDataRightWrap(cl, pt)/inputDataLeftWrap(cl, 0);
629 } // P-loop
630 } // C-loop
631 }
632 else {
633 for(int cl = 0; cl < numCells; cl++) {
634 for(int pt = 0; pt < numPoints; pt++) {
635 outputDataWrap(cl, pt) = inputDataRightWrap(cl, pt)*inputDataLeftWrap(cl, 0);
636 } // P-loop
637 } // C-loop
638 }
639 }// case 2
640 break;
641
642 case 3: {
643 if (reciprocal) {
644 for(int cl = 0; cl < numCells; cl++) {
645 for(int pt = 0; pt < numPoints; pt++) {
646 for( int iVec = 0; iVec < dim1Tens; iVec++) {
647 outputDataWrap(cl, pt, iVec) = inputDataRightWrap(cl, pt, iVec)/inputDataLeftWrap(cl, 0);
648 } // D1-loop
649 } // P-loop
650 } // C-loop
651 }
652 else {
653 for(int cl = 0; cl < numCells; cl++) {
654 for(int pt = 0; pt < numPoints; pt++) {
655 for( int iVec = 0; iVec < dim1Tens; iVec++) {
656 outputDataWrap(cl, pt, iVec) = inputDataRightWrap(cl, pt, iVec)*inputDataLeftWrap(cl, 0);
657 } // D1-loop
658 } // P-loop
659 } // C-loop
660 }
661 }// case 3
662 break;
663
664 case 4: {
665 if (reciprocal) {
666 for(int cl = 0; cl < numCells; cl++) {
667 for(int pt = 0; pt < numPoints; pt++) {
668 for( int iTens1 = 0; iTens1 < dim1Tens; iTens1++) {
669 for( int iTens2 = 0; iTens2 < dim2Tens; iTens2++) {
670 outputDataWrap(cl, pt, iTens1, iTens2) = inputDataRightWrap(cl, pt, iTens1, iTens2)/inputDataLeftWrap(cl, 0);
671 } // D2-loop
672 } // D1-loop
673 } // P-loop
674 } // C-loop
675 }
676 else {
677 for(int cl = 0; cl < numCells; cl++) {
678 for(int pt = 0; pt < numPoints; pt++) {
679 for( int iTens1 = 0; iTens1 < dim1Tens; iTens1++) {
680 for( int iTens2 = 0; iTens2 < dim2Tens; iTens2++) {
681 outputDataWrap(cl, pt, iTens1, iTens2) = inputDataRightWrap(cl, pt, iTens1, iTens2)*inputDataLeftWrap(cl, 0);
682 } // D2-loop
683 } // D1-loop
684 } // P-loop
685 } // C-loop
686 }
687 }// case 4
688 break;
689
690 default:
691 TEUCHOS_TEST_FOR_EXCEPTION( !( (invalRank == 2) || (invalRank == 3) || (invalRank == 4) ), std::invalid_argument,
692 ">>> ERROR (ArrayTools::scalarMultiplyDataData): This branch of the method is defined only for rank-2, 3 or 4 input containers.");
693
694 } // invalRank
695 } // numDataPoints
696
697 }
698 else {
699
700 if (numDataPoints != 1) { // nonconstant data
701
702 switch(invalRank) {
703 case 1: {
704 if (reciprocal) {
705 for(int cl = 0; cl < numCells; cl++) {
706 for(int pt = 0; pt < numPoints; pt++) {
707 outputDataWrap(cl, pt) = inputDataRightWrap(pt)/inputDataLeftWrap(cl, pt);
708 } // P-loop
709 } // C-loop
710 }
711 else {
712 for(int cl = 0; cl < numCells; cl++) {
713 for(int pt = 0; pt < numPoints; pt++) {
714 outputDataWrap(cl, pt) = inputDataRightWrap(pt)*inputDataLeftWrap(cl, pt);
715 } // P-loop
716 } // C-loop
717 }
718 }// case 1
719 break;
720
721 case 2: {
722 if (reciprocal) {
723 for(int cl = 0; cl < numCells; cl++) {
724 for(int pt = 0; pt < numPoints; pt++) {
725 for( int iVec = 0; iVec < dim1Tens; iVec++) {
726 outputDataWrap(cl, pt, iVec) = inputDataRightWrap(pt, iVec)/inputDataLeftWrap(cl, pt);
727 } // D1-loop
728 } // P-loop
729 } // C-loop
730 }
731 else {
732 for(int cl = 0; cl < numCells; cl++) {
733 for(int pt = 0; pt < numPoints; pt++) {
734 for( int iVec = 0; iVec < dim1Tens; iVec++) {
735 outputDataWrap(cl, pt, iVec) = inputDataRightWrap(pt, iVec)*inputDataLeftWrap(cl, pt);
736 } // D1-loop
737 } // P-loop
738 } // C-loop
739 }
740 }// case 2
741 break;
742
743 case 3: {
744 if (reciprocal) {
745 for(int cl = 0; cl < numCells; cl++) {
746 for(int pt = 0; pt < numPoints; pt++) {
747 for( int iTens1 = 0; iTens1 < dim1Tens; iTens1++) {
748 for( int iTens2 = 0; iTens2 < dim2Tens; iTens2++) {
749 outputDataWrap(cl, pt, iTens1, iTens2) = inputDataRightWrap(pt, iTens1, iTens2)/inputDataLeftWrap(cl, pt);
750 } // D2-loop
751 } // D1-loop
752 } // P-loop
753 } // C-loop
754 }
755 else {
756 for(int cl = 0; cl < numCells; cl++) {
757 for(int pt = 0; pt < numPoints; pt++) {
758 for( int iTens1 = 0; iTens1 < dim1Tens; iTens1++) {
759 for( int iTens2 = 0; iTens2 < dim2Tens; iTens2++) {
760 outputDataWrap(cl, pt, iTens1, iTens2) = inputDataRightWrap(pt, iTens1, iTens2)*inputDataLeftWrap(cl, pt);
761 } // D2-loop
762 } // D1-loop
763 } // P-loop
764 } // C-loop
765 }
766 }// case 3
767 break;
768
769 default:
770 TEUCHOS_TEST_FOR_EXCEPTION( !( (invalRank == 1) || (invalRank == 2) || (invalRank == 3) ), std::invalid_argument,
771 ">>> ERROR (ArrayTools::scalarMultiplyDataData): This branch of the method is defined only for rank-1, 2 or 3 input containers.");
772 }// invalRank
773
774 }
775 else { //constant data
776
777 switch(invalRank) {
778 case 1: {
779 if (reciprocal) {
780 for(int cl = 0; cl < numCells; cl++) {
781 for(int pt = 0; pt < numPoints; pt++) {
782 outputDataWrap(cl, pt) = inputDataRightWrap(pt)/inputDataLeftWrap(cl, 0);
783 } // P-loop
784 } // C-loop
785 }
786 else {
787 for(int cl = 0; cl < numCells; cl++) {
788 for(int pt = 0; pt < numPoints; pt++) {
789 outputDataWrap(cl, pt) = inputDataRightWrap(pt)*inputDataLeftWrap(cl, 0);
790 } // P-loop
791 } // C-loop
792 }
793 }// case 1
794 break;
795
796 case 2: {
797 if (reciprocal) {
798 for(int cl = 0; cl < numCells; cl++) {
799 for(int pt = 0; pt < numPoints; pt++) {
800 for( int iVec = 0; iVec < dim1Tens; iVec++) {
801 outputDataWrap(cl, pt, iVec) = inputDataRightWrap(pt, iVec)/inputDataLeftWrap(cl, 0);
802 } // D1-loop
803 } // P-loop
804 } // C-loop
805 }
806 else {
807 for(int cl = 0; cl < numCells; cl++) {
808 for(int pt = 0; pt < numPoints; pt++) {
809 for( int iVec = 0; iVec < dim1Tens; iVec++) {
810 outputDataWrap(cl, pt, iVec) = inputDataRightWrap(pt, iVec)*inputDataLeftWrap(cl, 0);
811 } // D1-loop
812 } // P-loop
813 } // C-loop
814 }
815 }// case 2
816 break;
817
818 case 3: {
819 if (reciprocal) {
820 for(int cl = 0; cl < numCells; cl++) {
821 for(int pt = 0; pt < numPoints; pt++) {
822 for( int iTens1 = 0; iTens1 < dim1Tens; iTens1++) {
823 for( int iTens2 = 0; iTens2 < dim2Tens; iTens2++) {
824 outputDataWrap(cl, pt, iTens1, iTens2) = inputDataRightWrap(pt, iTens1, iTens2)/inputDataLeftWrap(cl, 0);
825 } // D2-loop
826 } // D1-loop
827 } // P-loop
828 } // C-loop
829 }
830 else {
831 for(int cl = 0; cl < numCells; cl++) {
832 for(int pt = 0; pt < numPoints; pt++) {
833 for( int iTens1 = 0; iTens1 < dim1Tens; iTens1++) {
834 for( int iTens2 = 0; iTens2 < dim2Tens; iTens2++) {
835 outputDataWrap(cl, pt, iTens1, iTens2) = inputDataRightWrap(pt, iTens1, iTens2)*inputDataLeftWrap(cl, 0);
836 } // D2-loop
837 } // D1-loop
838 } // P-loop
839 } // C-loop
840 }
841 }// case 3
842 break;
843
844 default:
845 TEUCHOS_TEST_FOR_EXCEPTION( !( (invalRank == 1) || (invalRank == 2) || (invalRank == 3) ), std::invalid_argument,
846 ">>> ERROR (ArrayTools::scalarMultiplyDataData): This branch of the method is defined only for rank-1, 2 or 3 input containers.");
847
848 } // invalRank
849 } // numDataPoints
850
851 } // end if (outvalRank = invalRank)
852
853} // scalarMultiplyDataData
854} // end namespace Intrepid
855#endif
856
static void scalarMultiplyDataData(ArrayOutData &outputData, const ArrayInDataLeft &inputDataLeft, const ArrayInDataRight &inputDataRight, const bool reciprocal=false)
There are two use cases: (1) multiplies a rank-2, 3, or 4 container inputDataRight with dimensions (C...
static void scalarMultiplyDataField(ArrayOutFields &outputFields, const ArrayInData &inputData, const ArrayInFields &inputFields, const bool reciprocal=false)
There are two use cases: (1) multiplies a rank-3, 4, or 5 container inputFields with dimensions (C,...