#ifndef FILE_COEFFICIENT
#define FILE_COEFFICIENT

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

/** 
   coefficient functions
*/

class CoefficientFunction
{
public:
  ///
  CoefficientFunction ();
  ///
  virtual ~CoefficientFunction ();
  ///

  virtual double Evaluate (const BaseSpecificIntegrationPoint & ip) = 0;

  virtual double Evaluate (const BaseSpecificIntegrationPoint & ip, const double & t)
  {
    return Evaluate(ip);
  }
  virtual double EvaluateDeri (const BaseSpecificIntegrationPoint & ip, const double & t)
  {
    return 0;
  }

  // to be changed
  virtual double Evaluate (const BaseSpecificIntegrationPoint & ip,
			   const complex<double> & t)
  { return Evaluate(ip,t.real()); }
  // to be changed
    virtual double EvaluateDeri (const BaseSpecificIntegrationPoint & ip,
				 const complex<double> & t)
  { return EvaluateDeri(ip,t.real()); }


  virtual double EvaluateConst () 
  {
    throw Exception (string ("EvaluateConst called for non-const coefficient function ")+
		     typeid(*this).name());
  }

  virtual void PrintReport (ostream & ost)
  {
    ost << "Base-Class CoefficientFunction" << endl;
  }
};



template <int S, int R>
double Evaluate (CoefficientFunction & fun,
		 const SpecificIntegrationPoint<S,R> & ip) 
{ 
  return fun.Evaluate(ip); 
}

/*
template <>
inline double Evaluate (CoefficientFunction & fun,
			const SpecificIntegrationPoint<1,1> & ip)
{
  return fun.Evaluate11 (ip);
}

template <>
inline double Evaluate (CoefficientFunction & fun,
			const SpecificIntegrationPoint<1,2> & ip)
{
  return fun.Evaluate12 (ip);
}

template <>
inline double Evaluate (CoefficientFunction & fun,
			const SpecificIntegrationPoint<2,2> & ip)
{
  return fun.Evaluate22 (ip);
}

template <>
inline double Evaluate (CoefficientFunction & fun,
			const SpecificIntegrationPoint<2,3> & ip)
{
  return fun.Evaluate23 (ip);
}

template <>
inline double Evaluate (CoefficientFunction & fun,
			const SpecificIntegrationPoint<3,3> & ip)
{
  return fun.Evaluate33 (ip);
}





/// Base-class for template-polymorphismus
template <class T>
class SpecCoefficientFunction : public CoefficientFunction
{
public:
  SpecCoefficientFunction () { ; }
  virtual ~SpecCoefficientFunction () { ; }

  virtual double Evaluate11 (const SpecificIntegrationPoint<1,1> & ip) 
  {
    return static_cast<T*>(this) -> Evaluate (ip);
  }
  virtual double Evaluate12 (const SpecificIntegrationPoint<1,2> & ip)
  {
    return static_cast<T*>(this) -> Evaluate (ip);
  }
  virtual double Evaluate22 (const SpecificIntegrationPoint<2,2> & ip) 
  {
    return static_cast<T*>(this) -> Evaluate (ip);
  }
  virtual double Evaluate23 (const SpecificIntegrationPoint<2,3> & ip)
  {
    return static_cast<T*>(this) -> Evaluate (ip);
  }
  virtual double Evaluate33 (const SpecificIntegrationPoint<3,3> & ip) 
  {
    return static_cast<T*>(this) -> Evaluate (ip);
  }
};


*/


/// The coefficient is constant everywhere
class ConstantCoefficientFunction : public CoefficientFunction
// : public SpecCoefficientFunction<ConstantCoefficientFunction>
{
  ///
  double val;
public:
  ///
  ConstantCoefficientFunction (double aval);
  ///
  virtual ~ConstantCoefficientFunction ();
  ///

  template <int S, int R>
  double Evaluate (const SpecificIntegrationPoint<S,R> & ip)
  {
    return val;
  }

  virtual double Evaluate (const BaseSpecificIntegrationPoint & ip)
  {
    return val;
  }

  virtual double EvaluateConst () 
  {
    return val;
  }
};



