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