type_vector.h

Go to the documentation of this file.

// The libMesh Finite Element Library.
// Copyright (C) 2002-2026 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner

// This library is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.

// This library is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// Lesser General Public License for more details.

// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA



#ifndef LIBMESH_TYPE_VECTOR_H
#define LIBMESH_TYPE_VECTOR_H

// Local includes
#include "libmesh/libmesh_common.h"
#include "libmesh/compare_types.h"
#include "libmesh/tensor_tools.h"
#include "libmesh/int_range.h"
#include "libmesh/fuzzy_equals.h"

// C++ includes
#include <cstdlib> // *must* precede <cmath> for proper std:abs() on PGI, Sun Studio CC
#include <cmath>
#include <complex>
#include <iostream>
#include <iomanip> // for std::setw, std::setiosflags

namespace libMesh
{
// Forward declarations
template <typename T> class TypeTensor;
template <typename T> class VectorValue;
template <typename T> class TensorValue;
}

namespace std
{
template <typename T>
auto norm(const libMesh::TypeVector<T> & vector) -> decltype(std::norm(T()));
}

namespace libMesh
{
template <typename T>
class TypeVector
{
  template <typename T2>
  friend class TypeVector;

  friend class TypeTensor<T>;

public:
  TypeVector  ();

protected:
  TypeVector (const T & x,
              const T & y=0,
              const T & z=0);

  template <typename Scalar1, typename Scalar2, typename Scalar3>
  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);

  template <typename Scalar>
  TypeVector (const Scalar & x,
              typename
              std::enable_if<ScalarTraits<Scalar>::value,
              const Scalar>::type * sfinae = nullptr);

public:

  typedef T value_type;

  typedef unsigned int index_type;

  template <typename T2>
  TypeVector (const TypeVector<T2> & p);

  TypeVector (const TypeVector & p) = default;

  ~TypeVector () = default;

  template <typename T2>
  void assign (const TypeVector<T2> &);

  template <typename Scalar>
  typename std::enable_if<
    ScalarTraits<Scalar>::value,
    TypeVector &>::type
  operator = (const Scalar & libmesh_dbg_var(p))
  { libmesh_assert_equal_to (p, Scalar(0)); this->zero(); return *this; }

  const T & operator () (const unsigned int i) const;
  const T & slice (const unsigned int i) const { return (*this)(i); }

  T & operator () (const unsigned int i);
  T & slice (const unsigned int i) { return (*this)(i); }

  template <typename T2>
  TypeVector<typename CompareTypes<T, T2>::supertype>
  operator + (const TypeVector<T2> &) const;

  template <typename T2>
  const TypeVector<T> & operator += (const TypeVector<T2> &);

  template <typename T2>
  void add (const TypeVector<T2> &);

  template <typename T2>
  void add_scaled (const TypeVector<T2> &, const T &);

  template <typename T2>
  TypeVector<typename CompareTypes<T, T2>::supertype>
  operator - (const TypeVector<T2> &) const;

  template <typename T2>
  const TypeVector<T> & operator -= (const TypeVector<T2> &);

  template <typename T2>
  void subtract (const TypeVector<T2> &);

  template <typename T2>
  void subtract_scaled (const TypeVector<T2> &, const T &);

  TypeVector<T> operator - () const;

  template <typename Scalar>
  typename std::enable_if<
    ScalarTraits<Scalar>::value,
    TypeVector<typename CompareTypes<T, Scalar>::supertype>>::type
  operator * (const Scalar &) const;

  const TypeVector<T> & operator *= (const T &);

  template <typename Scalar>
  typename std::enable_if<
    ScalarTraits<Scalar>::value,
    TypeVector<typename CompareTypes<T, Scalar>::supertype>>::type
  operator / (const Scalar &) const;

  const TypeVector<T> & operator /= (const T &);

  template <typename T2>
  typename CompareTypes<T, T2>::supertype
  operator * (const TypeVector<T2> &) const;

  template <typename T2>
  typename CompareTypes<T, T2>::supertype
  contract (const TypeVector<T2> &) const;

  template <typename T2>
  TypeVector<typename CompareTypes<T, T2>::supertype>
  cross(const TypeVector<T2> & v) const;

  TypeVector<T> unit() const;

  auto norm() const;

  auto norm_sq() const;

  auto l1_norm() const;

  bool is_zero() const;

  void zero();

