//   
//   Project Name:        Kratos       
//   Last Modified by:    $Author: pooyan $
//   Date:                $Date: 2006-11-27 16:07:33 $
//   Revision:            $Revision: 1.1.1.1 $
//
//


#if !defined(KRATOS_FILENAME_H_INCLUDED )
#define  KRATOS_FILENAME_H_INCLUDED



// System includes


// External includes
#include <boost/array.hpp> 


// Project includes
#include "includes/define.h"
#include "geometries/geometry.h"


namespace Kratos
{

  ///@name Kratos Globals
  ///@{ 
  
  ///@} 
  ///@name Type Definitions
  ///@{ 
  
  ///@} 
  ///@name  Enum's
  ///@{
      
  ///@}
  ///@name  Functions 
  ///@{
      
  ///@}
  ///@name Kratos Classes
  ///@{
  
  /** 
  */ 
  template<class TPointType>
  class ClassName : public Geometry<TPointType>
    {
    public:
      ///@}
      ///@name Type Definitions
      ///@{
      
      /// Geometry as base class.
      typedef Geometry<TPointType> BaseType;

      /// Pointer definition of ClassName
      KRATOS_CLASS_POINTER_DEFINITION(ClassName);
  
      /** Integration methods implemented in geometry. 
      */
      typedef GeometryData::IntegrationMethod IntegrationMethod;

     /** A Vector of counted pointers to Geometries. Used for
	  returning edges of the geometry.
      */
      typedef BaseType::GeometryArrayType GeometryArrayType;
      
      /** Redefinition of template parameter TPointType.
       */
      typedef TPointType PointType;
      
      /** Type used for indexing in geometry class.std::size_t used for indexing
	  point or integration point access methods and also all other
	  methods which need point or integration point index.
      */
      typedef BaseType::IndexType IndexType;
      
      
      /** This typed used to return size or dimension in
	  geometry. Dimension, WorkingDimension, PointsNumber and
	  ... return this type as their results.
      */
      typedef BaseType::SizeType SizeType;
      
      /** Array of counted pointers to point. This type used to hold
	  geometry's points.
      */
      typedef  BaseType::PointsArrayType PointsArrayType;
      
      /** This type used for representing an integration point in
	  geometry. This integration point is a point with an
	  additional weight component.
      */
      typedef BaseType::IntegrationPointType IntegrationPointType;
      
      /** A Vector of IntegrationPointType which used to hold
	  integration points related to an integration
	  method. IntegrationPoints functions used this type to return
	  their results.
      */
      typedef BaseType::IntegrationPointsArrayType IntegrationPointsArrayType;
      
      /** A Vector of IntegrationPointsArrayType which used to hold
	  integration points related to different integration method
	  implemented in geometry.
      */
      typedef BaseType::IntegrationPointsContainerType IntegrationPointsContainerType;
      
      /** A third order tensor used as shape functions' values
	  continer.
      */
      typedef BaseType::ShapeFunctionsValuesContainerType ShapeFunctionsValuesContainerType;
      
      /** A fourth order tensor used as shape functions' local
	  gradients container in geometry.
      */
      typedef BaseType::ShapeFunctionsLocalGradientsContainerType ShapeFunctionsLocalGradientsContainerType;
      
      /** A third order tensor to hold jacobian matrices evaluated at
	  integration points. Jacobian and InverseOfJacobian functions
	  return this type as their result.
      */
      typedef BaseType::JacobiansType JacobiansType;
      
      /** A third order tensor to hold shape functions' local
	  gradients. ShapefunctionsLocalGradients function return this
	  type as its result.
      */
      typedef BaseType::ShapeFunctionsGradientsType ShapeFunctionsGradientsType;
      
      /** Type of the normal vector used for normal to edges in geomety.
       */
      typedef BaseType::NormalType NormalType;
      
      ///@}
      ///@name Life Cycle 
      ///@{ 
      
