Intrepid
http://trilinos.sandia.gov/packages/docs/r10.10/packages/intrepid/test/Shared/ArrayTools/test_04.cpp
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 #include "Intrepid_ArrayTools.hpp"
00050 #include "Intrepid_FieldContainer.hpp"
00051 #include "Intrepid_RealSpaceTools.hpp"
00052 #include "Teuchos_oblackholestream.hpp"
00053 #include "Teuchos_RCP.hpp"
00054 #include "Teuchos_ScalarTraits.hpp"
00055 
00056 using namespace std;
00057 using namespace Intrepid;
00058 
00059 #define INTREPID_TEST_COMMAND( S )                                                                                  \
00060 {                                                                                                                   \
00061   try {                                                                                                             \
00062     S ;                                                                                                             \
00063   }                                                                                                                 \
00064   catch (std::logic_error err) {                                                                                    \
00065       *outStream << "Expected Error ----------------------------------------------------------------\n";            \
00066       *outStream << err.what() << '\n';                                                                             \
00067       *outStream << "-------------------------------------------------------------------------------" << "\n\n";    \
00068   };                                                                                                                \
00069 }
00070 
00071 
00072 int main(int argc, char *argv[]) {
00073 
00074   // This little trick lets us print to std::cout only if
00075   // a (dummy) command-line argument is provided.
00076   int iprint     = argc - 1;
00077   Teuchos::RCP<std::ostream> outStream;
00078   Teuchos::oblackholestream bhs; // outputs nothing
00079   if (iprint > 0)
00080     outStream = Teuchos::rcp(&std::cout, false);
00081   else
00082     outStream = Teuchos::rcp(&bhs, false);
00083 
00084   // Save the format state of the original std::cout.
00085   Teuchos::oblackholestream oldFormatState;
00086   oldFormatState.copyfmt(std::cout);
00087 
00088   *outStream \
00089   << "===============================================================================\n" \
00090   << "|                                                                             |\n" \
00091   << "|                       Unit Test (ArrayTools)                                |\n" \
00092   << "|                                                                             |\n" \
00093   << "|     1) Array operations: multiplication, contractions                       |\n" \
00094   << "|                                                                             |\n" \
00095   << "|  Questions? Contact  Pavel Bochev (pbboche@sandia.gov) or                   |\n" \
00096   << "|                      Denis Ridzal (dridzal@sandia.gov).                     |\n" \
00097   << "|                                                                             |\n" \
00098   << "|  Intrepid's website: http://trilinos.sandia.gov/packages/intrepid           |\n" \
00099   << "|  Trilinos website:   http://trilinos.sandia.gov                             |\n" \
00100   << "|                                                                             |\n" \
00101   << "===============================================================================\n";
00102 
00103 
00104   int errorFlag = 0;
00105 #ifdef HAVE_INTREPID_DEBUG
00106   int beginThrowNumber = Teuchos::TestForException_getThrowNumber();
00107   int endThrowNumber = beginThrowNumber + 39 + 18 + 28 + 26 + 34 + 37 + 46 + 45;
00108 #endif
00109 
00110   typedef ArrayTools art; 
00111   typedef RealSpaceTools<double> rst; 
00112   
00113   
00114 #ifdef HAVE_INTREPID_DEBUG
00115   ArrayTools atools;
00116   /************************************************************************************************
00117     *                                                                                             *
00118     *  Exception tests: should only run when compiled in DEBUG mode                               *
00119     *                                                                                             *
00120     ***********************************************************************************************/
00121   int C  = 12,  C1 = 10;
00122   int P  = 8,   P1 = 16;
00123   int F  = 4,   F1 = 8; 
00124   int D1 = 1,   D2 = 2,   D3 = 3,   D4 = 4;
00125   
00126   FieldContainer<double> fc_P          (P);
00127   FieldContainer<double> fc_P_D1       (P, D1);
00128   FieldContainer<double> fc_P_D2       (P, D2);
00129   FieldContainer<double> fc_P_D3       (P, D3);
00130   FieldContainer<double> fc_P_D2_D2    (P, D2, D2);
00131   FieldContainer<double> fc_P1_D3      (P1, D3);
00132   FieldContainer<double> fc_P1_D2_D2   (P1, D2, D2);
00133   FieldContainer<double> fc_P1_D3_D3   (P1, D3, D3);
00134   
00135   FieldContainer<double> fc_C              (C);
00136   FieldContainer<double> fc_C1_P           (C1,P);
00137   FieldContainer<double> fc_C_P1           (C, P1);
00138   FieldContainer<double> fc_C_P            (C, P);
00139   FieldContainer<double> fc_C_P_D1         (C, P, D1);
00140   FieldContainer<double> fc_C_P_D2         (C, P, D2);
00141   FieldContainer<double> fc_C_P_D3         (C, P, D3);
00142   FieldContainer<double> fc_C_P_D4         (C, P, D4);
00143   FieldContainer<double> fc_C1_P_D2        (C1,P, D2);
00144   FieldContainer<double> fc_C1_P_D3        (C1,P, D3);
00145   FieldContainer<double> fc_C_P1_D1        (C, P1,D1);
00146   FieldContainer<double> fc_C_P1_D2        (C, P1,D2);
00147   FieldContainer<double> fc_C_P1_D3        (C, P1,D3);
00148   FieldContainer<double> fc_C_P_D1_D1      (C, P, D1, D1);
00149   FieldContainer<double> fc_C_P_D1_D2      (C, P, D1, D2);
00150   FieldContainer<double> fc_C_P_D1_D3      (C, P, D1, D3);
00151   FieldContainer<double> fc_C_P_D2_D2      (C, P, D2, D2);
00152   FieldContainer<double> fc_C_P_D3_D1      (C, P, D3, D1);
00153   FieldContainer<double> fc_C_P_D3_D2      (C, P, D3, D2);
00154   FieldContainer<double> fc_C_P_D2_D3      (C, P, D2, D3);
00155   FieldContainer<double> fc_C_P_D3_D3      (C, P, D3, D3);
00156   FieldContainer<double> fc_C1_P_D3_D3     (C1,P, D3, D3);
00157   FieldContainer<double> fc_C1_P_D2_D2     (C1,P, D2, D2);
00158   FieldContainer<double> fc_C_P1_D2_D2     (C, P1,D2, D2);
00159   FieldContainer<double> fc_C_P1_D3_D3     (C, P1,D3, D3);
00160   FieldContainer<double> fc_C_P_D3_D3_D3   (C, P, D3, D3, D3);
00161   
00162   FieldContainer<double> fc_F_P            (F, P);
00163   FieldContainer<double> fc_F_P_D1         (F, P, D1);
00164   FieldContainer<double> fc_F_P_D2         (F, P, D2);
00165   FieldContainer<double> fc_F_P1_D2        (F, P1, D2);
00166   FieldContainer<double> fc_F_P_D3         (F, P, D3);
00167   FieldContainer<double> fc_F_P_D3_D3      (F, P, D3, D3);
00168   FieldContainer<double> fc_F1_P_D2        (F1,P, D2);
00169   FieldContainer<double> fc_F1_P_D3        (F1,P, D3);
00170   FieldContainer<double> fc_F1_P_D3_D3     (F1,P, D3, D3);
00171   FieldContainer<double> fc_F_P1_D3        (F, P1,D3);
00172   FieldContainer<double> fc_F_P1_D3_D3     (F, P1,D3, D3);
00173   FieldContainer<double> fc_F_P_D2_D2      (F, P, D2, D2);
00174   FieldContainer<double> fc_F_P1_D2_D2     (F, P1,D2, D2);
00175   FieldContainer<double> fc_C_F_P          (C, F, P);
00176   FieldContainer<double> fc_C1_F_P         (C1, F, P);
00177   FieldContainer<double> fc_C_F1_P         (C, F1,P);
00178   FieldContainer<double> fc_C_F_P1         (C, F, P1);
00179   FieldContainer<double> fc_C_F_P_D1       (C, F, P, D1);
00180   FieldContainer<double> fc_C_F_P_D2       (C, F, P, D2);
00181   FieldContainer<double> fc_C_F_P_D3       (C, F, P, D3);
00182   FieldContainer<double> fc_C1_F_P_D2      (C1, F, P,D2);
00183   FieldContainer<double> fc_C1_F_P_D3      (C1, F, P,D3);
00184   FieldContainer<double> fc_C_F1_P_D2      (C, F1,P, D2);
00185   FieldContainer<double> fc_C_F1_P_D3      (C, F1,P, D3);
00186   FieldContainer<double> fc_C_F1_P_D3_D3   (C, F1,P, D3, D3);
00187   FieldContainer<double> fc_C_F_P1_D2      (C, F, P1,D2);
00188   FieldContainer<double> fc_C_F_P1_D3      (C, F, P1,D3);
00189   FieldContainer<double> fc_C_F_P_D1_D1    (C, F, P, D1, D1);
00190   FieldContainer<double> fc_C_F_P_D2_D2    (C, F, P, D2, D2);
00191   FieldContainer<double> fc_C_F_P_D1_D3    (C, F, P, D1, D3);
00192   FieldContainer<double> fc_C_F_P_D2_D3    (C, F, P, D2, D3);
00193   FieldContainer<double> fc_C_F_P_D3_D1    (C, F, P, D3, D1);
00194   FieldContainer<double> fc_C_F_P_D3_D2    (C, F, P, D3, D2);
00195   FieldContainer<double> fc_C_F_P_D3_D3    (C, F, P, D3, D3);
00196   FieldContainer<double> fc_C_F_P1_D2_D2   (C, F, P1,D2, D2);
00197   FieldContainer<double> fc_C_F_P1_D3_D3   (C, F, P1,D3, D3);
00198   FieldContainer<double> fc_C1_F_P_D2_D2   (C1,F, P, D2, D2);
00199   FieldContainer<double> fc_C1_F_P1_D2_D2  (C1,F, P1,D2, D2);
00200   FieldContainer<double> fc_C1_F_P_D3_D3   (C1,F, P, D3, D3);
00201   
00202   Teuchos::Array<int> dimensions(6);
00203   dimensions[0] = C; 
00204   dimensions[1] = F; 
00205   dimensions[2] = P;
00206   dimensions[3] = D3;
00207   dimensions[4] = D3;
00208   dimensions[5] = D3;
00209   FieldContainer<double> fc_C_F_P_D3_D3_D3 (dimensions);
00210   
00211   
00212   *outStream \
00213     << "\n"
00214     << "===============================================================================\n"\
00215     << "| TEST 1: crossProductDataField exceptions                                    |\n"\
00216     << "===============================================================================\n";
00217   
00218   try{
00219     // 39 exceptions
00220     // Test rank and D dimension of inputData
00221     INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C_F_P_D3, fc_C_P,       fc_C_F_P_D3) );
00222     INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C_F_P_D3, fc_C_P_D2_D2, fc_C_F_P_D3) );
00223     INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C_F_P_D3, fc_C_P_D1,    fc_C_F_P_D3) );
00224         
00225     // Test rank and D dimension of inputFields
00226     INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C_F_P_D3,  fc_C_P_D3,  fc_F_P) );
00227     INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C_F_P_D3,  fc_C_P_D3,  fc_C_F_P_D3_D3) );
00228     INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C_F_P_D3,  fc_C_P_D3,  fc_C_F_P_D1) );
00229     INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C_F_P_D3,  fc_C_P_D3,  fc_F_P_D1) );
00230 
00231     // Test rank of outputFields
00232     INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C_F_P_D3,  fc_C_P_D2,  fc_C_F_P_D2) );
00233     INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C_F_P,     fc_C_P_D3,  fc_C_F_P_D3) );
00234 
00235     // Dimension cross-check: (1) inputData    vs. inputFields
00236     INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C_F_P_D3,  fc_C_P_D3, fc_C1_F_P_D3) );
00237     INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C_F_P_D3,  fc_C_P_D3, fc_C_F_P1_D3) );
00238     INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C_F_P_D3,  fc_C_P_D3, fc_C_F_P_D2) );
00239     INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C_F_P,     fc_C_P_D2, fc_C1_F_P_D2) );
00240     INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C_F_P,     fc_C_P_D2, fc_C_F_P1_D2) );
00241     INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C_F_P,     fc_C_P_D2, fc_C_F_P_D3) );
00242     INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C_F_P_D3,  fc_C_P_D3, fc_F_P1_D3) );
00243     INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C_F_P_D3,  fc_C_P_D3, fc_F_P_D2) );
00244     INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C_F_P,     fc_C_P_D2, fc_F_P1_D2) );
00245     INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C_F_P,     fc_C_P_D2, fc_F_P_D3) );
00246     
00247     // Dimension cross-check: (2) outputFields vs. inputData
00248     INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C1_F_P,  fc_C_P_D2, fc_C_F_P_D2) );
00249     INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C_F_P1,  fc_C_P_D2, fc_C_F_P_D2) );
00250     INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C1_F_P,  fc_C_P_D2, fc_F_P_D2) );
00251     INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C_F_P1,  fc_C_P_D2, fc_F_P_D2) );
00252     INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C1_F_P_D3, fc_C_P_D3, fc_C_F_P_D3) );
00253     INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C_F_P1_D3, fc_C_P_D3, fc_C_F_P_D3) );
00254     INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C_F_P_D2,  fc_C_P_D3, fc_C_F_P_D3) );
00255     INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C1_F_P_D3, fc_C_P_D3, fc_F_P_D3) );
00256     INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C_F_P1_D3, fc_C_P_D3, fc_F_P_D3) );
00257     INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C_F_P_D2,  fc_C_P_D3, fc_F_P_D3) );
00258     
00259      // Dimension cross-check: (3) outputFields vs. inputFields
00260     INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C_F_P,  fc_C1_P_D2, fc_C1_F_P_D2) );
00261     INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C_F_P,  fc_C_P1_D2, fc_C_F_P1_D2) );
00262     INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C_F_P,  fc_C_P_D2, fc_C_F1_P_D2) );
00263     INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C_F_P,  fc_C_P1_D2, fc_F_P1_D2) );
00264     INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C_F_P,  fc_C_P_D2, fc_F1_P_D2) );
00265     INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C_F_P_D3, fc_C1_P_D3,  fc_C1_F_P_D3) );
00266     INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C_F_P_D3, fc_C_P1_D3,  fc_C_F_P1_D3) );
00267     INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C_F_P_D3, fc_C_P_D3,   fc_C_F1_P_D3) );
00268     INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C_F_P_D3, fc_C_P1_D3,  fc_F_P1_D3) );
00269     INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C_F_P_D3, fc_C_P_D3,   fc_F1_P_D3) );
00270      
00271      *outStream \
00272      << "\n"
00273      << "===============================================================================\n"\
00274      << "| TEST 2: crossProductDataData exceptions                                     |\n"\
00275      << "===============================================================================\n";
00276      
00277      // 18 exceptions 
00278      // inputDataL is (C, P, D) and 2 <= D <= 3 is required  
00279      INTREPID_TEST_COMMAND(atools.crossProductDataData<double>(fc_C_P_D3, fc_C_P,       fc_C_P_D3) );
00280      INTREPID_TEST_COMMAND(atools.crossProductDataData<double>(fc_C_P_D3, fc_C_P_D2_D2, fc_C_P_D3) );
00281      INTREPID_TEST_COMMAND(atools.crossProductDataData<double>(fc_C_P_D3, fc_C_P_D1,    fc_C_P_D3) );
00282      
00283      // inputDataRight is (C, P, D) or (P, D) and 2 <= (D=dimension(rank - 1)) <= 3 is required. 
00284      INTREPID_TEST_COMMAND(atools.crossProductDataData<double>(fc_C_P_D3, fc_C_P_D3,    fc_C) );
00285      INTREPID_TEST_COMMAND(atools.crossProductDataData<double>(fc_C_P_D3, fc_C_P_D3,    fc_C_P_D2_D2) );
00286      INTREPID_TEST_COMMAND(atools.crossProductDataData<double>(fc_C_P_D3, fc_C_P_D3,    fc_C_P_D1) );
00287 
00288      // outputData is (C,P,D) in 3D and (C,P) in 2D => rank = inputDataLeft.dimension(2)
00289      INTREPID_TEST_COMMAND(atools.crossProductDataData<double>(fc_C_P_D3, fc_C_P_D2,    fc_C_P_D2) ); 
00290      INTREPID_TEST_COMMAND(atools.crossProductDataData<double>(fc_C_P,    fc_C_P_D3,    fc_C_P_D3) );
00291    
00292      // Dimension cross-check (1): 
00293      // inputDataLeft(C,P,D) vs. inputDataRight(C,P,D): C, P, D must match 
00294      INTREPID_TEST_COMMAND(atools.crossProductDataData<double>(fc_C_P_D3, fc_C1_P_D3,   fc_C_P_D3) );
00295      INTREPID_TEST_COMMAND(atools.crossProductDataData<double>(fc_C_P_D3, fc_C_P1_D3,   fc_C_P_D3) );
00296      INTREPID_TEST_COMMAND(atools.crossProductDataData<double>(fc_C_P_D3, fc_C_P_D3,    fc_C_P_D2) );
00297      // inputDataLeft(C, P,D) vs. inputDataRight(P,D):  P, D must match
00298      INTREPID_TEST_COMMAND(atools.crossProductDataData<double>(fc_C_P_D3, fc_C_P1_D3,   fc_P_D3) );
00299      INTREPID_TEST_COMMAND(atools.crossProductDataData<double>(fc_C_P_D3, fc_C_P_D3,    fc_P_D2) );
00300      
00301      // Dimension cross-check (2):
00302      // in 2D: outputData(C,P) vs. inputDataLeft(C,P,D): dimensions C, P must match
00303      INTREPID_TEST_COMMAND(atools.crossProductDataData<double>(fc_C1_P,    fc_C_P_D2,   fc_P_D2) );
00304      INTREPID_TEST_COMMAND(atools.crossProductDataData<double>(fc_C_P1,    fc_C_P_D2,   fc_P_D2) ); 
00305      // in 3D: outputData(C,P,D) vs. inputDataLeft(C,P,D): all dimensions C, P, D must match
00306      INTREPID_TEST_COMMAND(atools.crossProductDataData<double>(fc_C1_P_D3, fc_C_P_D3,   fc_P_D3) );
00307      INTREPID_TEST_COMMAND(atools.crossProductDataData<double>(fc_C_P1_D3, fc_C_P_D3,   fc_P_D3) );
00308      INTREPID_TEST_COMMAND(atools.crossProductDataData<double>(fc_C_P_D2,  fc_C_P_D3,   fc_P_D3) );
00309       
00310     *outStream \
00311       << "\n"
00312       << "===============================================================================\n"\
00313       << "| TEST 3: outerProductDataField exceptions                                    |\n"\
00314       << "===============================================================================\n";
00315     // 28 exceptions
00316     // Test rank and D dimension: inputData(C, P, D) and 2 <= D <= 3 is required  
00317     INTREPID_TEST_COMMAND(atools.outerProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P,       fc_C_F_P_D3) );
00318     INTREPID_TEST_COMMAND(atools.outerProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D2_D2, fc_C_F_P_D3) );
00319     INTREPID_TEST_COMMAND(atools.outerProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D1,    fc_C_F_P_D3) );
00320     
00321     // Test rank and D dimension: inputFields(C,F,P,D)/(F,P,D) and 2 <= (D=dimension(rank - 1)) <= 3 is required. 
00322     INTREPID_TEST_COMMAND(atools.outerProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3,    fc_F_P) );
00323     INTREPID_TEST_COMMAND(atools.outerProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3,    fc_C_F_P_D3_D3) );
00324     INTREPID_TEST_COMMAND(atools.outerProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3,    fc_C_F_P_D1) );
00325     
00326     //  Test rank and D dimension: outputFields(C,F,P,D,D)
00327     INTREPID_TEST_COMMAND(atools.outerProductDataField<double>(fc_C_F_P,         fc_C_P_D3,    fc_C_F_P_D3) );
00328     INTREPID_TEST_COMMAND(atools.outerProductDataField<double>(fc_C_F_P_D3,      fc_C_P_D3,    fc_C_F_P_D3) );
00329     INTREPID_TEST_COMMAND(atools.outerProductDataField<double>(fc_C_F_P_D3_D3_D3,fc_C_P_D3,    fc_C_F_P_D3) );
00330     INTREPID_TEST_COMMAND(atools.outerProductDataField<double>(fc_C_F_P_D1_D3,   fc_C_P_D3,    fc_C_F_P_D3) );
00331     INTREPID_TEST_COMMAND(atools.outerProductDataField<double>(fc_C_F_P_D3_D1,   fc_C_P_D3,    fc_C_F_P_D3) );
00332     INTREPID_TEST_COMMAND(atools.outerProductDataField<double>(fc_C_F_P_D1_D1,   fc_C_P_D3,    fc_C_F_P_D3) );
00333     
00334     // Cross-check (2): outputFields(C,F,P,D,D) vs. inputData(C,P,D): dimensions C, P, D must match 
00335     INTREPID_TEST_COMMAND(atools.outerProductDataField<double>(fc_C_F_P_D3_D3, fc_C1_P_D3,   fc_C_F_P_D3) );
00336     INTREPID_TEST_COMMAND(atools.outerProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P1_D3,   fc_C_F_P_D3) );
00337     INTREPID_TEST_COMMAND(atools.outerProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D2,    fc_C_F_P_D2) );
00338     
00339     // Cross-check (1): inputData(C,P,D) vs. inputFields(C,F,P,D):  dimensions  C, P, D must match
00340     INTREPID_TEST_COMMAND(atools.outerProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3,   fc_C1_F_P_D3) );
00341     INTREPID_TEST_COMMAND(atools.outerProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3,   fc_C_F_P1_D3) );
00342     INTREPID_TEST_COMMAND(atools.outerProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3,   fc_C_F_P_D2) );
00343     // Cross-check (1): inputData(C,P,D) vs. inputFields(F,P,D): dimensions  P, D must match
00344     INTREPID_TEST_COMMAND(atools.outerProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3,   fc_F_P1_D3) );
00345     INTREPID_TEST_COMMAND(atools.outerProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3,   fc_F_P_D2) );
00346     
00347     // Cross-check (3): outputFields(C,F,P,D,D) vs. inputFields(C,F,P,D): dimensions C, F, P, D must match
00348     INTREPID_TEST_COMMAND(atools.outerProductDataField<double>(fc_C_F_P_D3_D3, fc_C1_P_D3,   fc_C1_F_P_D3) );
00349     INTREPID_TEST_COMMAND(atools.outerProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3,    fc_C_F1_P_D3) );
00350     INTREPID_TEST_COMMAND(atools.outerProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P1_D3,   fc_C_F_P1_D3) );
00351     INTREPID_TEST_COMMAND(atools.outerProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D2,    fc_C_F_P_D2) );
00352     INTREPID_TEST_COMMAND(atools.outerProductDataField<double>(fc_C_F_P_D2_D3, fc_C_P_D2,    fc_C_F_P_D2) );
00353     // Cross-check (3): outputFields(C,F,P,D,D) vs. inputFields(F,P,D): dimensions F, P, D must match
00354     INTREPID_TEST_COMMAND(atools.outerProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3,    fc_F1_P_D3) );
00355     INTREPID_TEST_COMMAND(atools.outerProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P1_D3,   fc_F_P1_D3) );
00356     INTREPID_TEST_COMMAND(atools.outerProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D2,    fc_F_P_D2) );  
00357    *outStream \
00358    << "\n"
00359    << "===============================================================================\n"\
00360    << "| TEST 4: outerProductDataData exceptions                                     |\n"\
00361    << "===============================================================================\n";
00362    // 26 exceptions
00363    // (1) inputDataLeft is (C, P, D) and 2 <= D <= 3 is required  
00364    INTREPID_TEST_COMMAND(atools.outerProductDataData<double>(fc_C_P_D3_D3, fc_C_P,         fc_C_P_D3) );
00365    INTREPID_TEST_COMMAND(atools.outerProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D3_D3,   fc_C_P_D3) );
00366    INTREPID_TEST_COMMAND(atools.outerProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D1,      fc_C_P_D3) );
00367 
00368    // (2) inputDataRight is (C, P, D) or (P, D) and 2 <= (D=dimension(rank - 1)) <= 3 is required. 
00369    INTREPID_TEST_COMMAND(atools.outerProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D3,      fc_C) );
00370    INTREPID_TEST_COMMAND(atools.outerProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D3,      fc_C_P_D3_D3) );
00371    INTREPID_TEST_COMMAND(atools.outerProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D3,      fc_C_P_D1) );
00372    INTREPID_TEST_COMMAND(atools.outerProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D3,      fc_P_D1) );
00373 
00374    // (3) outputData is (C,P,D,D)
00375    INTREPID_TEST_COMMAND(atools.outerProductDataData<double>(fc_C_P_D3,      fc_C_P_D3,      fc_C_P_D3) );
00376    INTREPID_TEST_COMMAND(atools.outerProductDataData<double>(fc_C_P_D3_D3_D3,fc_C_P_D3,      fc_C_P_D3) );
00377    INTREPID_TEST_COMMAND(atools.outerProductDataData<double>(fc_C_P_D3_D1,   fc_C_P_D3,      fc_C_P_D3) );
00378    INTREPID_TEST_COMMAND(atools.outerProductDataData<double>(fc_C_P_D1_D2,   fc_C_P_D3,      fc_C_P_D3) );
00379 
00380    // Cross-check (2): outputData(C,P,D,D) vs. inputDataLeft(C,P,D): dimensions C, P, D must match
00381    INTREPID_TEST_COMMAND(atools.outerProductDataData<double>(fc_C_P_D3_D3, fc_C1_P_D3,       fc_P_D3) );
00382    INTREPID_TEST_COMMAND(atools.outerProductDataData<double>(fc_C_P_D3_D3, fc_C_P1_D3,       fc_P_D3) );
00383    INTREPID_TEST_COMMAND(atools.outerProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D2,        fc_P_D3) );
00384    INTREPID_TEST_COMMAND(atools.outerProductDataData<double>(fc_C_P_D2_D3, fc_C_P_D2,        fc_P_D3) );
00385    INTREPID_TEST_COMMAND(atools.outerProductDataData<double>(fc_C_P_D3_D2, fc_C_P_D2,        fc_P_D3) );
00386 
00387    // Cross-check (1): inputDataLeft(C,P,D) vs. inputDataRight(C,P,D):  all dimensions  C, P, D must match
00388    INTREPID_TEST_COMMAND(atools.outerProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D3,       fc_C1_P_D3) );
00389    INTREPID_TEST_COMMAND(atools.outerProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D3,       fc_C_P1_D3) );
00390    INTREPID_TEST_COMMAND(atools.outerProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D3,       fc_C_P_D2) );
00391    // Cross-check (1): inputDataLeft(C,P,D) vs. inputDataRight(P,D): dimensions  P, D must match
00392    INTREPID_TEST_COMMAND(atools.outerProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D3,       fc_P1_D3) );
00393    INTREPID_TEST_COMMAND(atools.outerProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D3,       fc_P_D2) );
00394 
00395    // Cross-check (3): outputData(C,P,D,D) vs. inputDataRight(C,P,D): dimensions C, P, D must match
00396    INTREPID_TEST_COMMAND(atools.outerProductDataData<double>(fc_C_P_D3_D3, fc_C1_P_D3,       fc_C1_P_D3) );
00397    INTREPID_TEST_COMMAND(atools.outerProductDataData<double>(fc_C_P_D3_D3, fc_C_P1_D3,       fc_C_P1_D3) );
00398    INTREPID_TEST_COMMAND(atools.outerProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D2,        fc_C_P_D2) );
00399    // Cross-check (3): outputData(C,P,D,D) vs. inputDataRight(P,D): dimensions P, D must match
00400    INTREPID_TEST_COMMAND(atools.outerProductDataData<double>(fc_C_P_D3_D3, fc_C_P1_D3,      fc_P1_D3) );
00401    INTREPID_TEST_COMMAND(atools.outerProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D2,       fc_P_D2) );
00402    
00403    *outStream \
00404      << "\n"
00405      << "===============================================================================\n"\
00406      << "| TEST 5: matvecProductDataField exceptions                                   |\n"\
00407      << "===============================================================================\n";
00408    // 34 exceptions
00409    // (1) inputData is (C,P), (C,P,D) or (C, P, D, D) and 1 <= D <= 3 is required  
00410    INTREPID_TEST_COMMAND(atools.matvecProductDataField<double>(fc_C_F_P_D3, fc_C,              fc_C_F_P_D3) );
00411    INTREPID_TEST_COMMAND(atools.matvecProductDataField<double>(fc_C_F_P_D3, fc_C_P_D4,         fc_C_F_P_D3) );
00412    INTREPID_TEST_COMMAND(atools.matvecProductDataField<double>(fc_C_F_P_D3, fc_C_P_D3_D3_D3,   fc_C_F_P_D3) );
00413    INTREPID_TEST_COMMAND(atools.matvecProductDataField<double>(fc_C_F_P_D3, fc_C_P_D3_D1,      fc_C_F_P_D3) );   
00414    INTREPID_TEST_COMMAND(atools.matvecProductDataField<double>(fc_C_F_P_D3, fc_C_P_D1_D3,      fc_C_F_P_D3) ); 
00415    
00416    // (2) inputFields is (C, F, P, D) or (F, P, D) and 1 <= (D=dimension(rank - 1)) <= 3 is required. 
00417    INTREPID_TEST_COMMAND(atools.matvecProductDataField<double>(fc_C_F_P_D3, fc_C_P_D3_D3,    fc_P_D3) );
00418    INTREPID_TEST_COMMAND(atools.matvecProductDataField<double>(fc_C_F_P_D3, fc_C_P_D3_D3,    fc_C_F_P_D3_D3) );
00419    INTREPID_TEST_COMMAND(atools.matvecProductDataField<double>(fc_C_F_P_D3, fc_C_P_D3_D3,    fc_C_F_P_D1) );
00420    INTREPID_TEST_COMMAND(atools.matvecProductDataField<double>(fc_C_F_P_D3, fc_C_P_D3_D3,    fc_F_P_D1) );
00421    // (3) outputFields is (C,F,P,D)
00422    INTREPID_TEST_COMMAND(atools.matvecProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3_D3,    fc_C_F_P_D3) );
00423    INTREPID_TEST_COMMAND(atools.matvecProductDataField<double>(fc_C_F_P,       fc_C_P_D3_D3,    fc_C_F_P_D3) );
00424    INTREPID_TEST_COMMAND(atools.matvecProductDataField<double>(fc_C_F_P_D1,    fc_C_P_D3_D3,    fc_C_F_P_D3) );
00425 
00426    // Cross-check (2): outputFields(C,F,P,D) vs. inputData(C,P,D) and (C,P,D,D): dimensions C, P, D must match
00427    INTREPID_TEST_COMMAND(atools.matvecProductDataField<double>(fc_C_F_P_D3, fc_C1_P_D3_D3,    fc_C_F_P_D3) );
00428    INTREPID_TEST_COMMAND(atools.matvecProductDataField<double>(fc_C_F_P_D3, fc_C_P1_D3_D3,    fc_C_F_P_D3) );
00429    INTREPID_TEST_COMMAND(atools.matvecProductDataField<double>(fc_C_F_P_D3, fc_C_P_D2_D2,     fc_C_F_P_D3) );
00430    INTREPID_TEST_COMMAND(atools.matvecProductDataField<double>(fc_C_F_P_D3, fc_C1_P,          fc_C_F_P_D3) );
00431    INTREPID_TEST_COMMAND(atools.matvecProductDataField<double>(fc_C_F_P_D3, fc_C_P1,          fc_C_F_P_D3) );
00432    INTREPID_TEST_COMMAND(atools.matvecProductDataField<double>(fc_C_F_P_D3, fc_C1_P_D3,       fc_C_F_P_D3) );
00433    INTREPID_TEST_COMMAND(atools.matvecProductDataField<double>(fc_C_F_P_D3, fc_C_P1_D3,       fc_C_F_P_D3) );
00434    INTREPID_TEST_COMMAND(atools.matvecProductDataField<double>(fc_C_F_P_D3, fc_C_P_D2,        fc_C_F_P_D3) );
00435 
00436    // 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
00437    INTREPID_TEST_COMMAND(atools.matvecProductDataField<double>(fc_C_F_P_D3, fc_C_P_D3_D3,    fc_C1_F_P_D3) );
00438    INTREPID_TEST_COMMAND(atools.matvecProductDataField<double>(fc_C_F_P_D3, fc_C_P_D3_D3,    fc_C_F_P1_D3) );
00439    INTREPID_TEST_COMMAND(atools.matvecProductDataField<double>(fc_C_F_P_D3, fc_C_P_D3_D3,    fc_C_F_P_D2) );
00440    INTREPID_TEST_COMMAND(atools.matvecProductDataField<double>(fc_C_F_P_D3, fc_C_P_D3,       fc_C_F_P_D2) );   
00441    
00442    // Cross-check (1): inputData(C,P), (C,P,D) or (C,P,D,D) vs. inputFields(F,P,D): dimensions  P, D must match
00443    INTREPID_TEST_COMMAND(atools.matvecProductDataField<double>(fc_C_F_P_D3, fc_C_P_D3_D3,    fc_F_P1_D3) );
00444    INTREPID_TEST_COMMAND(atools.matvecProductDataField<double>(fc_C_F_P_D3, fc_C_P_D3_D3,    fc_F_P_D2) );
00445    INTREPID_TEST_COMMAND(atools.matvecProductDataField<double>(fc_C_F_P_D3, fc_C_P_D3,       fc_F_P_D2) );
00446 
00447    // Cross-check (3): outputFields(C,F,P,D) vs. inputFields(C,F,P,D): all dimensions C, F, P, D must match
00448    INTREPID_TEST_COMMAND(atools.matvecProductDataField<double>(fc_C_F_P_D3, fc_C_P_D3_D3,    fc_C1_F_P_D3) );
00449    INTREPID_TEST_COMMAND(atools.matvecProductDataField<double>(fc_C_F_P_D3, fc_C_P_D3_D3,    fc_C_F1_P_D3) );
00450    INTREPID_TEST_COMMAND(atools.matvecProductDataField<double>(fc_C_F_P_D3, fc_C_P1_D3_D3,   fc_C_F_P1_D3) );
00451    INTREPID_TEST_COMMAND(atools.matvecProductDataField<double>(fc_C_F_P_D2, fc_C_P_D2_D2,    fc_C_F_P_D3) );
00452 
00453    // Cross-check (3): outputFields(C,F,P,D) vs. inputFields(F,P,D): dimensions F, P, D must match
00454    INTREPID_TEST_COMMAND(atools.matvecProductDataField<double>(fc_C_F_P_D3, fc_C_P_D3_D3,    fc_C_F1_P_D3) );
00455    INTREPID_TEST_COMMAND(atools.matvecProductDataField<double>(fc_C_F_P_D3, fc_C_P1_D3_D3,   fc_C_F_P1_D3) );
00456    INTREPID_TEST_COMMAND(atools.matvecProductDataField<double>(fc_C_F_P_D3, fc_C_P_D3_D3,    fc_C_F_P_D2) );   
00457    
00458    *outStream \
00459      << "\n"
00460      << "===============================================================================\n"\
00461      << "| TEST 6: matvecProductDataData exceptions                                    |\n"\
00462      << "===============================================================================\n";
00463    // 37 exceptions
00464    // (1) inputDataLeft is (C,P), (C,P,D) or (C,P,D,D) and 1 <= D <= 3 is required  
00465    INTREPID_TEST_COMMAND(atools.matvecProductDataData<double>(fc_C_P_D3, fc_C,             fc_C_P_D3) );
00466    INTREPID_TEST_COMMAND(atools.matvecProductDataData<double>(fc_C_P_D3, fc_C_P_D3_D3_D3,  fc_C_P_D3) );
00467    INTREPID_TEST_COMMAND(atools.matvecProductDataData<double>(fc_C_P_D3, fc_C_P_D3_D1,     fc_C_P_D3) );
00468    INTREPID_TEST_COMMAND(atools.matvecProductDataData<double>(fc_C_P_D3, fc_C_P_D1_D3,     fc_C_P_D3) );
00469    
00470    // (2) inputDataRight is (C, P, D) or (P, D) and 1 <= (D=dimension(rank - 1)) <= 3 is required. 
00471    INTREPID_TEST_COMMAND(atools.matvecProductDataData<double>(fc_C_P_D3, fc_C_P_D3_D3,    fc_C_P_D3_D3) );
00472    INTREPID_TEST_COMMAND(atools.matvecProductDataData<double>(fc_C_P_D3, fc_C_P_D3_D3,    fc_P) );
00473    INTREPID_TEST_COMMAND(atools.matvecProductDataData<double>(fc_C_P_D3, fc_C_P_D3_D3,    fc_P_D1) );
00474    INTREPID_TEST_COMMAND(atools.matvecProductDataData<double>(fc_C_P_D3, fc_C_P_D3_D3,    fc_C_P_D1) );
00475    
00476    // (3) outputData is (C,P,D)
00477    INTREPID_TEST_COMMAND(atools.matvecProductDataData<double>(fc_C_P,       fc_C_P_D3_D3,    fc_C_P_D3) );
00478    INTREPID_TEST_COMMAND(atools.matvecProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D3_D3,    fc_C_P_D3) );
00479    INTREPID_TEST_COMMAND(atools.matvecProductDataData<double>(fc_C_P_D1,    fc_C_P_D3_D3,    fc_C_P_D3) );
00480 
00481    // Cross-check (2): outputData(C,P,D) vs. inputDataLeft(C,P), (C,P,D), (C,P,D,D): dimensions C, P, D must match
00482    INTREPID_TEST_COMMAND(atools.matvecProductDataData<double>(fc_C1_P_D3, fc_C_P_D3_D3,  fc_C_P_D3) );
00483    INTREPID_TEST_COMMAND(atools.matvecProductDataData<double>(fc_C_P1_D3, fc_C_P_D3_D3,  fc_C_P_D3) );
00484    INTREPID_TEST_COMMAND(atools.matvecProductDataData<double>(fc_C_P_D2,  fc_C_P_D3_D3,  fc_C_P_D3) );
00485    INTREPID_TEST_COMMAND(atools.matvecProductDataData<double>(fc_C1_P_D3, fc_C_P,        fc_C_P_D3) );
00486    INTREPID_TEST_COMMAND(atools.matvecProductDataData<double>(fc_C_P1_D3, fc_C_P,        fc_C_P_D3) );
00487    INTREPID_TEST_COMMAND(atools.matvecProductDataData<double>(fc_C1_P_D3, fc_C_P_D3,     fc_C_P_D3) );
00488    INTREPID_TEST_COMMAND(atools.matvecProductDataData<double>(fc_C_P1_D3, fc_C_P_D3,     fc_C_P_D3) );
00489    INTREPID_TEST_COMMAND(atools.matvecProductDataData<double>(fc_C_P_D2,  fc_C_P_D3,     fc_C_P_D3) );
00490    
00491    // 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
00492    INTREPID_TEST_COMMAND(atools.matvecProductDataData<double>(fc_C_P_D3, fc_C_P_D3_D3,    fc_C1_P_D3) );
00493    INTREPID_TEST_COMMAND(atools.matvecProductDataData<double>(fc_C_P_D3, fc_C_P_D3_D3,    fc_C_P1_D3) );
00494    INTREPID_TEST_COMMAND(atools.matvecProductDataData<double>(fc_C_P_D3, fc_C_P_D3_D3,    fc_C_P_D2) );   
00495    INTREPID_TEST_COMMAND(atools.matvecProductDataData<double>(fc_C_P_D3, fc_C_P,          fc_C1_P_D3) );
00496    INTREPID_TEST_COMMAND(atools.matvecProductDataData<double>(fc_C_P_D3, fc_C_P,          fc_C_P1_D3) );
00497    INTREPID_TEST_COMMAND(atools.matvecProductDataData<double>(fc_C_P_D3, fc_C_P_D3,       fc_C1_P_D3) );
00498    INTREPID_TEST_COMMAND(atools.matvecProductDataData<double>(fc_C_P_D3, fc_C_P_D3,       fc_C_P1_D3) );
00499    INTREPID_TEST_COMMAND(atools.matvecProductDataData<double>(fc_C_P_D3, fc_C_P_D3,       fc_C_P_D2) );
00500    
00501    // Cross-check (1): inputDataLeft(C,P), (C,P,D), or (C,P,D,D) vs. inputDataRight(P,D): dimensions  P, D must match
00502    INTREPID_TEST_COMMAND(atools.matvecProductDataData<double>(fc_C_P_D3, fc_C_P_D3_D3,    fc_P1_D3) );
00503    INTREPID_TEST_COMMAND(atools.matvecProductDataData<double>(fc_C_P_D3, fc_C_P_D3_D3,    fc_P_D2) );
00504    INTREPID_TEST_COMMAND(atools.matvecProductDataData<double>(fc_C_P_D3, fc_C_P,          fc_P1_D3) );
00505    INTREPID_TEST_COMMAND(atools.matvecProductDataData<double>(fc_C_P_D3, fc_C_P_D3,       fc_P1_D3) );
00506    INTREPID_TEST_COMMAND(atools.matvecProductDataData<double>(fc_C_P_D3, fc_C_P_D3,       fc_P_D2) );
00507    
00508    // Cross-check (3): outputData(C,P,D) vs. inputDataRight(C,P,D): all dimensions C, P, D must match
00509    INTREPID_TEST_COMMAND(atools.matvecProductDataData<double>(fc_C1_P_D3, fc_C1_P_D3_D3,    fc_C_P_D3) );
00510    INTREPID_TEST_COMMAND(atools.matvecProductDataData<double>(fc_C_P1_D3, fc_C_P1_D3_D3,    fc_C_P_D3) );
00511    INTREPID_TEST_COMMAND(atools.matvecProductDataData<double>(fc_C_P_D2,  fc_C_P_D3_D3,     fc_C_P_D3) );
00512 
00513    // Cross-check (3): outputData(C,P,D) vs. inputDataRight(P,D): dimensions P, D must match
00514    INTREPID_TEST_COMMAND(atools.matvecProductDataData<double>(fc_C_P1_D3, fc_C_P1_D3_D3,    fc_P_D3) );
00515    INTREPID_TEST_COMMAND(atools.matvecProductDataData<double>(fc_C_P_D3, fc_C_P_D3_D3,    fc_P_D2) );
00516    
00517    *outStream \
00518      << "\n"
00519      << "===============================================================================\n"\
00520      << "| TEST 7: matmatProductDataField exceptions                                   |\n"\
00521      << "===============================================================================\n";
00522    // 46 exceptions
00523    // (1) inputData is (C,P), (C,P,D), or (C,P,D,D) and 1 <= D <= 3 is required  
00524    INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C,             fc_C_F_P_D3_D3) );
00525    INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3_D3_D3,  fc_C_F_P_D3_D3) );
00526    INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D1_D3,     fc_C_F_P_D3_D3) );
00527    INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3_D1,     fc_C_F_P_D3_D3) );
00528 
00529    // (2) inputFields is (C,F,P,D,D) or (F,P,D,D) and 1 <= (dimension(rank-1), (rank-2)) <= 3 is required. 
00530    INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3_D3,  fc_C_F_P) );
00531    INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3_D3,  fc_C_F_P_D3_D3_D3) );
00532    INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3_D3,  fc_C_F_P_D1_D3) );
00533    INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3_D3,  fc_C_F_P_D3_D1) );
00534 
00535    // (3) outputFields is (C,F,P,D,D)
00536    INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3,       fc_C_P_D3_D3,  fc_C_F_P_D3_D3) );
00537    INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3_D3, fc_C_P_D3_D3,  fc_C_F_P_D3_D3) );
00538    INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D1_D3,    fc_C_P_D3_D3,  fc_C_F_P_D3_D3) );
00539    INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D1,    fc_C_P_D3_D3,  fc_C_F_P_D3_D3) );
00540 
00541    // Cross-check (2): outputFields(C,F,P,D,D) vs. inputData(C,P), (C,P,D), or (C,P,D,D): dimensions C, P, D must match
00542    INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C1_P_D3_D3,  fc_C_F_P_D3_D3) );
00543    INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P1_D3_D3,  fc_C_F_P_D3_D3) );
00544    INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D2_D2,   fc_C_F_P_D3_D3) );
00545    INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P1_D2_D2,  fc_C_F_P_D3_D3) );
00546    INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D2_D2, fc_C_P_D3_D3,   fc_C_F_P_D3_D3) );
00547    INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C1_P,  fc_C_F_P_D3_D3) );
00548    INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P1,  fc_C_F_P_D3_D3) );
00549    INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C1_P_D3,  fc_C_F_P_D3_D3) );
00550    INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P1_D3,  fc_C_F_P_D3_D3) );
00551    INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D1,   fc_C_F_P_D3_D3) );
00552    
00553    // 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
00554    INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3_D3,  fc_C1_F_P_D3_D3) );
00555    INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3_D3,  fc_C_F_P1_D3_D3) );
00556    INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3_D3,  fc_C_F_P_D2_D2) );
00557    INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3_D3,  fc_C_F_P1_D2_D2) );
00558    INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3_D3,  fc_C1_F_P_D2_D2) );
00559    INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3_D3,  fc_C1_F_P1_D2_D2) );
00560    INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P,  fc_C1_F_P_D3_D3) );
00561    INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P,  fc_C_F_P1_D3_D3) );
00562    INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3,  fc_C1_F_P_D3_D3) );
00563    INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3,  fc_C_F_P1_D3_D3) );
00564    INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3,  fc_C_F_P_D2_D2) );
00565       
00566    // 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
00567    INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3_D3,  fc_F_P1_D3_D3) );
00568    INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3_D3,  fc_F_P_D2_D2) );
00569    INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3_D3,  fc_F_P1_D2_D2) );
00570    INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P,    fc_F_P1_D3_D3) );
00571    INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3, fc_F_P1_D3_D3) );
00572    INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3, fc_F_P_D2_D2) );
00573    
00574    // Cross-check (3): outputFields(C,F,P,D,D) vs. inputFields(C,F,P,D,D): all dimensions C, F, P, D must match
00575    INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3_D3,  fc_C1_F_P_D3_D3) );
00576    INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3_D3,  fc_C_F1_P_D3_D3) );
00577    INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3_D3,  fc_C_F_P1_D3_D3) );
00578    INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3_D3,  fc_C_F_P_D2_D2) );
00579 
00580    // Cross-check (3): outputFields(C,F,P,D,D) vs. inputFields(F,P,D,D): dimensions F, P, D must match
00581    INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3_D3,  fc_F1_P_D3_D3) );
00582    INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3_D3,  fc_F_P1_D3_D3) );
00583    INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3_D3,  fc_F_P_D2_D2) );
00584    *outStream \
00585      << "\n"
00586      << "===============================================================================\n"\
00587      << "| TEST 8: matmatProductDataData exceptions                                    |\n"\
00588      << "===============================================================================\n";
00589    // 45 exceptions
00590    // (1) inputDataLeft is (C,P), (C,P,D), or (C,P,D,D) and 2 <= D <= 3 is required  
00591    INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3, fc_C,             fc_C_P_D3_D3) );
00592    INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D3_D3_D3,  fc_C_P_D3_D3) );
00593    INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D1_D3,     fc_C_P_D3_D3) );
00594    INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D3_D1,     fc_C_P_D3_D3) );
00595    
00596    // (2) inputDataRight is (C,P,D,D) or (P,D,D) and 1 <= (dimension(rank-1), (rank-2)) <= 3 is required. 
00597    INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3, fc_C_P,        fc_C_P) );
00598    INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D3,     fc_C_P_D3_D3_D3) );
00599    INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D3_D3,  fc_C_P_D1_D3) );
00600    INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D3_D3,  fc_C_P_D3_D1) );
00601    
00602    // (3) outputData is (C,P,D,D)
00603    INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3,       fc_C_P,        fc_C_P_D3_D3) );
00604    INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3_D3, fc_C_P_D3,     fc_C_P_D3_D3) );
00605    INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D1_D3,    fc_C_P_D3_D3,  fc_C_P_D3_D3) );
00606    INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D1,    fc_C_P_D3_D3,  fc_C_P_D3_D3) );
00607    
00608    // Cross-check (2): outputData(C,P,D,D) vs. inputDataLeft(C,P), (C,P,D), or (C,P,D,D): dimensions C, P, D must match
00609    INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3, fc_C1_P_D3_D3,  fc_C_P_D3_D3) );
00610    INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3, fc_C_P1_D3_D3,  fc_C_P_D3_D3) );
00611    INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D2_D2,   fc_C_P_D3_D3) );
00612    INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3, fc_C_P1_D2_D2,  fc_C_P_D3_D3) );
00613    INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D2_D2, fc_C_P_D3_D3,   fc_C_P_D3_D3) );
00614    INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3, fc_C1_P,        fc_C_P_D3_D3) );
00615    INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3, fc_C_P1,        fc_C_P_D3_D3) );
00616    INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3, fc_C1_P_D3,     fc_C_P_D3_D3) );
00617    INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3, fc_C_P1_D3,     fc_C_P_D3_D3) );
00618    INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D1,      fc_C_P_D3_D3) );
00619    
00620    // 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
00621    INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D3_D3,  fc_C1_P_D3_D3) );
00622    INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D3_D3,  fc_C_P1_D3_D3) );
00623    INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D3_D3,  fc_C_P_D2_D2) );
00624    INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D3_D3,  fc_C_P1_D2_D2) );
00625    INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D3_D3,  fc_C1_P_D2_D2) );
00626    INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D3_D3,  fc_C_P1_D2_D2) );
00627    INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3, fc_C_P,  fc_C1_P_D3_D3) );
00628    INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3, fc_C_P,  fc_C_P1_D3_D3) );
00629    INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D3,  fc_C1_P_D3_D3) );
00630    INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D3,  fc_C_P1_D3_D3) );
00631    INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D3,  fc_C_P_D2_D2) );
00632    
00633    // Cross-check (1): inputDataLeft(C,P), (C,P,D), or (C,P,D,D) vs. inputDataRight(P,D,D): dimensions  P, D must match
00634    INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D3_D3,  fc_P1_D3_D3) );
00635    INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D3_D3,  fc_P_D2_D2) );
00636    INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D3_D3,  fc_P1_D2_D2) );
00637    INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_P_D3_D3, fc_C_P,    fc_P1_D3_D3) );
00638    INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_P_D3_D3, fc_C_P_D3, fc_P1_D3_D3) );
00639    INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_P_D3_D3, fc_C_P_D3, fc_P_D2_D2) );
00640    
00641    // Cross-check (3): outputData(C,P,D,D) vs. inputDataRight(C,P,D,D): all dimensions C, P, D must match
00642    INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D3_D3,  fc_C1_P_D3_D3) );
00643    INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D3_D3,  fc_C1_P_D3_D3) );
00644    INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D3_D3,  fc_C_P1_D3_D3) );
00645    INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D3_D3,  fc_C_P_D2_D2) );
00646    
00647    // Cross-check (3): outputData(C,P,D,D) vs. inputDataRight(P,D,D): dimensions P, D must match
00648    INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D3_D3,  fc_P_D2_D2) );
00649    INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D3_D3,  fc_P1_D3_D3) );
00650   }  
00651 
00652   catch (std::logic_error err) {
00653     *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
00654     *outStream << err.what() << '\n';
00655     *outStream << "-------------------------------------------------------------------------------" << "\n\n";
00656     errorFlag = -1000;
00657   };
00658 
00659   if (Teuchos::TestForException_getThrowNumber() != endThrowNumber)
00660     errorFlag++;
00661 #endif
00662 
00663   /************************************************************************************************
00664     *                                                                                             *
00665     *                                 Operation corectness tests                                  *
00666     *                                                                                             *
00667     ***********************************************************************************************/
00668   
00669   try{
00670   *outStream \
00671     << "\n"
00672     << "===============================================================================\n"\
00673     << "| TEST 1.a: 3D crossProductDataField operations: (C,P,D) and (C,F,P,D)        |\n"\
00674     << "===============================================================================\n";
00675   /*
00676    *        ijkData_1a      ijkFields_1a       Expected result in (C,F,P,D) array
00677    *          i,i          i,i  j,j,  k,k           0, 0    k, k  -j,-j
00678    *          j,j          i,i  j,j,  k,k          -k,-k    0, 0   i, i
00679    *          k,k          i,i  j,j,  k,k           j, j   -i,-i   0, 0
00680    */
00681   
00682   // input data is (C,P,D)
00683   FieldContainer<double> ijkData_1a(3, 2, 3);
00684   // C=0 contains i
00685   ijkData_1a(0, 0, 0) = 1.0;   ijkData_1a(0, 0, 1) = 0.0;   ijkData_1a(0, 0, 2) = 0.0; 
00686   ijkData_1a(0, 1, 0) = 1.0;   ijkData_1a(0, 1, 1) = 0.0;   ijkData_1a(0, 1, 2) = 0.0; 
00687   // C=1 contains j
00688   ijkData_1a(1, 0, 0) = 0.0;   ijkData_1a(1, 0, 1) = 1.0;   ijkData_1a(1, 0, 2) = 0.0; 
00689   ijkData_1a(1, 1, 0) = 0.0;   ijkData_1a(1, 1, 1) = 1.0;   ijkData_1a(1, 1, 2) = 0.0; 
00690   // C=2 contains k
00691   ijkData_1a(2, 0, 0) = 0.0;   ijkData_1a(2, 0, 1) = 0.0;   ijkData_1a(2, 0, 2) = 1.0; 
00692   ijkData_1a(2, 1, 0) = 0.0;   ijkData_1a(2, 1, 1) = 0.0;   ijkData_1a(2, 1, 2) = 1.0; 
00693   
00694   
00695   FieldContainer<double> ijkFields_1a(3, 3, 2, 3);
00696   // C=0, F=0 is i
00697   ijkFields_1a(0, 0, 0, 0) = 1.0; ijkFields_1a(0, 0, 0, 1) = 0.0; ijkFields_1a(0, 0, 0, 2) = 0.0;
00698   ijkFields_1a(0, 0, 1, 0) = 1.0; ijkFields_1a(0, 0, 1, 1) = 0.0; ijkFields_1a(0, 0, 1, 2) = 0.0;
00699   // C=0, F=1 is j
00700   ijkFields_1a(0, 1, 0, 0) = 0.0; ijkFields_1a(0, 1, 0, 1) = 1.0; ijkFields_1a(0, 1, 0, 2) = 0.0;
00701   ijkFields_1a(0, 1, 1, 0) = 0.0; ijkFields_1a(0, 1, 1, 1) = 1.0; ijkFields_1a(0, 1, 1, 2) = 0.0;
00702   // C=0, F=2 is k
00703   ijkFields_1a(0, 2, 0, 0) = 0.0; ijkFields_1a(0, 2, 0, 1) = 0.0; ijkFields_1a(0, 2, 0, 2) = 1.0;
00704   ijkFields_1a(0, 2, 1, 0) = 0.0; ijkFields_1a(0, 2, 1, 1) = 0.0; ijkFields_1a(0, 2, 1, 2) = 1.0;
00705   
00706   // C=1, F=0 is i
00707   ijkFields_1a(1, 0, 0, 0) = 1.0; ijkFields_1a(1, 0, 0, 1) = 0.0; ijkFields_1a(1, 0, 0, 2) = 0.0;
00708   ijkFields_1a(1, 0, 1, 0) = 1.0; ijkFields_1a(1, 0, 1, 1) = 0.0; ijkFields_1a(1, 0, 1, 2) = 0.0;
00709   // C=1, F=1 is j
00710   ijkFields_1a(1, 1, 0, 0) = 0.0; ijkFields_1a(1, 1, 0, 1) = 1.0; ijkFields_1a(1, 1, 0, 2) = 0.0;
00711   ijkFields_1a(1, 1, 1, 0) = 0.0; ijkFields_1a(1, 1, 1, 1) = 1.0; ijkFields_1a(1, 1, 1, 2) = 0.0;
00712   // C=1, F=2 is k
00713   ijkFields_1a(1, 2, 0, 0) = 0.0; ijkFields_1a(1, 2, 0, 1) = 0.0; ijkFields_1a(1, 2, 0, 2) = 1.0;
00714   ijkFields_1a(1, 2, 1, 0) = 0.0; ijkFields_1a(1, 2, 1, 1) = 0.0; ijkFields_1a(1, 2, 1, 2) = 1.0;
00715   
00716   // C=2, F=0 is i
00717   ijkFields_1a(2, 0, 0, 0) = 1.0; ijkFields_1a(2, 0, 0, 1) = 0.0; ijkFields_1a(2, 0, 0, 2) = 0.0;
00718   ijkFields_1a(2, 0, 1, 0) = 1.0; ijkFields_1a(2, 0, 1, 1) = 0.0; ijkFields_1a(2, 0, 1, 2) = 0.0;
00719   // C=2, F=1 is j
00720   ijkFields_1a(2, 1, 0, 0) = 0.0; ijkFields_1a(2, 1, 0, 1) = 1.0; ijkFields_1a(2, 1, 0, 2) = 0.0;
00721   ijkFields_1a(2, 1, 1, 0) = 0.0; ijkFields_1a(2, 1, 1, 1) = 1.0; ijkFields_1a(2, 1, 1, 2) = 0.0;
00722   // C=2, F=2 is k
00723   ijkFields_1a(2, 2, 0, 0) = 0.0; ijkFields_1a(2, 2, 0, 1) = 0.0; ijkFields_1a(2, 2, 0, 2) = 1.0;
00724   ijkFields_1a(2, 2, 1, 0) = 0.0; ijkFields_1a(2, 2, 1, 1) = 0.0; ijkFields_1a(2, 2, 1, 2) = 1.0;
00725   
00726   
00727   FieldContainer<double> outFields(3, 3, 2, 3);
00728   art::crossProductDataField<double>(outFields, ijkData_1a, ijkFields_1a);
00729   
00730   // checks for C = 0
00731   if( !(outFields(0,0,0,0)==0.0 && outFields(0,0,0,1)==0.0 && outFields(0,0,0,2)==0.0 &&
00732         outFields(0,0,1,0)==0.0 && outFields(0,0,1,1)==0.0 && outFields(0,0,1,2)==0.0 ) ) {
00733     *outStream << "\n\nINCORRECT crossProductDataField (1): i x i != 0; ";
00734     errorFlag++;
00735   }    
00736   if( !(outFields(0,1,0,0)==0.0 && outFields(0,1,0,1)==0.0 && outFields(0,1,0,2)==1.0 &&
00737         outFields(0,1,1,0)==0.0 && outFields(0,1,1,1)==0.0 && outFields(0,1,1,2)==1.0 ) ) {
00738     *outStream << "\n\nINCORRECT crossProductDataField (2): i x j != k; ";
00739     errorFlag++;
00740   }
00741   if( !(outFields(0,2,0,0)==0.0 && outFields(0,2,0,1)==-1.0 && outFields(0,2,0,2)==0.0 &&
00742         outFields(0,2,1,0)==0.0 && outFields(0,2,1,1)==-1.0 && outFields(0,2,1,2)==0.0 ) ) {
00743     *outStream << "\n\nINCORRECT crossProductDataField (3): i x k != -j; ";
00744     errorFlag++;
00745   }
00746   
00747   // checks for C = 1
00748   if( !(outFields(1,0,0,0)==0.0 && outFields(1,0,0,1)==0.0 && outFields(1,0,0,2)==-1.0 &&
00749         outFields(1,0,1,0)==0.0 && outFields(1,0,1,1)==0.0 && outFields(1,0,1,2)==-1.0 ) ) {
00750     *outStream << "\n\nINCORRECT crossProductDataField (4): j x i != -k; ";
00751     errorFlag++;
00752   }    
00753   if( !(outFields(1,1,0,0)==0.0 && outFields(1,1,0,1)==0.0 && outFields(1,1,0,2)==0.0 &&
00754         outFields(1,1,1,0)==0.0 && outFields(1,1,1,1)==0.0 && outFields(1,1,1,2)==0.0 ) ) {
00755     *outStream << "\n\nINCORRECT crossProductDataField (5): j x j != 0; ";
00756     errorFlag++;
00757   }
00758   if( !(outFields(1,2,0,0)==1.0 && outFields(1,2,0,1)==0.0 && outFields(1,2,0,2)==0.0 &&
00759         outFields(1,2,1,0)==1.0 && outFields(1,2,1,1)==0.0 && outFields(1,2,1,2)==0.0 ) ) {
00760     *outStream << "\n\nINCORRECT crossProductDataField (6): j x k != i; ";
00761     errorFlag++;
00762   }
00763   
00764   // checks for C = 2
00765   if( !(outFields(2,0,0,0)==0.0 && outFields(2,0,0,1)==1.0 && outFields(2,0,0,2)==0.0 &&
00766         outFields(2,0,1,0)==0.0 && outFields(2,0,1,1)==1.0 && outFields(2,0,1,2)==0.0 ) ) {
00767     *outStream << "\n\nINCORRECT crossProductDataField (7): k x i != j; ";
00768     errorFlag++;
00769   }    
00770   if( !(outFields(2,1,0,0)==-1.0 && outFields(2,1,0,1)==0.0 && outFields(2,1,0,2)==0.0 &&
00771         outFields(2,1,1,0)==-1.0 && outFields(2,1,1,1)==0.0 && outFields(2,1,1,2)==0.0 ) ) {
00772     *outStream << "\n\nINCORRECT crossProductDataField (8): k x j != -i; ";
00773     errorFlag++;
00774   }
00775   if( !(outFields(2,2,0,0)==0.0 && outFields(2,2,0,1)==0.0 && outFields(2,2,0,2)==0.0 &&
00776         outFields(2,2,1,0)==0.0 && outFields(2,2,1,1)==0.0 && outFields(2,2,1,2)==0.0 ) ) {
00777     *outStream << "\n\nINCORRECT crossProductDataField (9): k x k != 0; ";
00778     errorFlag++;
00779   }
00780   
00781   *outStream \
00782     << "\n"
00783     << "===============================================================================\n"\
00784     << "| TEST 1.b: 3D crossProductDataField operations:  (C,P,D) and (F,P,D)         |\n"\
00785     << "===============================================================================\n";
00786   /*
00787    *        ijkData_1b         ijkFields_1b   expected result in (C,F,P,D) array
00788    *          i,i,i               i,j,k                 0, k,-j
00789    *          j,j,j                                    -k, 0, i
00790    *          k,k,k                                     j,-i, 0
00791    */
00792   
00793   // input data is (C,P,D)
00794   FieldContainer<double> ijkData_1b(3, 3, 3);
00795   // C=0 contains i
00796   ijkData_1b(0, 0, 0) = 1.0;   ijkData_1b(0, 0, 1) = 0.0;   ijkData_1b(0, 0, 2) = 0.0; 
00797   ijkData_1b(0, 1, 0) = 1.0;   ijkData_1b(0, 1, 1) = 0.0;   ijkData_1b(0, 1, 2) = 0.0; 
00798   ijkData_1b(0, 2, 0) = 1.0;   ijkData_1b(0, 2, 1) = 0.0;   ijkData_1b(0, 2, 2) = 0.0; 
00799   // C=1 contains j
00800   ijkData_1b(1, 0, 0) = 0.0;   ijkData_1b(1, 0, 1) = 1.0;   ijkData_1b(1, 0, 2) = 0.0; 
00801   ijkData_1b(1, 1, 0) = 0.0;   ijkData_1b(1, 1, 1) = 1.0;   ijkData_1b(1, 1, 2) = 0.0; 
00802   ijkData_1b(1, 2, 0) = 0.0;   ijkData_1b(1, 2, 1) = 1.0;   ijkData_1b(1, 2, 2) = 0.0; 
00803   // C=2 contains k
00804   ijkData_1b(2, 0, 0) = 0.0;   ijkData_1b(2, 0, 1) = 0.0;   ijkData_1b(2, 0, 2) = 1.0; 
00805   ijkData_1b(2, 1, 0) = 0.0;   ijkData_1b(2, 1, 1) = 0.0;   ijkData_1b(2, 1, 2) = 1.0; 
00806   ijkData_1b(2, 2, 0) = 0.0;   ijkData_1b(2, 2, 1) = 0.0;   ijkData_1b(2, 2, 2) = 1.0; 
00807   
00808   // input fields are (F,P,D)
00809   FieldContainer<double> ijkFields_1b(1, 3, 3);
00810   // F=0 at 3 points is (i,j,k)
00811   ijkFields_1b(0, 0, 0) = 1.0; ijkFields_1b(0, 0, 1) = 0.0; ijkFields_1b(0, 0, 2) = 0.0;
00812   ijkFields_1b(0, 1, 0) = 0.0; ijkFields_1b(0, 1, 1) = 1.0; ijkFields_1b(0, 1, 2) = 0.0;
00813   ijkFields_1b(0, 2, 0) = 0.0; ijkFields_1b(0, 2, 1) = 0.0; ijkFields_1b(0, 2, 2) = 1.0;
00814   
00815   // Output array is (C,F,P,D)
00816   outFields.resize(3, 1, 3, 3);
00817   art::crossProductDataField<double>(outFields, ijkData_1b, ijkFields_1b);
00818 
00819   // checks for C = 0
00820   if( !(outFields(0,0,0,0)==0.0 && outFields(0,0,0,1)==0.0 && outFields(0,0,0,2)==0.0) ) {
00821     *outStream << "\n\nINCORRECT crossProductDataField (10): i x i != 0; ";
00822     errorFlag++;
00823   }    
00824   if( !(outFields(0,0,1,0)==0.0 && outFields(0,0,1,1)==0.0 && outFields(0,0,1,2)==1.0) ) {
00825     *outStream << "\n\nINCORRECT crossProductDataField (11): i x j != k; ";
00826     errorFlag++;
00827   }    
00828   if( !(outFields(0,0,2,0)==0.0 && outFields(0,0,2,1)==-1.0 && outFields(0,0,2,2)==0.0) ) {
00829     *outStream << "\n\nINCORRECT crossProductDataField (12): i x k != -j; ";
00830     errorFlag++;
00831   }    
00832 
00833   // checks for C = 1
00834   if( !(outFields(1,0,0,0)==0.0 && outFields(1,0,0,1)==0.0 && outFields(1,0,0,2)==-1.0) ) {
00835     *outStream << "\n\nINCORRECT crossProductDataField (13): j x i != -k; ";
00836     errorFlag++;
00837   }    
00838   if( !(outFields(1,0,1,0)==0.0 && outFields(1,0,1,1)==0.0 && outFields(1,0,1,2)==0.0) ) {
00839     *outStream << "\n\nINCORRECT crossProductDataField (14): j x j != 0; ";
00840     errorFlag++;
00841   }    
00842   if( !(outFields(1,0,2,0)==1.0 && outFields(1,0,2,1)==0.0 && outFields(1,0,2,2)==0.0) ) {
00843     *outStream << "\n\nINCORRECT crossProductDataField (15): j x k != i; ";
00844     errorFlag++;
00845   }    
00846   
00847   // checks for C = 2
00848   if( !(outFields(2,0,0,0)==0.0 && outFields(2,0,0,1)==1.0 && outFields(2,0,0,2)==0.0) ) {
00849     *outStream << "\n\nINCORRECT crossProductDataField (16): k x i != j; ";
00850     errorFlag++;
00851   }    
00852   if( !(outFields(2,0,1,0)==-1.0 && outFields(2,0,1,1)==0.0 && outFields(2,0,1,2)==0.0) ) {
00853     *outStream << "\n\nINCORRECT crossProductDataField (17): k x j != -i; ";
00854     errorFlag++;
00855   }    
00856   if( !(outFields(2,0,2,0)==0.0 && outFields(2,0,2,1)==0.0 && outFields(2,0,2,2)==0.0) ) {
00857     *outStream << "\n\nINCORRECT crossProductDataField (18): k x k != 0; ";
00858     errorFlag++;
00859   }    
00860   
00861   *outStream \
00862     << "\n"
00863     << "===============================================================================\n"\
00864     << "| TEST 1.c: 2D crossProductDataField operations: (C,P,D) and (C,F,P,D)        |\n"\
00865     << "===============================================================================\n";
00866   
00867   /*
00868    *        ijData_1c        ijFields_1c    expected result in (C,F,P) array
00869    *          i,i               i,i  j,j             0, 0   1, 1
00870    *          j,j               i,i  j,j            -1,-1   0, 0
00871    */
00872   // input data is (C,P,D)
00873   FieldContainer<double> ijData_1c(2, 2, 2);
00874   // C=0 contains i
00875   ijData_1c(0, 0, 0) = 1.0;   ijData_1c(0, 0, 1) = 0.0;   
00876   ijData_1c(0, 1, 0) = 1.0;   ijData_1c(0, 1, 1) = 0.0;   
00877   // C=1 contains j
00878   ijData_1c(1, 0, 0) = 0.0;   ijData_1c(1, 0, 1) = 1.0;  
00879   ijData_1c(1, 1, 0) = 0.0;   ijData_1c(1, 1, 1) = 1.0;  
00880   
00881   
00882   FieldContainer<double> ijFields_1c(2, 2, 2, 2);
00883   // C=0, F=0 is i
00884   ijFields_1c(0, 0, 0, 0) = 1.0; ijFields_1c(0, 0, 0, 1) = 0.0; 
00885   ijFields_1c(0, 0, 1, 0) = 1.0; ijFields_1c(0, 0, 1, 1) = 0.0; 
00886   // C=0, F=1 is j
00887   ijFields_1c(0, 1, 0, 0) = 0.0; ijFields_1c(0, 1, 0, 1) = 1.0; 
00888   ijFields_1c(0, 1, 1, 0) = 0.0; ijFields_1c(0, 1, 1, 1) = 1.0; 
00889   
00890   // C=1, F=0 is i
00891   ijFields_1c(1, 0, 0, 0) = 1.0; ijFields_1c(1, 0, 0, 1) = 0.0; 
00892   ijFields_1c(1, 0, 1, 0) = 1.0; ijFields_1c(1, 0, 1, 1) = 0.0; 
00893   // C=1, F=1 is j
00894   ijFields_1c(1, 1, 0, 0) = 0.0; ijFields_1c(1, 1, 0, 1) = 1.0; 
00895   ijFields_1c(1, 1, 1, 0) = 0.0; ijFields_1c(1, 1, 1, 1) = 1.0; 
00896   
00897   // Output array is (C,F,P)
00898   outFields.resize(2, 2, 2);
00899   art::crossProductDataField<double>(outFields, ijData_1c, ijFields_1c);
00900   
00901   if( !(outFields(0,0,0)==0.0 && outFields(0,0,1)==0.0 ) ) {
00902     *outStream << "\n\nINCORRECT crossProductDataField (19): i x i != 0; ";
00903     errorFlag++;
00904   }    
00905   if( !(outFields(0,1,0)==1.0 && outFields(0,1,1)==1.0 ) ) {
00906     *outStream << "\n\nINCORRECT crossProductDataField (20): i x j != 1; ";
00907     errorFlag++;
00908   }    
00909   
00910   if( !(outFields(1,0,0)==-1.0 && outFields(1,0,1)==-1.0 ) ) {
00911     *outStream << "\n\nINCORRECT crossProductDataField (21): j x i != -1; ";
00912     errorFlag++;
00913   }    
00914   if( !(outFields(1,1,0)==0.0 && outFields(1,1,1)==0.0 ) ) {
00915     *outStream << "\n\nINCORRECT crossProductDataField (22): j x j != 0; ";
00916     errorFlag++;
00917   }    
00918   
00919   *outStream \
00920     << "\n"
00921     << "===============================================================================\n"\
00922     << "| TEST 1.d: 2D crossProductDataField operations: (C,P,D) and (F,P,D)          |\n"\
00923     << "===============================================================================\n";
00924   /*
00925    *        ijData_1c      ijFields_1d           expected result in (C,F,P) array
00926    *          i,i               i,j                      0, 1
00927    *          j,j                                       -1, 0
00928    */
00929   // inputFields is (F,P,D)
00930   FieldContainer<double> ijFields_1d(1, 2, 2);
00931   // F=0 at 2 points is i,j
00932   ijFields_1d(0, 0, 0) = 1.0; ijFields_1d(0, 0, 1) = 0.0; 
00933   ijFields_1d(0, 1, 0) = 0.0; ijFields_1d(0, 1, 1) = 1.0; 
00934 
00935   // Output array is (C,F,P)
00936   outFields.resize(2, 1, 2);
00937   art::crossProductDataField<double>(outFields, ijData_1c, ijFields_1d);
00938   
00939   if( !(outFields(0,0,0)==0.0 ) ) {
00940     *outStream << "\n\nINCORRECT crossProductDataField (23): i x i != 0; ";
00941     errorFlag++;
00942   }    
00943   if( !(outFields(0,0,1)==1.0 ) ) {
00944     *outStream << "\n\nINCORRECT crossProductDataField (24): i x j != 1; ";
00945     errorFlag++;
00946   }    
00947   if( !(outFields(1,0,0)==-1.0 ) ) {
00948     *outStream << "\n\nINCORRECT crossProductDataField (25): j x i != -1; ";
00949     errorFlag++;
00950   }    
00951   if( !(outFields(1,0,1)==0.0 ) ) {
00952     *outStream << "\n\nINCORRECT crossProductDataField (26): j x j != 0; ";
00953     errorFlag++;
00954   }    
00955   
00956   
00957   *outStream \
00958     << "\n"
00959     << "===============================================================================\n"\
00960     << "| TEST 2.a: 3D crossProductDataData operations: (C,P,D) and (C,P,D)           |\n"\
00961     << "===============================================================================\n";
00962   /*
00963    *     ijkData_1a    jkiData_2a   kijData_2a          
00964    *        i,i          j,j          k,k          
00965    *        j,j          k,k          i,i         
00966    *        k,k          i,i          j,j           
00967    */
00968   FieldContainer<double> jkiData_2a(3, 2, 3);
00969   // C=0 contains j
00970   jkiData_2a(0, 0, 0) = 0.0;   jkiData_2a(0, 0, 1) = 1.0;   jkiData_2a(0, 0, 2) = 0.0; 
00971   jkiData_2a(0, 1, 0) = 0.0;   jkiData_2a(0, 1, 1) = 1.0;   jkiData_2a(0, 1, 2) = 0.0; 
00972   // C=1 contains k
00973   jkiData_2a(1, 0, 0) = 0.0;   jkiData_2a(1, 0, 1) = 0.0;   jkiData_2a(1, 0, 2) = 1.0; 
00974   jkiData_2a(1, 1, 0) = 0.0;   jkiData_2a(1, 1, 1) = 0.0;   jkiData_2a(1, 1, 2) = 1.0; 
00975   // C=2 contains i
00976   jkiData_2a(2, 0, 0) = 1.0;   jkiData_2a(2, 0, 1) = 0.0;   jkiData_2a(2, 0, 2) = 0.0; 
00977   jkiData_2a(2, 1, 0) = 1.0;   jkiData_2a(2, 1, 1) = 0.0;   jkiData_2a(2, 1, 2) = 0.0; 
00978   
00979   FieldContainer<double> kijData_2a(3, 2, 3);
00980   // C=0 contains k
00981   kijData_2a(0, 0, 0) = 0.0;   kijData_2a(0, 0, 1) = 0.0;   kijData_2a(0, 0, 2) = 1.0; 
00982   kijData_2a(0, 1, 0) = 0.0;   kijData_2a(0, 1, 1) = 0.0;   kijData_2a(0, 1, 2) = 1.0; 
00983   // C=1 contains i
00984   kijData_2a(1, 0, 0) = 1.0;   kijData_2a(1, 0, 1) = 0.0;   kijData_2a(1, 0, 2) = 0.0; 
00985   kijData_2a(1, 1, 0) = 1.0;   kijData_2a(1, 1, 1) = 0.0;   kijData_2a(1, 1, 2) = 0.0; 
00986   // C=2 contains j
00987   kijData_2a(2, 0, 0) = 0.0;   kijData_2a(2, 0, 1) = 1.0;   kijData_2a(2, 0, 2) = 0.0; 
00988   kijData_2a(2, 1, 0) = 0.0;   kijData_2a(2, 1, 1) = 1.0;   kijData_2a(2, 1, 2) = 0.0; 
00989   
00990   
00991   // ijkData_1a x ijkData_1a: outData should contain ixi=0, jxj=0, kxk=0
00992   FieldContainer<double> outData(3,2,3);
00993   art::crossProductDataData<double>(outData, ijkData_1a, ijkData_1a);
00994   
00995   for(int i = 0; i < outData.size(); i++){
00996     if(outData[i] != 0) {
00997       *outStream << "\n\nINCORRECT crossProductDataData (1): i x i, j x j, or k x k != 0; "; 
00998       errorFlag++;
00999     }
01000   }
01001   
01002   
01003   // ijkData_1a x jkiData_2a
01004   art::crossProductDataData<double>(outData, ijkData_1a, jkiData_2a);
01005   
01006   // cell 0 should contain i x j = k
01007   if( !( outData(0,0,0)==0.0 && outData(0,0,1)==0.0 && outData(0,0,2)==1.0 &&
01008          outData(0,1,0)==0.0 && outData(0,1,1)==0.0 && outData(0,1,2)==1.0) ) {
01009     *outStream << "\n\nINCORRECT crossProductDataData (2): i x j != k; ";
01010     errorFlag++;
01011   }    
01012   
01013   // cell 1 should contain j x k = i
01014   if( !( outData(1,0,0)==1.0 && outData(1,0,1)==0.0 && outData(1,0,2)==0.0 &&
01015          outData(1,1,0)==1.0 && outData(1,1,1)==0.0 && outData(1,1,2)==0.0) ) {
01016     *outStream << "\n\nINCORRECT crossProductDataData (3): j x k != i; ";
01017     errorFlag++;
01018   }    
01019   
01020   // cell 2 should contain k x i = j
01021   if( !( outData(2,0,0)==0.0 && outData(2,0,1)==1.0 && outData(2,0,2)==0.0 &&
01022          outData(2,1,0)==0.0 && outData(2,1,1)==1.0 && outData(2,1,2)==0.0) ) {
01023     *outStream << "\n\nINCORRECT crossProductDataData (4): k x i != j; ";
01024     errorFlag++;
01025   }    
01026  
01027   
01028   // ijkData_1a x kijData_2a
01029   art::crossProductDataData<double>(outData, ijkData_1a, kijData_2a);
01030   
01031   // cell 0 should contain i x k = -j
01032   if( !( outData(0,0,0)==0.0 && outData(0,0,1)==-1.0 && outData(0,0,2)==0.0 &&
01033          outData(0,1,0)==0.0 && outData(0,1,1)==-1.0 && outData(0,1,2)==0.0) ) {
01034     *outStream << "\n\nINCORRECT crossProductDataData (5): i x k != -j; ";
01035     errorFlag++;
01036   }    
01037   
01038   // cell 1 should contain j x i = -k
01039   if( !( outData(1,0,0)==0.0 && outData(1,0,1)==0.0 && outData(1,0,2)==-1.0 &&
01040          outData(1,1,0)==0.0 && outData(1,1,1)==0.0 && outData(1,1,2)==-1.0) ) {
01041     *outStream << "\n\nINCORRECT crossProductDataData (6): j x i != -k; ";
01042     errorFlag++;
01043   }    
01044   
01045   // cell 2 should contain k x j = -i
01046   if( !( outData(2,0,0)==-1.0 && outData(2,0,1)==0.0 && outData(2,0,2)==0.0 &&
01047          outData(2,1,0)==-1.0 && outData(2,1,1)==0.0 && outData(2,1,2)==0.0) ) {
01048     *outStream << "\n\nINCORRECT crossProductDataData (7): k x j != -i; ";
01049     errorFlag++;
01050   }    
01051   
01052   
01053   *outStream \
01054     << "\n"
01055     << "===============================================================================\n"\
01056     << "| TEST 2.b: 3D crossProductDataData operations: (C,P,D) and (P,D)             |\n"\
01057     << "===============================================================================\n";
01058   /*
01059    *        ijkData_1b         ijkData_2b        expected result in (C,P,D) array
01060    *          i,i,i               i,j,k               0, k,-j
01061    *          j,j,j                                  -k, 0, i
01062    *          k,k,k                                   j,-i, 0
01063    */
01064   // input data is (P,D)
01065   FieldContainer<double> ijkData_2b(3, 3);
01066   // F=0 at 3 points is (i,j,k)
01067   ijkData_2b(0, 0) = 1.0;   ijkData_2b(0, 1) = 0.0;   ijkData_2b(0, 2) = 0.0;
01068   ijkData_2b(1, 0) = 0.0;   ijkData_2b(1, 1) = 1.0;   ijkData_2b(1, 2) = 0.0;
01069   ijkData_2b(2, 0) = 0.0;   ijkData_2b(2, 1) = 0.0;   ijkData_2b(2, 2) = 1.0;
01070   
01071   // Output array is (C,P,D)
01072   outData.resize(3, 3, 3);
01073   art::crossProductDataData<double>(outData, ijkData_1b, ijkData_2b);
01074   
01075   // checks for C = 0
01076   if( !(outData(0,0,0)==0.0 && outData(0,0,1)==0.0 && outData(0,0,2)==0.0) ) {
01077     *outStream << "\n\nINCORRECT crossProductDataData (5): i x i != 0; ";
01078     errorFlag++;
01079   }    
01080   if( !(outData(0,1,0)==0.0 && outData(0,1,1)==0.0 && outData(0,1,2)==1.0) ) {
01081     *outStream << "\n\nINCORRECT crossProductDataData (6): i x j != k; ";
01082     errorFlag++;
01083   }    
01084   if( !(outData(0,2,0)==0.0 && outData(0,2,1)==-1.0 && outData(0,2,2)==0.0) ) {
01085     *outStream << "\n\nINCORRECT crossProductDataData (7): i x k != -j; ";
01086     errorFlag++;
01087   }    
01088   
01089   // checks for C = 1
01090   if( !(outData(1,0,0)==0.0 && outData(1,0,1)==0.0 && outData(1,0,2)==-1.0) ) {
01091     *outStream << "\n\nINCORRECT crossProductDataData (8): j x i != -k; ";
01092     errorFlag++;
01093   }    
01094   if( !(outData(1,1,0)==0.0 && outData(1,1,1)==0.0 && outData(1,1,2)==0.0) ) {
01095     *outStream << "\n\nINCORRECT crossProductDataData (9): j x j != 0; ";
01096     errorFlag++;
01097   }    
01098   if( !(outData(1,2,0)==1.0 && outData(1,2,1)==0.0 && outData(1,2,2)==0.0) ) {
01099     *outStream << "\n\nINCORRECT crossProductDataData (10): j x k != i; ";
01100     errorFlag++;
01101   }    
01102   
01103   // checks for C = 2
01104   if( !(outData(2,0,0)==0.0 && outData(2,0,1)==1.0 && outData(2,0,2)==0.0) ) {
01105     *outStream << "\n\nINCORRECT crossProductDataData (11): k x i != j; ";
01106     errorFlag++;
01107   }    
01108   if( !(outData(2,1,0)==-1.0 && outData(2,1,1)==0.0 && outData(2,1,2)==0.0) ) {
01109     *outStream << "\n\nINCORRECT crossProductDataData (12): k x j != -i; ";
01110     errorFlag++;
01111   }    
01112   if( !(outData(2,2,0)==0.0 && outData(2,2,1)==0.0 && outData(2,2,2)==0.0) ) {
01113     *outStream << "\n\nINCORRECT crossProductDataData (13): k x k != 0; ";
01114     errorFlag++;
01115   }    
01116   
01117   
01118   *outStream \
01119     << "\n"
01120     << "===============================================================================\n"\
01121     << "| TEST 2.c: 2D crossProductDataData operations: (C,P,D) and (C,P,D)           |\n"\
01122     << "===============================================================================\n";
01123   /*
01124    *     ijData_1c    jiData_2c             
01125    *        i,i          j,j                   
01126    *        j,j          i,i                   
01127    */
01128   FieldContainer<double> jiData_2c(2, 2, 2);
01129   // C=0 contains j
01130   jiData_2c(0, 0, 0) = 0.0;   jiData_2c(0, 0, 1) = 1.0;   
01131   jiData_2c(0, 1, 0) = 0.0;   jiData_2c(0, 1, 1) = 1.0;   
01132   // C=1 contains i
01133   jiData_2c(1, 0, 0) = 1.0;   jiData_2c(1, 0, 1) = 0.0;  
01134   jiData_2c(1, 1, 0) = 1.0;   jiData_2c(1, 1, 1) = 0.0;  
01135   
01136   
01137   // ijData_1c x ijData_1c: outData should contain ixi=0, jxj=0
01138   outData.resize(2,2);
01139   art::crossProductDataData<double>(outData, ijData_1c, ijData_1c);
01140   
01141   for(int i = 0; i < outData.size(); i++){
01142     if(outData[i] != 0) {
01143       *outStream << "\n\nINCORRECT crossProductDataData (14): i x i or j x j != 0; "; 
01144       errorFlag++;
01145     }
01146   }
01147   
01148   // ijData_1c x jiData_1c: outData should contain ixi=0, jxj=0
01149   art::crossProductDataData<double>(outData, ijData_1c, jiData_2c);
01150  
01151   if( !(outData(0,0)==1.0 && outData(0,1)==1.0 ) ) {
01152     *outStream << "\n\nINCORRECT crossProductDataData (15): i x j != 1; ";
01153     errorFlag++;
01154   }    
01155   if( !(outData(1,0)==-1.0 && outData(1,1)==-1.0 ) ) {
01156     *outStream << "\n\nINCORRECT crossProductDataData (16): j x i != -1; ";
01157     errorFlag++;
01158   }    
01159   
01160   *outStream \
01161     << "\n"
01162     << "===============================================================================\n"\
01163     << "| TEST 2.d: 2D crossProductDataData operations: (C,P,D) and (P,D)             |\n"\
01164     << "===============================================================================\n";
01165   /*
01166    *        ijData_1c      ijData_2d        expected result in (C,P) array
01167    *          i,i             i,j                 0, 1
01168    *          j,j                                -1, 0
01169    */
01170   FieldContainer<double> ijData_2d(2, 2);
01171   ijData_2d(0, 0) = 1.0;   ijData_2d(0, 1) = 0.0; 
01172   ijData_2d(1, 0) = 0.0;   ijData_2d(1, 1) = 1.0; 
01173   
01174   outData.resize(2,2);
01175   art::crossProductDataData<double>(outData, ijData_1c, ijData_2d);
01176 
01177   if( !(outData(0,0)==0.0 ) ) {
01178     *outStream << "\n\nINCORRECT crossProductDataData (17): i x i != 0; ";
01179     errorFlag++;
01180   }    
01181   if( !(outData(0,1)==1.0 ) ) {
01182     *outStream << "\n\nINCORRECT crossProductDataData (18): i x j != 1; ";
01183     errorFlag++;
01184   }    
01185   if( !(outData(1,0)==-1.0 ) ) {
01186     *outStream << "\n\nINCORRECT crossProductDataData (19): j x i != -1; ";
01187     errorFlag++;
01188   }    
01189   if( !(outData(1,1)==0.0 ) ) {
01190     *outStream << "\n\nINCORRECT crossProductDataData (20): j x j != 0; ";
01191     errorFlag++;
01192   }    
01193   
01194   
01195   *outStream \
01196     << "\n"
01197     << "===============================================================================\n"\
01198     << "| TEST 3.a: outerProductDataField operations: (C,P,D) and (C,F,P,D)           |\n"\
01199     << "===============================================================================\n";
01200   /*
01201    *        ijkData_1a      ijkFields_1a       Expected result in (C,F,P,D,D) array:
01202    *          i,i          i,i  j,j,  k,k            (0,0) (0,1) (0,2)    
01203    *          j,j          i,i  j,j,  k,k            (1,0) (1,1) (1,2)     
01204    *          k,k          i,i  j,j,  k,k            (2,0) (2,1) (2,2)    
01205    *   Indicates the only non-zero element of (C,F,P,*,*), specifically, 
01206    *   element with row = cell and col = field should equal 1; all other should equal 0                        
01207    */
01208   
01209   outFields.resize(3, 3, 2, 3, 3);
01210   art::outerProductDataField<double>(outFields, ijkData_1a, ijkFields_1a);
01211 
01212   for(int cell = 0; cell < ijkData_1a.dimension(0); cell++){
01213     for(int field = 0; field < ijkFields_1a.dimension(1); field++){
01214       for(int point = 0; point < ijkData_1a.dimension(1); point++){
01215         for(int row = 0; row < ijkData_1a.dimension(2); row++){
01216           for(int col = 0; col < ijkData_1a.dimension(2); col++){
01217             
01218             // element with row = cell and col = field should equal 1; all other should equal 0
01219             if( (row == cell && col == field) ){
01220               if(outFields(cell, field, point, row, col) != 1.0) {
01221                 *outStream << "\n\nINCORRECT outerProductDataField (1): computed value is " 
01222                 << outFields(cell, field, point, row, col) << " whereas correct value is 1.0";
01223                 errorFlag++;
01224               }
01225             }
01226             else {
01227               if(outFields(cell, field, point, row, col) != 0.0) {
01228                 *outStream << "\n\nINCORRECT outerProductDataField (2): computed value is " 
01229                 << outFields(cell, field, point, row, col) << " whereas correct value is 0.0";
01230                 errorFlag++;
01231               }
01232             } // if
01233           }// col
01234         }// row
01235       }// point
01236     }// field
01237   }// cell
01238   
01239   *outStream \
01240     << "\n"
01241     << "===============================================================================\n"\
01242     << "| TEST 3.b: outerProductDataField operations: (C,P,D) and (F,P,D)             |\n"\
01243     << "===============================================================================\n";
01244   /*
01245    *        ijkData_1b         ijkFields_1b   expected result in (C,F,P,D,D) array
01246    *          i,i,i               i,j,k            (0,0) (0,1) (0,2)
01247    *          j,j,j                                (1,0) (1,1) (1,2)
01248    *          k,k,k                                (2,0) (2,1) (2,2)
01249    *   Indicates the only non-zero element of (C,F,P,*,*), specifically, 
01250    *   element with row = cell and col = point should equal 1; all other should equal 0                        
01251    */
01252   
01253   outFields.resize(3, 1, 3, 3, 3);
01254   art::outerProductDataField<double>(outFields, ijkData_1b, ijkFields_1b);
01255   
01256   for(int cell = 0; cell < ijkData_1b.dimension(0); cell++){
01257     for(int field = 0; field < ijkFields_1b.dimension(0); field++){
01258       for(int point = 0; point < ijkData_1b.dimension(1); point++){
01259         for(int row = 0; row < ijkData_1b.dimension(2); row++){
01260           for(int col = 0; col < ijkData_1b.dimension(2); col++){
01261             
01262             // element with row = cell and col = point should equal 1; all other should equal 0
01263             if( (row == cell && col == point) ){
01264               if(outFields(cell, field, point, row, col) != 1.0) {
01265                 *outStream << "\n\nINCORRECT outerProductDataField (3): computed value is " 
01266                 << outFields(cell, field, point, row, col) << " whereas correct value is 1.0";
01267                 errorFlag++;
01268               
01269               }
01270             }
01271             else {
01272               if(outFields(cell, field, point, row, col) != 0.0) {
01273                 *outStream << "\n\nINCORRECT outerProductDataField (4): computed value is " 
01274                 << outFields(cell, field, point, row, col) << " whereas correct value is 0.0";
01275                 errorFlag++;
01276               }
01277             } // if
01278           }// col
01279         }// row
01280       }// point
01281     }// field
01282   }// cell
01283 
01284   *outStream \
01285     << "\n"
01286     << "===============================================================================\n"\
01287     << "| TEST 4.a: outerProductDataData operations: (C,P,D) and (C,P,D)              |\n"\
01288     << "===============================================================================\n";
01289   /*
01290    *     ijkData_1a    jkiData_2a   kijData_2a          
01291    *        i,i          j,j          k,k          
01292    *        j,j          k,k          i,i         
01293    *        k,k          i,i          j,j    
01294    *     
01295    *     Expected results are stated with each test case.
01296    */
01297   outData.resize(3, 2, 3, 3);
01298  
01299   art::outerProductDataData<double>(outData, ijkData_1a, ijkData_1a);
01300   for(int cell = 0; cell < ijkData_1a.dimension(0); cell++){
01301       for(int point = 0; point < ijkData_1a.dimension(1); point++){
01302         for(int row = 0; row < ijkData_1a.dimension(2); row++){
01303           for(int col = 0; col < ijkData_1a.dimension(2); col++){
01304             
01305             // element with row = cell and col = cell should equal 1; all other should equal 0
01306             if( (row == cell && col == cell) ){
01307               if(outData(cell, point, row, col) != 1.0) {
01308                 *outStream << "\n\nINCORRECT outerProductDataData (1): computed value is " 
01309                 << outData(cell, point, row, col) << " whereas correct value is 1.0";
01310                 errorFlag++;
01311               }
01312             }
01313             else {
01314               if(outData(cell, point, row, col) != 0.0) {
01315                 *outStream << "\n\nINCORRECT outerProductDataData (2): computed value is " 
01316                 << outData(cell, point, row, col) << " whereas correct value is 0.0";
01317                 errorFlag++;
01318               }
01319             } // if
01320           }// col
01321         }// row
01322       }// point
01323   }// cell
01324   
01325   outData.initialize();
01326   art::outerProductDataData<double>(outData, ijkData_1a, jkiData_2a);
01327   for(int cell = 0; cell < ijkData_1a.dimension(0); cell++){
01328     for(int point = 0; point < ijkData_1a.dimension(1); point++){
01329       for(int row = 0; row < ijkData_1a.dimension(2); row++){
01330         for(int col = 0; col < ijkData_1a.dimension(2); col++){
01331           
01332           // element with row = cell and col = cell + 1 (mod 3) should equal 1; all other should equal 0
01333           if( (row == cell && col == (cell + 1) % 3) ){
01334             if(outData(cell, point, row, col) != 1.0) {
01335               *outStream << "\n\nINCORRECT outerProductDataData (3): computed value is " 
01336               << outData(cell, point, row, col) << " whereas correct value is 1.0";
01337               errorFlag++;
01338             }
01339           }
01340           else {
01341             if(outData(cell, point, row, col) != 0.0) {
01342               *outStream << "\n\nINCORRECT outerProductDataData (4): computed value is " 
01343               << outData(cell, point, row, col) << " whereas correct value is 0.0";
01344               errorFlag++;
01345             }
01346           } // if
01347         }// col
01348       }// row
01349     }// point
01350   }// cell
01351   
01352   
01353   outData.initialize();
01354   art::outerProductDataData<double>(outData, ijkData_1a, kijData_2a);
01355   for(int cell = 0; cell < ijkData_1a.dimension(0); cell++){
01356     for(int point = 0; point < ijkData_1a.dimension(1); point++){
01357       for(int row = 0; row < ijkData_1a.dimension(2); row++){
01358         for(int col = 0; col < ijkData_1a.dimension(2); col++){
01359           
01360           // element with row = cell and col = cell + 2 (mod 3) should equal 1; all other should equal 0
01361           if( (row == cell && col == (cell + 2) % 3) ){
01362             if(outData(cell, point, row, col) != 1.0) {
01363               *outStream << "\n\nINCORRECT outerProductDataData (5): computed value is " 
01364               << outData(cell, point, row, col) << " whereas correct value is 1.0";
01365               errorFlag++;
01366             }
01367           }
01368           else {
01369             if(outData(cell, point, row, col) != 0.0) {
01370               *outStream << "\n\nINCORRECT outerProductDataData (6): computed value is " 
01371               << outData(cell, point, row, col) << " whereas correct value is 0.0";
01372               errorFlag++;
01373             }
01374           } // if
01375         }// col
01376       }// row
01377     }// point
01378   }// cell
01379   
01380   
01381   *outStream \
01382     << "\n"
01383     << "===============================================================================\n"\
01384     << "| TEST 4.b: outerProductDataData operations: (C,P,D) and (P,D)                |\n"\
01385     << "===============================================================================\n";
01386   /*
01387    *        ijkData_1b         ijkData_2b   expected result in (C,P,D,D) array
01388    *          i,i,i               i,j,k            (0,0) (0,1) (0,2)
01389    *          j,j,j                                (1,0) (1,1) (1,2)
01390    *          k,k,k                                (2,0) (2,1) (2,2)
01391    *   Indicates the only non-zero element of (C,P,*,*), specifically, 
01392    *   element with row = cell and col = point should equal 1; all other should equal 0                        
01393    */
01394   outData.resize(3,3,3,3);
01395   art::outerProductDataData<double>(outData, ijkData_1b, ijkData_2b);
01396   for(int cell = 0; cell < ijkData_1b.dimension(0); cell++){
01397     for(int point = 0; point < ijkData_1b.dimension(1); point++){
01398       for(int row = 0; row < ijkData_1b.dimension(2); row++){
01399         for(int col = 0; col < ijkData_1b.dimension(2); col++){
01400           
01401           // element with row = cell and col = cell + 2 (mod 3) should equal 1; all other should equal 0
01402           if( (row == cell && col == point) ){
01403             if(outData(cell, point, row, col) != 1.0) {
01404               *outStream << "\n\nINCORRECT outerProductDataData (7): computed value is " 
01405               << outData(cell, point, row, col) << " whereas correct value is 1.0";
01406               errorFlag++;
01407             }
01408           }
01409           else {
01410             if(outData(cell, point, row, col) != 0.0) {
01411               *outStream << "\n\nINCORRECT outerProductDataData (8): computed value is " 
01412               << outData(cell, point, row, col) << " whereas correct value is 0.0";
01413               errorFlag++;
01414             }
01415           } // if
01416         }// col
01417       }// row
01418     }// point
01419   }// cell
01420   
01421   *outStream \
01422     << "\n"
01423     << "===============================================================================\n"\
01424     << "| TEST 5.a: matvecProductDataField operations: (C,P,D,D) and (C,F,P,D)        |\n"\
01425     << "===============================================================================\n";
01426   /*
01427    *         inputMat              inputVecFields           outFields
01428    *     1  1  1     0  0  0     0   1       -1  -1      0   3      0   0
01429    *    -1  2 -1    -1 -2 -3     0   1       -1   1      0   0      6   2
01430    *     1  2  3    -2  6 -4     0   1       -1  -1      0   6      0  12
01431    */
01432   
01433   // (C,P,D,D)
01434   FieldContainer<double> inputMat(2,1,3,3);
01435   // cell 0
01436   inputMat(0,0,0,0) = 1.0;  inputMat(0,0,0,1) = 1.0;  inputMat(0,0,0,2) = 1.0;
01437   inputMat(0,0,1,0) =-1.0;  inputMat(0,0,1,1) = 2.0;  inputMat(0,0,1,2) =-1.0;
01438   inputMat(0,0,2,0) = 1.0;  inputMat(0,0,2,1) = 2.0;  inputMat(0,0,2,2) = 3.0;
01439   // cell 1
01440   inputMat(1,0,0,0) = 0.0;  inputMat(1,0,0,1) = 0.0;  inputMat(1,0,0,2) = 0.0;
01441   inputMat(1,0,1,0) =-1.0;  inputMat(1,0,1,1) =-2.0;  inputMat(1,0,1,2) =-3.0;
01442   inputMat(1,0,2,0) =-2.0;  inputMat(1,0,2,1) = 6.0;  inputMat(1,0,2,2) =-4.0;
01443   
01444   // (C,F,P,D)
01445   FieldContainer<double> inputVecFields(2,2,1,3);
01446   // cell 0; fields 0,1
01447   inputVecFields(0,0,0,0) = 0.0;  inputVecFields(0,0,0,1) = 0.0;  inputVecFields(0,0,0,2) = 0.0;
01448   inputVecFields(0,1,0,0) = 1.0;  inputVecFields(0,1,0,1) = 1.0;  inputVecFields(0,1,0,2) = 1.0;
01449   // cell 1; fields 0,1
01450   inputVecFields(1,0,0,0) =-1.0;  inputVecFields(1,0,0,1) =-1.0;  inputVecFields(1,0,0,2) =-1.0;
01451   inputVecFields(1,1,0,0) =-1.0;  inputVecFields(1,1,0,1) = 1.0;  inputVecFields(1,1,0,2) =-1.0;
01452   
01453   // (C,F,P,D) - true
01454   FieldContainer<double> outFieldsCorrect(2,2,1,3);
01455   // cell 0; fields 0,1
01456   outFieldsCorrect(0,0,0,0) = 0.0;  outFieldsCorrect(0,0,0,1) = 0.0;  outFieldsCorrect(0,0,0,2) = 0.0;
01457   outFieldsCorrect(0,1,0,0) = 3.0;  outFieldsCorrect(0,1,0,1) = 0.0;  outFieldsCorrect(0,1,0,2) = 6.0;
01458   // cell 1; fields 0,1  
01459   outFieldsCorrect(1,0,0,0) = 0.0;  outFieldsCorrect(1,0,0,1) = 6.0;  outFieldsCorrect(1,0,0,2) = 0.0;
01460   outFieldsCorrect(1,1,0,0) = 0.0;  outFieldsCorrect(1,1,0,1) = 2.0;  outFieldsCorrect(1,1,0,2) = 12.0;
01461 
01462   // (C,F,P,D)
01463   outFields.resize(2,2,1,3);
01464   art::matvecProductDataField<double>(outFields, inputMat, inputVecFields);
01465   
01466   // test loop
01467   for(int cell = 0; cell < outFields.dimension(0); cell++){
01468     for(int field = 0; field < outFields.dimension(1); field++){
01469       for(int point = 0; point < outFields.dimension(2); point++){
01470         for(int row = 0; row < outFields.dimension(3); row++){
01471           if(outFields(cell, field, point, row) != outFieldsCorrect(cell, field, point, row)) {
01472             *outStream << "\n\nINCORRECT matvecProductDataField (1): \n value at multi-index ("
01473             << cell << "," << field << "," << point << "," << row << ") = " 
01474             << outFields(cell, field, point, row) << " but correct value is " 
01475             << outFieldsCorrect(cell, field, point, row) <<"\n";
01476             errorFlag++; 
01477           }
01478         }//row
01479       }// point
01480     }// field
01481   }// cell
01482   
01483 
01484   *outStream \
01485     << "\n"
01486     << "===============================================================================\n"\
01487     << "| TEST 5.b: matvecProductDataField operations: (C,P,D,D) and (F,P,D)          |\n"\
01488     << "===============================================================================\n";
01489   /*
01490    *         inputMat              inputVecFields        outFields
01491    *     1  1  1     0  0  0       0   1  -1  -1       0  3 -3 -1     0  0  0  0
01492    *    -1  2 -1    -1 -2 -3       0   1  -1   1       0  0  0  4     0 -6  6  2
01493    *     1  2  3    -2  6 -4       0   1  -1  -1       0  6 -6 -2     0  0  0 12
01494    *     Use the same 4 vector fields as above but formatted as (F,P,D) array
01495    */
01496   // (C,F,P,D)
01497   inputVecFields.resize(4,1,3);
01498   // fields 0,1,2,3
01499   inputVecFields(0,0,0) = 0.0;  inputVecFields(0,0,1) = 0.0;  inputVecFields(0,0,2) = 0.0;
01500   inputVecFields(1,0,0) = 1.0;  inputVecFields(1,0,1) = 1.0;  inputVecFields(1,0,2) = 1.0;
01501   inputVecFields(2,0,0) =-1.0;  inputVecFields(2,0,1) =-1.0;  inputVecFields(2,0,2) =-1.0;
01502   inputVecFields(3,0,0) =-1.0;  inputVecFields(3,0,1) = 1.0;  inputVecFields(3,0,2) =-1.0;
01503   
01504   // (C,F,P,D) - true
01505   outFieldsCorrect.resize(2,4,1,3);
01506   // cell 0; fields 0,1,2,3
01507   outFieldsCorrect(0,0,0,0) = 0.0;  outFieldsCorrect(0,0,0,1) = 0.0;  outFieldsCorrect(0,0,0,2) = 0.0;
01508   outFieldsCorrect(0,1,0,0) = 3.0;  outFieldsCorrect(0,1,0,1) = 0.0;  outFieldsCorrect(0,1,0,2) = 6.0;
01509   outFieldsCorrect(0,2,0,0) =-3.0;  outFieldsCorrect(0,2,0,1) = 0.0;  outFieldsCorrect(0,2,0,2) =-6.0;
01510   outFieldsCorrect(0,3,0,0) =-1.0;  outFieldsCorrect(0,3,0,1) = 4.0;  outFieldsCorrect(0,3,0,2) =-2.0;
01511   // cell 1; fields 0,1,2,3  
01512   outFieldsCorrect(1,0,0,0) = 0.0;  outFieldsCorrect(1,0,0,1) = 0.0;  outFieldsCorrect(1,0,0,2) = 0.0;
01513   outFieldsCorrect(1,1,0,0) = 0.0;  outFieldsCorrect(1,1,0,1) =-6.0;  outFieldsCorrect(1,1,0,2) = 0.0;
01514   outFieldsCorrect(1,2,0,0) = 0.0;  outFieldsCorrect(1,2,0,1) = 6.0;  outFieldsCorrect(1,2,0,2) = 0.0;
01515   outFieldsCorrect(1,3,0,0) = 0.0;  outFieldsCorrect(1,3,0,1) = 2.0;  outFieldsCorrect(1,3,0,2) =12.0;
01516 
01517   // (C,F,P,D)
01518   outFields.resize(2,4,1,3);
01519   art::matvecProductDataField<double>(outFields, inputMat, inputVecFields);
01520   
01521   // test loop
01522   for(int cell = 0; cell < outFields.dimension(0); cell++){
01523     for(int field = 0; field < outFields.dimension(1); field++){
01524       for(int point = 0; point < outFields.dimension(2); point++){
01525         for(int row = 0; row < outFields.dimension(3); row++){
01526           if(outFields(cell, field, point, row) != outFieldsCorrect(cell, field, point, row)) {
01527             *outStream << "\n\nINCORRECT matvecProductDataField (2): \n value at multi-index ("
01528             << cell << "," << field << "," << point << "," << row << ") = " 
01529             << outFields(cell, field, point, row) << " but correct value is " 
01530             << outFieldsCorrect(cell, field, point, row) <<"\n";
01531             errorFlag++; 
01532           }
01533         }//row
01534       }// point
01535     }// field
01536   }// cell
01537   
01538   
01539   *outStream \
01540     << "\n"
01541     << "===============================================================================\n"\
01542     << "| TEST 5.c: matvecProductDataField random tests: branch inputFields(C,F,P,D)  |\n"\
01543     << "===============================================================================\n";
01544   /*
01545    *  d1 is spatial dimension and should be 1, 2 or 3. If d1>3, the RealSpaceTools function 'inverse' will fail
01546    */
01547   {// test 5.c scope
01548     int c=5, p=9, f=7, d1=3;
01549     double zero = INTREPID_TOL*10000.0;
01550     
01551     FieldContainer<double> in_c_f_p_d(c, f, p, d1);
01552     FieldContainer<double> out_c_f_p_d(c, f, p, d1);
01553     FieldContainer<double> outi_c_f_p_d(c, f, p, d1);
01554     
01555     FieldContainer<double> data_c_p(c, p);
01556     FieldContainer<double> datainv_c_p(c, p);
01557     FieldContainer<double> data_c_1(c, 1);
01558     FieldContainer<double> datainv_c_1(c, 1);
01559     FieldContainer<double> data_c_p_d(c, p, d1);
01560     FieldContainer<double> datainv_c_p_d(c, p, d1);
01561     FieldContainer<double> data_c_1_d(c, 1, d1);
01562     FieldContainer<double> datainv_c_1_d(c, 1, d1);
01563     FieldContainer<double> data_c_p_d_d(c, p, d1, d1);
01564     FieldContainer<double> datainv_c_p_d_d(c, p, d1, d1);
01565     FieldContainer<double> datainvtrn_c_p_d_d(c, p, d1, d1);
01566     FieldContainer<double> data_c_1_d_d(c, 1, d1, d1);
01567     FieldContainer<double> datainv_c_1_d_d(c, 1, d1, d1);
01568     FieldContainer<double> datainvtrn_c_1_d_d(c, 1, d1, d1);
01569     /***********************************************************************************************
01570       *                          Constant diagonal tensor: inputData(C,P)                          *
01571       **********************************************************************************************/
01572     for (int i=0; i<in_c_f_p_d.size(); i++) {
01573       in_c_f_p_d[i] = Teuchos::ScalarTraits<double>::random();
01574     }
01575     for (int i=0; i<data_c_p.size(); i++) {
01576       data_c_p[i] = Teuchos::ScalarTraits<double>::random();
01577       datainv_c_p[i] = 1.0 / data_c_p[i];
01578     }
01579     for (int i=0; i<data_c_1.size(); i++) {
01580       data_c_1[i] = Teuchos::ScalarTraits<double>::random();
01581       datainv_c_1[i] = 1.0 / data_c_1[i];
01582     }
01583     // Tensor values vary by point:
01584     art::matvecProductDataField<double>(out_c_f_p_d, data_c_p, in_c_f_p_d);
01585     art::matvecProductDataField<double>(out_c_f_p_d, datainv_c_p, out_c_f_p_d);
01586     rst::subtract(&out_c_f_p_d[0], &in_c_f_p_d[0], out_c_f_p_d.size());
01587     if (rst::vectorNorm(&out_c_f_p_d[0], out_c_f_p_d.size(), NORM_ONE) > zero) {
01588       *outStream << "\n\nINCORRECT matvecProductDataField (3): check scalar inverse property\n\n";
01589       errorFlag = -1000;
01590     }
01591     // Tensor values do not vary by point:
01592     art::matvecProductDataField<double>(out_c_f_p_d, data_c_1, in_c_f_p_d);
01593     art::matvecProductDataField<double>(out_c_f_p_d, datainv_c_1, out_c_f_p_d);
01594     rst::subtract(&out_c_f_p_d[0], &in_c_f_p_d[0], out_c_f_p_d.size());
01595     if (rst::vectorNorm(&out_c_f_p_d[0], out_c_f_p_d.size(), NORM_ONE) > zero) {
01596       *outStream << "\n\nINCORRECT matvecProductDataField (4): check scalar inverse property\n\n";
01597       errorFlag = -1000;
01598     }
01599     /***********************************************************************************************
01600       *                     Non-onstant diagonal tensor: inputData(C,P,D)                          *
01601       **********************************************************************************************/
01602     for (int i=0; i<in_c_f_p_d.size(); i++) {
01603       in_c_f_p_d[i] = Teuchos::ScalarTraits<double>::random();
01604     }
01605     for (int i=0; i<data_c_p_d.size(); i++) {
01606       data_c_p_d[i] = Teuchos::ScalarTraits<double>::random();
01607       datainv_c_p_d[i] = 1.0 / data_c_p_d[i];
01608     }
01609     for (int i=0; i<data_c_1_d.size(); i++) {
01610       data_c_1_d[i] = Teuchos::ScalarTraits<double>::random();
01611       datainv_c_1_d[i] = 1.0 / data_c_1_d[i];
01612     }
01613     // Tensor values vary by point:
01614     art::matvecProductDataField<double>(out_c_f_p_d, data_c_p_d, in_c_f_p_d);
01615     art::matvecProductDataField<double>(out_c_f_p_d, datainv_c_p_d, out_c_f_p_d);
01616     rst::subtract(&out_c_f_p_d[0], &in_c_f_p_d[0], out_c_f_p_d.size());
01617     if (rst::vectorNorm(&out_c_f_p_d[0], out_c_f_p_d.size(), NORM_ONE) > zero) {
01618       *outStream << "\n\nINCORRECT matvecProductDataField (5): check scalar inverse property\n\n";
01619       errorFlag = -1000;
01620     }
01621     // Tensor values do not vary by point
01622     art::matvecProductDataField<double>(out_c_f_p_d, data_c_1_d, in_c_f_p_d);
01623     art::matvecProductDataField<double>(out_c_f_p_d, datainv_c_1_d, out_c_f_p_d);
01624     rst::subtract(&out_c_f_p_d[0], &in_c_f_p_d[0], out_c_f_p_d.size());
01625     if (rst::vectorNorm(&out_c_f_p_d[0], out_c_f_p_d.size(), NORM_ONE) > zero) {
01626       *outStream << "\n\nINCORRECT matvecProductDataField (6): check scalar inverse property\n\n";
01627       errorFlag = -1000;
01628     }
01629     /***********************************************************************************************
01630       *                               Full tensor: inputData(C,P,D,D)                              *
01631       **********************************************************************************************/
01632     for (int i=0; i<in_c_f_p_d.size(); i++) {
01633       in_c_f_p_d[i] = Teuchos::ScalarTraits<double>::random();
01634     }
01635     for (int ic=0; ic < c; ic++) {
01636       for (int ip=0; ip < p; ip++) {
01637         for (int i=0; i<d1*d1; i++) {
01638           (&data_c_p_d_d(ic, ip, 0, 0))[i] = Teuchos::ScalarTraits<double>::random();
01639         }
01640         RealSpaceTools<double>::inverse(&datainv_c_p_d_d(ic, ip, 0, 0), &data_c_p_d_d(ic, ip, 0, 0), d1);
01641       }
01642     }
01643     for (int ic=0; ic < c; ic++) {
01644       for (int ip=0; ip < 1; ip++) {
01645         for (int i=0; i<d1*d1; i++) {
01646           (&data_c_1_d_d(ic, 0, 0, 0))[i] = Teuchos::ScalarTraits<double>::random();
01647         }
01648         RealSpaceTools<double>::inverse(&datainv_c_1_d_d(ic, 0, 0, 0), &data_c_1_d_d(ic, 0, 0, 0), d1);
01649       }
01650     }
01651     // Tensor values vary by point: test "N" and "T" options (no transpose/transpose)
01652     art::matvecProductDataField<double>(out_c_f_p_d, data_c_p_d_d, in_c_f_p_d);
01653     art::matvecProductDataField<double>(outi_c_f_p_d, datainv_c_p_d_d, out_c_f_p_d);
01654     rst::subtract(&outi_c_f_p_d[0], &in_c_f_p_d[0], outi_c_f_p_d.size());
01655     if (rst::vectorNorm(&outi_c_f_p_d[0], outi_c_f_p_d.size(), NORM_ONE) > zero) {
01656       *outStream << "\n\nINCORRECT matvecProductDataField (7): check matrix inverse property\n\n";
01657       errorFlag = -1000;
01658     }
01659     art::matvecProductDataField<double>(out_c_f_p_d, data_c_p_d_d, in_c_f_p_d, 't');
01660     art::matvecProductDataField<double>(outi_c_f_p_d, datainv_c_p_d_d, out_c_f_p_d, 't');
01661     rst::subtract(&outi_c_f_p_d[0], &in_c_f_p_d[0], outi_c_f_p_d.size());
01662     if (rst::vectorNorm(&outi_c_f_p_d[0], outi_c_f_p_d.size(), NORM_ONE) > zero) {
01663       *outStream << "\n\nINCORRECT matvecProductDataField (8): check matrix inverse property, w/ double transpose\n\n";
01664       errorFlag = -1000;
01665     }
01666     // Tensor values do not vary by point: test "N" and "T" options (no transpose/transpose)
01667     art::matvecProductDataField<double>(out_c_f_p_d, data_c_1_d_d, in_c_f_p_d);
01668     art::matvecProductDataField<double>(outi_c_f_p_d, datainv_c_1_d_d, out_c_f_p_d);
01669     rst::subtract(&outi_c_f_p_d[0], &in_c_f_p_d[0], outi_c_f_p_d.size());
01670     if (rst::vectorNorm(&outi_c_f_p_d[0], outi_c_f_p_d.size(), NORM_ONE) > zero) {
01671       *outStream << "\n\nINCORRECT matvecProductDataField (9): check matrix inverse property\n\n";
01672       errorFlag = -1000;
01673     }
01674     art::matvecProductDataField<double>(out_c_f_p_d, data_c_1_d_d, in_c_f_p_d, 't');
01675     art::matvecProductDataField<double>(outi_c_f_p_d, datainv_c_1_d_d, out_c_f_p_d, 't');
01676     rst::subtract(&outi_c_f_p_d[0], &in_c_f_p_d[0], outi_c_f_p_d.size());
01677     if (rst::vectorNorm(&outi_c_f_p_d[0], outi_c_f_p_d.size(), NORM_ONE) > zero) {
01678       *outStream << "\n\nINCORRECT matvecProductDataField (10): check matrix inverse property, w/ double transpose\n\n";
01679       errorFlag = -1000;
01680     }
01681     /***********************************************************************************************
01682       *             Full tensor: inputData(C,P,D,D) test inverse transpose                         *
01683       **********************************************************************************************/
01684     for (int i=0; i<in_c_f_p_d.size(); i++) {
01685       in_c_f_p_d[i] = Teuchos::ScalarTraits<double>::random();
01686     }
01687     for (int ic=0; ic < c; ic++) {
01688       for (int ip=0; ip < p; ip++) {
01689         for (int i=0; i<d1*d1; i++) {
01690           (&data_c_p_d_d(ic, ip, 0, 0))[i] = Teuchos::ScalarTraits<double>::random();
01691         }
01692         RealSpaceTools<double>::inverse(&datainv_c_p_d_d(ic, ip, 0, 0), &data_c_p_d_d(ic, ip, 0, 0), d1);
01693         RealSpaceTools<double>::transpose(&datainvtrn_c_p_d_d(ic, ip, 0, 0), &datainv_c_p_d_d(ic, ip, 0, 0), d1);
01694       }
01695     }
01696     for (int ic=0; ic < c; ic++) {
01697       for (int ip=0; ip < 1; ip++) {
01698         for (int i=0; i<d1*d1; i++) {
01699           (&data_c_1_d_d(ic, 0, 0, 0))[i] = Teuchos::ScalarTraits<double>::random();
01700         }
01701         RealSpaceTools<double>::inverse(&datainv_c_1_d_d(ic, 0, 0, 0), &data_c_1_d_d(ic, 0, 0, 0), d1);
01702         RealSpaceTools<double>::transpose(&datainvtrn_c_1_d_d(ic, 0, 0, 0), &datainv_c_1_d_d(ic, 0, 0, 0), d1);
01703       }
01704     }
01705     // Tensor values vary by point:
01706     art::matvecProductDataField<double>(out_c_f_p_d, data_c_p_d_d, in_c_f_p_d);
01707     art::matvecProductDataField<double>(outi_c_f_p_d, datainvtrn_c_p_d_d, out_c_f_p_d, 't');
01708     rst::subtract(&outi_c_f_p_d[0], &in_c_f_p_d[0], outi_c_f_p_d.size());
01709     if (rst::vectorNorm(&outi_c_f_p_d[0], outi_c_f_p_d.size(), NORM_ONE) > zero) {
01710       *outStream << "\n\nINCORRECT matvecProductDataField (11): check matrix inverse transpose property\n\n";
01711       errorFlag = -1000;
01712     }
01713     // Tensor values do not vary by point
01714     art::matvecProductDataField<double>(out_c_f_p_d, data_c_1_d_d, in_c_f_p_d);
01715     art::matvecProductDataField<double>(outi_c_f_p_d, datainvtrn_c_1_d_d, out_c_f_p_d, 't');
01716     rst::subtract(&outi_c_f_p_d[0], &in_c_f_p_d[0], outi_c_f_p_d.size());
01717     if (rst::vectorNorm(&outi_c_f_p_d[0], outi_c_f_p_d.size(), NORM_ONE) > zero) {
01718       *outStream << "\n\nINCORRECT matvecProductDataField (12): check matrix inverse transpose property\n\n";
01719       errorFlag = -1000;
01720     }
01721   }// test 5.c scope
01722   
01723   
01724   *outStream \
01725     << "\n"
01726     << "===============================================================================\n"\
01727     << "| TEST 5.d: matvecProductDataField random tests: branch inputFields(F,P,D)    |\n"\
01728     << "===============================================================================\n";
01729   /*
01730    *  d1 is the spatial dimension and should be 1, 2 or 3. If d1>3, the RealSpaceTools function 'inverse' will fail
01731    */
01732   {// test 5.d scope
01733     int c=5, p=9, f=7, d1=3;
01734     double zero = INTREPID_TOL*10000.0;
01735     
01736     FieldContainer<double> in_f_p_d(f, p, d1);
01737     FieldContainer<double> in_c_f_p_d(c, f, p, d1);
01738     FieldContainer<double> data_c_p(c, p);
01739     FieldContainer<double> datainv_c_p(c, p);
01740     FieldContainer<double> data_c_1(c, 1);
01741     FieldContainer<double> datainv_c_1(c, 1);
01742     FieldContainer<double> data_c_p_d(c, p, d1);
01743     FieldContainer<double> datainv_c_p_d(c, p, d1);
01744     FieldContainer<double> data_c_1_d(c, 1, d1);
01745     FieldContainer<double> datainv_c_1_d(c, 1, d1);
01746     FieldContainer<double> data_c_p_d_d(c, p, d1, d1);
01747     FieldContainer<double> datainv_c_p_d_d(c, p, d1, d1);
01748     FieldContainer<double> datainvtrn_c_p_d_d(c, p, d1, d1);
01749     FieldContainer<double> data_c_1_d_d(c, 1, d1, d1);
01750     FieldContainer<double> datainv_c_1_d_d(c, 1, d1, d1);
01751     FieldContainer<double> datainvtrn_c_1_d_d(c, 1, d1, d1);
01752     FieldContainer<double> data_c_p_one(c, p);
01753     FieldContainer<double> data_c_1_one(c, 1);
01754     FieldContainer<double> out_c_f_p_d(c, f, p, d1);
01755     FieldContainer<double> outi_c_f_p_d(c, f, p, d1);
01756     /***********************************************************************************************
01757       *                          Constant diagonal tensor: inputData(C,P)                          *
01758       **********************************************************************************************/
01759     for (int i=0; i<in_f_p_d.size(); i++) {
01760       in_f_p_d[i] = Teuchos::ScalarTraits<double>::random();
01761     }
01762     for (int i=0; i<data_c_p.size(); i++) {
01763       data_c_p[i] = Teuchos::ScalarTraits<double>::random();
01764       datainv_c_p[i] = 1.0 / data_c_p[i];
01765       data_c_p_one[i] = 1.0;
01766     }
01767     for (int i=0; i<data_c_1.size(); i++) {
01768       data_c_1[i] = Teuchos::ScalarTraits<double>::random();
01769       datainv_c_1[i] = 1.0 / data_c_1[i];
01770     }
01771     // Tensor values vary by point
01772     art::matvecProductDataField<double>(in_c_f_p_d, data_c_p_one, in_f_p_d);
01773     art::matvecProductDataField<double>(out_c_f_p_d, data_c_p, in_f_p_d);
01774     art::matvecProductDataField<double>(outi_c_f_p_d, datainv_c_p, out_c_f_p_d);
01775     rst::subtract(&outi_c_f_p_d[0], &in_c_f_p_d[0], outi_c_f_p_d.size());
01776     if (rst::vectorNorm(&outi_c_f_p_d[0], outi_c_f_p_d.size(), NORM_ONE) > zero) {
01777       *outStream << "\n\nINCORRECT matvecProductDataField (13): check scalar inverse property\n\n";
01778       errorFlag = -1000;
01779     }
01780     // Tensor values do not vary by point
01781     art::matvecProductDataField<double>(in_c_f_p_d, data_c_p_one, in_f_p_d);
01782     art::matvecProductDataField<double>(out_c_f_p_d, data_c_1, in_f_p_d);
01783     art::matvecProductDataField<double>(outi_c_f_p_d, datainv_c_1, out_c_f_p_d);
01784     rst::subtract(&outi_c_f_p_d[0], &in_c_f_p_d[0], outi_c_f_p_d.size());
01785     if (rst::vectorNorm(&outi_c_f_p_d[0], outi_c_f_p_d.size(), NORM_ONE) > zero) {
01786       *outStream << "\n\nINCORRECT matvecProductDataField (14): check scalar inverse property\n\n";
01787       errorFlag = -1000;
01788     }
01789     /***********************************************************************************************
01790       *                       Non-constant diagonal tensor: inputData(C,P,D)                       *
01791       **********************************************************************************************/
01792     for (int i=0; i<in_f_p_d.size(); i++) {
01793       in_f_p_d[i] = Teuchos::ScalarTraits<double>::random();
01794     }
01795     for (int i=0; i<data_c_p_d.size(); i++) {
01796       data_c_p_d[i] = Teuchos::ScalarTraits<double>::random();
01797       datainv_c_p_d[i] = 1.0 / data_c_p_d[i];
01798     }
01799     for (int i=0; i<data_c_1_d.size(); i++) {
01800       data_c_1_d[i] = Teuchos::ScalarTraits<double>::random();
01801       datainv_c_1_d[i] = 1.0 / data_c_1_d[i];
01802     }
01803     // Tensor values vary by point:
01804     art::matvecProductDataField<double>(in_c_f_p_d, data_c_p_one, in_f_p_d);
01805     art::matvecProductDataField<double>(out_c_f_p_d, data_c_p_d, in_f_p_d);
01806     art::matvecProductDataField<double>(outi_c_f_p_d, datainv_c_p_d, out_c_f_p_d);
01807     rst::subtract(&outi_c_f_p_d[0], &in_c_f_p_d[0], outi_c_f_p_d.size());
01808     if (rst::vectorNorm(&outi_c_f_p_d[0], outi_c_f_p_d.size(), NORM_ONE) > zero) {
01809       *outStream << "\n\nINCORRECT matvecProductDataField (15): check scalar inverse property\n\n";
01810       errorFlag = -1000;
01811     }
01812     // Tensor values do not vary by point:
01813     art::matvecProductDataField<double>(in_c_f_p_d, data_c_p_one, in_f_p_d);
01814     art::matvecProductDataField<double>(out_c_f_p_d, data_c_1_d, in_f_p_d);
01815     art::matvecProductDataField<double>(outi_c_f_p_d, datainv_c_1_d, out_c_f_p_d);
01816     rst::subtract(&outi_c_f_p_d[0], &in_c_f_p_d[0], outi_c_f_p_d.size());
01817     if (rst::vectorNorm(&outi_c_f_p_d[0], outi_c_f_p_d.size(), NORM_ONE) > zero) {
01818       *outStream << "\n\nINCORRECT matvecProductDataField (16): check scalar inverse property\n\n";
01819       errorFlag = -1000;
01820     }
01821     /***********************************************************************************************
01822       *                              Full tensor: inputData(C,P,D,D)                               *
01823       **********************************************************************************************/
01824     for (int i=0; i<in_f_p_d.size(); i++) {
01825       in_f_p_d[i] = Teuchos::ScalarTraits<double>::random();
01826     }
01827     for (int ic=0; ic < c; ic++) {
01828       for (int ip=0; ip < p; ip++) {
01829         for (int i=0; i<d1*d1; i++) {
01830           (&data_c_p_d_d(ic, ip, 0, 0))[i] = Teuchos::ScalarTraits<double>::random();
01831         }
01832         RealSpaceTools<double>::inverse(&datainv_c_p_d_d(ic, ip, 0, 0), &data_c_p_d_d(ic, ip, 0, 0), d1);
01833       }
01834     }
01835     for (int ic=0; ic < c; ic++) {
01836       for (int ip=0; ip < 1; ip++) {
01837         for (int i=0; i<d1*d1; i++) {
01838           (&data_c_1_d_d(ic, 0, 0, 0))[i] = Teuchos::ScalarTraits<double>::random();
01839         }
01840         RealSpaceTools<double>::inverse(&datainv_c_1_d_d(ic, 0, 0, 0), &data_c_1_d_d(ic, 0, 0, 0), d1);
01841       }
01842     }
01843     // Tensor values vary by point: test "N" and "T" (no-transpose/transpose) options
01844     art::matvecProductDataField<double>(in_c_f_p_d, data_c_p_one, in_f_p_d);
01845     art::matvecProductDataField<double>(out_c_f_p_d, data_c_p_d_d, in_f_p_d);
01846     art::matvecProductDataField<double>(outi_c_f_p_d, datainv_c_p_d_d, out_c_f_p_d);
01847     rst::subtract(&outi_c_f_p_d[0], &in_c_f_p_d[0], outi_c_f_p_d.size());
01848     if (rst::vectorNorm(&outi_c_f_p_d[0], outi_c_f_p_d.size(), NORM_ONE) > zero) {
01849       *outStream << "\n\nINCORRECT matvecProductDataField (17): check matrix inverse property\n\n";
01850       errorFlag = -1000;
01851     }
01852     art::matvecProductDataField<double>(in_c_f_p_d, data_c_p_one, in_f_p_d);
01853     art::matvecProductDataField<double>(out_c_f_p_d, data_c_p_d_d, in_f_p_d, 't');
01854     art::matvecProductDataField<double>(outi_c_f_p_d, datainv_c_p_d_d, out_c_f_p_d, 't');
01855     rst::subtract(&outi_c_f_p_d[0], &in_c_f_p_d[0], outi_c_f_p_d.size());
01856     if (rst::vectorNorm(&outi_c_f_p_d[0], outi_c_f_p_d.size(), NORM_ONE) > zero) {
01857       *outStream << "\n\nINCORRECT matvecProductDataField (18): check matrix inverse property, w/ double transpose\n\n";
01858       errorFlag = -1000;
01859     }
01860     // Tensor values do not vary by point: test "N" and "T" (no-transpose/transpose) options
01861     art::matvecProductDataField<double>(in_c_f_p_d, data_c_p_one, in_f_p_d);
01862     art::matvecProductDataField<double>(out_c_f_p_d, data_c_1_d_d, in_f_p_d);
01863     art::matvecProductDataField<double>(outi_c_f_p_d, datainv_c_1_d_d, out_c_f_p_d);
01864     rst::subtract(&outi_c_f_p_d[0], &in_c_f_p_d[0], outi_c_f_p_d.size());
01865     if (rst::vectorNorm(&outi_c_f_p_d[0], outi_c_f_p_d.size(), NORM_ONE) > zero) {
01866       *outStream << "\n\nINCORRECT matvecProductDataField (19): check matrix inverse property\n\n";
01867       errorFlag = -1000;
01868     }
01869     art::matvecProductDataField<double>(in_c_f_p_d, data_c_p_one, in_f_p_d);
01870     art::matvecProductDataField<double>(out_c_f_p_d, data_c_1_d_d, in_f_p_d, 't');
01871     art::matvecProductDataField<double>(outi_c_f_p_d, datainv_c_1_d_d, out_c_f_p_d, 't');
01872     rst::subtract(&outi_c_f_p_d[0], &in_c_f_p_d[0], outi_c_f_p_d.size());
01873     if (rst::vectorNorm(&outi_c_f_p_d[0], outi_c_f_p_d.size(), NORM_ONE) > zero) {
01874       *outStream << "\n\nINCORRECT matvecProductDataField (20): check matrix inverse property, w/ double transpose\n\n";
01875       errorFlag = -1000;
01876     }
01877     /***********************************************************************************************
01878       *             Full tensor: inputData(C,P,D,D) test inverse transpose                         *
01879       **********************************************************************************************/
01880     for (int i=0; i<in_f_p_d.size(); i++) {
01881       in_f_p_d[i] = Teuchos::ScalarTraits<double>::random();
01882     }
01883     for (int ic=0; ic < c; ic++) {
01884       for (int ip=0; ip < p; ip++) {
01885         for (int i=0; i<d1*d1; i++) {
01886           (&data_c_p_d_d(ic, ip, 0, 0))[i] = Teuchos::ScalarTraits<double>::random();
01887         }
01888         RealSpaceTools<double>::inverse(&datainv_c_p_d_d(ic, ip, 0, 0), &data_c_p_d_d(ic, ip, 0, 0), d1);
01889         RealSpaceTools<double>::transpose(&datainvtrn_c_p_d_d(ic, ip, 0, 0), &datainv_c_p_d_d(ic, ip, 0, 0), d1);
01890       }
01891     }
01892     for (int ic=0; ic < c; ic++) {
01893       for (int ip=0; ip < 1; ip++) {
01894         for (int i=0; i<d1*d1; i++) {
01895           (&data_c_1_d_d(ic, 0, 0, 0))[i] = Teuchos::ScalarTraits<double>::random();
01896         }
01897         RealSpaceTools<double>::inverse(&datainv_c_1_d_d(ic, 0, 0, 0), &data_c_1_d_d(ic, 0, 0, 0), d1);
01898         RealSpaceTools<double>::transpose(&datainvtrn_c_1_d_d(ic, 0, 0, 0), &datainv_c_1_d_d(ic, 0, 0, 0), d1);
01899       }
01900     }
01901     // Tensor values vary by point:
01902     art::matvecProductDataField<double>(in_c_f_p_d, data_c_p_one, in_f_p_d);
01903     art::matvecProductDataField<double>(out_c_f_p_d, data_c_p_d_d, in_f_p_d);
01904     art::matvecProductDataField<double>(outi_c_f_p_d, datainvtrn_c_p_d_d, out_c_f_p_d, 't');
01905     rst::subtract(&outi_c_f_p_d[0], &in_c_f_p_d[0], outi_c_f_p_d.size());
01906     if (rst::vectorNorm(&outi_c_f_p_d[0], outi_c_f_p_d.size(), NORM_ONE) > zero) {
01907       *outStream << "\n\nINCORRECT matvecProductDataField (21): check matrix inverse transpose property\n\n";
01908       errorFlag = -1000;
01909     }
01910     // Tensor values do not vary by point:
01911     art::matvecProductDataField<double>(in_c_f_p_d, data_c_p_one, in_f_p_d);
01912     art::matvecProductDataField<double>(out_c_f_p_d, data_c_1_d_d, in_f_p_d);
01913     art::matvecProductDataField<double>(outi_c_f_p_d, datainvtrn_c_1_d_d, out_c_f_p_d, 't');
01914     rst::subtract(&outi_c_f_p_d[0], &in_c_f_p_d[0], outi_c_f_p_d.size());
01915     if (rst::vectorNorm(&outi_c_f_p_d[0], outi_c_f_p_d.size(), NORM_ONE) > zero) {
01916       *outStream << "\n\nINCORRECT matvecProductDataField (22): check matrix inverse transpose property\n\n";
01917       errorFlag = -1000;
01918     }    
01919   }// test 5.d scope
01920   
01921   
01922   
01923   *outStream \
01924     << "\n"
01925     << "===============================================================================\n"\
01926     << "| TEST 6.a: matvecProductDataData operations: (C,P,D,D) and (C,P,D)           |\n"\
01927     << "===============================================================================\n";
01928   /*
01929    *                     inputMat                     inputVecFields       outFields
01930    *      1  1  1   0  0  0    1  1  1    0  0  0     0   1  -1  -1      0   0  -3   0  
01931    *     -1  2 -1  -1 -2 -3   -1  2 -1   -1 -2 -3     0   1  -1   1      0  -6   0   2
01932    *      1  2  3  -2  6 -4    1  2  3   -2  6 -4     0   1  -1  -1      0   0  -6  12
01933    */
01934   
01935   // (C,P,D,D)
01936   inputMat.resize(4,1,3,3);
01937   // cell 0
01938   inputMat(0,0,0,0) = 1.0;  inputMat(0,0,0,1) = 1.0;  inputMat(0,0,0,2) = 1.0;
01939   inputMat(0,0,1,0) =-1.0;  inputMat(0,0,1,1) = 2.0;  inputMat(0,0,1,2) =-1.0;
01940   inputMat(0,0,2,0) = 1.0;  inputMat(0,0,2,1) = 2.0;  inputMat(0,0,2,2) = 3.0;
01941   // cell 1
01942   inputMat(1,0,0,0) = 0.0;  inputMat(1,0,0,1) = 0.0;  inputMat(1,0,0,2) = 0.0;
01943   inputMat(1,0,1,0) =-1.0;  inputMat(1,0,1,1) =-2.0;  inputMat(1,0,1,2) =-3.0;
01944   inputMat(1,0,2,0) =-2.0;  inputMat(1,0,2,1) = 6.0;  inputMat(1,0,2,2) =-4.0;
01945   // cell 2
01946   inputMat(2,0,0,0) = 1.0;  inputMat(2,0,0,1) = 1.0;  inputMat(2,0,0,2) = 1.0;
01947   inputMat(2,0,1,0) =-1.0;  inputMat(2,0,1,1) = 2.0;  inputMat(2,0,1,2) =-1.0;
01948   inputMat(2,0,2,0) = 1.0;  inputMat(2,0,2,1) = 2.0;  inputMat(2,0,2,2) = 3.0;
01949   // cell 3
01950   inputMat(3,0,0,0) = 0.0;  inputMat(3,0,0,1) = 0.0;  inputMat(3,0,0,2) = 0.0;
01951   inputMat(3,0,1,0) =-1.0;  inputMat(3,0,1,1) =-2.0;  inputMat(3,0,1,2) =-3.0;
01952   inputMat(3,0,2,0) =-2.0;  inputMat(3,0,2,1) = 6.0;  inputMat(3,0,2,2) =-4.0;
01953   
01954   // (C,P,D)
01955   inputVecFields.resize(4,1,3);
01956   inputVecFields(0,0,0) = 0.0;  inputVecFields(0,0,1) = 0.0;  inputVecFields(0,0,2) = 0.0;
01957   inputVecFields(1,0,0) = 1.0;  inputVecFields(1,0,1) = 1.0;  inputVecFields(1,0,2) = 1.0;
01958   inputVecFields(2,0,0) =-1.0;  inputVecFields(2,0,1) =-1.0;  inputVecFields(2,0,2) =-1.0;
01959   inputVecFields(3,0,0) =-1.0;  inputVecFields(3,0,1) = 1.0;  inputVecFields(3,0,2) =-1.0;
01960   
01961   // (C,P,D) - true
01962   outFieldsCorrect.resize(4,1,3);
01963   outFieldsCorrect(0,0,0) = 0.0;  outFieldsCorrect(0,0,1) = 0.0;  outFieldsCorrect(0,0,2) = 0.0;
01964   outFieldsCorrect(1,0,0) = 0.0;  outFieldsCorrect(1,0,1) =-6.0;  outFieldsCorrect(1,0,2) = 0.0;
01965   outFieldsCorrect(2,0,0) =-3.0;  outFieldsCorrect(2,0,1) = 0.0;  outFieldsCorrect(2,0,2) =-6.0;
01966   outFieldsCorrect(3,0,0) = 0.0;  outFieldsCorrect(3,0,1) = 2.0;  outFieldsCorrect(3,0,2) = 12.0;
01967 
01968   // (C,P,D)
01969   outFields.resize(4,1,3);
01970   art::matvecProductDataData<double>(outFields, inputMat, inputVecFields);
01971   
01972   // test loop
01973   for(int cell = 0; cell < outFields.dimension(0); cell++){
01974     for(int point = 0; point < outFields.dimension(1); point++){
01975       for(int row = 0; row < outFields.dimension(2); row++){
01976         if(outFields(cell, point, row) != outFieldsCorrect(cell, point, row)) {
01977           *outStream << "\n\nINCORRECT matvecProductDataData (1): \n value at multi-index ("
01978           << cell << "," << point << "," << row << ") = " 
01979           << outFields(cell, point, row) << " but correct value is " 
01980           << outFieldsCorrect(cell, point, row) <<"\n";
01981           errorFlag++; 
01982         }
01983       }//row
01984     }// point
01985   }// cell
01986   
01987   
01988   *outStream \
01989     << "\n"
01990     << "===============================================================================\n"\
01991     << "| TEST 6.b: matvecProductDataData operations: (C,P,D,D) and (P,D)             |\n"\
01992     << "===============================================================================\n";
01993   /*
01994    *                     inputMat                     inputVecFields       outFields
01995    *      1  1  1   0  0  0    1  1  1    0  0  0     0   1  -1  -1      0   0  -3   0  
01996    *     -1  2 -1  -1 -2 -3   -1  2 -1   -1 -2 -3     0   1  -1   1      0  -6   0   2
01997    *      1  2  3  -2  6 -4    1  2  3   -2  6 -4     0   1  -1  -1      0   0  -6  12
01998    */
01999   // (C,P,D,D)
02000   inputMat.resize(1,4,3,3);
02001   // point 0
02002   inputMat(0,0,0,0) = 1.0;  inputMat(0,0,0,1) = 1.0;  inputMat(0,0,0,2) = 1.0;
02003   inputMat(0,0,1,0) =-1.0;  inputMat(0,0,1,1) = 2.0;  inputMat(0,0,1,2) =-1.0;
02004   inputMat(0,0,2,0) = 1.0;  inputMat(0,0,2,1) = 2.0;  inputMat(0,0,2,2) = 3.0;
02005   // point 1
02006   inputMat(0,1,0,0) = 0.0;  inputMat(0,1,0,1) = 0.0;  inputMat(0,1,0,2) = 0.0;
02007   inputMat(0,1,1,0) =-1.0;  inputMat(0,1,1,1) =-2.0;  inputMat(0,1,1,2) =-3.0;
02008   inputMat(0,1,2,0) =-2.0;  inputMat(0,1,2,1) = 6.0;  inputMat(0,1,2,2) =-4.0;
02009   // point 2
02010   inputMat(0,2,0,0) = 1.0;  inputMat(0,2,0,1) = 1.0;  inputMat(0,2,0,2) = 1.0;
02011   inputMat(0,2,1,0) =-1.0;  inputMat(0,2,1,1) = 2.0;  inputMat(0,2,1,2) =-1.0;
02012   inputMat(0,2,2,0) = 1.0;  inputMat(0,2,2,1) = 2.0;  inputMat(0,2,2,2) = 3.0;
02013   // point 3
02014   inputMat(0,3,0,0) = 0.0;  inputMat(0,3,0,1) = 0.0;  inputMat(0,3,0,2) = 0.0;
02015   inputMat(0,3,1,0) =-1.0;  inputMat(0,3,1,1) =-2.0;  inputMat(0,3,1,2) =-3.0;
02016   inputMat(0,3,2,0) =-2.0;  inputMat(0,3,2,1) = 6.0;  inputMat(0,3,2,2) =-4.0;
02017   
02018   // (P,D)
02019   inputVecFields.resize(4,3);
02020   // 
02021   inputVecFields(0,0) = 0.0;  inputVecFields(0,1) = 0.0;  inputVecFields(0,2) = 0.0;
02022   inputVecFields(1,0) = 1.0;  inputVecFields(1,1) = 1.0;  inputVecFields(1,2) = 1.0;
02023   inputVecFields(2,0) =-1.0;  inputVecFields(2,1) =-1.0;  inputVecFields(2,2) =-1.0;
02024   inputVecFields(3,0) =-1.0;  inputVecFields(3,1) = 1.0;  inputVecFields(3,2) =-1.0;
02025   
02026   // (C,P,D) - true
02027   outFieldsCorrect.resize(1,4,3);
02028   outFieldsCorrect(0,0,0) = 0.0;  outFieldsCorrect(0,0,1) = 0.0;  outFieldsCorrect(0,0,2) = 0.0;
02029   outFieldsCorrect(0,1,0) = 0.0;  outFieldsCorrect(0,1,1) =-6.0;  outFieldsCorrect(0,1,2) = 0.0;
02030   outFieldsCorrect(0,2,0) =-3.0;  outFieldsCorrect(0,2,1) = 0.0;  outFieldsCorrect(0,2,2) =-6.0;
02031   outFieldsCorrect(0,3,0) = 0.0;  outFieldsCorrect(0,3,1) = 2.0;  outFieldsCorrect(0,3,2) = 12.0;
02032   
02033   // (C,P,D)
02034   outFields.resize(1,4,3);
02035   art::matvecProductDataData<double>(outFields, inputMat, inputVecFields);
02036   
02037   // test loop
02038   for(int cell = 0; cell < outFields.dimension(0); cell++){
02039     for(int point = 0; point < outFields.dimension(1); point++){
02040       for(int row = 0; row < outFields.dimension(2); row++){
02041         if(outFields(cell, point, row) != outFieldsCorrect(cell, point, row)) {
02042           *outStream << "\n\nINCORRECT matvecProductDataData (2): \n value at multi-index ("
02043           << cell << "," << point << "," << row << ") = " 
02044           << outFields(cell, point, row) << " but correct value is " 
02045           << outFieldsCorrect(cell, point, row) <<"\n";
02046           errorFlag++; 
02047         }
02048       }//row
02049     }// point
02050   }// cell
02051   
02052 
02053   
02054   *outStream \
02055     << "\n"
02056     << "===============================================================================\n"\
02057     << "| TEST 6.c: matvecProductDataData random tests: branch inputDataRight(C,P,D)  |\n"\
02058     << "===============================================================================\n";
02059   /*
02060    * Test derived from Test 5.c
02061    */
02062   {// test 6.c scope
02063     int c=5, p=9, d1=3;
02064     double zero = INTREPID_TOL*10000.0;
02065     
02066     FieldContainer<double> in_c_p_d(c, p, d1);
02067     FieldContainer<double> out_c_p_d(c, p, d1);
02068     FieldContainer<double> outi_c_p_d(c, p, d1);
02069     
02070     FieldContainer<double> data_c_p(c, p);
02071     FieldContainer<double> datainv_c_p(c, p);
02072     FieldContainer<double> data_c_1(c, 1);
02073     FieldContainer<double> datainv_c_1(c, 1);
02074     FieldContainer<double> data_c_p_d(c, p, d1);
02075     FieldContainer<double> datainv_c_p_d(c, p, d1);
02076     FieldContainer<double> data_c_1_d(c, 1, d1);
02077     FieldContainer<double> datainv_c_1_d(c, 1, d1);
02078     FieldContainer<double> data_c_p_d_d(c, p, d1, d1);
02079     FieldContainer<double> datainv_c_p_d_d(c, p, d1, d1);
02080     FieldContainer<double> datainvtrn_c_p_d_d(c, p, d1, d1);
02081     FieldContainer<double> data_c_1_d_d(c, 1, d1, d1);
02082     FieldContainer<double> datainv_c_1_d_d(c, 1, d1, d1);
02083     FieldContainer<double> datainvtrn_c_1_d_d(c, 1, d1, d1);
02084     /***********************************************************************************************
02085       *                          Constant diagonal tensor: inputDataLeft(C,P)                      *
02086       **********************************************************************************************/
02087     for (int i=0; i<in_c_p_d.size(); i++) {
02088       in_c_p_d[i] = Teuchos::ScalarTraits<double>::random();
02089     }
02090     for (int i=0; i<data_c_p.size(); i++) {
02091       data_c_p[i] = Teuchos::ScalarTraits<double>::random();
02092       datainv_c_p[i] = 1.0 / data_c_p[i];
02093     }
02094     for (int i=0; i<data_c_1.size(); i++) {
02095       data_c_1[i] = Teuchos::ScalarTraits<double>::random();
02096       datainv_c_1[i] = 1.0 / data_c_1[i];
02097     }
02098     // Tensor values vary by point:
02099     art::matvecProductDataData<double>(out_c_p_d, data_c_p, in_c_p_d);
02100     art::matvecProductDataData<double>(out_c_p_d, datainv_c_p, out_c_p_d);
02101     rst::subtract(&out_c_p_d[0], &in_c_p_d[0], out_c_p_d.size());
02102     if (rst::vectorNorm(&out_c_p_d[0], out_c_p_d.size(), NORM_ONE) > zero) {
02103       *outStream << "\n\nINCORRECT matvecProductDataData (3): check scalar inverse property\n\n";
02104       errorFlag = -1000;
02105     }
02106     // Tensor values do not vary by point:
02107     art::matvecProductDataData<double>(out_c_p_d, data_c_1, in_c_p_d);
02108     art::matvecProductDataData<double>(out_c_p_d, datainv_c_1, out_c_p_d);
02109     rst::subtract(&out_c_p_d[0], &in_c_p_d[0], out_c_p_d.size());
02110     if (rst::vectorNorm(&out_c_p_d[0], out_c_p_d.size(), NORM_ONE) > zero) {
02111       *outStream << "\n\nINCORRECT matvecProductDataData (4): check scalar inverse property\n\n";
02112       errorFlag = -1000;
02113     }
02114     /***********************************************************************************************
02115       *                     Non-onstant diagonal tensor: inputDataLeft(C,P,D)                      *
02116       **********************************************************************************************/
02117     for (int i=0; i<in_c_p_d.size(); i++) {
02118       in_c_p_d[i] = Teuchos::ScalarTraits<double>::random();
02119     }
02120     for (int i=0; i<data_c_p_d.size(); i++) {
02121       data_c_p_d[i] = Teuchos::ScalarTraits<double>::random();
02122       datainv_c_p_d[i] = 1.0 / data_c_p_d[i];
02123     }
02124     for (int i=0; i<data_c_1_d.size(); i++) {
02125       data_c_1_d[i] = Teuchos::ScalarTraits<double>::random();
02126       datainv_c_1_d[i] = 1.0 / data_c_1_d[i];
02127     }
02128     // Tensor values vary by point:
02129     art::matvecProductDataData<double>(out_c_p_d, data_c_p_d, in_c_p_d);
02130     art::matvecProductDataData<double>(out_c_p_d, datainv_c_p_d, out_c_p_d);
02131     rst::subtract(&out_c_p_d[0], &in_c_p_d[0], out_c_p_d.size());
02132     if (rst::vectorNorm(&out_c_p_d[0], out_c_p_d.size(), NORM_ONE) > zero) {
02133       *outStream << "\n\nINCORRECT matvecProductDataData (5): check scalar inverse property\n\n";
02134       errorFlag = -1000;
02135     }
02136     // Tensor values do not vary by point
02137     art::matvecProductDataData<double>(out_c_p_d, data_c_1_d, in_c_p_d);
02138     art::matvecProductDataData<double>(out_c_p_d, datainv_c_1_d, out_c_p_d);
02139     rst::subtract(&out_c_p_d[0], &in_c_p_d[0], out_c_p_d.size());
02140     if (rst::vectorNorm(&out_c_p_d[0], out_c_p_d.size(), NORM_ONE) > zero) {
02141       *outStream << "\n\nINCORRECT matvecProductDataData (6): check scalar inverse property\n\n";
02142       errorFlag = -1000;
02143     }
02144     /***********************************************************************************************
02145       *                            Full tensor: inputDataLeft(C,P,D,D)                             *
02146       **********************************************************************************************/
02147     for (int i=0; i<in_c_p_d.size(); i++) {
02148       in_c_p_d[i] = Teuchos::ScalarTraits<double>::random();
02149     }
02150     for (int ic=0; ic < c; ic++) {
02151       for (int ip=0; ip < p; ip++) {
02152         for (int i=0; i<d1*d1; i++) {
02153           (&data_c_p_d_d(ic, ip, 0, 0))[i] = Teuchos::ScalarTraits<double>::random();
02154         }
02155         RealSpaceTools<double>::inverse(&datainv_c_p_d_d(ic, ip, 0, 0), &data_c_p_d_d(ic, ip, 0, 0), d1);
02156       }
02157     }
02158     for (int ic=0; ic < c; ic++) {
02159       for (int ip=0; ip < 1; ip++) {
02160         for (int i=0; i<d1*d1; i++) {
02161           (&data_c_1_d_d(ic, 0, 0, 0))[i] = Teuchos::ScalarTraits<double>::random();
02162         }
02163         RealSpaceTools<double>::inverse(&datainv_c_1_d_d(ic, 0, 0, 0), &data_c_1_d_d(ic, 0, 0, 0), d1);
02164       }
02165     }
02166     // Tensor values vary by point: test "N" and "T" options (no transpose/transpose)
02167     art::matvecProductDataData<double>(out_c_p_d, data_c_p_d_d, in_c_p_d);
02168     art::matvecProductDataData<double>(outi_c_p_d, datainv_c_p_d_d, out_c_p_d);
02169     rst::subtract(&outi_c_p_d[0], &in_c_p_d[0], outi_c_p_d.size());
02170     if (rst::vectorNorm(&outi_c_p_d[0], outi_c_p_d.size(), NORM_ONE) > zero) {
02171       *outStream << "\n\nINCORRECT matvecProductDataData (7): check matrix inverse property\n\n";
02172       errorFlag = -1000;
02173     }
02174     art::matvecProductDataData<double>(out_c_p_d, data_c_p_d_d, in_c_p_d, 't');
02175     art::matvecProductDataData<double>(outi_c_p_d, datainv_c_p_d_d, out_c_p_d, 't');
02176     rst::subtract(&outi_c_p_d[0], &in_c_p_d[0], outi_c_p_d.size());
02177     if (rst::vectorNorm(&outi_c_p_d[0], outi_c_p_d.size(), NORM_ONE) > zero) {
02178       *outStream << "\n\nINCORRECT matvecProductDataData (8): check matrix inverse property, w/ double transpose\n\n";
02179       errorFlag = -1000;
02180     }
02181     // Tensor values do not vary by point: test "N" and "T" options (no transpose/transpose)
02182     art::matvecProductDataData<double>(out_c_p_d, data_c_1_d_d, in_c_p_d);
02183     art::matvecProductDataData<double>(outi_c_p_d, datainv_c_1_d_d, out_c_p_d);
02184     rst::subtract(&outi_c_p_d[0], &in_c_p_d[0], outi_c_p_d.size());
02185     if (rst::vectorNorm(&outi_c_p_d[0], outi_c_p_d.size(), NORM_ONE) > zero) {
02186       *outStream << "\n\nINCORRECT matvecProductDataData (9): check matrix inverse property\n\n";
02187       errorFlag = -1000;
02188     }
02189     art::matvecProductDataData<double>(out_c_p_d, data_c_1_d_d, in_c_p_d, 't');
02190     art::matvecProductDataData<double>(outi_c_p_d, datainv_c_1_d_d, out_c_p_d, 't');
02191     rst::subtract(&outi_c_p_d[0], &in_c_p_d[0], outi_c_p_d.size());
02192     if (rst::vectorNorm(&outi_c_p_d[0], outi_c_p_d.size(), NORM_ONE) > zero) {
02193       *outStream << "\n\nINCORRECT matvecProductDataData (10): check matrix inverse property, w/ double transpose\n\n";
02194       errorFlag = -1000;
02195     }
02196     /***********************************************************************************************
02197       *             Full tensor: inputDataLeft(C,P,D,D) test inverse transpose                     *
02198       **********************************************************************************************/
02199     for (int i=0; i<in_c_p_d.size(); i++) {
02200       in_c_p_d[i] = Teuchos::ScalarTraits<double>::random();
02201     }
02202     for (int ic=0; ic < c; ic++) {
02203       for (int ip=0; ip < p; ip++) {
02204         for (int i=0; i<d1*d1; i++) {
02205           (&data_c_p_d_d(ic, ip, 0, 0))[i] = Teuchos::ScalarTraits<double>::random();
02206         }
02207         RealSpaceTools<double>::inverse(&datainv_c_p_d_d(ic, ip, 0, 0), &data_c_p_d_d(ic, ip, 0, 0), d1);
02208         RealSpaceTools<double>::transpose(&datainvtrn_c_p_d_d(ic, ip, 0, 0), &datainv_c_p_d_d(ic, ip, 0, 0), d1);
02209       }
02210     }
02211     for (int ic=0; ic < c; ic++) {
02212       for (int ip=0; ip < 1; ip++) {
02213         for (int i=0; i<d1*d1; i++) {
02214           (&data_c_1_d_d(ic, 0, 0, 0))[i] = Teuchos::ScalarTraits<double>::random();
02215         }
02216         RealSpaceTools<double>::inverse(&datainv_c_1_d_d(ic, 0, 0, 0), &data_c_1_d_d(ic, 0, 0, 0), d1);
02217         RealSpaceTools<double>::transpose(&datainvtrn_c_1_d_d(ic, 0, 0, 0), &datainv_c_1_d_d(ic, 0, 0, 0), d1);
02218       }
02219     }
02220     // Tensor values vary by point:
02221     art::matvecProductDataData<double>(out_c_p_d, data_c_p_d_d, in_c_p_d);
02222     art::matvecProductDataData<double>(outi_c_p_d, datainvtrn_c_p_d_d, out_c_p_d, 't');
02223     rst::subtract(&outi_c_p_d[0], &in_c_p_d[0], outi_c_p_d.size());
02224     if (rst::vectorNorm(&outi_c_p_d[0], outi_c_p_d.size(), NORM_ONE) > zero) {
02225       *outStream << "\n\nINCORRECT matvecProductDataData (11): check matrix inverse transpose property\n\n";
02226       errorFlag = -1000;
02227     }
02228     // Tensor values do not vary by point
02229     art::matvecProductDataData<double>(out_c_p_d, data_c_1_d_d, in_c_p_d);
02230     art::matvecProductDataData<double>(outi_c_p_d, datainvtrn_c_1_d_d, out_c_p_d, 't');
02231     rst::subtract(&outi_c_p_d[0], &in_c_p_d[0], outi_c_p_d.size());
02232     if (rst::vectorNorm(&outi_c_p_d[0], outi_c_p_d.size(), NORM_ONE) > zero) {
02233       *outStream << "\n\nINCORRECT matvecProductDataData (12): check matrix inverse transpose property\n\n";
02234       errorFlag = -1000;
02235     }
02236   }// test 6.c scope
02237   
02238   
02239   *outStream \
02240     << "\n"
02241     << "===============================================================================\n"\
02242     << "| TEST 6.d: matvecProductDataData random tests: branch inputDataRight(P,D)    |\n"\
02243     << "===============================================================================\n";
02244   /*
02245    * Test derived from Test 5.d
02246    */
02247   {// test 6.d scope
02248     int c=5, p=9, d1=3;
02249     double zero = INTREPID_TOL*10000.0;
02250     
02251     FieldContainer<double> in_p_d(p, d1);
02252     FieldContainer<double> in_c_p_d(c, p, d1);
02253     FieldContainer<double> data_c_p(c, p);
02254     FieldContainer<double> datainv_c_p(c, p);
02255     FieldContainer<double> data_c_1(c, 1);
02256     FieldContainer<double> datainv_c_1(c, 1);
02257     FieldContainer<double> data_c_p_d(c, p, d1);
02258     FieldContainer<double> datainv_c_p_d(c, p, d1);
02259     FieldContainer<double> data_c_1_d(c, 1, d1);
02260     FieldContainer<double> datainv_c_1_d(c, 1, d1);
02261     FieldContainer<double> data_c_p_d_d(c, p, d1, d1);
02262     FieldContainer<double> datainv_c_p_d_d(c, p, d1, d1);
02263     FieldContainer<double> datainvtrn_c_p_d_d(c, p, d1, d1);
02264     FieldContainer<double> data_c_1_d_d(c, 1, d1, d1);
02265     FieldContainer<double> datainv_c_1_d_d(c, 1, d1, d1);
02266     FieldContainer<double> datainvtrn_c_1_d_d(c, 1, d1, d1);
02267     FieldContainer<double> data_c_p_one(c, p);
02268     FieldContainer<double> data_c_1_one(c, 1);
02269     FieldContainer<double> out_c_p_d(c, p, d1);
02270     FieldContainer<double> outi_c_p_d(c, p, d1);
02271     /***********************************************************************************************
02272       *                          Constant diagonal tensor: inputData(C,P)                          *
02273       **********************************************************************************************/
02274     for (int i=0; i<in_p_d.size(); i++) {
02275       in_p_d[i] = Teuchos::ScalarTraits<double>::random();
02276     }
02277     for (int i=0; i<data_c_p.size(); i++) {
02278       data_c_p[i] = Teuchos::ScalarTraits<double>::random();
02279       datainv_c_p[i] = 1.0 / data_c_p[i];
02280       data_c_p_one[i] = 1.0;
02281     }
02282     for (int i=0; i<data_c_1.size(); i++) {
02283       data_c_1[i] = Teuchos::ScalarTraits<double>::random();
02284       datainv_c_1[i] = 1.0 / data_c_1[i];
02285     }
02286     // Tensor values vary by point
02287     art::matvecProductDataData<double>(in_c_p_d, data_c_p_one, in_p_d);
02288     art::matvecProductDataData<double>(out_c_p_d, data_c_p, in_p_d);
02289     art::matvecProductDataData<double>(outi_c_p_d, datainv_c_p, out_c_p_d);
02290     rst::subtract(&outi_c_p_d[0], &in_c_p_d[0], outi_c_p_d.size());
02291     if (rst::vectorNorm(&outi_c_p_d[0], outi_c_p_d.size(), NORM_ONE) > zero) {
02292       *outStream << "\n\nINCORRECT matvecProductDataData (13): check scalar inverse property\n\n";
02293       errorFlag = -1000;
02294     }
02295     // Tensor values do not vary by point
02296     art::matvecProductDataData<double>(in_c_p_d, data_c_p_one, in_p_d);
02297     art::matvecProductDataData<double>(out_c_p_d, data_c_1, in_p_d);
02298     art::matvecProductDataData<double>(outi_c_p_d, datainv_c_1, out_c_p_d);
02299     rst::subtract(&outi_c_p_d[0], &in_c_p_d[0], outi_c_p_d.size());
02300     if (rst::vectorNorm(&outi_c_p_d[0], outi_c_p_d.size(), NORM_ONE) > zero) {
02301       *outStream << "\n\nINCORRECT matvecProductDataData (14): check scalar inverse property\n\n";
02302       errorFlag = -1000;
02303     }
02304     /***********************************************************************************************
02305       *                       Non-constant diagonal tensor: inputData(C,P,D)                       *
02306       **********************************************************************************************/
02307     for (int i=0; i<in_p_d.size(); i++) {
02308       in_p_d[i] = Teuchos::ScalarTraits<double>::random();
02309     }
02310     for (int i=0; i<data_c_p_d.size(); i++) {
02311       data_c_p_d[i] = Teuchos::ScalarTraits<double>::random();
02312       datainv_c_p_d[i] = 1.0 / data_c_p_d[i];
02313     }
02314     for (int i=0; i<data_c_1_d.size(); i++) {
02315       data_c_1_d[i] = Teuchos::ScalarTraits<double>::random();
02316       datainv_c_1_d[i] = 1.0 / data_c_1_d[i];
02317     }
02318     // Tensor values vary by point:
02319     art::matvecProductDataData<double>(in_c_p_d, data_c_p_one, in_p_d);
02320     art::matvecProductDataData<double>(out_c_p_d, data_c_p_d, in_p_d);
02321     art::matvecProductDataData<double>(outi_c_p_d, datainv_c_p_d, out_c_p_d);
02322     rst::subtract(&outi_c_p_d[0], &in_c_p_d[0], outi_c_p_d.size());
02323     if (rst::vectorNorm(&outi_c_p_d[0], outi_c_p_d.size(), NORM_ONE) > zero) {
02324       *outStream << "\n\nINCORRECT matvecProductDataData (15): check scalar inverse property\n\n";
02325       errorFlag = -1000;
02326     }
02327     // Tensor values do not vary by point:
02328     art::matvecProductDataData<double>(in_c_p_d, data_c_p_one, in_p_d);
02329     art::matvecProductDataData<double>(out_c_p_d, data_c_1_d, in_p_d);
02330     art::matvecProductDataData<double>(outi_c_p_d, datainv_c_1_d, out_c_p_d);
02331     rst::subtract(&outi_c_p_d[0], &in_c_p_d[0], outi_c_p_d.size());
02332     if (rst::vectorNorm(&outi_c_p_d[0], outi_c_p_d.size(), NORM_ONE) > zero) {
02333       *outStream << "\n\nINCORRECT matvecProductDataData (16): check scalar inverse property\n\n";
02334       errorFlag = -1000;
02335     }
02336     /***********************************************************************************************
02337       *                              Full tensor: inputData(C,P,D,D)                               *
02338       **********************************************************************************************/
02339     for (int i=0; i<in_p_d.size(); i++) {
02340       in_p_d[i] = Teuchos::ScalarTraits<double>::random();
02341     }
02342     for (int ic=0; ic < c; ic++) {
02343       for (int ip=0; ip < p; ip++) {
02344         for (int i=0; i<d1*d1; i++) {
02345           (&data_c_p_d_d(ic, ip, 0, 0))[i] = Teuchos::ScalarTraits<double>::random();
02346         }
02347         RealSpaceTools<double>::inverse(&datainv_c_p_d_d(ic, ip, 0, 0), &data_c_p_d_d(ic, ip, 0, 0), d1);
02348       }
02349     }
02350     for (int ic=0; ic < c; ic++) {
02351       for (int ip=0; ip < 1; ip++) {
02352         for (int i=0; i<d1*d1; i++) {
02353           (&data_c_1_d_d(ic, 0, 0, 0))[i] = Teuchos::ScalarTraits<double>::random();
02354         }
02355         RealSpaceTools<double>::inverse(&datainv_c_1_d_d(ic, 0, 0, 0), &data_c_1_d_d(ic, 0, 0, 0), d1);
02356       }
02357     }
02358     // Tensor values vary by point: test "N" and "T" (no-transpose/transpose) options
02359     art::matvecProductDataData<double>(in_c_p_d, data_c_p_one, in_p_d);
02360     art::matvecProductDataData<double>(out_c_p_d, data_c_p_d_d, in_p_d);
02361     art::matvecProductDataData<double>(outi_c_p_d, datainv_c_p_d_d, out_c_p_d);
02362     rst::subtract(&outi_c_p_d[0], &in_c_p_d[0], outi_c_p_d.size());
02363     if (rst::vectorNorm(&outi_c_p_d[0], outi_c_p_d.size(), NORM_ONE) > zero) {
02364       *outStream << "\n\nINCORRECT matvecProductDataData (17): check matrix inverse property\n\n";
02365       errorFlag = -1000;
02366     }
02367     art::matvecProductDataData<double>(in_c_p_d, data_c_p_one, in_p_d);
02368     art::matvecProductDataData<double>(out_c_p_d, data_c_p_d_d, in_p_d, 't');
02369     art::matvecProductDataData<double>(outi_c_p_d, datainv_c_p_d_d, out_c_p_d, 't');
02370     rst::subtract(&outi_c_p_d[0], &in_c_p_d[0], outi_c_p_d.size());
02371     if (rst::vectorNorm(&outi_c_p_d[0], outi_c_p_d.size(), NORM_ONE) > zero) {
02372       *outStream << "\n\nINCORRECT matvecProductDataData (18): check matrix inverse property, w/ double transpose\n\n";
02373       errorFlag = -1000;
02374     }
02375     // Tensor values do not vary by point: test "N" and "T" (no-transpose/transpose) options
02376     art::matvecProductDataData<double>(in_c_p_d, data_c_p_one, in_p_d);
02377     art::matvecProductDataData<double>(out_c_p_d, data_c_1_d_d, in_p_d);
02378     art::matvecProductDataData<double>(outi_c_p_d, datainv_c_1_d_d, out_c_p_d);
02379     rst::subtract(&outi_c_p_d[0], &in_c_p_d[0], outi_c_p_d.size());
02380     if (rst::vectorNorm(&outi_c_p_d[0], outi_c_p_d.size(), NORM_ONE) > zero) {
02381       *outStream << "\n\nINCORRECT matvecProductDataData (19): check matrix inverse property\n\n";
02382       errorFlag = -1000;
02383     }
02384     art::matvecProductDataData<double>(in_c_p_d, data_c_p_one, in_p_d);
02385     art::matvecProductDataData<double>(out_c_p_d, data_c_1_d_d, in_p_d, 't');
02386     art::matvecProductDataData<double>(outi_c_p_d, datainv_c_1_d_d, out_c_p_d, 't');
02387     rst::subtract(&outi_c_p_d[0], &in_c_p_d[0], outi_c_p_d.size());
02388     if (rst::vectorNorm(&outi_c_p_d[0], outi_c_p_d.size(), NORM_ONE) > zero) {
02389       *outStream << "\n\nINCORRECT matvecProductDataData (20): check matrix inverse property, w/ double transpose\n\n";
02390       errorFlag = -1000;
02391     }
02392     /***********************************************************************************************
02393       *             Full tensor: inputData(C,P,D,D) test inverse transpose                         *
02394       **********************************************************************************************/
02395     for (int i=0; i<in_p_d.size(); i++) {
02396       in_p_d[i] = Teuchos::ScalarTraits<double>::random();
02397     }
02398     for (int ic=0; ic < c; ic++) {
02399       for (int ip=0; ip < p; ip++) {
02400         for (int i=0; i<d1*d1; i++) {
02401           (&data_c_p_d_d(ic, ip, 0, 0))[i] = Teuchos::ScalarTraits<double>::random();
02402         }
02403         RealSpaceTools<double>::inverse(&datainv_c_p_d_d(ic, ip, 0, 0), &data_c_p_d_d(ic, ip, 0, 0), d1);
02404         RealSpaceTools<double>::transpose(&datainvtrn_c_p_d_d(ic, ip, 0, 0), &datainv_c_p_d_d(ic, ip, 0, 0), d1);
02405       }
02406     }
02407     for (int ic=0; ic < c; ic++) {
02408       for (int ip=0; ip < 1; ip++) {
02409         for (int i=0; i<d1*d1; i++) {
02410           (&data_c_1_d_d(ic, 0, 0, 0))[i] = Teuchos::ScalarTraits<double>::random();
02411         }
02412         RealSpaceTools<double>::inverse(&datainv_c_1_d_d(ic, 0, 0, 0), &data_c_1_d_d(ic, 0, 0, 0), d1);
02413         RealSpaceTools<double>::transpose(&datainvtrn_c_1_d_d(ic, 0, 0, 0), &datainv_c_1_d_d(ic, 0, 0, 0), d1);
02414       }
02415     }
02416     // Tensor values vary by point:
02417     art::matvecProductDataData<double>(in_c_p_d, data_c_p_one, in_p_d);
02418     art::matvecProductDataData<double>(out_c_p_d, data_c_p_d_d, in_p_d);
02419     art::matvecProductDataData<double>(outi_c_p_d, datainvtrn_c_p_d_d, out_c_p_d, 't');
02420     rst::subtract(&outi_c_p_d[0], &in_c_p_d[0], outi_c_p_d.size());
02421     if (rst::vectorNorm(&outi_c_p_d[0], outi_c_p_d.size(), NORM_ONE) > zero) {
02422       *outStream << "\n\nINCORRECT matvecProductDataData (21): check matrix inverse transpose property\n\n";
02423       errorFlag = -1000;
02424     }
02425     // Tensor values do not vary by point:
02426     art::matvecProductDataData<double>(in_c_p_d, data_c_p_one, in_p_d);
02427     art::matvecProductDataData<double>(out_c_p_d, data_c_1_d_d, in_p_d);
02428     art::matvecProductDataData<double>(outi_c_p_d, datainvtrn_c_1_d_d, out_c_p_d, 't');
02429     rst::subtract(&outi_c_p_d[0], &in_c_p_d[0], outi_c_p_d.size());
02430     if (rst::vectorNorm(&outi_c_p_d[0], outi_c_p_d.size(), NORM_ONE) > zero) {
02431       *outStream << "\n\nINCORRECT matvecProductDataData (22): check matrix inverse transpose property\n\n";
02432       errorFlag = -1000;
02433     }    
02434   }// test 6.d scope
02435   
02436 
02437   
02438   *outStream \
02439     << "\n"
02440     << "===============================================================================\n"\
02441     << "| TEST 7.a: matmatProductDataField random tests: branch inputFields(C,F,P,D,D)|\n"\
02442     << "===============================================================================\n";
02443   {// Test 7.a scope
02444     int c=5, p=9, f=7, d1=3;
02445     double zero = INTREPID_TOL*10000.0;
02446     
02447     FieldContainer<double> in_c_f_p_d_d(c, f, p, d1, d1);
02448     FieldContainer<double> data_c_p(c, p);
02449     FieldContainer<double> datainv_c_p(c, p);
02450     FieldContainer<double> data_c_1(c, 1);
02451     FieldContainer<double> datainv_c_1(c, 1);
02452     FieldContainer<double> data_c_p_d(c, p, d1);
02453     FieldContainer<double> datainv_c_p_d(c, p, d1);
02454     FieldContainer<double> data_c_1_d(c, 1, d1);
02455     FieldContainer<double> datainv_c_1_d(c, 1, d1);
02456     FieldContainer<double> data_c_p_d_d(c, p, d1, d1);
02457     FieldContainer<double> datainv_c_p_d_d(c, p, d1, d1);
02458     FieldContainer<double> datainvtrn_c_p_d_d(c, p, d1, d1);
02459     FieldContainer<double> data_c_1_d_d(c, 1, d1, d1);
02460     FieldContainer<double> datainv_c_1_d_d(c, 1, d1, d1);
02461     FieldContainer<double> datainvtrn_c_1_d_d(c, 1, d1, d1);
02462     FieldContainer<double> out_c_f_p_d_d(c, f, p, d1, d1);
02463     FieldContainer<double> outi_c_f_p_d_d(c, f, p, d1, d1);
02464     /***********************************************************************************************
02465      *                          Constant diagonal tensor: inputData(C,P)                          *
02466      **********************************************************************************************/
02467     for (int i=0; i<in_c_f_p_d_d.size(); i++) {
02468       in_c_f_p_d_d[i] = Teuchos::ScalarTraits<double>::random();
02469     }
02470     for (int i=0; i<data_c_p.size(); i++) {
02471       data_c_p[i] = Teuchos::ScalarTraits<double>::random();
02472       datainv_c_p[i] = 1.0 / data_c_p[i];
02473     }
02474     for (int i=0; i<data_c_1.size(); i++) {
02475       data_c_1[i] = Teuchos::ScalarTraits<double>::random();
02476       datainv_c_1[i] = 1.0 / data_c_1[i];
02477     }
02478     // Tensor values vary by point:
02479     art::matmatProductDataField<double>(out_c_f_p_d_d, data_c_p, in_c_f_p_d_d);
02480     art::matmatProductDataField<double>(out_c_f_p_d_d, datainv_c_p, out_c_f_p_d_d);
02481     rst::subtract(&out_c_f_p_d_d[0], &in_c_f_p_d_d[0], out_c_f_p_d_d.size());
02482     if (rst::vectorNorm(&out_c_f_p_d_d[0], out_c_f_p_d_d.size(), NORM_ONE) > zero) {
02483       *outStream << "\n\nINCORRECT matmatProductDataField (1): check scalar inverse property\n\n";
02484       errorFlag = -1000;
02485     }
02486     // Tensor values do not vary by point:
02487     art::matmatProductDataField<double>(out_c_f_p_d_d, data_c_1, in_c_f_p_d_d);
02488     art::matmatProductDataField<double>(out_c_f_p_d_d, datainv_c_1, out_c_f_p_d_d);
02489     rst::subtract(&out_c_f_p_d_d[0], &in_c_f_p_d_d[0], out_c_f_p_d_d.size());
02490     if (rst::vectorNorm(&out_c_f_p_d_d[0], out_c_f_p_d_d.size(), NORM_ONE) > zero) {
02491       *outStream << "\n\nINCORRECT matmatProductDataField (2): check scalar inverse property\n\n";
02492       errorFlag = -1000;
02493     }
02494     /***********************************************************************************************
02495      *                     Non-onstant diagonal tensor: inputData(C,P,D)                          *
02496      **********************************************************************************************/
02497     for (int i=0; i<in_c_f_p_d_d.size(); i++) {
02498       in_c_f_p_d_d[i] = Teuchos::ScalarTraits<double>::random();
02499     }
02500     for (int i=0; i<data_c_p_d.size(); i++) {
02501       data_c_p_d[i] = Teuchos::ScalarTraits<double>::random();
02502       datainv_c_p_d[i] = 1.0 / data_c_p_d[i];
02503     }
02504     for (int i=0; i<data_c_1_d.size(); i++) {
02505       data_c_1_d[i] = Teuchos::ScalarTraits<double>::random();
02506       datainv_c_1_d[i] = 1.0 / data_c_1_d[i];
02507     }
02508     // Tensor values vary by point:
02509     art::matmatProductDataField<double>(out_c_f_p_d_d, data_c_p_d, in_c_f_p_d_d);
02510     art::matmatProductDataField<double>(out_c_f_p_d_d, datainv_c_p_d, out_c_f_p_d_d);
02511     rst::subtract(&out_c_f_p_d_d[0], &in_c_f_p_d_d[0], out_c_f_p_d_d.size());
02512     if (rst::vectorNorm(&out_c_f_p_d_d[0], out_c_f_p_d_d.size(), NORM_ONE) > zero) {
02513       *outStream << "\n\nINCORRECT matmatProductDataField (3): check scalar inverse property\n\n";
02514       errorFlag = -1000;
02515     }
02516     // Tensor values do not vary by point
02517     art::matmatProductDataField<double>(out_c_f_p_d_d, data_c_1_d, in_c_f_p_d_d);
02518     art::matmatProductDataField<double>(out_c_f_p_d_d, datainv_c_1_d, out_c_f_p_d_d);
02519     rst::subtract(&out_c_f_p_d_d[0], &in_c_f_p_d_d[0], out_c_f_p_d_d.size());
02520     if (rst::vectorNorm(&out_c_f_p_d_d[0], out_c_f_p_d_d.size(), NORM_ONE) > zero) {
02521       *outStream << "\n\nINCORRECT matmatProductDataField (4): check scalar inverse property\n\n";
02522       errorFlag = -1000;
02523     }
02524     /***********************************************************************************************
02525      *                               Full tensor: inputData(C,P,D,D)                              *
02526      **********************************************************************************************/
02527     for (int i=0; i<in_c_f_p_d_d.size(); i++) {
02528       in_c_f_p_d_d[i] = Teuchos::ScalarTraits<double>::random();
02529     }
02530     for (int ic=0; ic < c; ic++) {
02531       for (int ip=0; ip < p; ip++) {
02532         for (int i=0; i<d1*d1; i++) {
02533           (&data_c_p_d_d(ic, ip, 0, 0))[i] = Teuchos::ScalarTraits<double>::random();
02534         }
02535         RealSpaceTools<double>::inverse(&datainv_c_p_d_d(ic, ip, 0, 0), &data_c_p_d_d(ic, ip, 0, 0), d1);
02536       }
02537     }
02538     for (int ic=0; ic < c; ic++) {
02539       for (int ip=0; ip < 1; ip++) {
02540         for (int i=0; i<d1*d1; i++) {
02541           (&data_c_1_d_d(ic, 0, 0, 0))[i] = Teuchos::ScalarTraits<double>::random();
02542         }
02543         RealSpaceTools<double>::inverse(&datainv_c_1_d_d(ic, 0, 0, 0), &data_c_1_d_d(ic, 0, 0, 0), d1);
02544       }
02545     }
02546     // Tensor values vary by point: test "N" and "T" options (no transpose/transpose)
02547     art::matmatProductDataField<double>(out_c_f_p_d_d, data_c_p_d_d, in_c_f_p_d_d);
02548     art::matmatProductDataField<double>(outi_c_f_p_d_d, datainv_c_p_d_d, out_c_f_p_d_d);
02549     rst::subtract(&outi_c_f_p_d_d[0], &in_c_f_p_d_d[0], outi_c_f_p_d_d.size());
02550     if (rst::vectorNorm(&outi_c_f_p_d_d[0], outi_c_f_p_d_d.size(), NORM_ONE) > zero) {
02551       *outStream << "\n\nINCORRECT matmatProductDataField (5): check matrix inverse property\n\n";
02552       errorFlag = -1000;
02553     }
02554     art::matmatProductDataField<double>(out_c_f_p_d_d, data_c_p_d_d, in_c_f_p_d_d, 't');
02555     art::matmatProductDataField<double>(outi_c_f_p_d_d, datainv_c_p_d_d, out_c_f_p_d_d, 't');
02556     rst::subtract(&outi_c_f_p_d_d[0], &in_c_f_p_d_d[0], outi_c_f_p_d_d.size());
02557     if (rst::vectorNorm(&outi_c_f_p_d_d[0], outi_c_f_p_d_d.size(), NORM_ONE) > zero) {
02558       *outStream << "\n\nINCORRECT matmatProductDataField (6): check matrix inverse property, w/ double transpose\n\n";
02559       errorFlag = -1000;
02560     }
02561     // Tensor values do not vary by point: test "N" and "T" options (no transpose/transpose)
02562     art::matmatProductDataField<double>(out_c_f_p_d_d, data_c_1_d_d, in_c_f_p_d_d);
02563     art::matmatProductDataField<double>(outi_c_f_p_d_d, datainv_c_1_d_d, out_c_f_p_d_d);
02564     rst::subtract(&outi_c_f_p_d_d[0], &in_c_f_p_d_d[0], outi_c_f_p_d_d.size());
02565     if (rst::vectorNorm(&outi_c_f_p_d_d[0], outi_c_f_p_d_d.size(), NORM_ONE) > zero) {
02566       *outStream << "\n\nINCORRECT matmatProductDataField (7): check matrix inverse property\n\n";
02567       errorFlag = -1000;
02568     }
02569     art::matmatProductDataField<double>(out_c_f_p_d_d, data_c_1_d_d, in_c_f_p_d_d,'t');
02570     art::matmatProductDataField<double>(outi_c_f_p_d_d, datainv_c_1_d_d, out_c_f_p_d_d, 't');
02571     rst::subtract(&outi_c_f_p_d_d[0], &in_c_f_p_d_d[0], outi_c_f_p_d_d.size());
02572     if (rst::vectorNorm(&outi_c_f_p_d_d[0], outi_c_f_p_d_d.size(), NORM_ONE) > zero) {
02573       *outStream << "\n\nINCORRECT matmatProductDataField (8): check matrix inverse property, w/ double transpose\n\n";
02574       errorFlag = -1000;
02575     }
02576     /***********************************************************************************************
02577      *             Full tensor: inputData(C,P,D,D) test inverse transpose                         *
02578      **********************************************************************************************/
02579     for (int i=0; i<in_c_f_p_d_d.size(); i++) {
02580       in_c_f_p_d_d[i] = Teuchos::ScalarTraits<double>::random();
02581     }
02582     for (int ic=0; ic < c; ic++) {
02583       for (int ip=0; ip < p; ip++) {
02584         for (int i=0; i<d1*d1; i++) {
02585           (&data_c_p_d_d(ic, ip, 0, 0))[i] = Teuchos::ScalarTraits<double>::random();
02586         }
02587         RealSpaceTools<double>::inverse(&datainv_c_p_d_d(ic, ip, 0, 0), &data_c_p_d_d(ic, ip, 0, 0), d1);
02588         RealSpaceTools<double>::transpose(&datainvtrn_c_p_d_d(ic, ip, 0, 0), &datainv_c_p_d_d(ic, ip, 0, 0), d1);
02589       }
02590     }
02591     for (int ic=0; ic < c; ic++) {
02592       for (int ip=0; ip < 1; ip++) {
02593         for (int i=0; i<d1*d1; i++) {
02594           (&data_c_1_d_d(ic, 0, 0, 0))[i] = Teuchos::ScalarTraits<double>::random();
02595         }
02596         RealSpaceTools<double>::inverse(&datainv_c_1_d_d(ic, 0, 0, 0), &data_c_1_d_d(ic, 0, 0, 0), d1);
02597         RealSpaceTools<double>::transpose(&datainvtrn_c_1_d_d(ic, 0, 0, 0), &datainv_c_1_d_d(ic, 0, 0, 0), d1);
02598       }
02599     }
02600     // Tensor values vary by point:
02601     art::matmatProductDataField<double>(out_c_f_p_d_d, data_c_p_d_d, in_c_f_p_d_d);
02602     art::matmatProductDataField<double>(outi_c_f_p_d_d, datainvtrn_c_p_d_d, out_c_f_p_d_d, 't');
02603     rst::subtract(&outi_c_f_p_d_d[0], &in_c_f_p_d_d[0], outi_c_f_p_d_d.size());
02604     if (rst::vectorNorm(&outi_c_f_p_d_d[0], outi_c_f_p_d_d.size(), NORM_ONE) > zero) {
02605       *outStream << "\n\nINCORRECT matmatProductDataField (9): check matrix inverse transpose property\n\n";
02606       errorFlag = -1000;
02607     }
02608     // Tensor values do not vary by point
02609     art::matmatProductDataField<double>(out_c_f_p_d_d, data_c_1_d_d, in_c_f_p_d_d);
02610     art::matmatProductDataField<double>(outi_c_f_p_d_d, datainvtrn_c_1_d_d, out_c_f_p_d_d, 't');
02611     rst::subtract(&outi_c_f_p_d_d[0], &in_c_f_p_d_d[0], outi_c_f_p_d_d.size());
02612     if (rst::vectorNorm(&outi_c_f_p_d_d[0], outi_c_f_p_d_d.size(), NORM_ONE) > zero) {
02613       *outStream << "\n\nINCORRECT matmatProductDataField (10): check matrix inverse transpose property\n\n";
02614       errorFlag = -1000;
02615     }
02616   }// test 7.a scope
02617   
02618   
02619   
02620   *outStream \
02621     << "\n"
02622     << "===============================================================================\n"\
02623     << "| TEST 7.b: matmatProductDataField random tests: branch inputFields(F,P,D,D)  |\n"\
02624     << "===============================================================================\n";
02625   {// Test 7.b scope
02626     int c=5, p=9, f=7, d1=3;
02627     double zero = INTREPID_TOL*10000.0;
02628     
02629     FieldContainer<double> in_f_p_d_d(f, p, d1, d1);
02630     FieldContainer<double> in_c_f_p_d_d(c, f, p, d1, d1);
02631     FieldContainer<double> data_c_p(c, p);
02632     FieldContainer<double> datainv_c_p(c, p);
02633     FieldContainer<double> data_c_1(c, 1);
02634     FieldContainer<double> datainv_c_1(c, 1);
02635     FieldContainer<double> data_c_p_d(c, p, d1);
02636     FieldContainer<double> datainv_c_p_d(c, p, d1);
02637     FieldContainer<double> data_c_1_d(c, 1, d1);
02638     FieldContainer<double> datainv_c_1_d(c, 1, d1);
02639     FieldContainer<double> data_c_p_d_d(c, p, d1, d1);
02640     FieldContainer<double> datainv_c_p_d_d(c, p, d1, d1);
02641     FieldContainer<double> datainvtrn_c_p_d_d(c, p, d1, d1);
02642     FieldContainer<double> data_c_1_d_d(c, 1, d1, d1);
02643     FieldContainer<double> datainv_c_1_d_d(c, 1, d1, d1);
02644     FieldContainer<double> datainvtrn_c_1_d_d(c, 1, d1, d1);
02645     FieldContainer<double> data_c_p_one(c, p);
02646     FieldContainer<double> data_c_1_one(c, 1);
02647     FieldContainer<double> out_c_f_p_d_d(c, f, p, d1, d1);
02648     FieldContainer<double> outi_c_f_p_d_d(c, f, p, d1, d1);
02649     /***********************************************************************************************
02650       *                          Constant diagonal tensor: inputData(C,P)                          *
02651       **********************************************************************************************/
02652     for (int i=0; i<in_f_p_d_d.size(); i++) {
02653       in_f_p_d_d[i] = Teuchos::ScalarTraits<double>::random();
02654     }
02655     for (int i=0; i<data_c_p.size(); i++) {
02656       data_c_p[i] = Teuchos::ScalarTraits<double>::random();
02657       datainv_c_p[i] = 1.0 / data_c_p[i];
02658       data_c_p_one[i] = 1.0;
02659     }
02660     for (int i=0; i<data_c_1.size(); i++) {
02661       data_c_1[i] = Teuchos::ScalarTraits<double>::random();
02662       datainv_c_1[i] = 1.0 / data_c_1[i];
02663     }
02664     
02665     // Tensor values vary by point
02666     art::matmatProductDataField<double>(in_c_f_p_d_d, data_c_p_one, in_f_p_d_d);
02667     art::matmatProductDataField<double>(out_c_f_p_d_d, data_c_p, in_f_p_d_d);
02668     art::matmatProductDataField<double>(outi_c_f_p_d_d, datainv_c_p, out_c_f_p_d_d);
02669     rst::subtract(&outi_c_f_p_d_d[0], &in_c_f_p_d_d[0], outi_c_f_p_d_d.size());
02670     if (rst::vectorNorm(&outi_c_f_p_d_d[0], outi_c_f_p_d_d.size(), NORM_ONE) > zero) {
02671       *outStream << "\n\nINCORRECT matmatProductDataField (11): check scalar inverse property\n\n";
02672       errorFlag = -1000;
02673     }
02674     // Tensor values do not vary by point
02675     art::matmatProductDataField<double>(in_c_f_p_d_d, data_c_p_one, in_f_p_d_d);
02676     art::matmatProductDataField<double>(out_c_f_p_d_d, data_c_1, in_f_p_d_d);
02677     art::matmatProductDataField<double>(outi_c_f_p_d_d, datainv_c_1, out_c_f_p_d_d);
02678     rst::subtract(&outi_c_f_p_d_d[0], &in_c_f_p_d_d[0], outi_c_f_p_d_d.size());
02679     if (rst::vectorNorm(&outi_c_f_p_d_d[0], outi_c_f_p_d_d.size(), NORM_ONE) > zero) {
02680       *outStream << "\n\nINCORRECT matmatProductDataField (12): check scalar inverse property\n\n";
02681       errorFlag = -1000;
02682     }
02683     /***********************************************************************************************
02684       *                       Non-constant diagonal tensor: inputData(C,P,D)                       *
02685       **********************************************************************************************/
02686     for (int i=0; i<in_f_p_d_d.size(); i++) {
02687       in_f_p_d_d[i] = Teuchos::ScalarTraits<double>::random();
02688     }
02689     for (int i=0; i<data_c_p_d.size(); i++) {
02690       data_c_p_d[i] = Teuchos::ScalarTraits<double>::random();
02691       datainv_c_p_d[i] = 1.0 / data_c_p_d[i];
02692     }
02693     for (int i=0; i<data_c_1_d.size(); i++) {
02694       data_c_1_d[i] = Teuchos::ScalarTraits<double>::random();
02695       datainv_c_1_d[i] = 1.0 / data_c_1_d[i];
02696     }
02697     // Tensor values vary by point:
02698     art::matmatProductDataField<double>(in_c_f_p_d_d, data_c_p_one, in_f_p_d_d);
02699     art::matmatProductDataField<double>(out_c_f_p_d_d, data_c_p_d, in_f_p_d_d);
02700     art::matmatProductDataField<double>(outi_c_f_p_d_d, datainv_c_p_d, out_c_f_p_d_d);
02701     rst::subtract(&outi_c_f_p_d_d[0], &in_c_f_p_d_d[0], outi_c_f_p_d_d.size());
02702     if (rst::vectorNorm(&outi_c_f_p_d_d[0], outi_c_f_p_d_d.size(), NORM_ONE) > zero) {
02703       *outStream << "\n\nINCORRECT matmatProductDataField (13): check scalar inverse property\n\n";
02704       errorFlag = -1000;
02705     }
02706     // Tensor values do not vary by point:
02707     art::matmatProductDataField<double>(in_c_f_p_d_d, data_c_p_one, in_f_p_d_d);
02708     art::matmatProductDataField<double>(out_c_f_p_d_d, data_c_1_d, in_f_p_d_d);
02709     art::matmatProductDataField<double>(outi_c_f_p_d_d, datainv_c_1_d, out_c_f_p_d_d);
02710     rst::subtract(&outi_c_f_p_d_d[0], &in_c_f_p_d_d[0], outi_c_f_p_d_d.size());
02711     if (rst::vectorNorm(&outi_c_f_p_d_d[0], outi_c_f_p_d_d.size(), NORM_ONE) > zero) {
02712       *outStream << "\n\nINCORRECT matmatProductDataField (14): check scalar inverse property\n\n";
02713       errorFlag = -1000;
02714     }
02715     /***********************************************************************************************
02716       *                              Full tensor: inputData(C,P,D,D)                               *
02717       **********************************************************************************************/
02718     for (int i=0; i<in_f_p_d_d.size(); i++) {
02719       in_f_p_d_d[i] = Teuchos::ScalarTraits<double>::random();
02720     }
02721     for (int ic=0; ic < c; ic++) {
02722       for (int ip=0; ip < p; ip++) {
02723         for (int i=0; i<d1*d1; i++) {
02724           (&data_c_p_d_d(ic, ip, 0, 0))[i] = Teuchos::ScalarTraits<double>::random();
02725         }
02726         RealSpaceTools<double>::inverse(&datainv_c_p_d_d(ic, ip, 0, 0), &data_c_p_d_d(ic, ip, 0, 0), d1);
02727       }
02728     }
02729     for (int ic=0; ic < c; ic++) {
02730       for (int ip=0; ip < 1; ip++) {
02731         for (int i=0; i<d1*d1; i++) {
02732           (&data_c_1_d_d(ic, 0, 0, 0))[i] = Teuchos::ScalarTraits<double>::random();
02733         }
02734         RealSpaceTools<double>::inverse(&datainv_c_1_d_d(ic, 0, 0, 0), &data_c_1_d_d(ic, 0, 0, 0), d1);
02735       }
02736     }
02737     // Tensor values vary by point: test "N" and "T" (no-transpose/transpose) options
02738     art::matmatProductDataField<double>(in_c_f_p_d_d, data_c_p_one, in_f_p_d_d);
02739     art::matmatProductDataField<double>(out_c_f_p_d_d, data_c_p_d_d, in_f_p_d_d);
02740     art::matmatProductDataField<double>(outi_c_f_p_d_d, datainv_c_p_d_d, out_c_f_p_d_d);
02741     rst::subtract(&outi_c_f_p_d_d[0], &in_c_f_p_d_d[0], outi_c_f_p_d_d.size());
02742     if (rst::vectorNorm(&outi_c_f_p_d_d[0], outi_c_f_p_d_d.size(), NORM_ONE) > zero) {
02743       *outStream << "\n\nINCORRECT matmatProductDataField (15): check matrix inverse property\n\n";
02744       errorFlag = -1000;
02745     }
02746     art::matmatProductDataField<double>(in_c_f_p_d_d, data_c_p_one, in_f_p_d_d);
02747     art::matmatProductDataField<double>(out_c_f_p_d_d, data_c_p_d_d, in_f_p_d_d, 't');
02748     art::matmatProductDataField<double>(outi_c_f_p_d_d, datainv_c_p_d_d, out_c_f_p_d_d, 't');
02749     rst::subtract(&outi_c_f_p_d_d[0], &in_c_f_p_d_d[0], outi_c_f_p_d_d.size());
02750     if (rst::vectorNorm(&outi_c_f_p_d_d[0], outi_c_f_p_d_d.size(), NORM_ONE) > zero) {
02751       *outStream << "\n\nINCORRECT matmatProductDataField (16): check matrix inverse property, w/ double transpose\n\n";
02752       errorFlag = -1000;
02753     }
02754     // Tensor values do not vary by point: test "N" and "T" (no-transpose/transpose) options
02755     art::matmatProductDataField<double>(in_c_f_p_d_d, data_c_p_one, in_f_p_d_d);
02756     art::matmatProductDataField<double>(out_c_f_p_d_d, data_c_1_d_d, in_f_p_d_d);
02757     art::matmatProductDataField<double>(outi_c_f_p_d_d, datainv_c_1_d_d, out_c_f_p_d_d);
02758     rst::subtract(&outi_c_f_p_d_d[0], &in_c_f_p_d_d[0], outi_c_f_p_d_d.size());
02759     if (rst::vectorNorm(&outi_c_f_p_d_d[0], outi_c_f_p_d_d.size(), NORM_ONE) > zero) {
02760       *outStream << "\n\nINCORRECT matmatProductDataField (17): check matrix inverse property\n\n";
02761       errorFlag = -1000;
02762     }
02763     art::matmatProductDataField<double>(in_c_f_p_d_d, data_c_p_one, in_f_p_d_d);
02764     art::matmatProductDataField<double>(out_c_f_p_d_d, data_c_1_d_d, in_f_p_d_d, 't');
02765     art::matmatProductDataField<double>(outi_c_f_p_d_d, datainv_c_1_d_d, out_c_f_p_d_d, 't');
02766     rst::subtract(&outi_c_f_p_d_d[0], &in_c_f_p_d_d[0], outi_c_f_p_d_d.size());
02767     if (rst::vectorNorm(&outi_c_f_p_d_d[0], outi_c_f_p_d_d.size(), NORM_ONE) > zero) {
02768       *outStream << "\n\nINCORRECT matmatProductDataField (18): check matrix inverse property, w/ double transpose\n\n";
02769       errorFlag = -1000;
02770     }
02771     /***********************************************************************************************
02772       *             Full tensor: inputData(C,P,D,D) test inverse transpose                         *
02773       **********************************************************************************************/
02774     for (int i=0; i<in_f_p_d_d.size(); i++) {
02775       in_f_p_d_d[i] = Teuchos::ScalarTraits<double>::random();
02776     }
02777     for (int ic=0; ic < c; ic++) {
02778       for (int ip=0; ip < p; ip++) {
02779         for (int i=0; i<d1*d1; i++) {
02780           (&data_c_p_d_d(ic, ip, 0, 0))[i] = Teuchos::ScalarTraits<double>::random();
02781         }
02782         RealSpaceTools<double>::inverse(&datainv_c_p_d_d(ic, ip, 0, 0), &data_c_p_d_d(ic, ip, 0, 0), d1);
02783         RealSpaceTools<double>::transpose(&datainvtrn_c_p_d_d(ic, ip, 0, 0), &datainv_c_p_d_d(ic, ip, 0, 0), d1);
02784       }
02785     }
02786     for (int ic=0; ic < c; ic++) {
02787       for (int ip=0; ip < 1; ip++) {
02788         for (int i=0; i<d1*d1; i++) {
02789           (&data_c_1_d_d(ic, 0, 0, 0))[i] = Teuchos::ScalarTraits<double>::random();
02790         }
02791         RealSpaceTools<double>::inverse(&datainv_c_1_d_d(ic, 0, 0, 0), &data_c_1_d_d(ic, 0, 0, 0), d1);
02792         RealSpaceTools<double>::transpose(&datainvtrn_c_1_d_d(ic, 0, 0, 0), &datainv_c_1_d_d(ic, 0, 0, 0), d1);
02793       }
02794     }
02795     // Tensor values vary by point:
02796     art::matmatProductDataField<double>(in_c_f_p_d_d, data_c_p_one, in_f_p_d_d);
02797     art::matmatProductDataField<double>(out_c_f_p_d_d, data_c_p_d_d, in_f_p_d_d);
02798     art::matmatProductDataField<double>(outi_c_f_p_d_d, datainvtrn_c_p_d_d, out_c_f_p_d_d, 't');
02799     rst::subtract(&outi_c_f_p_d_d[0], &in_c_f_p_d_d[0], outi_c_f_p_d_d.size());
02800     if (rst::vectorNorm(&outi_c_f_p_d_d[0], outi_c_f_p_d_d.size(), NORM_ONE) > zero) {
02801       *outStream << "\n\nINCORRECT matmatProductDataField (19): check matrix inverse transpose property\n\n";
02802       errorFlag = -1000;
02803     }
02804     // Tensor values do not vary by point:
02805     art::matmatProductDataField<double>(in_c_f_p_d_d, data_c_p_one, in_f_p_d_d);
02806     art::matmatProductDataField<double>(out_c_f_p_d_d, data_c_1_d_d, in_f_p_d_d);
02807     art::matmatProductDataField<double>(outi_c_f_p_d_d, datainvtrn_c_1_d_d, out_c_f_p_d_d, 't');
02808     rst::subtract(&outi_c_f_p_d_d[0], &in_c_f_p_d_d[0], outi_c_f_p_d_d.size());
02809     if (rst::vectorNorm(&outi_c_f_p_d_d[0], outi_c_f_p_d_d.size(), NORM_ONE) > zero) {
02810       *outStream << "\n\nINCORRECT matmatProductDataField (20): check matrix inverse transpose property\n\n";
02811       errorFlag = -1000;
02812     }    
02813   }// test 7.b scope
02814   
02815   
02816   
02817   *outStream \
02818     << "\n"
02819     << "===============================================================================\n"\
02820     << "| TEST 8.a: matmatProductDataData random tests: branch inputDataRight(C,P,D,D)|\n"\
02821     << "===============================================================================\n";
02822   /*
02823    * Test derived from Test 7.a
02824    */
02825   {// test 8.a scope
02826     int c=5, p=9, d1=3;
02827     double zero = INTREPID_TOL*10000.0;
02828     
02829     FieldContainer<double> in_c_p_d_d(c, p, d1, d1);
02830     FieldContainer<double> out_c_p_d_d(c, p, d1, d1);
02831     FieldContainer<double> outi_c_p_d_d(c, p, d1, d1);
02832     
02833     FieldContainer<double> data_c_p(c, p);
02834     FieldContainer<double> datainv_c_p(c, p);
02835     FieldContainer<double> data_c_1(c, 1);
02836     FieldContainer<double> datainv_c_1(c, 1);
02837     FieldContainer<double> data_c_p_d(c, p, d1);
02838     FieldContainer<double> datainv_c_p_d(c, p, d1);
02839     FieldContainer<double> data_c_1_d(c, 1, d1);
02840     FieldContainer<double> datainv_c_1_d(c, 1, d1);
02841     FieldContainer<double> data_c_p_d_d(c, p, d1, d1);
02842     FieldContainer<double> datainv_c_p_d_d(c, p, d1, d1);
02843     FieldContainer<double> datainvtrn_c_p_d_d(c, p, d1, d1);
02844     FieldContainer<double> data_c_1_d_d(c, 1, d1, d1);
02845     FieldContainer<double> datainv_c_1_d_d(c, 1, d1, d1);
02846     FieldContainer<double> datainvtrn_c_1_d_d(c, 1, d1, d1);
02847     /***********************************************************************************************
02848       *                          Constant diagonal tensor: inputDataLeft(C,P)                      *
02849       **********************************************************************************************/
02850     for (int i=0; i<in_c_p_d_d.size(); i++) {
02851       in_c_p_d_d[i] = Teuchos::ScalarTraits<double>::random();
02852     }
02853     for (int i=0; i<data_c_p.size(); i++) {
02854       data_c_p[i] = Teuchos::ScalarTraits<double>::random();
02855       datainv_c_p[i] = 1.0 / data_c_p[i];
02856     }
02857     for (int i=0; i<data_c_1.size(); i++) {
02858       data_c_1[i] = Teuchos::ScalarTraits<double>::random();
02859       datainv_c_1[i] = 1.0 / data_c_1[i];
02860     }
02861     // Tensor values vary by point:
02862     art::matmatProductDataData<double>(out_c_p_d_d, data_c_p, in_c_p_d_d);
02863     art::matmatProductDataData<double>(out_c_p_d_d, datainv_c_p, out_c_p_d_d);
02864     rst::subtract(&out_c_p_d_d[0], &in_c_p_d_d[0], out_c_p_d_d.size());
02865     if (rst::vectorNorm(&out_c_p_d_d[0], out_c_p_d_d.size(), NORM_ONE) > zero) {
02866       *outStream << "\n\nINCORRECT matmatProductDataData (1): check scalar inverse property\n\n";
02867       errorFlag = -1000;
02868     }
02869     // Tensor values do not vary by point:
02870     art::matmatProductDataData<double>(out_c_p_d_d, data_c_1, in_c_p_d_d);
02871     art::matmatProductDataData<double>(out_c_p_d_d, datainv_c_1, out_c_p_d_d);
02872     rst::subtract(&out_c_p_d_d[0], &in_c_p_d_d[0], out_c_p_d_d.size());
02873     if (rst::vectorNorm(&out_c_p_d_d[0], out_c_p_d_d.size(), NORM_ONE) > zero) {
02874       *outStream << "\n\nINCORRECT matmatProductDataData (2): check scalar inverse property\n\n";
02875       errorFlag = -1000;
02876     }
02877     /***********************************************************************************************
02878       *                     Non-onstant diagonal tensor: inputDataLeft(C,P,D)                      *
02879       **********************************************************************************************/
02880     for (int i=0; i<in_c_p_d_d.size(); i++) {
02881       in_c_p_d_d[i] = Teuchos::ScalarTraits<double>::random();
02882     }
02883     for (int i=0; i<data_c_p_d.size(); i++) {
02884       data_c_p_d[i] = Teuchos::ScalarTraits<double>::random();
02885       datainv_c_p_d[i] = 1.0 / data_c_p_d[i];
02886     }
02887     for (int i=0; i<data_c_1_d.size(); i++) {
02888       data_c_1_d[i] = Teuchos::ScalarTraits<double>::random();
02889       datainv_c_1_d[i] = 1.0 / data_c_1_d[i];
02890     }
02891     // Tensor values vary by point:
02892     art::matmatProductDataData<double>(out_c_p_d_d, data_c_p_d, in_c_p_d_d);
02893     art::matmatProductDataData<double>(out_c_p_d_d, datainv_c_p_d, out_c_p_d_d);
02894     rst::subtract(&out_c_p_d_d[0], &in_c_p_d_d[0], out_c_p_d_d.size());
02895     if (rst::vectorNorm(&out_c_p_d_d[0], out_c_p_d_d.size(), NORM_ONE) > zero) {
02896       *outStream << "\n\nINCORRECT matmatProductDataData (3): check scalar inverse property\n\n";
02897       errorFlag = -1000;
02898     }
02899     // Tensor values do not vary by point
02900     art::matmatProductDataData<double>(out_c_p_d_d, data_c_1_d, in_c_p_d_d);
02901     art::matmatProductDataData<double>(out_c_p_d_d, datainv_c_1_d, out_c_p_d_d);
02902     rst::subtract(&out_c_p_d_d[0], &in_c_p_d_d[0], out_c_p_d_d.size());
02903     if (rst::vectorNorm(&out_c_p_d_d[0], out_c_p_d_d.size(), NORM_ONE) > zero) {
02904       *outStream << "\n\nINCORRECT matmatProductDataData (4): check scalar inverse property\n\n";
02905       errorFlag = -1000;
02906     }
02907     /***********************************************************************************************
02908       *                            Full tensor: inputDataLeft(C,P,D,D)                             *
02909       **********************************************************************************************/
02910     for (int i=0; i<in_c_p_d_d.size(); i++) {
02911       in_c_p_d_d[i] = Teuchos::ScalarTraits<double>::random();
02912     }
02913     for (int ic=0; ic < c; ic++) {
02914       for (int ip=0; ip < p; ip++) {
02915         for (int i=0; i<d1*d1; i++) {
02916           (&data_c_p_d_d(ic, ip, 0, 0))[i] = Teuchos::ScalarTraits<double>::random();
02917         }
02918         RealSpaceTools<double>::inverse(&datainv_c_p_d_d(ic, ip, 0, 0), &data_c_p_d_d(ic, ip, 0, 0), d1);
02919       }
02920     }
02921     for (int ic=0; ic < c; ic++) {
02922       for (int ip=0; ip < 1; ip++) {
02923         for (int i=0; i<d1*d1; i++) {
02924           (&data_c_1_d_d(ic, 0, 0, 0))[i] = Teuchos::ScalarTraits<double>::random();
02925         }
02926         RealSpaceTools<double>::inverse(&datainv_c_1_d_d(ic, 0, 0, 0), &data_c_1_d_d(ic, 0, 0, 0), d1);
02927       }
02928     }
02929     // Tensor values vary by point: test "N" and "T" options (no transpose/transpose)
02930     art::matmatProductDataData<double>(out_c_p_d_d, data_c_p_d_d, in_c_p_d_d);
02931     art::matmatProductDataData<double>(outi_c_p_d_d, datainv_c_p_d_d, out_c_p_d_d);
02932     rst::subtract(&outi_c_p_d_d[0], &in_c_p_d_d[0], outi_c_p_d_d.size());
02933     if (rst::vectorNorm(&outi_c_p_d_d[0], outi_c_p_d_d.size(), NORM_ONE) > zero) {
02934       *outStream << "\n\nINCORRECT matmatProductDataData (5): check matrix inverse property\n\n";
02935       errorFlag = -1000;
02936     }
02937     art::matmatProductDataData<double>(out_c_p_d_d, data_c_p_d_d, in_c_p_d_d, 't');
02938     art::matmatProductDataData<double>(outi_c_p_d_d, datainv_c_p_d_d, out_c_p_d_d, 't');
02939     rst::subtract(&outi_c_p_d_d[0], &in_c_p_d_d[0], outi_c_p_d_d.size());
02940     if (rst::vectorNorm(&outi_c_p_d_d[0], outi_c_p_d_d.size(), NORM_ONE) > zero) {
02941       *outStream << "\n\nINCORRECT matmatProductDataData (6): check matrix inverse property, w/ double transpose\n\n";
02942       errorFlag = -1000;
02943     }
02944     // Tensor values do not vary by point: test "N" and "T" options (no transpose/transpose)
02945     art::matmatProductDataData<double>(out_c_p_d_d, data_c_1_d_d, in_c_p_d_d);
02946     art::matmatProductDataData<double>(outi_c_p_d_d, datainv_c_1_d_d, out_c_p_d_d);
02947     rst::subtract(&outi_c_p_d_d[0], &in_c_p_d_d[0], outi_c_p_d_d.size());
02948     if (rst::vectorNorm(&outi_c_p_d_d[0], outi_c_p_d_d.size(), NORM_ONE) > zero) {
02949       *outStream << "\n\nINCORRECT matmatProductDataData (7): check matrix inverse property\n\n";
02950       errorFlag = -1000;
02951     }
02952     art::matmatProductDataData<double>(out_c_p_d_d, data_c_1_d_d, in_c_p_d_d, 't');
02953     art::matmatProductDataData<double>(outi_c_p_d_d, datainv_c_1_d_d, out_c_p_d_d, 't');
02954     rst::subtract(&outi_c_p_d_d[0], &in_c_p_d_d[0], outi_c_p_d_d.size());
02955     if (rst::vectorNorm(&outi_c_p_d_d[0], outi_c_p_d_d.size(), NORM_ONE) > zero) {
02956       *outStream << "\n\nINCORRECT matmatProductDataData (8): check matrix inverse property, w/ double transpose\n\n";
02957       errorFlag = -1000;
02958     }
02959     /***********************************************************************************************
02960       *             Full tensor: inputDataLeft(C,P,D,D) test inverse transpose                     *
02961       **********************************************************************************************/
02962     for (int i=0; i<in_c_p_d_d.size(); i++) {
02963       in_c_p_d_d[i] = Teuchos::ScalarTraits<double>::random();
02964     }
02965     for (int ic=0; ic < c; ic++) {
02966       for (int ip=0; ip < p; ip++) {
02967         for (int i=0; i<d1*d1; i++) {
02968           (&data_c_p_d_d(ic, ip, 0, 0))[i] = Teuchos::ScalarTraits<double>::random();
02969         }
02970         RealSpaceTools<double>::inverse(&datainv_c_p_d_d(ic, ip, 0, 0), &data_c_p_d_d(ic, ip, 0, 0), d1);
02971         RealSpaceTools<double>::transpose(&datainvtrn_c_p_d_d(ic, ip, 0, 0), &datainv_c_p_d_d(ic, ip, 0, 0), d1);
02972       }
02973     }
02974     for (int ic=0; ic < c; ic++) {
02975       for (int ip=0; ip < 1; ip++) {
02976         for (int i=0; i<d1*d1; i++) {
02977           (&data_c_1_d_d(ic, 0, 0, 0))[i] = Teuchos::ScalarTraits<double>::random();
02978         }
02979         RealSpaceTools<double>::inverse(&datainv_c_1_d_d(ic, 0, 0, 0), &data_c_1_d_d(ic, 0, 0, 0), d1);
02980         RealSpaceTools<double>::transpose(&datainvtrn_c_1_d_d(ic, 0, 0, 0), &datainv_c_1_d_d(ic, 0, 0, 0), d1);
02981       }
02982     }
02983     // Tensor values vary by point:
02984     art::matmatProductDataData<double>(out_c_p_d_d, data_c_p_d_d, in_c_p_d_d);
02985     art::matmatProductDataData<double>(outi_c_p_d_d, datainvtrn_c_p_d_d, out_c_p_d_d, 't');
02986     rst::subtract(&outi_c_p_d_d[0], &in_c_p_d_d[0], outi_c_p_d_d.size());
02987     if (rst::vectorNorm(&outi_c_p_d_d[0], outi_c_p_d_d.size(), NORM_ONE) > zero) {
02988       *outStream << "\n\nINCORRECT matmatProductDataData (9): check matrix inverse transpose property\n\n";
02989       errorFlag = -1000;
02990     }
02991     // Tensor values do not vary by point
02992     art::matmatProductDataData<double>(out_c_p_d_d, data_c_1_d_d, in_c_p_d_d);
02993     art::matmatProductDataData<double>(outi_c_p_d_d, datainvtrn_c_1_d_d, out_c_p_d_d, 't');
02994     rst::subtract(&outi_c_p_d_d[0], &in_c_p_d_d[0], outi_c_p_d_d.size());
02995     if (rst::vectorNorm(&outi_c_p_d_d[0], outi_c_p_d_d.size(), NORM_ONE) > zero) {
02996       *outStream << "\n\nINCORRECT matmatProductDataData (10): check matrix inverse transpose property\n\n";
02997       errorFlag = -1000;
02998     }
02999   }// test 8.a scope
03000   *outStream \
03001     << "\n"
03002     << "===============================================================================\n"\
03003     << "| TEST 8.b: matmatProductDataData random tests: branch inputDataRight(P,D,D)  |\n"\
03004     << "===============================================================================\n";
03005   /*
03006    * Test derived from Test 7.b
03007    */
03008   {// test 8.b scope
03009     int c=5, p=9, d1=3;
03010     double zero = INTREPID_TOL*10000.0;
03011     
03012     FieldContainer<double> in_p_d_d(p, d1, d1);
03013     FieldContainer<double> in_c_p_d_d(c, p, d1, d1);
03014     FieldContainer<double> out_c_p_d_d(c, p, d1, d1);
03015     FieldContainer<double> outi_c_p_d_d(c, p, d1, d1);
03016     
03017     FieldContainer<double> data_c_p(c, p);
03018     FieldContainer<double> datainv_c_p(c, p);
03019     FieldContainer<double> data_c_1(c, 1);
03020     FieldContainer<double> datainv_c_1(c, 1);
03021     FieldContainer<double> data_c_p_d(c, p, d1);
03022     FieldContainer<double> datainv_c_p_d(c, p, d1);
03023     FieldContainer<double> data_c_1_d(c, 1, d1);
03024     FieldContainer<double> datainv_c_1_d(c, 1, d1);
03025     FieldContainer<double> data_c_p_d_d(c, p, d1, d1);
03026     FieldContainer<double> datainv_c_p_d_d(c, p, d1, d1);
03027     FieldContainer<double> datainvtrn_c_p_d_d(c, p, d1, d1);
03028     FieldContainer<double> data_c_1_d_d(c, 1, d1, d1);
03029     FieldContainer<double> datainv_c_1_d_d(c, 1, d1, d1);
03030     FieldContainer<double> datainvtrn_c_1_d_d(c, 1, d1, d1);
03031     FieldContainer<double> data_c_p_one(c, p);
03032     FieldContainer<double> data_c_1_one(c, 1);
03033     /***********************************************************************************************
03034       *                          Constant diagonal tensor: inputData(C,P)                          *
03035       **********************************************************************************************/
03036     for (int i=0; i<in_p_d_d.size(); i++) {
03037       in_p_d_d[i] = Teuchos::ScalarTraits<double>::random();
03038     }
03039     for (int i=0; i<data_c_p.size(); i++) {
03040       data_c_p[i] = Teuchos::ScalarTraits<double>::random();
03041       datainv_c_p[i] = 1.0 / data_c_p[i];
03042       data_c_p_one[i] = 1.0;
03043     }
03044     for (int i=0; i<data_c_1.size(); i++) {
03045       data_c_1[i] = Teuchos::ScalarTraits<double>::random();
03046       datainv_c_1[i] = 1.0 / data_c_1[i];
03047     }
03048     // Tensor values vary by point
03049     art::matmatProductDataData<double>(in_c_p_d_d, data_c_p_one, in_p_d_d);
03050     art::matmatProductDataData<double>(out_c_p_d_d, data_c_p, in_p_d_d);
03051     art::matmatProductDataData<double>(outi_c_p_d_d, datainv_c_p, out_c_p_d_d);
03052     rst::subtract(&outi_c_p_d_d[0], &in_c_p_d_d[0], outi_c_p_d_d.size());
03053     if (rst::vectorNorm(&outi_c_p_d_d[0], outi_c_p_d_d.size(), NORM_ONE) > zero) {
03054       *outStream << "\n\nINCORRECT matmatProductDataData (11): check scalar inverse property\n\n";
03055       errorFlag = -1000;
03056     }
03057     // Tensor values do not vary by point
03058     art::matmatProductDataData<double>(in_c_p_d_d, data_c_p_one, in_p_d_d);
03059     art::matmatProductDataData<double>(out_c_p_d_d, data_c_1, in_p_d_d);
03060     art::matmatProductDataData<double>(outi_c_p_d_d, datainv_c_1, out_c_p_d_d);
03061     rst::subtract(&outi_c_p_d_d[0], &in_c_p_d_d[0], outi_c_p_d_d.size());
03062     if (rst::vectorNorm(&outi_c_p_d_d[0], outi_c_p_d_d.size(), NORM_ONE) > zero) {
03063       *outStream << "\n\nINCORRECT matmatProductDataData (12): check scalar inverse property\n\n";
03064       errorFlag = -1000;
03065     }
03066     /***********************************************************************************************
03067       *                       Non-constant diagonal tensor: inputData(C,P,D)                       *
03068       **********************************************************************************************/
03069     for (int i=0; i<in_p_d_d.size(); i++) {
03070       in_p_d_d[i] = Teuchos::ScalarTraits<double>::random();
03071     }
03072     for (int i=0; i<data_c_p_d.size(); i++) {
03073       data_c_p_d[i] = Teuchos::ScalarTraits<double>::random();
03074       datainv_c_p_d[i] = 1.0 / data_c_p_d[i];
03075     }
03076     for (int i=0; i<data_c_1_d.size(); i++) {
03077       data_c_1_d[i] = Teuchos::ScalarTraits<double>::random();
03078       datainv_c_1_d[i] = 1.0 / data_c_1_d[i];
03079     }
03080     // Tensor values vary by point:
03081     art::matmatProductDataData<double>(in_c_p_d_d, data_c_p_one, in_p_d_d);
03082     art::matmatProductDataData<double>(out_c_p_d_d, data_c_p_d, in_p_d_d);
03083     art::matmatProductDataData<double>(outi_c_p_d_d, datainv_c_p_d, out_c_p_d_d);
03084     rst::subtract(&outi_c_p_d_d[0], &in_c_p_d_d[0], outi_c_p_d_d.size());
03085     if (rst::vectorNorm(&outi_c_p_d_d[0], outi_c_p_d_d.size(), NORM_ONE) > zero) {
03086       *outStream << "\n\nINCORRECT matmatProductDataData (13): check scalar inverse property\n\n";
03087       errorFlag = -1000;
03088     }
03089     // Tensor values do not vary by point:
03090     art::matmatProductDataData<double>(in_c_p_d_d, data_c_p_one, in_p_d_d);
03091     art::matmatProductDataData<double>(out_c_p_d_d, data_c_1_d, in_p_d_d);
03092     art::matmatProductDataData<double>(outi_c_p_d_d, datainv_c_1_d, out_c_p_d_d);
03093     rst::subtract(&outi_c_p_d_d[0], &in_c_p_d_d[0], outi_c_p_d_d.size());
03094     if (rst::vectorNorm(&outi_c_p_d_d[0], outi_c_p_d_d.size(), NORM_ONE) > zero) {
03095       *outStream << "\n\nINCORRECT matmatProductDataData (14): check scalar inverse property\n\n";
03096       errorFlag = -1000;
03097     }
03098     /***********************************************************************************************
03099       *                              Full tensor: inputData(C,P,D,D)                               *
03100       **********************************************************************************************/
03101     for (int i=0; i<in_p_d_d.size(); i++) {
03102       in_p_d_d[i] = Teuchos::ScalarTraits<double>::random();
03103     }
03104     for (int ic=0; ic < c; ic++) {
03105       for (int ip=0; ip < p; ip++) {
03106         for (int i=0; i<d1*d1; i++) {
03107           (&data_c_p_d_d(ic, ip, 0, 0))[i] = Teuchos::ScalarTraits<double>::random();
03108         }
03109         RealSpaceTools<double>::inverse(&datainv_c_p_d_d(ic, ip, 0, 0), &data_c_p_d_d(ic, ip, 0, 0), d1);
03110       }
03111     }
03112     for (int ic=0; ic < c; ic++) {
03113       for (int ip=0; ip < 1; ip++) {
03114         for (int i=0; i<d1*d1; i++) {
03115           (&data_c_1_d_d(ic, 0, 0, 0))[i] = Teuchos::ScalarTraits<double>::random();
03116         }
03117         RealSpaceTools<double>::inverse(&datainv_c_1_d_d(ic, 0, 0, 0), &data_c_1_d_d(ic, 0, 0, 0), d1);
03118       }
03119     }
03120     // Tensor values vary by point: test "N" and "T" (no-transpose/transpose) options
03121     art::matmatProductDataData<double>(in_c_p_d_d, data_c_p_one, in_p_d_d);
03122     art::matmatProductDataData<double>(out_c_p_d_d, data_c_p_d_d, in_p_d_d);
03123     art::matmatProductDataData<double>(outi_c_p_d_d, datainv_c_p_d_d, out_c_p_d_d);
03124     rst::subtract(&outi_c_p_d_d[0], &in_c_p_d_d[0], outi_c_p_d_d.size());
03125     if (rst::vectorNorm(&outi_c_p_d_d[0], outi_c_p_d_d.size(), NORM_ONE) > zero) {
03126       *outStream << "\n\nINCORRECT matmatProductDataData (15): check matrix inverse property\n\n";
03127       errorFlag = -1000;
03128     }
03129     art::matmatProductDataData<double>(in_c_p_d_d, data_c_p_one, in_p_d_d);
03130     art::matmatProductDataData<double>(out_c_p_d_d, data_c_p_d_d, in_p_d_d, 't');
03131     art::matmatProductDataData<double>(outi_c_p_d_d, datainv_c_p_d_d, out_c_p_d_d, 't');
03132     rst::subtract(&outi_c_p_d_d[0], &in_c_p_d_d[0], outi_c_p_d_d.size());
03133     if (rst::vectorNorm(&outi_c_p_d_d[0], outi_c_p_d_d.size(), NORM_ONE) > zero) {
03134       *outStream << "\n\nINCORRECT matmatProductDataData (16): check matrix inverse property, w/ double transpose\n\n";
03135       errorFlag = -1000;
03136     }
03137     // Tensor values do not vary by point: test "N" and "T" (no-transpose/transpose) options
03138     art::matmatProductDataData<double>(in_c_p_d_d, data_c_p_one, in_p_d_d);
03139     art::matmatProductDataData<double>(out_c_p_d_d, data_c_1_d_d, in_p_d_d);
03140     art::matmatProductDataData<double>(outi_c_p_d_d, datainv_c_1_d_d, out_c_p_d_d);
03141     rst::subtract(&outi_c_p_d_d[0], &in_c_p_d_d[0], outi_c_p_d_d.size());
03142     if (rst::vectorNorm(&outi_c_p_d_d[0], outi_c_p_d_d.size(), NORM_ONE) > zero) {
03143       *outStream << "\n\nINCORRECT matmatProductDataData (17): check matrix inverse property\n\n";
03144       errorFlag = -1000;
03145     }
03146     art::matmatProductDataData<double>(in_c_p_d_d, data_c_p_one, in_p_d_d);
03147     art::matmatProductDataData<double>(out_c_p_d_d, data_c_1_d_d, in_p_d_d, 't');
03148     art::matmatProductDataData<double>(outi_c_p_d_d, datainv_c_1_d_d, out_c_p_d_d, 't');
03149     rst::subtract(&outi_c_p_d_d[0], &in_c_p_d_d[0], outi_c_p_d_d.size());
03150     if (rst::vectorNorm(&outi_c_p_d_d[0], outi_c_p_d_d.size(), NORM_ONE) > zero) {
03151       *outStream << "\n\nINCORRECT matmatProductDataData (18): check matrix inverse property, w/ double transpose\n\n";
03152       errorFlag = -1000;
03153     }
03154     /***********************************************************************************************
03155       *             Full tensor: inputData(C,P,D,D) test inverse transpose                         *
03156       **********************************************************************************************/
03157     for (int i=0; i<in_p_d_d.size(); i++) {
03158       in_p_d_d[i] = Teuchos::ScalarTraits<double>::random();
03159     }
03160     for (int ic=0; ic < c; ic++) {
03161       for (int ip=0; ip < p; ip++) {
03162         for (int i=0; i<d1*d1; i++) {
03163           (&data_c_p_d_d(ic, ip, 0, 0))[i] = Teuchos::ScalarTraits<double>::random();
03164         }
03165         RealSpaceTools<double>::inverse(&datainv_c_p_d_d(ic, ip, 0, 0), &data_c_p_d_d(ic, ip, 0, 0), d1);
03166         RealSpaceTools<double>::transpose(&datainvtrn_c_p_d_d(ic, ip, 0, 0), &datainv_c_p_d_d(ic, ip, 0, 0), d1);
03167       }
03168     }
03169     for (int ic=0; ic < c; ic++) {
03170       for (int ip=0; ip < 1; ip++) {
03171         for (int i=0; i<d1*d1; i++) {
03172           (&data_c_1_d_d(ic, 0, 0, 0))[i] = Teuchos::ScalarTraits<double>::random();
03173         }
03174         RealSpaceTools<double>::inverse(&datainv_c_1_d_d(ic, 0, 0, 0), &data_c_1_d_d(ic, 0, 0, 0), d1);
03175         RealSpaceTools<double>::transpose(&datainvtrn_c_1_d_d(ic, 0, 0, 0), &datainv_c_1_d_d(ic, 0, 0, 0), d1);
03176       }
03177     }
03178     // Tensor values vary by point:
03179     art::matmatProductDataData<double>(in_c_p_d_d, data_c_p_one, in_p_d_d);
03180     art::matmatProductDataData<double>(out_c_p_d_d, data_c_p_d_d, in_p_d_d);
03181     art::matmatProductDataData<double>(outi_c_p_d_d, datainvtrn_c_p_d_d, out_c_p_d_d, 't');
03182     rst::subtract(&outi_c_p_d_d[0], &in_c_p_d_d[0], outi_c_p_d_d.size());
03183     if (rst::vectorNorm(&outi_c_p_d_d[0], outi_c_p_d_d.size(), NORM_ONE) > zero) {
03184       *outStream << "\n\nINCORRECT matmatProductDataData (19): check matrix inverse transpose property\n\n";
03185       errorFlag = -1000;
03186     }
03187     // Tensor values do not vary by point:
03188     art::matmatProductDataData<double>(in_c_p_d_d, data_c_p_one, in_p_d_d);
03189     art::matmatProductDataData<double>(out_c_p_d_d, data_c_1_d_d, in_p_d_d);
03190     art::matmatProductDataData<double>(outi_c_p_d_d, datainvtrn_c_1_d_d, out_c_p_d_d, 't');
03191     rst::subtract(&outi_c_p_d_d[0], &in_c_p_d_d[0], outi_c_p_d_d.size());
03192     if (rst::vectorNorm(&outi_c_p_d_d[0], outi_c_p_d_d.size(), NORM_ONE) > zero) {
03193       *outStream << "\n\nINCORRECT matmatProductDataData (20): check matrix inverse transpose property\n\n";
03194       errorFlag = -1000;
03195     }    
03196   }// test 8.b scope
03197   }//try
03198   
03199   /*************************************************************************************************
03200     *                                      Finish test                                             *
03201     ************************************************************************************************/
03202   
03203   catch (std::logic_error err) {
03204     *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
03205     *outStream << err.what() << '\n';
03206     *outStream << "-------------------------------------------------------------------------------" << "\n\n";
03207     errorFlag = -1000;
03208   };
03209   
03210   
03211   if (errorFlag != 0)
03212     std::cout << "End Result: TEST FAILED\n";
03213   else
03214     std::cout << "End Result: TEST PASSED\n";
03215 
03216   // reset format state of std::cout
03217   std::cout.copyfmt(oldFormatState);
03218 
03219   return errorFlag;
03220 }