type_tensor.h

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

// $Id: type_tensor.h 2789 2008-04-13 02:24:40Z roystgnr $

// The libMesh Finite Element Library.
// Copyright (C) 2002-2007  Benjamin S. Kirk, John W. Peterson
  
// 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 __type_tensor_h__
#define __type_tensor_h__

// C++ includes
#include <cmath>

// Local includes
#include "libmesh_common.h"
#include "type_vector.h"




/**
 * This class defines a tensor in \p DIM dimensional space of type T.
 * T may either be Real or Complex.  The default constructor for
 * this class is protected, suggesting that you should not instantiate
 * one of these directly.
 *
 * \author Roy Stogner, 2004.
 */

template <typename T>
class TypeTensor
{
template <typename T2>
friend class TypeTensor;

protected:

  /**
   * Constructor.  By default sets all entries to 0.  Gives the tensor 0 in
   * \p DIM dimensions.  This is a poor constructor for 2D tensors -
   * if the default arguments are to be overridden it requires that
   * the "z = 0." argument also be given explicitly.
   */
  TypeTensor  (const T xx=0.,
	       const T xy=0.,
	       const T xz=0.,
	       const T yx=0.,
	       const T yy=0.,
	       const T yz=0.,
	       const T zx=0.,
	       const T zy=0.,
	       const T zz=0.);
  
public:

  /**
   * Copy-constructor.
   */
  template<typename T2>
  TypeTensor(const TypeTensor<T2>& p);
  
  /**
   * Destructor.
   */ 
  ~TypeTensor();

  /**
   * Assign to a tensor without creating a temporary.
   */
  template<typename T2>
  void assign (const TypeTensor<T2> &);

  /**
   * Return the \f$ i,j^{th} \f$ element of the tensor.
   */
  T operator () (const unsigned int i, const unsigned int j) const;

  /**
   * Return a writeable reference to the \f$ i^{th} \f$ element of the
   * tensor.
   */
  T & operator () (const unsigned int i, const unsigned int j);

  /**
   * Add two tensors. 
   */
  template<typename T2>
  TypeTensor<typename CompareTypes<T, T2>::supertype>
  operator + (const TypeTensor<T2> &) const;

  /**
   * Add to this tensor. 
   */
  template<typename T2>
  const TypeTensor<T> & operator += (const TypeTensor<T2> &);
  
  /**
   * Add to this tensor without creating a temporary.
   */
  template<typename T2>
  void add (const TypeTensor<T2> &); 
  
  /**
   * Add a scaled tensor to this tensor without
   * creating a temporary.
   */
  template <typename T2>
  void add_scaled (const TypeTensor<T2> &, const T); 

  /**
   * Subtract two tensors.
   */
  template<typename T2>
  TypeTensor<typename CompareTypes<T, T2>::supertype>
  operator - (const TypeTensor<T2> &) const;

  /**
   * Subtract from this tensor. 
   */
  template<typename T2>
  const TypeTensor<T> & operator -= (const TypeTensor<T2> &);

  /**
   * Subtract from this tensor without creating a temporary.
   */
  template<typename T2>
  void subtract (const TypeTensor<T2> &); 
  
  /**
   * Subtract a scaled value from this tensor without
   * creating a temporary.
   */
  template <typename T2>
  void subtract_scaled (const TypeTensor<T2> &, const T); 
  
  /**
   * Return the opposite of a tensor 
   */
  TypeTensor<T> operator - () const;
  
  /**
   * Multiply a tensor by a number, i.e. scale.
   */
  template <typename Scalar>
  typename boostcopy::enable_if_c<
    ScalarTraits<Scalar>::value,
    TypeTensor<typename CompareTypes<T, Scalar>::supertype> >::type
  operator * (const Scalar) const;
  
  /**
   * Multiply this tensor by a number, i.e. scale.
   */
  template <typename Scalar>
  const TypeTensor<T> & operator *= (const Scalar);
  
  /**
   * Divide a tensor by a number, i.e. scale.
   */
  template <typename Scalar>
  typename boostcopy::enable_if_c<
    ScalarTraits<Scalar>::value,
    TypeTensor<typename CompareTypes<T, Scalar>::supertype> >::type
  operator / (const Scalar) const;