      ClassName(const PointsArrayType& ThisPoints) 
	: BaseType(ThisPoints, msGeometryData)
	{
	}

      /** Copy constructor.
	  Construct this geometry as a copy of given geometry.

	  @note This copy constructor don't copy the points and new
	  geometry shares points with given source geometry. It's
	  obvious that any change to this new geometry's point affect
	  source geometry's points too.
      */
      ClassName(ClassName const& rOther)
	: BaseType(rOther)
	{
	}


      /** Copy constructor from a geometry with other point type.
	  Construct this geometry as a copy of given geometry which
	  has different type of points. The given goemetry's
	  TOtherPointType* must be implicity convertible to this
	  geometry PointType.

	  @note This copy constructor don't copy the points and new
	  geometry shares points with given source geometry. It's
	  obvious that any change to this new geometry's point affect
	  source geometry's points too.
      */
      template<class TOtherPointType> ClassName(ClassName<TOtherPointType> const& rOther)
	: BaseType(rOther)
	{
	}

      /// Destructor. Do nothing!!!
      virtual ~ClassName(){}
      

      ///@}
      ///@name Operators 
      ///@{
      
      /** Assignment operator.

      @note This operator don't copy the points and this
      geometry shares points with given source geometry. It's
      obvious that any change to this geometry's point affect
      source geometry's points too.
      
      @see Clone
      @see ClonePoints
      */
      ClassName& operator=(const ClassName& rOther)
	{
	  BaseType::operator=(rOther);

	  return *this;
	}

      /** Assignment operator for geometries with different point type.

      @note This operator don't copy the points and this
      geometry shares points with given source geometry. It's
      obvious that any change to this geometry's point affect
      source geometry's points too.
      
      @see Clone
      @see ClonePoints
      */
      template<class TOtherPointType> 
	ClassName& operator=(ClassName<TOtherPointType> const & rOther)
	{
	  BaseType::operator=(rOther);

	  return *this;
	}
      
      ///@}
      ///@name Operations
      ///@{
      
      virtual boost::shared_ptr< Geometry< Point<3> > > Clone() const
	{
	  Geometry< Point<3> >::PointsArrayType NewPoints;

	  //making a copy of the nodes TO POINTS (not Nodes!!!)
	  for(IndexType i = 0 ; i < mPoints.size() ; i++)
	    NewPoints.push_back(mPoints[i]);

	  //creating a geometry with the new points
	  boost::shared_ptr< Geometry< Point<3> > > p_clone(new ClassName< Point<3> >(NewPoints));
	  p_clone->ClonePoints();

	  return p_clone;
	}
      
      //lumping factors for the calculation of the lumped mass matrix
      virtual void LumpingFactors(Vector& Result) 
	{ 
	  KRATOS_THROW_ERROR(std::logic_error, "Called the virtual function for LumpingFactors" , *this); 
	}

      ///@}
      ///@name Informations
      ///@{

      /** This method calculate and return Length or charactereistic
	  length of this geometry depending to it's dimension. For one
	  dimensional geometry for example Line it returns length of it
	  and for the other geometries it gives Characteristic length
	  otherwise.

	  @return double value contains length or Characteristic
	  length
	  @see Area()
	  @see Volume()
	  @see DomainSize()
      */
      virtual double Length()
	{
	}

      /** This method calculate and return area or surface area of
	  this geometry depending to it's dimension. For one dimensional
	  geometry it returns zero, for two dimensional it gives area
	  and for three dimensional geometries it gives surface area.

	  @return double value contains area or surface
	  area.
	  @see Length()
	  @see Volume()
	  @see DomainSize()
      */
      virtual double Area()
	{
	}

      /** This method calculate and return volume of this
	  geometry. For one and two dimensional geometry it returns
	  zero and for three dimensional it gives volume of geometry.

	  @return double value contains volume.
	  @see Length()
	  @see Area()
	  @see DomainSize()
      */
      virtual double Volume()
	{
	}


