75#ifdef HAVE_INTREPID_DEBUG
76 TEUCHOS_TEST_FOR_EXCEPTION( ( getrank(inArray) != getrank(absArray) ),
77 std::invalid_argument,
78 ">>> ERROR (RealSpaceTools::absval): Array arguments must have identical ranks!");
79 for (
size_t i=0; i<getrank(inArray); i++) {
80 TEUCHOS_TEST_FOR_EXCEPTION( (
static_cast<size_t>(inArray.dimension(i)) !=
static_cast<size_t>(absArray.dimension(i)) ),
81 std::invalid_argument,
82 ">>> ERROR (RealSpaceTools::absval): Dimensions of array arguments do not agree!");
86 ArrayWrapper<Scalar,ArrayAbs, Rank<ArrayAbs >::value,
false>absArrayWrap(absArray);
87 ArrayWrapper<Scalar,ArrayIn, Rank<ArrayIn >::value,
true>inArrayWrap(inArray);
89 int inArrayRank=getrank(inArray);
93 for (
size_t i=0; i<static_cast<size_t>(
static_cast<size_t>(inArray.dimension(0))); i++)
94 for (
size_t j=0; j<static_cast<size_t>(
static_cast<size_t>(inArray.dimension(1))); j++)
95 for (
size_t k=0; k<static_cast<size_t>(
static_cast<size_t>(inArray.dimension(2))); k++)
96 for (
size_t l=0; l<static_cast<size_t>(
static_cast<size_t>(inArray.dimension(3))); l++)
97 for (
size_t m=0; m<static_cast<size_t>(
static_cast<size_t>(inArray.dimension(4))); m++){
98 absArrayWrap(i,j,k,l,m) = std::abs(inArrayWrap(i,j,k,l,m));
100 }
else if(inArrayRank==4){
101 for (
size_t i=0; i<static_cast<size_t>(
static_cast<size_t>(inArray.dimension(0))); i++)
102 for (
size_t j=0; j<static_cast<size_t>(
static_cast<size_t>(inArray.dimension(1))); j++)
103 for (
size_t k=0; k<static_cast<size_t>(
static_cast<size_t>(inArray.dimension(2))); k++)
104 for (
size_t l=0; l<static_cast<size_t>(
static_cast<size_t>(inArray.dimension(3))); l++){
105 absArrayWrap(i,j,k,l) = std::abs(inArrayWrap(i,j,k,l));
107 }
else if(inArrayRank==3){
108 for (
size_t i=0; i<static_cast<size_t>(
static_cast<size_t>(inArray.dimension(0))); i++)
109 for (
size_t j=0; j<static_cast<size_t>(
static_cast<size_t>(inArray.dimension(1))); j++)
110 for (
size_t k=0; k<static_cast<size_t>(
static_cast<size_t>(inArray.dimension(2))); k++){
111 absArrayWrap(i,j,k) = std::abs(inArrayWrap(i,j,k));
113 }
else if(inArrayRank==2){
114 for (
size_t i=0; i<static_cast<size_t>(
static_cast<size_t>(inArray.dimension(0))); i++)
115 for (
size_t j=0; j<static_cast<size_t>(
static_cast<size_t>(inArray.dimension(1))); j++){
116 absArrayWrap(i,j) = std::abs(inArrayWrap(i,j));
118 }
else if(inArrayRank==1){
119 for (
size_t i=0; i<static_cast<size_t>(
static_cast<size_t>(inArray.dimension(0))); i++){
120 absArrayWrap(i) = std::abs(inArrayWrap(i));
171#ifdef HAVE_INTREPID_DEBUG
172 TEUCHOS_TEST_FOR_EXCEPTION( ( !(getrank(inVec) >= 1 && getrank(inVec) <= 5) ),
173 std::invalid_argument,
174 ">>> ERROR (RealSpaceTools::vectorNorm): Vector argument must have rank 1!");
176 ArrayWrapper<Scalar,ArrayIn, Rank<ArrayIn >::value,
true>inVecWrap(inVec);
177 int inVecRank=getrank(inVecWrap);
178 Scalar temp = (Scalar)0;
182 for (
size_t i=0; i<static_cast<size_t>(inVec.dimension(0)); i++)
183 for (
size_t j=0; j<static_cast<size_t>(inVec.dimension(1)); j++)
184 for (
size_t k=0; k<static_cast<size_t>(inVec.dimension(2)); k++)
185 for (
size_t l=0; l<static_cast<size_t>(inVec.dimension(3)); l++)
186 for (
size_t m=0; m<static_cast<size_t>(inVec.dimension(4)); m++)
187 temp += inVecWrap(i,j,k,l,m)*inVecWrap(i,j,k,l,m);
188 }
else if(inVecRank==4){
189 for (
size_t i=0; i<static_cast<size_t>(inVec.dimension(0)); i++)
190 for (
size_t j=0; j<static_cast<size_t>(inVec.dimension(1)); j++)
191 for (
size_t k=0; k<static_cast<size_t>(inVec.dimension(2)); k++)
192 for (
size_t l=0; l<static_cast<size_t>(inVec.dimension(3)); l++)
193 temp += inVecWrap(i,j,k,l)*inVecWrap(i,j,k,l);
194 }
else if(inVecRank==3){
195 for (
size_t i=0; i<static_cast<size_t>(inVec.dimension(0)); i++)
196 for (
size_t j=0; j<static_cast<size_t>(inVec.dimension(1)); j++)
197 for (
size_t k=0; k<static_cast<size_t>(inVec.dimension(2)); k++)
198 temp += inVecWrap(i,j,k)*inVecWrap(i,j,k);
199 }
else if(inVecRank==2){
200 for (
size_t i=0; i<static_cast<size_t>(inVec.dimension(0)); i++)
201 for (
size_t j=0; j<static_cast<size_t>(inVec.dimension(1)); j++)
202 temp += inVecWrap(i,j)*inVecWrap(i,j);
203 }
else if(inVecRank==1){
204 for (
size_t i=0; i<static_cast<size_t>(inVec.dimension(0)); i++)
205 temp += inVecWrap(i)*inVecWrap(i);
207 temp = std::sqrt(temp);
213 temp = std::abs(inVecWrap(0,0,0,0,0));
214 for (
size_t i=0; i<static_cast<size_t>(inVec.dimension(0)); i++)
215 for (
size_t j=0; j<static_cast<size_t>(inVec.dimension(1)); j++)
216 for (
size_t k=0; k<static_cast<size_t>(inVec.dimension(2)); k++)
217 for (
size_t l=0; l<static_cast<size_t>(inVec.dimension(3)); l++)
218 for (
size_t m=1; m<static_cast<size_t>(inVec.dimension(4)); m++){
219 Scalar absData = std::abs(inVecWrap(i,j,k,l,m));
220 if (temp < absData) temp = absData;
222 }
else if(inVecRank==4){
223 temp = std::abs(inVecWrap(0,0,0,0));
224 for (
size_t i=0; i<static_cast<size_t>(inVec.dimension(0)); i++)
225 for (
size_t j=0; j<static_cast<size_t>(inVec.dimension(1)); j++)
226 for (
size_t k=0; k<static_cast<size_t>(inVec.dimension(2)); k++)
227 for (
size_t l=1; l<static_cast<size_t>(inVec.dimension(3)); l++){
228 Scalar absData = std::abs(inVecWrap(i,j,k,l));
229 if (temp < absData) temp = absData;
231 }
else if(inVecRank==3){
232 temp = std::abs(inVecWrap(0,0,0));
233 for (
size_t i=0; i<static_cast<size_t>(inVec.dimension(0)); i++)
234 for (
size_t j=0; j<static_cast<size_t>(inVec.dimension(1)); j++)
235 for (
size_t k=0; k<static_cast<size_t>(inVec.dimension(2)); k++){
236 Scalar absData = std::abs(inVecWrap(i,j,k));
237 if (temp < absData) temp = absData;
239 }
else if(inVecRank==2){
240 temp = std::abs(inVecWrap(0,0));
241 for (
size_t i=0; i<static_cast<size_t>(inVec.dimension(0)); i++)
242 for (
size_t j=0; j<static_cast<size_t>(inVec.dimension(1)); j++){
243 Scalar absData = std::abs(inVecWrap(i,j));
244 if (temp < absData) temp = absData;
246 }
else if(inVecRank==1){
247 temp = std::abs(inVecWrap(0));
248 for (
size_t i=1; i<static_cast<size_t>(inVec.dimension(0)); i++){
249 Scalar absData = std::abs(inVecWrap(i));
250 if (temp < absData) temp = absData;
257 for (
size_t i=0; i<static_cast<size_t>(inVec.dimension(0)); i++)
258 for (
size_t j=0; j<static_cast<size_t>(inVec.dimension(1)); j++)
259 for (
size_t k=0; k<static_cast<size_t>(inVec.dimension(2)); k++)
260 for (
size_t l=0; l<static_cast<size_t>(inVec.dimension(3)); l++)
261 for (
size_t m=0; m<static_cast<size_t>(inVec.dimension(4)); m++){
262 temp += std::abs(inVecWrap(i,j,k,l,m));
264 }
else if(inVecRank==4){
265 for (
size_t i=0; i<static_cast<size_t>(inVec.dimension(0)); i++)
266 for (
size_t j=0; j<static_cast<size_t>(inVec.dimension(1)); j++)
267 for (
size_t k=0; k<static_cast<size_t>(inVec.dimension(2)); k++)
268 for (
size_t l=0; l<static_cast<size_t>(inVec.dimension(3)); l++){
269 temp += std::abs(inVecWrap(i,j,k,l));
271 }
else if(inVecRank==3){
272 for (
size_t i=0; i<static_cast<size_t>(inVec.dimension(0)); i++)
273 for (
size_t j=0; j<static_cast<size_t>(inVec.dimension(1)); j++)
274 for (
size_t k=0; k<static_cast<size_t>(inVec.dimension(2)); k++){
275 temp += std::abs(inVecWrap(i,j,k));
277 }
else if(inVecRank==2){
278 for (
size_t i=0; i<static_cast<size_t>(inVec.dimension(0)); i++)
279 for (
size_t j=0; j<static_cast<size_t>(inVec.dimension(1)); j++){
280 temp += std::abs(inVecWrap(i,j));
282 }
else if(inVecRank==1){
283 for (
size_t i=0; i<static_cast<size_t>(inVec.dimension(0)); i++){
284 temp += std::abs(inVecWrap(i));
290 TEUCHOS_TEST_FOR_EXCEPTION( ( (normType != NORM_TWO) && (normType != NORM_INF) && (normType != NORM_ONE) ),
291 std::invalid_argument,
292 ">>> ERROR (RealSpaceTools::vectorNorm): Invalid argument normType.");
340 ArrayWrapper<Scalar,ArrayNorm, Rank<ArrayNorm >::value,
false>normArrayWrap(normArray);
341 ArrayWrapper<Scalar,ArrayIn, Rank<ArrayIn >::value,
true>inVecsWrap(inVecs);
343 size_t arrayRank = getrank(inVecs);
344#ifdef HAVE_INTREPID_DEBUG
345 TEUCHOS_TEST_FOR_EXCEPTION( ( arrayRank != getrank(normArray)+1 ),
346 std::invalid_argument,
347 ">>> ERROR (RealSpaceTools::vectorNorm): Ranks of norm and vector array arguments are incompatible!");
348 TEUCHOS_TEST_FOR_EXCEPTION( ( (arrayRank < 2) || (arrayRank > 3) ),
349 std::invalid_argument,
350 ">>> ERROR (RealSpaceTools::vectorNorm): Rank of vector array must be 2 or 3!");
351 for (
size_t i=0; i<arrayRank-1; i++) {
352 TEUCHOS_TEST_FOR_EXCEPTION( ( inVecs.dimension(i) != normArray.dimension(i) ),
353 std::invalid_argument,
354 ">>> ERROR (RealSpaceTools::vectorNorm): Dimensions of norm and vector arguments do not agree!");
360 size_t dim =
static_cast<size_t>(inVecs.dimension(arrayRank-1));
365 dim_i0 =
static_cast<size_t>(inVecs.dimension(0));
366 dim_i1 =
static_cast<size_t>(inVecs.dimension(1));
369 for (
size_t i0=0; i0<dim_i0; i0++) {
370 for (
size_t i1=0; i1<dim_i1; i1++) {
371 Scalar temp = (Scalar)0;
372 for(
size_t i = 0; i < dim; i++){
373 temp += inVecsWrap(i0,i1,i)*inVecsWrap(i0,i1,i);
375 normArrayWrap(i0,i1) = std::sqrt(temp);
382 for (
size_t i0=0; i0<dim_i0; i0++) {
383 for (
size_t i1=0; i1<dim_i1; i1++) {
384 Scalar temp = (Scalar)0;
385 temp = std::abs(inVecsWrap(i0,i1,0));
386 for(
size_t i = 1; i < dim; i++){
387 Scalar absData = std::abs(inVecsWrap(i0,i1,i));
388 if (temp < absData) temp = absData;
390 normArrayWrap(i0,i1) = temp;
397 for (
size_t i0=0; i0<dim_i0; i0++) {
398 for (
size_t i1=0; i1<dim_i1; i1++) {
399 Scalar temp = (Scalar)0;
400 for(
size_t i = 0; i < dim; i++){
401 temp += std::abs(inVecsWrap(i0,i1,i));
403 normArrayWrap(i0,i1) = temp;
410 TEUCHOS_TEST_FOR_EXCEPTION( ( (normType != NORM_TWO) && (normType != NORM_INF) && (normType != NORM_ONE) ),
411 std::invalid_argument,
412 ">>> ERROR (RealSpaceTools::vectorNorm): Invalid argument normType.");
419 dim_i1 =
static_cast<size_t>(inVecs.dimension(0));
423 for (
size_t i1=0; i1<dim_i1; i1++) {
424 Scalar temp = (Scalar)0;
425 for(
size_t i = 0; i < dim; i++){
426 temp += inVecsWrap(i1,i)*inVecsWrap(i1,i);
428 normArrayWrap(i1) = std::sqrt(temp);
435 for (
size_t i1=0; i1<dim_i1; i1++) {
436 Scalar temp = (Scalar)0;
437 temp = std::abs(inVecsWrap(i1,0));
438 for(
size_t i = 1; i < dim; i++){
439 Scalar absData = std::abs(inVecsWrap(i1,i));
440 if (temp < absData) temp = absData;
442 normArrayWrap(i1) = temp;
448 for (
size_t i1=0; i1<dim_i1; i1++) {
449 Scalar temp = (Scalar)0;
450 for(
size_t i = 0; i < dim; i++){
451 temp += std::abs(inVecsWrap(i1,i));
453 normArrayWrap(i1) = temp;
459 TEUCHOS_TEST_FOR_EXCEPTION( ( (normType != NORM_TWO) && (normType != NORM_INF) && (normType != NORM_ONE) ),
460 std::invalid_argument,
461 ">>> ERROR (RealSpaceTools::vectorNorm): Invalid argument normType.");
591 size_t arrayRank = getrank(inMats);
592#ifdef HAVE_INTREPID_DEBUG
593 TEUCHOS_TEST_FOR_EXCEPTION( ( arrayRank != getrank(transposeMats) ),
594 std::invalid_argument,
595 ">>> ERROR (RealSpaceTools::transpose): Matrix array arguments do not have identical ranks!");
596 TEUCHOS_TEST_FOR_EXCEPTION( ( (arrayRank < 2) || (arrayRank > 4) ),
597 std::invalid_argument,
598 ">>> ERROR (RealSpaceTools::transpose): Rank of matrix array must be 2, 3, or 4!");
599 for (
size_t i=0; i<arrayRank; i++) {
600 TEUCHOS_TEST_FOR_EXCEPTION( (
static_cast<size_t>(inMats.dimension(i)) !=
static_cast<size_t>(transposeMats.dimension(i)) ),
601 std::invalid_argument,
602 ">>> ERROR (RealSpaceTools::transpose): Dimensions of matrix arguments do not agree!");
604 TEUCHOS_TEST_FOR_EXCEPTION( (
static_cast<size_t>(inMats.dimension(arrayRank-2)) !=
static_cast<size_t>(inMats.dimension(arrayRank-1)) ),
605 std::invalid_argument,
606 ">>> ERROR (RealSpaceTools::transpose): Matrices are not square!");
610 size_t dim =
static_cast<size_t>(inMats.dimension(arrayRank-2));
613 ArrayWrapper<Scalar,ArrayTranspose,Rank<ArrayTranspose>::value,
false>transposeArr(transposeMats);
614 ArrayWrapper<Scalar,ArrayIn,Rank<ArrayIn>::value,
true>inputArr(inMats);
618 dim_i0 =
static_cast<size_t>(inMats.dimension(0));
619 dim_i1 =
static_cast<size_t>(inMats.dimension(1));
621 for (
size_t i0=0; i0<dim_i0; i0++) {
622 for (
size_t i1=0; i1<dim_i1; i1++) {
623 for(
size_t i=0; i < dim; i++){
624 transposeArr(i0,i1,i,i)=inputArr(i0,i1,i,i);
625 for(
size_t j=i+1; j < dim; j++){
626 transposeArr(i0,i1,i,j)=inputArr(i0,i1,j,i);
627 transposeArr(i0,i1,j,i)=inputArr(i0,i1,i,j);
635 dim_i1 =
static_cast<size_t>(inMats.dimension(0));
636 for (
size_t i1=0; i1<dim_i1; i1++) {
637 for(
size_t i=0; i < dim; i++){
638 transposeArr(i1,i,i)=inputArr(i1,i,i);
639 for(
size_t j=i+1; j < dim; j++){
640 transposeArr(i1,i,j)=inputArr(i1,j,i);
641 transposeArr(i1,j,i)=inputArr(i1,i,j);
759 ArrayWrapper<Scalar,ArrayInverse, Rank<ArrayInverse >::value,
false>inverseMatsWrap(inverseMats);
760 ArrayWrapper<Scalar,ArrayIn, Rank<ArrayIn >::value,
true>inMatsWrap(inMats);
762 size_t arrayRank = getrank(inMats);
764#ifdef HAVE_INTREPID_DEBUG
765 TEUCHOS_TEST_FOR_EXCEPTION( ( arrayRank != getrank(inverseMats) ),
766 std::invalid_argument,
767 ">>> ERROR (RealSpaceTools::inverse): Matrix array arguments do not have identical ranks!");
768 TEUCHOS_TEST_FOR_EXCEPTION( ( (arrayRank < 2) || (arrayRank > 4) ),
769 std::invalid_argument,
770 ">>> ERROR (RealSpaceTools::inverse): Rank of matrix array must be 2, 3, or 4!");
771 for (
size_t i=0; i<arrayRank; i++) {
772 TEUCHOS_TEST_FOR_EXCEPTION( (
static_cast<size_t>(inMats.dimension(i)) !=
static_cast<size_t>(inverseMats.dimension(i)) ),
773 std::invalid_argument,
774 ">>> ERROR (RealSpaceTools::inverse): Dimensions of matrix arguments do not agree!");
776 TEUCHOS_TEST_FOR_EXCEPTION( (
static_cast<size_t>(inMats.dimension(arrayRank-2)) !=
static_cast<size_t>(inMats.dimension(arrayRank-1)) ),
777 std::invalid_argument,
778 ">>> ERROR (RealSpaceTools::inverse): Matrices are not square!");
779 TEUCHOS_TEST_FOR_EXCEPTION( ( (
static_cast<size_t>(inMats.dimension(arrayRank-2)) < 1) || (
static_cast<size_t>(inMats.dimension(arrayRank-2)) > 3) ),
780 std::invalid_argument,
781 ">>> ERROR (RealSpaceTools::inverse): Spatial dimension must be 1, 2, or 3!");
786 size_t dim =
static_cast<size_t>(inMats.dimension(arrayRank-2));
791 dim_i0 =
static_cast<size_t>(inMats.dimension(0));
792 dim_i1 =
static_cast<size_t>(inMats.dimension(1));
797 for (
size_t i0=0; i0<dim_i0; i0++) {
799 for (
size_t i1=0; i1<dim_i1; i1++) {
802 size_t i, j, rowID = 0, colID = 0;
803 int rowperm[3]={0,1,2};
804 int colperm[3]={0,1,2};
807 for(i=0; i < 3; ++i){
808 for(j=0; j < 3; ++j){
809 if( std::abs( inMatsWrap(i0,i1,i,j) ) > emax){
810 rowID = i; colID = j; emax = std::abs( inMatsWrap(i0,i1,i,j) );
814#ifdef HAVE_INTREPID_DEBUG
815#ifdef HAVE_INTREPID_DEBUG_INF_CHECK
816 TEUCHOS_TEST_FOR_EXCEPTION( ( emax == (Scalar)0 ),
817 std::invalid_argument,
818 ">>> ERROR (Matrix): Inverse of a zero matrix is undefined!");
829 Scalar B[3][3], S[2][2], Bi[3][3];
830 for(i=0; i < 3; ++i){
831 for(j=0; j < 3; ++j){
832 B[i][j] = inMatsWrap(i0,i1,rowperm[i],colperm[j]);
835 B[1][0] /= B[0][0]; B[2][0] /= B[0][0];
836 for(i=0; i < 2; ++i){
837 for(j=0; j < 2; ++j){
838 S[i][j] = B[i+1][j+1] - B[i+1][0] * B[0][j+1];
841 Scalar detS = S[0][0]*S[1][1]- S[0][1]*S[1][0], Si[2][2];
842#ifdef HAVE_INTREPID_DEBUG
843#ifdef HAVE_INTREPID_DEBUG_INF_CHECK
844 TEUCHOS_TEST_FOR_EXCEPTION( ( detS == (Scalar)0 ),
845 std::invalid_argument,
846 ">>> ERROR (Matrix): Inverse of a singular matrix is undefined!");
850 Si[0][0] = S[1][1]/detS; Si[0][1] = -S[0][1]/detS;
851 Si[1][0] = -S[1][0]/detS; Si[1][1] = S[0][0]/detS;
854 Bi[0][j+1] = -( B[0][1]*Si[0][j] + B[0][2]* Si[1][j])/B[0][0];
856 Bi[i+1][0] = -(Si[i][0]*B[1][0] + Si[i][1]*B[2][0]);
858 Bi[0][0] = ((Scalar)1/B[0][0])-Bi[0][1]*B[1][0]-Bi[0][2]*B[2][0];
863 for(i=0; i < 3; ++i){
864 for(j=0; j < 3; ++j){
865 inverseMatsWrap(i0,i1,i,j) = Bi[colperm[i]][rowperm[j]];
875 for (
size_t i0=0; i0<dim_i0; i0++) {
877 for (
size_t i1=0; i1<dim_i1; i1++) {
880 Scalar determinant = inMatsWrap(i0,i1,0,0)*inMatsWrap(i0,i1,1,1)-inMatsWrap(i0,i1,0,1)*inMatsWrap(i0,i1,1,0);
881#ifdef HAVE_INTREPID_DEBUG
882#ifdef HAVE_INTREPID_DEBUG_INF_CHECK
883 TEUCHOS_TEST_FOR_EXCEPTION( ( (inMatsWrap(i0,i1,0,0)==(Scalar)0) && (inMatsWrap(i0,i1,0,1)==(Scalar)0) &&
884 (inMatsWrap(i0,i1,1,0)==(Scalar)0) && (inMatsWrap(i0,i1,1,1)==(Scalar)0) ),
885 std::invalid_argument,
886 ">>> ERROR (Matrix): Inverse of a zero matrix is undefined!");
887 TEUCHOS_TEST_FOR_EXCEPTION( ( determinant == (Scalar)0 ),
888 std::invalid_argument,
889 ">>> ERROR (Matrix): Inverse of a singular matrix is undefined!");
892 inverseMatsWrap(i0,i1,0,0) = inMatsWrap(i0,i1,1,1) / determinant;
893 inverseMatsWrap(i0,i1,0,1) = - inMatsWrap(i0,i1,0,1) / determinant;
895 inverseMatsWrap(i0,i1,1,0) = - inMatsWrap(i0,i1,1,0) / determinant;
896 inverseMatsWrap(i0,i1,1,1) = inMatsWrap(i0,i1,0,0) / determinant;
903 for (
size_t i0=0; i0<dim_i0; i0++) {
904 for (
size_t i1=0; i1<dim_i1; i1++) {
905#ifdef HAVE_INTREPID_DEBUG
906#ifdef HAVE_INTREPID_DEBUG_INF_CHECK
907 TEUCHOS_TEST_FOR_EXCEPTION( ( inMatsWrap(i0,i1,0,0) == (Scalar)0 ),
908 std::invalid_argument,
909 ">>> ERROR (Matrix): Inverse of a zero matrix is undefined!");
913 inverseMatsWrap(i0,i1,0,0) = (Scalar)1 / inMatsWrap(i0,i1,0,0);
924 dim_i1 =
static_cast<size_t>(inMats.dimension(0));
928 for (
size_t i1=0; i1<dim_i1; i1++) {
930 size_t i, j, rowID = 0, colID = 0;
931 int rowperm[3]={0,1,2};
932 int colperm[3]={0,1,2};
935 for(i=0; i < 3; ++i){
936 for(j=0; j < 3; ++j){
937 if( std::abs( inMatsWrap(i1,i,j) ) > emax){
938 rowID = i; colID = j; emax = std::abs( inMatsWrap(i1,i,j) );
942#ifdef HAVE_INTREPID_DEBUG
943#ifdef HAVE_INTREPID_DEBUG_INF_CHECK
944 TEUCHOS_TEST_FOR_EXCEPTION( ( emax == (Scalar)0 ),
945 std::invalid_argument,
946 ">>> ERROR (Matrix): Inverse of a zero matrix is undefined!");
958 Scalar B[3][3], S[2][2], Bi[3][3];
959 for(i=0; i < 3; ++i){
960 for(j=0; j < 3; ++j){
961 B[i][j] = inMatsWrap(i1,rowperm[i],colperm[j]);
964 B[1][0] /= B[0][0]; B[2][0] /= B[0][0];
965 for(i=0; i < 2; ++i){
966 for(j=0; j < 2; ++j){
967 S[i][j] = B[i+1][j+1] - B[i+1][0] * B[0][j+1];
970 Scalar detS = S[0][0]*S[1][1]- S[0][1]*S[1][0], Si[2][2];
971#ifdef HAVE_INTREPID_DEBUG
972#ifdef HAVE_INTREPID_DEBUG_INF_CHECK
973 TEUCHOS_TEST_FOR_EXCEPTION( ( detS == (Scalar)0 ),
974 std::invalid_argument,
975 ">>> ERROR (Matrix): Inverse of a singular matrix is undefined!");
980 Si[0][0] = S[1][1]/detS; Si[0][1] = -S[0][1]/detS;
981 Si[1][0] = -S[1][0]/detS; Si[1][1] = S[0][0]/detS;
984 Bi[0][j+1] = -( B[0][1]*Si[0][j] + B[0][2]* Si[1][j])/B[0][0];
986 Bi[i+1][0] = -(Si[i][0]*B[1][0] + Si[i][1]*B[2][0]);
988 Bi[0][0] = ((Scalar)1/B[0][0])-Bi[0][1]*B[1][0]-Bi[0][2]*B[2][0];
993 for(i=0; i < 3; ++i){
994 for(j=0; j < 3; ++j){
995 inverseMatsWrap(i1,i,j) = Bi[colperm[i]][rowperm[j]];
1005 for (
size_t i1=0; i1<dim_i1; i1++) {
1008 Scalar determinant = inMatsWrap(i1,0,0)*inMatsWrap(i1,1,1)-inMatsWrap(i1,0,1)*inMatsWrap(i1,1,0);
1009#ifdef HAVE_INTREPID_DEBUG
1010#ifdef HAVE_INTREPID_DEBUG_INF_CHECK
1011 TEUCHOS_TEST_FOR_EXCEPTION( ( (inMatsWrap(i1,0,0)==(Scalar)0) && (inMatsWrap(i1,0,1)==(Scalar)0) &&
1012 (inMatsWrap(i1,1,0)==(Scalar)0) && (inMatsWrap(i1,1,1)==(Scalar)0) ),
1013 std::invalid_argument,
1014 ">>> ERROR (Matrix): Inverse of a zero matrix is undefined!");
1015 TEUCHOS_TEST_FOR_EXCEPTION( ( determinant == (Scalar)0 ),
1016 std::invalid_argument,
1017 ">>> ERROR (Matrix): Inverse of a singular matrix is undefined!");
1021 inverseMatsWrap(i1,0,0) = inMatsWrap(i1,1,1) / determinant;
1022 inverseMatsWrap(i1,0,1) = - inMatsWrap(i1,0,1) / determinant;
1024 inverseMatsWrap(i1,1,0) = - inMatsWrap(i1,1,0) / determinant;
1025 inverseMatsWrap(i1,1,1) = inMatsWrap(i1,0,0) / determinant;
1033 for (
size_t i1=0; i1<dim_i1; i1++) {
1034#ifdef HAVE_INTREPID_DEBUG
1035#ifdef HAVE_INTREPID_DEBUG_INF_CHECK
1036 TEUCHOS_TEST_FOR_EXCEPTION( ( inMatsWrap(i1,0,0) == (Scalar)0 ),
1037 std::invalid_argument,
1038 ">>> ERROR (Matrix): Inverse of a zero matrix is undefined!");
1043 inverseMatsWrap(i1,0,0) = (Scalar)1 / inMatsWrap(i1,0,0);
1056 size_t i, j, rowID = 0, colID = 0;
1057 int rowperm[3]={0,1,2};
1058 int colperm[3]={0,1,2};
1061 for(i=0; i < 3; ++i){
1062 for(j=0; j < 3; ++j){
1063 if( std::abs( inMatsWrap(i,j) ) > emax){
1064 rowID = i; colID = j; emax = std::abs( inMatsWrap(i,j) );
1068#ifdef HAVE_INTREPID_DEBUG
1069#ifdef HAVE_INTREPID_DEBUG_INF_CHECK
1070 TEUCHOS_TEST_FOR_EXCEPTION( ( emax == (Scalar)0 ),
1071 std::invalid_argument,
1072 ">>> ERROR (Matrix): Inverse of a zero matrix is undefined!");
1084 Scalar B[3][3], S[2][2], Bi[3][3];
1085 for(i=0; i < 3; ++i){
1086 for(j=0; j < 3; ++j){
1087 B[i][j] = inMatsWrap(rowperm[i],colperm[j]);
1090 B[1][0] /= B[0][0]; B[2][0] /= B[0][0];
1091 for(i=0; i < 2; ++i){
1092 for(j=0; j < 2; ++j){
1093 S[i][j] = B[i+1][j+1] - B[i+1][0] * B[0][j+1];
1096 Scalar detS = S[0][0]*S[1][1]- S[0][1]*S[1][0], Si[2][2];
1097#ifdef HAVE_INTREPID_DEBUG
1098#ifdef HAVE_INTREPID_DEBUG_INF_CHECK
1099 TEUCHOS_TEST_FOR_EXCEPTION( ( detS == (Scalar)0 ),
1100 std::invalid_argument,
1101 ">>> ERROR (Matrix): Inverse of a singular matrix is undefined!");
1105 Si[0][0] = S[1][1]/detS; Si[0][1] = -S[0][1]/detS;
1106 Si[1][0] = -S[1][0]/detS; Si[1][1] = S[0][0]/detS;
1109 Bi[0][j+1] = -( B[0][1]*Si[0][j] + B[0][2]* Si[1][j])/B[0][0];
1111 Bi[i+1][0] = -(Si[i][0]*B[1][0] + Si[i][1]*B[2][0]);
1113 Bi[0][0] = ((Scalar)1/B[0][0])-Bi[0][1]*B[1][0]-Bi[0][2]*B[2][0];
1114 Bi[1][1] = Si[0][0];
1115 Bi[1][2] = Si[0][1];
1116 Bi[2][1] = Si[1][0];
1117 Bi[2][2] = Si[1][1];
1118 for(i=0; i < 3; ++i){
1119 for(j=0; j < 3; ++j){
1120 inverseMatsWrap(i,j) = Bi[colperm[i]][rowperm[j]];
1128 Scalar determinant = inMatsWrap(0,0)*inMatsWrap(1,1)-inMatsWrap(0,1)*inMatsWrap(1,0);
1129#ifdef HAVE_INTREPID_DEBUG
1130#ifdef HAVE_INTREPID_DEBUG_INF_CHECK
1131 TEUCHOS_TEST_FOR_EXCEPTION( ( (inMatsWrap(0,0)==(Scalar)0) && (inMatsWrap(0,1)==(Scalar)0) &&
1132 (inMatsWrap(1,0)==(Scalar)0) && (inMatsWrap(1,1)==(Scalar)0) ),
1133 std::invalid_argument,
1134 ">>> ERROR (Matrix): Inverse of a zero matrix is undefined!");
1135 TEUCHOS_TEST_FOR_EXCEPTION( ( determinant == (Scalar)0 ),
1136 std::invalid_argument,
1137 ">>> ERROR (Matrix): Inverse of a singular matrix is undefined!");
1140 inverseMatsWrap(0,0) = inMatsWrap(1,1) / determinant;
1141 inverseMatsWrap(0,1) = - inMatsWrap(0,1) / determinant;
1143 inverseMatsWrap(1,0) = - inMatsWrap(1,0) / determinant;
1144 inverseMatsWrap(1,1) = inMatsWrap(0,0) / determinant;
1150#ifdef HAVE_INTREPID_DEBUG
1151#ifdef HAVE_INTREPID_DEBUG_INF_CHECK
1152 TEUCHOS_TEST_FOR_EXCEPTION( ( inMatsWrap(0,0) == (Scalar)0 ),
1153 std::invalid_argument,
1154 ">>> ERROR (Matrix): Inverse of a zero matrix is undefined!");
1157 inverseMatsWrap(0,0) = (Scalar)1 / inMatsWrap(0,0);
1322 ArrayWrapper<Scalar,ArrayDet, Rank<ArrayDet >::value,
false>detArrayWrap(detArray);
1323 ArrayWrapper<Scalar,ArrayIn, Rank<ArrayIn >::value,
true>inMatsWrap(inMats);
1325 size_t matArrayRank = getrank(inMats);
1326#ifdef HAVE_INTREPID_DEBUG
1327 TEUCHOS_TEST_FOR_EXCEPTION( ( matArrayRank != getrank(detArray)+2 ),
1328 std::invalid_argument,
1329 ">>> ERROR (RealSpaceTools::det): Determinant and matrix array arguments do not have compatible ranks!");
1330 TEUCHOS_TEST_FOR_EXCEPTION( ( (matArrayRank < 3) || (matArrayRank > 4) ),
1331 std::invalid_argument,
1332 ">>> ERROR (RealSpaceTools::det): Rank of matrix array must be 3 or 4!");
1333 for (
size_t i=0; i<matArrayRank-2; i++) {
1334 TEUCHOS_TEST_FOR_EXCEPTION( (
static_cast<size_t>(inMats.dimension(i)) !=
static_cast<size_t>(detArray.dimension(i)) ),
1335 std::invalid_argument,
1336 ">>> ERROR (RealSpaceTools::det): Dimensions of determinant and matrix array arguments do not agree!");
1338 TEUCHOS_TEST_FOR_EXCEPTION( (
static_cast<size_t>(inMats.dimension(matArrayRank-2)) !=
static_cast<size_t>(inMats.dimension(matArrayRank-1)) ),
1339 std::invalid_argument,
1340 ">>> ERROR (RealSpaceTools::det): Matrices are not square!");
1341 TEUCHOS_TEST_FOR_EXCEPTION( ( (
static_cast<size_t>(inMats.dimension(matArrayRank-2)) < 1) || (
static_cast<size_t>(inMats.dimension(matArrayRank-2)) > 3) ),
1342 std::invalid_argument,
1343 ">>> ERROR (RealSpaceTools::det): Spatial dimension must be 1, 2, or 3!");
1348 size_t dim = inMats.dimension(matArrayRank-2);
1351 switch(matArrayRank) {
1353 dim_i0 =
static_cast<size_t>(inMats.dimension(0));
1354 dim_i1 =
static_cast<size_t>(inMats.dimension(1));
1359 for (
size_t i0=0; i0<dim_i0; i0++) {
1361 for (
size_t i1=0; i1<dim_i1; i1++) {
1366 int rowperm[3]={0,1,2};
1367 int colperm[3]={0,1,2};
1368 Scalar emax(0), determinant(0);
1370 for(i=0; i < 3; ++i){
1371 for(j=0; j < 3; ++j){
1372 if( std::abs( inMatsWrap(i0,i1,i,j) ) > emax){
1373 rowID = i; colID = j; emax = std::abs( inMatsWrap(i0,i1,i,j) );
1386 Scalar B[3][3], S[2][2];
1387 for(i=0; i < 3; ++i){
1388 for(j=0; j < 3; ++j){
1389 B[i][j] = inMatsWrap(i0,i1,rowperm[i],colperm[j]);
1392 B[1][0] /= B[0][0]; B[2][0] /= B[0][0];
1393 for(i=0; i < 2; ++i){
1394 for(j=0; j < 2; ++j){
1395 S[i][j] = B[i+1][j+1] - B[i+1][0] * B[0][j+1];
1398 determinant = B[0][0] * (S[0][0] * S[1][1] - S[0][1] * S[1][0]);
1399 if( rowID ) determinant = -determinant;
1400 if( colID ) determinant = -determinant;
1402 detArrayWrap(i0,i1)= determinant;
1410 for (
size_t i0=0; i0<dim_i0; i0++) {
1412 for (
size_t i1=0; i1<dim_i1; i1++) {
1414 detArrayWrap(i0,i1) = inMatsWrap(i0,i1,0,0)*inMatsWrap(i0,i1,1,1)-inMatsWrap(i0,i1,0,1)*inMatsWrap(i0,i1,1,0);
1422 for (
size_t i0=0; i0<dim_i0; i0++) {
1424 for (
size_t i1=0; i1<dim_i1; i1++) {
1425 detArrayWrap(i0,i1) = inMatsWrap(i0,i1,0,0);
1434 dim_i1 =
static_cast<size_t>(inMats.dimension(0));
1437 for (
size_t i1=0; i1<dim_i1; i1++) {
1441 int rowperm[3]={0,1,2};
1442 int colperm[3]={0,1,2};
1443 Scalar emax(0), determinant(0);
1446 for(i=0; i < 3; ++i){
1447 for(j=0; j < 3; ++j){
1448 if( std::abs( inMatsWrap(i1,i,j) ) > emax){
1449 rowID = i; colID = j; emax = std::abs( inMatsWrap(i1,i,j) );
1462 Scalar B[3][3], S[2][2];
1463 for(i=0; i < 3; ++i){
1464 for(j=0; j < 3; ++j){
1465 B[i][j] = inMatsWrap(i1,rowperm[i],colperm[j]);
1468 B[1][0] /= B[0][0]; B[2][0] /= B[0][0];
1469 for(i=0; i < 2; ++i){
1470 for(j=0; j < 2; ++j){
1471 S[i][j] = B[i+1][j+1] - B[i+1][0] * B[0][j+1];
1474 determinant = B[0][0] * (S[0][0] * S[1][1] - S[0][1] * S[1][0]);
1475 if( rowID ) determinant = -determinant;
1476 if( colID ) determinant = -determinant;
1478 detArrayWrap(i1) = determinant;
1484 for (
size_t i1=0; i1<dim_i1; i1++) {
1485 detArrayWrap(i1) = inMatsWrap(i1,0,0)*inMatsWrap(i1,1,1)-inMatsWrap(i1,0,1)*inMatsWrap(i1,1,0);
1491 for (
size_t i1=0; i1<dim_i1; i1++) {
1492 detArrayWrap(i1) = inMatsWrap(i1,0,0);
1526#ifdef HAVE_INTREPID_DEBUG
1527 TEUCHOS_TEST_FOR_EXCEPTION( ( (getrank(inArray1) != getrank(inArray2)) || (getrank(inArray1) != getrank(sumArray)) ),
1528 std::invalid_argument,
1529 ">>> ERROR (RealSpaceTools::add): Array arguments must have identical ranks!");
1530 for (
size_t i=0; i<getrank(inArray1); i++) {
1531 TEUCHOS_TEST_FOR_EXCEPTION( ( (
static_cast<size_t>(inArray1.dimension(i)) !=
static_cast<size_t>(inArray2.dimension(i))) || (
static_cast<size_t>(inArray1.dimension(i)) !=
static_cast<size_t>(sumArray.dimension(i))) ),
1532 std::invalid_argument,
1533 ">>> ERROR (RealSpaceTools::add): Dimensions of array arguments do not agree!");
1540 ArrayWrapper<Scalar,ArraySum, Rank<ArraySum >::value,
false>sumArrayWrap(sumArray);
1541 ArrayWrapper<Scalar,ArrayIn1, Rank<ArrayIn1 >::value,
true>inArray1Wrap(inArray1);
1542 ArrayWrapper<Scalar,ArrayIn2, Rank<ArrayIn2 >::value,
true>inArray2Wrap(inArray2);
1543 int inArrayRank=getrank(inArray1);
1547 for (
size_t i=0; i<static_cast<size_t>(inArray1.dimension(0)); i++)
1548 for (
size_t j=0; j<static_cast<size_t>(inArray1.dimension(1)); j++)
1549 for (
size_t k=0; k<static_cast<size_t>(inArray1.dimension(2)); k++)
1550 for (
size_t l=0; l<static_cast<size_t>(inArray1.dimension(3)); l++)
1551 for (
size_t m=0; m<static_cast<size_t>(inArray1.dimension(4)); m++){
1552 sumArrayWrap(i,j,k,l,m) = inArray1Wrap(i,j,k,l,m)+inArray2Wrap(i,j,k,l,m);
1554 }
else if(inArrayRank==4){
1555 for (
size_t i=0; i<static_cast<size_t>(inArray1.dimension(0)); i++)
1556 for (
size_t j=0; j<static_cast<size_t>(inArray1.dimension(1)); j++)
1557 for (
size_t k=0; k<static_cast<size_t>(inArray1.dimension(2)); k++)
1558 for (
size_t l=0; l<static_cast<size_t>(inArray1.dimension(3)); l++){
1559 sumArrayWrap(i,j,k,l) = inArray1Wrap(i,j,k,l)+inArray2Wrap(i,j,k,l);
1561 }
else if(inArrayRank==3){
1562 for (
size_t i=0; i<static_cast<size_t>(inArray1.dimension(0)); i++)
1563 for (
size_t j=0; j<static_cast<size_t>(inArray1.dimension(1)); j++)
1564 for (
size_t k=0; k<static_cast<size_t>(inArray1.dimension(2)); k++){
1565 sumArrayWrap(i,j,k) = inArray1Wrap(i,j,k)+inArray2Wrap(i,j,k);
1567 }
else if(inArrayRank==2){
1568 for (
size_t i=0; i<static_cast<size_t>(inArray1.dimension(0)); i++)
1569 for (
size_t j=0; j<static_cast<size_t>(inArray1.dimension(1)); j++){
1570 sumArrayWrap(i,j) = inArray1Wrap(i,j)+inArray2Wrap(i,j);
1572 }
else if(inArrayRank==1){
1573 for (
size_t i=0; i<static_cast<size_t>(inArray1.dimension(0)); i++){
1574 sumArrayWrap(i) = inArray1Wrap(i)+inArray2Wrap(i);
1585#ifdef HAVE_INTREPID_DEBUG
1586 TEUCHOS_TEST_FOR_EXCEPTION( ( getrank(inArray) != getrank(inoutSumArray) ),
1587 std::invalid_argument,
1588 ">>> ERROR (RealSpaceTools::add): Array arguments must have identical ranks!");
1589 for (
size_t i=0; i<getrank(inArray); i++) {
1590 TEUCHOS_TEST_FOR_EXCEPTION( (
static_cast<size_t>(inArray.dimension(i)) !=
static_cast<size_t>(inoutSumArray.dimension(i)) ),
1591 std::invalid_argument,
1592 ">>> ERROR (RealSpaceTools::add): Dimensions of array arguments do not agree!");
1596 ArrayWrapper<Scalar,ArraySum, Rank<ArraySum >::value,
false>inoutSumArrayWrap(inoutSumArray);
1597 ArrayWrapper<Scalar,ArrayIn, Rank<ArrayIn >::value,
true>inArrayWrap(inArray);
1598 int inArrayRank=getrank(inArray);
1602 for (
size_t i=0; i<static_cast<size_t>(
static_cast<size_t>(inArray.dimension(0))); i++)
1603 for (
size_t j=0; j<static_cast<size_t>(
static_cast<size_t>(inArray.dimension(1))); j++)
1604 for (
size_t k=0; k<static_cast<size_t>(
static_cast<size_t>(inArray.dimension(2))); k++)
1605 for (
size_t l=0; l<static_cast<size_t>(
static_cast<size_t>(inArray.dimension(3))); l++)
1606 for (
size_t m=0; m<static_cast<size_t>(
static_cast<size_t>(inArray.dimension(4))); m++){
1607 inoutSumArrayWrap(i,j,k,l,m) += inArrayWrap(i,j,k,l,m);
1609 }
else if(inArrayRank==4){
1610 for (
size_t i=0; i<static_cast<size_t>(
static_cast<size_t>(inArray.dimension(0))); i++)
1611 for (
size_t j=0; j<static_cast<size_t>(
static_cast<size_t>(inArray.dimension(1))); j++)
1612 for (
size_t k=0; k<static_cast<size_t>(
static_cast<size_t>(inArray.dimension(2))); k++)
1613 for (
size_t l=0; l<static_cast<size_t>(
static_cast<size_t>(inArray.dimension(3))); l++){
1614 inoutSumArrayWrap(i,j,k,l) += inArrayWrap(i,j,k,l);
1616 }
else if(inArrayRank==3){
1617 for (
size_t i=0; i<static_cast<size_t>(
static_cast<size_t>(inArray.dimension(0))); i++)
1618 for (
size_t j=0; j<static_cast<size_t>(
static_cast<size_t>(inArray.dimension(1))); j++)
1619 for (
size_t k=0; k<static_cast<size_t>(inArray.dimension(2)); k++){
1620 inoutSumArrayWrap(i,j,k) += inArrayWrap(i,j,k);
1622 }
else if(inArrayRank==2){
1623 for (
size_t i=0; i<static_cast<size_t>(
static_cast<size_t>(inArray.dimension(0))); i++)
1624 for (
size_t j=0; j<static_cast<size_t>(
static_cast<size_t>(inArray.dimension(1))); j++){
1625 inoutSumArrayWrap(i,j) += inArrayWrap(i,j);
1627 }
else if(inArrayRank==1){
1628 for (
size_t i=0; i<static_cast<size_t>(
static_cast<size_t>(inArray.dimension(0))); i++){
1629 inoutSumArrayWrap(i) += inArrayWrap(i);
1658 ArrayWrapper<Scalar,ArrayDiff, Rank<ArrayDiff >::value,
false>diffArrayWrap(diffArray);
1659 ArrayWrapper<Scalar,ArrayIn1, Rank<ArrayIn1 >::value,
true>inArray1Wrap(inArray1);
1660 ArrayWrapper<Scalar,ArrayIn2, Rank<ArrayIn2 >::value,
true>inArray2Wrap(inArray2);
1661 size_t inArray1Rank=getrank(inArray1);
1662#ifdef HAVE_INTREPID_DEBUG
1663 TEUCHOS_TEST_FOR_EXCEPTION( ( (getrank(inArray1) != getrank(inArray2)) || (getrank(inArray1) != getrank(diffArray)) ),
1664 std::invalid_argument,
1665 ">>> ERROR (RealSpaceTools::subtract): Array arguments must have identical ranks!");
1666 for (
size_t i=0; i<getrank(inArray1); i++) {
1667 TEUCHOS_TEST_FOR_EXCEPTION( ( (
static_cast<size_t>(inArray1.dimension(i)) !=
static_cast<size_t>(inArray2.dimension(i))) || (
static_cast<size_t>(inArray1.dimension(i)) !=
static_cast<size_t>(diffArray.dimension(i))) ),
1668 std::invalid_argument,
1669 ">>> ERROR (RealSpaceTools::subtract): Dimensions of array arguments do not agree!");
1672 if(inArray1Rank==5){
1673 for (
size_t i=0; i<static_cast<size_t>(inArray1.dimension(0)); i++)
1674 for (
size_t j=0; j<static_cast<size_t>(inArray1.dimension(1)); j++)
1675 for (
size_t k=0; k<static_cast<size_t>(inArray1.dimension(2)); k++)
1676 for (
size_t l=0; l<static_cast<size_t>(inArray1.dimension(3)); l++)
1677 for (
size_t m=0; m<static_cast<size_t>(inArray1.dimension(4)); m++){
1678 diffArrayWrap(i,j,k,l,m) = inArray1Wrap(i,j,k,l,m)-inArray2Wrap(i,j,k,l,m);
1680 }
else if(inArray1Rank==4){
1681 for (
size_t i=0; i<static_cast<size_t>(inArray1.dimension(0)); i++)
1682 for (
size_t j=0; j<static_cast<size_t>(inArray1.dimension(1)); j++)
1683 for (
size_t k=0; k<static_cast<size_t>(inArray1.dimension(2)); k++)
1684 for (
size_t l=0; l<static_cast<size_t>(inArray1.dimension(3)); l++){
1685 diffArrayWrap(i,j,k,l) = inArray1Wrap(i,j,k,l)-inArray2Wrap(i,j,k,l);
1687 }
else if(inArray1Rank==3){
1688 for (
size_t i=0; i<static_cast<size_t>(inArray1.dimension(0)); i++)
1689 for (
size_t j=0; j<static_cast<size_t>(inArray1.dimension(1)); j++)
1690 for (
size_t k=0; k<static_cast<size_t>(inArray1.dimension(2)); k++){
1691 diffArrayWrap(i,j,k) = inArray1Wrap(i,j,k)-inArray2Wrap(i,j,k);
1693 }
else if(inArray1Rank==2){
1694 for (
size_t i=0; i<static_cast<size_t>(inArray1.dimension(0)); i++)
1695 for (
size_t j=0; j<static_cast<size_t>(inArray1.dimension(1)); j++){
1696 diffArrayWrap(i,j) = inArray1Wrap(i,j)-inArray2Wrap(i,j);
1698 }
else if(inArray1Rank==1){
1699 for (
size_t i=0; i<static_cast<size_t>(inArray1.dimension(0)); i++){
1700 diffArrayWrap(i) = inArray1Wrap(i)-inArray2Wrap(i);
1711 ArrayWrapper<Scalar,ArrayDiff, Rank<ArrayDiff >::value,
false>inoutDiffArrayWrap(inoutDiffArray);
1712 ArrayWrapper<Scalar,ArrayIn, Rank<ArrayIn >::value,
true>inArrayWrap(inArray);
1713 int inArrayRank=getrank(inArray);
1714#ifdef HAVE_INTREPID_DEBUG
1715 TEUCHOS_TEST_FOR_EXCEPTION( ( getrank(inArray) != getrank(inoutDiffArray) ),
1716 std::invalid_argument,
1717 ">>> ERROR (RealSpaceTools::subtract): Array arguments must have identical ranks!");
1718 for (
size_t i=0; i<getrank(inArray); i++) {
1719 TEUCHOS_TEST_FOR_EXCEPTION( (
static_cast<size_t>(inArray.dimension(i)) !=
static_cast<size_t>(inoutDiffArray.dimension(i)) ),
1720 std::invalid_argument,
1721 ">>> ERROR (RealSpaceTools::subtract): Dimensions of array arguments do not agree!");
1726 for (
size_t i=0; i<static_cast<size_t>(
static_cast<size_t>(inArray.dimension(0))); i++)
1727 for (
size_t j=0; j<static_cast<size_t>(
static_cast<size_t>(inArray.dimension(1))); j++)
1728 for (
size_t k=0; k<static_cast<size_t>(
static_cast<size_t>(inArray.dimension(2))); k++)
1729 for (
size_t l=0; l<static_cast<size_t>(
static_cast<size_t>(inArray.dimension(3))); l++)
1730 for (
size_t m=0; m<static_cast<size_t>(
static_cast<size_t>(inArray.dimension(4))); m++){
1731 inoutDiffArrayWrap(i,j,k,l,m) -= inArrayWrap(i,j,k,l,m);
1733 }
else if(inArrayRank==4){
1734 for (
size_t i=0; i<static_cast<size_t>(
static_cast<size_t>(inArray.dimension(0))); i++)
1735 for (
size_t j=0; j<static_cast<size_t>(
static_cast<size_t>(inArray.dimension(1))); j++)
1736 for (
size_t k=0; k<static_cast<size_t>(
static_cast<size_t>(inArray.dimension(2))); k++)
1737 for (
size_t l=0; l<static_cast<size_t>(
static_cast<size_t>(inArray.dimension(3))); l++){
1738 inoutDiffArrayWrap(i,j,k,l) -= inArrayWrap(i,j,k,l);
1740 }
else if(inArrayRank==3){
1741 for (
size_t i=0; i<static_cast<size_t>(
static_cast<size_t>(inArray.dimension(0))); i++)
1742 for (
size_t j=0; j<static_cast<size_t>(
static_cast<size_t>(inArray.dimension(1))); j++)
1743 for (
size_t k=0; k<static_cast<size_t>(
static_cast<size_t>(inArray.dimension(2))); k++){
1744 inoutDiffArrayWrap(i,j,k) -= inArrayWrap(i,j,k);
1746 }
else if(inArrayRank==2){
1747 for (
size_t i=0; i<static_cast<size_t>(
static_cast<size_t>(inArray.dimension(0))); i++)
1748 for (
size_t j=0; j<static_cast<size_t>(
static_cast<size_t>(inArray.dimension(1))); j++){
1749 inoutDiffArrayWrap(i,j) -= inArrayWrap(i,j);
1751 }
else if(inArrayRank==1){
1752 for (
size_t i=0; i<static_cast<size_t>(
static_cast<size_t>(inArray.dimension(0))); i++){
1753 inoutDiffArrayWrap(i) -= inArrayWrap(i);
1782#ifdef HAVE_INTREPID_DEBUG
1783 TEUCHOS_TEST_FOR_EXCEPTION( ( getrank(inArray) != getrank(scaledArray) ),
1784 std::invalid_argument,
1785 ">>> ERROR (RealSpaceTools::scale): Array arguments must have identical ranks!");
1786 for (
size_t i=0; i<getrank(inArray); i++) {
1787 TEUCHOS_TEST_FOR_EXCEPTION( (
static_cast<size_t>(inArray.dimension(i)) !=
static_cast<size_t>(scaledArray.dimension(i)) ),
1788 std::invalid_argument,
1789 ">>> ERROR (RealSpaceTools::scale): Dimensions of array arguments do not agree!");
1794 int inArrayRank=getrank(inArray);
1795 ArrayWrapper<Scalar,ArrayScaled, Rank<ArrayScaled >::value,
false>scaledArrayWrap(scaledArray);
1796 ArrayWrapper<Scalar,ArrayIn, Rank<ArrayIn >::value,
true>inArrayWrap(inArray);
1798 for (
size_t i=0; i<static_cast<size_t>(
static_cast<size_t>(inArray.dimension(0))); i++)
1799 for (
size_t j=0; j<static_cast<size_t>(
static_cast<size_t>(inArray.dimension(1))); j++)
1800 for (
size_t k=0; k<static_cast<size_t>(
static_cast<size_t>(inArray.dimension(2))); k++)
1801 for (
size_t l=0; l<static_cast<size_t>(
static_cast<size_t>(inArray.dimension(3))); l++)
1802 for (
size_t m=0; m<static_cast<size_t>(
static_cast<size_t>(inArray.dimension(4))); m++){
1803 scaledArrayWrap(i,j,k,l,m) = scalar*inArrayWrap(i,j,k,l,m);
1805 }
else if(inArrayRank==4){
1806 for (
size_t i=0; i<static_cast<size_t>(
static_cast<size_t>(inArray.dimension(0))); i++)
1807 for (
size_t j=0; j<static_cast<size_t>(
static_cast<size_t>(inArray.dimension(1))); j++)
1808 for (
size_t k=0; k<static_cast<size_t>(
static_cast<size_t>(inArray.dimension(2))); k++)
1809 for (
size_t l=0; l<static_cast<size_t>(
static_cast<size_t>(inArray.dimension(3))); l++){
1810 scaledArrayWrap(i,j,k,l) = scalar*inArrayWrap(i,j,k,l);
1812 }
else if(inArrayRank==3){
1813 for (
size_t i=0; i<static_cast<size_t>(
static_cast<size_t>(inArray.dimension(0))); i++)
1814 for (
size_t j=0; j<static_cast<size_t>(
static_cast<size_t>(inArray.dimension(1))); j++)
1815 for (
size_t k=0; k<static_cast<size_t>(
static_cast<size_t>(inArray.dimension(2))); k++){
1816 scaledArrayWrap(i,j,k) = scalar*inArrayWrap(i,j,k);
1818 }
else if(inArrayRank==2){
1819 for (
size_t i=0; i<static_cast<size_t>(
static_cast<size_t>(inArray.dimension(0))); i++)
1820 for (
size_t j=0; j<static_cast<size_t>(
static_cast<size_t>(inArray.dimension(1))); j++){
1821 scaledArrayWrap(i,j) = scalar*inArrayWrap(i,j);
1823 }
else if(inArrayRank==1){
1824 for (
size_t i=0; i<static_cast<size_t>(
static_cast<size_t>(inArray.dimension(0))); i++){
1825 scaledArrayWrap(i) = scalar*inArrayWrap(i);
2040 ArrayWrapper<Scalar,ArrayVecProd, Rank<ArrayVecProd >::value,
false>vecProdWrap(vecProd);
2041 ArrayWrapper<Scalar,ArrayIn1, Rank<ArrayIn1 >::value,
true>inLeftWrap(inLeft);
2042 ArrayWrapper<Scalar,ArrayIn2, Rank<ArrayIn2 >::value,
true>inRightWrap(inRight);
2043#ifdef HAVE_INTREPID_DEBUG
2049 std::string errmsg =
">>> ERROR (RealSpaceTools::vecprod):";
2052 TEUCHOS_TEST_FOR_EXCEPTION( !
requireRankRange(errmsg, inLeft, 1,3), std::invalid_argument, errmsg);
2053 TEUCHOS_TEST_FOR_EXCEPTION( !
requireDimensionMatch(errmsg, inLeft, inRight), std::invalid_argument, errmsg);
2054 TEUCHOS_TEST_FOR_EXCEPTION( !
requireDimensionMatch(errmsg, inLeft, vecProd), std::invalid_argument, errmsg);
2058 TEUCHOS_TEST_FOR_EXCEPTION( !
requireDimensionRange(errmsg, inLeft, getrank(inLeft) - 1, 2,3), std::invalid_argument, errmsg);
2062 int spaceDim =
static_cast<size_t>(inLeft.dimension(getrank(inLeft) - 1));
2064 switch(getrank(inLeft) ){
2068 vecProdWrap(0) = inLeftWrap(1)*inRightWrap(2) - inLeftWrap(2)*inRightWrap(1);
2069 vecProdWrap(1) = inLeftWrap(2)*inRightWrap(0) - inLeftWrap(0)*inRightWrap(2);
2070 vecProdWrap(2) = inLeftWrap(0)*inRightWrap(1) - inLeftWrap(1)*inRightWrap(0);
2076 size_t dim0 =
static_cast<size_t>(inLeft.dimension(0));
2078 for(
size_t i0 = 0; i0 < dim0; i0++){
2079 vecProdWrap(i0, 0) = inLeftWrap(i0, 1)*inRightWrap(i0, 2) - inLeftWrap(i0, 2)*inRightWrap(i0, 1);
2080 vecProdWrap(i0, 1) = inLeftWrap(i0, 2)*inRightWrap(i0, 0) - inLeftWrap(i0, 0)*inRightWrap(i0, 2);
2081 vecProdWrap(i0, 2) = inLeftWrap(i0, 0)*inRightWrap(i0, 1) - inLeftWrap(i0, 1)*inRightWrap(i0, 0);
2084 else if(spaceDim == 2){
2085 for(
size_t i0 = 0; i0 < dim0; i0++){
2087 vecProdWrap(i0, 0) = inLeftWrap(i0, 0)*inRightWrap(i0, 1) - inLeftWrap(i0, 1)*inRightWrap(i0, 0);
2095 size_t dim0 =
static_cast<size_t>(inLeft.dimension(0));
2096 size_t dim1 =
static_cast<size_t>(inLeft.dimension(1));
2098 for(
size_t i0 = 0; i0 < dim0; i0++){
2099 for(
size_t i1 = 0; i1 < dim1; i1++){
2100 vecProdWrap(i0, i1, 0) = inLeftWrap(i0, i1, 1)*inRightWrap(i0, i1, 2) - inLeftWrap(i0, i1, 2)*inRightWrap(i0, i1, 1);
2101 vecProdWrap(i0, i1, 1) = inLeftWrap(i0, i1, 2)*inRightWrap(i0, i1, 0) - inLeftWrap(i0, i1, 0)*inRightWrap(i0, i1, 2);
2102 vecProdWrap(i0, i1, 2) = inLeftWrap(i0, i1, 0)*inRightWrap(i0, i1, 1) - inLeftWrap(i0, i1, 1)*inRightWrap(i0, i1, 0);
2106 else if(spaceDim == 2){
2107 for(
size_t i0 = 0; i0 < dim0; i0++){
2108 for(
size_t i1 = 0; i1 < dim1; i1++){
2110 vecProdWrap(i0, i1, 0) = inLeftWrap(i0, i1, 0)*inRightWrap(i0, i1, 1) - inLeftWrap(i0, i1, 1)*inRightWrap(i0, i1, 0);
2118 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::invalid_argument,
2119 ">>> ERROR (RealSpaceTools::vecprod): rank-1,2,3 arrays required");