$extrastylesheet
type_vector.h
Go to the documentation of this file.
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