  bool relative_fuzzy_equals(const TypeVector<T> & rhs,
                             Real tol = TOLERANCE) const;

  bool absolute_fuzzy_equals(const TypeVector<T> & rhs,
                             Real tol = TOLERANCE) const;

  bool operator == (const TypeVector<T> & rhs) const;

  bool operator != (const TypeVector<T> & rhs) const;

  bool operator < (const TypeVector<T> & rhs) const;

  bool operator <= (const TypeVector<T> & rhs) const;

  bool operator > (const TypeVector<T> & rhs) const;

  bool operator >= (const TypeVector<T> & rhs) const;

  void print(std::ostream & os = libMesh::out) const;

  friend std::ostream & operator << (std::ostream & os, const TypeVector<T> & t)
  {
    t.print(os);
    return os;
  }

  void write_unformatted (std::ostream & out_stream,
                          const bool newline = true) const;

protected:

  T _coords[LIBMESH_DIM];
};





//------------------------------------------------------
// Inline functions

template <typename T>
inline
TypeVector<T>::TypeVector ()
{
  _coords[0] = {};

#if LIBMESH_DIM > 1
  _coords[1] = {};
#endif

#if LIBMESH_DIM > 2
  _coords[2] = {};
#endif
}



template <typename T>
inline
TypeVector<T>::TypeVector (const T & x,
                           const T & y,
                           const T & z)
{
  _coords[0] = x;

#if LIBMESH_DIM > 1
  _coords[1] = y;
#else
  libmesh_ignore(y);
  libmesh_assert_equal_to (y, 0);
#endif

#if LIBMESH_DIM > 2
  _coords[2] = z;
#else
  libmesh_ignore(z);
  libmesh_assert_equal_to (z, 0);
#endif
}


template <typename T>
template <typename Scalar1, typename Scalar2, typename Scalar3>
inline
TypeVector<T>::TypeVector (typename
                           std::enable_if<ScalarTraits<Scalar1>::value,
                           const Scalar1>::type & x,
                           typename
                           std::enable_if<ScalarTraits<Scalar2>::value,
                           const Scalar2>::type & y,
                           typename
                           std::enable_if<ScalarTraits<Scalar3>::value,
                           const Scalar3>::type & z)
{
  _coords[0] = x;

#if LIBMESH_DIM > 1
  _coords[1] = y;
#else
  libmesh_assert_equal_to (y, 0);
#endif

#if LIBMESH_DIM > 2
  _coords[2] = z;
#else
  libmesh_assert_equal_to (z, 0);
#endif
}



template <typename T>
template <typename Scalar>
inline
TypeVector<T>::TypeVector (const Scalar & x,
                           typename
                           std::enable_if<ScalarTraits<Scalar>::value,
                           const Scalar>::type * /*sfinae*/)
{
  _coords[0] = x;

#if LIBMESH_DIM > 1
  _coords[1] = 0;
#endif

#if LIBMESH_DIM > 2
  _coords[2] = 0;
#endif
}



template <typename T>
template <typename T2>
inline
TypeVector<T>::TypeVector (const TypeVector<T2> & p)
{
  // copy the nodes from vector p to me
  for (unsigned int i=0; i<LIBMESH_DIM; i++)
    _coords[i] = p._coords[i];
}



template <typename T>
template <typename T2>
inline
void TypeVector<T>::assign (const TypeVector<T2> & p)
{
  for (unsigned int i=0; i<LIBMESH_DIM; i++)
    _coords[i] = p._coords[i];
}



template <typename T>
inline
const T & TypeVector<T>::operator () (const unsigned int i) const
{
  libmesh_assert_less (i, LIBMESH_DIM);

  return _coords[i];
}



template <typename T>
inline
T & TypeVector<T>::operator () (const unsigned int i)
{
  libmesh_assert_less (i, LIBMESH_DIM);

  return _coords[i];
}



template <typename T>
template <typename T2>
inline
TypeVector<typename CompareTypes<T, T2>::supertype>
TypeVector<T>::operator + (const TypeVector<T2> & p) const
{
  typedef typename CompareTypes<T, T2>::supertype TS;
#if LIBMESH_DIM == 1
  return TypeVector<TS> (_coords[0] + p._coords[0]);
#endif

#if LIBMESH_DIM == 2
  return TypeVector<TS> (_coords[0] + p._coords[0],
                         _coords[1] + p._coords[1]);
#endif

#if LIBMESH_DIM == 3
  return TypeVector<TS> (_coords[0] + p._coords[0],
                         _coords[1] + p._coords[1],
                         _coords[2] + p._coords[2]);
#endif

}



