libMesh
type_vector.h
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2025 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_VECTOR_H
21 #define LIBMESH_TYPE_VECTOR_H
22 
23 // Local includes
24 #include "libmesh/libmesh_common.h"
25 #include "libmesh/compare_types.h"
26 #include "libmesh/tensor_tools.h"
27 #include "libmesh/int_range.h"
28 #include "libmesh/fuzzy_equals.h"
29 
30 // C++ includes
31 #include <cstdlib> // *must* precede <cmath> for proper std:abs() on PGI, Sun Studio CC
32 #include <cmath>
33 #include <complex>
34 #include <iostream>
35 #include <iomanip> // for std::setw, std::setiosflags
36 
37 namespace libMesh
38 {
39 // Forward declarations
40 template <typename T> class TypeTensor;
41 template <typename T> class VectorValue;
42 template <typename T> class TensorValue;
43 }
44 
45 namespace std
46 {
47 template <typename T>
48 auto norm(const libMesh::TypeVector<T> & vector) -> decltype(std::norm(T()));
49 }
50 
51 namespace libMesh
52 {
61 template <typename T>
62 class TypeVector
63 {
64  template <typename T2>
65  friend class TypeVector;
66 
67  friend class TypeTensor<T>;
68 
69 public:
74  TypeVector ();
75 
76 protected:
81  TypeVector (const T & x,
82  const T & y=0,
83  const T & z=0);
84 
89  template <typename Scalar1, typename Scalar2, typename Scalar3>
90  TypeVector (typename
91  boostcopy::enable_if_c<ScalarTraits<Scalar1>::value,
92  const Scalar1>::type & x,
93  typename
94  boostcopy::enable_if_c<ScalarTraits<Scalar2>::value,
95  const Scalar2>::type & y=0,
96  typename
97  boostcopy::enable_if_c<ScalarTraits<Scalar3>::value,
98  const Scalar3>::type & z=0);
99 
106  template <typename Scalar>
107  TypeVector (const Scalar & x,
108  typename
109  boostcopy::enable_if_c<ScalarTraits<Scalar>::value,
110  const Scalar>::type * sfinae = nullptr);
111 
112 public:
113 
117  typedef T value_type;
118 
122  typedef unsigned int index_type;
123 
127  template <typename T2>
128  TypeVector (const TypeVector<T2> & p);
129 
133  TypeVector (const TypeVector & p) = default;
134 
138  ~TypeVector () = default;
139 
143  template <typename T2>
144  void assign (const TypeVector<T2> &);
145 
149  template <typename Scalar>
150  typename boostcopy::enable_if_c<
151  ScalarTraits<Scalar>::value,
152  TypeVector &>::type
153  operator = (const Scalar & libmesh_dbg_var(p))
154  { libmesh_assert_equal_to (p, Scalar(0)); this->zero(); return *this; }
155 
159  const T & operator () (const unsigned int i) const;
160  const T & slice (const unsigned int i) const { return (*this)(i); }
161 
165  T & operator () (const unsigned int i);
166  T & slice (const unsigned int i) { return (*this)(i); }
167 
173  template <typename T2>
174  TypeVector<typename CompareTypes<T, T2>::supertype>
175  operator + (const TypeVector<T2> &) const;
176 
182  template <typename T2>
183  const TypeVector<T> & operator += (const TypeVector<T2> &);
184 
188  template <typename T2>
189  void add (const TypeVector<T2> &);
190 
194  template <typename T2>
195  void add_scaled (const TypeVector<T2> &, const T &);
196 
202  template <typename T2>
203  TypeVector<typename CompareTypes<T, T2>::supertype>
204  operator - (const TypeVector<T2> &) const;
205 
211  template <typename T2>
212  const TypeVector<T> & operator -= (const TypeVector<T2> &);
213 
217  template <typename T2>
218  void subtract (const TypeVector<T2> &);
219 
224  template <typename T2>
225  void subtract_scaled (const TypeVector<T2> &, const T &);
226 
230  TypeVector<T> operator - () const;
231 
237  template <typename Scalar>
238  typename boostcopy::enable_if_c<
239  ScalarTraits<Scalar>::value,
240  TypeVector<typename CompareTypes<T, Scalar>::supertype>>::type
241  operator * (const Scalar &) const;
242 
248  const TypeVector<T> & operator *= (const T &);
249 
255  template <typename Scalar>
256  typename boostcopy::enable_if_c<
257  ScalarTraits<Scalar>::value,
258  TypeVector<typename CompareTypes<T, Scalar>::supertype>>::type
259  operator / (const Scalar &) const;
260 
266  const TypeVector<T> & operator /= (const T &);
267 
274  template <typename T2>
275  typename CompareTypes<T, T2>::supertype
276  operator * (const TypeVector<T2> &) const;
277 
281  template <typename T2>
282  typename CompareTypes<T, T2>::supertype
283  contract (const TypeVector<T2> &) const;
284 
288  template <typename T2>
289  TypeVector<typename CompareTypes<T, T2>::supertype>
290  cross(const TypeVector<T2> & v) const;
291 
295  TypeVector<T> unit() const;
296 
301  auto norm() const -> decltype(std::norm(T()));
302 
307  auto norm_sq() const -> decltype(std::norm(T()));
308 
312  auto l1_norm() const;
313 
317  bool is_zero() const;
318 
322  void zero();
323 
328  bool relative_fuzzy_equals(const TypeVector<T> & rhs,
329  Real tol = TOLERANCE) const;
330 
335  bool absolute_fuzzy_equals(const TypeVector<T> & rhs,
336  Real tol = TOLERANCE) const;
337 
345  bool operator == (const TypeVector<T> & rhs) const;
346 
350  bool operator != (const TypeVector<T> & rhs) const;
351 
358  bool operator < (const TypeVector<T> & rhs) const;
359 
366  bool operator <= (const TypeVector<T> & rhs) const;
367 
374  bool operator > (const TypeVector<T> & rhs) const;
375 
382  bool operator >= (const TypeVector<T> & rhs) const;
383 
387  void print(std::ostream & os = libMesh::out) const;
388 
397  friend std::ostream & operator << (std::ostream & os, const TypeVector<T> & t)
398  {
399  t.print(os);
400  return os;
401  }
402 
409  void write_unformatted (std::ostream & out_stream,
410  const bool newline = true) const;
411 
412 protected:
413 
417  T _coords[LIBMESH_DIM];
418 };
419 
420 
421 
422 
423 
424 //------------------------------------------------------
425 // Inline functions
426 
427 template <typename T>
428 inline
429 TypeVector<T>::TypeVector ()
430 {
431  _coords[0] = {};
432 
433 #if LIBMESH_DIM > 1
434  _coords[1] = {};
435 #endif
436 
437 #if LIBMESH_DIM > 2
438  _coords[2] = {};
439 #endif
440 }
441 
442 
443 
444 template <typename T>
445 inline
446 TypeVector<T>::TypeVector (const T & x,
447  const T & y,
448  const T & z)
449 {
450  _coords[0] = x;
451 
452 #if LIBMESH_DIM > 1
453  _coords[1] = y;
454 #else
455  libmesh_ignore(y);
456  libmesh_assert_equal_to (y, 0);
457 #endif
458 
459 #if LIBMESH_DIM > 2
460  _coords[2] = z;
461 #else
462  libmesh_ignore(z);
463  libmesh_assert_equal_to (z, 0);
464 #endif
465 }
466 
467 
468 template <typename T>
469 template <typename Scalar1, typename Scalar2, typename Scalar3>
470 inline
471 TypeVector<T>::TypeVector (typename
472  boostcopy::enable_if_c<ScalarTraits<Scalar1>::value,
473  const Scalar1>::type & x,
474  typename
475  boostcopy::enable_if_c<ScalarTraits<Scalar2>::value,
476  const Scalar2>::type & y,
477  typename
478  boostcopy::enable_if_c<ScalarTraits<Scalar3>::value,
479  const Scalar3>::type & z)
480 {
481  _coords[0] = x;
482 
483 #if LIBMESH_DIM > 1
484  _coords[1] = y;
485 #else
486  libmesh_assert_equal_to (y, 0);
487 #endif
488 
489 #if LIBMESH_DIM > 2
490  _coords[2] = z;
491 #else
492  libmesh_assert_equal_to (z, 0);
493 #endif
494 }
495 
496 
497 
498 template <typename T>
499 template <typename Scalar>
500 inline
501 TypeVector<T>::TypeVector (const Scalar & x,
502  typename
503  boostcopy::enable_if_c<ScalarTraits<Scalar>::value,
504  const Scalar>::type * /*sfinae*/)
505 {
506  _coords[0] = x;
507 
508 #if LIBMESH_DIM > 1
509  _coords[1] = 0;
510 #endif
511 
512 #if LIBMESH_DIM > 2
513  _coords[2] = 0;
514 #endif
515 }
516 
517 
518 
519 template <typename T>
520 template <typename T2>
521 inline
522 TypeVector<T>::TypeVector (const TypeVector<T2> & p)
523 {
524  // copy the nodes from vector p to me
525  for (unsigned int i=0; i<LIBMESH_DIM; i++)
526  _coords[i] = p._coords[i];
527 }
528 
529 
530 
531 template <typename T>
532 template <typename T2>
533 inline
534 void TypeVector<T>::assign (const TypeVector<T2> & p)
535 {
536  for (unsigned int i=0; i<LIBMESH_DIM; i++)
537  _coords[i] = p._coords[i];
538 }
539 
540 
541 
542 template <typename T>
543 inline
544 const T & TypeVector<T>::operator () (const unsigned int i) const
545 {
546  libmesh_assert_less (i, LIBMESH_DIM);
547 
548  return _coords[i];
549 }
550 
551 
552 
553 template <typename T>
554 inline
555 T & TypeVector<T>::operator () (const unsigned int i)
556 {
557  libmesh_assert_less (i, LIBMESH_DIM);
558 
559  return _coords[i];
560 }
561 
562 
563 
564 template <typename T>
565 template <typename T2>
566 inline
567 TypeVector<typename CompareTypes<T, T2>::supertype>
568 TypeVector<T>::operator + (const TypeVector<T2> & p) const
569 {
570  typedef typename CompareTypes<T, T2>::supertype TS;
571 #if LIBMESH_DIM == 1
572  return TypeVector<TS> (_coords[0] + p._coords[0]);
573 #endif
574 
575 #if LIBMESH_DIM == 2
576  return TypeVector<TS> (_coords[0] + p._coords[0],
577  _coords[1] + p._coords[1]);
578 #endif
579 
580 #if LIBMESH_DIM == 3
581  return TypeVector<TS> (_coords[0] + p._coords[0],
582  _coords[1] + p._coords[1],
583  _coords[2] + p._coords[2]);
584 #endif
585 
586 }
587 
588 
589 
590 template <typename T>
591 template <typename T2>
592 inline
593 const TypeVector<T> & TypeVector<T>::operator += (const TypeVector<T2> & p)
594 {
595  this->add (p);
596 
597  return *this;
598 }
599 
600 
601 
602 template <typename T>
603 template <typename T2>
604 inline
605 void TypeVector<T>::add (const TypeVector<T2> & p)
606 {
607 #if LIBMESH_DIM == 1
608  _coords[0] += p._coords[0];
609 #endif
610 
611 #if LIBMESH_DIM == 2
612  _coords[0] += p._coords[0];
613  _coords[1] += p._coords[1];
614 #endif
615 
616 #if LIBMESH_DIM == 3
617  _coords[0] += p._coords[0];
618  _coords[1] += p._coords[1];
619  _coords[2] += p._coords[2];
620 #endif
621 
622 }
623 
624 
625 
626 template <typename T>
627 template <typename T2>
628 inline
629 void TypeVector<T>::add_scaled (const TypeVector<T2> & p, const T & factor)
630 {
631 #if LIBMESH_DIM == 1
632  _coords[0] += factor*p(0);
633 #endif
634 
635 #if LIBMESH_DIM == 2
636  _coords[0] += factor*p(0);
637  _coords[1] += factor*p(1);
638 #endif
639 
640 #if LIBMESH_DIM == 3
641  _coords[0] += factor*p(0);
642  _coords[1] += factor*p(1);
643  _coords[2] += factor*p(2);
644 #endif
645 
646 }
647 
648 
649 
650 template <typename T>
651 template <typename T2>
652 inline
653 TypeVector<typename CompareTypes<T, T2>::supertype>
654 TypeVector<T>::operator - (const TypeVector<T2> & p) const
655 {
656  typedef typename CompareTypes<T, T2>::supertype TS;
657 
658 #if LIBMESH_DIM == 1
659  return TypeVector<TS>(_coords[0] - p._coords[0]);
660 #endif
661 
662 #if LIBMESH_DIM == 2
663  return TypeVector<TS>(_coords[0] - p._coords[0],
664  _coords[1] - p._coords[1]);
665 #endif
666 
667 #if LIBMESH_DIM == 3
668  return TypeVector<TS>(_coords[0] - p._coords[0],
669  _coords[1] - p._coords[1],
670  _coords[2] - p._coords[2]);
671 #endif
672 
673 }
674 
675 
676 
677 template <typename T>
678 template <typename T2>
679 inline
680 const TypeVector<T> & TypeVector<T>::operator -= (const TypeVector<T2> & p)
681 {
682  this->subtract (p);
683 
684  return *this;
685 }
686 
687 
688 
689 template <typename T>
690 template <typename T2>
691 inline
692 void TypeVector<T>::subtract (const TypeVector<T2> & p)
693 {
694  for (unsigned int i=0; i<LIBMESH_DIM; i++)
695  _coords[i] -= p._coords[i];
696 }
697 
698 
699 
700 template <typename T>
701 template <typename T2>
702 inline
703 void TypeVector<T>::subtract_scaled (const TypeVector<T2> & p, const T & factor)
704 {
705  for (unsigned int i=0; i<LIBMESH_DIM; i++)
706  _coords[i] -= factor*p(i);
707 }
708 
709 
710 
711 template <typename T>
712 inline
713 TypeVector<T> TypeVector<T>::operator - () const
714 {
715 
716 #if LIBMESH_DIM == 1
717  return TypeVector(-_coords[0]);
718 #endif
719 
720 #if LIBMESH_DIM == 2
721  return TypeVector(-_coords[0],
722  -_coords[1]);
723 #endif
724 
725 #if LIBMESH_DIM == 3
726  return TypeVector(-_coords[0],
727  -_coords[1],
728  -_coords[2]);
729 #endif
730 
731 }
732 
733 
734 
735 template <typename T>
736 template <typename Scalar>
737 inline
738 typename boostcopy::enable_if_c<
739  ScalarTraits<Scalar>::value,
740  TypeVector<typename CompareTypes<T, Scalar>::supertype>>::type
741 TypeVector<T>::operator * (const Scalar & factor) const
742 {
743  typedef typename CompareTypes<T, Scalar>::supertype SuperType;
744 
745 #if LIBMESH_DIM == 1
746  return TypeVector<SuperType>(_coords[0]*factor);
747 #endif
748 
749 #if LIBMESH_DIM == 2
750  return TypeVector<SuperType>(_coords[0]*factor,
751  _coords[1]*factor);
752 #endif
753 
754 #if LIBMESH_DIM == 3
755  return TypeVector<SuperType>(_coords[0]*factor,
756  _coords[1]*factor,
757  _coords[2]*factor);
758 #endif
759 }
760 
761 
762 
763 template <typename T, typename Scalar>
764 inline
765 typename boostcopy::enable_if_c<
766  ScalarTraits<Scalar>::value,
767  TypeVector<typename CompareTypes<T, Scalar>::supertype>>::type
768 operator * (const Scalar & factor,
769  const TypeVector<T> & v)
770 {
771  return v * factor;
772 }
773 
774 
775 
776 template <typename T>
777 inline
778 const TypeVector<T> & TypeVector<T>::operator *= (const T & factor)
779 {
780 #if LIBMESH_DIM == 1
781  _coords[0] *= factor;
782 #endif
783 
784 #if LIBMESH_DIM == 2
785  _coords[0] *= factor;
786  _coords[1] *= factor;
787 #endif
788 
789 #if LIBMESH_DIM == 3
790  _coords[0] *= factor;
791  _coords[1] *= factor;
792  _coords[2] *= factor;
793 #endif
794 
795  return *this;
796 }
797 
798 
799 
800 template <typename T>
801 template <typename Scalar>
802 inline
803 typename boostcopy::enable_if_c<
804  ScalarTraits<Scalar>::value,
805  TypeVector<typename CompareTypes<T, Scalar>::supertype>>::type
806 TypeVector<T>::operator / (const Scalar & factor) const
807 {
808  libmesh_assert_not_equal_to (factor, static_cast<T>(0.));
809 
810  typedef typename CompareTypes<T, Scalar>::supertype TS;
811 
812 #if LIBMESH_DIM == 1
813  return TypeVector<TS>(_coords[0]/factor);
814 #endif
815 
816 #if LIBMESH_DIM == 2
817  return TypeVector<TS>(_coords[0]/factor,
818  _coords[1]/factor);
819 #endif
820 
821 #if LIBMESH_DIM == 3
822  return TypeVector<TS>(_coords[0]/factor,
823  _coords[1]/factor,
824  _coords[2]/factor);
825 #endif
826 
827 }
828 
829 
830 
831 
832 template <typename T>
833 inline
834 const TypeVector<T> &
835 TypeVector<T>::operator /= (const T & factor)
836 {
837  libmesh_assert_not_equal_to (factor, static_cast<T>(0.));
838 
839  for (unsigned int i=0; i<LIBMESH_DIM; i++)
840  _coords[i] /= factor;
841 
842  return *this;
843 }
844 
845 
846 
847 
848 template <typename T>
849 template <typename T2>
850 inline
851 typename CompareTypes<T, T2>::supertype
852 TypeVector<T>::operator * (const TypeVector<T2> & p) const
853 {
854 #if LIBMESH_DIM == 1
855  return _coords[0]*p._coords[0];
856 #endif
857 
858 #if LIBMESH_DIM == 2
859  return (_coords[0]*p._coords[0] +
860  _coords[1]*p._coords[1]);
861 #endif
862 
863 #if LIBMESH_DIM == 3
864  return (_coords[0]*p(0) +
865  _coords[1]*p(1) +
866  _coords[2]*p(2));
867 #endif
868 }
869 
870 template <typename T>
871 template <typename T2>
872 inline
873 typename CompareTypes<T, T2>::supertype
874 TypeVector<T>::contract(const TypeVector<T2> & p) const
875 {
876  return (*this)*(p);
877 }
878 
879 
880 
881 template <typename T>
882 template <typename T2>
883 TypeVector<typename CompareTypes<T, T2>::supertype>
884 TypeVector<T>::cross(const TypeVector<T2> & p) const
885 {
886  typedef typename CompareTypes<T, T2>::supertype TS;
887  libmesh_assert_equal_to (LIBMESH_DIM, 3);
888 
889  // | i j k |
890  // |(*this)(0) (*this)(1) (*this)(2)|
891  // | p(0) p(1) p(2) |
892 
893 #if LIBMESH_DIM == 3
894  return TypeVector<TS>( _coords[1]*p._coords[2] - _coords[2]*p._coords[1],
895  -_coords[0]*p._coords[2] + _coords[2]*p._coords[0],
896  _coords[0]*p._coords[1] - _coords[1]*p._coords[0]);
897 #else
898  libmesh_ignore(p);
899  return TypeVector<TS>(0);
900 #endif
901 }
902 
903 
904 
905 template <typename T>
906 inline
907 auto TypeVector<T>::norm() const -> decltype(std::norm(T()))
908 {
909  return std::sqrt(this->norm_sq());
910 }
911 
912 
913 
914 template <typename T>
915 inline
916 void TypeVector<T>::zero()
917 {
918  for (unsigned int i=0; i<LIBMESH_DIM; i++)
919  _coords[i] = 0.;
920 }
921 
922 
923 
924 template <typename T>
925 inline
926 auto TypeVector<T>::norm_sq() const -> decltype(std::norm(T()))
927 {
928 #if LIBMESH_DIM == 1
929  return (TensorTools::norm_sq(_coords[0]));
930 #endif
931 
932 #if LIBMESH_DIM == 2
933  return (TensorTools::norm_sq(_coords[0]) +
934  TensorTools::norm_sq(_coords[1]));
935 #endif
936 
937 #if LIBMESH_DIM == 3
938  return (TensorTools::norm_sq(_coords[0]) +
939  TensorTools::norm_sq(_coords[1]) +
940  TensorTools::norm_sq(_coords[2]));
941 #endif
942 }
943 
944 
945 template <typename T>
946 inline
947 bool TypeVector<T>::is_zero() const
948 {
949  for (const auto & val : _coords)
950  if (val != T(0))
951  return false;
952  return true;
953 }
954 
955 template <>
956 auto
957 TypeVector<bool>::l1_norm() const;
958 
959 template <typename T>
960 auto
961 TypeVector<T>::l1_norm() const
962 {
963  decltype(std::abs(T())) ret{};
964  for (const auto i : make_range(libmesh_dim))
965  ret += std::abs(_coords[i]);
966 
967  return ret;
968 }
969 
970 template <typename T>
971 inline
972 bool TypeVector<T>::absolute_fuzzy_equals(const TypeVector<T> & rhs, Real tol) const
973 {
974  return libMesh::absolute_fuzzy_equals(*this, rhs, tol);
975 }
976 
977 
978 
979 template <typename T>
980 inline
981 bool TypeVector<T>::relative_fuzzy_equals(const TypeVector<T> & rhs, Real tol) const
982 {
983  return libMesh::relative_fuzzy_equals(*this, rhs, tol);
984 }
985 
986 
987 
988 template <typename T>
989 inline
990 bool TypeVector<T>::operator == (const TypeVector<T> & rhs) const
991 {
992 #if LIBMESH_DIM == 1
993  return (_coords[0] == rhs._coords[0]);
994 #endif
995 
996 #if LIBMESH_DIM == 2
997  return (_coords[0] == rhs._coords[0] &&
998  _coords[1] == rhs._coords[1]);
999 #endif
1000 
1001 #if LIBMESH_DIM == 3
1002  return (_coords[0] == rhs._coords[0] &&
1003  _coords[1] == rhs._coords[1] &&
1004  _coords[2] == rhs._coords[2]);
1005 #endif
1006 }
1007 
1008 
1009 
1010 template <typename T>
1011 inline
1012 bool TypeVector<T>::operator != (const TypeVector<T> & rhs) const
1013 {
1014  return (!(*this == rhs));
1015 }
1016 
1017 
1018 //------------------------------------------------------
1019 // Non-member functions on TypeVectors
1020 
1021 // Compute a * (b.cross(c)) without creating a temporary
1022 // for the cross product. Equivalent to the determinant
1023 // of the 3x3 tensor:
1024 // [a0, a1, a2]
1025 // [b0, b1, b2]
1026 // [c0, c1, c2]
1027 template <typename T>
1028 inline
1029 T triple_product(const TypeVector<T> & a,
1030  const TypeVector<T> & b,
1031  const TypeVector<T> & c)
1032 {
1033 #if LIBMESH_DIM == 3
1034  return
1035  a(0)*(b(1)*c(2) - b(2)*c(1)) -
1036  a(1)*(b(0)*c(2) - b(2)*c(0)) +
1037  a(2)*(b(0)*c(1) - b(1)*c(0));
1038 #else
1039  libmesh_ignore(a, b, c);
1040  return 0;
1041 #endif
1042 }
1043 
1044 
1045 // compute the solid angle subtended by the tetrahedral vertex 0, as
1046 // defined by vectors v01, v02, and v03. The solid angle is defined
1047 // to be positive if the vectors are obey the right-hand rule, or
1048 // negative for a left-hand orientation.
1049 template <typename T>
1050 inline
1051 T solid_angle(const TypeVector<T> & v01,
1052  const TypeVector<T> & v02,
1053  const TypeVector<T> & v03)
1054 {
1055  const Real norm01 = v01.norm(),
1056  norm02 = v02.norm(),
1057  norm03 = v03.norm();
1058  const T tan_half_angle =
1059  triple_product(v01, v02, v03) /
1060  ((v01*v02)*norm03 + (v01*v03)*norm02 + (v02*v03)*norm01 +
1061  norm01*norm02*norm03);
1062 
1063  return Real(2)*std::atan(tan_half_angle);
1064 }
1065 
1066 
1067 
1068 
1073 template <typename T>
1074 inline
1075 T cross_norm_sq(const TypeVector<T> & b,
1076  const TypeVector<T> & c)
1077 {
1078  T z = b(0)*c(1) - b(1)*c(0);
1079 
1080 #if LIBMESH_DIM == 3
1081  T x = b(1)*c(2) - b(2)*c(1),
1082  y = b(0)*c(2) - b(2)*c(0);
1083  return x*x + y*y + z*z;
1084 #else
1085  return z*z;
1086 #endif
1087 }
1088 
1089 
1090 
1094 template <typename T>
1095 inline
1096 T cross_norm(const TypeVector<T> & b,
1097  const TypeVector<T> & c)
1098 {
1099  return std::sqrt(cross_norm_sq(b,c));
1100 }
1101 
1102 template <typename T>
1103 inline
1104 TypeVector<T> TypeVector<T>::unit() const
1105 {
1106 
1107  auto && length = norm();
1108 
1109  libmesh_assert_not_equal_to (length, static_cast<Real>(0.));
1110 
1111 #if LIBMESH_DIM == 1
1112  return TypeVector<T>(_coords[0]/length);
1113 #endif
1114 
1115 #if LIBMESH_DIM == 2
1116  return TypeVector<T>(_coords[0]/length,
1117  _coords[1]/length);
1118 #endif
1119 
1120 #if LIBMESH_DIM == 3
1121  return TypeVector<T>(_coords[0]/length,
1122  _coords[1]/length,
1123  _coords[2]/length);
1124 #endif
1125 
1126 }
1127 
1128 template <typename T>
1129 void TypeVector<T>::print(std::ostream & os) const
1130 {
1131 #if LIBMESH_DIM == 1
1132 
1133  os << "x=" << (*this)(0);
1134 
1135 #endif
1136 #if LIBMESH_DIM == 2
1137 
1138  os << "(x,y)=("
1139  << std::setw(8) << (*this)(0) << ", "
1140  << std::setw(8) << (*this)(1) << ")";
1141 
1142 #endif
1143 #if LIBMESH_DIM == 3
1144 
1145  os << "(x,y,z)=("
1146  << std::setw(8) << (*this)(0) << ", "
1147  << std::setw(8) << (*this)(1) << ", "
1148  << std::setw(8) << (*this)(2) << ")";
1149 #endif
1150 }
1151 
1152 template <typename T>
1153 struct CompareTypes<TypeVector<T>, TypeVector<T>>
1154 {
1155  typedef TypeVector<T> supertype;
1156 };
1157 
1158 template <typename T, typename T2>
1159 struct CompareTypes<TypeVector<T>, TypeVector<T2>>
1160 {
1161  typedef TypeVector<typename CompareTypes<T,T2>::supertype> supertype;
1162 };
1163 
1164 template <typename T, typename T2, typename std::enable_if<ScalarTraits<T>::value, int>::type = 0>
1165 TypeVector<typename CompareTypes<T, T2>::supertype>
1166 outer_product(const T & a, const TypeVector<T2> & b)
1167 {
1168  TypeVector<typename CompareTypes<T, T2>::supertype> ret;
1169  for (unsigned int i = 0; i < LIBMESH_DIM; i++)
1170  ret(i) = a * libmesh_conj(b(i));
1171 
1172  return ret;
1173 }
1174 
1175 template <typename T, typename T2, typename std::enable_if<ScalarTraits<T2>::value, int>::type = 0>
1176 TypeVector<typename CompareTypes<T, T2>::supertype>
1177 outer_product(const TypeVector<T> & a, const T2 & b)
1178 {
1179  TypeVector<typename CompareTypes<T, T2>::supertype> ret;
1180  const auto conj_b = libmesh_conj(b);
1181  for (unsigned int i = 0; i < LIBMESH_DIM; i++)
1182  ret(i) = a(i) * conj_b;
1183 
1184  return ret;
1185 }
1186 
1187 template <typename T>
1188 auto
1189 l1_norm(const TypeVector<T> & var)
1190 {
1191  return var.l1_norm();
1192 }
1193 
1194 template <typename T, typename T2>
1195 auto
1196 l1_norm_diff(const TypeVector<T> & vec1, const TypeVector<T2> & vec2)
1197 {
1198  return l1_norm(vec1 - vec2);
1199 }
1200 
1201 } // namespace libMesh
1202 
1203 namespace std
1204 {
1205 template <typename T>
1206 auto norm(const libMesh::TypeVector<T> & vector) -> decltype(std::norm(T()))
1207 {
1208  // Yea I agree it's dumb that the standard returns the square of the Euclidean norm
1209  return vector.norm_sq();
1210 }
1211 } // namespace std
1212 
1213 #ifdef LIBMESH_HAVE_METAPHYSICL
1214 namespace MetaPhysicL
1215 {
1216 template <typename T>
1217 struct RawType<libMesh::TypeVector<T>>
1218 {
1219  typedef libMesh::TypeVector<typename RawType<T>::value_type> value_type;
1220 
1221  static value_type value (const libMesh::TypeVector<T> & in)
1222  {
1223  value_type ret;
1224  for (unsigned int i = 0; i < LIBMESH_DIM; ++i)
1225  ret(i) = raw_value(in(i));
1226 
1227  return ret;
1228  }
1229 };
1230 
1231 template <typename T, typename U>
1232 struct ReplaceAlgebraicType<libMesh::TypeVector<T>, U>
1233 {
1234  typedef U type;
1235 };
1236 }
1237 #endif
1238 
1239 #endif // LIBMESH_TYPE_VECTOR_H
libMesh::TypeVector::_coords
T _coords[LIBMESH_DIM]
The coordinates of the TypeVector.
Definition: type_vector.h:417
libMesh::libmesh_dim
static constexpr std::size_t libmesh_dim
Definition: libmesh_common.h:640
libMesh::TypeVector::operator*=
const TypeVector< T > & operator*=(const T &)
Multiply this vector by a scalar value.
Definition: type_vector.h:778
libMesh::solid_angle
T solid_angle(const TypeVector< T > &v01, const TypeVector< T > &v02, const TypeVector< T > &v03)
Definition: type_vector.h:1051
libMesh::TypeVector::add_scaled
void add_scaled(const TypeVector< T2 > &, const T &)
Add a scaled value to this vector without creating a temporary.
Definition: type_vector.h:629
libMesh::outer_product
TypeTensor< typename CompareTypes< T, T2 >::supertype > outer_product(const TypeVector< T > &a, const TypeVector< T2 > &b)
Definition: type_tensor.h:1433
libMesh::TypeVector::norm
auto norm() const -> decltype(std::norm(T()))
Definition: type_vector.h:907
libMesh::libmesh_conj
T libmesh_conj(T a)
Definition: libmesh_common.h:188
libMesh::TypeVector::operator+
TypeVector< typename CompareTypes< T, T2 >::supertype > operator+(const TypeVector< T2 > &) const
Add two vectors.
Definition: type_vector.h:568
libMesh::TypeVector::subtract_scaled
void subtract_scaled(const TypeVector< T2 > &, const T &)
Subtract a scaled value from this vector without creating a temporary.
Definition: type_vector.h:703
MetaPhysicL::RawType< libMesh::TypeVector< T > >::value_type
libMesh::TypeVector< typename RawType< T >::value_type > value_type
Definition: type_vector.h:1219
libMesh::TypeVector::contract
CompareTypes< T, T2 >::supertype contract(const TypeVector< T2 > &) const
Definition: type_vector.h:874
libMesh::TOLERANCE
static constexpr Real TOLERANCE
Definition: libmesh_common.h:151
libMesh::cross_norm
T cross_norm(const TypeVector< T > &b, const TypeVector< T > &c)
Calls cross_norm_sq() and takes the square root of the result.
Definition: type_vector.h:1096
MetaPhysicL::RawType< libMesh::TypeVector< T > >::value
static value_type value(const libMesh::TypeVector< T > &in)
Definition: type_vector.h:1221
libMesh::CompareTypes< TypeVector< T >, TypeVector< T2 > >::supertype
TypeVector< typename CompareTypes< T, T2 >::supertype > supertype
Definition: type_vector.h:1161
libMesh::TypeVector::index_type
unsigned int index_type
Helper typedef for generic index programming.
Definition: type_vector.h:122
libMesh
The libMesh namespace provides an interface to certain functionality in the library.
libMesh::TypeVector::is_zero
bool is_zero() const
Definition: type_vector.h:947
libMesh::TypeVector::operator+=
const TypeVector< T > & operator+=(const TypeVector< T2 > &)
Add to this vector.
Definition: type_vector.h:593
libMesh::TypeVector::add
void add(const TypeVector< T2 > &)
Add to this vector without creating a temporary.
Definition: type_vector.h:605
libMesh::TypeTensor
This class defines a tensor in LIBMESH_DIM dimensional space of type T.
Definition: tensor_tools.h:36
libMesh::cross_norm_sq
T cross_norm_sq(const TypeVector< T > &b, const TypeVector< T > &c)
Compute |b x c|^2 without creating the extra temporary produced by calling b.cross(c).norm_sq().
Definition: type_vector.h:1075
libMesh::TypeVector::print
void print(std::ostream &os=libMesh::out) const
Formatted print, by default to libMesh::out.
Definition: type_vector.h:1129
libMesh::triple_product
T triple_product(const TypeVector< T > &a, const TypeVector< T > &b, const TypeVector< T > &c)
Definition: type_vector.h:1029
libMesh::TypeVector::unit
TypeVector< T > unit() const
Definition: type_vector.h:1104
libMesh::libmesh_ignore
void libmesh_ignore(const Args &...)
Definition: libmesh_common.h:515
libMesh::TypeVector::norm_sq
auto norm_sq() const -> decltype(std::norm(T()))
Definition: type_vector.h:926
libMesh::TypeVector::operator-
TypeVector< T > operator-() const
Definition: type_vector.h:713
libMesh::TypeVector::zero
void zero()
Set all entries of the vector to 0.
Definition: type_vector.h:916
libMesh::l1_norm
auto l1_norm(const NumericVector< T > &vec)
Definition: numeric_vector.h:1107
libMesh::TypeVector::subtract
void subtract(const TypeVector< T2 > &)
Subtract from this vector without creating a temporary.
Definition: type_vector.h:692
libMesh::operator*
boostcopy::enable_if_c< ScalarTraits< Scalar >::value, TypeNTensor< N, typename CompareTypes< Scalar, T >::supertype > >::type operator*(const Scalar &, const TypeNTensor< N, T > &)
Definition: type_n_tensor.h:323
libMesh::TypeVector::value_type
T value_type
Helper typedef for C++98 generic programming.
Definition: type_vector.h:117
libMesh::l1_norm_diff
auto l1_norm_diff(const NumericVector< T > &vec1, const NumericVector< T > &vec2)
Definition: numeric_vector.h:1114
libMesh::TypeVector
This class defines a vector in LIBMESH_DIM dimensional space of type T.
Definition: tensor_tools.h:34
libMesh::TypeVector::cross
TypeVector< typename CompareTypes< T, T2 >::supertype > cross(const TypeVector< T2 > &v) const
Definition: type_vector.h:884
libMesh::TypeVector::absolute_fuzzy_equals
bool absolute_fuzzy_equals(const TypeVector< T > &rhs, Real tol=TOLERANCE) const
Definition: type_vector.h:972
libMesh::boostcopy::enable_if_c
Definition: compare_types.h:38
libMesh::TypeVector::operator/=
const TypeVector< T > & operator/=(const T &)
Divide each entry of this vector by scalar value.
Definition: type_vector.h:835
libMesh::absolute_fuzzy_equals
bool absolute_fuzzy_equals(const T &var1, const T2 &var2, const Real tol=TOLERANCE *TOLERANCE)
Function to check whether two variables are equal within an absolute tolerance.
Definition: fuzzy_equals.h:64
libMesh::TypeVector::TypeVector
TypeVector()
Empty constructor.
Definition: type_vector.h:429
libMesh::TypeVector::operator-=
const TypeVector< T > & operator-=(const TypeVector< T2 > &)
Subtract from this vector.
Definition: type_vector.h:680
libMesh::Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Definition: libmesh_common.h:144
libMesh::TypeVector::operator*
boostcopy::enable_if_c< ScalarTraits< Scalar >::value, TypeVector< typename CompareTypes< T, Scalar >::supertype > >::type operator*(const Scalar &) const
Multiply this vector by a scalar value.
Definition: type_vector.h:741
libMesh::TensorTools::norm_sq
auto norm_sq(const T &a) -> decltype(std::norm(a))
Definition: tensor_tools.h:104
libMesh::out
OStreamProxy out
libMesh::TypeVector::assign
void assign(const TypeVector< T2 > &)
Assign to this vector without creating a temporary.
Definition: type_vector.h:534
value
static const bool value
Definition: xdr_io.C:54
libMesh::make_range
IntRange< T > make_range(T beg, T end)
The 2-parameter make_range() helper function returns an IntRange<T> when both input parameters are of...
Definition: int_range.h:140
libMesh::TypeVector::slice
T & slice(const unsigned int i)
Definition: type_vector.h:166
libMesh::TypeVector::operator()
const T & operator()(const unsigned int i) const
Definition: type_vector.h:544
libMesh::TypeVector::operator==
bool operator==(const TypeVector< T > &rhs) const
Definition: type_vector.h:990
libMesh::relative_fuzzy_equals
bool relative_fuzzy_equals(const T &var1, const T2 &var2, const Real tol=TOLERANCE *TOLERANCE)
Function to check whether two variables are equal within a relative tolerance.
Definition: fuzzy_equals.h:78
std::norm
auto norm(const libMesh::TypeVector< T > &vector) -> decltype(std::norm(T()))
Definition: type_vector.h:1206
libMesh::TypeVector::write_unformatted
void write_unformatted(std::ostream &out_stream, const bool newline=true) const
Unformatted print to the stream out.
Definition: type_vector.C:51
libMesh::TypeVector::operator/
boostcopy::enable_if_c< ScalarTraits< Scalar >::value, TypeVector< typename CompareTypes< T, Scalar >::supertype > >::type operator/(const Scalar &) const
Divide each entry of this vector by scalar value.
Definition: type_vector.h:806
libMesh::TypeVector::relative_fuzzy_equals
bool relative_fuzzy_equals(const TypeVector< T > &rhs, Real tol=TOLERANCE) const
Definition: type_vector.h:981
libMesh::TypeVector::operator=
boostcopy::enable_if_c< ScalarTraits< Scalar >::value, TypeVector & >::type operator=(const Scalar &libmesh_dbg_var(p))
Assignment-from-scalar operator.
Definition: type_vector.h:153
libMesh::TypeVector::l1_norm
auto l1_norm() const
Definition: type_vector.h:961
libMesh::TypeVector::operator!=
bool operator!=(const TypeVector< T > &rhs) const
Definition: type_vector.h:1012
libMesh::TypeVector::slice
const T & slice(const unsigned int i) const
Definition: type_vector.h:160
libMesh::TypeVector::~TypeVector
~TypeVector()=default
Destructor.
Generated on Sat Aug 2 2025 20:01:45 for libMesh by   doxygen 1.8.14