|
Intrepid
|
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 }
1.7.6.1