template <typename T>
template <typename T2>
inline
const TypeVector<T> & TypeVector<T>::operator += (const TypeVector<T2> & p)
{
  this->add (p);

  return *this;
}



template <typename T>
template <typename T2>
inline
void TypeVector<T>::add (const TypeVector<T2> & p)
{
#if LIBMESH_DIM == 1
  _coords[0] += p._coords[0];
#endif

#if LIBMESH_DIM == 2
  _coords[0] += p._coords[0];
  _coords[1] += p._coords[1];
#endif

#if LIBMESH_DIM == 3
  _coords[0] += p._coords[0];
  _coords[1] += p._coords[1];
  _coords[2] += p._coords[2];
#endif

}



template <typename T>
template <typename T2>
inline
void TypeVector<T>::add_scaled (const TypeVector<T2> & p, const T & factor)
{
#if LIBMESH_DIM == 1
  _coords[0] += factor*p(0);
#endif

#if LIBMESH_DIM == 2
  _coords[0] += factor*p(0);
  _coords[1] += factor*p(1);
#endif

#if LIBMESH_DIM == 3
  _coords[0] += factor*p(0);
  _coords[1] += factor*p(1);
  _coords[2] += factor*p(2);
#endif

}



template <typename T>
template <typename T2>
inline
TypeVector<typename CompareTypes<T, T2>::supertype>
TypeVector<T>::operator - (const TypeVector<T2> & p) const
{
  typedef typename CompareTypes<T, T2>::supertype TS;

#if LIBMESH_DIM == 1
  return TypeVector<TS>(_coords[0] - p._coords[0]);
#endif

#if LIBMESH_DIM == 2
  return TypeVector<TS>(_coords[0] - p._coords[0],
                        _coords[1] - p._coords[1]);
#endif

#if LIBMESH_DIM == 3
  return TypeVector<TS>(_coords[0] - p._coords[0],
                        _coords[1] - p._coords[1],
                        _coords[2] - p._coords[2]);
#endif

}



template <typename T>
template <typename T2>
inline
const TypeVector<T> & TypeVector<T>::operator -= (const TypeVector<T2> & p)
{
  this->subtract (p);

  return *this;
}



template <typename T>
template <typename T2>
inline
void TypeVector<T>::subtract (const TypeVector<T2> & p)
{
  for (unsigned int i=0; i<LIBMESH_DIM; i++)
    _coords[i] -= p._coords[i];
}



template <typename T>
template <typename T2>
inline
void TypeVector<T>::subtract_scaled (const TypeVector<T2> & p, const T & factor)
{
  for (unsigned int i=0; i<LIBMESH_DIM; i++)
    _coords[i] -= factor*p(i);
}



template <typename T>
inline
TypeVector<T> TypeVector<T>::operator - () const
{

#if LIBMESH_DIM == 1
  return TypeVector(-_coords[0]);
#endif

#if LIBMESH_DIM == 2
  return TypeVector(-_coords[0],
                    -_coords[1]);
#endif

#if LIBMESH_DIM == 3
  return TypeVector(-_coords[0],
                    -_coords[1],
                    -_coords[2]);
#endif

}



template <typename T>
template <typename Scalar>
inline
typename std::enable_if<
  ScalarTraits<Scalar>::value,
  TypeVector<typename CompareTypes<T, Scalar>::supertype>>::type
TypeVector<T>::operator * (const Scalar & factor) const
{
  typedef typename CompareTypes<T, Scalar>::supertype SuperType;

#if LIBMESH_DIM == 1
  return TypeVector<SuperType>(_coords[0]*factor);
#endif

#if LIBMESH_DIM == 2
  return TypeVector<SuperType>(_coords[0]*factor,
                               _coords[1]*factor);
#endif

#if LIBMESH_DIM == 3
  return TypeVector<SuperType>(_coords[0]*factor,
                               _coords[1]*factor,
                               _coords[2]*factor);
#endif
}



template <typename T, typename Scalar>
inline
typename std::enable_if<
  ScalarTraits<Scalar>::value,
  TypeVector<typename CompareTypes<T, Scalar>::supertype>>::type
operator * (const Scalar & factor,
            const TypeVector<T> & v)
{
  return v * factor;
}



