libMesh
Loading...
Searching...
No Matches
type_vector.h
Go to the documentation of this file.
1// The libMesh Finite Element Library.
2// Copyright (C) 2002-2026 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
37namespace libMesh
38{
39// Forward declarations
40template <typename T> class TypeTensor;
41template <typename T> class VectorValue;
42template <typename T> class TensorValue;
43}
44
45namespace std
46{
47template <typename T>
48auto norm(const libMesh::TypeVector<T> & vector) -> decltype(std::norm(T()));
49}
50
51namespace libMesh
52{
61template <typename T>
62class TypeVector
63{
64 template <typename T2>
65 friend class TypeVector;
66
67 friend class TypeTensor<T>;
68
69public:
74 TypeVector ();
75
76protected:
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 std::enable_if<ScalarTraits<Scalar1>::value,
92 const Scalar1>::type & x,
93 typename
94 std::enable_if<ScalarTraits<Scalar2>::value,
95 const Scalar2>::type & y=0,
96 typename
97 std::enable_if<ScalarTraits<Scalar3>::value,
98 const Scalar3>::type & z=0);
99
106 template <typename Scalar>
107 TypeVector (const Scalar & x,
108 typename
109 std::enable_if<ScalarTraits<Scalar>::value,
110 const Scalar>::type * sfinae = nullptr);
111
112public:
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 std::enable_if<
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 std::enable_if<
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 std::enable_if<
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;
302
307 auto norm_sq() const;
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
412protected:
413
417 T _coords[LIBMESH_DIM];
418};
419
420
421
422
423
424//------------------------------------------------------
425// Inline functions
426
427template <typename T>
428inline
429TypeVector<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
444template <typename T>
445inline
446TypeVector<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
468template <typename T>
469template <typename Scalar1, typename Scalar2, typename Scalar3>
470inline
471TypeVector<T>::TypeVector (typename
472 std::enable_if<ScalarTraits<Scalar1>::value,
473 const Scalar1>::type & x,
474 typename
475 std::enable_if<ScalarTraits<Scalar2>::value,
476 const Scalar2>::type & y,
477 typename
478 std::enable_if<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
498template <typename T>
499template <typename Scalar>
500inline
501TypeVector<T>::TypeVector (const Scalar & x,
502 typename
503 std::enable_if<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
519template <typename T>
520template <typename T2>
521inline
522TypeVector<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
531template <typename T>
532template <typename T2>
533inline
534void 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
542template <typename T>
543inline
544const 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
553template <typename T>
554inline
555T & TypeVector<T>::operator () (const unsigned int i)
556{
557 libmesh_assert_less (i, LIBMESH_DIM);
558
559 return _coords[i];
560}
561
562
563
564template <typename T>
565template <typename T2>
566inline
567TypeVector<typename CompareTypes<T, T2>::supertype>
568TypeVector<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
590template <typename T>
591template <typename T2>
592inline
593const TypeVector<T> & TypeVector<T>::operator += (const TypeVector<T2> & p)
594{
595 this->add (p);
596
597 return *this;
598}
599
600
601
602template <typename T>
603template <typename T2>
604inline
605void 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
626template <typename T>
627template <typename T2>
628inline
629void 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
650template <typename T>
651template <typename T2>
652inline
653TypeVector<typename CompareTypes<T, T2>::supertype>
654TypeVector<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
677template <typename T>
678template <typename T2>
679inline
680const TypeVector<T> & TypeVector<T>::operator -= (const TypeVector<T2> & p)
681{
682 this->subtract (p);
683
684 return *this;
685}
686
687
688
689template <typename T>
690template <typename T2>
691inline
692void 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
700template <typename T>
701template <typename T2>
702inline
703void 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
711template <typename T>
712inline
713TypeVector<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
735template <typename T>
736template <typename Scalar>
737inline
738typename std::enable_if<
739 ScalarTraits<Scalar>::value,
740 TypeVector<typename CompareTypes<T, Scalar>::supertype>>::type
741TypeVector<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
763template <typename T, typename Scalar>
764inline
765typename std::enable_if<
766 ScalarTraits<Scalar>::value,
767 TypeVector<typename CompareTypes<T, Scalar>::supertype>>::type
768operator * (const Scalar & factor,
769 const TypeVector<T> & v)
770{
771 return v * factor;
772}
773
774
775
776template <typename T>
777inline
778const 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
800template <typename T>
801template <typename Scalar>
802inline
803typename std::enable_if<
804 ScalarTraits<Scalar>::value,
805 TypeVector<typename CompareTypes<T, Scalar>::supertype>>::type
806TypeVector<T>::operator / (const Scalar & factor) const
807{
808 typedef typename CompareTypes<T, Scalar>::supertype TS;
809
810 libmesh_assert_not_equal_to (static_cast<TS>(factor),
811 static_cast<TS>(0.));
812
813#if LIBMESH_DIM == 1
814 return TypeVector<TS>(_coords[0]/factor);
815#endif
816
817#if LIBMESH_DIM == 2
818 return TypeVector<TS>(_coords[0]/factor,
819 _coords[1]/factor);
820#endif
821
822#if LIBMESH_DIM == 3
823 return TypeVector<TS>(_coords[0]/factor,
824 _coords[1]/factor,
825 _coords[2]/factor);
826#endif
827
828}
829
830
831
832
833template <typename T>
834inline
835const TypeVector<T> &
836TypeVector<T>::operator /= (const T & factor)
837{
838 libmesh_assert_not_equal_to (factor, static_cast<T>(0.));
839
840 for (unsigned int i=0; i<LIBMESH_DIM; i++)
841 _coords[i] /= factor;
842
843 return *this;
844}
845
846
847
848
849template <typename T>
850template <typename T2>
851inline
852typename CompareTypes<T, T2>::supertype
853TypeVector<T>::operator * (const TypeVector<T2> & p) const
854{
855#if LIBMESH_DIM == 1
856 return _coords[0]*p._coords[0];
857#endif
858
859#if LIBMESH_DIM == 2
860 return (_coords[0]*p._coords[0] +
861 _coords[1]*p._coords[1]);
862#endif
863
864#if LIBMESH_DIM == 3
865 return (_coords[0]*p(0) +
866 _coords[1]*p(1) +
867 _coords[2]*p(2));
868#endif
869}
870
871template <typename T>
872template <typename T2>
873inline
874typename CompareTypes<T, T2>::supertype
875TypeVector<T>::contract(const TypeVector<T2> & p) const
876{
877 return (*this)*(p);
878}
879
880
881
882template <typename T>
883template <typename T2>
884TypeVector<typename CompareTypes<T, T2>::supertype>
885TypeVector<T>::cross(const TypeVector<T2> & p) const
886{
887 typedef typename CompareTypes<T, T2>::supertype TS;
888 libmesh_assert_equal_to (LIBMESH_DIM, 3);
889
890 // | i j k |
891 // |(*this)(0) (*this)(1) (*this)(2)|
892 // | p(0) p(1) p(2) |
893
894#if LIBMESH_DIM == 3
895 return TypeVector<TS>( _coords[1]*p._coords[2] - _coords[2]*p._coords[1],
896 -_coords[0]*p._coords[2] + _coords[2]*p._coords[0],
897 _coords[0]*p._coords[1] - _coords[1]*p._coords[0]);
898#else
899 libmesh_ignore(p);
900 return TypeVector<TS>(0);
901#endif
902}
903
904
905
906template <typename T>
907inline
908auto TypeVector<T>::norm() const
909{
910 using std::sqrt;
911 return sqrt(this->norm_sq());
912}
913
914
915
916template <typename T>
917inline
918void TypeVector<T>::zero()
919{
920 for (unsigned int i=0; i<LIBMESH_DIM; i++)
921 _coords[i] = 0.;
922}
923
924
925
926template <typename T>
927inline
928auto TypeVector<T>::norm_sq() const
929{
930#if LIBMESH_DIM == 1
931 return (TensorTools::norm_sq(_coords[0]));
932#endif
933
934#if LIBMESH_DIM == 2
935 return (TensorTools::norm_sq(_coords[0]) +
936 TensorTools::norm_sq(_coords[1]));
937#endif
938
939#if LIBMESH_DIM == 3
940 return (TensorTools::norm_sq(_coords[0]) +
941 TensorTools::norm_sq(_coords[1]) +
942 TensorTools::norm_sq(_coords[2]));
943#endif
944}
945
946
947template <typename T>
948inline
949bool TypeVector<T>::is_zero() const
950{
951 for (const auto & val : _coords)
952 if (val != T(0))
953 return false;
954 return true;
955}
956
957template <>
958auto
959TypeVector<bool>::l1_norm() const;
960
961template <typename T>
962auto
963TypeVector<T>::l1_norm() const
964{
965 using std::abs;
966 decltype(abs(T())) ret{};
967 for (const auto i : make_range(libmesh_dim))
968 ret += abs(_coords[i]);
969
970 return ret;
971}
972
973template <typename T>
974inline
975bool TypeVector<T>::absolute_fuzzy_equals(const TypeVector<T> & rhs, Real tol) const
976{
977 return libMesh::absolute_fuzzy_equals(*this, rhs, tol);
978}
979
980
981
982template <typename T>
983inline
984bool TypeVector<T>::relative_fuzzy_equals(const TypeVector<T> & rhs, Real tol) const
985{
986 return libMesh::relative_fuzzy_equals(*this, rhs, tol);
987}
988
989
990
991template <typename T>
992inline
993bool TypeVector<T>::operator == (const TypeVector<T> & rhs) const
994{
995#if LIBMESH_DIM == 1
996 return (_coords[0] == rhs._coords[0]);
997#endif
998
999#if LIBMESH_DIM == 2
1000 return (_coords[0] == rhs._coords[0] &&
1001 _coords[1] == rhs._coords[1]);
1002#endif
1003
1004#if LIBMESH_DIM == 3
1005 return (_coords[0] == rhs._coords[0] &&
1006 _coords[1] == rhs._coords[1] &&
1007 _coords[2] == rhs._coords[2]);
1008#endif
1009}
1010
1011
1012
1013template <typename T>
1014inline
1015bool TypeVector<T>::operator != (const TypeVector<T> & rhs) const
1016{
1017 return (!(*this == rhs));
1018}
1019
1020
1021//------------------------------------------------------
1022// Non-member functions on TypeVectors
1023
1024// Compute a * (b.cross(c)) without creating a temporary
1025// for the cross product. Equivalent to the determinant
1026// of the 3x3 tensor:
1027// [a0, a1, a2]
1028// [b0, b1, b2]
1029// [c0, c1, c2]
1030template <typename T>
1031inline
1032T triple_product(const TypeVector<T> & a,
1033 const TypeVector<T> & b,
1034 const TypeVector<T> & c)
1035{
1036#if LIBMESH_DIM == 3
1037 return
1038 a(0)*(b(1)*c(2) - b(2)*c(1)) -
1039 a(1)*(b(0)*c(2) - b(2)*c(0)) +
1040 a(2)*(b(0)*c(1) - b(1)*c(0));
1041#else
1042 libmesh_ignore(a, b, c);
1043 return 0;
1044#endif
1045}
1046
1047
1048// compute the solid angle subtended by the tetrahedral vertex 0, as
1049// defined by vectors v01, v02, and v03. The solid angle is defined
1050// to be positive if the vectors are obey the right-hand rule, or
1051// negative for a left-hand orientation.
1052template <typename T>
1053inline
1054T solid_angle(const TypeVector<T> & v01,
1055 const TypeVector<T> & v02,
1056 const TypeVector<T> & v03)
1057{
1058 using std::atan;
1059
1060 const Real norm01 = v01.norm(),
1061 norm02 = v02.norm(),
1062 norm03 = v03.norm();
1063 const T tan_half_angle =
1064 triple_product(v01, v02, v03) /
1065 ((v01*v02)*norm03 + (v01*v03)*norm02 + (v02*v03)*norm01 +
1066 norm01*norm02*norm03);
1067
1068 return Real(2)*atan(tan_half_angle);
1069}
1070
1071
1072
1073// Compute the (possibly 3D) circumcenter of three non-collinear points
1074//
1075// If perfectly collinear points are provided, this algorithm computes
1076// NaN. If near-collinear points are provided, the algorithm may be
1077// numerically unstable.
1078template <typename T>
1079inline
1080TypeVector<T> circumcenter(const TypeVector<T> & p0,
1081 const TypeVector<T> & p1,
1082 const TypeVector<T> & p2)
1083{
1084 // smallest subchords in Edge3 order, but this works for any order
1085 const TypeVector<T> e02 = p2 - p0;
1086 const TypeVector<T> e21 = p1 - p2;
1087
1088 const TypeVector<T> scaled_vec = e02.norm_sq()*e21 + e21.norm_sq()*e02;
1089#if LIBMESH_DIM == 3
1090 const TypeVector<T> numerator =
1091 (e02.cross(e21)).cross(scaled_vec);
1092#else
1093 const T e02_cross_e21_z = e02(0)*e21(1)-e02(1)*e21(0);
1094 const TypeVector<T> numerator {-e02_cross_e21_z*scaled_vec(1),
1095 e02_cross_e21_z*scaled_vec(0)};
1096#endif
1097 const TypeVector<T> e2C =
1098 numerator*T(0.5)/cross_norm_sq(e02,e21);
1099
1100 return p2 + e2C;
1101}
1102
1103
1104
1109template <typename T>
1110inline
1111T cross_norm_sq(const TypeVector<T> & b,
1112 const TypeVector<T> & c)
1113{
1114 T z = b(0)*c(1) - b(1)*c(0);
1115
1116#if LIBMESH_DIM == 3
1117 T x = b(1)*c(2) - b(2)*c(1),
1118 y = b(0)*c(2) - b(2)*c(0);
1119 return x*x + y*y + z*z;
1120#else
1121 return z*z;
1122#endif
1123}
1124
1125
1126
1130template <typename T>
1131inline
1132T cross_norm(const TypeVector<T> & b,
1133 const TypeVector<T> & c)
1134{
1135 using std::sqrt;
1136 return sqrt(cross_norm_sq(b,c));
1137}
1138
1139template <typename T>
1140inline
1141TypeVector<T> TypeVector<T>::unit() const
1142{
1143
1144 auto && length = norm();
1145
1146 libmesh_assert_not_equal_to (length, static_cast<Real>(0.));
1147
1148#if LIBMESH_DIM == 1
1149 return TypeVector<T>(_coords[0]/length);
1150#endif
1151
1152#if LIBMESH_DIM == 2
1153 return TypeVector<T>(_coords[0]/length,
1154 _coords[1]/length);
1155#endif
1156
1157#if LIBMESH_DIM == 3
1158 return TypeVector<T>(_coords[0]/length,
1159 _coords[1]/length,
1160 _coords[2]/length);
1161#endif
1162
1163}
1164
1165template <typename T>
1166void TypeVector<T>::print(std::ostream & os) const
1167{
1168#if LIBMESH_DIM == 1
1169
1170 os << "x=" << (*this)(0);
1171
1172#endif
1173#if LIBMESH_DIM == 2
1174
1175 os << "(x,y)=("
1176 << std::setw(8) << (*this)(0) << ", "
1177 << std::setw(8) << (*this)(1) << ")";
1178
1179#endif
1180#if LIBMESH_DIM == 3
1181
1182 os << "(x,y,z)=("
1183 << std::setw(8) << (*this)(0) << ", "
1184 << std::setw(8) << (*this)(1) << ", "
1185 << std::setw(8) << (*this)(2) << ")";
1186#endif
1187}
1188
1189
1190// A vector is finite iff every component is
1191template <typename T>
1192bool isfinite (const TypeVector<T> & var)
1193{
1194 using std::isfinite;
1195 using libMesh::isfinite; // for T==complex
1196 for (auto i : make_range(LIBMESH_DIM))
1197 if (!isfinite(var(i)))
1198 return false;
1199 return true;
1200}
1201
1202
1203// A vector is infinite iff some component is infinite but no
1204// component is NaN.
1205//
1206// This is arguably inconsistent with our std::complex overload (and
1207// the C99 Annex G recommendations for _Complex, and C++ std::complex
1208// arithmetic), which treats mixed (inf,NaN) pairs as infinite, but
1209// this is probably safer for users.
1210template <typename T>
1211bool isinf (const TypeVector<T> & var)
1212{
1213 using std::isinf;
1214 using libMesh::isinf; // for T==complex
1215 using std::isnan;
1216 using libMesh::isnan;
1217 bool has_inf = false;
1218 for (auto i : make_range(LIBMESH_DIM))
1219 {
1220 // NaN anywhere makes us NaN, not inf
1221 if (isnan(var(i)))
1222 return false;
1223 has_inf = has_inf || isinf(var(i));
1224 }
1225 return has_inf;
1226}
1227
1228
1229// A vector is NaN iff some component is NaN
1230//
1231// This is arguably inconsistent with our std::complex overload (and
1232// the C99 Annex G recommendations for _Complex, and C++ std::complex
1233// arithmetic), which treats mixed (inf,NaN) pairs as infinite, but
1234// this is probably safer for users.
1235template <typename T>
1236bool isnan (const TypeVector<T> & var)
1237{
1238 using std::isnan;
1239 using libMesh::isnan; // for T==complex
1240 for (auto i : make_range(LIBMESH_DIM))
1241 if (isnan(var(i)))
1242 return true;
1243 return false;
1244}
1245
1246
1247template <typename T>
1252
1253template <typename T, typename T2>
1258
1259template <typename T, typename T2, typename std::enable_if<ScalarTraits<T>::value, int>::type = 0>
1260TypeVector<typename CompareTypes<T, T2>::supertype>
1261outer_product(const T & a, const TypeVector<T2> & b)
1262{
1263 TypeVector<typename CompareTypes<T, T2>::supertype> ret;
1264 for (unsigned int i = 0; i < LIBMESH_DIM; i++)
1265 ret(i) = a * libmesh_conj(b(i));
1266
1267 return ret;
1268}
1269
1270template <typename T, typename T2, typename std::enable_if<ScalarTraits<T2>::value, int>::type = 0>
1271TypeVector<typename CompareTypes<T, T2>::supertype>
1272outer_product(const TypeVector<T> & a, const T2 & b)
1273{
1274 TypeVector<typename CompareTypes<T, T2>::supertype> ret;
1275 const auto conj_b = libmesh_conj(b);
1276 for (unsigned int i = 0; i < LIBMESH_DIM; i++)
1277 ret(i) = a(i) * conj_b;
1278
1279 return ret;
1280}
1281
1282template <typename T>
1283auto
1284l1_norm(const TypeVector<T> & var)
1285{
1286 return var.l1_norm();
1287}
1288
1289template <typename T, typename T2>
1290auto
1291l1_norm_diff(const TypeVector<T> & vec1, const TypeVector<T2> & vec2)
1292{
1293 return l1_norm(vec1 - vec2);
1294}
1295
1296} // namespace libMesh
1297
1298namespace std
1299{
1300template <typename T>
1301auto norm(const libMesh::TypeVector<T> & vector) -> decltype(std::norm(T()))
1302{
1303 // Yea I agree it's dumb that the standard returns the square of the Euclidean norm
1304 return vector.norm_sq();
1305}
1306} // namespace std
1307
1308#ifdef LIBMESH_HAVE_METAPHYSICL
1309namespace MetaPhysicL
1310{
1311template <typename T>
1312struct RawType<libMesh::TypeVector<T>>
1313{
1314 typedef libMesh::TypeVector<typename RawType<T>::value_type> value_type;
1315
1316 static value_type value (const libMesh::TypeVector<T> & in)
1317 {
1318 value_type ret;
1319 for (unsigned int i = 0; i < LIBMESH_DIM; ++i)
1320 ret(i) = raw_value(in(i));
1321
1322 return ret;
1323 }
1324};
1325
1326template <typename T, typename U>
1327struct ReplaceAlgebraicType<libMesh::TypeVector<T>, U>
1328{
1329 typedef U type;
1330};
1331}
1332#endif
1333
1334#endif // LIBMESH_TYPE_VECTOR_H
libMesh::TypeTensor
This class defines a tensor in LIBMESH_DIM dimensional space of type T.
Definition type_tensor.h:53
libMesh::TypeVector
This class defines a vector in LIBMESH_DIM dimensional space of type T.
Definition type_vector.h:63
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
libMesh::TypeVector::slice
T & slice(const unsigned int i)
Definition type_vector.h:166
libMesh::TypeVector::operator>
bool operator>(const TypeVector< T > &rhs) const
Definition type_vector.C:95
libMesh::TypeVector::operator*=
const TypeVector< T > & operator*=(const T &)
Multiply this vector by a scalar value.
Definition type_vector.h:778
libMesh::TypeVector::_coords
T _coords[LIBMESH_DIM]
The coordinates of the TypeVector.
Definition type_vector.h:417
libMesh::TypeVector::operator=
std::enable_if< ScalarTraits< Scalar >::value, TypeVector & >::type operator=(const Scalar &libmesh_dbg_var(p))
Assignment-from-scalar operator.
Definition type_vector.h:153
libMesh::TypeVector::TypeVector
TypeVector(const T &x, const T &y=0, const T &z=0)
Constructor-from-T.
Definition type_vector.h:446
libMesh::TypeVector::absolute_fuzzy_equals
bool absolute_fuzzy_equals(const TypeVector< T > &rhs, Real tol=TOLERANCE) const
Definition type_vector.h:975
libMesh::TypeVector::assign
void assign(const TypeVector< T2 > &)
Assign to this vector without creating a temporary.
Definition type_vector.h:534
libMesh::TypeVector::operator>=
bool operator>=(const TypeVector< T > &rhs) const
Definition type_vector.C:109
libMesh::TypeVector::TypeVector
TypeVector(const Scalar &x, typename std::enable_if< ScalarTraits< Scalar >::value, const Scalar >::type *sfinae=nullptr)
Constructor-from-scalar.
Definition type_vector.h:501
libMesh::TypeVector::norm_sq
auto norm_sq() const
Definition type_vector.h:928
libMesh::TypeVector::cross
TypeVector< typename CompareTypes< T, T2 >::supertype > cross(const TypeVector< T2 > &v) const
Definition type_vector.h:885
libMesh::TypeVector::zero
void zero()
Set all entries of the vector to 0.
Definition type_vector.h:918
libMesh::TypeVector::TypeVector
TypeVector(const TypeVector &p)=default
Copy-constructor for the trivial case.
libMesh::TypeVector::operator+
TypeVector< typename CompareTypes< T, T2 >::supertype > operator+(const TypeVector< T2 > &) const
Add two vectors.
Definition type_vector.h:568
libMesh::TypeVector::slice
const T & slice(const unsigned int i) const
Definition type_vector.h:160
libMesh::TypeVector::operator!=
bool operator!=(const TypeVector< T > &rhs) const
Definition type_vector.h:1015
libMesh::TypeVector::operator-=
const TypeVector< T > & operator-=(const TypeVector< T2 > &)
Subtract from this vector.
Definition type_vector.h:680
libMesh::TypeVector::~TypeVector
~TypeVector()=default
Destructor.
libMesh::TypeVector::operator<=
bool operator<=(const TypeVector< T > &rhs) const
Definition type_vector.C:80
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::value_type
T value_type
Helper typedef for C++98 generic programming.
Definition type_vector.h:117
libMesh::TypeVector::is_zero
bool is_zero() const
Definition type_vector.h:949
libMesh::TypeVector::index_type
unsigned int index_type
Helper typedef for generic index programming.
Definition type_vector.h:122
libMesh::TypeVector::relative_fuzzy_equals
bool relative_fuzzy_equals(const TypeVector< T > &rhs, Real tol=TOLERANCE) const
Definition type_vector.h:984
libMesh::TypeVector::operator<<
friend std::ostream & operator<<(std::ostream &os, const TypeVector< T > &t)
Formatted print as above but supports the syntax:
Definition type_vector.h:397
libMesh::TypeVector::operator*
std::enable_if< ScalarTraits< Scalar >::value, TypeVector< typenameCompareTypes< T, Scalar >::supertype > >::type operator*(const Scalar &) const
Multiply this vector by a scalar value.
Definition type_vector.h:741
libMesh::TypeVector::TypeVector
TypeVector(typename std::enable_if< ScalarTraits< Scalar1 >::value, const Scalar1 >::type &x, typename std::enable_if< ScalarTraits< Scalar2 >::value, const Scalar2 >::type &y=0, typename std::enable_if< ScalarTraits< Scalar3 >::value, const Scalar3 >::type &z=0)
Constructor-from-scalars.
Definition type_vector.h:471
libMesh::TypeVector::contract
CompareTypes< T, T2 >::supertype contract(const TypeVector< T2 > &) const
Definition type_vector.h:875
libMesh::TypeVector::operator/=
const TypeVector< T > & operator/=(const T &)
Divide each entry of this vector by scalar value.
Definition type_vector.h:836
libMesh::TypeVector::TypeVector
TypeVector()
Empty constructor.
Definition type_vector.h:429
libMesh::TypeVector::add
void add(const TypeVector< T2 > &)
Add to this vector without creating a temporary.
Definition type_vector.h:605
libMesh::TypeVector::operator<
bool operator<(const TypeVector< T > &rhs) const
Definition type_vector.C:66
libMesh::TypeVector::subtract
void subtract(const TypeVector< T2 > &)
Subtract from this vector without creating a temporary.
Definition type_vector.h:692
libMesh::TypeVector::operator-
TypeVector< T > operator-() const
Definition type_vector.h:713
libMesh::TypeVector::norm
auto norm() const
Definition type_vector.h:908
libMesh::TypeVector::operator()
const T & operator()(const unsigned int i) const
Definition type_vector.h:544
libMesh::TypeVector::operator+=
const TypeVector< T > & operator+=(const TypeVector< T2 > &)
Add to this vector.
Definition type_vector.h:593
libMesh::TypeVector::TypeVector
friend class TypeVector
Definition type_vector.h:65
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::TypeVector::operator/
std::enable_if< ScalarTraits< Scalar >::value, TypeVector< typenameCompareTypes< T, Scalar >::supertype > >::type operator/(const Scalar &) const
Divide each entry of this vector by scalar value.
Definition type_vector.h:806
libMesh::TypeVector::TypeVector
TypeVector(const TypeVector< T2 > &p)
Copy-constructor.
Definition type_vector.h:522
libMesh::TypeVector::l1_norm
auto l1_norm() const
Definition type_vector.h:963
libMesh::TypeVector::unit
TypeVector< T > unit() const
Definition type_vector.h:1141
libMesh::TypeVector::operator==
bool operator==(const TypeVector< T > &rhs) const
Definition type_vector.h:993
libMesh::TypeVector::print
void print(std::ostream &os=libMesh::out) const
Formatted print, by default to libMesh::out.
Definition type_vector.h:1166
b
static const Real b
Definition disjoint_neighbor_test.C:24
libMesh::TensorTools::norm_sq
auto norm_sq(const T &a)
Definition tensor_tools.h:104
libMesh
The libMesh namespace provides an interface to certain functionality in the library.
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
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::l1_norm
auto l1_norm(const NumericVector< T > &vec)
Definition numeric_vector.h:1144
libMesh::libmesh_dim
static constexpr std::size_t libmesh_dim
Definition libmesh_common.h:720
libMesh::libmesh_ignore
void libmesh_ignore(const Args &...)
Definition libmesh_common.h:569
libMesh::triple_product
T triple_product(const TypeVector< T > &a, const TypeVector< T > &b, const TypeVector< T > &c)
Definition type_vector.h:1032
libMesh::isfinite
bool isfinite(std::complex< T > a)
Definition libmesh_common.h:231
libMesh::isnan
bool isnan(std::complex< T > a)
Definition libmesh_common.h:243
libMesh::out
OStreamProxy out
libMesh::operator*
std::enable_if< ScalarTraits< Scalar >::value, TypeNTensor< N, typenameCompareTypes< Scalar, T >::supertype > >::type operator*(const Scalar &, const TypeNTensor< N, T > &)
Definition type_n_tensor.h:323
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:1132
libMesh::outer_product
TypeTensor< typename CompareTypes< T, T2 >::supertype > outer_product(const TypeVector< T > &a, const TypeVector< T2 > &b)
Definition type_tensor.h:1504
libMesh::TOLERANCE
static constexpr Real TOLERANCE
Definition libmesh_common.h:151
libMesh::isinf
bool isinf(std::complex< T > a)
Definition libmesh_common.h:237
libMesh::solid_angle
T solid_angle(const TypeVector< T > &v01, const TypeVector< T > &v02, const TypeVector< T > &v03)
Definition type_vector.h:1054
libMesh::l1_norm_diff
auto l1_norm_diff(const NumericVector< T > &vec1, const NumericVector< T > &vec2)
Definition numeric_vector.h:1151
libMesh::libmesh_conj
T libmesh_conj(T a)
Definition libmesh_common.h:191
libMesh::circumcenter
TypeVector< T > circumcenter(const TypeVector< T > &p0, const TypeVector< T > &p1, const TypeVector< T > &p2)
Definition type_vector.h:1080
libMesh::Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Definition libmesh_common.h:144
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:176
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)....
Definition type_vector.h:1111
std::norm
auto norm(const libMesh::TypeVector< T > &vector) -> decltype(std::norm(T()))
Definition type_vector.h:1301
MetaPhysicL::RawType< libMesh::TypeVector< T > >::value_type
libMesh::TypeVector< typename RawType< T >::value_type > value_type
Definition type_vector.h:1314
MetaPhysicL::RawType< libMesh::TypeVector< T > >::value
static value_type value(const libMesh::TypeVector< T > &in)
Definition type_vector.h:1316
libMesh::CompareTypes< TypeVector< T >, TypeVector< T2 > >::supertype
TypeVector< typename CompareTypes< T, T2 >::supertype > supertype
Definition type_vector.h:1256
Generated on Mon Aug 24 2026 19:00:33 for libMesh by doxygen 1.9.8