/// The coefficient is constant in every sub-domain
class DomainConstantCoefficientFunction : public CoefficientFunction
//  : public SpecCoefficientFunction<DomainConstantCoefficientFunction>
{
  ///
  ARRAY<double> val;
public:
  ///
  DomainConstantCoefficientFunction (const ARRAY<double> & aval);
  ///
  virtual ~DomainConstantCoefficientFunction ();
  ///


  template <int S, int R>
  double Evaluate (const SpecificIntegrationPoint<S,R> & ip)
  {
    int elind = ip.GetTransformation().GetElementIndex();
    
    if (elind < 0 || elind >= val.Size())
      {
	ostringstream ost;
	ost << "DomainConstantCoefficientFunction: Element index "
	    << elind << " out of range 0 - " << val.Size()-1 << endl;
	throw Exception (ost.str());
      }
    
    return val[elind]; 
  }

  virtual double Evaluate (const BaseSpecificIntegrationPoint & ip)
  {
    int elind = ip.GetTransformation().GetElementIndex();
    
    if (elind < 0 || elind >= val.Size())
      {
	ostringstream ost;
	ost << "DomainConstantCoefficientFunction: Element index "
	    << elind << " out of range 0 - " << val.Size()-1 << endl;
	throw Exception (ost.str());
      }
    
    return val[elind]; 
  }

  virtual double EvaluateConst () 
  {
    return val[0];
  }
};




///
template <int DIM>
class DomainVariableCoefficientFunction : public CoefficientFunction
//  : public SpecCoefficientFunction<DomainVariableCoefficientFunction>
{
  ///
  ARRAY<EvalFunction*> fun;
public:
  ///
  DomainVariableCoefficientFunction (const ARRAY<EvalFunction*> & afun);
  ///
  virtual ~DomainVariableCoefficientFunction ();
  ///
  template <int S, int R>
  double Evaluate (const SpecificIntegrationPoint<S,R> & ip)
  {
    int elind = ip.GetTransformation().GetElementIndex();
    if (elind < 0 || elind >= fun.Size())
      {
	cerr << "element index " << elind << " out of range " 
	     << "0 - " << fun.Size()-1 << endl;
      return 0;
      }
    return fun[elind]->Eval ( &ip.GetPoint()(0));
  }

  virtual double Evaluate (const BaseSpecificIntegrationPoint & ip)
  {
    int elind = ip.GetTransformation().GetElementIndex();
    if (elind < 0 || elind >= fun.Size())
      {
	cerr << "element index " << elind << " out of range " 
	     << "0 - " << fun.Size()-1 << endl;
      return 0;
      }
    return fun[elind]->Eval (&static_cast<const DimSpecificIntegrationPoint<DIM>&>(ip).GetPoint()(0));
  }
};

///
template <int DIM>
class DomainInternalCoefficientFunction : public CoefficientFunction
{
  ///
  int matnr;
  ///
  double (*f)(const double*); 
public:
  ///
  DomainInternalCoefficientFunction (int amatnr, double (*af)(const double*))
    : matnr(amatnr), f(af) { ; }
  ///
  virtual ~DomainInternalCoefficientFunction () { ; }
  ///
  template <int S, int R>
  double Evaluate (const SpecificIntegrationPoint<S,R> & ip)
  {
    int elind = ip.GetTransformation().GetElementIndex();
    if (elind != matnr && matnr > 0) return 0;
  
    return f(&ip.GetPoint()(0));
  }
  
  virtual double Evaluate (const BaseSpecificIntegrationPoint & ip)
  {
    int elind = ip.GetTransformation().GetElementIndex();
    if (elind != matnr && matnr > 0) return 0;
  

    return f(&static_cast<const DimSpecificIntegrationPoint<DIM>&>(ip).GetPoint()(0));
    // return f(&ip.GetPoint().REval(0));
  }
  
};





/**
   coefficient function that is defined in every integration point.
   NOTE: for the constructor, the maximal number of integration
   points per element is required!
**/
class IntegrationPointCoefficientFunction : public CoefficientFunction
{
  int elems, ips_per_elem;
  ///
  ARRAY<double> values;
public:
  ///
  IntegrationPointCoefficientFunction (int aelems, int size)
    : elems(aelems), ips_per_elem(size), values(aelems*size) { ; }
  ///
  IntegrationPointCoefficientFunction (int aelems, int size, double val)
    : elems(aelems), ips_per_elem(size), values(aelems*size)
  {
    values = val;
  } 
  ///
  IntegrationPointCoefficientFunction (int aelems, int size, ARRAY<double> & avalues)
    : elems(aelems), ips_per_elem(size), values(avalues) 
  { 
    if ( avalues.Size() < aelems * size )
      {
	cout << "Warning: IntegrationPointCoefficientFunction, constructor: sizes don't match!" << endl;
	values.SetSize(aelems*size);
      }
  }