template <typename T>
inline
const TypeVector<T> & TypeVector<T>::operator *= (const T & factor)
{
#if LIBMESH_DIM == 1
  _coords[0] *= factor;
#endif

#if LIBMESH_DIM == 2
  _coords[0] *= factor;
  _coords[1] *= factor;
#endif

#if LIBMESH_DIM == 3
  _coords[0] *= factor;
  _coords[1] *= factor;
  _coords[2] *= factor;
#endif

  return *this;
}



template <typename T>
template <typename Scalar>
inline
typename std::enable_if<
  ScalarTraits<Scalar>::value,
  TypeVector<typename CompareTypes<T, Scalar>::supertype>>::type
TypeVector<T>::operator / (const Scalar & factor) const
{
  typedef typename CompareTypes<T, Scalar>::supertype TS;

  libmesh_assert_not_equal_to (static_cast<TS>(factor),
                               static_cast<TS>(0.));

#if LIBMESH_DIM == 1
  return TypeVector<TS>(_coords[0]/factor);
#endif

#if LIBMESH_DIM == 2
  return TypeVector<TS>(_coords[0]/factor,
                        _coords[1]/factor);
#endif

#if LIBMESH_DIM == 3
  return TypeVector<TS>(_coords[0]/factor,
                        _coords[1]/factor,
                        _coords[2]/factor);
#endif

}




template <typename T>
inline
const TypeVector<T> &
TypeVector<T>::operator /= (const T & factor)
{
  libmesh_assert_not_equal_to (factor, static_cast<T>(0.));

  for (unsigned int i=0; i<LIBMESH_DIM; i++)
    _coords[i] /= factor;

  return *this;
}




template <typename T>
template <typename T2>
inline
typename CompareTypes<T, T2>::supertype
TypeVector<T>::operator * (const TypeVector<T2> & p) const
{
#if LIBMESH_DIM == 1
  return _coords[0]*p._coords[0];
#endif

#if LIBMESH_DIM == 2
  return (_coords[0]*p._coords[0] +
          _coords[1]*p._coords[1]);
#endif

#if LIBMESH_DIM == 3
  return (_coords[0]*p(0) +
          _coords[1]*p(1) +
          _coords[2]*p(2));
#endif
}

template <typename T>
template <typename T2>
inline
typename CompareTypes<T, T2>::supertype
TypeVector<T>::contract(const TypeVector<T2> & p) const
{
  return (*this)*(p);
}



template <typename T>
template <typename T2>
TypeVector<typename CompareTypes<T, T2>::supertype>
TypeVector<T>::cross(const TypeVector<T2> & p) const
{
  typedef typename CompareTypes<T, T2>::supertype TS;
  libmesh_assert_equal_to (LIBMESH_DIM, 3);

  // |     i          j          k    |
  // |(*this)(0) (*this)(1) (*this)(2)|
  // |   p(0)       p(1)       p(2)   |

#if LIBMESH_DIM == 3
  return TypeVector<TS>( _coords[1]*p._coords[2] - _coords[2]*p._coords[1],
                         -_coords[0]*p._coords[2] + _coords[2]*p._coords[0],
                         _coords[0]*p._coords[1] - _coords[1]*p._coords[0]);
#else
  libmesh_ignore(p);
  return TypeVector<TS>(0);
#endif
}



template <typename T>
inline
auto TypeVector<T>::norm() const
{
  using std::sqrt;
  return sqrt(this->norm_sq());
}



template <typename T>
inline
void TypeVector<T>::zero()
{
  for (unsigned int i=0; i<LIBMESH_DIM; i++)
    _coords[i] = 0.;
}



template <typename T>
inline
auto TypeVector<T>::norm_sq() const
{
#if LIBMESH_DIM == 1
  return (TensorTools::norm_sq(_coords[0]));
#endif

#if LIBMESH_DIM == 2
  return (TensorTools::norm_sq(_coords[0]) +
          TensorTools::norm_sq(_coords[1]));
#endif

#if LIBMESH_DIM == 3
  return (TensorTools::norm_sq(_coords[0]) +
          TensorTools::norm_sq(_coords[1]) +
          TensorTools::norm_sq(_coords[2]));
#endif
}


template <typename T>
inline
bool TypeVector<T>::is_zero() const
{
  for (const auto & val : _coords)
    if (val != T(0))
      return false;
  return true;
}