      /** This method calculate and return length, area or volume of
	  this geometry depending to it's dimension. For one dimensional
	  geometry it returns its length, for two dimensional it gives area
	  and for three dimensional geometries it gives its volume.

	  @return double value contains length, area or volume.
	  @see Length()
	  @see Area()
	  @see Volume()
      */
      virtual double DomainSize()
	{
	}

      /** Calculates center of this geometry by a simple averaging algorithm.
	  Each center point component calculated using:
	  \f[
	  c_i = \sum_j^n(x_i^j) / n
	  \f]

	  where \f$ c_i \f$ is component i of center point and \f$
	  X_i^j \f$ is component i of j'th point of geometry and n is
	  number of the points in this geometry.

	  @return PointType which is the calculated center of this geometry.
      */
      virtual PointType Center() const
	{
	}

      ///@}
      ///@name Access
      ///@{ 
      

      ///@}
      ///@name Inquiry
      ///@{
      
      /** This method is to know if this geometry is symmetric or
	  not.

	  @todo Making some method related to symmetry axis and more...

	  @return bool true if this geometry is symmetric and false if
	  it's not.
      */
      virtual bool IsSymmetric() const
	{
	  return false;
	}


      ///@}      
      ///@name Jacobian 
      ///@{


      /** Jacobians for given  method. This method
	  calculate jacobians matrices in all integrations points of
	  given integration method.

	  @param ThisMethod integration method which jacobians has to
	  be calculated in its integration points.

	  @return JacobiansType a Vector of jacobian
	  matrices \f$ J_i \f$ where \f$ i=1,2,...,n \f$ is the integration
	  point index of given integration method.

	  @see DeterminantOfJacobian
	  @see InverseOfJacobian
      */
      virtual void Jacobian(JacobiansType& Result, IntegrationMethod ThisMethod)
	{
	  KRATOS_THROW_ERROR(std::logic_error,
		       "Calling base class Jacobian method instead of drived class one. Please check the definition of derived class." , *this);
	}

      /** Jacobian in specific integration point of given integration
	  method. This method calculate jacobian matrix in given
	  integration point of given integration method.

	  @param IntegrationPointIndex index of integration point which jacobians has to
	  be calculated in it.

	  @param ThisMethod integration method which jacobians has to
	  be calculated in its integration points.

	  @return Matrix<double> Jacobian matrix \f$ J_i \f$ where \f$
	  i \f$ is the given integration point index of given
	  integration method.

	  @see DeterminantOfJacobian
	  @see InverseOfJacobian
      */
      virtual void Jacobian(Matrix& Result, IndexType IntegrationPointIndex, IntegrationMethod ThisMethod)
	{
	  KRATOS_THROW_ERROR(std::logic_error,
		       "Calling base class Jacobian method instead of drived class one. Please check the definition of derived class." , *this);
	}

      /** Jacobian in given point. This method calculate jacobian
	  matrix in given point.

	  @param rPoint point which jacobians has to
	  be calculated in it.

	  @return Matrix of double which is jacobian matrix \f$ J \f$ in given point.

	  @see DeterminantOfJacobian
	  @see InverseOfJacobian
      */
      virtual void Jacobian(Matrix& Result, const PointType& rPoint)
	{
	  KRATOS_THROW_ERROR(std::logic_error,
		       "Calling base class Jacobian method instead of drived class one. Please check the definition of derived class." , *this);
	}

      /** Determinant of jacobians for given integration method. This
	  method calculate determinant of jacobian in all
	  integrations points of given integration method.

	  @return Vector of double which is vector of determinants of
	  jacobians \f$ |J|_i \f$ where \f$ i=1,2,...,n \f$ is the
	  integration point index of given integration method.

	  @see Jacobian
	  @see InverseOfJacobian
      */
      virtual Vector DeterminantOfJacobian(Vector& Result, IntegrationMethod ThisMethod)
	{
	  KRATOS_THROW_ERROR(std::logic_error,
		       "Calling base class DeterminantOfJacobian method instead of drived class one. Please check the definition of derived class." , *this);
	}

