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