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 
50 namespace Intrepid {
51 
52 template<class Scalar, class ArrayTypeOut, class ArrayTypeIn>
53 void FunctionSpaceTools::HGRADtransformVALUE(ArrayTypeOut & outVals,
54  const ArrayTypeIn & inVals) {
55 
56  ArrayTools::cloneFields<Scalar>(outVals, inVals);
57 
58 }
59 
60 template<class Scalar, class ArrayTypeOut, class ArrayTypeJac, class ArrayTypeIn>
61 void 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 
70 template<class Scalar, class ArrayTypeOut, class ArrayTypeJac, class ArrayTypeIn>
71 void 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 
80 template<class Scalar, class ArrayTypeOut, class ArrayTypeJac, class ArrayTypeDet, class ArrayTypeIn>
81 void 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 
92 template<class Scalar, class ArrayTypeOut, class ArrayTypeJac, class ArrayTypeDet, class ArrayTypeIn>
93 void 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 
104 template<class Scalar, class ArrayTypeOut, class ArrayTypeDet, class ArrayTypeIn>
105 void FunctionSpaceTools::HDIVtransformDIV(ArrayTypeOut & outVals,
106  const ArrayTypeDet & jacobianDet,
107  const ArrayTypeIn & inVals) {
108 
109  ArrayTools::scalarMultiplyDataField<Scalar>(outVals, jacobianDet, inVals, true);
110 
111 }
112 
113 template<class Scalar, class ArrayTypeOut, class ArrayTypeDet, class ArrayTypeIn>
114 void FunctionSpaceTools::HVOLtransformVALUE(ArrayTypeOut & outVals,
115  const ArrayTypeDet & jacobianDet,
116  const ArrayTypeIn & inVals) {
117 
118  ArrayTools::scalarMultiplyDataField<Scalar>(outVals, jacobianDet, inVals, true);
119 
120 }
121 template<class Scalar>
122 void FunctionSpaceTools::integrate(Intrepid::FieldContainer<Scalar> & outputValues,
123  const Intrepid::FieldContainer<Scalar> & leftValues,
124  const Intrepid::FieldContainer<Scalar> & rightValues,
125  const ECompEngine compEngine,
126  const bool sumInto) {
127 int outRank = getrank(outputValues);
128 int lRank = getrank(leftValues);
129 switch (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 
816 TEUCHOS_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;
1204 default:
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
1211 break;
1212 }
1213 }
1214 template<class Scalar, class ArrayOut, class ArrayInLeft, class ArrayInRight>
1215 void 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
1245 template<class Scalar, class ArrayOutFields, class ArrayInFieldsLeft, class ArrayInFieldsRight>
1246 void 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 
1276 template<class Scalar, class ArrayOutFields, class ArrayInData, class ArrayInFields>
1277 void 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 
1297 TEUCHOS_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 
1306 template<class Scalar, class ArrayOutData, class ArrayInDataLeft, class ArrayInDataRight>
1307 void 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 
1336 template<class Scalar, class ArrayOut, class ArrayDet, class ArrayWeights>
1337 inline 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 
1359 template<class Scalar, class ArrayOut, class ArrayJac, class ArrayWeights>
1360 void FunctionSpaceTools::computeFaceMeasure(ArrayOut & outVals,
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 
1388 template<class Scalar, class ArrayOut, class ArrayJac, class ArrayWeights>
1389 void FunctionSpaceTools::computeEdgeMeasure(ArrayOut & outVals,
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 
1418 template<class Scalar, class ArrayTypeOut, class ArrayTypeMeasure, class ArrayTypeIn>
1419 void 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 
1428 template<class Scalar, class ArrayOutFields, class ArrayInData, class ArrayInFields>
1429 void 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 
1439 template<class Scalar, class ArrayOutData, class ArrayInDataLeft, class ArrayInDataRight>
1440 void 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 
1450 template<class Scalar, class ArrayOutFields, class ArrayInData, class ArrayInFields>
1451 void 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 
1460 template<class Scalar, class ArrayOutData, class ArrayInDataLeft, class ArrayInDataRight>
1461 void 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 
1470 template<class Scalar, class ArrayOutFields, class ArrayInData, class ArrayInFields>
1471 void 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 
1493 template<class Scalar, class ArrayOutData, class ArrayInDataLeft, class ArrayInDataRight>
1494 void 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 
1516 template<class Scalar, class ArrayOutFields, class ArrayInData, class ArrayInFields>
1517 void 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 
1538 template<class Scalar, class ArrayOutData, class ArrayInDataLeft, class ArrayInDataRight>
1539 struct 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 };
1556 template<class Scalar, class ArrayOutData, class ArrayInDataLeft, class ArrayInDataRight>
1557 struct 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 };
1566 template<class Scalar, class ArrayOutData, class ArrayInDataLeft, class ArrayInDataRight>
1567 struct 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 };
1576 template<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){
1581  FunctionSpaceTools::tensorMultiplyDataDataTempSpec<Scalar,ArrayOutData,ArrayInDataLeft,ArrayInDataRight,Rank<ArrayOutData>::value>(outputData,inputDataLeft,inputDataRight,transpose);
1582 
1583  }
1584 template<class Scalar, class ArrayTypeInOut, class ArrayTypeSign>
1585 void 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 
1609 template<class Scalar, class ArrayTypeInOut, class ArrayTypeSign>
1610 void 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 
1635 template<class Scalar, class ArrayTypeInOut, class ArrayTypeSign>
1636 void 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 
1728 template<class Scalar, class ArrayOutPointVals, class ArrayInCoeffs, class ArrayInFields>
1729 void 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
Intrepid::FunctionSpaceTools::computeEdgeMeasure
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...
Definition: Intrepid_FunctionSpaceToolsDef.hpp:1389
Intrepid::FunctionSpaceTools::evaluate
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...
Definition: Intrepid_FunctionSpaceToolsDef.hpp:1729
Intrepid::FunctionSpaceTools::HDIVtransformDIV
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...
Definition: Intrepid_FunctionSpaceToolsDef.hpp:105
Intrepid::FunctionSpaceTools::dotMultiplyDataData
static void dotMultiplyDataData(ArrayOutData &outputData, const ArrayInDataLeft &inputDataLeft, const ArrayInDataRight &inputDataRight)
Dot product of data and data; please read the description below.
Definition: Intrepid_FunctionSpaceToolsDef.hpp:1461
Intrepid::CellTools::getPhysicalFaceNormals
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...
Definition: Intrepid_CellToolsDef.hpp:2388
Intrepid::FunctionSpaceTools::tensorMultiplyDataData
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...
Definition: Intrepid_FunctionSpaceToolsDef.hpp:1577
Intrepid::FunctionSpaceTools::integrate
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...
Definition: Intrepid_FunctionSpaceToolsDef.hpp:122
Intrepid::CellTools::getPhysicalEdgeTangents
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 ...
Definition: Intrepid_CellToolsDef.hpp:2239
Intrepid::FunctionSpaceTools::scalarMultiplyDataData
static void scalarMultiplyDataData(ArrayOutData &outputData, ArrayInDataLeft &inputDataLeft, ArrayInDataRight &inputDataRight, const bool reciprocal=false)
Scalar multiplication of data and data; please read the description below.
Definition: Intrepid_FunctionSpaceToolsDef.hpp:1440
Intrepid::FieldContainer::dimension
int dimension(const int whichDim) const
Returns the specified dimension.
Definition: Intrepid_FieldContainerDef.hpp:658
Intrepid::FunctionSpaceTools::multiplyMeasure
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 ...
Definition: Intrepid_FunctionSpaceToolsDef.hpp:1419
Intrepid::FunctionSpaceTools::HVOLtransformVALUE
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...
Definition: Intrepid_FunctionSpaceToolsDef.hpp:114
Intrepid::RealSpaceTools::vectorNorm
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.
Definition: Intrepid_RealSpaceToolsDef.hpp:139
Intrepid::FunctionSpaceTools::applyRightFieldSigns
static void applyRightFieldSigns(ArrayTypeInOut &inoutOperator, const ArrayTypeSign &fieldSigns)
Applies right (column) signs, stored in the user-provided container fieldSigns and indexed by (C...
Definition: Intrepid_FunctionSpaceToolsDef.hpp:1610
Intrepid::FunctionSpaceTools::functionalIntegral
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...
Definition: Intrepid_FunctionSpaceToolsDef.hpp:1277
Intrepid::FunctionSpaceTools::applyLeftFieldSigns
static void applyLeftFieldSigns(ArrayTypeInOut &inoutOperator, const ArrayTypeSign &fieldSigns)
Applies left (row) signs, stored in the user-provided container fieldSigns and indexed by (C...
Definition: Intrepid_FunctionSpaceToolsDef.hpp:1585
Intrepid::FieldContainer::rank
int rank() const
Return rank of the FieldContainer = number of indices used to tag the multi-indexed value...
Definition: Intrepid_FieldContainerDef.hpp:594
Intrepid::FunctionSpaceTools::HDIVtransformVALUE
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...
Definition: Intrepid_FunctionSpaceToolsDef.hpp:93
Intrepid::FunctionSpaceTools::HCURLtransformVALUE
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...
Definition: Intrepid_FunctionSpaceToolsDef.hpp:71
Intrepid::FunctionSpaceTools::dataIntegral
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...
Definition: Intrepid_FunctionSpaceToolsDef.hpp:1307
Intrepid::FunctionSpaceTools::operatorIntegral
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...
Definition: Intrepid_FunctionSpaceToolsDef.hpp:1246
Intrepid::FunctionSpaceTools::vectorMultiplyDataData
static void vectorMultiplyDataData(ArrayOutData &outputData, const ArrayInDataLeft &inputDataLeft, const ArrayInDataRight &inputDataRight)
Cross or outer product of data and data; please read the description below.
Definition: Intrepid_FunctionSpaceToolsDef.hpp:1494
Intrepid::FunctionSpaceTools::computeCellMeasure
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...
Definition: Intrepid_FunctionSpaceToolsDef.hpp:1337
Intrepid::FunctionSpaceTools::vectorMultiplyDataField
static void vectorMultiplyDataField(ArrayOutFields &outputFields, const ArrayInData &inputData, const ArrayInFields &inputFields)
Cross or outer product of data and fields; please read the description below.
Definition: Intrepid_FunctionSpaceToolsDef.hpp:1471
Intrepid::FunctionSpaceTools::dotMultiplyDataField
static void dotMultiplyDataField(ArrayOutFields &outputFields, const ArrayInData &inputData, const ArrayInFields &inputFields)
Dot product of data and fields; please read the description below.
Definition: Intrepid_FunctionSpaceToolsDef.hpp:1451
Intrepid::FunctionSpaceTools::applyFieldSigns
static void applyFieldSigns(ArrayTypeInOut &inoutFunction, const ArrayTypeSign &fieldSigns)
Applies field signs, stored in the user-provided container fieldSigns and indexed by (C...
Definition: Intrepid_FunctionSpaceToolsDef.hpp:1636
Intrepid::FunctionSpaceTools::tensorMultiplyDataField
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...
Definition: Intrepid_FunctionSpaceToolsDef.hpp:1517
Intrepid::FunctionSpaceTools::HGRADtransformGRAD
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...
Definition: Intrepid_FunctionSpaceToolsDef.hpp:61
Intrepid::FunctionSpaceTools::scalarMultiplyDataField
static void scalarMultiplyDataField(ArrayOutFields &outputFields, ArrayInData &inputData, ArrayInFields &inputFields, const bool reciprocal=false)
Scalar multiplication of data and fields; please read the description below.
Definition: Intrepid_FunctionSpaceToolsDef.hpp:1429
Intrepid::FunctionSpaceTools::HGRADtransformVALUE
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...
Definition: Intrepid_FunctionSpaceToolsDef.hpp:53
Intrepid::FunctionSpaceTools::computeFaceMeasure
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...
Definition: Intrepid_FunctionSpaceToolsDef.hpp:1360
Intrepid::FunctionSpaceTools::HCURLtransformCURL
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...
Definition: Intrepid_FunctionSpaceToolsDef.hpp:81
Generated by   doxygen 1.8.5