      /** Determinant of jacobian in specific integration point of
	  given integration method. This method calculate determinant
	  of jacobian in given integration point of given integration
	  method.

	  @param IntegrationPointIndex index of integration point which jacobians has to
	  be calculated in it.

	  @param IntegrationPointIndex index of integration point
	  which determinant of jacobians has to be calculated in it.

	  @return Determinamt of jacobian matrix \f$ |J|_i \f$ where \f$
	  i \f$ is the given integration point index of given
	  integration method.

	  @see Jacobian
	  @see InverseOfJacobian
      */
      virtual double DeterminantOfJacobian(IndexType IntegrationPointIndex, IntegrationMethod ThisMethod)
	{
	  KRATOS_THROW_ERROR(std::logic_error,
		       "Calling base class DeterminantOfJacobian method instead of drived class one. Please check the definition of derived class." , *this);

	  return 0;
	}


      /** Determinant of jacobian in given point. This method calculate determinant of jacobian
	  matrix in given point.

	  @param rPoint point which determinant of jacobians has to
	  be calculated in it.

	  @return Determinamt of jacobian matrix \f$ |J| \f$ in given
	  point.

	  @see DeterminantOfJacobian
	  @see InverseOfJacobian
      */
      virtual double DeterminantOfJacobian(const Point<3>& rPoint)
	{
	  KRATOS_THROW_ERROR(std::logic_error,
		       "Calling base class DeterminantOfJacobian method instead of drived class one. Please check the definition of derived class." , *this);

	  return 0;
	}


      /** Inverse of jacobians for given integration method. This method
	  calculate inverse of jacobians matrices in all integrations points of
	  given integration method.

	  @param ThisMethod integration method which inverse of jacobians has to
	  be calculated in its integration points.

	  @return Inverse of jacobian
	  matrices \f$ J^{-1}_i \f$ where \f$ i=1,2,...,n \f$ is the integration
	  point index of given integration method.

	  @see Jacobian
	  @see DeterminantOfJacobian
      */
      virtual void InverseOfJacobian(JacobiansType& Result, IntegrationMethod ThisMethod)
	{
	  KRATOS_THROW_ERROR(std::logic_error,
		       "Calling base class InverseOfJacobian method instead of drived class one. Please check the definition of derived class." , *this);
	}

      /** Inverse of jacobian in specific integration point of given integration
	  method. This method calculate Inverse of jacobian matrix in given
	  integration point of given integration method.

	  @param IntegrationPointIndex index of integration point which inverse of jacobians has to
	  be calculated in it.

	  @param ThisMethod integration method which inverse of jacobians has to
	  be calculated in its integration points.

	  @return Inverse of jacobian matrix \f$ J^{-1}_i \f$ where \f$
	  i \f$ is the given integration point index of given
	  integration method.

	  @see Jacobian
	  @see DeterminantOfJacobian
      */
      virtual void InverseOfJacobian(Matrix Result, IndexType IntegrationPointIndex, IntegrationMethod ThisMethod)
	{
	  KRATOS_THROW_ERROR(std::logic_error,
		       "Calling base class InverseOfJacobian method instead of drived class one. Please check the definition of derived class." , *this);
	}

      /** Inverse of jacobian in given point. This method calculate inverse of jacobian
	  matrix in given point.

	  @param rPoint point which inverse of jacobians has to
	  be calculated in it.

	  @return Inverse of jacobian matrix \f$ J^{-1} \f$ in given point.

	  @see DeterminantOfJacobian
	  @see InverseOfJacobian
      */
      virtual void InverseOfJacobian(Matrix Result, const PointType& rPoint)
	{
	  KRATOS_THROW_ERROR(std::logic_error,
		       "Calling base class InverseOfJacobian method instead of drived class one. Please check the definition of derived class." , *this);
	}