template <>
auto
TypeVector<bool>::l1_norm() const;

template <typename T>
auto
TypeVector<T>::l1_norm() const
{
  using std::abs;
  decltype(abs(T())) ret{};
  for (const auto i : make_range(libmesh_dim))
    ret += abs(_coords[i]);

  return ret;
}

template <typename T>
inline
bool TypeVector<T>::absolute_fuzzy_equals(const TypeVector<T> & rhs, Real tol) const
{
  return libMesh::absolute_fuzzy_equals(*this, rhs, tol);
}



template <typename T>
inline
bool TypeVector<T>::relative_fuzzy_equals(const TypeVector<T> & rhs, Real tol) const
{
  return libMesh::relative_fuzzy_equals(*this, rhs, tol);
}



template <typename T>
inline
bool TypeVector<T>::operator == (const TypeVector<T> & rhs) const
{
#if LIBMESH_DIM == 1
  return (_coords[0] == rhs._coords[0]);
#endif

#if LIBMESH_DIM == 2
  return (_coords[0] == rhs._coords[0] &&
          _coords[1] == rhs._coords[1]);
#endif

#if LIBMESH_DIM == 3
  return (_coords[0] == rhs._coords[0] &&
          _coords[1] == rhs._coords[1] &&
          _coords[2] == rhs._coords[2]);
#endif
}



template <typename T>
inline
bool TypeVector<T>::operator != (const TypeVector<T> & rhs) const
{
  return (!(*this == rhs));
}


//------------------------------------------------------
// Non-member functions on TypeVectors

// Compute a * (b.cross(c)) without creating a temporary
// for the cross product.  Equivalent to the determinant
// of the 3x3 tensor:
// [a0, a1, a2]
// [b0, b1, b2]
// [c0, c1, c2]
template <typename T>
inline
T triple_product(const TypeVector<T> & a,
                 const TypeVector<T> & b,
                 const TypeVector<T> & c)
{
#if LIBMESH_DIM == 3
  return
    a(0)*(b(1)*c(2) - b(2)*c(1)) -
    a(1)*(b(0)*c(2) - b(2)*c(0)) +
    a(2)*(b(0)*c(1) - b(1)*c(0));
#else
  libmesh_ignore(a, b, c);
  return 0;
#endif
}


// compute the solid angle subtended by the tetrahedral vertex 0, as
// defined by vectors v01, v02, and v03.  The solid angle is defined
// to be positive if the vectors are obey the right-hand rule, or
// negative for a left-hand orientation.
template <typename T>
inline
T solid_angle(const TypeVector<T> & v01,
              const TypeVector<T> & v02,
              const TypeVector<T> & v03)
{
  using std::atan;

  const Real norm01 = v01.norm(),
             norm02 = v02.norm(),
             norm03 = v03.norm();
  const T tan_half_angle =
    triple_product(v01, v02, v03) /
    ((v01*v02)*norm03 + (v01*v03)*norm02 + (v02*v03)*norm01 +
     norm01*norm02*norm03);

  return Real(2)*atan(tan_half_angle);
}



// Compute the (possibly 3D) circumcenter of three non-collinear points
//
// If perfectly collinear points are provided, this algorithm computes
// NaN.  If near-collinear points are provided, the algorithm may be
// numerically unstable.
template <typename T>
inline
TypeVector<T> circumcenter(const TypeVector<T> & p0,
                           const TypeVector<T> & p1,
                           const TypeVector<T> & p2)
{
  // smallest subchords in Edge3 order, but this works for any order
  const TypeVector<T> e02 = p2 - p0;
  const TypeVector<T> e21 = p1 - p2;

  const TypeVector<T> scaled_vec = e02.norm_sq()*e21 + e21.norm_sq()*e02;
#if LIBMESH_DIM == 3
  const TypeVector<T> numerator =
    (e02.cross(e21)).cross(scaled_vec);
#else
  const T e02_cross_e21_z = e02(0)*e21(1)-e02(1)*e21(0);
  const TypeVector<T> numerator {-e02_cross_e21_z*scaled_vec(1),
                         e02_cross_e21_z*scaled_vec(0)};
#endif
  const TypeVector<T> e2C =
    numerator*T(0.5)/cross_norm_sq(e02,e21);

  return p2 + e2C;
}



template <typename T>
inline
T cross_norm_sq(const TypeVector<T> & b,
                const TypeVector<T> & c)
{
  T z = b(0)*c(1) - b(1)*c(0);

#if LIBMESH_DIM == 3
  T x = b(1)*c(2) - b(2)*c(1),
    y = b(0)*c(2) - b(2)*c(0);
  return x*x + y*y + z*z;
#else
  return z*z;
#endif
}



template <typename T>
inline
T cross_norm(const TypeVector<T> & b,
             const TypeVector<T> & c)
{
  using std::sqrt;
  return sqrt(cross_norm_sq(b,c));
}

template <typename T>
inline
TypeVector<T> TypeVector<T>::unit() const
{

  auto && length = norm();

  libmesh_assert_not_equal_to (length, static_cast<Real>(0.));

#if LIBMESH_DIM == 1
  return TypeVector<T>(_coords[0]/length);
#endif

#if LIBMESH_DIM == 2
  return TypeVector<T>(_coords[0]/length,
                       _coords[1]/length);
#endif

#if LIBMESH_DIM == 3
  return TypeVector<T>(_coords[0]/length,
                       _coords[1]/length,
                       _coords[2]/length);
#endif

}

template <typename T>
void TypeVector<T>::print(std::ostream & os) const
{
#if LIBMESH_DIM == 1

  os << "x=" << (*this)(0);

#endif
#if LIBMESH_DIM == 2

  os << "(x,y)=("
     << std::setw(8) << (*this)(0) << ", "
     << std::setw(8) << (*this)(1) << ")";

#endif
#if LIBMESH_DIM == 3

  os <<  "(x,y,z)=("
     << std::setw(8) << (*this)(0) << ", "
     << std::setw(8) << (*this)(1) << ", "
     << std::setw(8) << (*this)(2) << ")";
#endif
}

template <typename T>
struct CompareTypes<TypeVector<T>, TypeVector<T>>
{
  typedef TypeVector<T> supertype;
};

template <typename T, typename T2>
struct CompareTypes<TypeVector<T>, TypeVector<T2>>
{
  typedef TypeVector<typename CompareTypes<T,T2>::supertype> supertype;
};

template <typename T, typename T2, typename std::enable_if<ScalarTraits<T>::value, int>::type = 0>
TypeVector<typename CompareTypes<T, T2>::supertype>
outer_product(const T & a, const TypeVector<T2> & b)
{
  TypeVector<typename CompareTypes<T, T2>::supertype> ret;
  for (unsigned int i = 0; i < LIBMESH_DIM; i++)
    ret(i) = a * libmesh_conj(b(i));

  return ret;
}

template <typename T, typename T2, typename std::enable_if<ScalarTraits<T2>::value, int>::type = 0>
TypeVector<typename CompareTypes<T, T2>::supertype>
outer_product(const TypeVector<T> & a, const T2 & b)
{
  TypeVector<typename CompareTypes<T, T2>::supertype> ret;
  const auto conj_b = libmesh_conj(b);
  for (unsigned int i = 0; i < LIBMESH_DIM; i++)
    ret(i) = a(i) * conj_b;

  return ret;
}

template <typename T>
auto
l1_norm(const TypeVector<T> & var)
{
  return var.l1_norm();
}

template <typename T, typename T2>
auto
l1_norm_diff(const TypeVector<T> & vec1, const TypeVector<T2> & vec2)
{
  return l1_norm(vec1 - vec2);
}

} // namespace libMesh

namespace std
{
template <typename T>
auto norm(const libMesh::TypeVector<T> & vector) -> decltype(std::norm(T()))
{
  // Yea I agree it's dumb that the standard returns the square of the Euclidean norm
  return vector.norm_sq();
}
} // namespace std

#ifdef LIBMESH_HAVE_METAPHYSICL
namespace MetaPhysicL
{
template <typename T>
struct RawType<libMesh::TypeVector<T>>
{
  typedef libMesh::TypeVector<typename RawType<T>::value_type> value_type;

  static value_type value (const libMesh::TypeVector<T> & in)
    {
      value_type ret;
      for (unsigned int i = 0; i < LIBMESH_DIM; ++i)
        ret(i) = raw_value(in(i));

      return ret;
    }
};

template <typename T, typename U>
struct ReplaceAlgebraicType<libMesh::TypeVector<T>, U>
{
  typedef U type;
};
}
#endif

#endif // LIBMESH_TYPE_VECTOR_H