$extrastylesheet
00001 // The libMesh Finite Element Library. 00002 // Copyright (C) 2002-2014 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_VECTOR_H 00021 #define LIBMESH_TYPE_VECTOR_H 00022 00023 // Local includes 00024 #include "libmesh/libmesh_common.h" 00025 #include "libmesh/compare_types.h" 00026 #include "libmesh/tensor_tools.h" 00027 00028 // C++ includes 00029 #include <cstdlib> // *must* precede <cmath> for proper std:abs() on PGI, Sun Studio CC 00030 #include <cmath> 00031 #include <complex> 00032 00033 namespace libMesh 00034 { 00035 00036 // Forward declarations 00037 template <typename T> class TypeTensor; 00038 template <typename T> class VectorValue; 00039 template <typename T> class TensorValue; 00040 00052 template <typename T> 00053 class TypeVector 00054 { 00055 template <typename T2> 00056 friend class TypeVector; 00057 00058 friend class TypeTensor<T>; 00059 00060 protected: 00061 00066 TypeVector (); 00067 00068 00073 TypeVector (const T x, 00074 const T y=0, 00075 const T z=0); 00076 00077 00082 template <typename Scalar> 00083 TypeVector (const Scalar x, 00084 const Scalar y=0, 00085 typename 00086 boostcopy::enable_if_c<ScalarTraits<Scalar>::value, 00087 const Scalar>::type z=0); 00088 00089 public: 00090 00094 template <typename T2> 00095 TypeVector (const TypeVector<T2>& p); 00096 00100 ~TypeVector (); 00101 00105 template <typename T2> 00106 void assign (const TypeVector<T2> &); 00107 00111 template <typename Scalar> 00112 typename boostcopy::enable_if_c< 00113 ScalarTraits<Scalar>::value, 00114 TypeVector&>::type 00115 operator = (const Scalar& p) 00116 { libmesh_assert_equal_to (p, Scalar(0)); this->zero(); return *this; } 00117 00121 const T & operator () (const unsigned int i) const; 00122 00123 const T & slice (const unsigned int i) const { return (*this)(i); } 00124 00128 T & operator () (const unsigned int i); 00129 00130 T & slice (const unsigned int i) { return (*this)(i); } 00131 00135 template <typename T2> 00136 TypeVector<typename CompareTypes<T, T2>::supertype> 00137 operator + (const TypeVector<T2> &) const; 00138 00142 template <typename T2> 00143 const TypeVector<T> & operator += (const TypeVector<T2> &); 00144 00148 template <typename T2> 00149 void add (const TypeVector<T2> &); 00150 00155 template <typename T2> 00156 void add_scaled (const TypeVector<T2> &, const T); 00157 00161 template <typename T2> 00162 TypeVector<typename CompareTypes<T, T2>::supertype> 00163 operator - (const TypeVector<T2> &) const; 00164 00168 template <typename T2> 00169 const TypeVector<T> & operator -= (const TypeVector<T2> &); 00170 00174 template <typename T2> 00175 void subtract (const TypeVector<T2> &); 00176 00181 template <typename T2> 00182 void subtract_scaled (const TypeVector<T2> &, const T); 00183 00187 TypeVector<T> operator - () const; 00188 00192 template <typename Scalar> 00193 typename boostcopy::enable_if_c< 00194 ScalarTraits<Scalar>::value, 00195 TypeVector<typename CompareTypes<T, Scalar>::supertype> >::type 00196 operator * (const Scalar) const; 00197 00201 const TypeVector<T> & operator *= (const T); 00202 00206 template <typename Scalar> 00207 typename boostcopy::enable_if_c< 00208 ScalarTraits<Scalar>::value, 00209 TypeVector<typename CompareTypes<T, Scalar>::supertype> >::type 00210 operator / (const Scalar) const; 00211 00215 const TypeVector<T> & operator /= (const T); 00216 00221 template <typename T2> 00222 typename CompareTypes<T, T2>::supertype 00223 operator * (const TypeVector<T2> &) const; 00224 00229 template <typename T2> 00230 typename CompareTypes<T, T2>::supertype 00231 contract (const TypeVector<T2> &) const; 00232 00236 template <typename T2> 00237 TypeVector<typename CompareTypes<T, T2>::supertype> 00238 cross(const TypeVector<T2> &) const; 00239 00244 TypeVector<T> unit() const; 00245 00250 Real size() const; 00251 00256 Real size_sq() const; 00257 00261 void zero(); 00262 00267 bool relative_fuzzy_equals(const TypeVector<T>& rhs, Real tol = TOLERANCE) const; 00268 00273 bool absolute_fuzzy_equals(const TypeVector<T>& rhs, Real tol = TOLERANCE) const; 00274 00279 bool operator == (const TypeVector<T>& rhs) const; 00280 00285 bool operator != (const TypeVector<T>& rhs) const; 00286 00293 bool operator < (const TypeVector<T>& rhs) const; 00294 00301 bool operator <= (const TypeVector<T>& rhs) const; 00302 00309 bool operator > (const TypeVector<T>& rhs) const; 00310 00317 bool operator >= (const TypeVector<T>& rhs) const; 00318 00322 void print(std::ostream& os = libMesh::out) const; 00323 00329 friend std::ostream& operator << (std::ostream& os, const TypeVector<T>& t) 00330 { 00331 t.print(os); 00332 return os; 00333 } 00334 00340 void write_unformatted (std::ostream &out, const bool newline = true) const; 00341 00342 protected: 00343 00347 T _coords[LIBMESH_DIM]; 00348 }; 00349 00350 00351 00352 00353 00354 //------------------------------------------------------ 00355 // Inline functions 00356 00357 template <typename T> 00358 inline 00359 TypeVector<T>::TypeVector () 00360 { 00361 _coords[0] = 0; 00362 00363 #if LIBMESH_DIM > 1 00364 _coords[1] = 0; 00365 #endif 00366 00367 #if LIBMESH_DIM > 2 00368 _coords[2] = 0; 00369 #endif 00370 } 00371 00372 00373 00374 template <typename T> 00375 inline 00376 TypeVector<T>::TypeVector (const T x, 00377 const T y, 00378 const T z) 00379 { 00380 _coords[0] = x; 00381 00382 #if LIBMESH_DIM > 1 00383 _coords[1] = y; 00384 #else 00385 libmesh_assert_equal_to (y, 0); 00386 #endif 00387 00388 #if LIBMESH_DIM > 2 00389 _coords[2] = z; 00390 #else 00391 libmesh_assert_equal_to (z, 0); 00392 #endif 00393 } 00394 00395 00396 template <typename T> 00397 template <typename Scalar> 00398 inline 00399 TypeVector<T>::TypeVector (const Scalar x, 00400 const Scalar y, 00401 typename 00402 boostcopy::enable_if_c<ScalarTraits<Scalar>::value, 00403 const Scalar>::type z) 00404 { 00405 _coords[0] = x; 00406 00407 #if LIBMESH_DIM > 1 00408 _coords[1] = y; 00409 #else 00410 libmesh_assert_equal_to (y, 0); 00411 #endif 00412 00413 #if LIBMESH_DIM > 2 00414 _coords[2] = z; 00415 #else 00416 libmesh_assert_equal_to (z, 0); 00417 #endif 00418 } 00419 00420 00421 00422 template <typename T> 00423 template <typename T2> 00424 inline 00425 TypeVector<T>::TypeVector (const TypeVector<T2> &p) 00426 { 00427 // copy the nodes from vector p to me 00428 for (unsigned int i=0; i<LIBMESH_DIM; i++) 00429 _coords[i] = p._coords[i]; 00430 } 00431 00432 00433 00434 template <typename T> 00435 inline 00436 TypeVector<T>::~TypeVector () 00437 { 00438 } 00439 00440 00441 00442 template <typename T> 00443 template <typename T2> 00444 inline 00445 void TypeVector<T>::assign (const TypeVector<T2> &p) 00446 { 00447 for (unsigned int i=0; i<LIBMESH_DIM; i++) 00448 _coords[i] = p._coords[i]; 00449 } 00450 00451 00452 00453 template <typename T> 00454 inline 00455 const T & TypeVector<T>::operator () (const unsigned int i) const 00456 { 00457 libmesh_assert_less (i, LIBMESH_DIM); 00458 00459 return _coords[i]; 00460 } 00461 00462 00463 00464 template <typename T> 00465 inline 00466 T & TypeVector<T>::operator () (const unsigned int i) 00467 { 00468 libmesh_assert_less (i, LIBMESH_DIM); 00469 00470 return _coords[i]; 00471 } 00472 00473 00474 00475 template <typename T> 00476 template <typename T2> 00477 inline 00478 TypeVector<typename CompareTypes<T, T2>::supertype> 00479 TypeVector<T>::operator + (const TypeVector<T2> &p) const 00480 { 00481 typedef typename CompareTypes<T, T2>::supertype TS; 00482 #if LIBMESH_DIM == 1 00483 return TypeVector<TS> (_coords[0] + p._coords[0]); 00484 #endif 00485 00486 #if LIBMESH_DIM == 2 00487 return TypeVector<TS> (_coords[0] + p._coords[0], 00488 _coords[1] + p._coords[1]); 00489 #endif 00490 00491 #if LIBMESH_DIM == 3 00492 return TypeVector<TS> (_coords[0] + p._coords[0], 00493 _coords[1] + p._coords[1], 00494 _coords[2] + p._coords[2]); 00495 #endif 00496 00497 } 00498 00499 00500 00501 template <typename T> 00502 template <typename T2> 00503 inline 00504 const TypeVector<T> & TypeVector<T>::operator += (const TypeVector<T2> &p) 00505 { 00506 this->add (p); 00507 00508 return *this; 00509 } 00510 00511 00512 00513 template <typename T> 00514 template <typename T2> 00515 inline 00516 void TypeVector<T>::add (const TypeVector<T2> &p) 00517 { 00518 #if LIBMESH_DIM == 1 00519 _coords[0] += p._coords[0]; 00520 #endif 00521 00522 #if LIBMESH_DIM == 2 00523 _coords[0] += p._coords[0]; 00524 _coords[1] += p._coords[1]; 00525 #endif 00526 00527 #if LIBMESH_DIM == 3 00528 _coords[0] += p._coords[0]; 00529 _coords[1] += p._coords[1]; 00530 _coords[2] += p._coords[2]; 00531 #endif 00532 00533 } 00534 00535 00536 00537 template <typename T> 00538 template <typename T2> 00539 inline 00540 void TypeVector<T>::add_scaled (const TypeVector<T2> &p, const T factor) 00541 { 00542 #if LIBMESH_DIM == 1 00543 _coords[0] += factor*p(0); 00544 #endif 00545 00546 #if LIBMESH_DIM == 2 00547 _coords[0] += factor*p(0); 00548 _coords[1] += factor*p(1); 00549 #endif 00550 00551 #if LIBMESH_DIM == 3 00552 _coords[0] += factor*p(0); 00553 _coords[1] += factor*p(1); 00554 _coords[2] += factor*p(2); 00555 #endif 00556 00557 } 00558 00559 00560 00561 template <typename T> 00562 template <typename T2> 00563 inline 00564 TypeVector<typename CompareTypes<T, T2>::supertype> 00565 TypeVector<T>::operator - (const TypeVector<T2> &p) const 00566 { 00567 typedef typename CompareTypes<T, T2>::supertype TS; 00568 00569 #if LIBMESH_DIM == 1 00570 return TypeVector<TS>(_coords[0] - p._coords[0]); 00571 #endif 00572 00573 #if LIBMESH_DIM == 2 00574 return TypeVector<TS>(_coords[0] - p._coords[0], 00575 _coords[1] - p._coords[1]); 00576 #endif 00577 00578 #if LIBMESH_DIM == 3 00579 return TypeVector<TS>(_coords[0] - p._coords[0], 00580 _coords[1] - p._coords[1], 00581 _coords[2] - p._coords[2]); 00582 #endif 00583 00584 } 00585 00586 00587 00588 template <typename T> 00589 template <typename T2> 00590 inline 00591 const TypeVector<T> & TypeVector<T>::operator -= (const TypeVector<T2> &p) 00592 { 00593 this->subtract (p); 00594 00595 return *this; 00596 } 00597 00598 00599 00600 template <typename T> 00601 template <typename T2> 00602 inline 00603 void TypeVector<T>::subtract (const TypeVector<T2>& p) 00604 { 00605 for (unsigned int i=0; i<LIBMESH_DIM; i++) 00606 _coords[i] -= p._coords[i]; 00607 } 00608 00609 00610 00611 template <typename T> 00612 template <typename T2> 00613 inline 00614 void TypeVector<T>::subtract_scaled (const TypeVector<T2> &p, const T factor) 00615 { 00616 for (unsigned int i=0; i<LIBMESH_DIM; i++) 00617 _coords[i] -= factor*p(i); 00618 } 00619 00620 00621 00622 template <typename T> 00623 inline 00624 TypeVector<T> TypeVector<T>::operator - () const 00625 { 00626 00627 #if LIBMESH_DIM == 1 00628 return TypeVector(-_coords[0]); 00629 #endif 00630 00631 #if LIBMESH_DIM == 2 00632 return TypeVector(-_coords[0], 00633 -_coords[1]); 00634 #endif 00635 00636 #if LIBMESH_DIM == 3 00637 return TypeVector(-_coords[0], 00638 -_coords[1], 00639 -_coords[2]); 00640 #endif 00641 00642 } 00643 00644 00645 00646 template <typename T> 00647 template <typename Scalar> 00648 inline 00649 typename boostcopy::enable_if_c< 00650 ScalarTraits<Scalar>::value, 00651 TypeVector<typename CompareTypes<T, Scalar>::supertype> >::type 00652 TypeVector<T>::operator * (const Scalar factor) const 00653 { 00654 typedef typename CompareTypes<T, Scalar>::supertype SuperType; 00655 00656 #if LIBMESH_DIM == 1 00657 return TypeVector<SuperType>(_coords[0]*factor); 00658 #endif 00659 00660 #if LIBMESH_DIM == 2 00661 return TypeVector<SuperType>(_coords[0]*factor, 00662 _coords[1]*factor); 00663 #endif 00664 00665 #if LIBMESH_DIM == 3 00666 return TypeVector<SuperType>(_coords[0]*factor, 00667 _coords[1]*factor, 00668 _coords[2]*factor); 00669 #endif 00670 } 00671 00672 00673 00674 template <typename T, typename Scalar> 00675 inline 00676 typename boostcopy::enable_if_c< 00677 ScalarTraits<Scalar>::value, 00678 TypeVector<typename CompareTypes<T, Scalar>::supertype> >::type 00679 operator * (const Scalar factor, 00680 const TypeVector<T> &v) 00681 { 00682 return v * factor; 00683 } 00684 00685 00686 00687 template <typename T> 00688 inline 00689 const TypeVector<T> & TypeVector<T>::operator *= (const T factor) 00690 { 00691 #if LIBMESH_DIM == 1 00692 _coords[0] *= factor; 00693 #endif 00694 00695 #if LIBMESH_DIM == 2 00696 _coords[0] *= factor; 00697 _coords[1] *= factor; 00698 #endif 00699 00700 #if LIBMESH_DIM == 3 00701 _coords[0] *= factor; 00702 _coords[1] *= factor; 00703 _coords[2] *= factor; 00704 #endif 00705 00706 return *this; 00707 } 00708 00709 00710 00711 template <typename T> 00712 template <typename Scalar> 00713 inline 00714 typename boostcopy::enable_if_c< 00715 ScalarTraits<Scalar>::value, 00716 TypeVector<typename CompareTypes<T, Scalar>::supertype> >::type 00717 TypeVector<T>::operator / (const Scalar factor) const 00718 { 00719 libmesh_assert_not_equal_to (factor, static_cast<T>(0.)); 00720 00721 typedef typename CompareTypes<T, Scalar>::supertype TS; 00722 00723 #if LIBMESH_DIM == 1 00724 return TypeVector<TS>(_coords[0]/factor); 00725 #endif 00726 00727 #if LIBMESH_DIM == 2 00728 return TypeVector<TS>(_coords[0]/factor, 00729 _coords[1]/factor); 00730 #endif 00731 00732 #if LIBMESH_DIM == 3 00733 return TypeVector<TS>(_coords[0]/factor, 00734 _coords[1]/factor, 00735 _coords[2]/factor); 00736 #endif 00737 00738 } 00739 00740 00741 00742 00743 template <typename T> 00744 inline 00745 const TypeVector<T> & 00746 TypeVector<T>::operator /= (const T factor) 00747 { 00748 libmesh_assert_not_equal_to (factor, static_cast<T>(0.)); 00749 00750 for (unsigned int i=0; i<LIBMESH_DIM; i++) 00751 _coords[i] /= factor; 00752 00753 return *this; 00754 } 00755 00756 00757 00758 00759 template <typename T> 00760 template <typename T2> 00761 inline 00762 typename CompareTypes<T, T2>::supertype 00763 TypeVector<T>::operator * (const TypeVector<T2> &p) const 00764 { 00765 #if LIBMESH_DIM == 1 00766 return _coords[0]*p._coords[0]; 00767 #endif 00768 00769 #if LIBMESH_DIM == 2 00770 return (_coords[0]*p._coords[0] + 00771 _coords[1]*p._coords[1]); 00772 #endif 00773 00774 #if LIBMESH_DIM == 3 00775 return (_coords[0]*p(0) + 00776 _coords[1]*p(1) + 00777 _coords[2]*p(2)); 00778 #endif 00779 } 00780 00781 template <typename T> 00782 template <typename T2> 00783 inline 00784 typename CompareTypes<T, T2>::supertype 00785 TypeVector<T>::contract(const TypeVector<T2> &p) const 00786 { 00787 return (*this)*(p); 00788 } 00789 00790 00791 00792 template <typename T> 00793 template <typename T2> 00794 TypeVector<typename CompareTypes<T, T2>::supertype> 00795 TypeVector<T>::cross(const TypeVector<T2>& p) const 00796 { 00797 typedef typename CompareTypes<T, T2>::supertype TS; 00798 libmesh_assert_equal_to (LIBMESH_DIM, 3); 00799 00800 // | i j k | 00801 // |(*this)(0) (*this)(1) (*this)(2)| 00802 // | p(0) p(1) p(2) | 00803 00804 return TypeVector<TS>( _coords[1]*p._coords[2] - _coords[2]*p._coords[1], 00805 -_coords[0]*p._coords[2] + _coords[2]*p._coords[0], 00806 _coords[0]*p._coords[1] - _coords[1]*p._coords[0]); 00807 } 00808 00809 00810 00811 template <typename T> 00812 inline 00813 Real TypeVector<T>::size() const 00814 { 00815 return std::sqrt(this->size_sq()); 00816 } 00817 00818 00819 00820 template <typename T> 00821 inline 00822 void TypeVector<T>::zero() 00823 { 00824 for (unsigned int i=0; i<LIBMESH_DIM; i++) 00825 _coords[i] = 0.; 00826 } 00827 00828 00829 00830 template <typename T> 00831 inline 00832 Real TypeVector<T>::size_sq() const 00833 { 00834 #if LIBMESH_DIM == 1 00835 return (TensorTools::norm_sq(_coords[0])); 00836 #endif 00837 00838 #if LIBMESH_DIM == 2 00839 return (TensorTools::norm_sq(_coords[0]) + 00840 TensorTools::norm_sq(_coords[1])); 00841 #endif 00842 00843 #if LIBMESH_DIM == 3 00844 return (TensorTools::norm_sq(_coords[0]) + 00845 TensorTools::norm_sq(_coords[1]) + 00846 TensorTools::norm_sq(_coords[2])); 00847 #endif 00848 } 00849 00850 00851 00852 template <typename T> 00853 inline 00854 bool TypeVector<T>::absolute_fuzzy_equals(const TypeVector<T>& rhs, Real tol) const 00855 { 00856 #if LIBMESH_DIM == 1 00857 return (std::abs(_coords[0] - rhs._coords[0]) 00858 <= tol); 00859 #endif 00860 00861 #if LIBMESH_DIM == 2 00862 return (std::abs(_coords[0] - rhs._coords[0]) + 00863 std::abs(_coords[1] - rhs._coords[1]) 00864 <= tol); 00865 #endif 00866 00867 #if LIBMESH_DIM == 3 00868 return (std::abs(_coords[0] - rhs._coords[0]) + 00869 std::abs(_coords[1] - rhs._coords[1]) + 00870 std::abs(_coords[2] - rhs._coords[2]) 00871 <= tol); 00872 #endif 00873 } 00874 00875 00876 00877 template <typename T> 00878 inline 00879 bool TypeVector<T>::relative_fuzzy_equals(const TypeVector<T>& rhs, Real tol) const 00880 { 00881 #if LIBMESH_DIM == 1 00882 return this->absolute_fuzzy_equals(rhs, tol * 00883 (std::abs(_coords[0]) + std::abs(rhs._coords[0]))); 00884 #endif 00885 00886 #if LIBMESH_DIM == 2 00887 return this->absolute_fuzzy_equals(rhs, tol * 00888 (std::abs(_coords[0]) + std::abs(rhs._coords[0]) + 00889 std::abs(_coords[1]) + std::abs(rhs._coords[1]))); 00890 #endif 00891 00892 #if LIBMESH_DIM == 3 00893 return this->absolute_fuzzy_equals(rhs, tol * 00894 (std::abs(_coords[0]) + std::abs(rhs._coords[0]) + 00895 std::abs(_coords[1]) + std::abs(rhs._coords[1]) + 00896 std::abs(_coords[2]) + std::abs(rhs._coords[2]))); 00897 #endif 00898 } 00899 00900 00901 00902 template <typename T> 00903 inline 00904 bool TypeVector<T>::operator == (const TypeVector<T>& rhs) const 00905 { 00906 #if LIBMESH_DIM == 1 00907 return (_coords[0] == rhs._coords[0]); 00908 #endif 00909 00910 #if LIBMESH_DIM == 2 00911 return (_coords[0] == rhs._coords[0] && 00912 _coords[1] == rhs._coords[1]); 00913 #endif 00914 00915 #if LIBMESH_DIM == 3 00916 return (_coords[0] == rhs._coords[0] && 00917 _coords[1] == rhs._coords[1] && 00918 _coords[2] == rhs._coords[2]); 00919 #endif 00920 } 00921 00922 00923 00924 template <> 00925 inline 00926 bool TypeVector<Real>::operator != (const TypeVector<Real>& rhs) const 00927 { 00928 return (!(*this == rhs)); 00929 } 00930 00931 } // namespace libMesh 00932 00933 #endif // LIBMESH_TYPE_VECTOR_H