libMesh
type_tensor.h
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2017 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 
4 // This library is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU Lesser General Public
6 // License as published by the Free Software Foundation; either
7 // version 2.1 of the License, or (at your option) any later version.
8 
9 // This library is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 // Lesser General Public License for more details.
13 
14 // You should have received a copy of the GNU Lesser General Public
15 // License along with this library; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 
18 
19 
20 #ifndef LIBMESH_TYPE_TENSOR_H
21 #define LIBMESH_TYPE_TENSOR_H
22 
23 // Local includes
24 #include "libmesh/libmesh_common.h"
25 #include "libmesh/type_vector.h"
26 
27 // C++ includes
28 #include <cstdlib> // *must* precede <cmath> for proper std:abs() on PGI, Sun Studio CC
29 #include <cmath>
30 
31 namespace libMesh
32 {
33 
34 // Forward declarations
35 template <typename T> class TypeTensorColumn;
36 template <typename T> class ConstTypeTensorColumn;
37 template <unsigned int N, typename T> class TypeNTensor;
38 
48 template <typename T>
49 class TypeTensor
50 {
51  template <typename T2>
52  friend class TypeTensor;
53 
54 protected:
55 
60  TypeTensor ();
61 
68  explicit TypeTensor (const T xx,
69  const T xy=0,
70  const T xz=0,
71  const T yx=0,
72  const T yy=0,
73  const T yz=0,
74  const T zx=0,
75  const T zy=0,
76  const T zz=0);
77 
78 
82  template <typename Scalar>
83  explicit TypeTensor (const Scalar xx,
84  const Scalar xy=0,
85  const Scalar xz=0,
86  const Scalar yx=0,
87  const Scalar yy=0,
88  const Scalar yz=0,
89  const Scalar zx=0,
90  const Scalar zy=0,
91  typename
93  const Scalar>::type zz=0);
94 
100  template<typename T2>
101  TypeTensor(const TypeVector<T2> & vx);
102 
103  template<typename T2>
104  TypeTensor(const TypeVector<T2> & vx,
105  const TypeVector<T2> & vy);
106 
107  template<typename T2>
108  TypeTensor(const TypeVector<T2> & vx,
109  const TypeVector<T2> & vy,
110  const TypeVector<T2> & vz);
111 
112 public:
113 
117  typedef T value_type;
118 
122  template<typename T2>
123  TypeTensor(const TypeTensor<T2> & p);
124 
128  ~TypeTensor();
129 
133  template<typename T2>
134  void assign (const TypeTensor<T2> &);
135 
141  template <typename Scalar>
142  typename boostcopy::enable_if_c<
144  TypeTensor &>::type
145  operator = (const Scalar & libmesh_dbg_var(p))
146  { libmesh_assert_equal_to (p, Scalar(0)); this->zero(); return *this; }
147 
151  const T & operator () (const unsigned int i, const unsigned int j) const;
152 
156  T & operator () (const unsigned int i, const unsigned int j);
157 
161  ConstTypeTensorColumn<T> slice (const unsigned int i) const;
162 
166  TypeTensorColumn<T> slice (const unsigned int i);
167 
171  TypeVector<T> row(const unsigned int r) const;
172 
178  template<typename T2>
180  operator + (const TypeTensor<T2> &) const;
181 
187  template<typename T2>
188  const TypeTensor<T> & operator += (const TypeTensor<T2> &);
189 
193  template<typename T2>
194  void add (const TypeTensor<T2> &);
195 
199  template <typename T2>
200  void add_scaled (const TypeTensor<T2> &, const T);
201 
207  template<typename T2>
209  operator - (const TypeTensor<T2> &) const;
210 
216  template<typename T2>
217  const TypeTensor<T> & operator -= (const TypeTensor<T2> &);
218 
222  template<typename T2>
223  void subtract (const TypeTensor<T2> &);
224 
229  template <typename T2>
230  void subtract_scaled (const TypeTensor<T2> &, const T);
231 
235  TypeTensor<T> operator - () const;
236 
242  template <typename Scalar>
243  typename boostcopy::enable_if_c<
244  ScalarTraits<Scalar>::value,
246  operator * (const Scalar) const;
247 
253  template <typename Scalar>
254  const TypeTensor<T> & operator *= (const Scalar);
255 
261  template <typename Scalar>
262  typename boostcopy::enable_if_c<
263  ScalarTraits<Scalar>::value,
265  operator / (const Scalar) const;
266 
272  const TypeTensor<T> & operator /= (const T);
273 
280  template <typename T2>
282 
290  template <typename T2>
292  contract (const TypeTensor<T2> &) const;
293 
300  template <typename T2>
302  operator * (const TypeVector<T2> &) const;
303 
307  TypeTensor<T> transpose() const;
308 
312  TypeTensor<T> inverse() const;
313 
320  void solve(const TypeVector<T> & b, TypeVector<T> & x) const;
321 
328 #ifdef LIBMESH_ENABLE_DEPRECATED
329  Real size() const;
330 #endif
331 
336  Real norm() const;
337 
344 #ifdef LIBMESH_ENABLE_DEPRECATED
345  Real size_sq() const;
346 #endif
347 
352  Real norm_sq() const;
353 
360  T det() const;
361 
365  T tr() const;
366 
370  void zero();
371 
375  bool operator == (const TypeTensor<T> & rhs) const;
376 
381  bool operator < (const TypeTensor<T> & rhs) const;
382 
386  bool operator > (const TypeTensor<T> & rhs) const;
387 
391  void print(std::ostream & os = libMesh::out) const;
392 
400  friend std::ostream & operator << (std::ostream & os, const TypeTensor<T> & t)
401  {
402  t.print(os);
403  return os;
404  }
405 
410  void write_unformatted (std::ostream & out, const bool newline = true) const;
411 
412 protected:
413 
417  T _coords[LIBMESH_DIM*LIBMESH_DIM];
418 };
419 
420 
421 template <typename T>
422 class TypeTensorColumn
423 {
424 public:
425 
427  unsigned int j) :
428  _tensor(&tensor), _j(j) {}
429 
434  T & operator () (const unsigned int i)
435  { return (*_tensor)(i,_j); }
436 
437  T & slice (const unsigned int i)
438  { return (*_tensor)(i,_j); }
439 
444  {
445  for (unsigned int i=0; i != LIBMESH_DIM; ++i)
446  (*this)(i) = rhs(i);
447  return *this;
448  }
449 
450 private:
452  const unsigned int _j;
453 };
454 
455 
456 template <typename T>
458 {
459 public:
460 
462  unsigned int j) :
463  _tensor(&tensor), _j(j) {}
464 
468  const T & operator () (const unsigned int i) const
469  { return (*_tensor)(i,_j); }
470 
471  const T & slice (const unsigned int i) const
472  { return (*_tensor)(i,_j); }
473 
474 private:
476  const unsigned int _j;
477 };
478 
479 //------------------------------------------------------
480 // Inline functions
481 template <typename T>
482 inline
484 {
485  _coords[0] = 0;
486 
487 #if LIBMESH_DIM > 1
488  _coords[1] = 0;
489  _coords[2] = 0;
490  _coords[3] = 0;
491 #endif
492 
493 #if LIBMESH_DIM > 2
494  _coords[4] = 0;
495  _coords[5] = 0;
496  _coords[6] = 0;
497  _coords[7] = 0;
498  _coords[8] = 0;
499 #endif
500 }
501 
502 
503 
504 template <typename T>
505 inline
507  const T xy,
508  const T xz,
509  const T yx,
510  const T yy,
511  const T yz,
512  const T zx,
513  const T zy,
514  const T zz)
515 {
516  _coords[0] = xx;
517 
518 #if LIBMESH_DIM == 2
519  _coords[1] = xy;
520  _coords[2] = yx;
521  _coords[3] = yy;
522 #endif
523 
524 #if LIBMESH_DIM == 3
525  _coords[1] = xy;
526  _coords[2] = xz;
527  _coords[3] = yx;
528  _coords[4] = yy;
529  _coords[5] = yz;
530  _coords[6] = zx;
531  _coords[7] = zy;
532  _coords[8] = zz;
533 #endif
534 }
535 
536 
537 template <typename T>
538 template <typename Scalar>
539 inline
540 TypeTensor<T>::TypeTensor (const Scalar xx,
541  const Scalar xy,
542  const Scalar xz,
543  const Scalar yx,
544  const Scalar yy,
545  const Scalar yz,
546  const Scalar zx,
547  const Scalar zy,
548  typename
549  boostcopy::enable_if_c<ScalarTraits<Scalar>::value,
550  const Scalar>::type zz)
551 {
552  _coords[0] = xx;
553 
554 #if LIBMESH_DIM == 2
555  _coords[1] = xy;
556  _coords[2] = yx;
557  _coords[3] = yy;
558 #endif
559 
560 #if LIBMESH_DIM == 3
561  _coords[1] = xy;
562  _coords[2] = xz;
563  _coords[3] = yx;
564  _coords[4] = yy;
565  _coords[5] = yz;
566  _coords[6] = zx;
567  _coords[7] = zy;
568  _coords[8] = zz;
569 #endif
570 }
571 
572 
573 
574 template <typename T>
575 template<typename T2>
576 inline
578 {
579  // copy the nodes from vector p to me
580  for (unsigned int i=0; i<LIBMESH_DIM*LIBMESH_DIM; i++)
581  _coords[i] = p._coords[i];
582 }
583 
584 
585 template <typename T>
586 template <typename T2>
588 {
589  libmesh_assert_equal_to (LIBMESH_DIM, 1);
590  _coords[0] = vx(0);
591 }
592 
593 template <typename T>
594 template <typename T2>
596  const TypeVector<T2> & vy)
597 {
598  libmesh_assert_equal_to (LIBMESH_DIM, 2);
599  _coords[0] = vx(0);
600  _coords[1] = vx(1);
601  _coords[2] = vy(0);
602  _coords[3] = vy(1);
603 }
604 
605 template <typename T>
606 template <typename T2>
608  const TypeVector<T2> & vy,
609  const TypeVector<T2> & vz)
610 {
611  libmesh_assert_equal_to (LIBMESH_DIM, 3);
612  _coords[0] = vx(0);
613  _coords[1] = vx(1);
614  _coords[2] = vx(2);
615  _coords[3] = vy(0);
616  _coords[4] = vy(1);
617  _coords[5] = vy(2);
618  _coords[6] = vz(0);
619  _coords[7] = vz(1);
620  _coords[8] = vz(2);
621 }
622 
623 
624 
625 
626 template <typename T>
627 inline
629 {
630 }
631 
632 
633 
634 template <typename T>
635 template<typename T2>
636 inline
638 {
639  for (unsigned int i=0; i<LIBMESH_DIM*LIBMESH_DIM; i++)
640  _coords[i] = p._coords[i];
641 }
642 
643 
644 
645 template <typename T>
646 inline
647 const T & TypeTensor<T>::operator () (const unsigned int i,
648  const unsigned int j) const
649 {
650  libmesh_assert_less (i, 3);
651  libmesh_assert_less (j, 3);
652 
653 #if LIBMESH_DIM < 3
654  const static T my_zero = 0;
655  if (i >= LIBMESH_DIM || j >= LIBMESH_DIM)
656  return my_zero;
657 #endif
658 
659  return _coords[i*LIBMESH_DIM+j];
660 }
661 
662 
663 
664 template <typename T>
665 inline
666 T & TypeTensor<T>::operator () (const unsigned int i,
667  const unsigned int j)
668 {
669 #if LIBMESH_DIM < 3
670 
671  if (i >= LIBMESH_DIM || j >= LIBMESH_DIM)
672  libmesh_error_msg("ERROR: You are assigning to a tensor component that is out of range for the compiled LIBMESH_DIM!");
673 
674 #endif
675 
676  libmesh_assert_less (i, LIBMESH_DIM);
677  libmesh_assert_less (j, LIBMESH_DIM);
678 
679  return _coords[i*LIBMESH_DIM+j];
680 }
681 
682 
683 template <typename T>
684 inline
686 TypeTensor<T>::slice (const unsigned int i) const
687 {
688  libmesh_assert_less (i, LIBMESH_DIM);
689  return ConstTypeTensorColumn<T>(*this, i);
690 }
691 
692 
693 template <typename T>
694 inline
696 TypeTensor<T>::slice (const unsigned int i)
697 {
698  libmesh_assert_less (i, LIBMESH_DIM);
699  return TypeTensorColumn<T>(*this, i);
700 }
701 
702 
703 template <typename T>
704 inline
706 TypeTensor<T>::row(const unsigned int r) const
707 {
708  TypeVector<T> return_vector;
709 
710  for (unsigned int j=0; j<LIBMESH_DIM; j++)
711  return_vector._coords[j] = _coords[r*LIBMESH_DIM + j];
712 
713  return return_vector;
714 }
715 
716 template <typename T>
717 template<typename T2>
718 inline
721 {
722 
723 #if LIBMESH_DIM == 1
724  return TypeTensor(_coords[0] + p._coords[0]);
725 #endif
726 
727 #if LIBMESH_DIM == 2
728  return TypeTensor(_coords[0] + p._coords[0],
729  _coords[1] + p._coords[1],
730  0.,
731  _coords[2] + p._coords[2],
732  _coords[3] + p._coords[3]);
733 #endif
734 
735 #if LIBMESH_DIM == 3
736  return TypeTensor(_coords[0] + p._coords[0],
737  _coords[1] + p._coords[1],
738  _coords[2] + p._coords[2],
739  _coords[3] + p._coords[3],
740  _coords[4] + p._coords[4],
741  _coords[5] + p._coords[5],
742  _coords[6] + p._coords[6],
743  _coords[7] + p._coords[7],
744  _coords[8] + p._coords[8]);
745 #endif
746 
747 }
748 
749 
750 
751 template <typename T>
752 template<typename T2>
753 inline
755 {
756  this->add (p);
757 
758  return *this;
759 }
760 
761 
762 
763 template <typename T>
764 template<typename T2>
765 inline
767 {
768  for (unsigned int i=0; i<LIBMESH_DIM*LIBMESH_DIM; i++)
769  _coords[i] += p._coords[i];
770 }
771 
772 
773 
774 template <typename T>
775 template <typename T2>
776 inline
777 void TypeTensor<T>::add_scaled (const TypeTensor<T2> & p, const T factor)
778 {
779  for (unsigned int i=0; i<LIBMESH_DIM*LIBMESH_DIM; i++)
780  _coords[i] += factor*p._coords[i];
781 
782 }
783 
784 
785 
786 template <typename T>
787 template<typename T2>
788 inline
791 {
792 
793 #if LIBMESH_DIM == 1
794  return TypeTensor(_coords[0] - p._coords[0]);
795 #endif
796 
797 #if LIBMESH_DIM == 2
798  return TypeTensor(_coords[0] - p._coords[0],
799  _coords[1] - p._coords[1],
800  0.,
801  _coords[2] - p._coords[2],
802  _coords[3] - p._coords[3]);
803 #endif
804 
805 #if LIBMESH_DIM == 3
806  return TypeTensor(_coords[0] - p._coords[0],
807  _coords[1] - p._coords[1],
808  _coords[2] - p._coords[2],
809  _coords[3] - p._coords[3],
810  _coords[4] - p._coords[4],
811  _coords[5] - p._coords[5],
812  _coords[6] - p._coords[6],
813  _coords[7] - p._coords[7],
814  _coords[8] - p._coords[8]);
815 #endif
816 
817 }
818 
819 
820 
821 template <typename T>
822 template <typename T2>
823 inline
825 {
826  this->subtract (p);
827 
828  return *this;
829 }
830 
831 
832 
833 template <typename T>
834 template <typename T2>
835 inline
837 {
838  for (unsigned int i=0; i<LIBMESH_DIM*LIBMESH_DIM; i++)
839  _coords[i] -= p._coords[i];
840 }
841 
842 
843 
844 template <typename T>
845 template <typename T2>
846 inline
847 void TypeTensor<T>::subtract_scaled (const TypeTensor<T2> & p, const T factor)
848 {
849  for (unsigned int i=0; i<LIBMESH_DIM*LIBMESH_DIM; i++)
850  _coords[i] -= factor*p._coords[i];
851 }
852 
853 
854 
855 template <typename T>
856 inline
858 {
859 
860 #if LIBMESH_DIM == 1
861  return TypeTensor(-_coords[0]);
862 #endif
863 
864 #if LIBMESH_DIM == 2
865  return TypeTensor(-_coords[0],
866  -_coords[1],
867  -_coords[2],
868  -_coords[3]);
869 #endif
870 
871 #if LIBMESH_DIM == 3
872  return TypeTensor(-_coords[0],
873  -_coords[1],
874  -_coords[2],
875  -_coords[3],
876  -_coords[4],
877  -_coords[5],
878  -_coords[6],
879  -_coords[7],
880  -_coords[8]);
881 #endif
882 
883 }
884 
885 
886 
887 template <typename T>
888 template <typename Scalar>
889 inline
890 typename boostcopy::enable_if_c<
891  ScalarTraits<Scalar>::value,
893 TypeTensor<T>::operator * (const Scalar factor) const
894 {
895  typedef typename CompareTypes<T, Scalar>::supertype TS;
896 
897 
898 #if LIBMESH_DIM == 1
899  return TypeTensor<TS>(_coords[0]*factor);
900 #endif
901 
902 #if LIBMESH_DIM == 2
903  return TypeTensor<TS>(_coords[0]*factor,
904  _coords[1]*factor,
905  _coords[2]*factor,
906  _coords[3]*factor);
907 #endif
908 
909 #if LIBMESH_DIM == 3
910  return TypeTensor<TS>(_coords[0]*factor,
911  _coords[1]*factor,
912  _coords[2]*factor,
913  _coords[3]*factor,
914  _coords[4]*factor,
915  _coords[5]*factor,
916  _coords[6]*factor,
917  _coords[7]*factor,
918  _coords[8]*factor);
919 #endif
920 }
921 
922 
923 
924 template <typename T, typename Scalar>
925 inline
926 typename boostcopy::enable_if_c<
927  ScalarTraits<Scalar>::value,
929 operator * (const Scalar factor,
930  const TypeTensor<T> & t)
931 {
932  return t * factor;
933 }
934 
935 
936 
937 template <typename T>
938 template <typename Scalar>
939 inline
940 const TypeTensor<T> & TypeTensor<T>::operator *= (const Scalar factor)
941 {
942  for (unsigned int i=0; i<LIBMESH_DIM*LIBMESH_DIM; i++)
943  _coords[i] *= factor;
944 
945  return *this;
946 }
947 
948 
949 
950 
951 template <typename T>
952 template <typename Scalar>
953 inline
954 typename boostcopy::enable_if_c<
955  ScalarTraits<Scalar>::value,
957 TypeTensor<T>::operator / (const Scalar factor) const
958 {
959  libmesh_assert_not_equal_to (factor, static_cast<T>(0.));
960 
961  typedef typename CompareTypes<T, Scalar>::supertype TS;
962 
963 #if LIBMESH_DIM == 1
964  return TypeTensor<TS>(_coords[0]/factor);
965 #endif
966 
967 #if LIBMESH_DIM == 2
968  return TypeTensor<TS>(_coords[0]/factor,
969  _coords[1]/factor,
970  _coords[2]/factor,
971  _coords[3]/factor);
972 #endif
973 
974 #if LIBMESH_DIM == 3
975  return TypeTensor<TS>(_coords[0]/factor,
976  _coords[1]/factor,
977  _coords[2]/factor,
978  _coords[3]/factor,
979  _coords[4]/factor,
980  _coords[5]/factor,
981  _coords[6]/factor,
982  _coords[7]/factor,
983  _coords[8]/factor);
984 #endif
985 
986 }
987 
988 
989 
990 template <typename T>
991 inline
993 {
994 #if LIBMESH_DIM == 1
995  return TypeTensor(_coords[0]);
996 #endif
997 
998 #if LIBMESH_DIM == 2
999  return TypeTensor(_coords[0],
1000  _coords[2],
1001  _coords[1],
1002  _coords[3]);
1003 #endif
1004 
1005 #if LIBMESH_DIM == 3
1006  return TypeTensor(_coords[0],
1007  _coords[3],
1008  _coords[6],
1009  _coords[1],
1010  _coords[4],
1011  _coords[7],
1012  _coords[2],
1013  _coords[5],
1014  _coords[8]);
1015 #endif
1016 }
1017 
1018 
1019 
1020 template <typename T>
1021 inline
1023 {
1024 #if LIBMESH_DIM == 1
1025  if (_coords[0] == static_cast<T>(0.))
1026  libmesh_convergence_failure();
1027  return TypeTensor(1. / _coords[0]);
1028 #endif
1029 
1030 #if LIBMESH_DIM == 2
1031  // Get convenient reference to this.
1032  const TypeTensor<T> & A = *this;
1033 
1034  // Use temporary variables, avoid multiple accesses
1035  T a = A(0,0), b = A(0,1),
1036  c = A(1,0), d = A(1,1);
1037 
1038  // Make sure det = ad - bc is not zero
1039  T my_det = a*d - b*c;
1040 
1041  if (my_det == static_cast<T>(0.))
1042  libmesh_convergence_failure();
1043 
1044  return TypeTensor(d/my_det, -b/my_det, -c/my_det, a/my_det);
1045 #endif
1046 
1047 #if LIBMESH_DIM == 3
1048  // Get convenient reference to this.
1049  const TypeTensor<T> & A = *this;
1050 
1051  T a11 = A(0,0), a12 = A(0,1), a13 = A(0,2),
1052  a21 = A(1,0), a22 = A(1,1), a23 = A(1,2),
1053  a31 = A(2,0), a32 = A(2,1), a33 = A(2,2);
1054 
1055  T my_det = a11*(a33*a22-a32*a23) - a21*(a33*a12-a32*a13) + a31*(a23*a12-a22*a13);
1056 
1057  if (my_det == static_cast<T>(0.))
1058  libmesh_convergence_failure();
1059 
1060  // Inline comment characters are for lining up columns.
1061  return TypeTensor( (a33*a22-a32*a23)/my_det, -(a33*a12-a32*a13)/my_det, (a23*a12-a22*a13)/my_det,
1062  -(a33*a21-a31*a23)/my_det, (a33*a11-a31*a13)/my_det, -(a23*a11-a21*a13)/my_det,
1063  (a32*a21-a31*a22)/my_det, -(a32*a11-a31*a12)/my_det, (a22*a11-a21*a12)/my_det);
1064 #endif
1065 }
1066 
1067 
1068 
1069 template <typename T>
1070 inline
1072 {
1073 #if LIBMESH_DIM == 1
1074  if (_coords[0] == static_cast<T>(0.))
1075  libmesh_convergence_failure();
1076  x(0) = b(0) / _coords[0];
1077 #endif
1078 
1079 #if LIBMESH_DIM == 2
1080  T my_det = _coords[0]*_coords[3] - _coords[1]*_coords[2];
1081 
1082  if (my_det == static_cast<T>(0.))
1083  libmesh_convergence_failure();
1084 
1085  T my_det_inv = 1./my_det;
1086 
1087  x(0) = my_det_inv*( _coords[3]*b(0) - _coords[1]*b(1));
1088  x(1) = my_det_inv*(-_coords[2]*b(0) + _coords[0]*b(1));
1089 #endif
1090 
1091 #if LIBMESH_DIM == 3
1092  T my_det =
1093  // a11*(a33 *a22 - a32 *a23)
1094  _coords[0]*(_coords[8]*_coords[4] - _coords[7]*_coords[5])
1095  // -a21*(a33 *a12 - a32 *a13)
1096  -_coords[3]*(_coords[8]*_coords[1] - _coords[7]*_coords[2]) +
1097  // +a31*(a23 *a12 - a22 *a13)
1098  +_coords[6]*(_coords[5]*_coords[1] - _coords[4]*_coords[2]);
1099 
1100  if (my_det == static_cast<T>(0.))
1101  libmesh_convergence_failure();
1102 
1103  T my_det_inv = 1./my_det;
1104  x(0) = my_det_inv*((_coords[8]*_coords[4] - _coords[7]*_coords[5])*b(0) -
1105  (_coords[8]*_coords[1] - _coords[7]*_coords[2])*b(1) +
1106  (_coords[5]*_coords[1] - _coords[4]*_coords[2])*b(2));
1107 
1108  x(1) = my_det_inv*((_coords[6]*_coords[5] - _coords[8]*_coords[3])*b(0) +
1109  (_coords[8]*_coords[0] - _coords[6]*_coords[2])*b(1) -
1110  (_coords[5]*_coords[0] - _coords[3]*_coords[2])*b(2));
1111 
1112  x(2) = my_det_inv*((_coords[7]*_coords[3] - _coords[6]*_coords[4])*b(0) -
1113  (_coords[7]*_coords[0] - _coords[6]*_coords[1])*b(1) +
1114  (_coords[4]*_coords[0] - _coords[3]*_coords[1])*b(2));
1115 #endif
1116 }
1117 
1118 
1119 
1120 template <typename T>
1121 inline
1123 {
1124  libmesh_assert_not_equal_to (factor, static_cast<T>(0.));
1125 
1126  for (unsigned int i=0; i<LIBMESH_DIM*LIBMESH_DIM; i++)
1127  _coords[i] /= factor;
1128 
1129  return *this;
1130 }
1131 
1132 
1133 
1134 
1135 template <typename T>
1136 template <typename T2>
1137 inline
1140 {
1142  for (unsigned int i=0; i<LIBMESH_DIM; i++)
1143  for (unsigned int j=0; j<LIBMESH_DIM; j++)
1144  returnval(i) += (*this)(i,j)*p(j);
1145 
1146  return returnval;
1147 }
1148 
1149 
1150 
1151 template <typename T>
1152 template <typename T2>
1153 inline
1155 {
1156  TypeTensor<T> returnval;
1157  for (unsigned int i=0; i<LIBMESH_DIM; i++)
1158  for (unsigned int j=0; j<LIBMESH_DIM; j++)
1159  for (unsigned int k=0; k<LIBMESH_DIM; k++)
1160  returnval(i,j) += (*this)(i,k)*p(k,j);
1161 
1162  return returnval;
1163 }
1164 
1165 
1166 
1171 template <typename T>
1172 template <typename T2>
1173 inline
1176 {
1177  typename CompareTypes<T,T2>::supertype sum = 0.;
1178  for (unsigned int i=0; i<LIBMESH_DIM*LIBMESH_DIM; i++)
1179  sum += _coords[i]*t._coords[i];
1180  return sum;
1181 }
1182 
1183 
1184 
1185 #ifdef LIBMESH_ENABLE_DEPRECATED
1186 template <typename T>
1187 inline
1189 {
1190  libmesh_deprecated();
1191  return this->norm();
1192 }
1193 #endif
1194 
1195 
1196 
1197 template <typename T>
1198 inline
1200 {
1201  return std::sqrt(this->norm_sq());
1202 }
1203 
1204 
1205 
1206 template <typename T>
1207 inline
1209 {
1210 #if LIBMESH_DIM == 1
1211  return _coords[0];
1212 #endif
1213 
1214 #if LIBMESH_DIM == 2
1215  return (_coords[0] * _coords[3]
1216  - _coords[1] * _coords[2]);
1217 #endif
1218 
1219 #if LIBMESH_DIM == 3
1220  return (_coords[0] * _coords[4] * _coords[8]
1221  + _coords[1] * _coords[5] * _coords[6]
1222  + _coords[2] * _coords[3] * _coords[7]
1223  - _coords[0] * _coords[5] * _coords[7]
1224  - _coords[1] * _coords[3] * _coords[8]
1225  - _coords[2] * _coords[4] * _coords[6]);
1226 #endif
1227 }
1228 
1229 template <typename T>
1230 inline
1232 {
1233 #if LIBMESH_DIM == 1
1234  return _coords[0];
1235 #endif
1236 
1237 #if LIBMESH_DIM == 2
1238  return _coords[0] + _coords[3];
1239 #endif
1240 
1241 #if LIBMESH_DIM == 3
1242  return _coords[0] + _coords[4] + _coords[8];
1243 #endif
1244 }
1245 
1246 template <typename T>
1247 inline
1249 {
1250  for (unsigned int i=0; i<LIBMESH_DIM*LIBMESH_DIM; i++)
1251  _coords[i] = 0.;
1252 }
1253 
1254 
1255 
1256 #ifdef LIBMESH_ENABLE_DEPRECATED
1257 template <typename T>
1258 inline
1260 {
1261  libmesh_deprecated();
1262  return this->norm_sq();
1263 }
1264 #endif
1265 
1266 
1267 
1268 template <typename T>
1269 inline
1271 {
1272  Real sum = 0.;
1273  for (unsigned int i=0; i<LIBMESH_DIM*LIBMESH_DIM; i++)
1274  sum += TensorTools::norm_sq(_coords[i]);
1275  return sum;
1276 }
1277 
1278 
1279 
1280 template <typename T>
1281 inline
1283 {
1284 #if LIBMESH_DIM == 1
1285  return (std::abs(_coords[0] - rhs._coords[0])
1286  < TOLERANCE);
1287 #endif
1288 
1289 #if LIBMESH_DIM == 2
1290  return ((std::abs(_coords[0] - rhs._coords[0]) +
1291  std::abs(_coords[1] - rhs._coords[1]) +
1292  std::abs(_coords[2] - rhs._coords[2]) +
1293  std::abs(_coords[3] - rhs._coords[3]))
1294  < 4.*TOLERANCE);
1295 #endif
1296 
1297 #if LIBMESH_DIM == 3
1298  return ((std::abs(_coords[0] - rhs._coords[0]) +
1299  std::abs(_coords[1] - rhs._coords[1]) +
1300  std::abs(_coords[2] - rhs._coords[2]) +
1301  std::abs(_coords[3] - rhs._coords[3]) +
1302  std::abs(_coords[4] - rhs._coords[4]) +
1303  std::abs(_coords[5] - rhs._coords[5]) +
1304  std::abs(_coords[6] - rhs._coords[6]) +
1305  std::abs(_coords[7] - rhs._coords[7]) +
1306  std::abs(_coords[8] - rhs._coords[8]))
1307  < 9.*TOLERANCE);
1308 #endif
1309 
1310 }
1311 
1312 } // namespace libMesh
1313 
1314 #endif // LIBMESH_TYPE_TENSOR_H
T _coords[LIBMESH_DIM]
The coordinates of the TypeVector.
Definition: type_vector.h:406
const TypeTensor< T > & operator*=(const Scalar)
Multiply this tensor by a scalar value in place.
Definition: type_tensor.h:940
bool operator==(const TypeTensor< T > &rhs) const
Definition: type_tensor.h:1282
double abs(double a)
const unsigned int _j
Definition: type_tensor.h:452
TypeTensor< T > * _tensor
Definition: type_tensor.h:451
Real norm_sq() const
Definition: type_tensor.h:1270
T value_type
Helper typedef for C++98 generic programming.
Definition: type_tensor.h:117
void add_scaled(const TypeTensor< T2 > &, const T)
Add a scaled tensor to this tensor without creating a temporary.
Definition: type_tensor.h:777
TypeTensor< typename CompareTypes< T, T2 >::supertype > operator+(const TypeTensor< T2 > &) const
Add another tensor to this tensor.
Definition: type_tensor.h:720
void zero()
Set all entries of the tensor to 0.
Definition: type_tensor.h:1248
TypeTensor()
Empty constructor.
Definition: type_tensor.h:483
const T & slice(const unsigned int i) const
Definition: type_tensor.h:471
ConstTypeTensorColumn(const TypeTensor< T > &tensor, unsigned int j)
Definition: type_tensor.h:461
static const Real TOLERANCE
The libMesh namespace provides an interface to certain functionality in the library.
void sum(T &r, const Communicator &comm=Communicator_World)
T _coords[LIBMESH_DIM *LIBMESH_DIM]
The coordinates of the TypeTensor.
Definition: type_tensor.h:417
const TypeTensor< T > * _tensor
Definition: type_tensor.h:475
void print(std::ostream &os=libMesh::out) const
Formatted print to a stream which defaults to libMesh::out.
Definition: type_tensor.C:39
ConstTypeTensorColumn< T > slice(const unsigned int i) const
Definition: type_tensor.h:686
This class defines a tensor in LIBMESH_DIM dimensional space of type T.
Definition: tensor_tools.h:32
TypeTensor< T > transpose() const
Definition: type_tensor.h:992
Real norm() const
Definition: type_tensor.h:1199
TypeTensorColumn(TypeTensor< T > &tensor, unsigned int j)
Definition: type_tensor.h:426
PetscErrorCode Vec x
Real size() const
Definition: type_tensor.h:1188
Real size_sq() const
Definition: type_tensor.h:1259
void write_unformatted(std::ostream &out, const bool newline=true) const
Unformatted print to the stream out.
Definition: type_tensor.C:83
TypeTensor< T > operator-() const
Definition: type_tensor.h:857
const TypeTensor< T > & operator/=(const T)
Divide each entry of this tensor by a scalar value.
Definition: type_tensor.h:1122
void subtract(const TypeTensor< T2 > &)
Subtract from this tensor without creating a temporary.
Definition: type_tensor.h:836
void assign(const TypeTensor< T2 > &)
Assign to this tensor without creating a temporary.
Definition: type_tensor.h:637
This class defines a vector in LIBMESH_DIM dimensional space of type T.
Definition: tensor_tools.h:30
PetscErrorCode Vec Mat libmesh_dbg_var(j)
void add(const TypeTensor< T2 > &)
Add to this tensor without creating a temporary.
Definition: type_tensor.h:766
~TypeTensor()
Destructor.
Definition: type_tensor.h:628
TypeVector< T > row(const unsigned int r) const
Definition: type_tensor.h:706
TypeTensor< T > inverse() const
Definition: type_tensor.h:1022
void solve(const TypeVector< T > &b, TypeVector< T > &x) const
Solve the 2x2 or 3x3 system of equations A*x = b for x by directly inverting A and multiplying it by ...
Definition: type_tensor.h:1071
const T & operator()(const unsigned int i, const unsigned int j) const
Definition: type_tensor.h:647
This class will eventually define a rank-N tensor in LIBMESH_DIM dimensional space of type T...
Definition: tensor_tools.h:34
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
static PetscErrorCode Mat * A
CompareTypes< T, T2 >::supertype contract(const TypeTensor< T2 > &) const
Multiply 2 tensors together, i.e.
Definition: type_tensor.h:1175
OStreamProxy out
boostcopy::enable_if_c< ScalarTraits< Scalar >::value, TypeTensor & >::type operator=(const Scalar &libmesh_dbg_var(p))
Assignment-from-scalar operator.
Definition: type_tensor.h:145
const TypeTensor< T > & operator-=(const TypeTensor< T2 > &)
Subtract from this tensor.
Definition: type_tensor.h:824
const TypeTensor< T > & operator+=(const TypeTensor< T2 > &)
Add to this tensor.
Definition: type_tensor.h:754
boostcopy::enable_if_c< ScalarTraits< Scalar >::value, TypeTensor< typename CompareTypes< T, Scalar >::supertype > >::type operator*(const Scalar) const
Multiply this tensor by a scalar value.
Definition: type_tensor.h:893
boostcopy::enable_if_c< ScalarTraits< Scalar >::value, TypeTensor< typename CompareTypes< T, Scalar >::supertype > >::type operator/(const Scalar) const
Divide each entry of this tensor by a scalar value.
Definition: type_tensor.h:957
T & slice(const unsigned int i)
Definition: type_tensor.h:437
bool operator>(const TypeTensor< T > &rhs) const
void subtract_scaled(const TypeTensor< T2 > &, const T)
Subtract a scaled value from this tensor without creating a temporary.
Definition: type_tensor.h:847