#ifndef FILE_INTRULE
#define FILE_INTRULE

/*********************************************************************/
/* File:   intrule.hh                                                */
/* Author: Joachim Schoeberl                                         */
/* Date:   25. Mar. 2000                                             */
/*********************************************************************/



// class IntegrationPoint;


/// An integration point 
class IntegrationPoint
{
private:
  /// number within intergration Rule
  int nr; 
  /// enumerates ALL intergration points for element type
  int glob_nr; 
  /// coordinates (empty values for 1D and 2D)
  double pi[3];
  /// weight of integration point
  double weight;

public:
  ///
  IntegrationPoint ()
    : nr(-1), glob_nr(-1) { ; }
  ///
  IntegrationPoint (double api[3], double aw);
  ///
  IntegrationPoint (double p1, double p2 = 0, double p3 = 0, double aw = 0);
  ///
  IntegrationPoint (const FlatVector<double> & ap, double aw);

  ///
  void SetNr (int anr) { nr = anr; }
  ///
  void SetGlobNr (int aipnr) { glob_nr = aipnr; }
  ///
  const double * Point () const { return pi; }
  ///
  int Size() const { return 3; }
  ///
  const double & operator() (int i) const { return pi[i]; }
  ///
  double & operator() (int i) { return pi[i]; }

  /*
  ///
  double X() const { return pi[0]; }
  ///
  double Y() const { return pi[1]; }
  ///
  double Z() const { return pi[2]; }
  */

  ///
  double Weight () const { return weight; }
  /// number within integration rule
  int Nr () const { return nr; }

  /// global number of ip
  int IPNr () const { return glob_nr; }

  /*
  ///
  bool operator== (const IntegrationPoint & ip2) const
  {
    return (pi[0] == ip2.pi[0] &&
	    pi[1] == ip2.pi[1] &&
	    pi[2] == ip2.pi[2] &&
	    weight == ip2.weight);
  }
  */
  ///
  friend ostream & operator<< (ostream & ost, const IntegrationPoint & ip);
};





// class BaseSpecificIntegrationPoint;
// template <int R, typename SCAL> class DimSpecificIntegrationPoint;
// template <int S, int R, typename SCAL> class SpecificIntegrationPoint;
// class IntegrationRule;

class ElementTransformation;


/**
   Bsae class for SpecificIntegrationPoint.
   A specific integration point is the mapped point, and stores
   point coordinates, the Jacobimatrix, and the determinant of the Jacobimatrix.
*/
class BaseSpecificIntegrationPoint
{
protected:
  /// IP on the reference element
  const IntegrationPoint & ip;
  /// 
  const ElementTransformation & eltrans;
  ///
  // int ipnr;

public:
  BaseSpecificIntegrationPoint (const IntegrationPoint & aip,
				const ElementTransformation & aeltrans)
    : ip(aip), eltrans(aeltrans) /* , ipnr(-1) */ { ; }

  ///
  const IntegrationPoint & IP () const { return ip; }
  ///
  const ElementTransformation & GetTransformation () const { return eltrans; }
  ///
  int GetIPNr() const { return ip.Nr(); }
};

template <int R, typename SCAL = double>
class DimSpecificIntegrationPoint : public BaseSpecificIntegrationPoint
{
protected:
  ///
  Vec<R,SCAL> point;
public:
  ///
  DimSpecificIntegrationPoint (const IntegrationPoint & aip,
			       const ElementTransformation & aeltrans)
    : BaseSpecificIntegrationPoint (aip, aeltrans)
  { ; }
  ///
  const Vec<R,SCAL> GetPoint () const { return point; }
  ///
  SCAL operator() (int i) const 
  { return point(i); }
};


/// ip, dimension source, dimension range
template <int DIMS = 2, int DIMR = 2, typename SCAL = double> 
class SpecificIntegrationPoint : public DimSpecificIntegrationPoint<DIMR,SCAL>
{
private:
  ///
  Mat<DIMR,DIMS,SCAL> dxdxi;
  ///
  Mat<DIMS,DIMR,SCAL> dxidx;
  /// 
  SCAL det;
  ///
  Vec<DIMR,SCAL> normalvec;

public:
  ///
  SpecificIntegrationPoint (const IntegrationPoint & aip,
			    const ElementTransformation & aeltrans,
			    LocalHeap & lh);

