Intrepid
http://trilinos.sandia.gov/packages/docs/r10.12/packages/intrepid/src/Shared/Intrepid_ArrayToolsDefTensor.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, 
00052            class ArrayOutFields, 
00053            class ArrayInData, 
00054            class ArrayInFields>
00055   void ArrayTools::crossProductDataField(ArrayOutFields &       outputFields,
00056                                          const ArrayInData &    inputData,
00057                                          const ArrayInFields &  inputFields){
00058     
00059 #ifdef HAVE_INTREPID_DEBUG
00060     std::string errmsg = ">>> ERROR (ArrayTools::crossProductDataField):";
00061     /*
00062      *   Check array rank and spatial dimension range (if applicable)
00063      *      (1) inputData(C,P,D);    
00064      *      (2) inputFields(C,F,P,D) or (F,P,D);   
00065      *      (3) outputFields(C,F,P,D) in 3D, or (C,F,P) in 2D
00066      */
00067     // (1) inputData is (C, P, D) and 2 <= D <= 3 is required  
00068     TEUCHOS_TEST_FOR_EXCEPTION( !requireRankRange(errmsg, inputData, 3, 3), 
00069                         std::invalid_argument, errmsg);
00070     TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, inputData, 2, 2,3), 
00071                         std::invalid_argument, errmsg);
00072     // (2) inputFields is (C, F, P, D) or (F, P, D) and 2 <= (D=dimension(rank - 1)) <= 3 is required. 
00073     TEUCHOS_TEST_FOR_EXCEPTION( !requireRankRange(errmsg, inputFields, 3,4), std::invalid_argument, errmsg);
00074     TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, inputFields, inputFields.rank()-1, 2,3), 
00075                         std::invalid_argument, errmsg);
00076     // (3) outputFields is (C,F,P,D) in 3D and (C,F,P) in 2D => rank = inputData.dimension(2) + 1
00077     TEUCHOS_TEST_FOR_EXCEPTION( !requireRankRange(errmsg, outputFields, inputData.dimension(2)+1, inputData.dimension(2)+1), 
00078                         std::invalid_argument, errmsg); 
00079     /*
00080      *   Dimension cross-checks:
00081      *      (1) inputData    vs. inputFields
00082      *      (2) outputFields vs. inputData
00083      *      (3) outputFields vs. inputFields 
00084      *
00085      *   Cross-check (1):
00086      */
00087     if( inputFields.rank() == 4) {
00088       // inputData(C,P,D) vs. inputFields(C,F,P,D): dimensions C, P, D must match 
00089       TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, inputData, 0,1,2,  inputFields, 0,2,3),
00090                           std::invalid_argument, errmsg);
00091     }
00092     else{
00093       // inputData(C,P,D) vs. inputFields(F,P,D): dimensions P, D must match 
00094       TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, inputData, 1,2,  inputFields, 1,2),
00095                           std::invalid_argument, errmsg);      
00096     }
00097     /* 
00098      *  Cross-check (2): 
00099      */
00100     if(inputData.dimension(2) == 2) {
00101       //  in 2D: outputFields(C,F,P) vs. inputData(C,P,D): dimensions C,P must match
00102       TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, outputFields, 0,2,  inputData, 0,1),
00103                           std::invalid_argument, errmsg);
00104     }
00105     else{
00106       // in 3D: outputFields(C,F,P,D) vs. inputData(C,P,D): dimensions C,P,D must match
00107       TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, outputFields, 0,2,3,  inputData, 0,1,2),
00108                           std::invalid_argument, errmsg);
00109     }
00110     /* 
00111      *  Cross-check (3): 
00112      */
00113     if(inputData.dimension(2) == 2) {
00114       // In 2D:
00115       if(inputFields.rank() == 4){
00116         //  and rank-4 inputFields: outputFields(C,F,P) vs. inputFields(C,F,P,D): dimensions C,F,P must match
00117         TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, outputFields, 0,1,2,  inputFields, 0,1,2),
00118                             std::invalid_argument, errmsg);
00119       }
00120       else{
00121         //  and rank-3 inputFields: outputFields(C,F,P) vs. inputFields(F,P,D): dimensions F,P must match
00122         TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, outputFields, 1,2,  inputFields, 0,1),
00123                             std::invalid_argument, errmsg);
00124       }
00125     }
00126     else{
00127       // In 3D:
00128       if(inputFields.rank() == 4){
00129         //  and rank-4 inputFields: outputFields(C,F,P,D) vs. inputFields(C,F,P,D): all dimensions C,F,P,D must match
00130         TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, outputFields,  inputFields),
00131                             std::invalid_argument, errmsg);
00132       }
00133       else{
00134         // and rank-3 inputFields: outputFields(C,F,P,D) vs. inputFields(F,P,D): dimensions F,P,D must match
00135         TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, outputFields, 1,2,3,  inputFields, 0,1,2),
00136                             std::invalid_argument, errmsg);
00137       }
00138     }
00139 #endif  
00140     // 3D cross product
00141     if(inputData.dimension(2) == 3) {
00142       
00143       // inputFields is (C,F,P,D)
00144       if(inputFields.rank() == 4){
00145         
00146         for(int cell = 0; cell < outputFields.dimension(0); cell++){
00147           for(int field = 0; field < outputFields.dimension(1); field++){
00148             for(int point = 0; point < outputFields.dimension(2); point++){
00149               // 
00150               outputFields(cell, field, point, 0) = \
00151                 inputData(cell, point, 1)*inputFields(cell, field, point, 2) - 
00152                 inputData(cell, point, 2)*inputFields(cell, field, point, 1); 
00153               // 
00154               outputFields(cell, field, point, 1) = \
00155                 inputData(cell, point, 2)*inputFields(cell, field, point, 0) - 
00156                 inputData(cell, point, 0)*inputFields(cell, field, point, 2); 
00157               // 
00158               outputFields(cell, field, point, 2) = \
00159                 inputData(cell, point, 0)*inputFields(cell, field, point, 1) - 
00160                 inputData(cell, point, 1)*inputFields(cell, field, point, 0); 
00161             }// point
00162           }// field
00163         } // cell
00164       }// rank = 4
00165       // inputFields is (F,P,D)
00166       else if(inputFields.rank() == 3){
00167         
00168         for(int cell = 0; cell < outputFields.dimension(0); cell++){
00169           for(int field = 0; field < outputFields.dimension(1); field++){
00170             for(int point = 0; point < outputFields.dimension(2); point++){
00171               // 
00172               outputFields(cell, field, point, 0) = \
00173               inputData(cell, point, 1)*inputFields(field, point, 2) - 
00174               inputData(cell, point, 2)*inputFields(field, point, 1); 
00175               // 
00176               outputFields(cell, field, point, 1) = \
00177                 inputData(cell, point, 2)*inputFields(field, point, 0) - 
00178                 inputData(cell, point, 0)*inputFields(field, point, 2); 
00179               // 
00180               outputFields(cell, field, point, 2) = \
00181                 inputData(cell, point, 0)*inputFields(field, point, 1) - 
00182                 inputData(cell, point, 1)*inputFields(field, point, 0); 
00183             }// point
00184           }// field
00185         } // cell
00186       }// rank = 3
00187       else{
00188         TEUCHOS_TEST_FOR_EXCEPTION( true, std::invalid_argument, 
00189                             ">>> ERROR (ArrayTools::crossProductDataField): inputFields rank 3 or 4 required.")
00190       }    
00191     }
00192     // 2D cross product
00193     else if(inputData.dimension(2) == 2){
00194       
00195       // inputFields is (C,F,P,D)
00196       if(inputFields.rank() == 4){
00197         
00198         for(int cell = 0; cell < outputFields.dimension(0); cell++){
00199           for(int field = 0; field < outputFields.dimension(1); field++){
00200             for(int point = 0; point < outputFields.dimension(2); point++){
00201               outputFields(cell, field, point) = \
00202                 inputData(cell, point, 0)*inputFields(cell, field, point, 1) - 
00203                 inputData(cell, point, 1)*inputFields(cell, field, point, 0); 
00204             }// point
00205           }// field
00206         } // cell
00207       }// rank = 4
00208       // inputFields is (F,P,D)
00209       else if(inputFields.rank() == 3) {
00210         
00211         for(int cell = 0; cell < outputFields.dimension(0); cell++){
00212           for(int field = 0; field < outputFields.dimension(1); field++){
00213             for(int point = 0; point < outputFields.dimension(2); point++){
00214               outputFields(cell, field, point) = \
00215                 inputData(cell, point, 0)*inputFields(field, point, 1) - 
00216                 inputData(cell, point, 1)*inputFields(field, point, 0); 
00217             }// point
00218           }// field
00219         } // cell
00220       }// rank = 3
00221     }
00222     // Error: wrong dimension
00223     else {
00224       TEUCHOS_TEST_FOR_EXCEPTION( true, std::invalid_argument, 
00225                           ">>> ERROR (ArrayTools::crossProductDataField): spatial dimension 2 or 3 required.")
00226     }
00227   }
00228   
00229   
00230   
00231   template<class Scalar, 
00232            class ArrayOutData, 
00233            class ArrayInDataLeft, 
00234            class ArrayInDataRight>
00235   void ArrayTools::crossProductDataData(ArrayOutData &            outputData,
00236                                         const ArrayInDataLeft &   inputDataLeft,
00237                                         const ArrayInDataRight &  inputDataRight){
00238 #ifdef HAVE_INTREPID_DEBUG
00239     std::string errmsg = ">>> ERROR (ArrayTools::crossProductDataData):";
00240     /*
00241      *   Check array rank and spatial dimension range (if applicable)
00242      *      (1) inputDataLeft(C,P,D);    
00243      *      (2) inputDataRight(C,P,D) or (P,D);   
00244      *      (3) outputData(C,P,D) in 3D, or (C,P) in 2D
00245      */
00246     // (1) inputDataLeft is (C, P, D) and 2 <= D <= 3 is required  
00247     TEUCHOS_TEST_FOR_EXCEPTION( !requireRankRange(errmsg, inputDataLeft, 3,3), 
00248                         std::invalid_argument, errmsg);
00249     TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, inputDataLeft, 2, 2,3), 
00250                         std::invalid_argument, errmsg);
00251     // (2) inputDataRight is (C, P, D) or (P, D) and 2 <= (D=dimension(rank - 1)) <= 3 is required. 
00252     TEUCHOS_TEST_FOR_EXCEPTION( !requireRankRange(errmsg, inputDataRight, 2,3), 
00253                         std::invalid_argument, errmsg);
00254     TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, inputDataRight, inputDataRight.rank()-1,  2,3), 
00255                         std::invalid_argument, errmsg);
00256     // (3) outputData is (C,P,D) in 3D and (C,P) in 2D => rank = inputDataLeft.dimension(2)
00257     TEUCHOS_TEST_FOR_EXCEPTION( !requireRankRange(errmsg, outputData, inputDataLeft.dimension(2), inputDataLeft.dimension(2)), 
00258                         std::invalid_argument, errmsg); 
00259     /*
00260      *   Dimension cross-checks:
00261      *      (1) inputDataLeft vs. inputDataRight
00262      *      (2) outputData    vs. inputDataLeft
00263      *      (3) outputData    vs. inputDataRight 
00264      *
00265      *   Cross-check (1):
00266      */
00267     if( inputDataRight.rank() == 3) {
00268       // inputDataLeft(C,P,D) vs. inputDataRight(C,P,D): all dimensions C, P, D must match 
00269       TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, inputDataLeft, inputDataRight),
00270                           std::invalid_argument, errmsg);
00271     }
00272     // inputDataLeft(C, P,D) vs. inputDataRight(P,D): dimensions P, D must match
00273     else{
00274       TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, inputDataLeft, 1,2, inputDataRight, 0,1),
00275                           std::invalid_argument, errmsg);      
00276     }
00277     /* 
00278      *  Cross-check (2): 
00279      */
00280     if(inputDataLeft.dimension(2) == 2){
00281       // in 2D: outputData(C,P) vs. inputDataLeft(C,P,D): dimensions C, P must match
00282       TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, inputDataLeft, 0,1,  outputData, 0,1),
00283                           std::invalid_argument, errmsg);
00284     }
00285     else{
00286       // in 3D: outputData(C,P,D) vs. inputDataLeft(C,P,D): all dimensions C, P, D must match
00287       TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, inputDataLeft, outputData),
00288                           std::invalid_argument, errmsg);
00289     }
00290     /* 
00291      *  Cross-check (3): 
00292      */
00293     if(inputDataLeft.dimension(2) == 2) {
00294       // In 2D:
00295       if(inputDataRight.rank() == 3){
00296         //  and rank-3 inputDataRight: outputData(C,P) vs. inputDataRight(C,P,D): dimensions C,P must match
00297         TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, outputData, 0,1,  inputDataRight, 0,1),
00298                             std::invalid_argument, errmsg);
00299       }
00300       else{
00301         //  and rank-2 inputDataRight: outputData(C,P) vs. inputDataRight(P,D): dimension P must match
00302         TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, outputData, 1,  inputDataRight, 0),
00303                             std::invalid_argument, errmsg);
00304       }
00305     }
00306     else{
00307       // In 3D:
00308       if(inputDataRight.rank() == 3){
00309         //  and rank-3 inputDataRight: outputData(C,P,D) vs. inputDataRight(C,P,D): all dimensions C,P,D must match
00310         TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, outputData,  inputDataRight),
00311                             std::invalid_argument, errmsg);
00312       }
00313       else{
00314         //  and rank-2 inputDataRight: outputData(C,P,D) vs. inputDataRight(P,D): dimensions P, D must match
00315         TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, outputData, 1,2,  inputDataRight, 0,1),
00316                             std::invalid_argument, errmsg);
00317       }
00318     }
00319 #endif  
00320     // 3D cross product
00321     if(inputDataLeft.dimension(2) == 3) {
00322       
00323       // inputDataRight is (C,P,D)
00324       if(inputDataRight.rank() == 3){
00325         
00326         for(int cell = 0; cell < inputDataLeft.dimension(0); cell++){
00327           for(int point = 0; point < inputDataLeft.dimension(1); point++){
00328             // 
00329             outputData(cell, point, 0) = \
00330             inputDataLeft(cell, point, 1)*inputDataRight(cell, point, 2) - 
00331             inputDataLeft(cell, point, 2)*inputDataRight(cell, point, 1); 
00332             // 
00333             outputData(cell, point, 1) = \
00334               inputDataLeft(cell, point, 2)*inputDataRight(cell, point, 0) - 
00335               inputDataLeft(cell, point, 0)*inputDataRight(cell, point, 2); 
00336             // 
00337             outputData(cell, point, 2) = \
00338               inputDataLeft(cell, point, 0)*inputDataRight(cell, point, 1) - 
00339               inputDataLeft(cell, point, 1)*inputDataRight(cell, point, 0); 
00340           }// point
00341         } // cell
00342       }// rank = 3
00343        // inputDataRight is (P,D)
00344       else if(inputDataRight.rank() == 2){
00345         
00346         for(int cell = 0; cell < inputDataLeft.dimension(0); cell++){
00347           for(int point = 0; point < inputDataLeft.dimension(1); point++){
00348             // 
00349             outputData(cell, point, 0) = \
00350             inputDataLeft(cell, point, 1)*inputDataRight(point, 2) - 
00351             inputDataLeft(cell, point, 2)*inputDataRight(point, 1); 
00352             // 
00353             outputData(cell, point, 1) = \
00354               inputDataLeft(cell, point, 2)*inputDataRight(point, 0) - 
00355               inputDataLeft(cell, point, 0)*inputDataRight(point, 2); 
00356             // 
00357             outputData(cell, point, 2) = \
00358               inputDataLeft(cell, point, 0)*inputDataRight(point, 1) - 
00359               inputDataLeft(cell, point, 1)*inputDataRight(point, 0); 
00360           }// point
00361         } // cell
00362       }// rank = 2
00363       else{
00364         TEUCHOS_TEST_FOR_EXCEPTION( true, std::invalid_argument, 
00365                             ">>> ERROR (ArrayTools::crossProductDataData): inputDataRight rank 2 or 3 required.")
00366       }    
00367     }
00368     // 2D cross product
00369     else if(inputDataLeft.dimension(2) == 2){
00370       
00371       // inputDataRight is (C,P,D)
00372       if(inputDataRight.rank() == 3){
00373         
00374         for(int cell = 0; cell < inputDataLeft.dimension(0); cell++){
00375             for(int point = 0; point < inputDataLeft.dimension(1); point++){
00376               outputData(cell, point) = \
00377                 inputDataLeft(cell, point, 0)*inputDataRight(cell, point, 1) - 
00378                 inputDataLeft(cell, point, 1)*inputDataRight(cell, point, 0); 
00379             }// point
00380         } // cell
00381       }// rank = 3
00382        // inputDataRight is (P,D)
00383       else if(inputDataRight.rank() == 2) {
00384         
00385         for(int cell = 0; cell < inputDataLeft.dimension(0); cell++){
00386             for(int point = 0; point < inputDataLeft.dimension(1); point++){
00387               outputData(cell, point) = \
00388                 inputDataLeft(cell, point, 0)*inputDataRight(point, 1) - 
00389                 inputDataLeft(cell, point, 1)*inputDataRight(point, 0); 
00390             }// point
00391         } // cell
00392       }// rank = 2
00393     }
00394     // Error: wrong dimension
00395     else {
00396       TEUCHOS_TEST_FOR_EXCEPTION( true, std::invalid_argument, 
00397                           ">>> ERROR (ArrayTools::crossProductDataData): spatial dimension 2 or 3 required.")
00398     }
00399   }
00400   
00401   
00402   
00403   template<class Scalar, 
00404            class ArrayOutFields, 
00405            class ArrayInData, 
00406            class ArrayInFields>
00407   void ArrayTools::outerProductDataField(ArrayOutFields &       outputFields,
00408                                          const ArrayInData &    inputData,
00409                                          const ArrayInFields &  inputFields){
00410 #ifdef HAVE_INTREPID_DEBUG
00411     std::string errmsg = ">>> ERROR (ArrayTools::outerProductDataField):";
00412     /*
00413      *   Check array rank and spatial dimension range (if applicable)
00414      *      (1) inputData(C,P,D);    
00415      *      (2) inputFields(C,F,P,D) or (F,P,D);   
00416      *      (3) outputFields(C,F,P,D,D)
00417      */
00418      // (1) inputData is (C, P, D) and 2 <= D <= 3 is required  
00419     TEUCHOS_TEST_FOR_EXCEPTION( !requireRankRange(errmsg, inputData,  3,3), 
00420                         std::invalid_argument, errmsg);
00421     TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, inputData, 2,  2,3), 
00422                         std::invalid_argument, errmsg);
00423     // (2) inputFields is (C, F, P, D) or (F, P, D) and 2 <= (D=dimension(rank - 1)) <= 3 is required. 
00424     TEUCHOS_TEST_FOR_EXCEPTION( !requireRankRange(errmsg, inputFields, 3,4), std::invalid_argument, errmsg);
00425     TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, inputFields, inputFields.rank()-1,  2,3), 
00426                         std::invalid_argument, errmsg);
00427     // (3) outputFields is (C,F,P,D,D)
00428     TEUCHOS_TEST_FOR_EXCEPTION( !requireRankRange(errmsg, outputFields,  5,5), std::invalid_argument, errmsg);      
00429     TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, outputFields, 3,  2,3), 
00430                         std::invalid_argument, errmsg);
00431     TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, outputFields, 4,  2,3), 
00432                         std::invalid_argument, errmsg);
00433     /*
00434      *   Dimension cross-checks:
00435      *      (1) inputData    vs. inputFields
00436      *      (2) outputFields vs. inputData
00437      *      (3) outputFields vs. inputFields 
00438      *
00439      *   Cross-check (2): outputFields(C,F,P,D,D) vs. inputData(C,P,D): dimensions C, P, D must match
00440      */
00441     TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg,
00442                                                outputFields, 0,2,3,4,
00443                                                inputData,    0,1,2,2),
00444                         std::invalid_argument, errmsg);    
00445     /*
00446      *   Cross-checks (1,3):
00447      */
00448     if( inputFields.rank() == 4) {
00449       // Cross-check (1): inputData(C,P,D) vs. inputFields(C,F,P,D):  dimensions  C, P, D must match
00450       TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, 
00451                                                  inputData,    0,1,2, 
00452                                                  inputFields,  0,2,3),
00453                           std::invalid_argument, errmsg);  
00454       // Cross-check (3): outputFields(C,F,P,D,D) vs. inputFields(C,F,P,D): dimensions C, F, P, D must match
00455       TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, 
00456                                                  outputFields, 0,1,2,3,4,
00457                                                  inputFields,  0,1,2,3,3),
00458                           std::invalid_argument, errmsg);
00459     }
00460     else{
00461       // Cross-check (1): inputData(C,P,D) vs. inputFields(F,P,D): dimensions  P, D must match
00462       TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, 
00463                                                  inputData,    1,2, 
00464                                                  inputFields,  1,2),
00465                           std::invalid_argument, errmsg);      
00466       // Cross-check (3): outputFields(C,F,P,D,D) vs. inputFields(F,P,D): dimensions F, P, D must match
00467       TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, 
00468                                                  outputFields, 1,2,3,4, 
00469                                                  inputFields,  0,1,2,2),
00470                           std::invalid_argument, errmsg);
00471     }
00472 #endif  
00473     
00474     // inputFields is (C,F,P,D)
00475     if(inputFields.rank() == 4){
00476       
00477       for(int cell = 0; cell < outputFields.dimension(0); cell++){
00478         for(int field = 0; field < outputFields.dimension(1); field++){
00479           for(int point = 0; point < outputFields.dimension(2); point++){
00480             for(int row = 0; row < outputFields.dimension(3); row++){
00481               for(int col = 0; col < outputFields.dimension(4); col++){
00482                 outputFields(cell, field, point, row, col) = \
00483                   inputData(cell, point, row)*inputFields(cell, field, point, col);
00484               }// col
00485             }// row
00486           }// point
00487         }// field
00488       } // cell
00489     }// rank = 4
00490      // inputFields is (F,P,D)
00491     else if(inputFields.rank() == 3){
00492       
00493       for(int cell = 0; cell < outputFields.dimension(0); cell++){
00494         for(int field = 0; field < outputFields.dimension(1); field++){
00495           for(int point = 0; point < outputFields.dimension(2); point++){
00496             for(int row = 0; row < outputFields.dimension(3); row++){
00497               for(int col = 0; col < outputFields.dimension(4); col++){
00498                 outputFields(cell, field, point, row, col) = \
00499                   inputData(cell, point, row)*inputFields(field, point, col);
00500               }// col
00501             }// row
00502           }// point
00503         }// field
00504       } // cell
00505     }// rank = 3
00506     else{
00507       TEUCHOS_TEST_FOR_EXCEPTION( true, std::invalid_argument, 
00508                           ">>> ERROR (ArrayTools::outerProductDataField): inputFields rank 3 or 4 required.")
00509     }    
00510   }
00511   
00512   
00513   
00514   template<class Scalar, 
00515            class ArrayOutData, 
00516            class ArrayInDataLeft, 
00517            class ArrayInDataRight>
00518   void ArrayTools::outerProductDataData(ArrayOutData &            outputData,
00519                                         const ArrayInDataLeft &   inputDataLeft,
00520                                         const ArrayInDataRight &  inputDataRight){
00521 #ifdef HAVE_INTREPID_DEBUG
00522     std::string errmsg = ">>> ERROR (ArrayTools::outerProductDataData):";
00523     /*
00524      *   Check array rank and spatial dimension range (if applicable)
00525      *      (1) inputDataLeft(C,P,D);    
00526      *      (2) inputDataRight(C,P,D) or (P,D);   
00527      *      (3) outputData(C,P,D,D)
00528      */
00529     // (1) inputDataLeft is (C, P, D) and 2 <= D <= 3 is required  
00530     TEUCHOS_TEST_FOR_EXCEPTION( !requireRankRange(errmsg, inputDataLeft,  3,3), 
00531                         std::invalid_argument, errmsg);
00532     TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, inputDataLeft, 2,  2,3), 
00533                         std::invalid_argument, errmsg);
00534     // (2) inputDataRight is (C, P, D) or (P, D) and 2 <= (D=dimension(rank - 1)) <= 3 is required. 
00535     TEUCHOS_TEST_FOR_EXCEPTION( !requireRankRange(errmsg, inputDataRight,  2,3), 
00536                         std::invalid_argument, errmsg);
00537     TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, inputDataRight, inputDataRight.rank()-1,  2,3), 
00538                         std::invalid_argument, errmsg);
00539     // (3) outputData is (C,P,D,D)
00540     TEUCHOS_TEST_FOR_EXCEPTION( !requireRankRange(errmsg, outputData, 4, 4), std::invalid_argument, errmsg);      
00541     TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, outputData, 2,  2,3), 
00542                         std::invalid_argument, errmsg);
00543     TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, outputData, 3,  2,3), 
00544                         std::invalid_argument, errmsg);
00545     /*
00546      *   Dimension cross-checks:
00547      *      (1) inputDataLeft vs. inputDataRight
00548      *      (2) outputData    vs. inputDataLeft
00549      *      (3) outputData    vs. inputDataRight 
00550      *
00551      *   Cross-check (2): outputData(C,P,D,D) vs. inputDataLeft(C,P,D): dimensions C, P, D must match
00552      */
00553     TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg,
00554                                                outputData,    0,1,2,3,
00555                                                inputDataLeft, 0,1,2,2),
00556                         std::invalid_argument, errmsg);    
00557     /*
00558      *   Cross-checks (1,3):
00559      */
00560     if( inputDataRight.rank() == 3) {
00561       // Cross-check (1): inputDataLeft(C,P,D) vs. inputDataRight(C,P,D):  all dimensions  C, P, D must match
00562       TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, inputDataLeft, inputDataRight),
00563                           std::invalid_argument, errmsg);  
00564       // Cross-check (3): outputData(C,P,D,D) vs. inputDataRight(C,P,D): dimensions C, P, D must match
00565       TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, 
00566                                                  outputData,     0,1,2,3,
00567                                                  inputDataRight, 0,1,2,2),
00568                           std::invalid_argument, errmsg);
00569     }
00570     else{
00571       // Cross-check (1): inputDataLeft(C,P,D) vs. inputDataRight(P,D): dimensions  P, D must match
00572       TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, 
00573                                                  inputDataLeft,  1,2, 
00574                                                  inputDataRight, 0,1),
00575                           std::invalid_argument, errmsg);      
00576       // Cross-check (3): outputData(C,P,D,D) vs. inputDataRight(P,D): dimensions P, D must match
00577       TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, 
00578                                                  outputData,     1,2,3, 
00579                                                  inputDataRight, 0,1,1),
00580                           std::invalid_argument, errmsg);
00581     }
00582 #endif
00583     // inputDataRight is (C,P,D)
00584     if(inputDataRight.rank() == 3){
00585       
00586       for(int cell = 0; cell < inputDataLeft.dimension(0); cell++){
00587         for(int point = 0; point < inputDataLeft.dimension(1); point++){
00588           for(int row = 0; row < inputDataLeft.dimension(2); row++){
00589             for(int col = 0; col < inputDataLeft.dimension(2); col++){
00590               
00591               outputData(cell, point, row, col) = \
00592                 inputDataLeft(cell, point, row)*inputDataRight(cell, point, col); 
00593             }// col
00594           }// row
00595         }// point
00596       } // cell
00597     }// rank = 3
00598      // inputDataRight is (P,D)
00599     else if(inputDataRight.rank() == 2){
00600       
00601       for(int cell = 0; cell < inputDataLeft.dimension(0); cell++){
00602         for(int point = 0; point < inputDataLeft.dimension(1); point++){
00603           for(int row = 0; row < inputDataLeft.dimension(2); row++){
00604             for(int col = 0; col < inputDataLeft.dimension(2); col++){
00605               // 
00606               outputData(cell, point, row, col) = \
00607               inputDataLeft(cell, point, row)*inputDataRight(point, col); 
00608             } // col
00609           } // row
00610         } // point
00611       } // cell
00612     }// rank = 2
00613     else{
00614       TEUCHOS_TEST_FOR_EXCEPTION( true, std::invalid_argument, 
00615                           ">>> ERROR (ArrayTools::crossProductDataData): inputDataRight rank 2 or 3 required.")
00616     }    
00617   }
00618   
00619   
00620   
00621   template<class Scalar, 
00622            class ArrayOutFields, 
00623            class ArrayInData, 
00624            class ArrayInFields>
00625   void ArrayTools::matvecProductDataField(ArrayOutFields &       outputFields,
00626                                           const ArrayInData &    inputData,
00627                                           const ArrayInFields &  inputFields,
00628                                           const char              transpose){
00629 #ifdef HAVE_INTREPID_DEBUG
00630     std::string errmsg = ">>> ERROR (ArrayTools::matvecProductDataField):";
00631     /*
00632      *   Check array rank and spatial dimension range (if applicable)
00633      *      (1) inputData(C,P), (C,P,D) or (C,P,D,D);   P=1 is admissible to allow multiply by const. data
00634      *      (2) inputFields(C,F,P,D) or (F,P,D);   
00635      *      (3) outputFields(C,F,P,D)
00636      */
00637     // (1) inputData is (C,P), (C, P, D) or (C, P, D, D) and 1 <= D <= 3 is required  
00638     TEUCHOS_TEST_FOR_EXCEPTION( !requireRankRange(errmsg, inputData,  2,4), 
00639                         std::invalid_argument, errmsg);
00640     if(inputData.rank() > 2) {
00641       TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, inputData, 2,  1,3), 
00642                           std::invalid_argument, errmsg);
00643     }
00644     if(inputData.rank() == 4) {
00645       TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, inputData, 3,  1,3), 
00646                         std::invalid_argument, errmsg);
00647     }
00648     // (2) inputFields is (C, F, P, D) or (F, P, D) and 1 <= (D=dimension(rank - 1)) <= 3 is required. 
00649     TEUCHOS_TEST_FOR_EXCEPTION( !requireRankRange(errmsg, inputFields, 3,4), std::invalid_argument, errmsg);
00650     TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, inputFields, inputFields.rank()-1,  1,3), 
00651                         std::invalid_argument, errmsg);
00652     // (3) outputFields is (C,F,P,D)
00653     TEUCHOS_TEST_FOR_EXCEPTION( !requireRankRange(errmsg, outputFields,  4,4), std::invalid_argument, errmsg);      
00654     TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, outputFields, 3,  1,3), 
00655                         std::invalid_argument, errmsg);
00656     /*
00657      *   Dimension cross-checks:
00658      *      (1) inputData    vs. inputFields
00659      *      (2) outputFields vs. inputData
00660      *      (3) outputFields vs. inputFields
00661      *     
00662      *   Cross-check (2): outputFields(C,F,P,D) vs. inputData(C,P), (C,P,D) or (C,P,D,D): 
00663      *   dimensions C, and D must match in all cases, dimension P must match only when non-constant
00664      *   data is specified (P>1). Do not check P dimensions with constant data, i.e., when P=1 in
00665      *   inputData(C,1,...)
00666      */
00667     if(inputData.dimension(1) > 1){ // check P dimension if P>1 in inputData
00668       TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg,
00669                                                  outputFields, 2,
00670                                                  inputData,    1),
00671                           std::invalid_argument, errmsg);    
00672     }
00673     if(inputData.rank() == 2) { // inputData(C,P) -> C match
00674       TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg,
00675                                                  outputFields, 0,
00676                                                  inputData,    0),
00677                           std::invalid_argument, errmsg);    
00678     }
00679     if(inputData.rank() == 3){ // inputData(C,P,D) -> C, D match
00680       TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg,
00681                                                  outputFields, 0,3,
00682                                                  inputData,    0,2),
00683                           std::invalid_argument, errmsg);    
00684       
00685     }
00686     if(inputData.rank() == 4){ // inputData(C,P,D,D) -> C, D, D match
00687       TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg,
00688                                                  outputFields, 0,3,3,
00689                                                  inputData,    0,2,3),
00690                           std::invalid_argument, errmsg);
00691     }
00692     /*
00693      *   Cross-checks (1,3):
00694      */
00695     if(inputFields.rank() == 4) {      
00696       // Cross-check (1): inputData(C,P), (C,P,D) or (C,P,D,D) vs. inputFields(C,F,P,D):  dimensions  C, P, D must match      
00697       if(inputData.dimension(1) > 1){ // check P dimension if P>1 in inputData
00698         TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg,
00699                                                    inputFields,  2,
00700                                                    inputData,    1),
00701                             std::invalid_argument, errmsg);    
00702       }      
00703       if(inputData.rank() == 2){ // inputData(C,P) -> C match
00704         TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, 
00705                                                    inputData,    0, 
00706                                                    inputFields,  0),
00707                             std::invalid_argument, errmsg);  
00708       }
00709       if(inputData.rank() == 3){  // inputData(C,P,D) -> C, D match
00710         TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, 
00711                                                    inputData,    0,2, 
00712                                                    inputFields,  0,3),
00713                             std::invalid_argument, errmsg);  
00714       }
00715       if(inputData.rank() == 4){   // inputData(C,P,D,D) -> C, D, D match
00716         TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, 
00717                                                    inputData,    0,2,3, 
00718                                                    inputFields,  0,3,3),
00719                             std::invalid_argument, errmsg);  
00720       }
00721       // Cross-check (3): outputFields(C,F,P,D) vs. inputFields(C,F,P,D): all dimensions C, F, P, D must match
00722       TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, outputFields, inputFields),
00723                           std::invalid_argument, errmsg);
00724     }
00725     else{
00726       // Cross-check (1): inputData(C,P), (C,P,D) or (C,P,D,D) vs. inputFields(F,P,D): dimensions  P, D must match
00727       if(inputData.dimension(1) > 1){ // check P if P>1 in inputData 
00728         TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, 
00729                                                    inputData,    1, 
00730                                                    inputFields,  1),
00731                             std::invalid_argument, errmsg);    
00732       }
00733       if(inputData.rank() == 3){
00734         TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, 
00735                                                    inputData,    2, 
00736                                                    inputFields,  2),
00737                             std::invalid_argument, errmsg);    
00738       }
00739       if(inputData.rank() == 4){
00740         TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, 
00741                                                    inputData,    2,3, 
00742                                                    inputFields,  2,2),
00743                             std::invalid_argument, errmsg);            
00744       }
00745       // Cross-check (3): outputFields(C,F,P,D) vs. inputFields(F,P,D): dimensions F, P, D must match
00746       TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, 
00747                                                  outputFields, 1,2,3, 
00748                                                  inputFields,  0,1,2),
00749                           std::invalid_argument, errmsg);
00750     }
00751 #endif
00752     int dataRank   = inputData.rank();
00753     int numDataPts = inputData.dimension(1);
00754     int inRank     = inputFields.rank();    
00755     int numCells   = outputFields.dimension(0);
00756     int numFields  = outputFields.dimension(1);
00757     int numPoints  = outputFields.dimension(2);
00758     int matDim     = outputFields.dimension(3);
00759     /*********************************************************************************************
00760      *                              inputFields is (C,F,P,D)                                     *
00761      *********************************************************************************************/
00762     if(inRank == 4){
00763       if(numDataPts != 1){  // non-constant data
00764         
00765         switch(dataRank){
00766           case 2:
00767             for(int cell = 0; cell < numCells; cell++) {
00768               for(int field = 0; field < numFields; field++) {
00769                 for(int point = 0; point < numPoints; point++) {
00770                   for( int row = 0; row < matDim; row++) {
00771                     outputFields(cell, field, point, row) = \
00772                       inputData(cell, point)*inputFields(cell, field, point, row);
00773                   } // Row-loop
00774                 } // P-loop
00775               } // F-loop
00776             }// C-loop
00777             break;
00778             
00779           case 3:
00780             for(int cell = 0; cell < numCells; cell++) {
00781               for(int field = 0; field < numFields; field++) {
00782                 for(int point = 0; point < numPoints; point++) {
00783                   for( int row = 0; row < matDim; row++) {
00784                     outputFields(cell, field, point, row) = \
00785                       inputData(cell, point, row)*inputFields(cell, field, point, row);
00786                   } // Row-loop
00787                 } // P-loop
00788               } // F-loop
00789             }// C-loop
00790             break;
00791             
00792           case 4:
00793             if ((transpose == 'n') || (transpose == 'N')) {
00794               for(int cell = 0; cell < numCells; cell++){
00795                 for(int field = 0; field < numFields; field++){
00796                   for(int point = 0; point < numPoints; point++){
00797                     for(int row = 0; row < matDim; row++){
00798                       outputFields(cell, field, point, row) = 0.0;
00799                       for(int col = 0; col < matDim; col++){
00800                         outputFields(cell, field, point, row) += \
00801                           inputData(cell, point, row, col)*inputFields(cell, field, point, col);
00802                       }// col
00803                     } //row
00804                   }// point
00805                 }// field
00806               }// cell
00807             } // no transpose
00808             else if ((transpose == 't') || (transpose == 'T')) {
00809               for(int cell = 0; cell < numCells; cell++){
00810                 for(int field = 0; field < numFields; field++){
00811                   for(int point = 0; point < numPoints; point++){
00812                     for(int row = 0; row < matDim; row++){
00813                       outputFields(cell, field, point, row) = 0.0;
00814                       for(int col = 0; col < matDim; col++){
00815                         outputFields(cell, field, point, row) += \
00816                           inputData(cell, point, col, row)*inputFields(cell, field, point, col);
00817                       }// col
00818                     } //row
00819                   }// point
00820                 }// field
00821               }// cell
00822             } //transpose
00823             else {
00824               TEUCHOS_TEST_FOR_EXCEPTION( !( (transpose == 'n') || (transpose == 'N') || (transpose == 't') || (transpose == 'T') ), std::invalid_argument,
00825                                   ">>> ERROR (ArrayTools::matvecProductDataField): The transpose flag must be 'n', 'N', 't' or 'T'.");
00826             }
00827             break;
00828             
00829           default:
00830             TEUCHOS_TEST_FOR_EXCEPTION( !( (dataRank == 2) || (dataRank == 3) || (dataRank == 4) ), std::invalid_argument,
00831                                 ">>> ERROR (ArrayTools::matvecProductDataField): inputData rank 2, 3 or 4 required.")      
00832         } // switch inputData rank
00833       }
00834       else{  // constant data case
00835         switch(dataRank){
00836           case 2:
00837             for(int cell = 0; cell < numCells; cell++) {
00838               for(int field = 0; field < numFields; field++) {
00839                 for(int point = 0; point < numPoints; point++) {
00840                   for( int row = 0; row < matDim; row++) {
00841                     outputFields(cell, field, point, row) = \
00842                       inputData(cell, 0)*inputFields(cell, field, point, row);
00843                   } // Row-loop
00844                 } // P-loop
00845               } // F-loop
00846             }// C-loop
00847             break;
00848             
00849           case 3:
00850             for(int cell = 0; cell < numCells; cell++) {
00851               for(int field = 0; field < numFields; field++) {
00852                 for(int point = 0; point < numPoints; point++) {
00853                   for( int row = 0; row < matDim; row++) {
00854                     outputFields(cell, field, point, row) = \
00855                       inputData(cell, 0, row)*inputFields(cell, field, point, row);
00856                   } // Row-loop
00857                 } // P-loop
00858               } // F-loop
00859             }// C-loop
00860             break;
00861             
00862           case 4:
00863             if ((transpose == 'n') || (transpose == 'N')) {
00864               for(int cell = 0; cell < numCells; cell++){
00865                 for(int field = 0; field < numFields; field++){
00866                   for(int point = 0; point < numPoints; point++){
00867                     for(int row = 0; row < matDim; row++){
00868                       outputFields(cell, field, point, row) = 0.0;
00869                       for(int col = 0; col < matDim; col++){
00870                         outputFields(cell, field, point, row) += \
00871                           inputData(cell, 0, row, col)*inputFields(cell, field, point, col);
00872                       }// col
00873                     } //row
00874                   }// point
00875                 }// field
00876               }// cell
00877             } // no transpose
00878             else if ((transpose == 't') || (transpose == 'T')) {
00879               for(int cell = 0; cell < numCells; cell++){
00880                 for(int field = 0; field < numFields; field++){
00881                   for(int point = 0; point < numPoints; point++){
00882                     for(int row = 0; row < matDim; row++){
00883                       outputFields(cell, field, point, row) = 0.0;
00884                       for(int col = 0; col < matDim; col++){
00885                         outputFields(cell, field, point, row) += \
00886                           inputData(cell, 0, col, row)*inputFields(cell, field, point, col);
00887                       }// col
00888                     } //row
00889                   }// point
00890                 }// field
00891               }// cell
00892             } //transpose
00893             else {
00894               TEUCHOS_TEST_FOR_EXCEPTION( !( (transpose == 'n') || (transpose == 'N') || (transpose == 't') || (transpose == 'T') ), std::invalid_argument,
00895                                   ">>> ERROR (ArrayTools::matvecProductDataField): The transpose flag must be 'n', 'N', 't' or 'T'.");
00896             }
00897             break;
00898             
00899           default:
00900             TEUCHOS_TEST_FOR_EXCEPTION( !( (dataRank == 2) || (dataRank == 3) || (dataRank == 4) ), std::invalid_argument,
00901                                 ">>> ERROR (ArrayTools::matvecProductDataField): inputData rank 2, 3 or 4 required.")      
00902         } // switch inputData rank
00903       } // end constant data case
00904     } // inputFields rank 4
00905     /*********************************************************************************************
00906      *                              inputFields is (F,P,D)                                       *
00907      *********************************************************************************************/
00908     else if(inRank == 3) {
00909       if(numDataPts != 1){  // non-constant data
00910         
00911         switch(dataRank){
00912           case 2:
00913             for(int cell = 0; cell < numCells; cell++) {
00914               for(int field = 0; field < numFields; field++) {
00915                 for(int point = 0; point < numPoints; point++) {
00916                   for( int row = 0; row < matDim; row++) {
00917                     outputFields(cell, field, point, row) = \
00918                       inputData(cell, point)*inputFields(field, point, row);
00919                   } // Row-loop
00920                 } // P-loop
00921               } // F-loop
00922             }// C-loop
00923             break;
00924             
00925           case 3:
00926             for(int cell = 0; cell < numCells; cell++) {
00927               for(int field = 0; field < numFields; field++) {
00928                 for(int point = 0; point < numPoints; point++) {
00929                   for( int row = 0; row < matDim; row++) {
00930                     outputFields(cell, field, point, row) = \
00931                       inputData(cell, point, row)*inputFields(field, point, row);
00932                   } // Row-loop
00933                 } // P-loop
00934               } // F-loop
00935             }// C-loop
00936             break;
00937             
00938           case 4:
00939             if ((transpose == 'n') || (transpose == 'N')) {
00940               for(int cell = 0; cell < numCells; cell++){
00941                 for(int field = 0; field < numFields; field++){
00942                   for(int point = 0; point < numPoints; point++){
00943                     for(int row = 0; row < matDim; row++){
00944                       outputFields(cell, field, point, row) = 0.0;
00945                       for(int col = 0; col < matDim; col++){
00946                         outputFields(cell, field, point, row) += \
00947                           inputData(cell, point, row, col)*inputFields(field, point, col);
00948                       }// col
00949                     } //row
00950                   }// point
00951                 }// field
00952               }// cell
00953             } // no transpose
00954             else if ((transpose == 't') || (transpose == 'T')) {
00955               for(int cell = 0; cell < numCells; cell++){
00956                 for(int field = 0; field < numFields; field++){
00957                   for(int point = 0; point < numPoints; point++){
00958                     for(int row = 0; row < matDim; row++){
00959                       outputFields(cell, field, point, row) = 0.0;
00960                       for(int col = 0; col < matDim; col++){
00961                         outputFields(cell, field, point, row) += \
00962                           inputData(cell, point, col, row)*inputFields(field, point, col);
00963                       }// col
00964                     } //row
00965                   }// point
00966                 }// field
00967               }// cell
00968             } //transpose
00969             else {
00970               TEUCHOS_TEST_FOR_EXCEPTION( !( (transpose == 'n') || (transpose == 'N') || (transpose == 't') || (transpose == 'T') ), std::invalid_argument,
00971                                   ">>> ERROR (ArrayTools::matvecProductDataField): The transpose flag must be 'n', 'N', 't' or 'T'.");
00972             }
00973             break;
00974             
00975           default:
00976             TEUCHOS_TEST_FOR_EXCEPTION( !( (dataRank == 2) || (dataRank == 3) || (dataRank == 4) ), std::invalid_argument,
00977                                 ">>> ERROR (ArrayTools::matvecProductDataField): inputData rank 2, 3 or 4 required.")      
00978         } // switch inputData rank
00979       }
00980       else{  // constant data case
00981         switch(dataRank){
00982           case 2:
00983             for(int cell = 0; cell < numCells; cell++) {
00984               for(int field = 0; field < numFields; field++) {
00985                 for(int point = 0; point < numPoints; point++) {
00986                   for( int row = 0; row < matDim; row++) {
00987                     outputFields(cell, field, point, row) = \
00988                     inputData(cell, 0)*inputFields(field, point, row);
00989                   } // Row-loop
00990                 } // P-loop
00991               } // F-loop
00992             }// C-loop
00993             break;
00994             
00995           case 3:
00996             for(int cell = 0; cell < numCells; cell++) {
00997               for(int field = 0; field < numFields; field++) {
00998                 for(int point = 0; point < numPoints; point++) {
00999                   for( int row = 0; row < matDim; row++) {
01000                     outputFields(cell, field, point, row) = \
01001                     inputData(cell, 0, row)*inputFields(field, point, row);
01002                   } // Row-loop
01003                 } // P-loop
01004               } // F-loop
01005             }// C-loop
01006             break;
01007             
01008           case 4:
01009             if ((transpose == 'n') || (transpose == 'N')) {
01010               for(int cell = 0; cell < numCells; cell++){
01011                 for(int field = 0; field < numFields; field++){
01012                   for(int point = 0; point < numPoints; point++){
01013                     for(int row = 0; row < matDim; row++){
01014                       outputFields(cell, field, point, row) = 0.0;
01015                       for(int col = 0; col < matDim; col++){
01016                         outputFields(cell, field, point, row) += \
01017                         inputData(cell, 0, row, col)*inputFields(field, point, col);
01018                       }// col
01019                     } //row
01020                   }// point
01021                 }// field
01022               }// cell
01023             } // no transpose
01024             else if ((transpose == 't') || (transpose == 'T')) {
01025               for(int cell = 0; cell < numCells; cell++){
01026                 for(int field = 0; field < numFields; field++){
01027                   for(int point = 0; point < numPoints; point++){
01028                     for(int row = 0; row < matDim; row++){
01029                       outputFields(cell, field, point, row) = 0.0;
01030                       for(int col = 0; col < matDim; col++){
01031                         outputFields(cell, field, point, row) += \
01032                         inputData(cell, 0, col, row)*inputFields(field, point, col);
01033                       }// col
01034                     } //row
01035                   }// point
01036                 }// field
01037               }// cell
01038             } //transpose
01039             else {
01040               TEUCHOS_TEST_FOR_EXCEPTION( !( (transpose == 'n') || (transpose == 'N') || (transpose == 't') || (transpose == 'T') ), std::invalid_argument,
01041                                   ">>> ERROR (ArrayTools::matvecProductDataField): The transpose flag must be 'n', 'N', 't' or 'T'.");
01042             }
01043             break;
01044             
01045           default:
01046             TEUCHOS_TEST_FOR_EXCEPTION( !( (dataRank == 2) || (dataRank == 3) || (dataRank == 4) ), std::invalid_argument,
01047                                 ">>> ERROR (ArrayTools::matvecProductDataField): inputData rank 2, 3 or 4 required.")      
01048         } // switch inputData rank
01049       } // end constant data case
01050     } // inputFields rank 3
01051     else {
01052       TEUCHOS_TEST_FOR_EXCEPTION( true, std::invalid_argument, 
01053                           ">>> ERROR (ArrayTools::matvecProductDataField): inputFields rank 3 or 4 required.")      
01054     }// rank error
01055   }
01056 
01057   
01058   
01059   template<class Scalar, 
01060            class ArrayOutData, 
01061            class ArrayInDataLeft, 
01062            class ArrayInDataRight>
01063   void ArrayTools::matvecProductDataData(ArrayOutData &            outputData,
01064                                          const ArrayInDataLeft &   inputDataLeft,
01065                                          const ArrayInDataRight &  inputDataRight,
01066                                          const char                transpose){
01067 #ifdef HAVE_INTREPID_DEBUG
01068     std::string errmsg = ">>> ERROR (ArrayTools::matvecProductDataData):";
01069     /*
01070      *   Check array rank and spatial dimension range (if applicable)
01071      *      (1) inputDataLeft(C,P), (C,P,D) or (C,P,D,D); P=1 is admissible to allow multiply by const. left data   
01072      *      (2) inputDataRight(C,P,D) or (P,D);   
01073      *      (3) outputData(C,P,D)
01074      */
01075     // (1) inputDataLeft is (C,P), (C,P,D) or (C,P,D,D) and 1 <= D <= 3 is required  
01076     TEUCHOS_TEST_FOR_EXCEPTION( !requireRankRange(errmsg, inputDataLeft,  2,4), 
01077                         std::invalid_argument, errmsg);
01078     if(inputDataLeft.rank() > 2) {
01079       TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, inputDataLeft, 2,  1,3), 
01080                           std::invalid_argument, errmsg);
01081     }
01082     if(inputDataLeft.rank() == 4) {
01083       TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, inputDataLeft, 3,  1,3), 
01084                           std::invalid_argument, errmsg);
01085     }
01086     // (2) inputDataRight is (C, P, D) or (P, D) and 1 <= (D=dimension(rank - 1)) <= 3 is required. 
01087     TEUCHOS_TEST_FOR_EXCEPTION( !requireRankRange(errmsg, inputDataRight,  2,3), 
01088                         std::invalid_argument, errmsg);
01089     TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, inputDataRight, inputDataRight.rank()-1,  1,3), 
01090                         std::invalid_argument, errmsg);
01091     // (3) outputData is (C,P,D)
01092     TEUCHOS_TEST_FOR_EXCEPTION( !requireRankRange(errmsg, outputData, 3,3), std::invalid_argument, errmsg);
01093     TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, outputData, 2,  1,3), 
01094                         std::invalid_argument, errmsg);
01095     /*
01096      *   Dimension cross-checks:
01097      *      (1) inputDataLeft vs. inputDataRight
01098      *      (2) outputData    vs. inputDataLeft
01099      *      (3) outputData    vs. inputDataRight 
01100      *
01101      *   Cross-check (2): outputData(C,P,D) vs. inputDataLeft(C,P), (C,P,D) or (C,P,D,D):
01102      *   dimensions C, and D must match in all cases, dimension P must match only when non-constant
01103      *   data is specified (P>1). Do not check P dimensions with constant left data, i.e., when P=1 in
01104      *   inputDataLeft(C,1,...)
01105      */
01106     if(inputDataLeft.dimension(1) > 1){ // check P dimension if P>1 in inputDataLeft
01107       TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg,
01108                                                  outputData,     1,
01109                                                  inputDataLeft,  1),
01110                           std::invalid_argument, errmsg);    
01111     }
01112     if(inputDataLeft.rank() == 2){  // inputDataLeft(C,P): check C
01113       TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg,
01114                                                  outputData,    0,
01115                                                  inputDataLeft, 0),
01116                           std::invalid_argument, errmsg);    
01117     }
01118     if(inputDataLeft.rank() == 3){   // inputDataLeft(C,P,D): check C and D
01119       TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg,
01120                                                  outputData,    0,2,
01121                                                  inputDataLeft, 0,2),
01122                           std::invalid_argument, errmsg);    
01123     }
01124     if(inputDataLeft.rank() == 4){   // inputDataLeft(C,P,D,D): check C and D
01125       TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg,
01126                                                  outputData,    0,2,2,
01127                                                  inputDataLeft, 0,2,3),
01128                           std::invalid_argument, errmsg);    
01129     }
01130     /*
01131      *   Cross-checks (1,3):
01132      */
01133     if( inputDataRight.rank() == 3) {
01134       // Cross-check (1): inputDataLeft(C,P), (C,P,D), or (C,P,D,D) vs. inputDataRight(C,P,D):  dimensions  C, P, D must match
01135       if(inputDataLeft.dimension(1) > 1){ // check P dimension if P>1 in inputDataLeft
01136         TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg,
01137                                                    inputDataLeft,  1,
01138                                                    inputDataRight, 1),
01139                             std::invalid_argument, errmsg);    
01140       }      
01141       if(inputDataLeft.rank() == 2){  // inputDataLeft(C,P): check C
01142         TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, 
01143                                                    inputDataLeft,  0, 
01144                                                    inputDataRight, 0),
01145                             std::invalid_argument, errmsg);  
01146       }      
01147       if(inputDataLeft.rank() == 3){   // inputDataLeft(C,P,D): check C and D
01148         TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, 
01149                                                    inputDataLeft,  0,2, 
01150                                                    inputDataRight, 0,2),
01151                             std::invalid_argument, errmsg);  
01152       }      
01153       if(inputDataLeft.rank() == 4){   // inputDataLeft(C,P,D,D): check C and D
01154         TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, 
01155                                                    inputDataLeft,  0,2,3, 
01156                                                    inputDataRight, 0,2,2),
01157                             std::invalid_argument, errmsg);  
01158       }
01159       
01160       // Cross-check (3): outputData(C,P,D) vs. inputDataRight(C,P,D): all dimensions C, P, D must match
01161       TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, outputData, inputDataRight),
01162                           std::invalid_argument, errmsg);
01163     }
01164     else{
01165       // Cross-check (1): inputDataLeft(C,P), (C,P,D), or (C,P,D,D) vs. inputDataRight(P,D): dimensions  P, D must match
01166       if(inputDataLeft.dimension(1) > 1){ // check P if P>1 in inputData 
01167         TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg,
01168                                                    inputDataLeft,  1,
01169                                                    inputDataRight, 0),
01170                             std::invalid_argument, errmsg);    
01171       }
01172       if(inputDataLeft.rank() == 3){   // inputDataLeft(C,P,D): check D
01173         TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, 
01174                                                    inputDataLeft,  2, 
01175                                                    inputDataRight, 1),
01176                             std::invalid_argument, errmsg); 
01177       }
01178       if(inputDataLeft.rank() == 4){   // inputDataLeft(C,P,D,D): check D      
01179         TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, 
01180                                                    inputDataLeft,  2,3, 
01181                                                    inputDataRight, 1,1),
01182                             std::invalid_argument, errmsg); 
01183       }
01184       // Cross-check (3): outputData(C,P,D) vs. inputDataRight(P,D): dimensions P, D must match
01185       TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, 
01186                                                  outputData,     1,2, 
01187                                                  inputDataRight, 0,1),
01188                           std::invalid_argument, errmsg);
01189     }
01190 #endif
01191     int dataLeftRank   = inputDataLeft.rank();
01192     int numDataLeftPts = inputDataLeft.dimension(1);
01193     int dataRightRank  = inputDataRight.rank();    
01194     int numCells       = outputData.dimension(0);
01195     int numPoints      = outputData.dimension(1);
01196     int matDim         = outputData.dimension(2);
01197     
01198     /*********************************************************************************************
01199       *                              inputDataRight is (C,P,D)                                   *
01200       *********************************************************************************************/
01201     if(dataRightRank == 3){
01202       if(numDataLeftPts != 1){  // non-constant left data
01203         
01204         switch(dataLeftRank){
01205           case 2:
01206             for(int cell = 0; cell < numCells; cell++) {
01207               for(int point = 0; point < numPoints; point++) {
01208                 for( int row = 0; row < matDim; row++) {
01209                   outputData(cell, point, row) = \
01210                   inputDataLeft(cell, point)*inputDataRight(cell, point, row);
01211                 } // Row-loop
01212               } // P-loop
01213             }// C-loop
01214             break;
01215             
01216           case 3:
01217             for(int cell = 0; cell < numCells; cell++) {
01218               for(int point = 0; point < numPoints; point++) {
01219                 for( int row = 0; row < matDim; row++) {
01220                   outputData(cell, point, row) = \
01221                   inputDataLeft(cell, point, row)*inputDataRight(cell, point, row);
01222                 } // Row-loop
01223               } // P-loop
01224             }// C-loop
01225             break;
01226             
01227           case 4:
01228             if ((transpose == 'n') || (transpose == 'N')) {
01229               for(int cell = 0; cell < numCells; cell++){
01230                 for(int point = 0; point < numPoints; point++){
01231                   for(int row = 0; row < matDim; row++){
01232                     outputData(cell, point, row) = 0.0;
01233                     for(int col = 0; col < matDim; col++){
01234                       outputData(cell, point, row) += \
01235                       inputDataLeft(cell, point, row, col)*inputDataRight(cell, point, col);
01236                     }// col
01237                   } //row
01238                 }// point
01239               }// cell
01240             } // no transpose
01241             else if ((transpose == 't') || (transpose == 'T')) {
01242               for(int cell = 0; cell < numCells; cell++){
01243                 for(int point = 0; point < numPoints; point++){
01244                   for(int row = 0; row < matDim; row++){
01245                     outputData(cell, point, row) = 0.0;
01246                     for(int col = 0; col < matDim; col++){
01247                       outputData(cell, point, row) += \
01248                       inputDataLeft(cell, point, col, row)*inputDataRight(cell, point, col);
01249                     }// col
01250                   } //row
01251                 }// point
01252               }// cell
01253             } //transpose
01254             else {
01255               TEUCHOS_TEST_FOR_EXCEPTION( !( (transpose == 'n') || (transpose == 'N') || (transpose == 't') || (transpose == 'T') ), std::invalid_argument,
01256                                   ">>> ERROR (ArrayTools::matvecProductDataData): The transpose flag must be 'n', 'N', 't' or 'T'.");
01257             }
01258             break;
01259             
01260           default:
01261             TEUCHOS_TEST_FOR_EXCEPTION( !( (dataLeftRank == 2) || (dataLeftRank == 3) || (dataLeftRank == 4) ), std::invalid_argument,
01262                                 ">>> ERROR (ArrayTools::matvecProductDataData): inputDataLeft rank 2, 3 or 4 required.")      
01263         } // switch inputDataLeft rank
01264       }
01265       else{  // constant data case
01266         switch(dataLeftRank){
01267           case 2:
01268             for(int cell = 0; cell < numCells; cell++) {
01269               for(int point = 0; point < numPoints; point++) {
01270                 for( int row = 0; row < matDim; row++) {
01271                   outputData(cell, point, row) = \
01272                   inputDataLeft(cell, 0)*inputDataRight(cell, point, row);
01273                 } // Row-loop
01274               } // F-loop
01275             }// C-loop
01276             break;
01277             
01278           case 3:
01279             for(int cell = 0; cell < numCells; cell++) {
01280               for(int point = 0; point < numPoints; point++) {
01281                 for( int row = 0; row < matDim; row++) {
01282                   outputData(cell, point, row) = \
01283                   inputDataLeft(cell, 0, row)*inputDataRight(cell, point, row);
01284                 } // Row-loop
01285               } // P-loop
01286             }// C-loop
01287             break;
01288             
01289           case 4:
01290             if ((transpose == 'n') || (transpose == 'N')) {
01291               for(int cell = 0; cell < numCells; cell++){
01292                 for(int point = 0; point < numPoints; point++){
01293                   for(int row = 0; row < matDim; row++){
01294                     outputData(cell, point, row) = 0.0;
01295                     for(int col = 0; col < matDim; col++){
01296                       outputData(cell, point, row) += \
01297                       inputDataLeft(cell, 0, row, col)*inputDataRight(cell, point, col);
01298                     }// col
01299                   } //row
01300                 }// point
01301               }// cell
01302             } // no transpose
01303             else if ((transpose == 't') || (transpose == 'T')) {
01304               for(int cell = 0; cell < numCells; cell++){
01305                 for(int point = 0; point < numPoints; point++){
01306                   for(int row = 0; row < matDim; row++){
01307                     outputData(cell, point, row) = 0.0;
01308                     for(int col = 0; col < matDim; col++){
01309                       outputData(cell, point, row) += \
01310                       inputDataLeft(cell, 0, col, row)*inputDataRight(cell, point, col);
01311                     }// col
01312                   } //row
01313                 }// point
01314               }// cell
01315             } //transpose
01316             else {
01317               TEUCHOS_TEST_FOR_EXCEPTION( !( (transpose == 'n') || (transpose == 'N') || (transpose == 't') || (transpose == 'T') ), std::invalid_argument,
01318                                   ">>> ERROR (ArrayTools::matvecProductDataData): The transpose flag must be 'n', 'N', 't' or 'T'.");
01319             }
01320             break;
01321             
01322           default:
01323             TEUCHOS_TEST_FOR_EXCEPTION( !( (dataLeftRank == 2) || (dataLeftRank == 3) || (dataLeftRank == 4) ), std::invalid_argument,
01324                                 ">>> ERROR (ArrayTools::matvecProductDataData): inputDataLeft rank 2, 3 or 4 required.")      
01325         } // switch inputDataLeft rank
01326       } // end constant data case
01327     } // inputDataRight rank 4
01328     /*********************************************************************************************
01329       *                              inputDataRight is (P,D)                                     *
01330       *********************************************************************************************/
01331     else if(dataRightRank == 2) {
01332       if(numDataLeftPts != 1){  // non-constant data
01333         
01334         switch(dataLeftRank){
01335           case 2:
01336             for(int cell = 0; cell < numCells; cell++) {
01337               for(int point = 0; point < numPoints; point++) {
01338                 for( int row = 0; row < matDim; row++) {
01339                   outputData(cell, point, row) = \
01340                   inputDataLeft(cell, point)*inputDataRight(point, row);
01341                 } // Row-loop
01342               } // P-loop
01343             }// C-loop
01344             break;
01345             
01346           case 3:
01347             for(int cell = 0; cell < numCells; cell++) {
01348               for(int point = 0; point < numPoints; point++) {
01349                 for( int row = 0; row < matDim; row++) {
01350                   outputData(cell, point, row) = \
01351                   inputDataLeft(cell, point, row)*inputDataRight(point, row);
01352                 } // Row-loop
01353               } // P-loop
01354             }// C-loop
01355             break;
01356             
01357           case 4:
01358             if ((transpose == 'n') || (transpose == 'N')) {
01359               for(int cell = 0; cell < numCells; cell++){
01360                 for(int point = 0; point < numPoints; point++){
01361                   for(int row = 0; row < matDim; row++){
01362                     outputData(cell, point, row) = 0.0;
01363                     for(int col = 0; col < matDim; col++){
01364                       outputData(cell, point, row) += \
01365                       inputDataLeft(cell, point, row, col)*inputDataRight(point, col);
01366                     }// col
01367                   } //row
01368                 }// point
01369               }// cell
01370             } // no transpose
01371             else if ((transpose == 't') || (transpose == 'T')) {
01372               for(int cell = 0; cell < numCells; cell++){
01373                 for(int point = 0; point < numPoints; point++){
01374                   for(int row = 0; row < matDim; row++){
01375                     outputData(cell, point, row) = 0.0;
01376                     for(int col = 0; col < matDim; col++){
01377                       outputData(cell, point, row) += \
01378                       inputDataLeft(cell, point, col, row)*inputDataRight(point, col);
01379                     }// col
01380                   } //row
01381                 }// point
01382               }// cell
01383             } //transpose
01384             else {
01385               TEUCHOS_TEST_FOR_EXCEPTION( !( (transpose == 'n') || (transpose == 'N') || (transpose == 't') || (transpose == 'T') ), std::invalid_argument,
01386                                   ">>> ERROR (ArrayTools::matvecProductDataData): The transpose flag must be 'n', 'N', 't' or 'T'.");
01387             }
01388             break;
01389             
01390           default:
01391             TEUCHOS_TEST_FOR_EXCEPTION( !( (dataLeftRank == 2) || (dataLeftRank == 3) || (dataLeftRank == 4) ), std::invalid_argument,
01392                                 ">>> ERROR (ArrayTools::matvecProductDataData): inputDataLeft rank 2, 3 or 4 required.")      
01393         } // switch inputDataLeft rank
01394       }
01395       else{  // constant data case
01396         switch(dataLeftRank){
01397           case 2:
01398             for(int cell = 0; cell < numCells; cell++) {
01399               for(int point = 0; point < numPoints; point++) {
01400                 for( int row = 0; row < matDim; row++) {
01401                   outputData(cell, point, row) = \
01402                   inputDataLeft(cell, 0)*inputDataRight(point, row);
01403                 } // Row-loop
01404               } // P-loop
01405             }// C-loop
01406             break;
01407             
01408           case 3:
01409             for(int cell = 0; cell < numCells; cell++) {
01410               for(int point = 0; point < numPoints; point++) {
01411                 for( int row = 0; row < matDim; row++) {
01412                   outputData(cell, point, row) = \
01413                   inputDataLeft(cell, 0, row)*inputDataRight(point, row);
01414                 } // Row-loop
01415               } // P-loop
01416             }// C-loop
01417             break;
01418             
01419           case 4:
01420             if ((transpose == 'n') || (transpose == 'N')) {
01421               for(int cell = 0; cell < numCells; cell++){
01422                 for(int point = 0; point < numPoints; point++){
01423                   for(int row = 0; row < matDim; row++){
01424                     outputData(cell, point, row) = 0.0;
01425                     for(int col = 0; col < matDim; col++){
01426                       outputData(cell, point, row) += \
01427                       inputDataLeft(cell, 0, row, col)*inputDataRight(point, col);
01428                     }// col
01429                   } //row
01430                 }// point
01431               }// cell
01432             } // no transpose
01433             else if ((transpose == 't') || (transpose == 'T')) {
01434               for(int cell = 0; cell < numCells; cell++){
01435                 for(int point = 0; point < numPoints; point++){
01436                   for(int row = 0; row < matDim; row++){
01437                     outputData(cell, point, row) = 0.0;
01438                     for(int col = 0; col < matDim; col++){
01439                       outputData(cell, point, row) += \
01440                       inputDataLeft(cell, 0, col, row)*inputDataRight(point, col);
01441                     }// col
01442                   } //row
01443                 }// point
01444               }// cell
01445             } //transpose
01446             else {
01447               TEUCHOS_TEST_FOR_EXCEPTION( !( (transpose == 'n') || (transpose == 'N') || (transpose == 't') || (transpose == 'T') ), std::invalid_argument,
01448                                   ">>> ERROR (ArrayTools::matvecProductDataData): The transpose flag must be 'n', 'N', 't' or 'T'.");
01449             }
01450             break;
01451             
01452           default:
01453             TEUCHOS_TEST_FOR_EXCEPTION( !( (dataLeftRank == 2) || (dataLeftRank == 3) || (dataLeftRank == 4) ), std::invalid_argument,
01454                                 ">>> ERROR (ArrayTools::matvecProductDataData): inputDataLeft rank 2, 3 or 4 required.")      
01455         } // switch inputDataLeft rank
01456       } // end constant inputDataLeft case
01457     } // inputDataRight rank 2
01458     else {
01459       TEUCHOS_TEST_FOR_EXCEPTION( true, std::invalid_argument, 
01460                           ">>> ERROR (ArrayTools::matvecProductDataData): inputDataRight rank 2 or 3 required.")      
01461     }// rank error
01462   }
01463   
01464   
01465   
01466   template<class Scalar, 
01467            class ArrayOutFields, 
01468            class ArrayInData, 
01469            class ArrayInFields>
01470   void ArrayTools::matmatProductDataField(ArrayOutFields &       outputFields,
01471                                           const ArrayInData &    inputData,
01472                                           const ArrayInFields &  inputFields,
01473                                           const char             transpose){
01474 #ifdef HAVE_INTREPID_DEBUG
01475     std::string errmsg = ">>> ERROR (ArrayTools::matmatProductDataField):";
01476     /*
01477      *   Check array rank and spatial dimension range (if applicable)
01478      *      (1) inputData(C,P), (C,P,D) or (C,P,D,D);   P=1 is admissible to allow multiply by const. data
01479      *      (2) inputFields(C,F,P,D,D) or (F,P,D,D);   
01480      *      (3) outputFields(C,F,P,D,D)
01481      */
01482     // (1) inputData is (C,P), (C, P, D) or (C, P, D, D) and 1 <= D <= 3 is required  
01483     TEUCHOS_TEST_FOR_EXCEPTION( !requireRankRange(errmsg, inputData,  2,4), 
01484                         std::invalid_argument, errmsg);
01485     if(inputData.rank() > 2) {
01486       TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, inputData, 2,  1,3), 
01487                           std::invalid_argument, errmsg);
01488     }
01489     if(inputData.rank() == 4) {
01490       TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, inputData, 3,  1,3), 
01491                           std::invalid_argument, errmsg);
01492     }
01493     // (2) inputFields is (C,F,P,D,D) or (F,P,D,D) and 1 <= (dimension(rank-1), (rank-2)) <= 3 is required. 
01494     TEUCHOS_TEST_FOR_EXCEPTION( !requireRankRange(errmsg, inputFields, 4,5), std::invalid_argument, errmsg);
01495     TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, inputFields, inputFields.rank()-1,  1,3), 
01496                         std::invalid_argument, errmsg);
01497     TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, inputFields, inputFields.rank()-2,  1,3), 
01498                         std::invalid_argument, errmsg);
01499     // (3) outputFields is (C,F,P,D,D)
01500     TEUCHOS_TEST_FOR_EXCEPTION( !requireRankRange(errmsg, outputFields,  5,5), std::invalid_argument, errmsg);      
01501     TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, outputFields, 3,  1,3), 
01502                         std::invalid_argument, errmsg);
01503     TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, outputFields, 4,  1,3), 
01504                         std::invalid_argument, errmsg);
01505     /*
01506      *   Dimension cross-checks:
01507      *      (1) inputData    vs. inputFields
01508      *      (2) outputFields vs. inputData
01509      *      (3) outputFields vs. inputFields 
01510      *
01511      *   Cross-check (2): outputFields(C,F,P,D,D) vs. inputData(C,P), (C,P,D) or (C,P,D,D): 
01512      *   dimensions C, and D must match in all cases, dimension P must match only when non-constant
01513      *   data is specified (P>1). Do not check P dimensions with constant data, i.e., when P=1 in
01514      *   inputData(C,1,...)
01515      */
01516     if(inputData.dimension(1) > 1){ // check P dimension if P>1 in inputData
01517       TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg,
01518                                                  outputFields, 2,
01519                                                  inputData,    1),
01520                           std::invalid_argument, errmsg);    
01521     }
01522     if(inputData.rank() == 2) { // inputData(C,P) -> C match
01523       TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg,
01524                                                  outputFields, 0,
01525                                                  inputData,    0),
01526                           std::invalid_argument, errmsg);    
01527     }
01528     if(inputData.rank() == 3){ // inputData(C,P,D) -> C, D match
01529       TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg,
01530                                                  outputFields, 0,3,4,
01531                                                  inputData,    0,2,2),
01532                           std::invalid_argument, errmsg);    
01533       
01534     }
01535     if(inputData.rank() == 4){ // inputData(C,P,D,D) -> C, D, D match
01536       TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg,
01537                                                  outputFields, 0,3,4,
01538                                                  inputData,    0,2,3),
01539                           std::invalid_argument, errmsg);
01540     }
01541     /*
01542      *   Cross-checks (1,3):
01543      */
01544     if( inputFields.rank() == 5) {
01545       // Cross-check (1): inputData(C,P), (C,P,D) or (C,P,D,D) vs. inputFields(C,F,P,D,D):  dimensions  C, P, D must match
01546       if(inputData.dimension(1) > 1){ // check P dimension if P>1 in inputData
01547         TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg,
01548                                                    inputData,    1,
01549                                                    inputFields,  2),
01550                             std::invalid_argument, errmsg);    
01551       }      
01552       if(inputData.rank() == 2){ // inputData(C,P) -> C match
01553         TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, 
01554                                                    inputData,    0, 
01555                                                    inputFields,  0),
01556                             std::invalid_argument, errmsg);  
01557       }
01558       if(inputData.rank() == 3){  // inputData(C,P,D) -> C, D match
01559         TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, 
01560                                                    inputData,    0,2,2, 
01561                                                    inputFields,  0,3,4),
01562                             std::invalid_argument, errmsg);  
01563       }
01564       if(inputData.rank() == 4){   // inputData(C,P,D,D) -> C, D, D match
01565         TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, 
01566                                                    inputData,    0,2,3, 
01567                                                    inputFields,  0,3,4),
01568                             std::invalid_argument, errmsg);  
01569       }
01570       // Cross-check (3): outputFields(C,F,P,D,D) vs. inputFields(C,F,P,D,D): all dimensions C, F, P, D must match
01571       TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, outputFields, inputFields),
01572                           std::invalid_argument, errmsg);
01573     }
01574     else{
01575       // Cross-check (1): inputData(C,P), (C,P,D) or (C,P,D,D) vs. inputFields(F,P,D,D): dimensions  P, D must match
01576       if(inputData.dimension(1) > 1){ // check P if P>1 in inputData 
01577         TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, 
01578                                                    inputData,    1, 
01579                                                    inputFields,  1),
01580                             std::invalid_argument, errmsg);    
01581       }
01582       if(inputData.rank() == 3){
01583         TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, 
01584                                                    inputData,    2,2, 
01585                                                    inputFields,  2,3),
01586                             std::invalid_argument, errmsg);    
01587       }
01588       if(inputData.rank() == 4){
01589         TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, 
01590                                                    inputData,    2,3, 
01591                                                    inputFields,  2,3),
01592                             std::invalid_argument, errmsg);            
01593       }
01594       // Cross-check (3): outputFields(C,F,P,D,D) vs. inputFields(F,P,D,D): dimensions F, P, D must match
01595       TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, 
01596                                                  outputFields, 1,2,3,4, 
01597                                                  inputFields,  0,1,2,3),
01598                           std::invalid_argument, errmsg);
01599     }
01600 #endif
01601     int dataRank   = inputData.rank();
01602     int numDataPts = inputData.dimension(1);
01603     int inRank     = inputFields.rank();    
01604     int numCells   = outputFields.dimension(0);
01605     int numFields  = outputFields.dimension(1);
01606     int numPoints  = outputFields.dimension(2);
01607     int matDim     = outputFields.dimension(3);
01608     
01609     /*********************************************************************************************
01610       *                              inputFields is (C,F,P,D,D)                                     *
01611       *********************************************************************************************/
01612     if(inRank == 5){
01613       if(numDataPts != 1){  // non-constant data
01614         
01615         switch(dataRank){
01616           case 2:
01617             for(int cell = 0; cell < numCells; cell++) {
01618               for(int field = 0; field < numFields; field++) {
01619                 for(int point = 0; point < numPoints; point++) {
01620                   for( int row = 0; row < matDim; row++) {
01621                     for( int col = 0; col < matDim; col++) {
01622                       outputFields(cell, field, point, row, col) = \
01623                       inputData(cell, point)*inputFields(cell, field, point, row, col);
01624                     }// Col-loop
01625                   } // Row-loop
01626                 } // P-loop
01627               } // F-loop
01628             }// C-loop
01629             break;
01630             
01631           case 3:
01632             for(int cell = 0; cell < numCells; cell++) {
01633               for(int field = 0; field < numFields; field++) {
01634                 for(int point = 0; point < numPoints; point++) {
01635                   for( int row = 0; row < matDim; row++) {
01636                     for( int col = 0; col < matDim; col++) {
01637                       outputFields(cell, field, point, row, col) = \
01638                       inputData(cell, point, row)*inputFields(cell, field, point, row, col);
01639                     }// Col-loop
01640                   } // Row-loop
01641                 } // P-loop
01642               } // F-loop
01643             }// C-loop
01644             break;
01645             
01646           case 4:
01647             if ((transpose == 'n') || (transpose == 'N')) {
01648               for(int cell = 0; cell < numCells; cell++){
01649                 for(int field = 0; field < numFields; field++){
01650                   for(int point = 0; point < numPoints; point++){
01651                     for(int row = 0; row < matDim; row++){
01652                       for(int col = 0; col < matDim; col++){
01653                         outputFields(cell, field, point, row, col) = 0.0;
01654                         for(int i = 0; i < matDim; i++){
01655                           outputFields(cell, field, point, row, col) += \
01656                           inputData(cell, point, row, i)*inputFields(cell, field, point, i, col);
01657                         }// i
01658                       } // col
01659                     } //row
01660                   }// point
01661                 }// field
01662               }// cell
01663             } // no transpose
01664             else if ((transpose == 't') || (transpose == 'T')) {
01665               for(int cell = 0; cell < numCells; cell++){
01666                 for(int field = 0; field < numFields; field++){
01667                   for(int point = 0; point < numPoints; point++){
01668                     for(int row = 0; row < matDim; row++){
01669                       for(int col = 0; col < matDim; col++){
01670                         outputFields(cell, field, point, row, col) = 0.0;
01671                         for(int i = 0; i < matDim; i++){
01672                           outputFields(cell, field, point, row, col) += \
01673                           inputData(cell, point, i, row)*inputFields(cell, field, point, i, col);
01674                         }// i
01675                       } // col
01676                     } //row
01677                   }// point
01678                 }// field
01679               }// cell
01680             } //transpose
01681             else {
01682               TEUCHOS_TEST_FOR_EXCEPTION( !( (transpose == 'n') || (transpose == 'N') || (transpose == 't') || (transpose == 'T') ), std::invalid_argument,
01683                                   ">>> ERROR (ArrayTools::matmatProductDataField): The transpose flag must be 'n', 'N', 't' or 'T'.");
01684             }
01685             break;
01686             
01687           default:
01688             TEUCHOS_TEST_FOR_EXCEPTION( !( (dataRank == 2) || (dataRank == 3) || (dataRank == 4) ), std::invalid_argument,
01689                                 ">>> ERROR (ArrayTools::matmatProductDataField): inputData rank 2, 3 or 4 required.")      
01690         } // switch inputData rank
01691       }
01692       else{  // constant data case
01693         switch(dataRank){
01694           case 2:
01695             for(int cell = 0; cell < numCells; cell++) {
01696               for(int field = 0; field < numFields; field++) {
01697                 for(int point = 0; point < numPoints; point++) {
01698                   for( int row = 0; row < matDim; row++) {
01699                     for( int col = 0; col < matDim; col++) {
01700                       outputFields(cell, field, point, row, col) = \
01701                       inputData(cell, 0)*inputFields(cell, field, point, row, col);
01702                     }// Col-loop
01703                   } // Row-loop
01704                 } // P-loop
01705               } // F-loop
01706             }// C-loop
01707             break;
01708             
01709           case 3:
01710             for(int cell = 0; cell < numCells; cell++) {
01711               for(int field = 0; field < numFields; field++) {
01712                 for(int point = 0; point < numPoints; point++) {
01713                   for( int row = 0; row < matDim; row++) {
01714                     for( int col = 0; col < matDim; col++) {
01715                       outputFields(cell, field, point, row, col) = \
01716                       inputData(cell, 0, row)*inputFields(cell, field, point, row, col);
01717                     }// Col-loop
01718                   } // Row-loop
01719                 } // P-loop
01720               } // F-loop
01721             }// C-loop
01722             break;
01723             
01724           case 4:
01725             if ((transpose == 'n') || (transpose == 'N')) {
01726               for(int cell = 0; cell < numCells; cell++){
01727                 for(int field = 0; field < numFields; field++){
01728                   for(int point = 0; point < numPoints; point++){
01729                     for(int row = 0; row < matDim; row++){
01730                       for(int col = 0; col < matDim; col++){
01731                         outputFields(cell, field, point, row, col) = 0.0;
01732                         for(int i = 0; i < matDim; i++){
01733                           outputFields(cell, field, point, row, col) += \
01734                           inputData(cell, 0, row, i)*inputFields(cell, field, point, i, col);
01735                         }// i
01736                       } // col
01737                     } //row
01738                   }// point
01739                 }// field
01740               }// cell
01741             } // no transpose
01742             else if ((transpose == 't') || (transpose == 'T')) {
01743               for(int cell = 0; cell < numCells; cell++){
01744                 for(int field = 0; field < numFields; field++){
01745                   for(int point = 0; point < numPoints; point++){
01746                     for(int row = 0; row < matDim; row++){
01747                       for(int col = 0; col < matDim; col++){
01748                         outputFields(cell, field, point, row, col) = 0.0;
01749                         for(int i = 0; i < matDim; i++){
01750                           outputFields(cell, field, point, row, col) += \
01751                           inputData(cell, 0, i, row)*inputFields(cell, field, point, i, col);
01752                         }// i
01753                       } // col
01754                     } //row
01755                   }// point
01756                 }// field
01757               }// cell
01758             } //transpose
01759             else {
01760               TEUCHOS_TEST_FOR_EXCEPTION( !( (transpose == 'n') || (transpose == 'N') || (transpose == 't') || (transpose == 'T') ), std::invalid_argument,
01761                                   ">>> ERROR (ArrayTools::matmatProductDataField): The transpose flag must be 'n', 'N', 't' or 'T'.");
01762             }
01763             break;
01764             
01765           default:
01766             TEUCHOS_TEST_FOR_EXCEPTION( !( (dataRank == 2) || (dataRank == 3) || (dataRank == 4) ), std::invalid_argument,
01767                                 ">>> ERROR (ArrayTools::matmatProductDataField): inputData rank 2, 3 or 4 required.")      
01768         } // switch inputData rank
01769       } // end constant data case
01770     } // inputFields rank 5
01771     /**********************************************************************************************
01772       *                              inputFields is (F,P,D,D)                                     *
01773       *********************************************************************************************/
01774     else if(inRank == 4) {
01775       if(numDataPts != 1){  // non-constant data
01776         
01777         switch(dataRank){
01778           case 2:
01779             for(int cell = 0; cell < numCells; cell++) {
01780               for(int field = 0; field < numFields; field++) {
01781                 for(int point = 0; point < numPoints; point++) {
01782                   for( int row = 0; row < matDim; row++) {
01783                     for( int col = 0; col < matDim; col++) {
01784                       outputFields(cell, field, point, row, col) = \
01785                       inputData(cell, point)*inputFields(field, point, row, col);
01786                     }// Col-loop
01787                   } // Row-loop
01788                 } // P-loop
01789               } // F-loop
01790             }// C-loop
01791             break;
01792             
01793           case 3:
01794             for(int cell = 0; cell < numCells; cell++) {
01795               for(int field = 0; field < numFields; field++) {
01796                 for(int point = 0; point < numPoints; point++) {
01797                   for( int row = 0; row < matDim; row++) {
01798                     for( int col = 0; col < matDim; col++) {
01799                       outputFields(cell, field, point, row, col) = \
01800                       inputData(cell, point, row)*inputFields(field, point, row, col);
01801                     }// Col-loop
01802                   } // Row-loop
01803                 } // P-loop
01804               } // F-loop
01805             }// C-loop
01806             break;
01807             
01808           case 4:
01809             if ((transpose == 'n') || (transpose == 'N')) {
01810               for(int cell = 0; cell < numCells; cell++){
01811                 for(int field = 0; field < numFields; field++){
01812                   for(int point = 0; point < numPoints; point++){
01813                     for(int row = 0; row < matDim; row++){
01814                       for(int col = 0; col < matDim; col++){
01815                         outputFields(cell, field, point, row, col) = 0.0;
01816                         for(int i = 0; i < matDim; i++){
01817                           outputFields(cell, field, point, row, col) += \
01818                           inputData(cell, point, row, i)*inputFields(field, point, i, col);
01819                         }// i
01820                       } // col
01821                     } //row
01822                   }// point
01823                 }// field
01824               }// cell
01825             } // no transpose
01826             else if ((transpose == 't') || (transpose == 'T')) {
01827               for(int cell = 0; cell < numCells; cell++){
01828                 for(int field = 0; field < numFields; field++){
01829                   for(int point = 0; point < numPoints; point++){
01830                     for(int row = 0; row < matDim; row++){
01831                       for(int col = 0; col < matDim; col++){
01832                         outputFields(cell, field, point, row, col) = 0.0;
01833                         for(int i = 0; i < matDim; i++){
01834                           outputFields(cell, field, point, row, col) += \
01835                           inputData(cell, point, i, row)*inputFields(field, point, i, col);
01836                         }// i
01837                       } // col
01838                     } //row
01839                   }// point
01840                 }// field
01841               }// cell
01842             } //transpose
01843             else {
01844               TEUCHOS_TEST_FOR_EXCEPTION( !( (transpose == 'n') || (transpose == 'N') || (transpose == 't') || (transpose == 'T') ), std::invalid_argument,
01845                                   ">>> ERROR (ArrayTools::matmatProductDataField): The transpose flag must be 'n', 'N', 't' or 'T'.");
01846             }
01847             break;
01848             
01849           default:
01850             TEUCHOS_TEST_FOR_EXCEPTION( !( (dataRank == 2) || (dataRank == 3) || (dataRank == 4) ), std::invalid_argument,
01851                                 ">>> ERROR (ArrayTools::matmatProductDataField): inputData rank 2, 3 or 4 required.")      
01852         } // switch inputData rank
01853       }
01854       else{  // constant data case
01855         switch(dataRank){
01856           case 2:
01857             for(int cell = 0; cell < numCells; cell++) {
01858               for(int field = 0; field < numFields; field++) {
01859                 for(int point = 0; point < numPoints; point++) {
01860                   for( int row = 0; row < matDim; row++) {
01861                     for( int col = 0; col < matDim; col++) {
01862                       outputFields(cell, field, point, row, col) = \
01863                       inputData(cell, 0)*inputFields(field, point, row, col);
01864                     }// Col-loop
01865                   } // Row-loop
01866                 } // P-loop
01867               } // F-loop
01868             }// C-loop
01869             break;
01870             
01871           case 3:
01872             for(int cell = 0; cell < numCells; cell++) {
01873               for(int field = 0; field < numFields; field++) {
01874                 for(int point = 0; point < numPoints; point++) {
01875                   for( int row = 0; row < matDim; row++) {
01876                     for( int col = 0; col < matDim; col++) {
01877                       outputFields(cell, field, point, row, col) = \
01878                       inputData(cell, 0, row)*inputFields(field, point, row, col);
01879                     }// Col-loop
01880                   } // Row-loop
01881                 } // P-loop
01882               } // F-loop
01883             }// C-loop
01884             break;
01885             
01886           case 4:
01887             if ((transpose == 'n') || (transpose == 'N')) {
01888               for(int cell = 0; cell < numCells; cell++){
01889                 for(int field = 0; field < numFields; field++){
01890                   for(int point = 0; point < numPoints; point++){
01891                     for(int row = 0; row < matDim; row++){
01892                       for(int col = 0; col < matDim; col++){
01893                         outputFields(cell, field, point, row, col) = 0.0;
01894                         for(int i = 0; i < matDim; i++){
01895                           outputFields(cell, field, point, row, col) += \
01896                           inputData(cell, 0, row, i)*inputFields(field, point, i, col);
01897                         }// i
01898                       } // col
01899                     } //row
01900                   }// point
01901                 }// field
01902               }// cell
01903             } // no transpose
01904             else if ((transpose == 't') || (transpose == 'T')) {
01905               for(int cell = 0; cell < numCells; cell++){
01906                 for(int field = 0; field < numFields; field++){
01907                   for(int point = 0; point < numPoints; point++){
01908                     for(int row = 0; row < matDim; row++){
01909                       for(int col = 0; col < matDim; col++){
01910                         outputFields(cell, field, point, row, col) = 0.0;
01911                         for(int i = 0; i < matDim; i++){
01912                           outputFields(cell, field, point, row, col) += \
01913                           inputData(cell, 0, i, row)*inputFields(field, point, i, col);
01914                         }// i
01915                       } // col
01916                     } //row
01917                   }// point
01918                 }// field
01919               }// cell
01920             } //transpose
01921             else {
01922               TEUCHOS_TEST_FOR_EXCEPTION( !( (transpose == 'n') || (transpose == 'N') || (transpose == 't') || (transpose == 'T') ), std::invalid_argument,
01923                                   ">>> ERROR (ArrayTools::matmatProductDataField): The transpose flag must be 'n', 'N', 't' or 'T'.");
01924             }
01925             break;
01926             
01927           default:
01928             TEUCHOS_TEST_FOR_EXCEPTION( !( (dataRank == 2) || (dataRank == 3) || (dataRank == 4) ), std::invalid_argument,
01929                                 ">>> ERROR (ArrayTools::matmatProductDataField): inputData rank 2, 3 or 4 required.")      
01930         } // switch inputData rank
01931       } // end constant data case
01932     } // inputFields rank 4
01933     else {
01934       TEUCHOS_TEST_FOR_EXCEPTION( true, std::invalid_argument, 
01935                           ">>> ERROR (ArrayTools::matmatProductDataField): inputFields rank 4 or 5 required.")      
01936     }// rank error
01937   }
01938   
01939   
01940   
01941   template<class Scalar, 
01942            class ArrayOutData, 
01943            class ArrayInDataLeft, 
01944            class ArrayInDataRight>
01945   void ArrayTools::matmatProductDataData(ArrayOutData &            outputData,
01946                                          const ArrayInDataLeft &   inputDataLeft,
01947                                          const ArrayInDataRight &  inputDataRight,
01948                                          const char                transpose){
01949 #ifdef HAVE_INTREPID_DEBUG
01950     std::string errmsg = ">>> ERROR (ArrayTools::matmatProductDataData):";
01951     /*
01952      *   Check array rank and spatial dimension range (if applicable)
01953      *      (1) inputDataLeft(C,P), (C,P,D) or (C,P,D,D); P=1 is admissible to allow multiply by const. left data   
01954      *      (2) inputDataRight(C,P,D,D) or (P,D,D);   
01955      *      (3) outputData(C,P,D,D)
01956      */
01957     // (1) inputDataLeft is (C,P), (C,P,D) or (C,P,D,D) and 1 <= D <= 3 is required  
01958     TEUCHOS_TEST_FOR_EXCEPTION( !requireRankRange(errmsg, inputDataLeft,  2,4), 
01959                         std::invalid_argument, errmsg);
01960     if(inputDataLeft.rank() > 2) {
01961       TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, inputDataLeft, 2,  1,3), 
01962                           std::invalid_argument, errmsg);
01963     }
01964     if(inputDataLeft.rank() == 4) {
01965       TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, inputDataLeft, 3,  1,3), 
01966                           std::invalid_argument, errmsg);
01967     }
01968     // (2) inputDataRight is (C,P,D,D) or (P,D,D) and 1 <= (D=dimension(rank-1),(rank-2)) <= 3 is required. 
01969     TEUCHOS_TEST_FOR_EXCEPTION( !requireRankRange(errmsg, inputDataRight,  3,4), 
01970                         std::invalid_argument, errmsg);
01971     TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, inputDataRight, inputDataRight.rank()-1,  1,3), 
01972                         std::invalid_argument, errmsg);
01973     TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, inputDataRight, inputDataRight.rank()-2,  1,3), 
01974                         std::invalid_argument, errmsg);
01975     // (3) outputData is (C,P,D,D)
01976     TEUCHOS_TEST_FOR_EXCEPTION( !requireRankRange(errmsg, outputData, 4, 4), std::invalid_argument, errmsg);      
01977     TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, outputData, 2,  1,3), 
01978                         std::invalid_argument, errmsg);
01979     TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, outputData, 3,  1,3), 
01980                         std::invalid_argument, errmsg);
01981     /*
01982      *   Dimension cross-checks:
01983      *      (1) inputDataLeft vs. inputDataRight
01984      *      (2) outputData    vs. inputDataLeft
01985      *      (3) outputData    vs. inputDataRight 
01986      *
01987      *   Cross-check (2): outputData(C,P,D,D) vs. inputDataLeft(C,P), (C,P,D) or (C,P,D,D):
01988      *   dimensions C, and D must match in all cases, dimension P must match only when non-constant
01989      *   data is specified (P>1). Do not check P dimensions with constant left data, i.e., when P=1 in
01990      *   inputDataLeft(C,1,...)
01991      */
01992     if(inputDataLeft.dimension(1) > 1){ // check P dimension if P>1 in inputDataLeft
01993       TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg,
01994                                                  outputData,     1,
01995                                                  inputDataLeft,  1),
01996                           std::invalid_argument, errmsg);    
01997     }
01998     if(inputDataLeft.rank() == 2){  // inputDataLeft(C,P): check C
01999       TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg,
02000                                                  outputData,    0,
02001                                                  inputDataLeft, 0),
02002                           std::invalid_argument, errmsg);    
02003     }
02004     if(inputDataLeft.rank() == 3){   // inputDataLeft(C,P,D): check C and D
02005       TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg,
02006                                                  outputData,    0,2,3,
02007                                                  inputDataLeft, 0,2,2),
02008                           std::invalid_argument, errmsg);    
02009     }
02010     if(inputDataLeft.rank() == 4){   // inputDataLeft(C,P,D,D): check C and D
02011       TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg,
02012                                                  outputData,    0,2,3,
02013                                                  inputDataLeft, 0,2,3),
02014                           std::invalid_argument, errmsg);    
02015     }
02016     /*
02017      *   Cross-checks (1,3):
02018      */
02019     if( inputDataRight.rank() == 4) {
02020       // Cross-check (1): inputDataLeft(C,P), (C,P,D), or (C,P,D,D) vs. inputDataRight(C,P,D,D):  dimensions  C, P, D must match
02021       if(inputDataLeft.dimension(1) > 1){ // check P dimension if P>1 in inputDataLeft
02022         TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg,
02023                                                    inputDataLeft,  1,
02024                                                    inputDataRight, 1),
02025                             std::invalid_argument, errmsg);    
02026       }      
02027       if(inputDataLeft.rank() == 2){  // inputDataLeft(C,P): check C
02028         TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, 
02029                                                    inputDataLeft,  0, 
02030                                                    inputDataRight, 0),
02031                             std::invalid_argument, errmsg);  
02032       }      
02033       if(inputDataLeft.rank() == 3){   // inputDataLeft(C,P,D): check C and D
02034         TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, 
02035                                                    inputDataLeft,  0,2,2, 
02036                                                    inputDataRight, 0,2,3),
02037                             std::invalid_argument, errmsg);  
02038       }      
02039       if(inputDataLeft.rank() == 4){   // inputDataLeft(C,P,D,D): check C and D
02040         TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, 
02041                                                    inputDataLeft,  0,2,3, 
02042                                                    inputDataRight, 0,2,3),
02043                             std::invalid_argument, errmsg);  
02044       }
02045       // Cross-check (3): outputData(C,P,D,D) vs. inputDataRight(C,P,D,D): all dimensions C, P, D must match
02046       TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, outputData, inputDataRight),
02047                           std::invalid_argument, errmsg);
02048     }
02049     else{
02050       // Cross-check (1): inputDataLeft(C,P), (C,P,D), or (C,P,D,D) vs. inputDataRight(P,D,D): dimensions  P, D must match
02051       if(inputDataLeft.dimension(1) > 1){ // check P if P>1 in inputData 
02052         TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg,
02053                                                    inputDataLeft,  1,
02054                                                    inputDataRight, 0),
02055                             std::invalid_argument, errmsg);    
02056       }
02057       if(inputDataLeft.rank() == 3){   // inputDataLeft(C,P,D): check D
02058         TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, 
02059                                                    inputDataLeft,  2,2, 
02060                                                    inputDataRight, 1,2),
02061                             std::invalid_argument, errmsg); 
02062       }
02063       if(inputDataLeft.rank() == 4){   // inputDataLeft(C,P,D,D): check D      
02064         TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, 
02065                                                    inputDataLeft,  2,3, 
02066                                                    inputDataRight, 1,2),
02067                             std::invalid_argument, errmsg); 
02068       }
02069       // Cross-check (3): outputData(C,P,D,D) vs. inputDataRight(P,D,D): dimensions P, D must match
02070       TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, 
02071                                                  outputData,     1,2,3, 
02072                                                  inputDataRight, 0,1,2),
02073                           std::invalid_argument, errmsg);
02074     }
02075 #endif
02076     int dataLeftRank   = inputDataLeft.rank();
02077     int numDataLeftPts = inputDataLeft.dimension(1);
02078     int dataRightRank  = inputDataRight.rank();    
02079     int numCells       = outputData.dimension(0);
02080     int numPoints      = outputData.dimension(1);
02081     int matDim         = outputData.dimension(2);
02082     
02083     /*********************************************************************************************
02084       *                              inputDataRight is (C,P,D,D)                                 *
02085       *********************************************************************************************/
02086     if(dataRightRank == 4){
02087       if(numDataLeftPts != 1){  // non-constant data
02088         
02089         switch(dataLeftRank){
02090           case 2:
02091             for(int cell = 0; cell < numCells; cell++) {
02092                 for(int point = 0; point < numPoints; point++) {
02093                   for( int row = 0; row < matDim; row++) {
02094                     for( int col = 0; col < matDim; col++) {
02095                       outputData(cell, point, row, col) = \
02096                       inputDataLeft(cell, point)*inputDataRight(cell, point, row, col);
02097                     }// Col-loop
02098                   } // Row-loop
02099                 } // P-loop
02100             }// C-loop
02101             break;
02102             
02103           case 3:
02104             for(int cell = 0; cell < numCells; cell++) {
02105                 for(int point = 0; point < numPoints; point++) {
02106                   for( int row = 0; row < matDim; row++) {
02107                     for( int col = 0; col < matDim; col++) {
02108                       outputData(cell, point, row, col) = \
02109                       inputDataLeft(cell, point, row)*inputDataRight(cell, point, row, col);
02110                     }// Col-loop
02111                   } // Row-loop
02112                 } // P-loop
02113             }// C-loop
02114             break;
02115             
02116           case 4:
02117             if ((transpose == 'n') || (transpose == 'N')) {
02118               for(int cell = 0; cell < numCells; cell++){
02119                   for(int point = 0; point < numPoints; point++){
02120                     for(int row = 0; row < matDim; row++){
02121                       for(int col = 0; col < matDim; col++){
02122                         outputData(cell, point, row, col) = 0.0;
02123                         for(int i = 0; i < matDim; i++){
02124                           outputData(cell, point, row, col) += \
02125                           inputDataLeft(cell, point, row, i)*inputDataRight(cell, point, i, col);
02126                         }// i
02127                       } // col
02128                     } //row
02129                   }// point
02130               }// cell
02131             } // no transpose
02132             else if ((transpose == 't') || (transpose == 'T')) {
02133               for(int cell = 0; cell < numCells; cell++){
02134                   for(int point = 0; point < numPoints; point++){
02135                     for(int row = 0; row < matDim; row++){
02136                       for(int col = 0; col < matDim; col++){
02137                         outputData(cell, point, row, col) = 0.0;
02138                         for(int i = 0; i < matDim; i++){
02139                           outputData(cell, point, row, col) += \
02140                           inputDataLeft(cell, point, i, row)*inputDataRight(cell, point, i, col);
02141                         }// i
02142                       } // col
02143                     } //row
02144                   }// point
02145               }// cell
02146             } //transpose
02147             else {
02148               TEUCHOS_TEST_FOR_EXCEPTION( !( (transpose == 'n') || (transpose == 'N') || (transpose == 't') || (transpose == 'T') ), std::invalid_argument,
02149                                   ">>> ERROR (ArrayTools::matmatProductDataData): The transpose flag must be 'n', 'N', 't' or 'T'.");
02150             }
02151             break;
02152             
02153           default:
02154             TEUCHOS_TEST_FOR_EXCEPTION( !( (dataLeftRank == 2) || (dataLeftRank == 3) || (dataLeftRank == 4) ), std::invalid_argument,
02155                                 ">>> ERROR (ArrayTools::matmatProductDataData): inputDataLeft rank 2, 3 or 4 required.")      
02156         } // switch inputData rank
02157       }
02158       else{  // constant data case
02159         switch(dataLeftRank){
02160           case 2:
02161             for(int cell = 0; cell < numCells; cell++) {
02162                 for(int point = 0; point < numPoints; point++) {
02163                   for( int row = 0; row < matDim; row++) {
02164                     for( int col = 0; col < matDim; col++) {
02165                       outputData(cell, point, row, col) = \
02166                       inputDataLeft(cell, 0)*inputDataRight(cell, point, row, col);
02167                     }// Col-loop
02168                   } // Row-loop
02169                 } // P-loop
02170             }// C-loop
02171             break;
02172             
02173           case 3:
02174             for(int cell = 0; cell < numCells; cell++) {
02175                 for(int point = 0; point < numPoints; point++) {
02176                   for( int row = 0; row < matDim; row++) {
02177                     for( int col = 0; col < matDim; col++) {
02178                       outputData(cell, point, row, col) = \
02179                       inputDataLeft(cell, 0, row)*inputDataRight(cell, point, row, col);
02180                     }// Col-loop
02181                   } // Row-loop
02182                 } // P-loop
02183             }// C-loop
02184             break;
02185             
02186           case 4:
02187             if ((transpose == 'n') || (transpose == 'N')) {
02188               for(int cell = 0; cell < numCells; cell++){
02189                   for(int point = 0; point < numPoints; point++){
02190                     for(int row = 0; row < matDim; row++){
02191                       for(int col = 0; col < matDim; col++){
02192                         outputData(cell, point, row, col) = 0.0;
02193                         for(int i = 0; i < matDim; i++){
02194                           outputData(cell, point, row, col) += \
02195                           inputDataLeft(cell, 0, row, i)*inputDataRight(cell, point, i, col);
02196                         }// i
02197                       } // col
02198                     } //row
02199                   }// point
02200               }// cell
02201             } // no transpose
02202             else if ((transpose == 't') || (transpose == 'T')) {
02203               for(int cell = 0; cell < numCells; cell++){
02204                   for(int point = 0; point < numPoints; point++){
02205                     for(int row = 0; row < matDim; row++){
02206                       for(int col = 0; col < matDim; col++){
02207                         outputData(cell, point, row, col) = 0.0;
02208                         for(int i = 0; i < matDim; i++){
02209                           outputData(cell, point, row, col) += \
02210                           inputDataLeft(cell, 0, i, row)*inputDataRight(cell, point, i, col);
02211                         }// i
02212                       } // col
02213                     } //row
02214                   }// point
02215               }// cell
02216             } //transpose
02217             else {
02218               TEUCHOS_TEST_FOR_EXCEPTION( !( (transpose == 'n') || (transpose == 'N') || (transpose == 't') || (transpose == 'T') ), std::invalid_argument,
02219                                   ">>> ERROR (ArrayTools::matmatProductDataData): The transpose flag must be 'n', 'N', 't' or 'T'.");
02220             }
02221             break;
02222             
02223           default:
02224             TEUCHOS_TEST_FOR_EXCEPTION( !( (dataLeftRank == 2) || (dataLeftRank == 3) || (dataLeftRank == 4) ), std::invalid_argument,
02225                                 ">>> ERROR (ArrayTools::matmatProductDataData): inputDataLeft rank 2, 3 or 4 required.")      
02226         } // switch inputDataLeft rank
02227       } // end constant data case
02228     } // inputDataRight rank 4
02229     /**********************************************************************************************
02230       *                              inputDataRight is (P,D,D)                                    *
02231       *********************************************************************************************/
02232     else if(dataRightRank == 3) {
02233       if(numDataLeftPts != 1){  // non-constant data
02234         
02235         switch(dataLeftRank){
02236           case 2:
02237             for(int cell = 0; cell < numCells; cell++) {
02238                 for(int point = 0; point < numPoints; point++) {
02239                   for( int row = 0; row < matDim; row++) {
02240                     for( int col = 0; col < matDim; col++) {
02241                       outputData(cell, point, row, col) = \
02242                       inputDataLeft(cell, point)*inputDataRight(point, row, col);
02243                     }// Col-loop
02244                   } // Row-loop
02245                 } // P-loop
02246             }// C-loop
02247             break;
02248             
02249           case 3:
02250             for(int cell = 0; cell < numCells; cell++) {
02251                 for(int point = 0; point < numPoints; point++) {
02252                   for( int row = 0; row < matDim; row++) {
02253                     for( int col = 0; col < matDim; col++) {
02254                       outputData(cell, point, row, col) = \
02255                       inputDataLeft(cell, point, row)*inputDataRight(point, row, col);
02256                     }// Col-loop
02257                   } // Row-loop
02258                 } // P-loop
02259             }// C-loop
02260             break;
02261             
02262           case 4:
02263             if ((transpose == 'n') || (transpose == 'N')) {
02264               for(int cell = 0; cell < numCells; cell++){
02265                   for(int point = 0; point < numPoints; point++){
02266                     for(int row = 0; row < matDim; row++){
02267                       for(int col = 0; col < matDim; col++){
02268                         outputData(cell, point, row, col) = 0.0;
02269                         for(int i = 0; i < matDim; i++){
02270                           outputData(cell, point, row, col) += \
02271                           inputDataLeft(cell, point, row, i)*inputDataRight(point, i, col);
02272                         }// i
02273                       } // col
02274                     } //row
02275                   }// point
02276               }// cell
02277             } // no transpose
02278             else if ((transpose == 't') || (transpose == 'T')) {
02279               for(int cell = 0; cell < numCells; cell++){
02280                   for(int point = 0; point < numPoints; point++){
02281                     for(int row = 0; row < matDim; row++){
02282                       for(int col = 0; col < matDim; col++){
02283                         outputData(cell, point, row, col) = 0.0;
02284                         for(int i = 0; i < matDim; i++){
02285                           outputData(cell, point, row, col) += \
02286                           inputDataLeft(cell, point, i, row)*inputDataRight(point, i, col);
02287                         }// i
02288                       } // col
02289                     } //row
02290                   }// point
02291               }// cell
02292             } //transpose
02293             else {
02294               TEUCHOS_TEST_FOR_EXCEPTION( !( (transpose == 'n') || (transpose == 'N') || (transpose == 't') || (transpose == 'T') ), std::invalid_argument,
02295                                   ">>> ERROR (ArrayTools::matmatProductDataData): The transpose flag must be 'n', 'N', 't' or 'T'.");
02296             }
02297             break;
02298             
02299           default:
02300             TEUCHOS_TEST_FOR_EXCEPTION( !( (dataLeftRank == 2) || (dataLeftRank == 3) || (dataLeftRank == 4) ), std::invalid_argument,
02301                                 ">>> ERROR (ArrayTools::matmatProductDataData): inputDataLeft rank 2, 3 or 4 required.")      
02302         } // switch inputDataLeft rank
02303       }
02304       else{  // constant data case
02305         switch(dataLeftRank){
02306           case 2:
02307             for(int cell = 0; cell < numCells; cell++) {
02308                 for(int point = 0; point < numPoints; point++) {
02309                   for( int row = 0; row < matDim; row++) {
02310                     for( int col = 0; col < matDim; col++) {
02311                       outputData(cell, point, row, col) = \
02312                       inputDataLeft(cell, 0)*inputDataRight(point, row, col);
02313                     }// Col-loop
02314                   } // Row-loop
02315                 } // P-loop
02316             }// C-loop
02317             break;
02318             
02319           case 3:
02320             for(int cell = 0; cell < numCells; cell++) {
02321                 for(int point = 0; point < numPoints; point++) {
02322                   for( int row = 0; row < matDim; row++) {
02323                     for( int col = 0; col < matDim; col++) {
02324                       outputData(cell, point, row, col) = \
02325                       inputDataLeft(cell, 0, row)*inputDataRight(point, row, col);
02326                     }// Col-loop
02327                   } // Row-loop
02328                 } // P-loop
02329             }// C-loop
02330             break;
02331             
02332           case 4:
02333             if ((transpose == 'n') || (transpose == 'N')) {
02334               for(int cell = 0; cell < numCells; cell++){
02335                   for(int point = 0; point < numPoints; point++){
02336                     for(int row = 0; row < matDim; row++){
02337                       for(int col = 0; col < matDim; col++){
02338                         outputData(cell, point, row, col) = 0.0;
02339                         for(int i = 0; i < matDim; i++){
02340                           outputData(cell, point, row, col) += \
02341                           inputDataLeft(cell, 0, row, i)*inputDataRight(point, i, col);
02342                         }// i
02343                       } // col
02344                     } //row
02345                   }// point
02346               }// cell
02347             } // no transpose
02348             else if ((transpose == 't') || (transpose == 'T')) {
02349               for(int cell = 0; cell < numCells; cell++){
02350                   for(int point = 0; point < numPoints; point++){
02351                     for(int row = 0; row < matDim; row++){
02352                       for(int col = 0; col < matDim; col++){
02353                         outputData(cell, point, row, col) = 0.0;
02354                         for(int i = 0; i < matDim; i++){
02355                           outputData(cell, point, row, col) += \
02356                           inputDataLeft(cell, 0, i, row)*inputDataRight(point, i, col);
02357                         }// i
02358                       } // col
02359                     } //row
02360                   }// point
02361               }// cell
02362             } //transpose
02363             else {
02364               TEUCHOS_TEST_FOR_EXCEPTION( !( (transpose == 'n') || (transpose == 'N') || (transpose == 't') || (transpose == 'T') ), std::invalid_argument,
02365                                   ">>> ERROR (ArrayTools::matmatProductDataData): The transpose flag must be 'n', 'N', 't' or 'T'.");
02366             }
02367             break;
02368             
02369           default:
02370             TEUCHOS_TEST_FOR_EXCEPTION( !( (dataLeftRank == 2) || (dataLeftRank == 3) || (dataLeftRank == 4) ), std::invalid_argument,
02371                                 ">>> ERROR (ArrayTools::matmatProductDataData): inputDataLeft rank 2, 3 or 4 required.")      
02372         } // switch inputDataLeft rank
02373       } // end constant data case
02374     } // inputDataRight rank 3
02375     else {
02376       TEUCHOS_TEST_FOR_EXCEPTION( true, std::invalid_argument, 
02377                           ">>> ERROR (ArrayTools::matmatProductDataData): inputDataRight rank 3 or 4 required.")      
02378     }// rank error
02379   }
02380   
02381   
02382   
02383   
02384   
02385 } // end namespace Intrepid
02386 
02387 
02388 
02389 
02390 
02391 
02392 
02393 
02394 
02395 
02396 
02397 
02398 
02399 
02400 
02401 
02402 
02403 
02404 
02405 
02406