type_tensor.h
Go to the documentation of this file.00001 // The libMesh Finite Element Library. 00002 // Copyright (C) 2002-2012 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner 00003 00004 // This library is free software; you can redistribute it and/or 00005 // modify it under the terms of the GNU Lesser General Public 00006 // License as published by the Free Software Foundation; either 00007 // version 2.1 of the License, or (at your option) any later version. 00008 00009 // This library is distributed in the hope that it will be useful, 00010 // but WITHOUT ANY WARRANTY; without even the implied warranty of 00011 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 00012 // Lesser General Public License for more details. 00013 00014 // You should have received a copy of the GNU Lesser General Public 00015 // License along with this library; if not, write to the Free Software 00016 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA 00017 00018 00019 00020 #ifndef LIBMESH_TYPE_TENSOR_H 00021 #define LIBMESH_TYPE_TENSOR_H 00022 00023 // Local includes 00024 #include "libmesh/libmesh_common.h" 00025 #include "libmesh/type_vector.h" 00026 00027 // C++ includes 00028 #include <cstdlib> // *must* precede <cmath> for proper std:abs() on PGI, Sun Studio CC 00029 #include <cmath> 00030 00031 namespace libMesh 00032 { 00033 00034 // Forward declarations 00035 template <typename T> class TypeTensorColumn; 00036 template <typename T> class ConstTypeTensorColumn; 00037 template <unsigned int N, typename T> class TypeNTensor; 00038 00039 00040 00050 template <typename T> 00051 class TypeTensor 00052 { 00053 template <typename T2> 00054 friend class TypeTensor; 00055 00056 protected: 00057 00062 TypeTensor (); 00063 00070 explicit TypeTensor (const T xx, 00071 const T xy=0, 00072 const T xz=0, 00073 const T yx=0, 00074 const T yy=0, 00075 const T yz=0, 00076 const T zx=0, 00077 const T zy=0, 00078 const T zz=0); 00079 00080 00084 template <typename Scalar> 00085 explicit TypeTensor (const Scalar xx, 00086 const Scalar xy=0, 00087 const Scalar xz=0, 00088 const Scalar yx=0, 00089 const Scalar yy=0, 00090 const Scalar yz=0, 00091 const Scalar zx=0, 00092 const Scalar zy=0, 00093 typename 00094 boostcopy::enable_if_c<ScalarTraits<Scalar>::value, 00095 const Scalar>::type zz=0); 00096 00102 template<typename T2> 00103 TypeTensor(const TypeVector<T2>& vx); 00104 00105 template<typename T2> 00106 TypeTensor(const TypeVector<T2>& vx, const TypeVector<T2> &vy); 00107 00108 template<typename T2> 00109 TypeTensor(const TypeVector<T2>& vx, const TypeVector<T2> &vy, const TypeVector<T2> &vz); 00110 00111 public: 00112 00116 template<typename T2> 00117 TypeTensor(const TypeTensor<T2>& p); 00118 00122 ~TypeTensor(); 00123 00127 template<typename T2> 00128 void assign (const TypeTensor<T2> &); 00129 00133 template <typename Scalar> 00134 typename boostcopy::enable_if_c< 00135 ScalarTraits<Scalar>::value, 00136 TypeTensor&>::type 00137 operator = (const Scalar& p) 00138 { libmesh_assert_equal_to (p, Scalar(0)); this->zero(); return *this; } 00139 00143 const T & operator () (const unsigned int i, const unsigned int j) const; 00144 00149 T & operator () (const unsigned int i, const unsigned int j); 00150 00154 ConstTypeTensorColumn<T> slice (const unsigned int i) const; 00155 00159 TypeTensorColumn<T> slice (const unsigned int i); 00160 00164 TypeVector<T> row(const unsigned int r); 00165 00169 template<typename T2> 00170 TypeTensor<typename CompareTypes<T, T2>::supertype> 00171 operator + (const TypeTensor<T2> &) const; 00172 00176 template<typename T2> 00177 const TypeTensor<T> & operator += (const TypeTensor<T2> &); 00178 00182 template<typename T2> 00183 void add (const TypeTensor<T2> &); 00184 00189 template <typename T2> 00190 void add_scaled (const TypeTensor<T2> &, const T); 00191 00195 template<typename T2> 00196 TypeTensor<typename CompareTypes<T, T2>::supertype> 00197 operator - (const TypeTensor<T2> &) const; 00198 00202 template<typename T2> 00203 const TypeTensor<T> & operator -= (const TypeTensor<T2> &); 00204 00208 template<typename T2> 00209 void subtract (const TypeTensor<T2> &); 00210 00215 template <typename T2> 00216 void subtract_scaled (const TypeTensor<T2> &, const T); 00217 00221 TypeTensor<T> operator - () const; 00222 00226 template <typename Scalar> 00227 typename boostcopy::enable_if_c< 00228 ScalarTraits<Scalar>::value, 00229 TypeTensor<typename CompareTypes<T, Scalar>::supertype> >::type 00230 operator * (const Scalar) const; 00231 00235 template <typename Scalar> 00236 const TypeTensor<T> & operator *= (const Scalar); 00237 00241 template <typename Scalar> 00242 typename boostcopy::enable_if_c< 00243 ScalarTraits<Scalar>::value, 00244 TypeTensor<typename CompareTypes<T, Scalar>::supertype> >::type 00245 operator / (const Scalar) const; 00246 00250 const TypeTensor<T> & operator /= (const T); 00251 00256 template <typename T2> 00257 TypeTensor<T> operator * (const TypeTensor<T2> &) const; 00258 00264 template <typename T2> 00265 typename CompareTypes<T,T2>::supertype 00266 contract (const TypeTensor<T2> &) const; 00267 00272 template <typename T2> 00273 TypeVector<typename CompareTypes<T,T2>::supertype> 00274 operator * (const TypeVector<T2> &) const; 00275 00279 TypeTensor<T> transpose() const; 00280 00285 Real size() const; 00286 00291 Real size_sq() const; 00292 00298 T det() const; 00299 00303 T tr() const; 00304 00308 void zero(); 00309 00313 bool operator == (const TypeTensor<T>& rhs) const; 00314 00319 bool operator < (const TypeTensor<T>& rhs) const; 00320 00325 bool operator > (const TypeTensor<T>& rhs) const; 00326 00330 void print(std::ostream& os = libMesh::out) const; 00331 00336 friend std::ostream& operator << (std::ostream& os, const TypeTensor<T>& t) 00337 { 00338 t.print(os); 00339 return os; 00340 } 00341 00346 void write_unformatted (std::ostream &out, const bool newline = true) const; 00347 00348 // protected: 00349 00353 T _coords[LIBMESH_DIM*LIBMESH_DIM]; 00354 }; 00355 00356 00357 template <typename T> 00358 class TypeTensorColumn 00359 { 00360 public: 00361 00362 TypeTensorColumn(TypeTensor<T> &tensor, 00363 unsigned int j) : 00364 _tensor(&tensor), _j(j) {} 00365 00370 T & operator () (const unsigned int i) 00371 { return (*_tensor)(i,_j); } 00372 00373 T & slice (const unsigned int i) 00374 { return (*_tensor)(i,_j); } 00375 00379 TypeTensorColumn<T>& operator = (const TypeVector<T>& rhs) 00380 { 00381 for (unsigned int i=0; i != LIBMESH_DIM; ++i) 00382 (*this)(i) = rhs(i); 00383 return *this; 00384 } 00385 00386 private: 00387 TypeTensor<T> *_tensor; 00388 const unsigned int _j; 00389 }; 00390 00391 00392 template <typename T> 00393 class ConstTypeTensorColumn 00394 { 00395 public: 00396 00397 ConstTypeTensorColumn(const TypeTensor<T> &tensor, 00398 unsigned int j) : 00399 _tensor(&tensor), _j(j) {} 00400 00404 const T & operator () (const unsigned int i) const 00405 { return (*_tensor)(i,_j); } 00406 00407 const T & slice (const unsigned int i) const 00408 { return (*_tensor)(i,_j); } 00409 00410 private: 00411 const TypeTensor<T> *_tensor; 00412 const unsigned int _j; 00413 }; 00414 00415 //------------------------------------------------------ 00416 // Inline functions 00417 template <typename T> 00418 inline 00419 TypeTensor<T>::TypeTensor () 00420 { 00421 _coords[0] = 0; 00422 00423 #if LIBMESH_DIM > 1 00424 _coords[1] = 0; 00425 _coords[2] = 0; 00426 _coords[3] = 0; 00427 #endif 00428 00429 #if LIBMESH_DIM > 2 00430 _coords[4] = 0; 00431 _coords[5] = 0; 00432 _coords[6] = 0; 00433 _coords[7] = 0; 00434 _coords[8] = 0; 00435 #endif 00436 } 00437 00438 00439 00440 template <typename T> 00441 inline 00442 TypeTensor<T>::TypeTensor 00443 (const T xx, 00444 const T xy, 00445 const T xz, 00446 const T yx, 00447 const T yy, 00448 const T yz, 00449 const T zx, 00450 const T zy, 00451 const T zz) 00452 { 00453 _coords[0] = xx; 00454 00455 #if LIBMESH_DIM == 2 00456 _coords[1] = xy; 00457 _coords[2] = yx; 00458 _coords[3] = yy; 00459 #endif 00460 00461 #if LIBMESH_DIM == 3 00462 _coords[1] = xy; 00463 _coords[2] = xz; 00464 _coords[3] = yx; 00465 _coords[4] = yy; 00466 _coords[5] = yz; 00467 _coords[6] = zx; 00468 _coords[7] = zy; 00469 _coords[8] = zz; 00470 #endif 00471 } 00472 00473 00474 template <typename T> 00475 template <typename Scalar> 00476 inline 00477 TypeTensor<T>::TypeTensor 00478 (const Scalar xx, 00479 const Scalar xy, 00480 const Scalar xz, 00481 const Scalar yx, 00482 const Scalar yy, 00483 const Scalar yz, 00484 const Scalar zx, 00485 const Scalar zy, 00486 typename 00487 boostcopy::enable_if_c<ScalarTraits<Scalar>::value, 00488 const Scalar>::type zz) 00489 { 00490 _coords[0] = xx; 00491 00492 #if LIBMESH_DIM == 2 00493 _coords[1] = xy; 00494 _coords[2] = yx; 00495 _coords[3] = yy; 00496 #endif 00497 00498 #if LIBMESH_DIM == 3 00499 _coords[1] = xy; 00500 _coords[2] = xz; 00501 _coords[3] = yx; 00502 _coords[4] = yy; 00503 _coords[5] = yz; 00504 _coords[6] = zx; 00505 _coords[7] = zy; 00506 _coords[8] = zz; 00507 #endif 00508 } 00509 00510 00511 00512 template <typename T> 00513 template<typename T2> 00514 inline 00515 TypeTensor<T>::TypeTensor (const TypeTensor <T2> &p) 00516 { 00517 // copy the nodes from vector p to me 00518 for (unsigned int i=0; i<LIBMESH_DIM*LIBMESH_DIM; i++) 00519 _coords[i] = p._coords[i]; 00520 } 00521 00522 00523 template <typename T> 00524 template <typename T2> 00525 TypeTensor<T>::TypeTensor(const TypeVector<T2>& vx) 00526 { 00527 libmesh_assert_equal_to (LIBMESH_DIM, 1); 00528 _coords[0] = vx(0); 00529 } 00530 00531 template <typename T> 00532 template <typename T2> 00533 TypeTensor<T>::TypeTensor(const TypeVector<T2>& vx, const TypeVector<T2> &vy) 00534 { 00535 libmesh_assert_equal_to (LIBMESH_DIM, 2); 00536 _coords[0] = vx(0); 00537 _coords[1] = vx(1); 00538 _coords[2] = vy(0); 00539 _coords[3] = vy(1); 00540 } 00541 00542 template <typename T> 00543 template <typename T2> 00544 TypeTensor<T>::TypeTensor(const TypeVector<T2>& vx, const TypeVector<T2> &vy, const TypeVector<T2> &vz) 00545 { 00546 libmesh_assert_equal_to (LIBMESH_DIM, 3); 00547 _coords[0] = vx(0); 00548 _coords[1] = vx(1); 00549 _coords[2] = vx(2); 00550 _coords[3] = vy(0); 00551 _coords[4] = vy(1); 00552 _coords[5] = vy(2); 00553 _coords[6] = vz(0); 00554 _coords[7] = vz(1); 00555 _coords[8] = vz(2); 00556 } 00557 00558 00559 00560 00561 template <typename T> 00562 inline 00563 TypeTensor<T>::~TypeTensor () 00564 { 00565 } 00566 00567 00568 00569 template <typename T> 00570 template<typename T2> 00571 inline 00572 void TypeTensor<T>::assign (const TypeTensor<T2> &p) 00573 { 00574 for (unsigned int i=0; i<LIBMESH_DIM*LIBMESH_DIM; i++) 00575 _coords[i] = p._coords[i]; 00576 } 00577 00578 00579 00580 template <typename T> 00581 inline 00582 const T & TypeTensor<T>::operator () (const unsigned int i, 00583 const unsigned int j) const 00584 { 00585 libmesh_assert_less (i, 3); 00586 libmesh_assert_less (j, 3); 00587 00588 #if LIBMESH_DIM < 3 00589 const static T my_zero = 0; 00590 if (i >= LIBMESH_DIM || j >= LIBMESH_DIM) 00591 return my_zero; 00592 #endif 00593 00594 return _coords[i*LIBMESH_DIM+j]; 00595 } 00596 00597 00598 00599 template <typename T> 00600 inline 00601 T & TypeTensor<T>::operator () (const unsigned int i, 00602 const unsigned int j) 00603 { 00604 #if LIBMESH_DIM < 3 00605 00606 if (i >= LIBMESH_DIM || j >= LIBMESH_DIM) 00607 { 00608 // libMesh::err << "ERROR: You are assigning to a tensor component" << std::endl 00609 // << "that is out of range for the compiled LIBMESH_DIM!" << std::endl 00610 // << " LIBMESH_DIM=" << LIBMESH_DIM << " , i=" << i << " , j=" << j << std::endl; 00611 libmesh_error(); 00612 } 00613 00614 #endif 00615 00616 libmesh_assert_less (i, LIBMESH_DIM); 00617 libmesh_assert_less (j, LIBMESH_DIM); 00618 00619 return _coords[i*LIBMESH_DIM+j]; 00620 } 00621 00622 00623 template <typename T> 00624 inline 00625 ConstTypeTensorColumn<T> 00626 TypeTensor<T>::slice (const unsigned int i) const 00627 { 00628 libmesh_assert_less (i, LIBMESH_DIM); 00629 return ConstTypeTensorColumn<T>(*this, i); 00630 } 00631 00632 00633 template <typename T> 00634 inline 00635 TypeTensorColumn<T> 00636 TypeTensor<T>::slice (const unsigned int i) 00637 { 00638 libmesh_assert_less (i, LIBMESH_DIM); 00639 return TypeTensorColumn<T>(*this, i); 00640 } 00641 00642 00643 template <typename T> 00644 inline 00645 TypeVector<T> 00646 TypeTensor<T>::row(const unsigned int r) 00647 { 00648 TypeVector<T> return_vector; 00649 00650 for(unsigned int j=0; j<LIBMESH_DIM; j++) 00651 return_vector._coords[j] = _coords[r*LIBMESH_DIM + j]; 00652 00653 return return_vector; 00654 } 00655 00656 template <typename T> 00657 template<typename T2> 00658 inline 00659 TypeTensor<typename CompareTypes<T, T2>::supertype> 00660 TypeTensor<T>::operator + (const TypeTensor<T2> &p) const 00661 { 00662 00663 #if LIBMESH_DIM == 1 00664 return TypeTensor(_coords[0] + p._coords[0]); 00665 #endif 00666 00667 #if LIBMESH_DIM == 2 00668 return TypeTensor(_coords[0] + p._coords[0], 00669 _coords[1] + p._coords[1], 00670 0., 00671 _coords[2] + p._coords[2], 00672 _coords[3] + p._coords[3]); 00673 #endif 00674 00675 #if LIBMESH_DIM == 3 00676 return TypeTensor(_coords[0] + p._coords[0], 00677 _coords[1] + p._coords[1], 00678 _coords[2] + p._coords[2], 00679 _coords[3] + p._coords[3], 00680 _coords[4] + p._coords[4], 00681 _coords[5] + p._coords[5], 00682 _coords[6] + p._coords[6], 00683 _coords[7] + p._coords[7], 00684 _coords[8] + p._coords[8]); 00685 #endif 00686 00687 } 00688 00689 00690 00691 template <typename T> 00692 template<typename T2> 00693 inline 00694 const TypeTensor<T> & TypeTensor<T>::operator += (const TypeTensor<T2> &p) 00695 { 00696 this->add (p); 00697 00698 return *this; 00699 } 00700 00701 00702 00703 template <typename T> 00704 template<typename T2> 00705 inline 00706 void TypeTensor<T>::add (const TypeTensor<T2> &p) 00707 { 00708 for (unsigned int i=0; i<LIBMESH_DIM*LIBMESH_DIM; i++) 00709 _coords[i] += p._coords[i]; 00710 } 00711 00712 00713 00714 template <typename T> 00715 template <typename T2> 00716 inline 00717 void TypeTensor<T>::add_scaled (const TypeTensor<T2> &p, const T factor) 00718 { 00719 for (unsigned int i=0; i<LIBMESH_DIM*LIBMESH_DIM; i++) 00720 _coords[i] += factor*p._coords[i]; 00721 00722 } 00723 00724 00725 00726 template <typename T> 00727 template<typename T2> 00728 inline 00729 TypeTensor<typename CompareTypes<T, T2>::supertype> 00730 TypeTensor<T>::operator - (const TypeTensor<T2> &p) const 00731 { 00732 00733 #if LIBMESH_DIM == 1 00734 return TypeTensor(_coords[0] - p._coords[0]); 00735 #endif 00736 00737 #if LIBMESH_DIM == 2 00738 return TypeTensor(_coords[0] - p._coords[0], 00739 _coords[1] - p._coords[1], 00740 0., 00741 _coords[2] - p._coords[2], 00742 _coords[3] - p._coords[3]); 00743 #endif 00744 00745 #if LIBMESH_DIM == 3 00746 return TypeTensor(_coords[0] - p._coords[0], 00747 _coords[1] - p._coords[1], 00748 _coords[2] - p._coords[2], 00749 _coords[3] - p._coords[3], 00750 _coords[4] - p._coords[4], 00751 _coords[5] - p._coords[5], 00752 _coords[6] - p._coords[6], 00753 _coords[7] - p._coords[7], 00754 _coords[8] - p._coords[8]); 00755 #endif 00756 00757 } 00758 00759 00760 00761 template <typename T> 00762 template <typename T2> 00763 inline 00764 const TypeTensor<T> & TypeTensor<T>::operator -= (const TypeTensor<T2> &p) 00765 { 00766 this->subtract (p); 00767 00768 return *this; 00769 } 00770 00771 00772 00773 template <typename T> 00774 template <typename T2> 00775 inline 00776 void TypeTensor<T>::subtract (const TypeTensor<T2>& p) 00777 { 00778 for (unsigned int i=0; i<LIBMESH_DIM*LIBMESH_DIM; i++) 00779 _coords[i] -= p._coords[i]; 00780 } 00781 00782 00783 00784 template <typename T> 00785 template <typename T2> 00786 inline 00787 void TypeTensor<T>::subtract_scaled (const TypeTensor<T2> &p, const T factor) 00788 { 00789 for (unsigned int i=0; i<LIBMESH_DIM*LIBMESH_DIM; i++) 00790 _coords[i] -= factor*p._coords[i]; 00791 } 00792 00793 00794 00795 template <typename T> 00796 inline 00797 TypeTensor<T> TypeTensor<T>::operator - () const 00798 { 00799 00800 #if LIBMESH_DIM == 1 00801 return TypeTensor(-_coords[0]); 00802 #endif 00803 00804 #if LIBMESH_DIM == 2 00805 return TypeTensor(-_coords[0], 00806 -_coords[1], 00807 -_coords[2], 00808 -_coords[3]); 00809 #endif 00810 00811 #if LIBMESH_DIM == 3 00812 return TypeTensor(-_coords[0], 00813 -_coords[1], 00814 -_coords[2], 00815 -_coords[3], 00816 -_coords[4], 00817 -_coords[5], 00818 -_coords[6], 00819 -_coords[7], 00820 -_coords[8]); 00821 #endif 00822 00823 } 00824 00825 00826 00827 template <typename T> 00828 template <typename Scalar> 00829 inline 00830 typename boostcopy::enable_if_c< 00831 ScalarTraits<Scalar>::value, 00832 TypeTensor<typename CompareTypes<T, Scalar>::supertype> >::type 00833 TypeTensor<T>::operator * (const Scalar factor) const 00834 { 00835 typedef typename CompareTypes<T, Scalar>::supertype TS; 00836 00837 00838 #if LIBMESH_DIM == 1 00839 return TypeTensor<TS>(_coords[0]*factor); 00840 #endif 00841 00842 #if LIBMESH_DIM == 2 00843 return TypeTensor<TS>(_coords[0]*factor, 00844 _coords[1]*factor, 00845 _coords[2]*factor, 00846 _coords[3]*factor); 00847 #endif 00848 00849 #if LIBMESH_DIM == 3 00850 return TypeTensor<TS>(_coords[0]*factor, 00851 _coords[1]*factor, 00852 _coords[2]*factor, 00853 _coords[3]*factor, 00854 _coords[4]*factor, 00855 _coords[5]*factor, 00856 _coords[6]*factor, 00857 _coords[7]*factor, 00858 _coords[8]*factor); 00859 #endif 00860 } 00861 00862 00863 00864 template <typename T, typename Scalar> 00865 inline 00866 typename boostcopy::enable_if_c< 00867 ScalarTraits<Scalar>::value, 00868 TypeTensor<typename CompareTypes<T, Scalar>::supertype> >::type 00869 operator * (const Scalar factor, 00870 const TypeTensor<T> &t) 00871 { 00872 return t * factor; 00873 } 00874 00875 00876 00877 template <typename T> 00878 template <typename Scalar> 00879 inline 00880 const TypeTensor<T> & TypeTensor<T>::operator *= (const Scalar factor) 00881 { 00882 for (unsigned int i=0; i<LIBMESH_DIM*LIBMESH_DIM; i++) 00883 _coords[i] *= factor; 00884 00885 return *this; 00886 } 00887 00888 00889 00890 00891 template <typename T> 00892 template <typename Scalar> 00893 inline 00894 typename boostcopy::enable_if_c< 00895 ScalarTraits<Scalar>::value, 00896 TypeTensor<typename CompareTypes<T, Scalar>::supertype> >::type 00897 TypeTensor<T>::operator / (const Scalar factor) const 00898 { 00899 libmesh_assert_not_equal_to (factor, static_cast<T>(0.)); 00900 00901 typedef typename CompareTypes<T, Scalar>::supertype TS; 00902 00903 #if LIBMESH_DIM == 1 00904 return TypeTensor<TS>(_coords[0]/factor); 00905 #endif 00906 00907 #if LIBMESH_DIM == 2 00908 return TypeTensor<TS>(_coords[0]/factor, 00909 _coords[1]/factor, 00910 _coords[2]/factor, 00911 _coords[3]/factor); 00912 #endif 00913 00914 #if LIBMESH_DIM == 3 00915 return TypeTensor<TS>(_coords[0]/factor, 00916 _coords[1]/factor, 00917 _coords[2]/factor, 00918 _coords[3]/factor, 00919 _coords[4]/factor, 00920 _coords[5]/factor, 00921 _coords[6]/factor, 00922 _coords[7]/factor, 00923 _coords[8]/factor); 00924 #endif 00925 00926 } 00927 00928 00929 00930 template <typename T> 00931 inline 00932 TypeTensor<T> TypeTensor<T>::transpose() const 00933 { 00934 #if LIBMESH_DIM == 1 00935 return TypeTensor(_coords[0]); 00936 #endif 00937 00938 #if LIBMESH_DIM == 2 00939 return TypeTensor(_coords[0], 00940 _coords[2], 00941 _coords[1], 00942 _coords[3]); 00943 #endif 00944 00945 #if LIBMESH_DIM == 3 00946 return TypeTensor(_coords[0], 00947 _coords[3], 00948 _coords[6], 00949 _coords[1], 00950 _coords[4], 00951 _coords[7], 00952 _coords[2], 00953 _coords[5], 00954 _coords[8]); 00955 #endif 00956 } 00957 00958 00959 00960 template <typename T> 00961 inline 00962 const TypeTensor<T> & TypeTensor<T>::operator /= (const T factor) 00963 { 00964 libmesh_assert_not_equal_to (factor, static_cast<T>(0.)); 00965 00966 for (unsigned int i=0; i<LIBMESH_DIM*LIBMESH_DIM; i++) 00967 _coords[i] /= factor; 00968 00969 return *this; 00970 } 00971 00972 00973 00974 00975 template <typename T> 00976 template <typename T2> 00977 inline 00978 TypeVector<typename CompareTypes<T,T2>::supertype> 00979 TypeTensor<T>::operator * (const TypeVector<T2> &p) const 00980 { 00981 TypeVector<typename CompareTypes<T,T2>::supertype> returnval; 00982 for (unsigned int i=0; i<LIBMESH_DIM; i++) 00983 for (unsigned int j=0; j<LIBMESH_DIM; j++) 00984 returnval(i) += (*this)(i,j)*p(j); 00985 00986 return returnval; 00987 } 00988 00989 00990 00991 template <typename T> 00992 template <typename T2> 00993 inline 00994 TypeTensor<T> TypeTensor<T>::operator * (const TypeTensor<T2> &p) const 00995 { 00996 TypeTensor<T> returnval; 00997 for (unsigned int i=0; i<LIBMESH_DIM; i++) 00998 for (unsigned int j=0; j<LIBMESH_DIM; j++) 00999 for (unsigned int k=0; k<LIBMESH_DIM; k++) 01000 returnval(i,j) += (*this)(i,k)*p(k,j); 01001 01002 return returnval; 01003 } 01004 01005 01006 01011 template <typename T> 01012 template <typename T2> 01013 inline 01014 typename CompareTypes<T,T2>::supertype 01015 TypeTensor<T>::contract (const TypeTensor<T2> &t) const 01016 { 01017 typename CompareTypes<T,T2>::supertype sum = 0.; 01018 for (unsigned int i=0; i<LIBMESH_DIM*LIBMESH_DIM; i++) 01019 sum += _coords[i]*t._coords[i]; 01020 return sum; 01021 } 01022 01023 01024 01025 template <typename T> 01026 inline 01027 Real TypeTensor<T>::size() const 01028 { 01029 return std::sqrt(this->size_sq()); 01030 } 01031 01032 01033 01034 template <typename T> 01035 inline 01036 T TypeTensor<T>::det() const 01037 { 01038 #if LIBMESH_DIM == 1 01039 return _coords[0]; 01040 #endif 01041 01042 #if LIBMESH_DIM == 2 01043 return (_coords[0] * _coords[3] 01044 - _coords[1] * _coords[2]); 01045 #endif 01046 01047 #if LIBMESH_DIM == 3 01048 return (_coords[0] * _coords[4] * _coords[8] 01049 + _coords[1] * _coords[5] * _coords[6] 01050 + _coords[2] * _coords[3] * _coords[7] 01051 - _coords[0] * _coords[5] * _coords[7] 01052 - _coords[1] * _coords[3] * _coords[8] 01053 - _coords[2] * _coords[4] * _coords[6]); 01054 #endif 01055 } 01056 01057 template <typename T> 01058 inline 01059 T TypeTensor<T>::tr() const 01060 { 01061 #if LIBMESH_DIM == 1 01062 return _coords[0]; 01063 #endif 01064 01065 #if LIBMESH_DIM == 2 01066 return _coords[0] + _coords[3]; 01067 #endif 01068 01069 #if LIBMESH_DIM == 3 01070 return _coords[0] + _coords[4] + _coords[8]; 01071 #endif 01072 } 01073 01074 template <typename T> 01075 inline 01076 void TypeTensor<T>::zero() 01077 { 01078 for (unsigned int i=0; i<LIBMESH_DIM*LIBMESH_DIM; i++) 01079 _coords[i] = 0.; 01080 } 01081 01082 01083 01084 template <typename T> 01085 inline 01086 Real TypeTensor<T>::size_sq () const 01087 { 01088 Real sum = 0.; 01089 for (unsigned int i=0; i<LIBMESH_DIM*LIBMESH_DIM; i++) 01090 sum += TensorTools::norm_sq(_coords[i]); 01091 return sum; 01092 } 01093 01094 01095 01096 template <typename T> 01097 inline 01098 bool TypeTensor<T>::operator == (const TypeTensor<T>& rhs) const 01099 { 01100 #if LIBMESH_DIM == 1 01101 return (std::abs(_coords[0] - rhs._coords[0]) 01102 < TOLERANCE); 01103 #endif 01104 01105 #if LIBMESH_DIM == 2 01106 return ((std::abs(_coords[0] - rhs._coords[0]) + 01107 std::abs(_coords[1] - rhs._coords[1]) + 01108 std::abs(_coords[2] - rhs._coords[2]) + 01109 std::abs(_coords[3] - rhs._coords[3])) 01110 < 4.*TOLERANCE); 01111 #endif 01112 01113 #if LIBMESH_DIM == 3 01114 return ((std::abs(_coords[0] - rhs._coords[0]) + 01115 std::abs(_coords[1] - rhs._coords[1]) + 01116 std::abs(_coords[2] - rhs._coords[2]) + 01117 std::abs(_coords[3] - rhs._coords[3]) + 01118 std::abs(_coords[4] - rhs._coords[4]) + 01119 std::abs(_coords[5] - rhs._coords[5]) + 01120 std::abs(_coords[6] - rhs._coords[6]) + 01121 std::abs(_coords[7] - rhs._coords[7]) + 01122 std::abs(_coords[8] - rhs._coords[8])) 01123 < 9.*TOLERANCE); 01124 #endif 01125 01126 } 01127 01128 } // namespace libMesh 01129 01130 #endif // LIBMESH_TYPE_TENSOR_H
Site Created By: libMesh Developers
Last modified: February 05 2013 19:54:49 UTC
Hosted By: