type_vector.h

一个用来实现偏微分方程中网格的计算库

  return *this;
}



template <typename T>
template <typename T2>
inline
void TypeVector<T>::add (const TypeVector<T2> &p)
{
#if DIM == 1
  _coords[0] += p._coords[0];
#endif
  
#if DIM == 2
  _coords[0] += p._coords[0];
  _coords[1] += p._coords[1];
#endif
  
#if 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 DIM == 1
  _coords[0] += factor*p(0);
#endif
  
#if DIM == 2
  _coords[0] += factor*p(0);
  _coords[1] += factor*p(1);
#endif
  
#if 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 DIM == 1
  return TypeVector<TS>(_coords[0] - p._coords[0]);
#endif

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

#if 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<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<DIM; i++)
    _coords[i] -= factor*p(i);
}



template <typename T>
inline
TypeVector<T> TypeVector<T>::operator - () const
{
  
#if DIM == 1
  return TypeVector(-_coords[0]);
#endif

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

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



template <typename T>
template <typename Scalar>
inline
typename boostcopy::enable_if_c<
  ScalarTraits<Scalar>::value,
  TypeVector<typename CompareTypes<T, Scalar>::supertype> >::type
TypeVector<T>::operator * (const Scalar factor) const
{
  typedef typename CompareTypes<T, Scalar>::supertype TS;

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



template <typename T, typename Scalar>
inline
typename boostcopy::enable_if_c<
  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 DIM == 1
  _coords[0] *= factor;
#endif
  
#if DIM == 2
  _coords[0] *= factor;
  _coords[1] *= factor;
#endif
  
#if DIM == 3
  _coords[0] *= factor;
  _coords[1] *= factor;
  _coords[2] *= factor;
#endif

  return *this;
}



template <typename T>
template <typename Scalar>
inline
typename boostcopy::enable_if_c<
  ScalarTraits<Scalar>::value,
  TypeVector<typename CompareTypes<T, Scalar>::supertype> >::type
TypeVector<T>::operator / (const Scalar factor) const
{
  libmesh_assert (factor != static_cast<T>(0.));

  typedef typename CompareTypes<T, Scalar>::supertype TS;
  
#if DIM == 1
  return TypeVector<TS>(_coords[0]/factor);
#endif
  
#if DIM == 2 
  return TypeVector<TS>(_coords[0]/factor,
		        _coords[1]/factor);
#endif
  
#if 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 (factor != static_cast<T>(0.));
  
  for (unsigned int i=0; i<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 DIM == 1
  return _coords[0]*p._coords[0];
#endif
  
#if DIM == 2
  return (_coords[0]*p._coords[0] +
	  _coords[1]*p._coords[1]);
#endif
  
#if DIM == 3
  return (_coords[0]*p(0) +
	  _coords[1]*p(1) +
	  _coords[2]*p(2));
#endif
}



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 (DIM == 3);

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

  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]);
}



template <typename T>
inline
Real TypeVector<T>::size() const
{
  return std::sqrt(this->size_sq());  
}



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



template <typename T>
inline
Real TypeVector<T>::size_sq() const
{
#if DIM == 1
  return (libmesh_norm(_coords[0]));
#endif
  
#if DIM == 2
  return (libmesh_norm(_coords[0]) +
	  libmesh_norm(_coords[1]));
#endif
  
#if DIM == 3
  return (libmesh_norm(_coords[0]) +
	  libmesh_norm(_coords[1]) + 
	  libmesh_norm(_coords[2]));
#endif
}



template <typename T>
inline
bool TypeVector<T>::absolute_fuzzy_equals(const TypeVector<T>& rhs, Real tol) const
{
#if DIM == 1
  return (std::abs(_coords[0] - rhs._coords[0])
	  <= tol);
#endif

#if DIM == 2
  return (std::abs(_coords[0] - rhs._coords[0]) +
	  std::abs(_coords[1] - rhs._coords[1])
	  <= tol);
#endif

#if DIM == 3
  return (std::abs(_coords[0] - rhs._coords[0]) +
	  std::abs(_coords[1] - rhs._coords[1]) +
	  std::abs(_coords[2] - rhs._coords[2])
	  <= tol);
#endif
}



template <typename T>
inline
bool TypeVector<T>::relative_fuzzy_equals(const TypeVector<T>& rhs, Real tol) const
{
#if DIM == 1
  return this->absolute_fuzzy_equals(rhs, tol *
				     (std::abs(_coords[0]) + std::abs(rhs._coords[0])));
#endif

#if DIM == 2
  return this->absolute_fuzzy_equals(rhs, tol * 
				     (std::abs(_coords[0]) + std::abs(rhs._coords[0]) +
				      std::abs(_coords[1]) + std::abs(rhs._coords[1])));
#endif

#if DIM == 3
  return this->absolute_fuzzy_equals(rhs, tol *
				     (std::abs(_coords[0]) + std::abs(rhs._coords[0]) +
				      std::abs(_coords[1]) + std::abs(rhs._coords[1]) +
				      std::abs(_coords[2]) + std::abs(rhs._coords[2])));
#endif
}



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

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

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



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

#endif // #define __type_vector_h__