/*********************************************************************/
/* File:   coefficient.cpp                                           */
/* Author: Joachim Schoeberl                                         */
/* Date:   24. Jan. 2002                                             */
/*********************************************************************/

/* 
   Finite Element Coefficient Function
*/


#include <fem.hpp>

namespace ngfem
{
  using namespace ngfem;

  
  CoefficientFunction :: CoefficientFunction ()
  { ; }

  CoefficientFunction :: ~CoefficientFunction ()
  { ; }

    ///
  ConstantCoefficientFunction ::   
  ConstantCoefficientFunction (double aval) 
    : val(aval) 
  { ; }

  ConstantCoefficientFunction ::
  ~ConstantCoefficientFunction ()
  { ; }


  DomainConstantCoefficientFunction :: 
  DomainConstantCoefficientFunction (const ARRAY<double> & aval)
    : val(aval) { ; }
  /*
  .Size())
  {
    for (int i = 0; i < val.Size(); i++)
      val[i] = aval[i];
  }
  */

  DomainConstantCoefficientFunction :: 
  ~DomainConstantCoefficientFunction ()
  { ; }




  template <int DIM>
  DomainVariableCoefficientFunction<DIM> ::
  DomainVariableCoefficientFunction (const ARRAY<EvalFunction*> & afun)
    : fun(afun.Size())
  {
    for (int i = 0; i < fun.Size(); i++)
      fun[i] = new EvalFunction (*afun[i]);
  }

  template <int DIM>
  DomainVariableCoefficientFunction<DIM> ::
  ~DomainVariableCoefficientFunction ()
  {
    for (int i = 0; i < fun.Size(); i++)
      delete fun[i];
  }


  template class DomainVariableCoefficientFunction<2>;
  template class DomainVariableCoefficientFunction<3>;


  PolynomialCoefficientFunction::PolynomialCoefficientFunction(const ARRAY < ARRAY< ARRAY<double>* >* > & polycoeffs_in,
							       const ARRAY < ARRAY<double>* > & polybounds_in)
    : polycoeffs(polycoeffs_in), polybounds(polybounds_in)
  {}

  PolynomialCoefficientFunction::PolynomialCoefficientFunction(const ARRAY < ARRAY<double>* > & polycoeffs_in)
  {
    polycoeffs.SetSize(polycoeffs_in.Size());
    polybounds.SetSize(polycoeffs_in.Size());
    
    for(int i=0; i<polycoeffs_in.Size(); i++)
      {
	polycoeffs[i] = new ARRAY< ARRAY<double>* >(1);
	(*polycoeffs[i])[0] = polycoeffs_in[i];
	polybounds[i] = new ARRAY<double>(0);
      } 
  }


  PolynomialCoefficientFunction::~PolynomialCoefficientFunction()
  {
    for(int i=0; i<polycoeffs.Size(); i++)
      {
	delete polybounds[i];
	for(int j=0; j<polycoeffs[i]->Size(); j++)
	  {
	    delete (*polycoeffs[i])[j];
	  }
	delete polycoeffs[i];
      }
    polycoeffs.DeleteAll();
    polybounds.DeleteAll();
  }
    
  
  
  double PolynomialCoefficientFunction::Evaluate (const BaseSpecificIntegrationPoint & ip)
  {
    return Evaluate(ip,0);
  }



  double PolynomialCoefficientFunction::EvalPoly(const double t, const ARRAY<double> & coeffs) const
  {
    const int last = coeffs.Size()-1;
    
    double retval = coeffs[last];
    for(int i=last-1; i>=0; i--)
      {
	retval *= t;
	retval += coeffs[i];
      }

    return retval;    
  }


  double PolynomialCoefficientFunction::EvalPolyDeri(const double t, const ARRAY<double> & coeffs) const
  {
    const int last = coeffs.Size()-1;

    double retval = last*coeffs[last];
    for(int i=last-1; i>=1; i--)
      {
	retval *= t;
	retval += i*coeffs[i];
      }

    return retval;    
  }


  double PolynomialCoefficientFunction::Evaluate (const BaseSpecificIntegrationPoint & ip, const double & t)
  {
    const int elind = ip.GetTransformation().GetElementIndex();
    
    if (elind < 0 || elind >= polycoeffs.Size())
      {
	ostringstream ost;
	ost << "PolynomialCoefficientFunction: Element index "
	    << elind << " out of range 0 - " << polycoeffs.Size()-1 << endl;
	throw Exception (ost.str());
      }
 
    int pos;
    for(pos=0; pos < polybounds[elind]->Size() && t > (*polybounds[elind])[pos]; pos++){}
   
    return EvalPoly(t,*((*(polycoeffs[elind]))[pos]));

    
  }


 
  double PolynomialCoefficientFunction::EvaluateDeri (const BaseSpecificIntegrationPoint & ip, const double & t)
  {
    const int elind = ip.GetTransformation().GetElementIndex();
    
    if (elind < 0 || elind >= polycoeffs.Size())
      {
	ostringstream ost;
	ost << "PolynomialCoefficientFunction: Element index "
	    << elind << " out of range 0 - " << polycoeffs.Size()-1 << endl;
	throw Exception (ost.str());
      }

    int pos;
    for(pos=0; pos < polybounds[elind]->Size() && t > (*polybounds[elind])[pos]; pos++){}

    return EvalPolyDeri(t,*((*(polycoeffs[elind]))[pos]));
  }


  double PolynomialCoefficientFunction::EvaluateConst () 
  {
    return (*(*polycoeffs[0])[0])[0];
  }

}