  /**
   * Divide this tensor by a number, i.e. scale.
   */
  const TypeTensor<T> & operator /= (const T);

  /**
   * Multiply 2 tensors together, i.e. matrix product.
   * The tensors may be of different types.
   */
  template <typename T2>
  TypeTensor<T> operator * (const TypeTensor<T2> &) const;

  /**
   * Multiply 2 tensors together, i.e. sum Aij*Bij.
   * The tensors may be of different types.
   */
  template <typename T2>
  typename CompareTypes<T,T2>::supertype
  contract (const TypeTensor<T2> &) const;

  /**
   * Multiply a tensor and vector together, i.e. matrix-vector product.
   * The tensor and vector may be of different types.
   */
  template <typename T2>
  TypeVector<typename CompareTypes<T,T2>::supertype>
  operator * (const TypeVector<T2> &) const;

  /**
   * The transpose (with complex numbers not conjugated) of the tensor.
   */
  TypeTensor<T> transpose() const;

  /**
   * Returns the Frobenius norm of the tensor, i.e. the square-root of
   * the sum of the elements squared.
   */
  Real size() const;

  /**
   * Returns the Frobenius norm of the tensor squared, i.e.  sum of the
   * element magnitudes squared.
   */
  Real size_sq() const;

  /**
   * Zero the tensor in any dimension.
   */
  void zero();

  /**
   * @returns \p true if two tensors are equal valued.
   */
  bool operator == (const TypeTensor<T>& rhs) const;
  
  /**
   * @returns \p true if this tensor is "less"
   * than another.  Useful for sorting.
   */
  bool operator < (const TypeTensor<T>& rhs) const;
  
  /**
   * @returns \p true if this tensor is "greater"
   * than another.
   */
  bool operator > (const TypeTensor<T>& rhs) const;
  
  /**
   * Formatted print to \p std::cout.
   */
  void print(std::ostream& os) const;

  /**
   * Formatted print as above but allows you to do
   * std::cout << t << std::endl;
   */
  friend std::ostream& operator << (std::ostream& os, const TypeTensor<T>& t)
  {
    t.print(os);
    return os;
  }
  
  /**
   * Unformatted print to the stream \p out.  Simply prints the elements
   * of the tensor separated by spaces and newlines.
   */ 
  void write_unformatted (std::ostream &out, const bool newline = true) const;
    
  // protected:

  /**
   * The coordinates of the \p TypeTensor
   */
  T _coords[DIM*DIM];
};



//------------------------------------------------------
// Inline functions
template <typename T>
inline
TypeTensor<T>::TypeTensor (const T xx,
	                   const T xy,
	                   const T xz,
	                   const T yx,
	                   const T yy,
	                   const T yz,
	                   const T zx,
	                   const T zy,
	                   const T zz)
{
  _coords[0] = xx;

  if (DIM == 2)
    {
      _coords[1] = xy;
      _coords[2] = yx;
      _coords[3] = yy;
    }

  if (DIM == 3)
    {
      _coords[1] = xy;
      _coords[2] = xz;
      _coords[3] = yx;
      _coords[4] = yy;
      _coords[5] = yz;
      _coords[6] = zx;
      _coords[7] = zy;
      _coords[8] = zz;
    }
}



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



template <typename T>
inline
TypeTensor<T>::~TypeTensor ()
{
}



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



template <typename T>
inline
T TypeTensor<T>::operator () (const unsigned int i,
			      const unsigned int j) const
{
  libmesh_assert (i<3);
  libmesh_assert (j<3);

#if DIM < 3
  if (i >= DIM || j >= DIM)
    return 0.;
#endif
  
  return _coords[i*DIM+j];
}



template <typename T>
inline
T & TypeTensor<T>::operator () (const unsigned int i,
				const unsigned int j)
{
#if DIM < 3

  if (i >= DIM || j >= DIM)
    {
//       std::cerr << "ERROR:  You are assigning to a tensor component" << std::endl
// 		<< "that is out of range for the compiled DIM!"      << std::endl
// 		<< " DIM=" << DIM << " , i=" << i << " , j=" << j << std::endl;
      libmesh_error();
    }
  
#endif
  
  libmesh_assert (i<DIM);
  libmesh_assert (j<DIM);
  
  return _coords[i*DIM+j];
}



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],