      ///@}
      ///@name Shape Function
      ///@{ 

      /** This method gives value of given shape function evaluated in given
	  point.

	  @param rPoint Point of evaluation of the shape
	  function. This point must be in local coordinate.

	  @param ShapeFunctionIndex index of node which correspounding
	  shape function evaluated in given integration point.

	  @return Value of given shape function in given point.

	  @see ShapeFunctionsValues
	  @see ShapeFunctionsLocalGradients
	  @see ShapeFunctionLocalGradient
      */ 
      virtual double ShapeFunctionValue(IndexType ShapeFunctionIndex, const PointType& rPoint)
	{
	  KRATOS_THROW_ERROR(std::logic_error,
		       "Calling base class DeterminantOfJacobian method instead of drived class one. Please check the definition of derived class." , *this);

	  return 0;
	}



      virtual void CalculateShapeFunctionsIntegrationPointsGradients(ShapeFunctionsGradientsType& rResult)
	{
	  ShapeFunctionsIntegrationPointsGradients(mGeometryData.DefaultIntegrationMethod(), rResult);
	}

      virtual void ShapeFunctionsIntegrationPointsGradients(IntegrationMethod ThisMethod, ShapeFunctionsGradientsType& rResult)
	{
	  KRATOS_THROW_ERROR(std::logic_error,
		       "Calling base class ShapeFunctionsGaussPointsGradients method instead of derived class one. Please check the definition of derived class." , *this);
	}

      
      ///@}      
      ///@name Input and output
      ///@{

      /** Turn back information as a string.

      @return String contains information about this geometry.
      @see PrintData()
      @see PrintInfo()
      */
      virtual std::string Info() const
	{
	}

      /** Print information about this object.

      @param rOStream Stream to print into it.
      @see PrintData()
      @see Info()
      */
      virtual void PrintInfo(std::ostream& rOStream) const
	{
	}

      /** Print geometry's data into given stream. Prints it's points
	  by the order they stored in the geometry and then center
	  point of geometry.

	  @param rOStream Stream to print into it.
	  @see PrintInfo()
	  @see Info()
      */
      virtual void PrintData(std::ostream& rOStream) const
	{
	}      
      
            
      ///@}      
      ///@name Friends
      ///@{
      
            
      ///@}
      
    protected:
      ///@name Protected static Member Variables 
      ///@{ 
        
        
      ///@} 
      ///@name Protected member Variables 
      ///@{ 
        
        
      ///@} 
      ///@name Protected Operators
      ///@{ 
        
        
      ///@} 
      ///@name Protected Operations
      ///@{ 
        
        
      ///@} 
      ///@name Protected  Access 
      ///@{ 
        
        
      ///@}      
      ///@name Protected Inquiry 
      ///@{ 
        
        
      ///@}    
      ///@name Protected LifeCycle 
      ///@{ 
      

      ///@}
      
    private:
      ///@name Static Member Variables 
      ///@{ 
        
      static const GeometryData msGeometryData;
        
      ///@} 
      ///@name Member Variables 
      ///@{ 
        

      ///@} 
      ///@name Private Operators
      ///@{ 
        
        
      ///@} 
      ///@name Private Operations
      ///@{

      static ShapeFunctionsValuesContainerType CalculateShapeFunctionsIntegrationPointsValues()
      {
      }
        
      static ShapeFunctionsLocalGradientsContainerType CalculateShapeFunctionsIntegrationPointsLocalGradients()
      {
      }

      ///@} 
      ///@name Private  Access 
      ///@{ 
        
        
      ///@}    
      ///@name Private Inquiry 
      ///@{ 
        
        
      ///@}    
      ///@name Private Friends
      ///@{
      