  ///
  virtual ~IntegrationPointCoefficientFunction () { ; }
  ///
  template <int S, int R>
  double Evaluate (const SpecificIntegrationPoint<S,R> & ip)
  {
    int ipnr = ip.GetIPNr();
    int elnr = ip.GetTransformation().GetElementNr();

    if ( ipnr < 0 || ipnr >= ips_per_elem || elnr < 0 || elnr >= elems ) 
      {
	ostringstream ost;
	ost << "IntegrationPointCoefficientFunction: ip = "
	    << ipnr << " / elem = " << elnr << ". Ranges: 0 - " 
	    << ips_per_elem << "/ 0 - " << elems << "!" << endl;
	throw Exception (ost.str());
      }
  
    return values[elnr*ips_per_elem+ipnr];
  }
  
  virtual double Evaluate (const BaseSpecificIntegrationPoint & ip)
  {
    int ipnr = ip.GetIPNr();
    int elnr = ip.GetTransformation().GetElementNr();

    if ( ipnr < 0 || ipnr >= ips_per_elem || elnr < 0 || elnr >= elems ) 
      {
	ostringstream ost;
	ost << "IntegrationPointCoefficientFunction: ip = "
	    << ipnr << " / elem = " << elnr << ". Ranges: 0 - " 
	    << ips_per_elem << "/ 0 - " << elems << "!" << endl;
	throw Exception (ost.str());
      }
  
    return values[elnr*ips_per_elem+ipnr];
  }


  // direct access to the values at the integration points
  double & operator() (int elnr, int ipnr)
  {
    if ( ipnr < 0 || ipnr >= ips_per_elem || elnr < 0 || elnr >= elems ) 
      {
	ostringstream ost;
	ost << "IntegrationPointCoefficientFunction: ip = "
	    << ipnr << " / elem = " << elnr << ". Ranges: 0 - " 
	    << ips_per_elem-1 << " / 0 - " << elems-1 << "!" << endl;
	throw Exception (ost.str());
      }
  
    return values[elnr*ips_per_elem+ipnr];
  }

  double operator() (int elnr, int ipnr) const
  {
    if ( ipnr < 0 || ipnr >= ips_per_elem || elnr < 0 || elnr >= elems ) 
      {
	ostringstream ost;
	ost << "IntegrationPointCoefficientFunction: ip = "
	    << ipnr << " / elem = " << elnr << ". Ranges: 0 - " 
	    << ips_per_elem-1 << " / 0 - " << elems-1 << "!" << endl;
	throw Exception (ost.str());
      }
  
    return values[elnr*ips_per_elem+ipnr];
  }



  int GetNumIPs() const { return ips_per_elem; }
  int GetNumElems() const { return elems; }



  void ReSetValues( ARRAY<double> & avalues )
  {
    if ( avalues.Size() < values.Size() )
      {
	throw Exception("IntegrationPointCoefficientFunction::ReSetValues - sizes don't match!");
      }
    values = avalues;
  }
  
};


/// Coefficient function that depends (piecewise polynomially) on a parameter
class PolynomialCoefficientFunction : public CoefficientFunction
{
private:
  ARRAY < ARRAY< ARRAY<double>* >* > polycoeffs;
  ARRAY < ARRAY<double>* > polybounds;

private:
  double EvalPoly(const double t, const ARRAY<double> & coeffs) const;
  double EvalPolyDeri(const double t, const ARRAY<double> & coeffs) const;

public:
  PolynomialCoefficientFunction(const ARRAY < ARRAY<double>* > & polycoeffs_in);
  PolynomialCoefficientFunction(const ARRAY < ARRAY< ARRAY<double>* >* > & polycoeffs_in, const ARRAY < ARRAY<double>* > & polybounds_in);
  
  virtual ~PolynomialCoefficientFunction();

  virtual double Evaluate (const BaseSpecificIntegrationPoint & ip);

  virtual double Evaluate (const BaseSpecificIntegrationPoint & ip, const double & t);

  virtual double EvaluateDeri (const BaseSpecificIntegrationPoint & ip, const double & t);

  virtual double EvaluateConst ();
};


#endif