type_tensor.h

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

		    _coords[3] + p._coords[3],
		    _coords[4] + p._coords[4],
		    _coords[5] + p._coords[5],
		    _coords[6] + p._coords[6],
		    _coords[7] + p._coords[7],
		    _coords[8] + p._coords[8]);
#endif
	       
}



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

  return *this;
}



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



template <typename T>
template <typename T2>
inline
void TypeTensor<T>::add_scaled (const TypeTensor<T2> &p, const T factor)
{
  for (unsigned int i=0; i<DIM*DIM; i++)
    _coords[i] += factor*p._coords[i];

}



template <typename T>
template<typename T2>
inline
TypeTensor<typename CompareTypes<T, T2>::supertype>
TypeTensor<T>::operator - (const TypeTensor<T2> &p) const
{

#if DIM == 1
  return TypeTensor(_coords[0] - p._coords[0]);
#endif

#if DIM == 2 
  return TypeTensor(_coords[0] - p._coords[0],
		    _coords[1] - p._coords[1],
		    0.,
		    _coords[2] - p._coords[2],
		    _coords[3] - p._coords[3]);
#endif

#if DIM == 3
  return TypeTensor(_coords[0] - p._coords[0],
		    _coords[1] - p._coords[1],
		    _coords[2] - p._coords[2],
		    _coords[3] - p._coords[3],
		    _coords[4] - p._coords[4],
		    _coords[5] - p._coords[5],
		    _coords[6] - p._coords[6],
		    _coords[7] - p._coords[7],
		    _coords[8] - p._coords[8]);
#endif

}



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

  return *this;
}



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



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



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

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

#if DIM == 3
  return TypeTensor(-_coords[0],
		    -_coords[1], 
		    -_coords[2], 
		    -_coords[3], 
		    -_coords[4], 
		    -_coords[5], 
		    -_coords[6], 
		    -_coords[7], 
		    -_coords[8]);
#endif
  
}



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


#if DIM == 1
  return TypeTensor<TS>(_coords[0]*factor);
#endif
  
#if DIM == 2 
  return TypeTensor<TS>(_coords[0]*factor,
		        _coords[1]*factor,
		        _coords[2]*factor,
		        _coords[3]*factor);
#endif
  
#if DIM == 3
  return TypeTensor<TS>(_coords[0]*factor,
		        _coords[1]*factor, 
		        _coords[2]*factor, 
		        _coords[3]*factor, 
		        _coords[4]*factor, 
		        _coords[5]*factor, 
		        _coords[6]*factor, 
		        _coords[7]*factor, 
		        _coords[8]*factor);
#endif  
}



template <typename T, typename Scalar>
inline
typename boostcopy::enable_if_c<
  ScalarTraits<Scalar>::value,
  TypeTensor<typename CompareTypes<T, Scalar>::supertype> >::type
operator * (const Scalar factor,
	    const TypeTensor<T> &t)
{
  return t * factor;
}



template <typename T>
template <typename Scalar>
inline
const TypeTensor<T> & TypeTensor<T>::operator *= (const Scalar factor)
{
  for (unsigned int i=0; i<DIM*DIM; i++)
    _coords[i] *= factor;

  return *this;
}




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

  typedef typename CompareTypes<T, Scalar>::supertype TS;
  
#if DIM == 1
  return TypeTensor<TS>(_coords[0]/factor);
#endif
  
#if DIM == 2 
  return TypeTensor<TS>(_coords[0]/factor,
		        _coords[1]/factor,
		        _coords[2]/factor,
		        _coords[3]/factor);
#endif
  
#if DIM == 3
  return TypeTensor<TS>(_coords[0]/factor,
		        _coords[1]/factor, 
		        _coords[2]/factor, 
		        _coords[3]/factor, 
		        _coords[4]/factor, 
		        _coords[5]/factor, 
		        _coords[6]/factor, 
		        _coords[7]/factor, 
		        _coords[8]/factor);
#endif
  
}



template <typename T>
inline
TypeTensor<T> TypeTensor<T>::transpose() const
{
#if DIM == 1
  return TypeTensor(_coords[0]);
#endif
  
#if DIM == 2 
  return TypeTensor(_coords[0],
		    _coords[2],
		    _coords[1],
		    _coords[3]);
#endif
  
#if DIM == 3
  return TypeTensor(_coords[0],
		    _coords[3], 
		    _coords[6], 
		    _coords[1], 
		    _coords[4], 
		    _coords[7], 
		    _coords[2], 
		    _coords[5], 
		    _coords[8]);
#endif
}



template <typename T>
inline
const TypeTensor<T> & TypeTensor<T>::operator /= (const T factor)
{
  libmesh_assert (factor != static_cast<T>(0.));
  
  for (unsigned int i=0; i<DIM*DIM; i++)
    _coords[i] /= factor;

  return *this;
}




template <typename T>
template <typename T2>
inline
TypeVector<typename CompareTypes<T,T2>::supertype>
TypeTensor<T>::operator * (const TypeVector<T2> &p) const
{
  TypeVector<typename CompareTypes<T,T2>::supertype> returnval;
  for (unsigned int i=0; i<DIM; i++)
    for (unsigned int j=0; j<DIM; j++)
      returnval(i) += (*this)(i,j)*p(j);

  return returnval;
}



template <typename T>
template <typename T2>
inline
TypeTensor<T> TypeTensor<T>::operator * (const TypeTensor<T2> &p) const
{
  TypeTensor<T> returnval;
  for (unsigned int i=0; i<DIM; i++)
    for (unsigned int j=0; j<DIM; j++)
      for (unsigned int k=0; k<DIM; k++)
        returnval(i,j) += (*this)(i,k)*p(k,j);

  return returnval;
}



  /**
   * Multiply 2 tensors together, i.e. sum Aij*Bij.
   * The tensors may be of different types.
   */
template <typename T>
template <typename T2>
inline
typename CompareTypes<T,T2>::supertype
TypeTensor<T>::contract (const TypeTensor<T2> &t) const
{
  typename CompareTypes<T,T2>::supertype sum = 0.;
  for (unsigned int i=0; i<DIM*DIM; i++)
    sum += _coords[i]*t._coords[i];
  return sum;
}



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



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



template <typename T>
inline
Real TypeTensor<T>::size_sq () const
{
  Real sum = 0.;
  for (unsigned int i=0; i<DIM*DIM; i++)
    sum += libmesh_norm(_coords[i]);
  return sum;
}



template <typename T>
inline
bool TypeTensor<T>::operator == (const TypeTensor<T>& rhs) const
{
#if DIM == 1
  return (std::abs(_coords[0] - rhs._coords[0])
	  < TOLERANCE);
#endif
  
#if DIM == 2
  return ((std::abs(_coords[0] - rhs._coords[0]) +
	   std::abs(_coords[1] - rhs._coords[1]) +
	   std::abs(_coords[2] - rhs._coords[2]) +
	   std::abs(_coords[3] - rhs._coords[3]))
	  < 4.*TOLERANCE);
#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]) +
	   std::abs(_coords[3] - rhs._coords[3]) +
	   std::abs(_coords[4] - rhs._coords[4]) +
	   std::abs(_coords[5] - rhs._coords[5]) +
	   std::abs(_coords[6] - rhs._coords[6]) +
	   std::abs(_coords[7] - rhs._coords[7]) +
	   std::abs(_coords[8] - rhs._coords[8]))
	  < 9.*TOLERANCE);
#endif
  
}

#endif // #define __type_tensor_h__