      template<class TOtherPointType> friend class ClassName;

      ///@}    
      ///@name Un accessible methods 
      ///@{ 
      
      ClassName();

            
        
      ///@}    
        
    }; // Class Geometry 

  ///@} 
  
  ///@name Type Definitions       
  ///@{ 
  
  
  ///@} 
  ///@name Input and output 
  ///@{ 
        
 
      /// input stream function
      template<class TPointType>
	inline std::istream& operator >> (std::istream& rIStream, 
					  ClassName<TPointType>& rThis);

      /// output stream function
      template<class TPointType>
	inline std::ostream& operator << (std::ostream& rOStream, 
					  const ClassName<TPointType>& rThis)
	{
	  rThis.PrintInfo(rOStream);
	  rOStream << std::endl;
	  rThis.PrintData(rOStream);

	  return rOStream;
	}
      ///@} 

  template<class TPointType>  
  typename Triangle2D<TPointType>::IntegrationPointsContainerType Triangle2D<TPointType>::msIntegrationPoints = {
	  Quadrature<TriangleGaussLegendreIntegrationPoints<1>, 2, IntegrationPoint<3> >::IntegrationPoints(),
	  Quadrature<TriangleGaussLegendreIntegrationPoints<2>, 2, IntegrationPoint<3> >::IntegrationPoints()
  };
	

  template<class TPointType> 
  typename Triangle2D<TPointType>::ShapeFunctionsValuesContainerType 
  Triangle2D<TPointType>::msShapeFunctionsValues = {
	  Triangle2D<TPointType>::CalculateShapeFunctionsIntegrationPointsValues(GeometryType::GI_GAUSS_1),
	  Triangle2D<TPointType>::CalculateShapeFunctionsIntegrationPointsValues(GeometryType::GI_GAUSS_2)
  };

	
  template<class TPointType> 
  typename Triangle2D<TPointType>::ShapeFunctionsLocalGradientsContainerType 
  Triangle2D<TPointType>::msShapeFunctionsLocalGradients = {
	  Triangle2D<TPointType>::CalculateShapeFunctionsIntegrationPointsLocalGradients(GeometryType::GI_GAUSS_1),  
	  Triangle2D<TPointType>::CalculateShapeFunctionsIntegrationPointsLocalGradients(GeometryType::GI_GAUSS_2)
  };
  template<class TPointType>  
  typename ClassName<TPointType>::IntegrationPointsContainerType ClassName<TPointType>::msIntegrationPoints = {}	
	

  template<class TPointType> 
  typename ClassName<TPointType>::ShapeFunctionsValuesContainerType 
  ClassName<TPointType>::msShapeFunctionsValues = {
//	  ClassName<TPointType>::CalculateShapeFunctionsIntegrationPointsValues(GeometryType::GI_GAUSS_1),
//	  ClassName<TPointType>::CalculateShapeFunctionsIntegrationPointsValues(GeometryType::GI_GAUSS_2)
  };
	
  template<class TPointType> 
  typename ClassName<TPointType>::ShapeFunctionsLocalGradientsContainerType 
  ClassName<TPointType>::msShapeFunctionsLocalGradients( = {
//	  ClassName<TPointType>::CalculateShapeFunctionsIntegrationPointsLocalGradients(GeometryType::GI_GAUSS_1),  
//	  ClassName<TPointType>::CalculateShapeFunctionsIntegrationPointsLocalGradients(GeometryType::GI_GAUSS_2)
  };
      
  template<class TPointType>  
  const GeometryData ClassName<TPointType>::msGeometryData(/*Dimension*/,
							  /*WorkingSpaceDimension*/,
							  /*LocalSpaceDimension*/, 
							  msIntegrationPoints,
							  msShapeFunctionsValues,
							  msShapeFunctionsLocalGradients);

  
}  // namespace Kratos.

#endif // KRATOS_FILENAME_H_INCLUDED  defined 


