Intrepid
Intrepid_FunctionSpaceToolsDef.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
50namespace Intrepid {
51
52template<class Scalar, class ArrayTypeOut, class ArrayTypeIn>
53void FunctionSpaceTools::HGRADtransformVALUE(ArrayTypeOut & outVals,
54 const ArrayTypeIn & inVals) {
55
56 ArrayTools::cloneFields<Scalar>(outVals, inVals);
57
58}
59
60template<class Scalar, class ArrayTypeOut, class ArrayTypeJac, class ArrayTypeIn>
61void FunctionSpaceTools::HGRADtransformGRAD(ArrayTypeOut & outVals,
62 const ArrayTypeJac & jacobianInverse,
63 const ArrayTypeIn & inVals,
64 const char transpose) {
65
66 ArrayTools::matvecProductDataField<Scalar>(outVals, jacobianInverse, inVals, transpose);
67
68}
69
70template<class Scalar, class ArrayTypeOut, class ArrayTypeJac, class ArrayTypeIn>
71void FunctionSpaceTools::HCURLtransformVALUE(ArrayTypeOut & outVals,
72 const ArrayTypeJac & jacobianInverse,
73 const ArrayTypeIn & inVals,
74 const char transpose) {
75
76 ArrayTools::matvecProductDataField<Scalar>(outVals, jacobianInverse, inVals, transpose);
77
78}
79
80template<class Scalar, class ArrayTypeOut, class ArrayTypeJac, class ArrayTypeDet, class ArrayTypeIn>
81void FunctionSpaceTools::HCURLtransformCURL(ArrayTypeOut & outVals,
82 const ArrayTypeJac & jacobian,
83 const ArrayTypeDet & jacobianDet,
84 const ArrayTypeIn & inVals,
85 const char transpose) {
86
87 ArrayTools::matvecProductDataField<Scalar>(outVals, jacobian, inVals, transpose);
88 ArrayTools::scalarMultiplyDataField<Scalar>(outVals, jacobianDet, outVals, true);
89
90}
91
92template<class Scalar, class ArrayTypeOut, class ArrayTypeJac, class ArrayTypeDet, class ArrayTypeIn>
93void FunctionSpaceTools::HDIVtransformVALUE(ArrayTypeOut & outVals,
94 const ArrayTypeJac & jacobian,
95 const ArrayTypeDet & jacobianDet,
96 const ArrayTypeIn & inVals,
97 const char transpose) {
98
99 ArrayTools::matvecProductDataField<Scalar>(outVals, jacobian, inVals, transpose);
100 ArrayTools::scalarMultiplyDataField<Scalar>(outVals, jacobianDet, outVals, true);
101
102}
103
104template<class Scalar, class ArrayTypeOut, class ArrayTypeDet, class ArrayTypeIn>
105void FunctionSpaceTools::HDIVtransformDIV(ArrayTypeOut & outVals,
106 const ArrayTypeDet & jacobianDet,
107 const ArrayTypeIn & inVals) {
108
109 ArrayTools::scalarMultiplyDataField<Scalar>(outVals, jacobianDet, inVals, true);
110
111}
112
113template<class Scalar, class ArrayTypeOut, class ArrayTypeDet, class ArrayTypeIn>
114void FunctionSpaceTools::HVOLtransformVALUE(ArrayTypeOut & outVals,
115 const ArrayTypeDet & jacobianDet,
116 const ArrayTypeIn & inVals) {
117
118 ArrayTools::scalarMultiplyDataField<Scalar>(outVals, jacobianDet, inVals, true);
119
120}
121template<class Scalar>
123 const Intrepid::FieldContainer<Scalar> & leftValues,
124 const Intrepid::FieldContainer<Scalar> & rightValues,
125 const ECompEngine compEngine,
126 const bool sumInto) {
127int outRank = getrank(outputValues);
128int lRank = getrank(leftValues);
129switch (outRank) {
130 case 1:{
131 switch (lRank){
132 case 2:{
133 #ifdef HAVE_INTREPID_DEBUG
134 TEUCHOS_TEST_FOR_EXCEPTION( (leftValues.rank() != 2 ), std::invalid_argument,
135 ">>> ERROR (ArrayTools::contractDataDataScalar): Rank of the left input argument must equal 2!");
136 TEUCHOS_TEST_FOR_EXCEPTION( (rightValues.rank() != 2 ), std::invalid_argument,
137 ">>> ERROR (ArrayTools::contractDataDataScalar): Rank of right input argument must equal 2!");
138 TEUCHOS_TEST_FOR_EXCEPTION( (outputValues.rank() != 1 ), std::invalid_argument,
139 ">>> ERROR (ArrayTools::contractDataDataScalar): Rank of output argument must equal 1!");
140 TEUCHOS_TEST_FOR_EXCEPTION( (leftValues.dimension(0) != rightValues.dimension(0) ), std::invalid_argument,
141 ">>> ERROR (ArrayTools::contractDataDataScalar): Zeroth dimensions (number of integration domains) of the left and right input containers must agree!");
142 TEUCHOS_TEST_FOR_EXCEPTION( (leftValues.dimension(1) != rightValues.dimension(1) ), std::invalid_argument,
143 ">>> ERROR (ArrayTools::contractDataDataScalar): First dimensions (numbers of integration points) of the left and right input containers must agree!");
144 TEUCHOS_TEST_FOR_EXCEPTION( (outputValues.dimension(0) != rightValues.dimension(0) ), std::invalid_argument,
145 ">>> ERROR (ArrayTools::contractDataDataScalar): Zeroth dimensions (numbers of integration domains) of the input and output containers must agree!");
146#endif
147
148 // get sizes
149 int numCells = leftValues.dimension(0);
150 int numPoints = leftValues.dimension(1);
151
152 switch(compEngine) {
153 case COMP_CPP: {
154 if (sumInto) {
155 for (int cl = 0; cl < numCells; cl++) {
156 Scalar tmpVal(0);
157 for (int qp = 0; qp < numPoints; qp++) {
158 tmpVal += leftValues(cl, qp)*rightValues(cl, qp);
159 } // P-loop
160 outputValues(cl) += tmpVal;
161 } // C-loop
162 }
163 else {
164 for (int cl = 0; cl < numCells; cl++) {
165 Scalar tmpVal(0);
166 for (int qp = 0; qp < numPoints; qp++) {
167 tmpVal += leftValues(cl, qp)*rightValues(cl, qp);
168 } // P-loop
169 outputValues(cl) = tmpVal;
170 } // C-loop
171 }
172 }
173 break;
174
175 case COMP_BLAS: {
176 int incr = 1; // increment
177 if (sumInto) {
178 for (int cl=0; cl < numCells; cl++) {
179 Teuchos::BLAS<int, Scalar> myblas;
180 outputValues(cl) += myblas.DOT(numPoints, &leftValues[cl*numPoints], incr, &rightValues[cl*numPoints], incr);
181 }
182 }
183 else {
184 for (int cl=0; cl < numCells; cl++) {
185 Teuchos::BLAS<int, Scalar> myblas;
186 outputValues(cl) = myblas.DOT(numPoints, &leftValues[cl*numPoints], incr, &rightValues[cl*numPoints], incr);
187 }
188 }
189 }
190 break;
191
192 default:
193 TEUCHOS_TEST_FOR_EXCEPTION( ( ~isValidCompEngine(compEngine) ), std::invalid_argument,
194 ">>> ERROR (ArrayTools::contractDataDataScalar): Computational engine not defined!");
195 } // switch(compEngine)
196 }
197 break;
198 case 3:{
199#ifdef HAVE_INTREPID_DEBUG
200 TEUCHOS_TEST_FOR_EXCEPTION( (leftValues.rank() != 3 ), std::invalid_argument,
201 ">>> ERROR (ArrayTools::contractDataDataVector): Rank of the left input argument must equal 3!");
202 TEUCHOS_TEST_FOR_EXCEPTION( (rightValues.rank() != 3 ), std::invalid_argument,
203 ">>> ERROR (ArrayTools::contractDataDataVector): Rank of right input argument must equal 3!");
204 TEUCHOS_TEST_FOR_EXCEPTION( (outputValues.rank() != 1 ), std::invalid_argument,
205 ">>> ERROR (ArrayTools::contractDataDataVector): Rank of output argument must equal 1!");
206 TEUCHOS_TEST_FOR_EXCEPTION( (leftValues.dimension(0) != rightValues.dimension(0) ), std::invalid_argument,
207 ">>> ERROR (ArrayTools::contractDataDataVector): Zeroth dimensions (number of integration domains) of the left and right input containers must agree!");
208 TEUCHOS_TEST_FOR_EXCEPTION( (leftValues.dimension(1) != rightValues.dimension(1) ), std::invalid_argument,
209 ">>> ERROR (ArrayTools::contractDataDataVector): First dimensions (numbers of integration points) of the left and right input containers must agree!");
210 TEUCHOS_TEST_FOR_EXCEPTION( (leftValues.dimension(2) != rightValues.dimension(2) ), std::invalid_argument,
211 ">>> ERROR (ArrayTools::contractDataDataVector): Second dimensions (numbers of vector components) of the left and right input containers must agree!");
212 TEUCHOS_TEST_FOR_EXCEPTION( (outputValues.dimension(0) != rightValues.dimension(0) ), std::invalid_argument,
213 ">>> ERROR (ArrayTools::contractDataDataVector): Zeroth dimensions (numbers of integration domains) of the input and output containers must agree!");
214#endif
215
216 // get sizes
217 int numCells = leftValues.dimension(0);
218 int numPoints = leftValues.dimension(1);
219 int dimVec = leftValues.dimension(2);
220
221 switch(compEngine) {
222 case COMP_CPP: {
223 if (sumInto) {
224 for (int cl = 0; cl < numCells; cl++) {
225 Scalar tmpVal(0);
226 for (int qp = 0; qp < numPoints; qp++) {
227 for (int iVec = 0; iVec < dimVec; iVec++) {
228 tmpVal += leftValues(cl, qp, iVec)*rightValues(cl, qp, iVec);
229 } // D-loop
230 } // P-loop
231 outputValues(cl) += tmpVal;
232 } // C-loop
233 }
234 else {
235 for (int cl = 0; cl < numCells; cl++) {
236 Scalar tmpVal(0);
237 for (int qp = 0; qp < numPoints; qp++) {
238 for (int iVec = 0; iVec < dimVec; iVec++) {
239 tmpVal += leftValues(cl, qp, iVec)*rightValues(cl, qp, iVec);
240 } // D-loop
241 } // P-loop
242 outputValues(cl) = tmpVal;
243 } // C-loop
244 }
245 }
246 break;
247
248 case COMP_BLAS: {
249 int skip = numPoints*dimVec; // size of the left data chunk per cell
250 int incr = 1; // increment
251 if (sumInto) {
252 for (int cl=0; cl < numCells; cl++) {
253 Teuchos::BLAS<int, Scalar> myblas;
254 outputValues(cl) += myblas.DOT(skip, &leftValues[cl*skip], incr, &rightValues[cl*skip], incr);
255 }
256 }
257 else {
258 for (int cl=0; cl < numCells; cl++) {
259 Teuchos::BLAS<int, Scalar> myblas;
260 outputValues(cl) = myblas.DOT(skip, &leftValues[cl*skip], incr, &rightValues[cl*skip], incr);
261 }
262 }
263 }
264 break;
265
266 default:
267 TEUCHOS_TEST_FOR_EXCEPTION( ( ~isValidCompEngine(compEngine) ), std::invalid_argument,
268 ">>> ERROR (ArrayTools::contractDataDataVector): Computational engine not defined!");
269 } // switch(compEngine)
270
271
272}
273
274 break;
275 case 4:{
276#ifdef HAVE_INTREPID_DEBUG
277 TEUCHOS_TEST_FOR_EXCEPTION( (leftValues.rank() != 4 ), std::invalid_argument,
278 ">>> ERROR (ArrayTools::contractDataDataTensor): Rank of the left input argument must equal 4");
279 TEUCHOS_TEST_FOR_EXCEPTION( (rightValues.rank() != 4 ), std::invalid_argument,
280 ">>> ERROR (ArrayTools::contractDataDataTensor): Rank of right input argument must equal 4!");
281 TEUCHOS_TEST_FOR_EXCEPTION( (outputValues.rank() != 1 ), std::invalid_argument,
282 ">>> ERROR (ArrayTools::contractDataDataTensor): Rank of output argument must equal 1!");
283 TEUCHOS_TEST_FOR_EXCEPTION( (leftValues.dimension(0) != rightValues.dimension(0) ), std::invalid_argument,
284 ">>> ERROR (ArrayTools::contractDataDataTensor): Zeroth dimensions (number of integration domains) of the left and right input containers must agree!");
285 TEUCHOS_TEST_FOR_EXCEPTION( (leftValues.dimension(1) != rightValues.dimension(1) ), std::invalid_argument,
286 ">>> ERROR (ArrayTools::contractDataDataTensor): First dimensions (numbers of integration points) of the left and right input containers must agree!");
287 TEUCHOS_TEST_FOR_EXCEPTION( (leftValues.dimension(2) != rightValues.dimension(2) ), std::invalid_argument,
288 ">>> ERROR (ArrayTools::contractDataDataTensor): Second dimensions (first tensor dimensions) of the left and right input containers must agree!");
289 TEUCHOS_TEST_FOR_EXCEPTION( (leftValues.dimension(3) != rightValues.dimension(3) ), std::invalid_argument,
290 ">>> ERROR (ArrayTools::contractDataDataTensor): Third dimensions (second tensor dimensions) of the left and right input containers must agree!");
291 TEUCHOS_TEST_FOR_EXCEPTION( (outputValues.dimension(0) != rightValues.dimension(0) ), std::invalid_argument,
292 ">>> ERROR (ArrayTools::contractDataDataTensor): Zeroth dimensions (numbers of integration domains) of the input and output containers must agree!");
293#endif
294
295 // get sizes
296 int numCells = leftValues.dimension(0);
297 int numPoints = leftValues.dimension(1);
298 int dim1Tensor = leftValues.dimension(2);
299 int dim2Tensor = leftValues.dimension(3);
300
301 switch(compEngine) {
302 case COMP_CPP: {
303 if (sumInto) {
304 for (int cl = 0; cl < numCells; cl++) {
305 Scalar tmpVal(0);
306 for (int qp = 0; qp < numPoints; qp++) {
307 for (int iTens1 = 0; iTens1 < dim1Tensor; iTens1++) {
308 for (int iTens2 = 0; iTens2 < dim2Tensor; iTens2++) {
309 tmpVal += leftValues(cl, qp, iTens1, iTens2)*rightValues(cl, qp, iTens1, iTens2);
310 } // D2-loop
311 } // D1-loop
312 } // P-loop
313 outputValues(cl) += tmpVal;
314 } // C-loop
315 }
316 else {
317 for (int cl = 0; cl < numCells; cl++) {
318 Scalar tmpVal(0);
319 for (int qp = 0; qp < numPoints; qp++) {
320 for (int iTens1 = 0; iTens1 < dim1Tensor; iTens1++) {
321 for (int iTens2 = 0; iTens2 < dim2Tensor; iTens2++) {
322 tmpVal += leftValues(cl, qp, iTens1, iTens2)*rightValues(cl, qp, iTens1, iTens2);
323 } // D2-loop
324 } // D1-loop
325 } // P-loop
326 outputValues(cl) = tmpVal;
327 } // C-loop
328 }
329 }
330 break;
331
332 case COMP_BLAS: {
333 int skip = numPoints*dim1Tensor*dim2Tensor; // size of the left data chunk per cell
334 int incr = 1; // increment
335 if (sumInto) {
336 for (int cl=0; cl < numCells; cl++) {
337 Teuchos::BLAS<int, Scalar> myblas;
338 outputValues(cl) += myblas.DOT(skip, &leftValues[cl*skip], incr, &rightValues[cl*skip], incr);
339 }
340 }
341 else {
342 for (int cl=0; cl < numCells; cl++) {
343 Teuchos::BLAS<int, Scalar> myblas;
344 outputValues(cl) = myblas.DOT(skip, &leftValues[cl*skip], incr, &rightValues[cl*skip], incr);
345 }
346 }
347 }
348 break;
349
350 default:
351 TEUCHOS_TEST_FOR_EXCEPTION( ( ~isValidCompEngine(compEngine) ), std::invalid_argument,
352 ">>> ERROR (ArrayTools::contractDataDataTensor): Computational engine not defined!");
353 } // switch(compEngine)
354
355}
356 break;
357 default:
358#ifdef HAVE_INTREPID_DEBUG
359
360 TEUCHOS_TEST_FOR_EXCEPTION( ((lRank != 2) && (lRank != 3) && (lRank != 4)), std::invalid_argument,
361 ">>> ERROR (FunctionSpaceTools::dataIntegral): Left data input container must have rank 2, 3 or 4.");
362
363#endif
364 break;
365
366}
367
368}
369 break;
370 case 2:{
371 switch (lRank) {
372 case 2:{
373#ifdef HAVE_INTREPID_DEBUG
374 TEUCHOS_TEST_FOR_EXCEPTION( (rightValues.rank() != 3 ), std::invalid_argument,
375 ">>> ERROR (ArrayTools::contractDataFieldScalar): Rank of the fields input argument must equal 3!");
376 TEUCHOS_TEST_FOR_EXCEPTION( (leftValues.rank() != 2 ), std::invalid_argument,
377 ">>> ERROR (ArrayTools::contractDataFieldScalar): Rank of the data input argument must equal 2!");
378 TEUCHOS_TEST_FOR_EXCEPTION( (outputValues.rank() != 2 ), std::invalid_argument,
379 ">>> ERROR (ArrayTools::contractDataFieldScalar): Rank of output argument must equal 2!");
380 TEUCHOS_TEST_FOR_EXCEPTION( (rightValues.dimension(0) != leftValues.dimension(0) ), std::invalid_argument,
381 ">>> ERROR (ArrayTools::contractDataFieldScalar): Zeroth dimensions (number of integration domains) of the fields and data input containers must agree!");
382 TEUCHOS_TEST_FOR_EXCEPTION( ( (rightValues.dimension(2) != leftValues.dimension(1)) && (leftValues.dimension(1) != 1) ), std::invalid_argument,
383 ">>> 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!");
384 TEUCHOS_TEST_FOR_EXCEPTION( (outputValues.dimension(0) != rightValues.dimension(0) ), std::invalid_argument,
385 ">>> ERROR (ArrayTools::contractDataFieldScalar): Zeroth dimensions (numbers of integration domains) of the fields input and output containers must agree!");
386 TEUCHOS_TEST_FOR_EXCEPTION( (outputValues.dimension(1) != rightValues.dimension(1) ), std::invalid_argument,
387 ">>> ERROR (ArrayTools::contractDataFieldScalar): First dimensions (number of fields) of the fields input and output containers must agree!");
388#endif
389
390 // get sizes
391 int numCells = rightValues.dimension(0);
392 int numFields = rightValues.dimension(1);
393 int numPoints = rightValues.dimension(2);
394 int numDataPoints = leftValues.dimension(1);
395
396 ECompEngine myCompEngine = (numDataPoints == 1 ? COMP_CPP : compEngine);
397
398 switch(myCompEngine) {
399 case COMP_CPP: {
400 if (sumInto) {
401 if (numDataPoints != 1) { // nonconstant data
402 for (int cl = 0; cl < numCells; cl++) {
403 for (int lbf = 0; lbf < numFields; lbf++) {
404 Scalar tmpVal(0);
405 for (int qp = 0; qp < numPoints; qp++) {
406 tmpVal += rightValues(cl, lbf, qp)*leftValues(cl, qp);
407 } // P-loop
408 outputValues(cl, lbf) += tmpVal;
409 } // F-loop
410 } // C-loop
411 }
412 else { // constant data
413 for (int cl = 0; cl < numCells; cl++) {
414 for (int lbf = 0; lbf < numFields; lbf++) {
415 Scalar tmpVal(0);
416 for (int qp = 0; qp < numPoints; qp++) {
417 tmpVal += rightValues(cl, lbf, qp)*leftValues(cl, 0);
418 } // P-loop
419 outputValues(cl, lbf) += tmpVal;
420 } // F-loop
421 } // C-loop
422 } // numDataPoints
423 }
424 else {
425 if (numDataPoints != 1) { // nonconstant data
426 for (int cl = 0; cl < numCells; cl++) {
427 for (int lbf = 0; lbf < numFields; lbf++) {
428 Scalar tmpVal(0);
429 for (int qp = 0; qp < numPoints; qp++) {
430 tmpVal += rightValues(cl, lbf, qp)*leftValues(cl, qp);
431 } // P-loop
432 outputValues(cl, lbf) = tmpVal;
433 } // F-loop
434 } // C-loop
435 }
436 else { // constant data
437 for (int cl = 0; cl < numCells; cl++) {
438 for (int lbf = 0; lbf < numFields; lbf++) {
439 Scalar tmpVal(0);
440 for (int qp = 0; qp < numPoints; qp++) {
441 tmpVal += rightValues(cl, lbf, qp)*leftValues(cl, 0);
442 } // P-loop
443 outputValues(cl, lbf) = tmpVal;
444 } // F-loop
445 } // C-loop
446 } // numDataPoints
447 }
448 }
449 break;
450
451 case COMP_BLAS: {
452 /*
453 GEMM parameters and their values.
454 (Note: It is assumed that the result needs to be transposed into row-major format.
455 Think of left and right input matrices as A(p x f) and B(p x 1), respectively,
456 even though the indexing is ((C),F,P) and ((C),P). Due to BLAS formatting
457 assumptions, we are computing (A^T*B)^T = B^T*A.)
458 TRANSA TRANS
459 TRANSB NO_TRANS
460 M #rows(B^T) = 1
461 N #cols(A) = number of input fields
462 K #cols(B^T) = number of integration points * size of data
463 ALPHA 1.0
464 A right data for cell cl = &rightFields[cl*skipR]
465 LDA #rows(B) = number of integration points * size of data
466 B left data for cell cl = &leftFields[cl*skipL]
467 LDB #rows(A) = number of integration points * size of data
468 BETA 0.0
469 C result for cell cl = &outputFields[cl*skipOp]
470 LDC #rows(C) = 1
471 */
472 int numData = numPoints;
473 int skipL = numFields*numPoints; // size of the left data chunk per cell
474 int skipR = numPoints; // size of the right data chunk per cell
475 int skipOp = numFields; // size of the output data chunk per cell
476 Scalar alpha(1.0); // these are left unchanged by GEMM
477 Scalar beta(0.0);
478 if (sumInto) {
479 beta = 1.0;
480 }
481
482 for (int cl=0; cl < numCells; cl++) {
483 /* Use this if data is used in row-major format */
484 Teuchos::BLAS<int, Scalar> myblas;
485 myblas.GEMM(Teuchos::TRANS, Teuchos::NO_TRANS,
486 1, numFields, numData,
487 alpha, &leftValues[cl*skipR], numData,
488 &rightValues[cl*skipL], numData,
489 beta, &outputValues[cl*skipOp], 1);
490 /* Use this if data is used in column-major format */
491 /*
492 myblas.GEMM(Teuchos::TRANS, Teuchos::NO_TRANS,
493 numFields, 1, numData,
494 alpha, &inputFields[cl*skipL], numData,
495 &inputData[cl*skipR], numData,
496 beta, &outputFields[cl*skipOp], numFields);
497 */
498 }
499 }
500 break;
501
502 default:
503 TEUCHOS_TEST_FOR_EXCEPTION( ( ~isValidCompEngine(compEngine) ), std::invalid_argument,
504 ">>> ERROR (ArrayTools::contractDataFieldScalar): Computational engine not defined!");
505 } // switch(compEngine)
506 }
507 break;
508 case 3:{
509 #ifdef HAVE_INTREPID_DEBUG
510 TEUCHOS_TEST_FOR_EXCEPTION( (rightValues.rank() != 4 ), std::invalid_argument,
511 ">>> ERROR (ArrayTools::contractDataFieldVector): Rank of the fields input argument must equal 4!");
512 TEUCHOS_TEST_FOR_EXCEPTION( (leftValues.rank() != 3 ), std::invalid_argument,
513 ">>> ERROR (ArrayTools::contractDataFieldVector): Rank of the data input argument must equal 3!");
514 TEUCHOS_TEST_FOR_EXCEPTION( (outputValues.rank() != 2 ), std::invalid_argument,
515 ">>> ERROR (ArrayTools::contractDataFieldVector): Rank of output argument must equal 2!");
516 TEUCHOS_TEST_FOR_EXCEPTION( (rightValues.dimension(0) != leftValues.dimension(0) ), std::invalid_argument,
517 ">>> ERROR (ArrayTools::contractDataFieldVector): Zeroth dimensions (number of integration domains) of the fields and data input containers must agree!");
518 TEUCHOS_TEST_FOR_EXCEPTION( ( (rightValues.dimension(2) != leftValues.dimension(1)) && (leftValues.dimension(1) != 1) ), std::invalid_argument,
519 ">>> 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!");
520 TEUCHOS_TEST_FOR_EXCEPTION( (rightValues.dimension(3) != leftValues.dimension(2) ), std::invalid_argument,
521 ">>> ERROR (ArrayTools::contractDataFieldVector): Third dimension of the fields input container and second dimension of data input container (vector index) must agree!");
522 TEUCHOS_TEST_FOR_EXCEPTION( (outputValues.dimension(0) != rightValues.dimension(0) ), std::invalid_argument,
523 ">>> ERROR (ArrayTools::contractDataFieldVector): Zeroth dimensions (numbers of integration domains) of the fields input and output containers must agree!");
524 TEUCHOS_TEST_FOR_EXCEPTION( (outputValues.dimension(1) != rightValues.dimension(1) ), std::invalid_argument,
525 ">>> ERROR (ArrayTools::contractDataFieldVector): First dimensions of output container and fields input container (number of fields) must agree!");
526#endif
527
528 // get sizes
529 int numCells = rightValues.dimension(0);
530 int numFields = rightValues.dimension(1);
531 int numPoints = rightValues.dimension(2);
532 int dimVec = rightValues.dimension(3);
533 int numDataPoints = leftValues.dimension(1);
534
535 ECompEngine myCompEngine = (numDataPoints == 1 ? COMP_CPP : compEngine);
536
537 switch(myCompEngine) {
538 case COMP_CPP: {
539 if (sumInto) {
540 if (numDataPoints != 1) { // nonconstant data
541 for (int cl = 0; cl < numCells; cl++) {
542 for (int lbf = 0; lbf < numFields; lbf++) {
543 Scalar tmpVal(0);
544 for (int qp = 0; qp < numPoints; qp++) {
545 for (int iVec = 0; iVec < dimVec; iVec++) {
546 tmpVal += rightValues(cl, lbf, qp, iVec)*leftValues(cl, qp, iVec);
547 } // D-loop
548 } // P-loop
549 outputValues(cl, lbf) += tmpVal;
550 } // F-loop
551 } // C-loop
552 }
553 else { // constant data
554 for (int cl = 0; cl < numCells; cl++) {
555 for (int lbf = 0; lbf < numFields; lbf++) {
556 Scalar tmpVal(0);
557 for (int qp = 0; qp < numPoints; qp++) {
558 for (int iVec = 0; iVec < dimVec; iVec++) {
559 tmpVal += rightValues(cl, lbf, qp, iVec)*leftValues(cl, 0, iVec);
560 } //D-loop
561 } // P-loop
562 outputValues(cl, lbf) += tmpVal;
563 } // F-loop
564 } // C-loop
565 } // numDataPoints
566 }
567 else {
568 if (numDataPoints != 1) { // nonconstant data
569 for (int cl = 0; cl < numCells; cl++) {
570 for (int lbf = 0; lbf < numFields; lbf++) {
571 Scalar tmpVal(0);
572 for (int qp = 0; qp < numPoints; qp++) {
573 for (int iVec = 0; iVec < dimVec; iVec++) {
574 tmpVal += rightValues(cl, lbf, qp, iVec)*leftValues(cl, qp, iVec);
575 } // D-loop
576 } // P-loop
577 outputValues(cl, lbf) = tmpVal;
578 } // F-loop
579 } // C-loop
580 }
581 else { // constant data
582 for (int cl = 0; cl < numCells; cl++) {
583 for (int lbf = 0; lbf < numFields; lbf++) {
584 Scalar tmpVal(0);
585 for (int qp = 0; qp < numPoints; qp++) {
586 for (int iVec = 0; iVec < dimVec; iVec++) {
587 tmpVal += rightValues(cl, lbf, qp, iVec)*leftValues(cl, 0, iVec);
588 } //D-loop
589 } // P-loop
590 outputValues(cl, lbf) = tmpVal;
591 } // F-loop
592 } // C-loop
593 } // numDataPoints
594 }
595 }
596 break;
597
598 case COMP_BLAS: {
599 /*
600 GEMM parameters and their values.
601 (Note: It is assumed that the result needs to be transposed into row-major format.
602 Think of left and right input matrices as A(p x f) and B(p x 1), respectively,
603 even though the indexing is ((C),F,P) and ((C),P). Due to BLAS formatting
604 assumptions, we are computing (A^T*B)^T = B^T*A.)
605 TRANSA TRANS
606 TRANSB NO_TRANS
607 M #rows(B^T) = 1
608 N #cols(A) = number of input fields
609 K #cols(B^T) = number of integration points * size of data
610 ALPHA 1.0
611 A right data for cell cl = &rightFields[cl*skipR]
612 LDA #rows(B) = number of integration points * size of data
613 B left data for cell cl = &leftFields[cl*skipL]
614 LDB #rows(A) = number of integration points * size of data
615 BETA 0.0
616 C result for cell cl = &outputFields[cl*skipOp]
617 LDC #rows(C) = 1
618 */
619 int numData = numPoints*dimVec;
620 int skipL = numFields*numData; // size of the left data chunk per cell
621 int skipR = numData; // size of the right data chunk per cell
622 int skipOp = numFields; // size of the output data chunk per cell
623 Scalar alpha(1.0); // these are left unchanged by GEMM
624 Scalar beta(0.0);
625 if (sumInto) {
626 beta = 1.0;
627 }
628
629 for (int cl=0; cl < numCells; cl++) {
630 /* Use this if data is used in row-major format */
631 Teuchos::BLAS<int, Scalar> myblas;
632 myblas.GEMM(Teuchos::TRANS, Teuchos::NO_TRANS,
633 1, numFields, numData,
634 alpha, &leftValues[cl*skipR], numData,
635 &rightValues[cl*skipL], numData,
636 beta, &outputValues[cl*skipOp], 1);
637 /* Use this if data is used in column-major format */
638 /*
639 myblas.GEMM(Teuchos::TRANS, Teuchos::NO_TRANS,
640 numFields, 1, numData,
641 alpha, &inputFields[cl*skipL], numData,
642 &inputData[cl*skipR], numData,
643 beta, &outputFields[cl*skipOp], numFields);
644 */
645 }
646 }
647 break;
648
649 default:
650 TEUCHOS_TEST_FOR_EXCEPTION( ( ~isValidCompEngine(compEngine) ), std::invalid_argument,
651 ">>> ERROR (ArrayTools::contractDataFieldVector): Computational engine not defined!");
652 } // switch(compEngine)
653 }
654 break;
655 case 4:{
656#ifdef HAVE_INTREPID_DEBUG
657 TEUCHOS_TEST_FOR_EXCEPTION( (rightValues.rank() != 5 ), std::invalid_argument,
658 ">>> ERROR (ArrayTools::contractDataFieldTensor): Rank of the fields input argument must equal 5!");
659 TEUCHOS_TEST_FOR_EXCEPTION( (leftValues.rank() != 4 ), std::invalid_argument,
660 ">>> ERROR (ArrayTools::contractDataFieldTensor): Rank of the data input argument must equal 4!");
661 TEUCHOS_TEST_FOR_EXCEPTION( (outputValues.rank() != 2 ), std::invalid_argument,
662 ">>> ERROR (ArrayTools::contractDataFieldTensor): Rank of output argument must equal 2!");
663 TEUCHOS_TEST_FOR_EXCEPTION( (rightValues.dimension(0) != leftValues.dimension(0) ), std::invalid_argument,
664 ">>> ERROR (ArrayTools::contractDataFieldTensor): Zeroth dimensions (number of integration domains) of the fields and data input containers must agree!");
665 TEUCHOS_TEST_FOR_EXCEPTION( ( (rightValues.dimension(2) != leftValues.dimension(1)) && (leftValues.dimension(1) != 1) ), std::invalid_argument,
666 ">>> 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!");
667 TEUCHOS_TEST_FOR_EXCEPTION( (rightValues.dimension(3) != leftValues.dimension(2) ), std::invalid_argument,
668 ">>> ERROR (ArrayTools::contractDataFieldTensor): Third dimension of the fields input container and second dimension of data input container (first tensor dimension) must agree!");
669 TEUCHOS_TEST_FOR_EXCEPTION( (rightValues.dimension(4) != leftValues.dimension(3) ), std::invalid_argument,
670 ">>> ERROR (ArrayTools::contractDataFieldTensor): Fourth dimension of the fields input container and third dimension of data input container (second tensor dimension) must agree!");
671 TEUCHOS_TEST_FOR_EXCEPTION( (outputValues.dimension(0) != rightValues.dimension(0) ), std::invalid_argument,
672 ">>> ERROR (ArrayTools::contractDataFieldTensor): Zeroth dimensions (numbers of integration domains) of the fields input and output containers must agree!");
673 TEUCHOS_TEST_FOR_EXCEPTION( (outputValues.dimension(1) != rightValues.dimension(1) ), std::invalid_argument,
674 ">>> ERROR (ArrayTools::contractDataFieldTensor): First dimensions (number of fields) of output container and fields input container must agree!");
675#endif
676
677 // get sizes
678 int numCells = rightValues.dimension(0);
679 int numFields = rightValues.dimension(1);
680 int numPoints = rightValues.dimension(2);
681 int dim1Tens = rightValues.dimension(3);
682 int dim2Tens = rightValues.dimension(4);
683 int numDataPoints = leftValues.dimension(1);
684
685 ECompEngine myCompEngine = (numDataPoints == 1 ? COMP_CPP : compEngine);
686
687 switch(myCompEngine) {
688 case COMP_CPP: {
689 if (sumInto) {
690 if (numDataPoints != 1) { // nonconstant data
691 for (int cl = 0; cl < numCells; cl++) {
692 for (int lbf = 0; lbf < numFields; lbf++) {
693 Scalar tmpVal(0);
694 for (int qp = 0; qp < numPoints; qp++) {
695 for (int iTens1 = 0; iTens1 < dim1Tens; iTens1++) {
696 for (int iTens2 =0; iTens2 < dim2Tens; iTens2++) {
697 tmpVal += rightValues(cl, lbf, qp, iTens1, iTens2)*leftValues(cl, qp, iTens1, iTens2);
698 } // D2-loop
699 } // D1-loop
700 } // P-loop
701 outputValues(cl, lbf) += tmpVal;
702 } // F-loop
703 } // C-loop
704 }
705 else { // constant data
706 for (int cl = 0; cl < numCells; cl++) {
707 for (int lbf = 0; lbf < numFields; lbf++) {
708 Scalar tmpVal(0);
709 for (int qp = 0; qp < numPoints; qp++) {
710 for (int iTens1 = 0; iTens1 < dim1Tens; iTens1++) {
711 for (int iTens2 = 0; iTens2 < dim2Tens; iTens2++) {
712 tmpVal += rightValues(cl, lbf, qp, iTens1, iTens2)*leftValues(cl, 0, iTens1, iTens2);
713 } // D2-loop
714 } // D1-loop
715 } // P-loop
716 outputValues(cl, lbf) += tmpVal;
717 } // F-loop
718 } // C-loop
719 } // numDataPoints
720 }
721 else {
722 if (numDataPoints != 1) { // nonconstant data
723 for (int cl = 0; cl < numCells; cl++) {
724 for (int lbf = 0; lbf < numFields; lbf++) {
725 Scalar tmpVal(0);
726 for (int qp = 0; qp < numPoints; qp++) {
727 for (int iTens1 = 0; iTens1 < dim1Tens; iTens1++) {
728 for (int iTens2 =0; iTens2 < dim2Tens; iTens2++) {
729 tmpVal += rightValues(cl, lbf, qp, iTens1, iTens2)*leftValues(cl, qp, iTens1, iTens2);
730 } // D2-loop
731 } // D1-loop
732 } // P-loop
733 outputValues(cl, lbf) = tmpVal;
734 } // F-loop
735 } // C-loop
736 }
737 else { // constant data
738 for (int cl = 0; cl < numCells; cl++) {
739 for (int lbf = 0; lbf < numFields; lbf++) {
740 Scalar tmpVal(0);
741 for (int qp = 0; qp < numPoints; qp++) {
742 for (int iTens1 = 0; iTens1 < dim1Tens; iTens1++) {
743 for (int iTens2 = 0; iTens2 < dim2Tens; iTens2++) {
744 tmpVal += rightValues(cl, lbf, qp, iTens1, iTens2)*leftValues(cl, 0, iTens1, iTens2);
745 } // D2-loop
746 } // D1-loop
747 } // P-loop
748 outputValues(cl, lbf) = tmpVal;
749 } // F-loop
750 } // C-loop
751 } // numDataPoints
752 }
753 }
754 break;
755
756 case COMP_BLAS: {
757 /*
758 GEMM parameters and their values.
759 (Note: It is assumed that the result needs to be transposed into row-major format.
760 Think of left and right input matrices as A(p x f) and B(p x 1), respectively,
761 even though the indexing is ((C),F,P) and ((C),P). Due to BLAS formatting
762 assumptions, we are computing (A^T*B)^T = B^T*A.)
763 TRANSA TRANS
764 TRANSB NO_TRANS
765 M #rows(B^T) = 1
766 N #cols(A) = number of input fields
767 K #cols(B^T) = number of integration points * size of data
768 ALPHA 1.0
769 A right data for cell cl = &rightFields[cl*skipR]
770 LDA #rows(B) = number of integration points * size of data
771 B left data for cell cl = &leftFields[cl*skipL]
772 LDB #rows(A) = number of integration points * size of data
773 BETA 0.0
774 C result for cell cl = &outputFields[cl*skipOp]
775 LDC #rows(C) = 1
776 */
777 int numData = numPoints*dim1Tens*dim2Tens;
778 int skipL = numFields*numData; // size of the left data chunk per cell
779 int skipR = numData; // size of the right data chunk per cell
780 int skipOp = numFields; // size of the output data chunk per cell
781 Scalar alpha(1.0); // these are left unchanged by GEMM
782 Scalar beta(0.0);
783 if (sumInto) {
784 beta = 1.0;
785 }
786
787 for (int cl=0; cl < numCells; cl++) {
788 /* Use this if data is used in row-major format */
789 Teuchos::BLAS<int, Scalar> myblas;
790 myblas.GEMM(Teuchos::TRANS, Teuchos::NO_TRANS,
791 1, numFields, numData,
792 alpha, &leftValues[cl*skipR], numData,
793 &rightValues[cl*skipL], numData,
794 beta, &outputValues[cl*skipOp], 1);
795 /* Use this if data is used in column-major format */
796 /*
797 myblas.GEMM(Teuchos::TRANS, Teuchos::NO_TRANS,
798 numFields, 1, numData,
799 alpha, &inputFields[cl*skipL], numData,
800 &inputData[cl*skipR], numData,
801 beta, &outputFields[cl*skipOp], numFields);
802 */
803 }
804 }
805 break;
806
807 default:
808 TEUCHOS_TEST_FOR_EXCEPTION( ( ~isValidCompEngine(compEngine) ), std::invalid_argument,
809 ">>> ERROR (ArrayTools::contractDataFieldTensor): Computational engine not defined!");
810 } // switch(compEngine)
811}
812 break;
813 default:
814#ifdef HAVE_INTREPID_DEBUG
815
816TEUCHOS_TEST_FOR_EXCEPTION( ((lRank != 2) && (lRank != 3) && (lRank != 4)), std::invalid_argument,
817 ">>> ERROR (FunctionSpaceTools::functionalIntegral): Data input container must have rank 2, 3 or 4.");
818
819#endif
820 break;
821}
822
823}
824 break;
825 case 3:{
826 switch (lRank) {
827 case 3:{
828#ifdef HAVE_INTREPID_DEBUG
829 TEUCHOS_TEST_FOR_EXCEPTION( (leftValues.rank() != 3 ), std::invalid_argument,
830 ">>> ERROR (ArrayTools::contractFieldFieldScalar): Rank of the left input argument must equal 3!");
831 TEUCHOS_TEST_FOR_EXCEPTION( (rightValues.rank() != 3 ), std::invalid_argument,
832 ">>> ERROR (ArrayTools::contractFieldFieldScalar): Rank of right input argument must equal 3!");
833 TEUCHOS_TEST_FOR_EXCEPTION( (outputValues.rank() != 3 ), std::invalid_argument,
834 ">>> ERROR (ArrayTools::contractFieldFieldScalar): Rank of output argument must equal 3!");
835 TEUCHOS_TEST_FOR_EXCEPTION( (leftValues.dimension(0) != rightValues.dimension(0) ), std::invalid_argument,
836 ">>> ERROR (ArrayTools::contractFieldFieldScalar): Zeroth dimensions (number of integration domains) of the left and right input containers must agree!");
837 TEUCHOS_TEST_FOR_EXCEPTION( (leftValues.dimension(2) != rightValues.dimension(2) ), std::invalid_argument,
838 ">>> ERROR (ArrayTools::contractFieldFieldScalar): Second dimensions (numbers of integration points) of the left and right input containers must agree!");
839 TEUCHOS_TEST_FOR_EXCEPTION( (outputValues.dimension(0) != rightValues.dimension(0) ), std::invalid_argument,
840 ">>> ERROR (ArrayTools::contractFieldFieldScalar): Zeroth dimensions (numbers of integration domains) of the input and output containers must agree!");
841 TEUCHOS_TEST_FOR_EXCEPTION( (outputValues.dimension(1) != leftValues.dimension(1) ), std::invalid_argument,
842 ">>> ERROR (ArrayTools::contractFieldFieldScalar): First dimension of output container and first dimension of left input container must agree!");
843 TEUCHOS_TEST_FOR_EXCEPTION( (outputValues.dimension(2) != rightValues.dimension(1) ), std::invalid_argument,
844 ">>> ERROR (ArrayTools::contractFieldFieldScalar): Second dimension of output container and first dimension of right input container must agree!");
845#endif
846
847 // get sizes
848 int numCells = leftValues.dimension(0);
849 int numLeftFields = leftValues.dimension(1);
850 int numRightFields = rightValues.dimension(1);
851 int numPoints = leftValues.dimension(2);
852
853 switch(compEngine) {
854 case COMP_CPP: {
855 if (sumInto) {
856 for (int cl = 0; cl < numCells; cl++) {
857 for (int lbf = 0; lbf < numLeftFields; lbf++) {
858 for (int rbf = 0; rbf < numRightFields; rbf++) {
859 Scalar tmpVal(0);
860 for (int qp = 0; qp < numPoints; qp++) {
861 tmpVal += leftValues(cl, lbf, qp)*rightValues(cl, rbf, qp);
862 } // P-loop
863 outputValues(cl, lbf, rbf) += tmpVal;
864 } // R-loop
865 } // L-loop
866 } // C-loop
867 }
868 else {
869 for (int cl = 0; cl < numCells; cl++) {
870 for (int lbf = 0; lbf < numLeftFields; lbf++) {
871 for (int rbf = 0; rbf < numRightFields; rbf++) {
872 Scalar tmpVal(0);
873 for (int qp = 0; qp < numPoints; qp++) {
874 tmpVal += leftValues(cl, lbf, qp)*rightValues(cl, rbf, qp);
875 } // P-loop
876 outputValues(cl, lbf, rbf) = tmpVal;
877 } // R-loop
878 } // L-loop
879 } // C-loop
880 }
881 }
882 break;
883
884 case COMP_BLAS: {
885 /*
886 GEMM parameters and their values.
887 (Note: It is assumed that the result needs to be transposed into row-major format.
888 Think of left and right input matrices as A(p x l) and B(p x r), respectively,
889 even though the indexing is ((C),L,P) and ((C),R,P). Due to BLAS formatting
890 assumptions, we are computing (A^T*B)^T = B^T*A.)
891 TRANSA TRANS
892 TRANSB NO_TRANS
893 M #rows(B^T) = number of right fields
894 N #cols(A) = number of left fields
895 K #cols(B^T) = number of integration points
896 ALPHA 1.0
897 A right data for cell cl = &rightFields[cl*skipR]
898 LDA #rows(B) = number of integration points
899 B left data for cell cl = &leftFields[cl*skipL]
900 LDB #rows(A) = number of integration points
901 BETA 0.0
902 C result for cell cl = &outputFields[cl*skipOp]
903 LDC #rows(C) = number of right fields
904 */
905 int skipL = numLeftFields*numPoints; // size of the left data chunk per cell
906 int skipR = numRightFields*numPoints; // size of the right data chunk per cell
907 int skipOp = numLeftFields*numRightFields; // size of the output data chunk per cell
908 Scalar alpha(1.0); // these are left unchanged by GEMM
909 Scalar beta(0.0);
910 if (sumInto) {
911 beta = 1.0;
912 }
913
914 for (int cl=0; cl < numCells; cl++) {
915 /* Use this if data is used in row-major format */
916 Teuchos::BLAS<int, Scalar> myblas;
917 myblas.GEMM(Teuchos::TRANS, Teuchos::NO_TRANS,
918 numRightFields, numLeftFields, numPoints,
919 alpha, &rightValues[cl*skipR], numPoints,
920 &leftValues[cl*skipL], numPoints,
921 beta, &outputValues[cl*skipOp], numRightFields);
922 /* Use this if data is used in column-major format */
923 /*
924 myblas.GEMM(Teuchos::TRANS, Teuchos::NO_TRANS,
925 numLeftFields, numRightFields, numPoints,
926 alpha, &leftFields[cl*skipL], numPoints,
927 &rightFields[cl*skipR], numPoints,
928 beta, &outputFields[cl*skipOp], numLeftFields);
929 */
930 }
931 }
932 break;
933
934 default:
935 TEUCHOS_TEST_FOR_EXCEPTION( ( ~isValidCompEngine(compEngine) ), std::invalid_argument,
936 ">>> ERROR (ArrayTools::contractFieldFieldScalar): Computational engine not defined!");
937 } // switch(compEngine)
938
939}
940 break;
941 case 4:{
942#ifdef HAVE_INTREPID_DEBUG
943 TEUCHOS_TEST_FOR_EXCEPTION( (leftValues.rank() != 4 ), std::invalid_argument,
944 ">>> ERROR (ArrayTools::contractFieldFieldVector): Rank of the left input argument must equal 4!");
945 TEUCHOS_TEST_FOR_EXCEPTION( (rightValues.rank() != 4 ), std::invalid_argument,
946 ">>> ERROR (ArrayTools::contractFieldFieldVector): Rank of right input argument must equal 4!");
947 TEUCHOS_TEST_FOR_EXCEPTION( (outputValues.rank() != 3 ), std::invalid_argument,
948 ">>> ERROR (ArrayTools::contractFieldFieldVector): Rank of output argument must equal 3!");
949 TEUCHOS_TEST_FOR_EXCEPTION( (leftValues.dimension(0) != rightValues.dimension(0) ), std::invalid_argument,
950 ">>> ERROR (ArrayTools::contractFieldFieldVector): Zeroth dimensions (number of integration domains) of the left and right input containers must agree!");
951 TEUCHOS_TEST_FOR_EXCEPTION( (leftValues.dimension(2) != rightValues.dimension(2) ), std::invalid_argument,
952 ">>> ERROR (ArrayTools::contractFieldFieldVector): Second dimensions (numbers of integration points) of the left and right input containers must agree!");
953 TEUCHOS_TEST_FOR_EXCEPTION( (leftValues.dimension(3) != rightValues.dimension(3) ), std::invalid_argument,
954 ">>> ERROR (ArrayTools::contractFieldFieldVector): Third dimensions (numbers of vector components) of the left and right input containers must agree!");
955 TEUCHOS_TEST_FOR_EXCEPTION( (outputValues.dimension(0) != rightValues.dimension(0) ), std::invalid_argument,
956 ">>> ERROR (ArrayTools::contractFieldFieldVector): Zeroth dimensions (numbers of integration domains) of the input and output containers must agree!");
957 TEUCHOS_TEST_FOR_EXCEPTION( (outputValues.dimension(1) != leftValues.dimension(1) ), std::invalid_argument,
958 ">>> ERROR (ArrayTools::contractFieldFieldVector): First dimension of output container and first dimension of left input container must agree!");
959 TEUCHOS_TEST_FOR_EXCEPTION( (outputValues.dimension(2) != rightValues.dimension(1) ), std::invalid_argument,
960 ">>> ERROR (ArrayTools::contractFieldFieldVector): Second dimension of output container and first dimension of right input container must agree!");
961#endif
962
963 // get sizes
964 int numCells = leftValues.dimension(0);
965 int numLeftFields = leftValues.dimension(1);
966 int numRightFields = rightValues.dimension(1);
967 int numPoints = leftValues.dimension(2);
968 int dimVec = leftValues.dimension(3);
969
970 switch(compEngine) {
971 case COMP_CPP: {
972 if (sumInto) {
973 for (int cl = 0; cl < numCells; cl++) {
974 for (int lbf = 0; lbf < numLeftFields; lbf++) {
975 for (int rbf = 0; rbf < numRightFields; rbf++) {
976 Scalar tmpVal(0);
977 for (int qp = 0; qp < numPoints; qp++) {
978 for (int iVec = 0; iVec < dimVec; iVec++) {
979 tmpVal += leftValues(cl, lbf, qp, iVec)*rightValues(cl, rbf, qp, iVec);
980 } //D-loop
981 } // P-loop
982 outputValues(cl, lbf, rbf) += tmpVal;
983 } // R-loop
984 } // L-loop
985 } // C-loop
986 }
987 else {
988 for (int cl = 0; cl < numCells; cl++) {
989 for (int lbf = 0; lbf < numLeftFields; lbf++) {
990 for (int rbf = 0; rbf < numRightFields; rbf++) {
991 Scalar tmpVal(0);
992 for (int qp = 0; qp < numPoints; qp++) {
993 for (int iVec = 0; iVec < dimVec; iVec++) {
994 tmpVal += leftValues(cl, lbf, qp, iVec)*rightValues(cl, rbf, qp, iVec);
995 } //D-loop
996 } // P-loop
997 outputValues(cl, lbf, rbf) = tmpVal;
998 } // R-loop
999 } // L-loop
1000 } // C-loop
1001 }
1002 }
1003 break;
1004
1005 case COMP_BLAS: {
1006 /*
1007 GEMM parameters and their values.
1008 (Note: It is assumed that the result needs to be transposed into row-major format.
1009 Think of left and right input matrices as A(p x l) and B(p x r), respectively,
1010 even though the indexing is ((C),L,P) and ((C),R,P). Due to BLAS formatting
1011 assumptions, we are computing (A^T*B)^T = B^T*A.)
1012 TRANSA TRANS
1013 TRANSB NO_TRANS
1014 M #rows(B^T) = number of right fields
1015 N #cols(A) = number of left fields
1016 K #cols(B^T) = number of integration points * size of vector
1017 ALPHA 1.0
1018 A right data for cell cl = &rightFields[cl*skipR]
1019 LDA #rows(B) = number of integration points * size of vector
1020 B left data for cell cl = &leftFields[cl*skipL]
1021 LDB #rows(A) = number of integration points * size of vector
1022 BETA 0.0
1023 C result for cell cl = &outputFields[cl*skipOp]
1024 LDC #rows(C) = number of right fields
1025 */
1026 int numData = numPoints*dimVec;
1027 int skipL = numLeftFields*numData; // size of the left data chunk per cell
1028 int skipR = numRightFields*numData; // size of the right data chunk per cell
1029 int skipOp = numLeftFields*numRightFields; // size of the output data chunk per cell
1030 Scalar alpha(1.0); // these are left unchanged by GEMM
1031 Scalar beta(0.0);
1032 if (sumInto) {
1033 beta = 1.0;
1034 }
1035
1036 for (int cl=0; cl < numCells; cl++) {
1037 /* Use this if data is used in row-major format */
1038 Teuchos::BLAS<int, Scalar> myblas;
1039 myblas.GEMM(Teuchos::TRANS, Teuchos::NO_TRANS,
1040 numRightFields, numLeftFields, numData,
1041 alpha, &rightValues[cl*skipR], numData,
1042 &leftValues[cl*skipL], numData,
1043 beta, &outputValues[cl*skipOp], numRightFields);
1044 /* Use this if data is used in column-major format */
1045 /*
1046 myblas.GEMM(Teuchos::TRANS, Teuchos::NO_TRANS,
1047 numLeftFields, numRightFields, numData,
1048 alpha, &leftFields[cl*skipL], numData,
1049 &rightFields[cl*skipR], numData,
1050 beta, &outputFields[cl*skipOp], numLeftFields);
1051 */
1052 }
1053 }
1054 break;
1055
1056 default:
1057 TEUCHOS_TEST_FOR_EXCEPTION( ( ~isValidCompEngine(compEngine) ), std::invalid_argument,
1058 ">>> ERROR (ArrayTools::contractFieldFieldVector): Computational engine not defined!");
1059 } // switch(compEngine)
1060}
1061 break;
1062 case 5:{
1063#ifdef HAVE_INTREPID_DEBUG
1064 TEUCHOS_TEST_FOR_EXCEPTION( (leftValues.rank() != 5 ), std::invalid_argument,
1065 ">>> ERROR (ArrayTools::contractFieldFieldTensor): Rank of the left input argument must equal 5!");
1066 TEUCHOS_TEST_FOR_EXCEPTION( (rightValues.rank() != 5 ), std::invalid_argument,
1067 ">>> ERROR (ArrayTools::contractFieldFieldTensor): Rank of right input argument must equal 5!");
1068 TEUCHOS_TEST_FOR_EXCEPTION( (outputValues.rank() != 3 ), std::invalid_argument,
1069 ">>> ERROR (ArrayTools::contractFieldFieldTensor): Rank of output argument must equal 3!");
1070 TEUCHOS_TEST_FOR_EXCEPTION( (leftValues.dimension(0) != rightValues.dimension(0) ), std::invalid_argument,
1071 ">>> ERROR (ArrayTools::contractFieldFieldTensor): Zeroth dimensions (number of integration domains) of the left and right input containers must agree!");
1072 TEUCHOS_TEST_FOR_EXCEPTION( (leftValues.dimension(2) != rightValues.dimension(2) ), std::invalid_argument,
1073 ">>> ERROR (ArrayTools::contractFieldFieldTensor): Second dimensions (numbers of integration points) of the left and right input containers must agree!");
1074 TEUCHOS_TEST_FOR_EXCEPTION( (leftValues.dimension(3) != rightValues.dimension(3) ), std::invalid_argument,
1075 ">>> ERROR (ArrayTools::contractFieldFieldTensor): Third dimensions (first tensor dimensions) of the left and right input containers must agree!");
1076 TEUCHOS_TEST_FOR_EXCEPTION( (leftValues.dimension(4) != rightValues.dimension(4) ), std::invalid_argument,
1077 ">>> ERROR (ArrayTools::contractFieldFieldTensor): Fourth dimensions (second tensor dimensions) of the left and right input containers must agree!");
1078 TEUCHOS_TEST_FOR_EXCEPTION( (outputValues.dimension(0) != rightValues.dimension(0) ), std::invalid_argument,
1079 ">>> ERROR (ArrayTools::contractFieldFieldTensor): Zeroth dimensions (numbers of integration domains) of the input and output containers must agree!");
1080 TEUCHOS_TEST_FOR_EXCEPTION( (outputValues.dimension(1) != leftValues.dimension(1) ), std::invalid_argument,
1081 ">>> ERROR (ArrayTools::contractFieldFieldTensor): First dimension of output container and first dimension of left input container must agree!");
1082 TEUCHOS_TEST_FOR_EXCEPTION( (outputValues.dimension(2) != rightValues.dimension(1) ), std::invalid_argument,
1083 ">>> ERROR (ArrayTools::contractFieldFieldTensor): Second dimension of output container and first dimension of right input container must agree!");
1084#endif
1085
1086 // get sizes
1087 int numCells = leftValues.dimension(0);
1088 int numLeftFields = leftValues.dimension(1);
1089 int numRightFields = rightValues.dimension(1);
1090 int numPoints = leftValues.dimension(2);
1091 int dim1Tensor = leftValues.dimension(3);
1092 int dim2Tensor = leftValues.dimension(4);
1093
1094 switch(compEngine) {
1095 case COMP_CPP: {
1096 if (sumInto) {
1097 for (int cl = 0; cl < numCells; cl++) {
1098 for (int lbf = 0; lbf < numLeftFields; lbf++) {
1099 for (int rbf = 0; rbf < numRightFields; rbf++) {
1100 Scalar tmpVal(0);
1101 for (int qp = 0; qp < numPoints; qp++) {
1102 for (int iTens1 = 0; iTens1 < dim1Tensor; iTens1++) {
1103 for (int iTens2 = 0; iTens2 < dim2Tensor; iTens2++) {
1104 tmpVal += leftValues(cl, lbf, qp, iTens1, iTens2)*rightValues(cl, rbf, qp, iTens1, iTens2);
1105 } // D2-loop
1106 } // D1-loop
1107 } // P-loop
1108 outputValues(cl, lbf, rbf) += tmpVal;
1109 } // R-loop
1110 } // L-loop
1111 } // C-loop
1112 }
1113 else {
1114 for (int cl = 0; cl < numCells; cl++) {
1115 for (int lbf = 0; lbf < numLeftFields; lbf++) {
1116 for (int rbf = 0; rbf < numRightFields; rbf++) {
1117 Scalar tmpVal(0);
1118 for (int qp = 0; qp < numPoints; qp++) {
1119 for (int iTens1 = 0; iTens1 < dim1Tensor; iTens1++) {
1120 for (int iTens2 = 0; iTens2 < dim2Tensor; iTens2++) {
1121 tmpVal += leftValues(cl, lbf, qp, iTens1, iTens2)*rightValues(cl, rbf, qp, iTens1, iTens2);
1122 } // D2-loop
1123 } // D1-loop
1124 } // P-loop
1125 outputValues(cl, lbf, rbf) = tmpVal;
1126 } // R-loop
1127 } // L-loop
1128 } // C-loop
1129 }
1130 }
1131 break;
1132
1133 case COMP_BLAS: {
1134 /*
1135 GEMM parameters and their values.
1136 (Note: It is assumed that the result needs to be transposed into row-major format.
1137 Think of left and right input matrices as A(p x l) and B(p x r), respectively,
1138 even though the indexing is ((C),L,P) and ((C),R,P). Due to BLAS formatting
1139 assumptions, we are computing (A^T*B)^T = B^T*A.)
1140 TRANSA TRANS
1141 TRANSB NO_TRANS
1142 M #rows(B^T) = number of right fields
1143 N #cols(A) = number of left fields
1144 K #cols(B^T) = number of integration points * size of tensor
1145 ALPHA 1.0
1146 A right data for cell cl = &rightFields[cl*skipR]
1147 LDA #rows(B) = number of integration points * size of tensor
1148 B left data for cell cl = &leftFields[cl*skipL]
1149 LDB #rows(A) = number of integration points * size of tensor
1150 BETA 0.0
1151 C result for cell cl = &outputFields[cl*skipOp]
1152 LDC #rows(C) = number of right fields
1153 */
1154 int numData = numPoints*dim1Tensor*dim2Tensor;
1155 int skipL = numLeftFields*numData; // size of the left data chunk per cell
1156 int skipR = numRightFields*numData; // size of the right data chunk per cell
1157 int skipOp = numLeftFields*numRightFields; // size of the output data chunk per cell
1158 Scalar alpha(1.0); // these are left unchanged by GEMM
1159 Scalar beta(0.0);
1160 if (sumInto) {
1161 beta = 1.0;
1162 }
1163
1164 for (int cl=0; cl < numCells; cl++) {
1165 /* Use this if data is used in row-major format */
1166 Teuchos::BLAS<int, Scalar> myblas;
1167 myblas.GEMM(Teuchos::TRANS, Teuchos::NO_TRANS,
1168 numRightFields, numLeftFields, numData,
1169 alpha, &rightValues[cl*skipR], numData,
1170 &leftValues[cl*skipL], numData,
1171 beta, &outputValues[cl*skipOp], numRightFields);
1172 /* Use this if data is used in column-major format */
1173 /*
1174 myblas.GEMM(Teuchos::TRANS, Teuchos::NO_TRANS,
1175 numLeftFields, numRightFields, numData,
1176 alpha, &leftFields[cl*skipL], numData,
1177 &rightFields[cl*skipR], numData,
1178 beta, &outputFields[cl*skipOp], numLeftFields);
1179 */
1180 }
1181 }
1182 break;
1183
1184 default:
1185 TEUCHOS_TEST_FOR_EXCEPTION( ( ~isValidCompEngine(compEngine) ), std::invalid_argument,
1186 ">>> ERROR (ArrayTools::contractFieldFieldTensor): Computational engine not defined!");
1187 } // switch(compEngine)
1188}
1189 break;
1190 default:
1191#ifdef HAVE_INTREPID_DEBUG
1192
1193 TEUCHOS_TEST_FOR_EXCEPTION( ((lRank != 3) && (lRank != 4) && (lRank != 5)), std::invalid_argument,
1194 ">>> ERROR (FunctionSpaceTools::operatorIntegral): Left fields input container must have rank 3, 4 or 5.");
1195
1196#endif
1197 break;
1198}
1199
1200
1201
1202}
1203 break;
1204default:
1205#ifdef HAVE_INTREPID_DEBUG
1206
1207 TEUCHOS_TEST_FOR_EXCEPTION( ((outRank != 1) && (outRank != 2) && (outRank != 3)), std::invalid_argument,
1208 ">>> ERROR (FunctionSpaceTools::integrate): Output container must have rank 1, 2 or 3.");
1209
1210#endif
1211break;
1212}
1213}
1214template<class Scalar, class ArrayOut, class ArrayInLeft, class ArrayInRight>
1215void FunctionSpaceTools::integrate(ArrayOut & outputValues,
1216 const ArrayInLeft & leftValues,
1217 const ArrayInRight & rightValues,
1218 const ECompEngine compEngine,
1219 const bool sumInto) {
1220 ArrayWrapper<Scalar,ArrayOut, Rank<ArrayOut >::value, false>outputValuesWrap(outputValues);
1221 ArrayWrapper<Scalar,ArrayInLeft, Rank<ArrayInLeft >::value, true>leftValuesWrap(leftValues);
1222 ArrayWrapper<Scalar,ArrayInRight, Rank<ArrayInRight >::value, true>rightValuesWrap(rightValues);
1223 int outRank = getrank(outputValues);
1224 switch (outRank) {
1225 case 1:
1226 dataIntegral<Scalar>(outputValuesWrap, leftValuesWrap, rightValuesWrap, compEngine, sumInto);
1227 break;
1228 case 2:
1229 functionalIntegral<Scalar>(outputValuesWrap, leftValuesWrap, rightValuesWrap, compEngine, sumInto);
1230 break;
1231 case 3:
1232 operatorIntegral<Scalar>(outputValuesWrap, leftValuesWrap, rightValuesWrap, compEngine, sumInto);
1233 break;
1234 default:
1235#ifdef HAVE_INTREPID_DEBUG
1236
1237 TEUCHOS_TEST_FOR_EXCEPTION( ((outRank != 1) && (outRank != 2) && (outRank != 3)), std::invalid_argument,
1238 ">>> ERROR (FunctionSpaceTools::integrate): Output container must have rank 1, 2 or 3.");
1239
1240#endif
1241 break;
1242 }
1243
1244} // integrate
1245template<class Scalar, class ArrayOutFields, class ArrayInFieldsLeft, class ArrayInFieldsRight>
1246void FunctionSpaceTools::operatorIntegral(ArrayOutFields & outputFields,
1247 const ArrayInFieldsLeft & leftFields,
1248 const ArrayInFieldsRight & rightFields,
1249 const ECompEngine compEngine,
1250 const bool sumInto) {
1251 int lRank = getrank(leftFields);
1252
1253 switch (lRank) {
1254 case 3:
1255 ArrayTools::contractFieldFieldScalar<Scalar>(outputFields, leftFields, rightFields, compEngine, sumInto);
1256 break;
1257 case 4:
1258 ArrayTools::contractFieldFieldVector<Scalar>(outputFields, leftFields, rightFields, compEngine, sumInto);
1259 break;
1260 case 5:
1261 ArrayTools::contractFieldFieldTensor<Scalar>(outputFields, leftFields, rightFields, compEngine, sumInto);
1262 break;
1263 default:
1264#ifdef HAVE_INTREPID_DEBUG
1265
1266 TEUCHOS_TEST_FOR_EXCEPTION( ((lRank != 3) && (lRank != 4) && (lRank != 5)), std::invalid_argument,
1267 ">>> ERROR (FunctionSpaceTools::operatorIntegral): Left fields input container must have rank 3, 4 or 5.");
1268
1269#endif
1270 break;
1271}
1272
1273} // operatorIntegral
1274
1275
1276template<class Scalar, class ArrayOutFields, class ArrayInData, class ArrayInFields>
1277void FunctionSpaceTools::functionalIntegral(ArrayOutFields & outputFields,
1278 const ArrayInData & inputData,
1279 const ArrayInFields & inputFields,
1280 const ECompEngine compEngine,
1281 const bool sumInto) {
1282 int dRank = getrank(inputData);
1283
1284 switch (dRank) {
1285 case 2:
1286 ArrayTools::contractDataFieldScalar<Scalar>(outputFields, inputData, inputFields, compEngine, sumInto);
1287 break;
1288 case 3:
1289 ArrayTools::contractDataFieldVector<Scalar>(outputFields, inputData, inputFields, compEngine, sumInto);
1290 break;
1291 case 4:
1292 ArrayTools::contractDataFieldTensor<Scalar>(outputFields, inputData, inputFields, compEngine, sumInto);
1293 break;
1294 default:
1295#ifdef HAVE_INTREPID_DEBUG
1296
1297TEUCHOS_TEST_FOR_EXCEPTION( ((dRank != 2) && (dRank != 3) && (dRank != 4)), std::invalid_argument,
1298 ">>> ERROR (FunctionSpaceTools::functionalIntegral): Data input container must have rank 2, 3 or 4.");
1299
1300#endif
1301 break;
1302}
1303
1304} // functionalIntegral
1305
1306template<class Scalar, class ArrayOutData, class ArrayInDataLeft, class ArrayInDataRight>
1307void FunctionSpaceTools::dataIntegral(ArrayOutData & outputData,
1308 const ArrayInDataLeft & inputDataLeft,
1309 const ArrayInDataRight & inputDataRight,
1310 const ECompEngine compEngine,
1311 const bool sumInto) {
1312 int lRank = getrank(inputDataLeft);
1313
1314 switch (lRank) {
1315 case 2:
1316 ArrayTools::contractDataDataScalar<Scalar>(outputData, inputDataLeft, inputDataRight, compEngine, sumInto);
1317 break;
1318 case 3:
1319 ArrayTools::contractDataDataVector<Scalar>(outputData, inputDataLeft, inputDataRight, compEngine, sumInto);
1320 break;
1321 case 4:
1322 ArrayTools::contractDataDataTensor<Scalar>(outputData, inputDataLeft, inputDataRight, compEngine, sumInto);
1323 break;
1324 default:
1325#ifdef HAVE_INTREPID_DEBUG
1326
1327 TEUCHOS_TEST_FOR_EXCEPTION( ((lRank != 2) && (lRank != 3) && (lRank != 4)), std::invalid_argument,
1328 ">>> ERROR (FunctionSpaceTools::dataIntegral): Left data input container must have rank 2, 3 or 4.");
1329
1330#endif
1331 break;
1332}
1333
1334} // dataIntegral
1335
1336template<class Scalar, class ArrayOut, class ArrayDet, class ArrayWeights>
1337inline void FunctionSpaceTools::computeCellMeasure(ArrayOut & outVals,
1338 const ArrayDet & inDet,
1339 const ArrayWeights & inWeights) {
1340#ifdef HAVE_INTREPID_DEBUG
1341
1342 TEUCHOS_TEST_FOR_EXCEPTION( (inDet.rank() != 2), std::invalid_argument,
1343 ">>> ERROR (FunctionSpaceTools::computeCellMeasure): Input determinants container must have rank 2.");
1344
1345#endif
1346
1347 ArrayTools::scalarMultiplyDataData<Scalar>(outVals, inDet, inWeights);
1348 // must use absolute value of inDet, so flip sign where needed
1349 for (int cell=0; cell<outVals.dimension(0); cell++) {
1350 if (inDet(cell,0) < 0.0) {
1351 for (int point=0; point<outVals.dimension(1); point++) {
1352 outVals(cell, point) *= -1.0;
1353 }
1354 }
1355 }
1356
1357} // computeCellMeasure
1358
1359template<class Scalar, class ArrayOut, class ArrayJac, class ArrayWeights>
1361 const ArrayJac & inJac,
1362 const ArrayWeights & inWeights,
1363 const int whichFace,
1364 const shards::CellTopology & parentCell) {
1365
1366#ifdef HAVE_INTREPID_DEBUG
1367
1368 TEUCHOS_TEST_FOR_EXCEPTION( (inJac.rank() != 4), std::invalid_argument,
1369 ">>> ERROR (FunctionSpaceTools::computeFaceMeasure): Input Jacobian container must have rank 4.");
1370
1371#endif
1372
1373 // temporary storage for face normals
1374 FieldContainer<Scalar> faceNormals(inJac.dimension(0), inJac.dimension(1), inJac.dimension(2));
1375
1376 // compute normals
1377 CellTools<Scalar>::getPhysicalFaceNormals(faceNormals, inJac, whichFace, parentCell);
1378
1379 // compute lenghts of normals
1380 RealSpaceTools<Scalar>::vectorNorm(outVals, faceNormals, NORM_TWO);
1381
1382 // multiply with weights
1383 ArrayTools::scalarMultiplyDataData<Scalar>(outVals, outVals, inWeights);
1384
1385}
1386
1387
1388template<class Scalar, class ArrayOut, class ArrayJac, class ArrayWeights>
1390 const ArrayJac & inJac,
1391 const ArrayWeights & inWeights,
1392 const int whichEdge,
1393 const shards::CellTopology & parentCell) {
1394
1395#ifdef HAVE_INTREPID_DEBUG
1396
1397 TEUCHOS_TEST_FOR_EXCEPTION( (inJac.rank() != 4), std::invalid_argument,
1398 ">>> ERROR (FunctionSpaceTools::computeEdgeMeasure): Input Jacobian container must have rank 4.");
1399
1400#endif
1401
1402 // temporary storage for edge tangents
1403 FieldContainer<Scalar> edgeTangents(inJac.dimension(0), inJac.dimension(1), inJac.dimension(2));
1404
1405 // compute normals
1406 CellTools<Scalar>::getPhysicalEdgeTangents(edgeTangents, inJac, whichEdge, parentCell);
1407
1408 // compute lenghts of tangents
1409 RealSpaceTools<Scalar>::vectorNorm(outVals, edgeTangents, NORM_TWO);
1410
1411 // multiply with weights
1412 ArrayTools::scalarMultiplyDataData<Scalar>(outVals, outVals, inWeights);
1413
1414}
1415
1416
1417
1418template<class Scalar, class ArrayTypeOut, class ArrayTypeMeasure, class ArrayTypeIn>
1419void FunctionSpaceTools::multiplyMeasure(ArrayTypeOut & outVals,
1420 const ArrayTypeMeasure & inMeasure,
1421 const ArrayTypeIn & inVals) {
1422
1423 ArrayTools::scalarMultiplyDataField<Scalar>(outVals, inMeasure, inVals);
1424
1425} // multiplyMeasure
1426
1427
1428template<class Scalar, class ArrayOutFields, class ArrayInData, class ArrayInFields>
1429void FunctionSpaceTools::scalarMultiplyDataField(ArrayOutFields & outputFields,
1430 ArrayInData & inputData,
1431 ArrayInFields & inputFields,
1432 const bool reciprocal) {
1433
1434 ArrayTools::scalarMultiplyDataField<Scalar>(outputFields, inputData, inputFields, reciprocal);
1435
1436} // scalarMultiplyDataField
1437
1438
1439template<class Scalar, class ArrayOutData, class ArrayInDataLeft, class ArrayInDataRight>
1440void FunctionSpaceTools::scalarMultiplyDataData(ArrayOutData & outputData,
1441 ArrayInDataLeft & inputDataLeft,
1442 ArrayInDataRight & inputDataRight,
1443 const bool reciprocal) {
1444
1445 ArrayTools::scalarMultiplyDataData<Scalar>(outputData, inputDataLeft, inputDataRight, reciprocal);
1446
1447} // scalarMultiplyDataData
1448
1449
1450template<class Scalar, class ArrayOutFields, class ArrayInData, class ArrayInFields>
1451void FunctionSpaceTools::dotMultiplyDataField(ArrayOutFields & outputFields,
1452 const ArrayInData & inputData,
1453 const ArrayInFields & inputFields) {
1454
1455 ArrayTools::dotMultiplyDataField<Scalar>(outputFields, inputData, inputFields);
1456
1457} // dotMultiplyDataField
1458
1459
1460template<class Scalar, class ArrayOutData, class ArrayInDataLeft, class ArrayInDataRight>
1461void FunctionSpaceTools::dotMultiplyDataData(ArrayOutData & outputData,
1462 const ArrayInDataLeft & inputDataLeft,
1463 const ArrayInDataRight & inputDataRight) {
1464
1465 ArrayTools::dotMultiplyDataData<Scalar>(outputData, inputDataLeft, inputDataRight);
1466
1467} // dotMultiplyDataData
1468
1469
1470template<class Scalar, class ArrayOutFields, class ArrayInData, class ArrayInFields>
1471void FunctionSpaceTools::vectorMultiplyDataField(ArrayOutFields & outputFields,
1472 const ArrayInData & inputData,
1473 const ArrayInFields & inputFields) {
1474
1475 int outRank = getrank(outputFields);
1476
1477 switch (outRank) {
1478 case 3:
1479 case 4:
1480 ArrayTools::crossProductDataField<Scalar>(outputFields, inputData, inputFields);
1481 break;
1482 case 5:
1483 ArrayTools::outerProductDataField<Scalar>(outputFields, inputData, inputFields);
1484 break;
1485 default:
1486 TEUCHOS_TEST_FOR_EXCEPTION( ((outRank != 3) && (outRank != 4) && (outRank != 5)), std::invalid_argument,
1487 ">>> ERROR (FunctionSpaceTools::vectorMultiplyDataField): Output container must have rank 3, 4 or 5.");
1488 }
1489
1490} // vectorMultiplyDataField
1491
1492
1493template<class Scalar, class ArrayOutData, class ArrayInDataLeft, class ArrayInDataRight>
1494void FunctionSpaceTools::vectorMultiplyDataData(ArrayOutData & outputData,
1495 const ArrayInDataLeft & inputDataLeft,
1496 const ArrayInDataRight & inputDataRight) {
1497
1498 int outRank = getrank(outputData);
1499
1500 switch (outRank) {
1501 case 2:
1502 case 3:
1503 ArrayTools::crossProductDataData<Scalar>(outputData, inputDataLeft, inputDataRight);
1504 break;
1505 case 4:
1506 ArrayTools::outerProductDataData<Scalar>(outputData, inputDataLeft, inputDataRight);
1507 break;
1508 default:
1509 TEUCHOS_TEST_FOR_EXCEPTION( ((outRank != 2) && (outRank != 3) && (outRank != 4)), std::invalid_argument,
1510 ">>> ERROR (FunctionSpaceTools::vectorMultiplyDataData): Output container must have rank 2, 3 or 4.");
1511 }
1512
1513} // vectorMultiplyDataData
1514
1515
1516template<class Scalar, class ArrayOutFields, class ArrayInData, class ArrayInFields>
1517void FunctionSpaceTools::tensorMultiplyDataField(ArrayOutFields & outputFields,
1518 const ArrayInData & inputData,
1519 const ArrayInFields & inputFields,
1520 const char transpose) {
1521
1522 int outRank = outputFields.rank();
1523
1524 switch (outRank) {
1525 case 4:
1526 ArrayTools::matvecProductDataField<Scalar>(outputFields, inputData, inputFields, transpose);
1527 break;
1528 case 5:
1529 ArrayTools::matmatProductDataField<Scalar>(outputFields, inputData, inputFields, transpose);
1530 break;
1531 default:
1532 TEUCHOS_TEST_FOR_EXCEPTION( ((outRank != 4) && (outRank != 5)), std::invalid_argument,
1533 ">>> ERROR (FunctionSpaceTools::tensorMultiplyDataField): Output container must have rank 4 or 5.");
1534 }
1535
1536} // tensorMultiplyDataField
1537
1538template<class Scalar, class ArrayOutData, class ArrayInDataLeft, class ArrayInDataRight>
1539struct FunctionSpaceTools::tensorMultiplyDataDataTempSpec<Scalar, ArrayOutData, ArrayInDataLeft, ArrayInDataRight,-1>{
1540 tensorMultiplyDataDataTempSpec(ArrayOutData & outputData,
1541 const ArrayInDataLeft & inputDataLeft,
1542 const ArrayInDataRight & inputDataRight,
1543 const char transpose) {
1544 int outRank = getrank(outputData);
1545
1546 switch (outRank) {
1547 case 3:
1548 ArrayTools::matvecProductDataData<Scalar>(outputData, inputDataLeft, inputDataRight, transpose);
1549 break;
1550 case 4:
1551 ArrayTools::matmatProductDataData<Scalar>(outputData, inputDataLeft, inputDataRight, transpose);
1552 break;
1553 }
1554}
1555};
1556template<class Scalar, class ArrayOutData, class ArrayInDataLeft, class ArrayInDataRight>
1557struct FunctionSpaceTools::tensorMultiplyDataDataTempSpec<Scalar, ArrayOutData, ArrayInDataLeft, ArrayInDataRight,3>{
1558 tensorMultiplyDataDataTempSpec(ArrayOutData & outputData,
1559 const ArrayInDataLeft & inputDataLeft,
1560 const ArrayInDataRight & inputDataRight,
1561 const char transpose) {
1562 ArrayTools::matvecProductDataData<Scalar>(outputData, inputDataLeft, inputDataRight, transpose);
1563
1564}
1565};
1566template<class Scalar, class ArrayOutData, class ArrayInDataLeft, class ArrayInDataRight>
1567struct FunctionSpaceTools::tensorMultiplyDataDataTempSpec<Scalar, ArrayOutData, ArrayInDataLeft, ArrayInDataRight,4>{
1568 tensorMultiplyDataDataTempSpec(ArrayOutData & outputData,
1569 const ArrayInDataLeft & inputDataLeft,
1570 const ArrayInDataRight & inputDataRight,
1571 const char transpose) {
1572 ArrayTools::matmatProductDataData<Scalar>(outputData, inputDataLeft, inputDataRight, transpose);
1573}
1574
1575};
1576template<class Scalar, class ArrayOutData, class ArrayInDataLeft, class ArrayInDataRight>
1577 void FunctionSpaceTools::tensorMultiplyDataData(ArrayOutData & outputData,
1578 const ArrayInDataLeft & inputDataLeft,
1579 const ArrayInDataRight & inputDataRight,
1580 const char transpose){
1582
1583 }
1584template<class Scalar, class ArrayTypeInOut, class ArrayTypeSign>
1585void FunctionSpaceTools::applyLeftFieldSigns(ArrayTypeInOut & inoutOperator,
1586 const ArrayTypeSign & fieldSigns) {
1587#ifdef HAVE_INTREPID_DEBUG
1588 TEUCHOS_TEST_FOR_EXCEPTION( (inoutOperator.rank() != 3), std::invalid_argument,
1589 ">>> ERROR (FunctionSpaceTools::applyLeftFieldSigns): Input operator container must have rank 3.");
1590 TEUCHOS_TEST_FOR_EXCEPTION( (fieldSigns.rank() != 2), std::invalid_argument,
1591 ">>> ERROR (FunctionSpaceTools::applyLeftFieldSigns): Input field signs container must have rank 2.");
1592 TEUCHOS_TEST_FOR_EXCEPTION( (inoutOperator.dimension(0) != fieldSigns.dimension(0) ), std::invalid_argument,
1593 ">>> ERROR (FunctionSpaceTools::applyLeftFieldSigns): Zeroth dimensions (number of cells) of the operator and field signs containers must agree!");
1594 TEUCHOS_TEST_FOR_EXCEPTION( (inoutOperator.dimension(1) != fieldSigns.dimension(1) ), std::invalid_argument,
1595 ">>> ERROR (FunctionSpaceTools::applyLeftFieldSigns): First dimensions (number of left fields) of the operator and field signs containers must agree!");
1596#endif
1597
1598 for (int cell=0; cell<inoutOperator.dimension(0); cell++) {
1599 for (int lbf=0; lbf<inoutOperator.dimension(1); lbf++) {
1600 for (int rbf=0; rbf<inoutOperator.dimension(2); rbf++) {
1601 inoutOperator(cell, lbf, rbf) *= fieldSigns(cell, lbf);
1602 }
1603 }
1604 }
1605
1606} // applyLeftFieldSigns
1607
1608
1609template<class Scalar, class ArrayTypeInOut, class ArrayTypeSign>
1610void FunctionSpaceTools::applyRightFieldSigns(ArrayTypeInOut & inoutOperator,
1611 const ArrayTypeSign & fieldSigns) {
1612#ifdef HAVE_INTREPID_DEBUG
1613 TEUCHOS_TEST_FOR_EXCEPTION( (getrank(inoutOperator) != 3), std::invalid_argument,
1614 ">>> ERROR (FunctionSpaceTools::applyRightFieldSigns): Input operator container must have rank 3.");
1615 TEUCHOS_TEST_FOR_EXCEPTION( (getrank(fieldSigns) != 2), std::invalid_argument,
1616 ">>> ERROR (FunctionSpaceTools::applyRightFieldSigns): Input field signs container must have rank 2.");
1617 TEUCHOS_TEST_FOR_EXCEPTION( (inoutOperator.dimension(0) != fieldSigns.dimension(0) ), std::invalid_argument,
1618 ">>> ERROR (FunctionSpaceTools::applyRightFieldSigns): Zeroth dimensions (number of cells) of the operator and field signs containers must agree!");
1619 TEUCHOS_TEST_FOR_EXCEPTION( (inoutOperator.dimension(2) != fieldSigns.dimension(1) ), std::invalid_argument,
1620 ">>> ERROR (FunctionSpaceTools::applyRightFieldSigns): Second dimension of the operator container and first dimension of the field signs container (number of right fields) must agree!");
1621#endif
1622
1623 for (int cell=0; cell<inoutOperator.dimension(0); cell++) {
1624 for (int lbf=0; lbf<inoutOperator.dimension(1); lbf++) {
1625 for (int rbf=0; rbf<inoutOperator.dimension(2); rbf++) {
1626 inoutOperator(cell, lbf, rbf) *= fieldSigns(cell, rbf);
1627 }
1628 }
1629 }
1630
1631} // applyRightFieldSigns
1632
1633
1634
1635template<class Scalar, class ArrayTypeInOut, class ArrayTypeSign>
1636void FunctionSpaceTools::applyFieldSigns(ArrayTypeInOut & inoutFunction,
1637 const ArrayTypeSign & fieldSigns) {
1638
1639#ifdef HAVE_INTREPID_DEBUG
1640
1641 TEUCHOS_TEST_FOR_EXCEPTION( ((inoutFunction.rank() < 2) || (inoutFunction.rank() > 5)), std::invalid_argument,
1642 ">>> ERROR (FunctionSpaceTools::applyFieldSigns): Input function container must have rank 2, 3, 4, or 5.");
1643 TEUCHOS_TEST_FOR_EXCEPTION( (fieldSigns.rank() != 2), std::invalid_argument,
1644 ">>> ERROR (FunctionSpaceTools::applyFieldSigns): Input field signs container must have rank 2.");
1645 TEUCHOS_TEST_FOR_EXCEPTION( (inoutFunction.dimension(0) != fieldSigns.dimension(0) ), std::invalid_argument,
1646 ">>> ERROR (FunctionSpaceTools::applyFieldSigns): Zeroth dimensions (number of integration domains) of the function and field signs containers must agree!");
1647 TEUCHOS_TEST_FOR_EXCEPTION( (inoutFunction.dimension(1) != fieldSigns.dimension(1) ), std::invalid_argument,
1648 ">>> ERROR (FunctionSpaceTools::applyFieldSigns): First dimensions (number of fields) of the function and field signs containers must agree!");
1649
1650#endif
1651
1652 ArrayWrapper<Scalar,ArrayTypeInOut, Rank<ArrayTypeInOut >::value, false>inoutFunctionWrap(inoutFunction);
1653 ArrayWrapper<Scalar,ArrayTypeSign, Rank<ArrayTypeSign >::value, true>fieldSignsWrap(fieldSigns);
1654
1655
1656 int numCells = inoutFunction.dimension(0);
1657 int numFields = inoutFunction.dimension(1);
1658 int fRank = getrank(inoutFunction);
1659
1660 switch (fRank) {
1661 case 2: {
1662 for (int cell=0; cell<numCells; cell++) {
1663 for (int bf=0; bf<numFields; bf++) {
1664 inoutFunctionWrap(cell, bf) *= fieldSignsWrap(cell, bf);
1665 }
1666 }
1667 }
1668 break;
1669
1670 case 3: {
1671 int numPoints = inoutFunction.dimension(2);
1672 for (int cell=0; cell<numCells; cell++) {
1673 for (int bf=0; bf<numFields; bf++) {
1674 for (int pt=0; pt<numPoints; pt++) {
1675 inoutFunctionWrap(cell, bf, pt) *= fieldSignsWrap(cell, bf);
1676 }
1677 }
1678 }
1679 }
1680 break;
1681
1682 case 4: {
1683 int numPoints = inoutFunction.dimension(2);
1684 int spaceDim1 = inoutFunction.dimension(3);
1685 for (int cell=0; cell<numCells; cell++) {
1686 for (int bf=0; bf<numFields; bf++) {
1687 for (int pt=0; pt<numPoints; pt++) {
1688 for (int d1=0; d1<spaceDim1; d1++) {
1689 inoutFunctionWrap(cell, bf, pt, d1) *= fieldSignsWrap(cell, bf);
1690 }
1691 }
1692 }
1693 }
1694 }
1695 break;
1696
1697 case 5: {
1698 int numPoints = inoutFunction.dimension(2);
1699 int spaceDim1 = inoutFunction.dimension(3);
1700 int spaceDim2 = inoutFunction.dimension(4);
1701 for (int cell=0; cell<numCells; cell++) {
1702 for (int bf=0; bf<numFields; bf++) {
1703 for (int pt=0; pt<numPoints; pt++) {
1704 for (int d1=0; d1<spaceDim1; d1++) {
1705 for (int d2=0; d2<spaceDim2; d2++) {
1706 inoutFunctionWrap(cell, bf, pt, d1, d2) *= fieldSignsWrap(cell, bf);
1707 }
1708 }
1709 }
1710 }
1711 }
1712 }
1713 break;
1714
1715 default:
1716#ifdef HAVE_INTREPID_DEBUG
1717
1718 TEUCHOS_TEST_FOR_EXCEPTION( !( (inoutFunction.rank() == 2) || (inoutFunction.rank() == 3) || (inoutFunction.rank() == 4) || (inoutFunction.rank() == 5)), std::invalid_argument,
1719 ">>> ERROR (FunctionSpaceTools::applyFieldSigns): Method defined only for rank-2, 3, 4, or 5 input function containers.");
1720#endif
1721 break;
1722 } // end switch fRank
1723
1724
1725
1726} // applyFieldSigns
1727
1728template<class Scalar, class ArrayOutPointVals, class ArrayInCoeffs, class ArrayInFields>
1729void FunctionSpaceTools::evaluate(ArrayOutPointVals & outPointVals,
1730 const ArrayInCoeffs & inCoeffs,
1731 const ArrayInFields & inFields) {
1732
1733#ifdef HAVE_INTREPID_DEBUG
1734
1735 TEUCHOS_TEST_FOR_EXCEPTION( ((inFields.rank() < 3) || (inFields.rank() > 5)), std::invalid_argument,
1736 ">>> ERROR (FunctionSpaceTools::evaluate): Input fields container must have rank 3, 4, or 5.");
1737 TEUCHOS_TEST_FOR_EXCEPTION( (inCoeffs.rank() != 2), std::invalid_argument,
1738 ">>> ERROR (FunctionSpaceTools::evaluate): Input coefficient container must have rank 2.");
1739 TEUCHOS_TEST_FOR_EXCEPTION( (outPointVals.rank() != inFields.rank()-1), std::invalid_argument,
1740 ">>> ERROR (FunctionSpaceTools::evaluate): Output values container must have rank one less than the rank of the input fields container.");
1741 TEUCHOS_TEST_FOR_EXCEPTION( (inCoeffs.dimension(0) != inFields.dimension(0) ), std::invalid_argument,
1742 ">>> ERROR (FunctionSpaceTools::evaluate): Zeroth dimensions (number of cells) of the coefficient and fields input containers must agree!");
1743 TEUCHOS_TEST_FOR_EXCEPTION( (inCoeffs.dimension(1) != inFields.dimension(1) ), std::invalid_argument,
1744 ">>> ERROR (FunctionSpaceTools::evaluate): First dimensions (number of fields) of the coefficient and fields input containers must agree!");
1745 TEUCHOS_TEST_FOR_EXCEPTION( (outPointVals.dimension(0) != inFields.dimension(0) ), std::invalid_argument,
1746 ">>> ERROR (FunctionSpaceTools::evaluate): Zeroth dimensions (number of cells) of the input fields container and the output values container must agree!");
1747 for (int i=1; i<outPointVals.rank(); i++) {
1748 std::string errmsg = ">>> ERROR (FunctionSpaceTools::evaluate): Dimensions ";
1749 errmsg += (char)(48+i);
1750 errmsg += " and ";
1751 errmsg += (char)(48+i+1);
1752 errmsg += " of the output values and input fields containers must agree!";
1753 TEUCHOS_TEST_FOR_EXCEPTION( (outPointVals.dimension(i) != inFields.dimension(i+1)), std::invalid_argument, errmsg );
1754 }
1755
1756#endif
1757 ArrayWrapper<Scalar,ArrayOutPointVals, Rank<ArrayOutPointVals >::value, false>outPointValsWrap(outPointVals);
1758 ArrayWrapper<Scalar,ArrayInCoeffs, Rank<ArrayInCoeffs>::value, true>inCoeffsWrap(inCoeffs);
1759 ArrayWrapper<Scalar,ArrayInFields, Rank<ArrayInFields>::value, true>inFieldsWrap(inFields);
1760
1761 int numCells = inFields.dimension(0);
1762 int numFields = inFields.dimension(1);
1763 int numPoints = inFields.dimension(2);
1764 int fRank = getrank(inFields);
1765
1766 switch (fRank) {
1767 case 3: {
1768 for (int cell=0; cell<numCells; cell++) {
1769 for (int pt=0; pt<numPoints; pt++) {
1770 for (int bf=0; bf<numFields; bf++) {
1771 outPointValsWrap(cell, pt) += inCoeffsWrap(cell, bf) * inFieldsWrap(cell, bf, pt);
1772 }
1773 }
1774 }
1775 }
1776 break;
1777
1778 case 4: {
1779 int spaceDim1 = inFields.dimension(3);
1780 for (int cell=0; cell<numCells; cell++) {
1781 for (int pt=0; pt<numPoints; pt++) {
1782 for (int d1=0; d1<spaceDim1; d1++) {
1783 for (int bf=0; bf<numFields; bf++) {
1784 outPointValsWrap(cell, pt, d1) += inCoeffsWrap(cell, bf) * inFieldsWrap(cell, bf, pt, d1);
1785 }
1786 }
1787 }
1788 }
1789 }
1790 break;
1791
1792 case 5: {
1793 int spaceDim1 = inFields.dimension(3);
1794 int spaceDim2 = inFields.dimension(4);
1795 for (int cell=0; cell<numCells; cell++) {
1796 for (int pt=0; pt<numPoints; pt++) {
1797 for (int d1=0; d1<spaceDim1; d1++) {
1798 for (int d2=0; d2<spaceDim2; d2++) {
1799 for (int bf=0; bf<numFields; bf++) {
1800 outPointValsWrap(cell, pt, d1, d2) += inCoeffsWrap(cell, bf) * inFieldsWrap(cell, bf, pt, d1, d2);
1801 }
1802 }
1803 }
1804 }
1805 }
1806 }
1807 break;
1808
1809 default:
1810#ifdef HAVE_INTREPID_DEBUG
1811
1812 TEUCHOS_TEST_FOR_EXCEPTION( !( (fRank == 3) || (fRank == 4) || (fRank == 5)), std::invalid_argument,
1813 ">>> ERROR (FunctionSpaceTools::evaluate): Method defined only for rank-3, 4, or 5 input fields containers.");
1814
1815#endif
1816 break;
1817 } // end switch fRank
1818
1819} // evaluate
1820
1821} // end namespace Intrepid
int isValidCompEngine(const ECompEngine compEngType)
Verifies validity of a computational engine enum.
static void getPhysicalFaceNormals(ArrayFaceNormal &faceNormals, const ArrayJac &worksetJacobians, const int &worksetFaceOrd, const shards::CellTopology &parentCell)
Computes non-normalized normal vectors to physical faces in a face workset ; (see Subcell worksets fo...
static void getPhysicalEdgeTangents(ArrayEdgeTangent &edgeTangents, const ArrayJac &worksetJacobians, const int &worksetEdgeOrd, const shards::CellTopology &parentCell)
Computes non-normalized tangent vectors to physical edges in an edge workset ; (see Subcell worksets ...
Implementation of a templated lexicographical container for a multi-indexed scalar quantity....
static void computeFaceMeasure(ArrayOut &outVals, const ArrayJac &inJac, const ArrayWeights &inWeights, const int whichFace, const shards::CellTopology &parentCell)
Returns the weighted integration measures outVals with dimensions (C,P) used for the computation of f...
static void vectorMultiplyDataField(ArrayOutFields &outputFields, const ArrayInData &inputData, const ArrayInFields &inputFields)
Cross or outer product of data and fields; please read the description below.
static void integrate(Intrepid::FieldContainer< Scalar > &outputValues, const Intrepid::FieldContainer< Scalar > &leftValues, const Intrepid::FieldContainer< Scalar > &rightValues, const ECompEngine compEngine, const bool sumInto=false)
Contracts leftValues and rightValues arrays on the point and possibly space dimensions and stores the...
static void HDIVtransformDIV(ArrayTypeOut &outVals, const ArrayTypeDet &jacobianDet, const ArrayTypeIn &inVals)
Transformation of a divergence field in the H-div space, defined at points on a reference cell,...
static void HGRADtransformVALUE(ArrayTypeOut &outVals, const ArrayTypeIn &inVals)
Transformation of a (scalar) value field in the H-grad space, defined at points on a reference cell,...
static void tensorMultiplyDataData(ArrayOutData &outputData, const ArrayInDataLeft &inputDataLeft, const ArrayInDataRight &inputDataRight, const char transpose='N')
Matrix-vector or matrix-matrix product of data and data; please read the description below.
static void HVOLtransformVALUE(ArrayTypeOut &outVals, const ArrayTypeDet &jacobianDet, const ArrayTypeIn &inVals)
Transformation of a (scalar) value field in the H-vol space, defined at points on a reference cell,...
static void scalarMultiplyDataField(ArrayOutFields &outputFields, ArrayInData &inputData, ArrayInFields &inputFields, const bool reciprocal=false)
Scalar multiplication of data and fields; please read the description below.
static void HCURLtransformVALUE(ArrayTypeOut &outVals, const ArrayTypeJac &jacobianInverse, const ArrayTypeIn &inVals, const char transpose='T')
Transformation of a (vector) value field in the H-curl space, defined at points on a reference cell,...
static void applyRightFieldSigns(ArrayTypeInOut &inoutOperator, const ArrayTypeSign &fieldSigns)
Applies right (column) signs, stored in the user-provided container fieldSigns and indexed by (C,...
static void scalarMultiplyDataData(ArrayOutData &outputData, ArrayInDataLeft &inputDataLeft, ArrayInDataRight &inputDataRight, const bool reciprocal=false)
Scalar multiplication of data and data; please read the description below.
static void dotMultiplyDataField(ArrayOutFields &outputFields, const ArrayInData &inputData, const ArrayInFields &inputFields)
Dot product of data and fields; please read the description below.
static void evaluate(ArrayOutPointVals &outPointVals, const ArrayInCoeffs &inCoeffs, const ArrayInFields &inFields)
Computes point values outPointVals of a discrete function specified by the basis inFields and coeffic...
static void dotMultiplyDataData(ArrayOutData &outputData, const ArrayInDataLeft &inputDataLeft, const ArrayInDataRight &inputDataRight)
Dot product of data and data; please read the description below.
static void multiplyMeasure(ArrayTypeOut &outVals, const ArrayTypeMeasure &inMeasure, const ArrayTypeIn &inVals)
Multiplies fields inVals by weighted measures inMeasure and returns the field array outVals; this is ...
static void applyLeftFieldSigns(ArrayTypeInOut &inoutOperator, const ArrayTypeSign &fieldSigns)
Applies left (row) signs, stored in the user-provided container fieldSigns and indexed by (C,...
static void vectorMultiplyDataData(ArrayOutData &outputData, const ArrayInDataLeft &inputDataLeft, const ArrayInDataRight &inputDataRight)
Cross or outer product of data and data; please read the description below.
static void HGRADtransformGRAD(ArrayTypeOut &outVals, const ArrayTypeJac &jacobianInverse, const ArrayTypeIn &inVals, const char transpose='T')
Transformation of a gradient field in the H-grad space, defined at points on a reference cell,...
static void computeCellMeasure(ArrayOut &outVals, const ArrayDet &inDet, const ArrayWeights &inWeights)
Returns the weighted integration measures outVals with dimensions (C,P) used for the computation of c...
static void applyFieldSigns(ArrayTypeInOut &inoutFunction, const ArrayTypeSign &fieldSigns)
Applies field signs, stored in the user-provided container fieldSigns and indexed by (C,...
static void operatorIntegral(ArrayOutFields &outputFields, const ArrayInFieldsLeft &leftFields, const ArrayInFieldsRight &rightFields, const ECompEngine compEngine, const bool sumInto=false)
Contracts the point (and space) dimensions P (and D1 and D2) of two rank-3, 4, or 5 containers with d...
static void functionalIntegral(ArrayOutFields &outputFields, const ArrayInData &inputData, const ArrayInFields &inputFields, const ECompEngine compEngine, const bool sumInto=false)
Contracts the point (and space) dimensions P (and D1 and D2) of a rank-3, 4, or 5 container and a ran...
static void computeEdgeMeasure(ArrayOut &outVals, const ArrayJac &inJac, const ArrayWeights &inWeights, const int whichEdge, const shards::CellTopology &parentCell)
Returns the weighted integration measures outVals with dimensions (C,P) used for the computation of e...
static void dataIntegral(ArrayOutData &outputData, const ArrayInDataLeft &inputDataLeft, const ArrayInDataRight &inputDataRight, const ECompEngine compEngine, const bool sumInto=false)
Contracts the point (and space) dimensions P (and D1 and D2) of two rank-2, 3, or 4 containers with d...
static void tensorMultiplyDataField(ArrayOutFields &outputFields, const ArrayInData &inputData, const ArrayInFields &inputFields, const char transpose='N')
Matrix-vector or matrix-matrix product of data and fields; please read the description below.
static void HCURLtransformCURL(ArrayTypeOut &outVals, const ArrayTypeJac &jacobian, const ArrayTypeDet &jacobianDet, const ArrayTypeIn &inVals, const char transpose='N')
Transformation of a curl field in the H-curl space, defined at points on a reference cell,...
static void HDIVtransformVALUE(ArrayTypeOut &outVals, const ArrayTypeJac &jacobian, const ArrayTypeDet &jacobianDet, const ArrayTypeIn &inVals, const char transpose='N')
Transformation of a (vector) value field in the H-div space, defined at points on a reference cell,...
static Scalar vectorNorm(const Scalar *inVec, const size_t dim, const ENorm normType)
Computes norm (1, 2, infinity) of the vector inVec of size dim.