|
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 00050 namespace Intrepid { 00051 00052 00053 00054 template<class Scalar> 00055 void RealSpaceTools<Scalar>::absval(Scalar* absArray, const Scalar* inArray, const int size) { 00056 for (int i=0; i<size; i++) { 00057 absArray[i] = std::abs(inArray[i]); 00058 } 00059 } 00060 00061 00062 00063 template<class Scalar> 00064 void RealSpaceTools<Scalar>::absval(Scalar* inoutAbsArray, const int size) { 00065 for (int i=0; i<size; i++) { 00066 inoutAbsArray[i] = std::abs(inoutAbsArray[i]); 00067 } 00068 } 00069 00070 00071 00072 template<class Scalar> 00073 template<class ArrayAbs, class ArrayIn> 00074 void RealSpaceTools<Scalar>::absval(ArrayAbs & absArray, const ArrayIn & inArray) { 00075 #ifdef HAVE_INTREPID_DEBUG 00076 TEUCHOS_TEST_FOR_EXCEPTION( ( inArray.rank() != absArray.rank() ), 00077 std::invalid_argument, 00078 ">>> ERROR (RealSpaceTools::absval): Array arguments must have identical ranks!"); 00079 for (int i=0; i<inArray.rank(); i++) { 00080 TEUCHOS_TEST_FOR_EXCEPTION( ( inArray.dimension(i) != absArray.dimension(i) ), 00081 std::invalid_argument, 00082 ">>> ERROR (RealSpaceTools::absval): Dimensions of array arguments do not agree!"); 00083 } 00084 #endif 00085 00086 for (int i=0; i<inArray.size(); i++) { 00087 absArray[i] = std::abs(inArray[i]); 00088 } 00089 } 00090 00091 00092 00093 template<class Scalar> 00094 template<class ArrayInOut> 00095 void RealSpaceTools<Scalar>::absval(ArrayInOut & inoutAbsArray) { 00096 for (int i=0; i<inoutAbsArray.size(); i++) { 00097 inoutAbsArray[i] = std::abs(inoutAbsArray[i]); 00098 } 00099 } 00100 00101 00102 00103 template<class Scalar> 00104 Scalar RealSpaceTools<Scalar>::vectorNorm(const Scalar* inVec, const int dim, const ENorm normType) { 00105 Scalar temp = (Scalar)0; 00106 switch(normType) { 00107 case NORM_TWO: 00108 for(int i = 0; i < dim; i++){ 00109 temp += inVec[i]*inVec[i]; 00110 } 00111 temp = std::sqrt(temp); 00112 break; 00113 case NORM_INF: 00114 temp = std::abs(inVec[0]); 00115 for(int i = 1; i < dim; i++){ 00116 Scalar absData = std::abs(inVec[i]); 00117 if (temp < absData) temp = absData; 00118 } 00119 break; 00120 case NORM_ONE: 00121 for(int i = 0; i < dim; i++){ 00122 temp += std::abs(inVec[i]); 00123 } 00124 break; 00125 default: 00126 TEUCHOS_TEST_FOR_EXCEPTION( ( (normType != NORM_TWO) && (normType != NORM_INF) && (normType != NORM_ONE) ), 00127 std::invalid_argument, 00128 ">>> ERROR (RealSpaceTools::vectorNorm): Invalid argument normType."); 00129 } 00130 return temp; 00131 } 00132 00133 00134 00135 template<class Scalar> 00136 template<class ArrayIn> 00137 Scalar RealSpaceTools<Scalar>::vectorNorm(const ArrayIn & inVec, const ENorm normType) { 00138 00139 #ifdef HAVE_INTREPID_DEBUG 00140 TEUCHOS_TEST_FOR_EXCEPTION( ( inVec.rank() != 1 ), 00141 std::invalid_argument, 00142 ">>> ERROR (RealSpaceTools::vectorNorm): Vector argument must have rank 1!"); 00143 #endif 00144 00145 int size = inVec.size(); 00146 00147 Scalar temp = (Scalar)0; 00148 switch(normType) { 00149 case NORM_TWO: 00150 for(int i = 0; i < size; i++){ 00151 temp += inVec[i]*inVec[i]; 00152 } 00153 temp = std::sqrt(temp); 00154 break; 00155 case NORM_INF: 00156 temp = std::abs(inVec[0]); 00157 for(int i = 1; i < size; i++){ 00158 Scalar absData = std::abs(inVec[i]); 00159 if (temp < absData) temp = absData; 00160 } 00161 break; 00162 case NORM_ONE: 00163 for(int i = 0; i < size; i++){ 00164 temp += std::abs(inVec[i]); 00165 } 00166 break; 00167 default: 00168 TEUCHOS_TEST_FOR_EXCEPTION( ( (normType != NORM_TWO) && (normType != NORM_INF) && (normType != NORM_ONE) ), 00169 std::invalid_argument, 00170 ">>> ERROR (RealSpaceTools::vectorNorm): Invalid argument normType."); 00171 } 00172 return temp; 00173 } 00174 00175 00176 00177 template<class Scalar> 00178 template<class ArrayNorm, class ArrayIn> 00179 void RealSpaceTools<Scalar>::vectorNorm(ArrayNorm & normArray, const ArrayIn & inVecs, const ENorm normType) { 00180 00181 int arrayRank = inVecs.rank(); 00182 00183 #ifdef HAVE_INTREPID_DEBUG 00184 TEUCHOS_TEST_FOR_EXCEPTION( ( arrayRank != normArray.rank()+1 ), 00185 std::invalid_argument, 00186 ">>> ERROR (RealSpaceTools::vectorNorm): Ranks of norm and vector array arguments are incompatible!"); 00187 TEUCHOS_TEST_FOR_EXCEPTION( ( (arrayRank < 2) || (arrayRank > 3) ), 00188 std::invalid_argument, 00189 ">>> ERROR (RealSpaceTools::vectorNorm): Rank of vector array must be 2 or 3!"); 00190 for (int i=0; i<arrayRank-1; i++) { 00191 TEUCHOS_TEST_FOR_EXCEPTION( ( inVecs.dimension(i) != normArray.dimension(i) ), 00192 std::invalid_argument, 00193 ">>> ERROR (RealSpaceTools::vectorNorm): Dimensions of norm and vector arguments do not agree!"); 00194 } 00195 #endif 00196 00197 int dim_i0 = 1; // first index dimension (e.g. cell index) 00198 int dim_i1 = 1; // second index dimension (e.g. point index) 00199 int dim = inVecs.dimension(arrayRank-1); // spatial dimension 00200 00201 // determine i0 and i1 dimensions 00202 switch(arrayRank) { 00203 case 3: 00204 dim_i0 = inVecs.dimension(0); 00205 dim_i1 = inVecs.dimension(1); 00206 break; 00207 case 2: 00208 dim_i1 = inVecs.dimension(0); 00209 break; 00210 } 00211 00212 switch(normType) { 00213 case NORM_TWO: { 00214 int offset_i0, offset, normOffset; 00215 for (int i0=0; i0<dim_i0; i0++) { 00216 offset_i0 = i0*dim_i1; 00217 for (int i1=0; i1<dim_i1; i1++) { 00218 offset = offset_i0 + i1; 00219 normOffset = offset; 00220 offset *= dim; 00221 Scalar temp = (Scalar)0; 00222 for(int i = 0; i < dim; i++){ 00223 temp += inVecs[offset+i]*inVecs[offset+i]; 00224 } 00225 normArray[normOffset] = std::sqrt(temp); 00226 } 00227 } 00228 break; 00229 } // case NORM_TWO 00230 00231 case NORM_INF: { 00232 int offset_i0, offset, normOffset; 00233 for (int i0=0; i0<dim_i0; i0++) { 00234 offset_i0 = i0*dim_i1; 00235 for (int i1=0; i1<dim_i1; i1++) { 00236 offset = offset_i0 + i1; 00237 normOffset = offset; 00238 offset *= dim; 00239 Scalar temp = (Scalar)0; 00240 temp = std::abs(inVecs[offset]); 00241 for(int i = 1; i < dim; i++){ 00242 Scalar absData = std::abs(inVecs[offset+i]); 00243 if (temp < absData) temp = absData; 00244 } 00245 normArray[normOffset] = temp; 00246 } 00247 } 00248 break; 00249 } // case NORM_INF 00250 00251 case NORM_ONE: { 00252 int offset_i0, offset, normOffset; 00253 for (int i0=0; i0<dim_i0; i0++) { 00254 offset_i0 = i0*dim_i1; 00255 for (int i1=0; i1<dim_i1; i1++) { 00256 offset = offset_i0 + i1; 00257 normOffset = offset; 00258 offset *= dim; 00259 Scalar temp = (Scalar)0; 00260 for(int i = 0; i < dim; i++){ 00261 temp += std::abs(inVecs[offset+i]); 00262 } 00263 normArray[normOffset] = temp; 00264 } 00265 } 00266 break; 00267 } // case NORM_ONE 00268 00269 default: 00270 TEUCHOS_TEST_FOR_EXCEPTION( ( (normType != NORM_TWO) && (normType != NORM_INF) && (normType != NORM_ONE) ), 00271 std::invalid_argument, 00272 ">>> ERROR (RealSpaceTools::vectorNorm): Invalid argument normType."); 00273 } 00274 } 00275 00276 00277 00278 template<class Scalar> 00279 void RealSpaceTools<Scalar>::transpose(Scalar* transposeMat, const Scalar* inMat, const int dim) { 00280 for(int i=0; i < dim; i++){ 00281 transposeMat[i*dim+i]=inMat[i*dim+i]; // Set diagonal elements 00282 for(int j=i+1; j < dim; j++){ 00283 transposeMat[i*dim+j]=inMat[j*dim+i]; // Set off-diagonal elements 00284 transposeMat[j*dim+i]=inMat[i*dim+j]; 00285 } 00286 } 00287 } 00288 00289 00290 00291 template<class Scalar> 00292 template<class ArrayTranspose, class ArrayIn> 00293 void RealSpaceTools<Scalar>::transpose(ArrayTranspose & transposeMats, const ArrayIn & inMats) { 00294 int arrayRank = inMats.rank(); 00295 00296 #ifdef HAVE_INTREPID_DEBUG 00297 TEUCHOS_TEST_FOR_EXCEPTION( ( arrayRank != transposeMats.rank() ), 00298 std::invalid_argument, 00299 ">>> ERROR (RealSpaceTools::transpose): Matrix array arguments do not have identical ranks!"); 00300 TEUCHOS_TEST_FOR_EXCEPTION( ( (arrayRank < 2) || (arrayRank > 4) ), 00301 std::invalid_argument, 00302 ">>> ERROR (RealSpaceTools::transpose): Rank of matrix array must be 2, 3, or 4!"); 00303 for (int i=0; i<arrayRank; i++) { 00304 TEUCHOS_TEST_FOR_EXCEPTION( ( inMats.dimension(i) != transposeMats.dimension(i) ), 00305 std::invalid_argument, 00306 ">>> ERROR (RealSpaceTools::transpose): Dimensions of matrix arguments do not agree!"); 00307 } 00308 TEUCHOS_TEST_FOR_EXCEPTION( ( inMats.dimension(arrayRank-2) != inMats.dimension(arrayRank-1) ), 00309 std::invalid_argument, 00310 ">>> ERROR (RealSpaceTools::transpose): Matrices are not square!"); 00311 #endif 00312 00313 int dim_i0 = 1; // first index dimension (e.g. cell index) 00314 int dim_i1 = 1; // second index dimension (e.g. point index) 00315 int dim = inMats.dimension(arrayRank-2); // spatial dimension 00316 00317 // determine i0 and i1 dimensions 00318 switch(arrayRank) { 00319 case 4: 00320 dim_i0 = inMats.dimension(0); 00321 dim_i1 = inMats.dimension(1); 00322 break; 00323 case 3: 00324 dim_i1 = inMats.dimension(0); 00325 break; 00326 } 00327 00328 int offset_i0, offset; 00329 00330 for (int i0=0; i0<dim_i0; i0++) { 00331 offset_i0 = i0*dim_i1; 00332 for (int i1=0; i1<dim_i1; i1++) { 00333 offset = offset_i0 + i1; 00334 offset *= (dim*dim); 00335 00336 for(int i=0; i < dim; i++){ 00337 transposeMats[offset+i*dim+i]=inMats[offset+i*dim+i]; // Set diagonal elements 00338 for(int j=i+1; j < dim; j++){ 00339 transposeMats[offset+i*dim+j]=inMats[offset+j*dim+i]; // Set off-diagonal elements 00340 transposeMats[offset+j*dim+i]=inMats[offset+i*dim+j]; 00341 } 00342 } 00343 00344 } // i1 00345 } // i0 00346 00347 } 00348 00349 00350 00351 template<class Scalar> 00352 void RealSpaceTools<Scalar>::inverse(Scalar* inverseMat, const Scalar* inMat, const int dim) { 00353 00354 switch(dim) { 00355 case 3: { 00356 int i, j, rowID = 0, colID = 0; 00357 int rowperm[3]={0,1,2}; 00358 int colperm[3]={0,1,2}; // Complete pivoting 00359 Scalar emax(0); 00360 00361 for(i=0; i < 3; ++i){ 00362 for(j=0; j < 3; ++j){ 00363 if( std::abs( inMat[i*3+j] ) > emax){ 00364 rowID = i; colID = j; emax = std::abs( inMat[i*3+j] ); 00365 } 00366 } 00367 } 00368 #ifdef HAVE_INTREPID_DEBUG 00369 TEUCHOS_TEST_FOR_EXCEPTION( ( emax == (Scalar)0 ), 00370 std::invalid_argument, 00371 ">>> ERROR (Matrix): Inverse of a zero matrix is undefined!"); 00372 #endif 00373 if( rowID ){ 00374 rowperm[0] = rowID; 00375 rowperm[rowID] = 0; 00376 } 00377 if( colID ){ 00378 colperm[0] = colID; 00379 colperm[colID] = 0; 00380 } 00381 Scalar B[3][3], S[2][2], Bi[3][3]; // B=rowperm inMat colperm, S=Schur complement(Boo) 00382 for(i=0; i < 3; ++i){ 00383 for(j=0; j < 3; ++j){ 00384 B[i][j] = inMat[rowperm[i]*3+colperm[j]]; 00385 } 00386 } 00387 B[1][0] /= B[0][0]; B[2][0] /= B[0][0];// B(:,0)/=pivot 00388 for(i=0; i < 2; ++i){ 00389 for(j=0; j < 2; ++j){ 00390 S[i][j] = B[i+1][j+1] - B[i+1][0] * B[0][j+1]; // S = B -z*y' 00391 } 00392 } 00393 Scalar detS = S[0][0]*S[1][1]- S[0][1]*S[1][0], Si[2][2]; 00394 #ifdef HAVE_INTREPID_DEBUG 00395 TEUCHOS_TEST_FOR_EXCEPTION( ( detS == (Scalar)0 ), 00396 std::invalid_argument, 00397 ">>> ERROR (Matrix): Inverse of a singular matrix is undefined!"); 00398 #endif 00399 00400 Si[0][0] = S[1][1]/detS; Si[0][1] = -S[0][1]/detS; 00401 Si[1][0] = -S[1][0]/detS; Si[1][1] = S[0][0]/detS; 00402 00403 for(j=0; j<2;j++) 00404 Bi[0][j+1] = -( B[0][1]*Si[0][j] + B[0][2]* Si[1][j])/B[0][0]; 00405 for(i=0; i<2;i++) 00406 Bi[i+1][0] = -(Si[i][0]*B[1][0] + Si[i][1]*B[2][0]); 00407 00408 Bi[0][0] = ((Scalar)1/B[0][0])-Bi[0][1]*B[1][0]-Bi[0][2]*B[2][0]; 00409 Bi[1][1] = Si[0][0]; 00410 Bi[1][2] = Si[0][1]; 00411 Bi[2][1] = Si[1][0]; 00412 Bi[2][2] = Si[1][1]; 00413 for(i=0; i < 3; ++i){ 00414 for(j=0; j < 3; ++j){ 00415 inverseMat[i*3+j] = Bi[colperm[i]][rowperm[j]]; // set inverse 00416 } 00417 } 00418 break; 00419 } // case 3 00420 00421 case 2: { 00422 00423 Scalar determinant = inMat[0]*inMat[3]-inMat[1]*inMat[2];; 00424 #ifdef HAVE_INTREPID_DEBUG 00425 TEUCHOS_TEST_FOR_EXCEPTION( ( (inMat[0]==(Scalar)0) && (inMat[1]==(Scalar)0) && 00426 (inMat[2]==(Scalar)0) && (inMat[3]==(Scalar)0) ), 00427 std::invalid_argument, 00428 ">>> ERROR (Matrix): Inverse of a zero matrix is undefined!"); 00429 TEUCHOS_TEST_FOR_EXCEPTION( ( determinant == (Scalar)0 ), 00430 std::invalid_argument, 00431 ">>> ERROR (Matrix): Inverse of a singular matrix is undefined!"); 00432 #endif 00433 inverseMat[0] = inMat[3] / determinant; 00434 inverseMat[1] = - inMat[1] / determinant; 00435 // 00436 inverseMat[2] = - inMat[2] / determinant; 00437 inverseMat[3] = inMat[0] / determinant; 00438 break; 00439 } // case 2 00440 00441 case 1: { 00442 #ifdef HAVE_INTREPID_DEBUG 00443 TEUCHOS_TEST_FOR_EXCEPTION( ( inMat[0] == (Scalar)0 ), 00444 std::invalid_argument, 00445 ">>> ERROR (Matrix): Inverse of a zero matrix is undefined!"); 00446 #endif 00447 inverseMat[0] = (Scalar)1 / inMat[0]; 00448 break; 00449 } // case 1 00450 00451 } // switch (dim) 00452 } 00453 00454 00455 00456 template<class Scalar> 00457 template<class ArrayInverse, class ArrayIn> 00458 void RealSpaceTools<Scalar>::inverse(ArrayInverse & inverseMats, const ArrayIn & inMats) { 00459 00460 int arrayRank = inMats.rank(); 00461 00462 #ifdef HAVE_INTREPID_DEBUG 00463 TEUCHOS_TEST_FOR_EXCEPTION( ( arrayRank != inverseMats.rank() ), 00464 std::invalid_argument, 00465 ">>> ERROR (RealSpaceTools::inverse): Matrix array arguments do not have identical ranks!"); 00466 TEUCHOS_TEST_FOR_EXCEPTION( ( (arrayRank < 2) || (arrayRank > 4) ), 00467 std::invalid_argument, 00468 ">>> ERROR (RealSpaceTools::inverse): Rank of matrix array must be 2, 3, or 4!"); 00469 for (int i=0; i<arrayRank; i++) { 00470 TEUCHOS_TEST_FOR_EXCEPTION( ( inMats.dimension(i) != inverseMats.dimension(i) ), 00471 std::invalid_argument, 00472 ">>> ERROR (RealSpaceTools::inverse): Dimensions of matrix arguments do not agree!"); 00473 } 00474 TEUCHOS_TEST_FOR_EXCEPTION( ( inMats.dimension(arrayRank-2) != inMats.dimension(arrayRank-1) ), 00475 std::invalid_argument, 00476 ">>> ERROR (RealSpaceTools::inverse): Matrices are not square!"); 00477 TEUCHOS_TEST_FOR_EXCEPTION( ( (inMats.dimension(arrayRank-2) < 1) || (inMats.dimension(arrayRank-2) > 3) ), 00478 std::invalid_argument, 00479 ">>> ERROR (RealSpaceTools::inverse): Spatial dimension must be 1, 2, or 3!"); 00480 #endif 00481 00482 int dim_i0 = 1; // first index dimension (e.g. cell index) 00483 int dim_i1 = 1; // second index dimension (e.g. point index) 00484 int dim = inMats.dimension(arrayRank-2); // spatial dimension 00485 00486 // determine i0 and i1 dimensions 00487 switch(arrayRank) { 00488 case 4: 00489 dim_i0 = inMats.dimension(0); 00490 dim_i1 = inMats.dimension(1); 00491 break; 00492 case 3: 00493 dim_i1 = inMats.dimension(0); 00494 break; 00495 } 00496 00497 switch(dim) { 00498 case 3: { 00499 int offset_i0, offset; 00500 00501 for (int i0=0; i0<dim_i0; i0++) { 00502 offset_i0 = i0*dim_i1; 00503 for (int i1=0; i1<dim_i1; i1++) { 00504 offset = offset_i0 + i1; 00505 offset *= 9; 00506 00507 int i, j, rowID = 0, colID = 0; 00508 int rowperm[3]={0,1,2}; 00509 int colperm[3]={0,1,2}; // Complete pivoting 00510 Scalar emax(0); 00511 00512 for(i=0; i < 3; ++i){ 00513 for(j=0; j < 3; ++j){ 00514 if( std::abs( inMats[offset+i*3+j] ) > emax){ 00515 rowID = i; colID = j; emax = std::abs( inMats[offset+i*3+j] ); 00516 } 00517 } 00518 } 00519 #ifdef HAVE_INTREPID_DEBUG 00520 #ifdef HAVE_INTREPID_DEBUG_INF_CHECK 00521 TEUCHOS_TEST_FOR_EXCEPTION( ( emax == (Scalar)0 ), 00522 std::invalid_argument, 00523 ">>> ERROR (Matrix): Inverse of a zero matrix is undefined!"); 00524 #endif 00525 #endif 00526 if( rowID ){ 00527 rowperm[0] = rowID; 00528 rowperm[rowID] = 0; 00529 } 00530 if( colID ){ 00531 colperm[0] = colID; 00532 colperm[colID] = 0; 00533 } 00534 Scalar B[3][3], S[2][2], Bi[3][3]; // B=rowperm inMat colperm, S=Schur complement(Boo) 00535 for(i=0; i < 3; ++i){ 00536 for(j=0; j < 3; ++j){ 00537 B[i][j] = inMats[offset+rowperm[i]*3+colperm[j]]; 00538 } 00539 } 00540 B[1][0] /= B[0][0]; B[2][0] /= B[0][0];// B(:,0)/=pivot 00541 for(i=0; i < 2; ++i){ 00542 for(j=0; j < 2; ++j){ 00543 S[i][j] = B[i+1][j+1] - B[i+1][0] * B[0][j+1]; // S = B -z*y' 00544 } 00545 } 00546 Scalar detS = S[0][0]*S[1][1]- S[0][1]*S[1][0], Si[2][2]; 00547 #ifdef HAVE_INTREPID_DEBUG 00548 #ifdef HAVE_INTREPID_DEBUG_INF_CHECK 00549 TEUCHOS_TEST_FOR_EXCEPTION( ( detS == (Scalar)0 ), 00550 std::invalid_argument, 00551 ">>> ERROR (Matrix): Inverse of a singular matrix is undefined!"); 00552 #endif 00553 #endif 00554 00555 Si[0][0] = S[1][1]/detS; Si[0][1] = -S[0][1]/detS; 00556 Si[1][0] = -S[1][0]/detS; Si[1][1] = S[0][0]/detS; 00557 00558 for(j=0; j<2;j++) 00559 Bi[0][j+1] = -( B[0][1]*Si[0][j] + B[0][2]* Si[1][j])/B[0][0]; 00560 for(i=0; i<2;i++) 00561 Bi[i+1][0] = -(Si[i][0]*B[1][0] + Si[i][1]*B[2][0]); 00562 00563 Bi[0][0] = ((Scalar)1/B[0][0])-Bi[0][1]*B[1][0]-Bi[0][2]*B[2][0]; 00564 Bi[1][1] = Si[0][0]; 00565 Bi[1][2] = Si[0][1]; 00566 Bi[2][1] = Si[1][0]; 00567 Bi[2][2] = Si[1][1]; 00568 for(i=0; i < 3; ++i){ 00569 for(j=0; j < 3; ++j){ 00570 inverseMats[offset+i*3+j] = Bi[colperm[i]][rowperm[j]]; // set inverse 00571 } 00572 } 00573 } // for i1 00574 } // for i0 00575 break; 00576 } // case 3 00577 00578 case 2: { 00579 int offset_i0, offset; 00580 00581 for (int i0=0; i0<dim_i0; i0++) { 00582 offset_i0 = i0*dim_i1; 00583 for (int i1=0; i1<dim_i1; i1++) { 00584 offset = offset_i0 + i1;; 00585 offset *= 4; 00586 00587 Scalar determinant = inMats[offset]*inMats[offset+3]-inMats[offset+1]*inMats[offset+2];; 00588 #ifdef HAVE_INTREPID_DEBUG 00589 #ifdef HAVE_INTREPID_DEBUG_INF_CHECK 00590 TEUCHOS_TEST_FOR_EXCEPTION( ( (inMats[offset]==(Scalar)0) && (inMats[offset+1]==(Scalar)0) && 00591 (inMats[offset+2]==(Scalar)0) && (inMats[offset+3]==(Scalar)0) ), 00592 std::invalid_argument, 00593 ">>> ERROR (Matrix): Inverse of a zero matrix is undefined!"); 00594 TEUCHOS_TEST_FOR_EXCEPTION( ( determinant == (Scalar)0 ), 00595 std::invalid_argument, 00596 ">>> ERROR (Matrix): Inverse of a singular matrix is undefined!"); 00597 #endif 00598 #endif 00599 inverseMats[offset] = inMats[offset+3] / determinant; 00600 inverseMats[offset+1] = - inMats[offset+1] / determinant; 00601 // 00602 inverseMats[offset+2] = - inMats[offset+2] / determinant; 00603 inverseMats[offset+3] = inMats[offset] / determinant; 00604 } // for i1 00605 } // for i0 00606 break; 00607 } // case 2 00608 00609 case 1: { 00610 int offset_i0, offset; 00611 00612 for (int i0=0; i0<dim_i0; i0++) { 00613 offset_i0 = i0*dim_i1; 00614 for (int i1=0; i1<dim_i1; i1++) { 00615 offset = offset_i0 + i1;; 00616 #ifdef HAVE_INTREPID_DEBUG 00617 #ifdef HAVE_INTREPID_DEBUG_INF_CHECK 00618 TEUCHOS_TEST_FOR_EXCEPTION( ( inMats[offset] == (Scalar)0 ), 00619 std::invalid_argument, 00620 ">>> ERROR (Matrix): Inverse of a zero matrix is undefined!"); 00621 #endif 00622 #endif 00623 inverseMats[offset] = (Scalar)1 / inMats[offset]; 00624 } // for i1 00625 } // for i2 00626 break; 00627 } // case 1 00628 00629 } // switch (dim) 00630 } 00631 00632 00633 00634 template<class Scalar> 00635 Scalar RealSpaceTools<Scalar>::det(const Scalar* inMat, const int dim) { 00636 Scalar determinant(0); 00637 00638 switch (dim) { 00639 case 3: { 00640 int i,j,rowID = 0; 00641 int colID = 0; 00642 int rowperm[3]={0,1,2}; 00643 int colperm[3]={0,1,2}; // Complete pivoting 00644 Scalar emax(0); 00645 00646 for(i=0; i < 3; ++i){ 00647 for(j=0; j < 3; ++j){ 00648 if( std::abs( inMat[i*dim+j] ) > emax){ 00649 rowID = i; colID = j; emax = std::abs( inMat[i*dim+j] ); 00650 } 00651 } 00652 } 00653 if( emax > 0 ){ 00654 if( rowID ){ 00655 rowperm[0] = rowID; 00656 rowperm[rowID] = 0; 00657 } 00658 if( colID ){ 00659 colperm[0] = colID; 00660 colperm[colID] = 0; 00661 } 00662 Scalar B[3][3], S[2][2]; // B=rowperm inMat colperm, S=Schur complement(Boo) 00663 for(i=0; i < 3; ++i){ 00664 for(j=0; j < 3; ++j){ 00665 B[i][j] = inMat[rowperm[i]*dim+colperm[j]]; 00666 } 00667 } 00668 B[1][0] /= B[0][0]; B[2][0] /= B[0][0];// B(:,0)/=pivot 00669 for(i=0; i < 2; ++i){ 00670 for(j=0; j < 2; ++j){ 00671 S[i][j] = B[i+1][j+1] - B[i+1][0] * B[0][j+1]; // S = B -z*y' 00672 } 00673 } 00674 determinant = B[0][0] * (S[0][0] * S[1][1] - S[0][1] * S[1][0]); // det(B) 00675 if( rowID ) determinant = -determinant; 00676 if( colID ) determinant = -determinant; 00677 } 00678 break; 00679 } // case 3 00680 00681 case 2: 00682 determinant = inMat[0]*inMat[3]- 00683 inMat[1]*inMat[2]; 00684 break; 00685 00686 case 1: 00687 determinant = inMat[0]; 00688 break; 00689 00690 default: 00691 TEUCHOS_TEST_FOR_EXCEPTION( ( (dim != 1) && (dim != 2) && (dim != 3) ), 00692 std::invalid_argument, 00693 ">>> ERROR (Matrix): Invalid matrix dimension."); 00694 } // switch (dim) 00695 00696 return determinant; 00697 } 00698 00699 00700 00701 template<class Scalar> 00702 template<class ArrayIn> 00703 Scalar RealSpaceTools<Scalar>::det(const ArrayIn & inMat) { 00704 00705 #ifdef HAVE_INTREPID_DEBUG 00706 TEUCHOS_TEST_FOR_EXCEPTION( (inMat.rank() != 2), 00707 std::invalid_argument, 00708 ">>> ERROR (RealSpaceTools::det): Rank of matrix argument must be 2!"); 00709 TEUCHOS_TEST_FOR_EXCEPTION( ( inMat.dimension(0) != inMat.dimension(1) ), 00710 std::invalid_argument, 00711 ">>> ERROR (RealSpaceTools::det): Matrix is not square!"); 00712 TEUCHOS_TEST_FOR_EXCEPTION( ( (inMat.dimension(0) < 1) || (inMat.dimension(0) > 3) ), 00713 std::invalid_argument, 00714 ">>> ERROR (RealSpaceTools::det): Spatial dimension must be 1, 2, or 3!"); 00715 #endif 00716 00717 int dim = inMat.dimension(0); 00718 Scalar determinant(0); 00719 00720 switch (dim) { 00721 case 3: { 00722 int i,j,rowID = 0; 00723 int colID = 0; 00724 int rowperm[3]={0,1,2}; 00725 int colperm[3]={0,1,2}; // Complete pivoting 00726 Scalar emax(0); 00727 00728 for(i=0; i < 3; ++i){ 00729 for(j=0; j < 3; ++j){ 00730 if( std::abs( inMat[i*dim+j] ) > emax){ 00731 rowID = i; colID = j; emax = std::abs( inMat[i*dim+j] ); 00732 } 00733 } 00734 } 00735 if( emax > 0 ){ 00736 if( rowID ){ 00737 rowperm[0] = rowID; 00738 rowperm[rowID] = 0; 00739 } 00740 if( colID ){ 00741 colperm[0] = colID; 00742 colperm[colID] = 0; 00743 } 00744 Scalar B[3][3], S[2][2]; // B=rowperm inMat colperm, S=Schur complement(Boo) 00745 for(i=0; i < 3; ++i){ 00746 for(j=0; j < 3; ++j){ 00747 B[i][j] = inMat[rowperm[i]*dim+colperm[j]]; 00748 } 00749 } 00750 B[1][0] /= B[0][0]; B[2][0] /= B[0][0];// B(:,0)/=pivot 00751 for(i=0; i < 2; ++i){ 00752 for(j=0; j < 2; ++j){ 00753 S[i][j] = B[i+1][j+1] - B[i+1][0] * B[0][j+1]; // S = B -z*y' 00754 } 00755 } 00756 determinant = B[0][0] * (S[0][0] * S[1][1] - S[0][1] * S[1][0]); // det(B) 00757 if( rowID ) determinant = -determinant; 00758 if( colID ) determinant = -determinant; 00759 } 00760 break; 00761 } // case 3 00762 00763 case 2: 00764 determinant = inMat[0]*inMat[3]- 00765 inMat[1]*inMat[2]; 00766 break; 00767 00768 case 1: 00769 determinant = inMat[0]; 00770 break; 00771 00772 default: 00773 TEUCHOS_TEST_FOR_EXCEPTION( ( (dim != 1) && (dim != 2) && (dim != 3) ), 00774 std::invalid_argument, 00775 ">>> ERROR (Matrix): Invalid matrix dimension."); 00776 } // switch (dim) 00777 00778 return determinant; 00779 } 00780 00781 00782 00783 00784 template<class Scalar> 00785 template<class ArrayDet, class ArrayIn> 00786 void RealSpaceTools<Scalar>::det(ArrayDet & detArray, const ArrayIn & inMats) { 00787 00788 int matArrayRank = inMats.rank(); 00789 00790 #ifdef HAVE_INTREPID_DEBUG 00791 TEUCHOS_TEST_FOR_EXCEPTION( ( matArrayRank != detArray.rank()+2 ), 00792 std::invalid_argument, 00793 ">>> ERROR (RealSpaceTools::det): Determinant and matrix array arguments do not have compatible ranks!"); 00794 TEUCHOS_TEST_FOR_EXCEPTION( ( (matArrayRank < 3) || (matArrayRank > 4) ), 00795 std::invalid_argument, 00796 ">>> ERROR (RealSpaceTools::det): Rank of matrix array must be 3 or 4!"); 00797 for (int i=0; i<matArrayRank-2; i++) { 00798 TEUCHOS_TEST_FOR_EXCEPTION( ( inMats.dimension(i) != detArray.dimension(i) ), 00799 std::invalid_argument, 00800 ">>> ERROR (RealSpaceTools::det): Dimensions of determinant and matrix array arguments do not agree!"); 00801 } 00802 TEUCHOS_TEST_FOR_EXCEPTION( ( inMats.dimension(matArrayRank-2) != inMats.dimension(matArrayRank-1) ), 00803 std::invalid_argument, 00804 ">>> ERROR (RealSpaceTools::det): Matrices are not square!"); 00805 TEUCHOS_TEST_FOR_EXCEPTION( ( (inMats.dimension(matArrayRank-2) < 1) || (inMats.dimension(matArrayRank-2) > 3) ), 00806 std::invalid_argument, 00807 ">>> ERROR (RealSpaceTools::det): Spatial dimension must be 1, 2, or 3!"); 00808 #endif 00809 00810 int dim_i0 = 1; // first index dimension (e.g. cell index) 00811 int dim_i1 = 1; // second index dimension (e.g. point index) 00812 int dim = inMats.dimension(matArrayRank-2); // spatial dimension 00813 00814 // determine i0 and i1 dimensions 00815 switch(matArrayRank) { 00816 case 4: 00817 dim_i0 = inMats.dimension(0); 00818 dim_i1 = inMats.dimension(1); 00819 break; 00820 case 3: 00821 dim_i1 = inMats.dimension(0); 00822 break; 00823 } 00824 00825 switch(dim) { 00826 case 3: { 00827 int offset_i0, offset, detOffset; 00828 00829 for (int i0=0; i0<dim_i0; i0++) { 00830 offset_i0 = i0*dim_i1; 00831 for (int i1=0; i1<dim_i1; i1++) { 00832 offset = offset_i0 + i1; 00833 detOffset = offset; 00834 offset *= 9; 00835 00836 int i,j,rowID = 0; 00837 int colID = 0; 00838 int rowperm[3]={0,1,2}; 00839 int colperm[3]={0,1,2}; // Complete pivoting 00840 Scalar emax(0), determinant(0); 00841 00842 for(i=0; i < 3; ++i){ 00843 for(j=0; j < 3; ++j){ 00844 if( std::abs( inMats[offset+i*3+j] ) > emax){ 00845 rowID = i; colID = j; emax = std::abs( inMats[offset+i*3+j] ); 00846 } 00847 } 00848 } 00849 if( emax > 0 ){ 00850 if( rowID ){ 00851 rowperm[0] = rowID; 00852 rowperm[rowID] = 0; 00853 } 00854 if( colID ){ 00855 colperm[0] = colID; 00856 colperm[colID] = 0; 00857 } 00858 Scalar B[3][3], S[2][2]; // B=rowperm inMat colperm, S=Schur complement(Boo) 00859 for(i=0; i < 3; ++i){ 00860 for(j=0; j < 3; ++j){ 00861 B[i][j] = inMats[offset+rowperm[i]*3+colperm[j]]; 00862 } 00863 } 00864 B[1][0] /= B[0][0]; B[2][0] /= B[0][0];// B(:,0)/=pivot 00865 for(i=0; i < 2; ++i){ 00866 for(j=0; j < 2; ++j){ 00867 S[i][j] = B[i+1][j+1] - B[i+1][0] * B[0][j+1]; // S = B -z*y' 00868 } 00869 } 00870 determinant = B[0][0] * (S[0][0] * S[1][1] - S[0][1] * S[1][0]); // det(B) 00871 if( rowID ) determinant = -determinant; 00872 if( colID ) determinant = -determinant; 00873 } 00874 detArray[detOffset] = determinant; 00875 } // for i1 00876 } // for i0 00877 break; 00878 } // case 3 00879 00880 case 2: { 00881 int offset_i0, offset, detOffset; 00882 00883 for (int i0=0; i0<dim_i0; i0++) { 00884 offset_i0 = i0*dim_i1; 00885 for (int i1=0; i1<dim_i1; i1++) { 00886 offset = offset_i0 + i1; 00887 detOffset = offset; 00888 offset *= 4; 00889 00890 detArray[detOffset] = inMats[offset]*inMats[offset+3]-inMats[offset+1]*inMats[offset+2];; 00891 } // for i1 00892 } // for i0 00893 break; 00894 } // case 2 00895 00896 case 1: { 00897 int offset_i0, offset; 00898 00899 for (int i0=0; i0<dim_i0; i0++) { 00900 offset_i0 = i0*dim_i1; 00901 for (int i1=0; i1<dim_i1; i1++) { 00902 offset = offset_i0 + i1;; 00903 detArray[offset] = inMats[offset]; 00904 } // for i1 00905 } // for i2 00906 break; 00907 } // case 1 00908 00909 } // switch (dim) 00910 } 00911 00912 00913 00914 template<class Scalar> 00915 void RealSpaceTools<Scalar>::add(Scalar* sumArray, const Scalar* inArray1, const Scalar* inArray2, const int size) { 00916 for (int i=0; i<size; i++) { 00917 sumArray[i] = inArray1[i] + inArray2[i]; 00918 } 00919 } 00920 00921 00922 00923 template<class Scalar> 00924 void RealSpaceTools<Scalar>::add(Scalar* inoutSumArray, const Scalar* inArray, const int size) { 00925 for (int i=0; i<size; i++) { 00926 inoutSumArray[i] += inArray[i]; 00927 } 00928 } 00929 00930 00931 00932 template<class Scalar> 00933 template<class ArraySum, class ArrayIn1, class ArrayIn2> 00934 void RealSpaceTools<Scalar>::add(ArraySum & sumArray, const ArrayIn1 & inArray1, const ArrayIn2 & inArray2) { 00935 #ifdef HAVE_INTREPID_DEBUG 00936 TEUCHOS_TEST_FOR_EXCEPTION( ( (inArray1.rank() != inArray2.rank()) || (inArray1.rank() != sumArray.rank()) ), 00937 std::invalid_argument, 00938 ">>> ERROR (RealSpaceTools::add): Array arguments must have identical ranks!"); 00939 for (int i=0; i<inArray1.rank(); i++) { 00940 TEUCHOS_TEST_FOR_EXCEPTION( ( (inArray1.dimension(i) != inArray2.dimension(i)) || (inArray1.dimension(i) != sumArray.dimension(i)) ), 00941 std::invalid_argument, 00942 ">>> ERROR (RealSpaceTools::add): Dimensions of array arguments do not agree!"); 00943 } 00944 #endif 00945 00946 for (int i=0; i<inArray1.size(); i++) { 00947 sumArray[i] = inArray1[i] + inArray2[i]; 00948 } 00949 } 00950 00951 00952 00953 template<class Scalar> 00954 template<class ArraySum, class ArrayIn> 00955 void RealSpaceTools<Scalar>::add(ArraySum & inoutSumArray, const ArrayIn & inArray) { 00956 #ifdef HAVE_INTREPID_DEBUG 00957 TEUCHOS_TEST_FOR_EXCEPTION( ( inArray.rank() != inoutSumArray.rank() ), 00958 std::invalid_argument, 00959 ">>> ERROR (RealSpaceTools::add): Array arguments must have identical ranks!"); 00960 for (int i=0; i<inArray.rank(); i++) { 00961 TEUCHOS_TEST_FOR_EXCEPTION( ( inArray.dimension(i) != inoutSumArray.dimension(i) ), 00962 std::invalid_argument, 00963 ">>> ERROR (RealSpaceTools::add): Dimensions of array arguments do not agree!"); 00964 } 00965 #endif 00966 00967 for (int i=0; i<inArray.size(); i++) { 00968 inoutSumArray[i] += inArray[i]; 00969 } 00970 } 00971 00972 00973 00974 template<class Scalar> 00975 void RealSpaceTools<Scalar>::subtract(Scalar* diffArray, const Scalar* inArray1, const Scalar* inArray2, const int size) { 00976 for (int i=0; i<size; i++) { 00977 diffArray[i] = inArray1[i] - inArray2[i]; 00978 } 00979 } 00980 00981 00982 00983 template<class Scalar> 00984 void RealSpaceTools<Scalar>::subtract(Scalar* inoutDiffArray, const Scalar* inArray, const int size) { 00985 for (int i=0; i<size; i++) { 00986 inoutDiffArray[i] -= inArray[i]; 00987 } 00988 } 00989 00990 00991 00992 template<class Scalar> 00993 template<class ArrayDiff, class ArrayIn1, class ArrayIn2> 00994 void RealSpaceTools<Scalar>::subtract(ArrayDiff & diffArray, const ArrayIn1 & inArray1, const ArrayIn2 & inArray2) { 00995 #ifdef HAVE_INTREPID_DEBUG 00996 TEUCHOS_TEST_FOR_EXCEPTION( ( (inArray1.rank() != inArray2.rank()) || (inArray1.rank() != diffArray.rank()) ), 00997 std::invalid_argument, 00998 ">>> ERROR (RealSpaceTools::subtract): Array arguments must have identical ranks!"); 00999 for (int i=0; i<inArray1.rank(); i++) { 01000 TEUCHOS_TEST_FOR_EXCEPTION( ( (inArray1.dimension(i) != inArray2.dimension(i)) || (inArray1.dimension(i) != diffArray.dimension(i)) ), 01001 std::invalid_argument, 01002 ">>> ERROR (RealSpaceTools::subtract): Dimensions of array arguments do not agree!"); 01003 } 01004 #endif 01005 01006 for (int i=0; i<inArray1.size(); i++) { 01007 diffArray[i] = inArray1[i] - inArray2[i]; 01008 } 01009 } 01010 01011 01012 01013 template<class Scalar> 01014 template<class ArrayDiff, class ArrayIn> 01015 void RealSpaceTools<Scalar>::subtract(ArrayDiff & inoutDiffArray, const ArrayIn & inArray) { 01016 #ifdef HAVE_INTREPID_DEBUG 01017 TEUCHOS_TEST_FOR_EXCEPTION( ( inArray.rank() != inoutDiffArray.rank() ), 01018 std::invalid_argument, 01019 ">>> ERROR (RealSpaceTools::subtract): Array arguments must have identical ranks!"); 01020 for (int i=0; i<inArray.rank(); i++) { 01021 TEUCHOS_TEST_FOR_EXCEPTION( ( inArray.dimension(i) != inoutDiffArray.dimension(i) ), 01022 std::invalid_argument, 01023 ">>> ERROR (RealSpaceTools::subtract): Dimensions of array arguments do not agree!"); 01024 } 01025 #endif 01026 01027 for (int i=0; i<inArray.size(); i++) { 01028 inoutDiffArray[i] -= inArray[i]; 01029 } 01030 } 01031 01032 01033 01034 01035 template<class Scalar> 01036 void RealSpaceTools<Scalar>::scale(Scalar* scaledArray, const Scalar* inArray, const int size, const Scalar scalar) { 01037 for (int i=0; i<size; i++) { 01038 scaledArray[i] = scalar*inArray[i]; 01039 } 01040 } 01041 01042 01043 01044 template<class Scalar> 01045 void RealSpaceTools<Scalar>::scale(Scalar* inoutScaledArray, const int size, const Scalar scalar) { 01046 for (int i=0; i<size; i++) { 01047 inoutScaledArray[i] *= scalar; 01048 } 01049 } 01050 01051 01052 01053 template<class Scalar> 01054 template<class ArrayScaled, class ArrayIn> 01055 void RealSpaceTools<Scalar>::scale(ArrayScaled & scaledArray, const ArrayIn & inArray, const Scalar scalar) { 01056 #ifdef HAVE_INTREPID_DEBUG 01057 TEUCHOS_TEST_FOR_EXCEPTION( ( inArray.rank() != scaledArray.rank() ), 01058 std::invalid_argument, 01059 ">>> ERROR (RealSpaceTools::scale): Array arguments must have identical ranks!"); 01060 for (int i=0; i<inArray.rank(); i++) { 01061 TEUCHOS_TEST_FOR_EXCEPTION( ( inArray.dimension(i) != scaledArray.dimension(i) ), 01062 std::invalid_argument, 01063 ">>> ERROR (RealSpaceTools::scale): Dimensions of array arguments do not agree!"); 01064 } 01065 #endif 01066 01067 for (int i=0; i<inArray.size(); i++) { 01068 scaledArray[i] = scalar*inArray[i]; 01069 } 01070 } 01071 01072 01073 01074 template<class Scalar> 01075 template<class ArrayScaled> 01076 void RealSpaceTools<Scalar>::scale(ArrayScaled & inoutScaledArray, const Scalar scalar) { 01077 for (int i=0; i<inoutScaledArray.size(); i++) { 01078 inoutScaledArray[i] *= scalar; 01079 } 01080 } 01081 01082 01083 01084 01085 template<class Scalar> 01086 Scalar RealSpaceTools<Scalar>::dot(const Scalar* inArray1, const Scalar* inArray2, const int size) { 01087 Scalar dot(0); 01088 for (int i=0; i<size; i++) { 01089 dot += inArray1[i]*inArray2[i]; 01090 } 01091 return dot; 01092 } 01093 01094 01095 01096 template<class Scalar> 01097 template<class ArrayVec1, class ArrayVec2> 01098 Scalar RealSpaceTools<Scalar>::dot(const ArrayVec1 & inVec1, const ArrayVec2 & inVec2) { 01099 #ifdef HAVE_INTREPID_DEBUG 01100 TEUCHOS_TEST_FOR_EXCEPTION( ( (inVec1.rank() != 1) || (inVec2.rank() != 1) ), 01101 std::invalid_argument, 01102 ">>> ERROR (RealSpaceTools::dot): Vector arguments must have rank 1!"); 01103 TEUCHOS_TEST_FOR_EXCEPTION( ( inVec1.dimension(0) != inVec2.dimension(0) ), 01104 std::invalid_argument, 01105 ">>> ERROR (RealSpaceTools::dot): Dimensions of vector arguments must agree!"); 01106 #endif 01107 01108 Scalar dot(0); 01109 for (int i=0; i<inVec1.size(); i++) { 01110 dot += inVec1[i]*inVec2[i]; 01111 } 01112 return dot; 01113 01114 } 01115 01116 01117 01118 template<class Scalar> 01119 template<class ArrayDot, class ArrayVec1, class ArrayVec2> 01120 void RealSpaceTools<Scalar>::dot(ArrayDot & dotArray, const ArrayVec1 & inVecs1, const ArrayVec2 & inVecs2) { 01121 01122 int arrayRank = inVecs1.rank(); 01123 01124 #ifdef HAVE_INTREPID_DEBUG 01125 TEUCHOS_TEST_FOR_EXCEPTION( ( arrayRank != dotArray.rank()+1 ), 01126 std::invalid_argument, 01127 ">>> ERROR (RealSpaceTools::dot): Ranks of norm and vector array arguments are incompatible!"); 01128 TEUCHOS_TEST_FOR_EXCEPTION( ( arrayRank != inVecs2.rank() ), 01129 std::invalid_argument, 01130 ">>> ERROR (RealSpaceTools::dot): Ranks of input vector arguments must be identical!"); 01131 TEUCHOS_TEST_FOR_EXCEPTION( ( (arrayRank < 2) || (arrayRank > 3) ), 01132 std::invalid_argument, 01133 ">>> ERROR (RealSpaceTools::dot): Rank of input vector arguments must be 2 or 3!"); 01134 for (int i=0; i<arrayRank; i++) { 01135 TEUCHOS_TEST_FOR_EXCEPTION( ( inVecs1.dimension(i) != inVecs2.dimension(i) ), 01136 std::invalid_argument, 01137 ">>> ERROR (RealSpaceTools::dot): Dimensions of input vector arguments do not agree!"); 01138 } 01139 for (int i=0; i<arrayRank-1; i++) { 01140 TEUCHOS_TEST_FOR_EXCEPTION( ( inVecs1.dimension(i) != dotArray.dimension(i) ), 01141 std::invalid_argument, 01142 ">>> ERROR (RealSpaceTools::dot): Dimensions of dot-product and vector arrays do not agree!"); 01143 } 01144 #endif 01145 01146 int dim_i0 = 1; // first index dimension (e.g. cell index) 01147 int dim_i1 = 1; // second index dimension (e.g. point index) 01148 int dim = inVecs1.dimension(arrayRank-1); // spatial dimension 01149 01150 // determine i0 and i1 dimensions 01151 switch(arrayRank) { 01152 case 3: 01153 dim_i0 = inVecs1.dimension(0); 01154 dim_i1 = inVecs1.dimension(1); 01155 break; 01156 case 2: 01157 dim_i1 = inVecs1.dimension(0); 01158 break; 01159 } 01160 01161 int offset_i0, offset, dotOffset; 01162 for (int i0=0; i0<dim_i0; i0++) { 01163 offset_i0 = i0*dim_i1; 01164 for (int i1=0; i1<dim_i1; i1++) { 01165 offset = offset_i0 + i1; 01166 dotOffset = offset; 01167 offset *= dim; 01168 Scalar dot(0); 01169 for (int i=0; i<dim; i++) { 01170 dot += inVecs1[offset+i]*inVecs2[offset+i]; 01171 } 01172 dotArray[dotOffset] = dot; 01173 } 01174 } 01175 } 01176 01177 01178 01179 template<class Scalar> 01180 void RealSpaceTools<Scalar>::matvec(Scalar* matVec, const Scalar* inMat, const Scalar* inVec, const int dim) { 01181 for (int i=0; i<dim; i++) { 01182 Scalar sumdot(0); 01183 for (int j=0; j<dim; j++) { 01184 sumdot += inMat[i*dim+j]*inVec[j]; 01185 } 01186 matVec[i] = sumdot; 01187 } 01188 } 01189 01190 01191 01192 template<class Scalar> 01193 template<class ArrayMatVec, class ArrayMat, class ArrayVec> 01194 void RealSpaceTools<Scalar>::matvec(ArrayMatVec & matVecs, const ArrayMat & inMats, const ArrayVec & inVecs) { 01195 int matArrayRank = inMats.rank(); 01196 01197 #ifdef HAVE_INTREPID_DEBUG 01198 TEUCHOS_TEST_FOR_EXCEPTION( ( matArrayRank != inVecs.rank()+1 ), 01199 std::invalid_argument, 01200 ">>> ERROR (RealSpaceTools::matvec): Vector and matrix array arguments do not have compatible ranks!"); 01201 TEUCHOS_TEST_FOR_EXCEPTION( ( (matArrayRank < 3) || (matArrayRank > 4) ), 01202 std::invalid_argument, 01203 ">>> ERROR (RealSpaceTools::matvec): Rank of matrix array must be 3 or 4!"); 01204 TEUCHOS_TEST_FOR_EXCEPTION( ( matVecs.rank() != inVecs.rank() ), 01205 std::invalid_argument, 01206 ">>> ERROR (RealSpaceTools::matvec): Vector arrays must be have the same rank!"); 01207 for (int i=0; i<matArrayRank-1; i++) { 01208 TEUCHOS_TEST_FOR_EXCEPTION( ( inMats.dimension(i) != inVecs.dimension(i) ), 01209 std::invalid_argument, 01210 ">>> ERROR (RealSpaceTools::matvec): Dimensions of vector and matrix array arguments do not agree!"); 01211 } 01212 for (int i=0; i<inVecs.rank(); i++) { 01213 TEUCHOS_TEST_FOR_EXCEPTION( ( matVecs.dimension(i) != inVecs.dimension(i) ), 01214 std::invalid_argument, 01215 ">>> ERROR (RealSpaceTools::matvec): Dimensions of vector array arguments do not agree!"); 01216 } 01217 TEUCHOS_TEST_FOR_EXCEPTION( ( inMats.dimension(matArrayRank-2) != inMats.dimension(matArrayRank-1) ), 01218 std::invalid_argument, 01219 ">>> ERROR (RealSpaceTools::matvec): Matrices are not square!"); 01220 #endif 01221 01222 int dim_i0 = 1; // first index dimension (e.g. cell index) 01223 int dim_i1 = 1; // second index dimension (e.g. point index) 01224 int dim = inMats.dimension(matArrayRank-2); // spatial dimension 01225 01226 // determine i0 and i1 dimensions 01227 switch(matArrayRank) { 01228 case 4: 01229 dim_i0 = inMats.dimension(0); 01230 dim_i1 = inMats.dimension(1); 01231 break; 01232 case 3: 01233 dim_i1 = inMats.dimension(0); 01234 break; 01235 } 01236 01237 int offset_i0, offset, vecOffset; 01238 01239 for (int i0=0; i0<dim_i0; i0++) { 01240 offset_i0 = i0*dim_i1; 01241 for (int i1=0; i1<dim_i1; i1++) { 01242 offset = offset_i0 + i1; 01243 vecOffset = offset*dim; 01244 offset = vecOffset*dim; 01245 01246 for (int i=0; i<dim; i++) { 01247 Scalar sumdot(0); 01248 for (int j=0; j<dim; j++) { 01249 sumdot += inMats[offset+i*dim+j]*inVecs[vecOffset+j]; 01250 } 01251 matVecs[vecOffset+i] = sumdot; 01252 } 01253 } 01254 } 01255 } 01256 01257 01258 template<class Scalar> 01259 template<class ArrayVecProd, class ArrayIn1, class ArrayIn2> 01260 void RealSpaceTools<Scalar>::vecprod(ArrayVecProd & vecProd, const ArrayIn1 & inLeft, const ArrayIn2 & inRight) { 01261 01262 #ifdef HAVE_INTREPID_DEBUG 01263 /* 01264 * Check array rank and spatial dimension range (if applicable) 01265 * (1) all array arguments are required to have matching dimensions and rank: (D), (I0,D) or (I0,I1,D) 01266 * (2) spatial dimension should be 2 or 3 01267 */ 01268 std::string errmsg = ">>> ERROR (RealSpaceTools::vecprod):"; 01269 01270 // (1) check rank range on inLeft and then compare the other arrays with inLeft 01271 TEUCHOS_TEST_FOR_EXCEPTION( !requireRankRange(errmsg, inLeft, 1,3), std::invalid_argument, errmsg); 01272 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, inLeft, inRight), std::invalid_argument, errmsg); 01273 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, inLeft, vecProd), std::invalid_argument, errmsg); 01274 01275 // (2) spatial dimension ordinal = array rank - 1. Suffices to check only one array because we just 01276 // checked whether or not the arrays have matching dimensions. 01277 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, inLeft, inLeft.rank() - 1, 2,3), std::invalid_argument, errmsg); 01278 01279 #endif 01280 01281 int spaceDim = inLeft.dimension(inLeft.rank() - 1); 01282 01283 switch(inLeft.rank() ){ 01284 01285 case 1: 01286 { 01287 vecProd(0) = inLeft(1)*inRight(2) - inLeft(2)*inRight(1); 01288 vecProd(1) = inLeft(2)*inRight(0) - inLeft(0)*inRight(2); 01289 vecProd(2) = inLeft(0)*inRight(1) - inLeft(1)*inRight(0); 01290 } 01291 break; 01292 01293 case 2: 01294 { 01295 int dim0 = inLeft.dimension(0); 01296 if(spaceDim == 3) { 01297 for(int i0 = 0; i0 < dim0; i0++){ 01298 vecProd(i0, 0) = inLeft(i0, 1)*inRight(i0, 2) - inLeft(i0, 2)*inRight(i0, 1); 01299 vecProd(i0, 1) = inLeft(i0, 2)*inRight(i0, 0) - inLeft(i0, 0)*inRight(i0, 2); 01300 vecProd(i0, 2) = inLeft(i0, 0)*inRight(i0, 1) - inLeft(i0, 1)*inRight(i0, 0); 01301 }// i0 01302 } //spaceDim == 3 01303 else if(spaceDim == 2){ 01304 for(int i0 = 0; i0 < dim0; i0++){ 01305 // vecprod is scalar - do we still want result to be (i0,i1,D)? 01306 vecProd(i0, 0) = inLeft(i0, 0)*inRight(i0, 1) - inLeft(i0, 1)*inRight(i0, 0); 01307 }// i0 01308 }// spaceDim == 2 01309 }// case 2 01310 break; 01311 01312 case 3: 01313 { 01314 int dim0 = inLeft.dimension(0); 01315 int dim1 = inLeft.dimension(1); 01316 if(spaceDim == 3) { 01317 for(int i0 = 0; i0 < dim0; i0++){ 01318 for(int i1 = 0; i1 < dim1; i1++){ 01319 vecProd(i0, i1, 0) = inLeft(i0, i1, 1)*inRight(i0, i1, 2) - inLeft(i0, i1, 2)*inRight(i0, i1, 1); 01320 vecProd(i0, i1, 1) = inLeft(i0, i1, 2)*inRight(i0, i1, 0) - inLeft(i0, i1, 0)*inRight(i0, i1, 2); 01321 vecProd(i0, i1, 2) = inLeft(i0, i1, 0)*inRight(i0, i1, 1) - inLeft(i0, i1, 1)*inRight(i0, i1, 0); 01322 }// i1 01323 }// i0 01324 } //spaceDim == 3 01325 else if(spaceDim == 2){ 01326 for(int i0 = 0; i0 < dim0; i0++){ 01327 for(int i1 = 0; i1 < dim1; i1++){ 01328 // vecprod is scalar - do we still want result to be (i0,i1,D)? 01329 vecProd(i0, i1, 0) = inLeft(i0, i1, 0)*inRight(i0, i1, 1) - inLeft(i0, i1, 1)*inRight(i0, i1, 0); 01330 }// i1 01331 }// i0 01332 }// spaceDim == 2 01333 } // case 3 01334 break; 01335 01336 default: 01337 TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, 01338 ">>> ERROR (RealSpaceTools::vecprod): rank-1,2,3 arrays required"); 01339 } 01340 01341 } 01342 01343 } // namespace Intrepid
1.7.6.1