Intrepid
http://trilinos.sandia.gov/packages/docs/r10.10/packages/intrepid/src/Shared/Intrepid_ArrayToolsDefScalar.hpp
Go to the documentation of this file.
00001 // @HEADER
00002 // ************************************************************************
00003 //
00004 //                           Intrepid Package
00005 //                 Copyright (2007) Sandia Corporation
00006 //
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 //
00010 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Pavel Bochev  (pbboche@sandia.gov)
00038 //                    Denis Ridzal  (dridzal@sandia.gov), or
00039 //                    Kara Peterson (kjpeter@sandia.gov)
00040 //
00041 // ************************************************************************
00042 // @HEADER
00043 
00049 namespace Intrepid {
00050 
00051 template<class Scalar, class ArrayOutFields, class ArrayInData, class ArrayInFields>
00052 void ArrayTools::scalarMultiplyDataField(ArrayOutFields &     outputFields,
00053                                          const ArrayInData &  inputData,
00054                                          ArrayInFields &      inputFields,
00055                                          const bool           reciprocal) {
00056 
00057 #ifdef HAVE_INTREPID_DEBUG
00058   TEUCHOS_TEST_FOR_EXCEPTION( (inputData.rank() != 2), std::invalid_argument,
00059                       ">>> ERROR (ArrayTools::scalarMultiplyDataField): Input data container must have rank 2.");
00060   if (outputFields.rank() <= inputFields.rank()) {
00061     TEUCHOS_TEST_FOR_EXCEPTION( ( (inputFields.rank() < 3) || (inputFields.rank() > 5) ), std::invalid_argument,
00062                         ">>> ERROR (ArrayTools::scalarMultiplyDataField): Input fields container must have rank 3, 4, or 5.");
00063     TEUCHOS_TEST_FOR_EXCEPTION( (outputFields.rank() != inputFields.rank()), std::invalid_argument,
00064                         ">>> ERROR (ArrayTools::scalarMultiplyDataField): Input and output fields containers must have the same rank.");
00065     TEUCHOS_TEST_FOR_EXCEPTION( (inputFields.dimension(0) != inputData.dimension(0) ), std::invalid_argument,
00066                         ">>> ERROR (ArrayTools::scalarMultiplyDataField): Zeroth dimensions (number of integration domains) of the fields and data input containers must agree!");
00067     TEUCHOS_TEST_FOR_EXCEPTION( ( (inputFields.dimension(2) != inputData.dimension(1)) && (inputData.dimension(1) != 1) ), std::invalid_argument,
00068                         ">>> ERROR (ArrayTools::scalarMultiplyDataField): Second dimension of the fields input container and first dimension of data input container (number of integration points) must agree or first data dimension must be 1!");
00069     for (int i=0; i<inputFields.rank(); i++) {
00070       std::string errmsg  = ">>> ERROR (ArrayTools::scalarMultiplyDataField): Dimension ";
00071       errmsg += (char)(48+i);
00072       errmsg += " of the input and output fields containers must agree!";
00073       TEUCHOS_TEST_FOR_EXCEPTION( (inputFields.dimension(i) != outputFields.dimension(i)), std::invalid_argument, errmsg );
00074     }
00075   }
00076   else {
00077     TEUCHOS_TEST_FOR_EXCEPTION( ( (inputFields.rank() < 2) || (inputFields.rank() > 4) ), std::invalid_argument,
00078                         ">>> ERROR (ArrayTools::scalarMultiplyDataField): Input fields container must have rank 2, 3, or 4.");
00079     TEUCHOS_TEST_FOR_EXCEPTION( (outputFields.rank() != inputFields.rank()+1), std::invalid_argument,
00080                         ">>> ERROR (ArrayTools::scalarMultiplyDataField): The rank of the input fields container must be one less than the rank of the output fields container.");
00081     TEUCHOS_TEST_FOR_EXCEPTION( ( (inputFields.dimension(1) != inputData.dimension(1)) && (inputData.dimension(1) != 1) ), std::invalid_argument,
00082                         ">>> ERROR (ArrayTools::scalarMultiplyDataField): First dimensions of fields input container and data input container (number of integration points) must agree or first data dimension must be 1!");
00083     TEUCHOS_TEST_FOR_EXCEPTION( ( inputData.dimension(0) != outputFields.dimension(0) ), std::invalid_argument,
00084                         ">>> ERROR (ArrayTools::scalarMultiplyDataField): Zeroth dimensions of fields output container and data input containers (number of integration domains) must agree!");
00085     for (int i=0; i<inputFields.rank(); i++) {
00086       std::string errmsg  = ">>> ERROR (ArrayTools::scalarMultiplyDataField): Dimensions ";
00087       errmsg += (char)(48+i);
00088       errmsg += " and ";
00089       errmsg += (char)(48+i+1);
00090       errmsg += " of the input and output fields containers must agree!";
00091       TEUCHOS_TEST_FOR_EXCEPTION( (inputFields.dimension(i) != outputFields.dimension(i+1)), std::invalid_argument, errmsg );
00092     }
00093   }
00094 #endif
00095 
00096   // get sizes
00097   int invalRank      = inputFields.rank();
00098   int outvalRank     = outputFields.rank();
00099   int numCells       = outputFields.dimension(0);
00100   int numFields      = outputFields.dimension(1);
00101   int numPoints      = outputFields.dimension(2);
00102   int numDataPoints  = inputData.dimension(1);
00103   int dim1Tens       = 0;
00104   int dim2Tens       = 0;
00105   if (outvalRank > 3) {
00106     dim1Tens = outputFields.dimension(3);
00107     if (outvalRank > 4) {
00108       dim2Tens = outputFields.dimension(4);
00109     }
00110   }
00111 
00112   if (outvalRank == invalRank) {
00113 
00114     if (numDataPoints != 1) { // nonconstant data
00115 
00116       switch(invalRank) {
00117         case 3: {
00118           if (reciprocal) {
00119             for(int cl = 0; cl < numCells; cl++) {
00120               for(int bf = 0; bf < numFields; bf++) {
00121                 for(int pt = 0; pt < numPoints; pt++) {
00122                   outputFields(cl, bf, pt) = inputFields(cl, bf, pt)/inputData(cl, pt);
00123                 } // P-loop
00124               } // F-loop
00125             } // C-loop
00126           }
00127           else {
00128             for(int cl = 0; cl < numCells; cl++) {
00129               for(int bf = 0; bf < numFields; bf++) {
00130                 for(int pt = 0; pt < numPoints; pt++) {
00131                   outputFields(cl, bf, pt) = inputFields(cl, bf, pt)*inputData(cl, pt);
00132                 } // P-loop
00133               } // F-loop
00134             } // C-loop
00135           }
00136         }// case 3
00137         break;
00138 
00139         case 4: {
00140           if (reciprocal) {
00141             for(int cl = 0; cl < numCells; cl++) {
00142               for(int bf = 0; bf < numFields; bf++) {
00143                 for(int pt = 0; pt < numPoints; pt++) {
00144                   for( int iVec = 0; iVec < dim1Tens; iVec++) {
00145                     outputFields(cl, bf, pt, iVec) = inputFields(cl, bf, pt, iVec)/inputData(cl, pt);
00146                   } // D1-loop
00147                 } // P-loop
00148               } // F-loop
00149             } // C-loop
00150           }
00151           else {
00152             for(int cl = 0; cl < numCells; cl++) {
00153               for(int bf = 0; bf < numFields; bf++) {
00154                 for(int pt = 0; pt < numPoints; pt++) {
00155                   for( int iVec = 0; iVec < dim1Tens; iVec++) {
00156                     outputFields(cl, bf, pt, iVec) = inputFields(cl, bf, pt, iVec)*inputData(cl, pt);
00157                   } // D1-loop
00158                 } // P-loop
00159               } // F-loop
00160             } // C-loop
00161           }
00162         }// case 4
00163         break;
00164 
00165         case 5: {
00166           if (reciprocal) {
00167             for(int cl = 0; cl < numCells; cl++) {
00168               for(int bf = 0; bf < numFields; bf++) {
00169                 for(int pt = 0; pt < numPoints; pt++) {
00170                   for( int iTens1 = 0; iTens1 < dim1Tens; iTens1++) {
00171                     for( int iTens2 = 0; iTens2 < dim2Tens; iTens2++) {
00172                       outputFields(cl, bf, pt, iTens1, iTens2) = inputFields(cl, bf, pt, iTens1, iTens2)/inputData(cl, pt);
00173                     } // D2-loop
00174                   } // D1-loop
00175                 } // F-loop
00176               } // P-loop
00177             } // C-loop
00178           }
00179           else {
00180             for(int cl = 0; cl < numCells; cl++) {
00181               for(int bf = 0; bf < numFields; bf++) {
00182                 for(int pt = 0; pt < numPoints; pt++) {
00183                   for( int iTens1 = 0; iTens1 < dim1Tens; iTens1++) {
00184                     for( int iTens2 = 0; iTens2 < dim2Tens; iTens2++) {
00185                       outputFields(cl, bf, pt, iTens1, iTens2) = inputFields(cl, bf, pt, iTens1, iTens2)*inputData(cl, pt);
00186                     } // D2-loop
00187                   } // D1-loop
00188                 } // F-loop
00189               } // P-loop
00190             } // C-loop
00191           }
00192         }// case 5
00193         break;
00194 
00195         default:
00196               TEUCHOS_TEST_FOR_EXCEPTION( !( (invalRank == 3) || (invalRank == 4) || (invalRank == 5) ), std::invalid_argument,
00197                                   ">>> ERROR (ArrayTools::scalarMultiplyDataField): This branch of the method is defined only for rank-3,4 or 5 containers.");
00198       }// invalRank
00199 
00200     }
00201     else { //constant data
00202 
00203       switch(invalRank) {
00204         case 3: {
00205           if (reciprocal) {
00206             for(int cl = 0; cl < numCells; cl++) {
00207               for(int bf = 0; bf < numFields; bf++) {
00208                 for(int pt = 0; pt < numPoints; pt++) {
00209                   outputFields(cl, bf, pt) = inputFields(cl, bf, pt)/inputData(cl, 0);
00210                 } // P-loop
00211               } // F-loop
00212             } // C-loop
00213           }
00214           else {
00215             for(int cl = 0; cl < numCells; cl++) {
00216               for(int bf = 0; bf < numFields; bf++) {
00217                 for(int pt = 0; pt < numPoints; pt++) {
00218                   outputFields(cl, bf, pt) = inputFields(cl, bf, pt)*inputData(cl, 0);
00219                 } // P-loop
00220               } // F-loop
00221             } // C-loop
00222           }
00223         }// case 3
00224         break;
00225 
00226         case 4: {
00227           if (reciprocal) {
00228             for(int cl = 0; cl < numCells; cl++) {
00229               for(int bf = 0; bf < numFields; bf++) {
00230                 for(int pt = 0; pt < numPoints; pt++) {
00231                   for( int iVec = 0; iVec < dim1Tens; iVec++) {
00232                     outputFields(cl, bf, pt, iVec) = inputFields(cl, bf, pt, iVec)/inputData(cl, 0);
00233                   } // D1-loop
00234                 } // P-loop
00235               } // F-loop
00236             } // C-loop
00237           }
00238           else {
00239             for(int cl = 0; cl < numCells; cl++) {
00240               for(int bf = 0; bf < numFields; bf++) {
00241                 for(int pt = 0; pt < numPoints; pt++) {
00242                   for( int iVec = 0; iVec < dim1Tens; iVec++) {
00243                     outputFields(cl, bf, pt, iVec) = inputFields(cl, bf, pt, iVec)*inputData(cl, 0);
00244                   } // D1-loop
00245                 } // P-loop
00246               } // F-loop
00247             } // C-loop
00248           }
00249         }// case 4
00250         break;
00251 
00252         case 5: {
00253           if (reciprocal) {
00254             for(int cl = 0; cl < numCells; cl++) {
00255               for(int bf = 0; bf < numFields; bf++) {
00256                 for(int pt = 0; pt < numPoints; pt++) {
00257                   for( int iTens1 = 0; iTens1 < dim1Tens; iTens1++) {
00258                     for( int iTens2 = 0; iTens2 < dim2Tens; iTens2++) {
00259                       outputFields(cl, bf, pt, iTens1, iTens2) = inputFields(cl, bf, pt, iTens1, iTens2)/inputData(cl, 0);
00260                     } // D2-loop
00261                   } // D1-loop
00262                 } // F-loop
00263               } // P-loop
00264             } // C-loop
00265           }
00266           else {
00267             for(int cl = 0; cl < numCells; cl++) {
00268               for(int bf = 0; bf < numFields; bf++) {
00269                 for(int pt = 0; pt < numPoints; pt++) {
00270                   for( int iTens1 = 0; iTens1 < dim1Tens; iTens1++) {
00271                     for( int iTens2 = 0; iTens2 < dim2Tens; iTens2++) {
00272                       outputFields(cl, bf, pt, iTens1, iTens2) = inputFields(cl, bf, pt, iTens1, iTens2)*inputData(cl, 0);
00273                     } // D2-loop
00274                   } // D1-loop
00275                 } // F-loop
00276               } // P-loop
00277             } // C-loop
00278           }
00279         }// case 5
00280         break;
00281 
00282         default:
00283               TEUCHOS_TEST_FOR_EXCEPTION( !( (invalRank == 3) || (invalRank == 4) || (invalRank == 5) ), std::invalid_argument,
00284                                   ">>> ERROR (ArrayTools::scalarMultiplyDataField): This branch of the method is defined only for rank-3, 4 or 5 input containers.");
00285 
00286       } // invalRank
00287     } // numDataPoints
00288 
00289   }
00290   else {
00291 
00292     if (numDataPoints != 1) { // nonconstant data
00293 
00294       switch(invalRank) {
00295         case 2: {
00296           if (reciprocal) {
00297             for(int cl = 0; cl < numCells; cl++) {
00298               for(int bf = 0; bf < numFields; bf++) {
00299                 for(int pt = 0; pt < numPoints; pt++) {
00300                   outputFields(cl, bf, pt) = inputFields(bf, pt)/inputData(cl, pt);
00301                 } // P-loop
00302               } // F-loop
00303             } // C-loop
00304           }
00305           else {
00306             for(int cl = 0; cl < numCells; cl++) {
00307               for(int bf = 0; bf < numFields; bf++) {
00308                 for(int pt = 0; pt < numPoints; pt++) {
00309                   outputFields(cl, bf, pt) = inputFields(bf, pt)*inputData(cl, pt);
00310                 } // P-loop
00311               } // F-loop
00312             } // C-loop
00313           }
00314         }// case 2
00315         break;
00316 
00317         case 3: {
00318           if (reciprocal) {
00319             for(int cl = 0; cl < numCells; cl++) {
00320               for(int bf = 0; bf < numFields; bf++) {
00321                 for(int pt = 0; pt < numPoints; pt++) {
00322                   for( int iVec = 0; iVec < dim1Tens; iVec++) {
00323                     outputFields(cl, bf, pt, iVec) = inputFields(bf, pt, iVec)/inputData(cl, pt);
00324                   } // D1-loop
00325                 } // P-loop
00326               } // F-loop
00327             } // C-loop
00328           }
00329           else {
00330             for(int cl = 0; cl < numCells; cl++) {
00331               for(int bf = 0; bf < numFields; bf++) {
00332                 for(int pt = 0; pt < numPoints; pt++) {
00333                   for( int iVec = 0; iVec < dim1Tens; iVec++) {
00334                     outputFields(cl, bf, pt, iVec) = inputFields(bf, pt, iVec)*inputData(cl, pt);
00335                   } // D1-loop
00336                 } // P-loop
00337               } // F-loop
00338             } // C-loop
00339           }
00340         }// case 3
00341         break;
00342 
00343         case 4: {
00344           if (reciprocal) {
00345             for(int cl = 0; cl < numCells; cl++) {
00346               for(int bf = 0; bf < numFields; bf++) {
00347                 for(int pt = 0; pt < numPoints; pt++) {
00348                   for( int iTens1 = 0; iTens1 < dim1Tens; iTens1++) {
00349                     for( int iTens2 = 0; iTens2 < dim2Tens; iTens2++) {
00350                       outputFields(cl, bf, pt, iTens1, iTens2) = inputFields(bf, pt, iTens1, iTens2)/inputData(cl, pt);
00351                     } // D2-loop
00352                   } // D1-loop
00353                 } // F-loop
00354               } // P-loop
00355             } // C-loop
00356           }
00357           else {
00358             for(int cl = 0; cl < numCells; cl++) {
00359               for(int bf = 0; bf < numFields; bf++) {
00360                 for(int pt = 0; pt < numPoints; pt++) {
00361                   for( int iTens1 = 0; iTens1 < dim1Tens; iTens1++) {
00362                     for( int iTens2 = 0; iTens2 < dim2Tens; iTens2++) {
00363                       outputFields(cl, bf, pt, iTens1, iTens2) = inputFields(bf, pt, iTens1, iTens2)*inputData(cl, pt);
00364                     } // D2-loop
00365                   } // D1-loop
00366                 } // F-loop
00367               } // P-loop
00368             } // C-loop
00369           }
00370         }// case 4
00371         break;
00372 
00373         default:
00374               TEUCHOS_TEST_FOR_EXCEPTION( !( (invalRank == 2) || (invalRank == 3) || (invalRank == 4) ), std::invalid_argument,
00375                                   ">>> ERROR (ArrayTools::scalarMultiplyDataField): This branch of the method is defined only for rank-2, 3 or 4 input containers.");
00376       }// invalRank
00377 
00378     }
00379     else { //constant data
00380 
00381       switch(invalRank) {
00382         case 2: {
00383           if (reciprocal) {
00384             for(int cl = 0; cl < numCells; cl++) {
00385               for(int bf = 0; bf < numFields; bf++) {
00386                 for(int pt = 0; pt < numPoints; pt++) {
00387                   outputFields(cl, bf, pt) = inputFields(bf, pt)/inputData(cl, 0);
00388                 } // P-loop
00389               } // F-loop
00390             } // C-loop
00391           }
00392           else {
00393             for(int cl = 0; cl < numCells; cl++) {
00394               for(int bf = 0; bf < numFields; bf++) {
00395                 for(int pt = 0; pt < numPoints; pt++) {
00396                   outputFields(cl, bf, pt) = inputFields(bf, pt)*inputData(cl, 0);
00397                 } // P-loop
00398               } // F-loop
00399             } // C-loop
00400           }
00401         }// case 2
00402         break;
00403 
00404         case 3: {
00405           if (reciprocal) {
00406             for(int cl = 0; cl < numCells; cl++) {
00407               for(int bf = 0; bf < numFields; bf++) {
00408                 for(int pt = 0; pt < numPoints; pt++) {
00409                   for( int iVec = 0; iVec < dim1Tens; iVec++) {
00410                     outputFields(cl, bf, pt, iVec) = inputFields(bf, pt, iVec)/inputData(cl, 0);
00411                   } // D1-loop
00412                 } // P-loop
00413               } // F-loop
00414             } // C-loop
00415           }
00416           else {
00417             for(int cl = 0; cl < numCells; cl++) {
00418               for(int bf = 0; bf < numFields; bf++) {
00419                 for(int pt = 0; pt < numPoints; pt++) {
00420                   for( int iVec = 0; iVec < dim1Tens; iVec++) {
00421                     outputFields(cl, bf, pt, iVec) = inputFields(bf, pt, iVec)*inputData(cl, 0);
00422                   } // D1-loop
00423                 } // P-loop
00424               } // F-loop
00425             } // C-loop
00426           }
00427         }// case 3
00428         break;
00429 
00430         case 4: {
00431           if (reciprocal) {
00432             for(int cl = 0; cl < numCells; cl++) {
00433               for(int bf = 0; bf < numFields; bf++) {
00434                 for(int pt = 0; pt < numPoints; pt++) {
00435                   for( int iTens1 = 0; iTens1 < dim1Tens; iTens1++) {
00436                     for( int iTens2 = 0; iTens2 < dim2Tens; iTens2++) {
00437                       outputFields(cl, bf, pt, iTens1, iTens2) = inputFields(bf, pt, iTens1, iTens2)/inputData(cl, 0);
00438                     } // D2-loop
00439                   } // D1-loop
00440                 } // F-loop
00441               } // P-loop
00442             } // C-loop
00443           }
00444           else {
00445             for(int cl = 0; cl < numCells; cl++) {
00446               for(int bf = 0; bf < numFields; bf++) {
00447                 for(int pt = 0; pt < numPoints; pt++) {
00448                   for( int iTens1 = 0; iTens1 < dim1Tens; iTens1++) {
00449                     for( int iTens2 = 0; iTens2 < dim2Tens; iTens2++) {
00450                       outputFields(cl, bf, pt, iTens1, iTens2) = inputFields(bf, pt, iTens1, iTens2)*inputData(cl, 0);
00451                     } // D2-loop
00452                   } // D1-loop
00453                 } // F-loop
00454               } // P-loop
00455             } // C-loop
00456           }
00457         }// case 4
00458         break;
00459 
00460         default:
00461               TEUCHOS_TEST_FOR_EXCEPTION( !( (invalRank == 2) || (invalRank == 3) || (invalRank == 3) ), std::invalid_argument,
00462                                   ">>> ERROR (ArrayTools::scalarMultiplyDataField): This branch of the method is defined only for rank-2, 3 or 4 input containers.");
00463 
00464       } // invalRank
00465     } // numDataPoints
00466 
00467   } // end if (outvalRank = invalRank)
00468 
00469 } // scalarMultiplyDataField
00470 
00471 
00472 
00473 template<class Scalar, class ArrayOutData, class ArrayInDataLeft, class ArrayInDataRight>
00474 void ArrayTools::scalarMultiplyDataData(ArrayOutData &           outputData,
00475                                         ArrayInDataLeft &        inputDataLeft,
00476                                         ArrayInDataRight &       inputDataRight,
00477                                         const bool               reciprocal) {
00478 
00479 #ifdef HAVE_INTREPID_DEBUG
00480   TEUCHOS_TEST_FOR_EXCEPTION( (inputDataLeft.rank() != 2), std::invalid_argument,
00481                       ">>> ERROR (ArrayTools::scalarMultiplyDataData): Left input data container must have rank 2.");
00482   if (outputData.rank() <= inputDataRight.rank()) {
00483     TEUCHOS_TEST_FOR_EXCEPTION( ( (inputDataRight.rank() < 2) || (inputDataRight.rank() > 4) ), std::invalid_argument,
00484                         ">>> ERROR (ArrayTools::scalarMultiplyDataData): Right input data container must have rank 2, 3, or 4.");
00485     TEUCHOS_TEST_FOR_EXCEPTION( (outputData.rank() != inputDataRight.rank()), std::invalid_argument,
00486                         ">>> ERROR (ArrayTools::scalarMultiplyDataData): Right input and output data containers must have the same rank.");
00487     TEUCHOS_TEST_FOR_EXCEPTION( (inputDataRight.dimension(0) != inputDataLeft.dimension(0) ), std::invalid_argument,
00488                         ">>> ERROR (ArrayTools::scalarMultiplyDataData): Zeroth dimensions (number of integration domains) of the left and right data input containers must agree!");
00489     TEUCHOS_TEST_FOR_EXCEPTION( ( (inputDataRight.dimension(1) != inputDataLeft.dimension(1)) && (inputDataLeft.dimension(1) != 1) ), std::invalid_argument,
00490                         ">>> ERROR (ArrayTools::scalarMultiplyDataData): First dimensions of the left and right data input containers (number of integration points) must agree or first dimension of the left data input container must be 1!");
00491     for (int i=0; i<inputDataRight.rank(); i++) {
00492       std::string errmsg  = ">>> ERROR (ArrayTools::scalarMultiplyDataData): Dimension ";
00493       errmsg += (char)(48+i);
00494       errmsg += " of the right input and output data containers must agree!";
00495       TEUCHOS_TEST_FOR_EXCEPTION( (inputDataRight.dimension(i) != outputData.dimension(i)), std::invalid_argument, errmsg );
00496     }
00497   }
00498   else {
00499     TEUCHOS_TEST_FOR_EXCEPTION( ( (inputDataRight.rank() < 1) || (inputDataRight.rank() > 3) ), std::invalid_argument,
00500                         ">>> ERROR (ArrayTools::scalarMultiplyDataData): Right input data container must have rank 1, 2, or 3.");
00501     TEUCHOS_TEST_FOR_EXCEPTION( (outputData.rank() != inputDataRight.rank()+1), std::invalid_argument,
00502                         ">>> ERROR (ArrayTools::scalarMultiplyDataData): The rank of the right input data container must be one less than the rank of the output data container.");
00503     TEUCHOS_TEST_FOR_EXCEPTION( ( (inputDataRight.dimension(0) != inputDataLeft.dimension(1)) && (inputDataLeft.dimension(1) != 1) ), std::invalid_argument,
00504                         ">>> ERROR (ArrayTools::scalarMultiplyDataData): Zeroth dimension of the right input data container and first dimension of the left data input container (number of integration points) must agree or first dimension of the left data input container must be 1!");
00505     TEUCHOS_TEST_FOR_EXCEPTION( ( inputDataLeft.dimension(0) != outputData.dimension(0) ), std::invalid_argument,
00506                         ">>> ERROR (ArrayTools::scalarMultiplyDataData): Zeroth dimensions of data output and left data input containers (number of integration domains) must agree!");
00507     for (int i=0; i<inputDataRight.rank(); i++) {
00508       std::string errmsg  = ">>> ERROR (ArrayTools::scalarMultiplyDataData): Dimensions ";
00509       errmsg += (char)(48+i);
00510       errmsg += " and ";
00511       errmsg += (char)(48+i+1);
00512       errmsg += " of the right input and output data containers must agree!";
00513       TEUCHOS_TEST_FOR_EXCEPTION( (inputDataRight.dimension(i) != outputData.dimension(i+1)), std::invalid_argument, errmsg );
00514     }
00515   }
00516 #endif
00517 
00518   // get sizes
00519   int invalRank      = inputDataRight.rank();
00520   int outvalRank     = outputData.rank();
00521   int numCells       = outputData.dimension(0);
00522   int numPoints      = outputData.dimension(1);
00523   int numDataPoints  = inputDataLeft.dimension(1);
00524   int dim1Tens       = 0;
00525   int dim2Tens       = 0;
00526   if (outvalRank > 2) {
00527     dim1Tens = outputData.dimension(2);
00528     if (outvalRank > 3) {
00529       dim2Tens = outputData.dimension(3);
00530     }
00531   }
00532 
00533   if (outvalRank == invalRank) {
00534 
00535     if (numDataPoints != 1) { // nonconstant data
00536 
00537       switch(invalRank) {
00538         case 2: {
00539           if (reciprocal) {
00540             for(int cl = 0; cl < numCells; cl++) {
00541               for(int pt = 0; pt < numPoints; pt++) {
00542                   outputData(cl, pt) = inputDataRight(cl, pt)/inputDataLeft(cl, pt);
00543               } // P-loop
00544             } // C-loop
00545           }
00546           else {
00547             for(int cl = 0; cl < numCells; cl++) {
00548               for(int pt = 0; pt < numPoints; pt++) {
00549                   outputData(cl, pt) = inputDataRight(cl, pt)*inputDataLeft(cl, pt);
00550               } // P-loop
00551             } // C-loop
00552           }
00553         }// case 2
00554         break;
00555 
00556         case 3: {
00557           if (reciprocal) {
00558             for(int cl = 0; cl < numCells; cl++) {
00559               for(int pt = 0; pt < numPoints; pt++) {
00560                 for( int iVec = 0; iVec < dim1Tens; iVec++) {
00561                     outputData(cl, pt, iVec) = inputDataRight(cl, pt, iVec)/inputDataLeft(cl, pt);
00562                 } // D1-loop
00563               } // P-loop
00564             } // C-loop
00565           }
00566           else {
00567             for(int cl = 0; cl < numCells; cl++) {
00568               for(int pt = 0; pt < numPoints; pt++) {
00569                 for( int iVec = 0; iVec < dim1Tens; iVec++) {
00570                     outputData(cl, pt, iVec) = inputDataRight(cl, pt, iVec)*inputDataLeft(cl, pt);
00571                 } // D1-loop
00572               } // P-loop
00573             } // C-loop
00574           }
00575         }// case 3
00576         break;
00577 
00578         case 4: {
00579           if (reciprocal) {
00580             for(int cl = 0; cl < numCells; cl++) {
00581               for(int pt = 0; pt < numPoints; pt++) {
00582                 for( int iTens1 = 0; iTens1 < dim1Tens; iTens1++) {
00583                   for( int iTens2 = 0; iTens2 < dim2Tens; iTens2++) {
00584                       outputData(cl, pt, iTens1, iTens2) = inputDataRight(cl, pt, iTens1, iTens2)/inputDataLeft(cl, pt);
00585                   } // D2-loop
00586                 } // D1-loop
00587               } // P-loop
00588             } // C-loop
00589           }
00590           else {
00591             for(int cl = 0; cl < numCells; cl++) {
00592               for(int pt = 0; pt < numPoints; pt++) {
00593                 for( int iTens1 = 0; iTens1 < dim1Tens; iTens1++) {
00594                   for( int iTens2 = 0; iTens2 < dim2Tens; iTens2++) {
00595                       outputData(cl, pt, iTens1, iTens2) = inputDataRight(cl, pt, iTens1, iTens2)*inputDataLeft(cl, pt);
00596                   } // D2-loop
00597                 } // D1-loop
00598               } // P-loop
00599             } // C-loop
00600           }
00601         }// case 4
00602         break;
00603 
00604         default:
00605               TEUCHOS_TEST_FOR_EXCEPTION( !( (invalRank == 2) || (invalRank == 3) || (invalRank == 4) ), std::invalid_argument,
00606                                   ">>> ERROR (ArrayTools::scalarMultiplyDataData): This branch of the method is defined only for rank-2, 3 or 4 containers.");
00607       }// invalRank
00608 
00609     }
00610     else { // constant left data
00611 
00612       switch(invalRank) {
00613         case 2: {
00614           if (reciprocal) {
00615             for(int cl = 0; cl < numCells; cl++) {
00616               for(int pt = 0; pt < numPoints; pt++) {
00617                   outputData(cl, pt) = inputDataRight(cl, pt)/inputDataLeft(cl, 0);
00618               } // P-loop
00619             } // C-loop
00620           }
00621           else {
00622             for(int cl = 0; cl < numCells; cl++) {
00623               for(int pt = 0; pt < numPoints; pt++) {
00624                   outputData(cl, pt) = inputDataRight(cl, pt)*inputDataLeft(cl, 0);
00625               } // P-loop
00626             } // C-loop
00627           }
00628         }// case 2
00629         break;
00630 
00631         case 3: {
00632           if (reciprocal) {
00633             for(int cl = 0; cl < numCells; cl++) {
00634               for(int pt = 0; pt < numPoints; pt++) {
00635                 for( int iVec = 0; iVec < dim1Tens; iVec++) {
00636                     outputData(cl, pt, iVec) = inputDataRight(cl, pt, iVec)/inputDataLeft(cl, 0);
00637                 } // D1-loop
00638               } // P-loop
00639             } // C-loop
00640           }
00641           else {
00642             for(int cl = 0; cl < numCells; cl++) {
00643               for(int pt = 0; pt < numPoints; pt++) {
00644                 for( int iVec = 0; iVec < dim1Tens; iVec++) {
00645                     outputData(cl, pt, iVec) = inputDataRight(cl, pt, iVec)*inputDataLeft(cl, 0);
00646                 } // D1-loop
00647               } // P-loop
00648             } // C-loop
00649           }
00650         }// case 3
00651         break;
00652 
00653         case 4: {
00654           if (reciprocal) {
00655             for(int cl = 0; cl < numCells; cl++) {
00656               for(int pt = 0; pt < numPoints; pt++) {
00657                 for( int iTens1 = 0; iTens1 < dim1Tens; iTens1++) {
00658                   for( int iTens2 = 0; iTens2 < dim2Tens; iTens2++) {
00659                       outputData(cl, pt, iTens1, iTens2) = inputDataRight(cl, pt, iTens1, iTens2)/inputDataLeft(cl, 0);
00660                   } // D2-loop
00661                 } // D1-loop
00662               } // P-loop
00663             } // C-loop
00664           }
00665           else {
00666             for(int cl = 0; cl < numCells; cl++) {
00667               for(int pt = 0; pt < numPoints; pt++) {
00668                 for( int iTens1 = 0; iTens1 < dim1Tens; iTens1++) {
00669                   for( int iTens2 = 0; iTens2 < dim2Tens; iTens2++) {
00670                       outputData(cl, pt, iTens1, iTens2) = inputDataRight(cl, pt, iTens1, iTens2)*inputDataLeft(cl, 0);
00671                   } // D2-loop
00672                 } // D1-loop
00673               } // P-loop
00674             } // C-loop
00675           }
00676         }// case 4
00677         break;
00678 
00679         default:
00680               TEUCHOS_TEST_FOR_EXCEPTION( !( (invalRank == 2) || (invalRank == 3) || (invalRank == 4) ), std::invalid_argument,
00681                                   ">>> ERROR (ArrayTools::scalarMultiplyDataData): This branch of the method is defined only for rank-2, 3 or 4 input containers.");
00682 
00683       } // invalRank
00684     } // numDataPoints
00685 
00686   }
00687   else {
00688 
00689     if (numDataPoints != 1) { // nonconstant data
00690 
00691       switch(invalRank) {
00692         case 1: {
00693           if (reciprocal) {
00694             for(int cl = 0; cl < numCells; cl++) {
00695               for(int pt = 0; pt < numPoints; pt++) {
00696                   outputData(cl, pt) = inputDataRight(pt)/inputDataLeft(cl, pt);
00697               } // P-loop
00698             } // C-loop
00699           }
00700           else {
00701             for(int cl = 0; cl < numCells; cl++) {
00702               for(int pt = 0; pt < numPoints; pt++) {
00703                   outputData(cl, pt) = inputDataRight(pt)*inputDataLeft(cl, pt);
00704               } // P-loop
00705             } // C-loop
00706           }
00707         }// case 1
00708         break;
00709 
00710         case 2: {
00711           if (reciprocal) {
00712             for(int cl = 0; cl < numCells; cl++) {
00713               for(int pt = 0; pt < numPoints; pt++) {
00714                 for( int iVec = 0; iVec < dim1Tens; iVec++) {
00715                     outputData(cl, pt, iVec) = inputDataRight(pt, iVec)/inputDataLeft(cl, pt);
00716                 } // D1-loop
00717               } // P-loop
00718             } // C-loop
00719           }
00720           else {
00721             for(int cl = 0; cl < numCells; cl++) {
00722               for(int pt = 0; pt < numPoints; pt++) {
00723                 for( int iVec = 0; iVec < dim1Tens; iVec++) {
00724                     outputData(cl, pt, iVec) = inputDataRight(pt, iVec)*inputDataLeft(cl, pt);
00725                 } // D1-loop
00726               } // P-loop
00727             } // C-loop
00728           }
00729         }// case 2
00730         break;
00731 
00732         case 3: {
00733           if (reciprocal) {
00734             for(int cl = 0; cl < numCells; cl++) {
00735               for(int pt = 0; pt < numPoints; pt++) {
00736                 for( int iTens1 = 0; iTens1 < dim1Tens; iTens1++) {
00737                   for( int iTens2 = 0; iTens2 < dim2Tens; iTens2++) {
00738                       outputData(cl, pt, iTens1, iTens2) = inputDataRight(pt, iTens1, iTens2)/inputDataLeft(cl, pt);
00739                   } // D2-loop
00740                 } // D1-loop
00741               } // P-loop
00742             } // C-loop
00743           }
00744           else {
00745             for(int cl = 0; cl < numCells; cl++) {
00746               for(int pt = 0; pt < numPoints; pt++) {
00747                 for( int iTens1 = 0; iTens1 < dim1Tens; iTens1++) {
00748                   for( int iTens2 = 0; iTens2 < dim2Tens; iTens2++) {
00749                       outputData(cl, pt, iTens1, iTens2) = inputDataRight(pt, iTens1, iTens2)*inputDataLeft(cl, pt);
00750                   } // D2-loop
00751                 } // D1-loop
00752               } // P-loop
00753             } // C-loop
00754           }
00755         }// case 3
00756         break;
00757 
00758         default:
00759               TEUCHOS_TEST_FOR_EXCEPTION( !( (invalRank == 1) || (invalRank == 2) || (invalRank == 3) ), std::invalid_argument,
00760                                   ">>> ERROR (ArrayTools::scalarMultiplyDataData): This branch of the method is defined only for rank-1, 2 or 3 input containers.");
00761       }// invalRank
00762 
00763     }
00764     else { //constant data
00765 
00766       switch(invalRank) {
00767         case 1: {
00768           if (reciprocal) {
00769             for(int cl = 0; cl < numCells; cl++) {
00770               for(int pt = 0; pt < numPoints; pt++) {
00771                   outputData(cl, pt) = inputDataRight(pt)/inputDataLeft(cl, 0);
00772               } // P-loop
00773             } // C-loop
00774           }
00775           else {
00776             for(int cl = 0; cl < numCells; cl++) {
00777               for(int pt = 0; pt < numPoints; pt++) {
00778                   outputData(cl, pt) = inputDataRight(pt)*inputDataLeft(cl, 0);
00779               } // P-loop
00780             } // C-loop
00781           }
00782         }// case 1
00783         break;
00784 
00785         case 2: {
00786           if (reciprocal) {
00787             for(int cl = 0; cl < numCells; cl++) {
00788               for(int pt = 0; pt < numPoints; pt++) {
00789                   for( int iVec = 0; iVec < dim1Tens; iVec++) {
00790                     outputData(cl, pt, iVec) = inputDataRight(pt, iVec)/inputDataLeft(cl, 0);
00791                 } // D1-loop
00792               } // P-loop
00793             } // C-loop
00794           }
00795           else {
00796             for(int cl = 0; cl < numCells; cl++) {
00797               for(int pt = 0; pt < numPoints; pt++) {
00798                 for( int iVec = 0; iVec < dim1Tens; iVec++) {
00799                     outputData(cl, pt, iVec) = inputDataRight(pt, iVec)*inputDataLeft(cl, 0);
00800                 } // D1-loop
00801               } // P-loop
00802             } // C-loop
00803           }
00804         }// case 2
00805         break;
00806 
00807         case 3: {
00808           if (reciprocal) {
00809             for(int cl = 0; cl < numCells; cl++) {
00810               for(int pt = 0; pt < numPoints; pt++) {
00811                 for( int iTens1 = 0; iTens1 < dim1Tens; iTens1++) {
00812                   for( int iTens2 = 0; iTens2 < dim2Tens; iTens2++) {
00813                       outputData(cl, pt, iTens1, iTens2) = inputDataRight(pt, iTens1, iTens2)/inputDataLeft(cl, 0);
00814                   } // D2-loop
00815                 } // D1-loop
00816               } // P-loop
00817             } // C-loop
00818           }
00819           else {
00820             for(int cl = 0; cl < numCells; cl++) {
00821               for(int pt = 0; pt < numPoints; pt++) {
00822                 for( int iTens1 = 0; iTens1 < dim1Tens; iTens1++) {
00823                   for( int iTens2 = 0; iTens2 < dim2Tens; iTens2++) {
00824                       outputData(cl, pt, iTens1, iTens2) = inputDataRight(pt, iTens1, iTens2)*inputDataLeft(cl, 0);
00825                   } // D2-loop
00826                 } // D1-loop
00827               } // P-loop
00828             } // C-loop
00829           }
00830         }// case 3
00831         break;
00832 
00833         default:
00834               TEUCHOS_TEST_FOR_EXCEPTION( !( (invalRank == 1) || (invalRank == 2) || (invalRank == 3) ), std::invalid_argument,
00835                                   ">>> ERROR (ArrayTools::scalarMultiplyDataData): This branch of the method is defined only for rank-1, 2 or 3 input containers.");
00836 
00837       } // invalRank
00838     } // numDataPoints
00839 
00840   } // end if (outvalRank = invalRank)
00841 
00842 } // scalarMultiplyDataData
00843 
00844 } // end namespace Intrepid