  ///
  const Mat<DIMR,DIMS,SCAL> & GetJacobian() const { return dxdxi; }
  ///
  SCAL GetJacobiDet() const { return det; }
  ///
  const Mat<DIMS,DIMR,SCAL> & GetJacobianInverse() const { return dxidx; }
  ///
  const Vec<DIMR,SCAL> GetNV () const { return normalvec; }
  ///
  int IsBoundary () const { return DIMS != DIMR; }
};



/**
   An integration rule.
   Contains array of integration points
*/
class IntegrationRule
{
  ///
  ARRAY<IntegrationPoint*> ipts;
public:
  ///
  IntegrationRule ();
  ///
  IntegrationRule (int nips);
  ///
  ~IntegrationRule ();
  ///
  void AddIntegrationPoint (IntegrationPoint * ip);

  /// number of integration points
  int GetNIP() const { return ipts.Size(); }

  /// returns i-th integration point
  const IntegrationPoint & operator[] (int i) const
  {
#ifdef CHECK_RANGE
    if (i < 0 || i >= ipts.Size())
      {
	stringstream st;
	st << "Integration Point out of range " << i << endl;
	throw Exception (st.str());
      }
#endif
    return *ipts[i]; 
  }


  // old style, should be replaced by []
  const IntegrationPoint & GetIP(int i) const
  {
#ifdef CHECK_RANGE
    if (i < 0 || i >= ipts.Size())
      {
	stringstream st;
	st << "Integration Point out of range " << i << endl;
	throw Exception (st.str());
      }
#endif
    return *ipts[i]; 
  }
};


/** 
   Integration Rules.
   A global class maintaining integration rules. If a rule of specific
   order is requested for the first time, than the rule is generated.
*/
class IntegrationRules
{
  ///
  ARRAY<IntegrationRule*> segmentrules;
  ///
  IntegrationRule segmentlumping;
  ///
  ARRAY<IntegrationPoint*> segmentpoints;

  ///
  ARRAY<IntegrationRule*> trigrules;
  ///
  ARRAY<IntegrationRule*> trignodalrules;
  ///
  IntegrationRule triglumping;
  ///
  IntegrationRule triglumping2;
  ///
  ARRAY<IntegrationPoint*> trigpoints;

  ///
  ARRAY<IntegrationRule*> quadrules;
  ///
  IntegrationRule quadlumping;
  ///
  ARRAY<IntegrationPoint*> quadpoints;

  ///
  ARRAY<IntegrationRule*> tetrules;
  ///
  IntegrationRule tetlumping;
  ///
  ARRAY<IntegrationRule*> tetnodalrules;
  ///
  ARRAY<IntegrationPoint*> tetpoints;

  ///
  ARRAY<IntegrationRule*> prismrules;
  ///
  IntegrationRule prismlumping;
  ///
  IntegrationRule prismfacemidpoint;
  ///
  ARRAY<IntegrationPoint*> prismpoints;

  ///
  ARRAY<IntegrationRule*> pyramidrules;
  ///
  IntegrationRule pyramidlumping;
  ///
  ARRAY<IntegrationPoint*> pyramidpoints;

  ///
  ARRAY<IntegrationRule*> hexrules;
  ///
  ARRAY<IntegrationPoint*> hexpoints;

public:
  ///
  IntegrationRules ();
  ///
  ~IntegrationRules ();
  /// returns the integration rule for the element type and integration order
  const IntegrationRule & SelectIntegrationRule (ELEMENT_TYPE eltyp, int order) const;
  ///
  const IntegrationRule & SelectLumpingIntegrationRule (ELEMENT_TYPE eltyp) const;
  ///
  const IntegrationRule & SelectNodalIntegrationRule (ELEMENT_TYPE eltyp, int order) const;
  ///
  const IntegrationRule & PrismFaceMidPoint () const
  { return prismfacemidpoint; }
  ///
  const ARRAY<IntegrationPoint*> & GetIntegrationPoints (ELEMENT_TYPE eltyp) const;
  ///
  const IntegrationRule & GenerateIntegrationRule (ELEMENT_TYPE eltyp, int order);
};


/**
   Computes Gaussean integration.
   Integration rule on (0,1), contains n points,
   Exact for polynomials up to order 2n-1
*/ 
extern void ComputeGaussRule (int n, 
			      ARRAY<double> & xi, 
			      ARRAY<double> & wi);

/// Get a reference to the integration-rules container
extern const IntegrationRules & GetIntegrationRules ();


#endif