#ifndef FILE_GRIDFUNCTION
#define FILE_GRIDFUNCTION

/*********************************************************************/
/* File:   gridfunction.hpp                                           */
/* Author: Joachim Schoeberl                                         */
/* Date:   10. Jul. 2000                                             */
/*********************************************************************/


/** 
    Grid-functions
*/

class GridFunction : public NGS_Object
{
protected:
  const FESpace & fespace;
  bool nested;
  bool visual;
  int multidim;
  int level_updated;
  NgMutex mutex;
public:
  ///
  GridFunction (const FESpace & afespace, const string & name, const Flags & flags);
  ///
  virtual ~GridFunction () { ; }
  
  ///
  virtual void Update () = 0;
  ///  
  virtual BaseVector & GetVector (int comp = 0) = 0;

  ///
  void SetNested (int anested = 1)
  { nested = anested; }

  ///
  void SetVisual (bool avisual = 1)
  { visual = avisual; }

  int GetMultiDim () const { return multidim; }
  
  int GetLevelUpdated() const { return level_updated; }
  ///
  const FESpace & GetFESpace() const
  { return fespace; }

  
  NgMutex & Mutex () { return mutex; }

  ///
  virtual string GetClassName () const
  {
    return "GridFunction";
  }

  ///
  virtual void PrintReport (ostream & ost);
  ///
  virtual void MemoryUsage (ARRAY<MemoryUsageStruct*> & mu) const;

  ///
  ///
  void Visualize(const string & name);
};



template <class SCAL>
class S_GridFunction : public GridFunction
{
public:
  S_GridFunction (const FESpace & afespace, const string & aname, const Flags & flags)
    : GridFunction (afespace, aname, flags) { ; }
  
  ///
  virtual void GetElementVector (const ARRAY<int> & dnums,
				 FlatVector<SCAL> & elvec) const = 0;
  ///
  virtual void SetElementVector (const ARRAY<int> & dnums,
				 const FlatVector<SCAL> & elvec) = 0;


  ///
  virtual void GetElementVector (int comp,
				 const ARRAY<int> & dnums,
				 FlatVector<SCAL> & elvec) const = 0;
  ///
  virtual void SetElementVector (int comp,
				 const ARRAY<int> & dnums,
				 const FlatVector<SCAL> & elvec) = 0;
};



template <class TV>
class T_GridFunction : public S_GridFunction<typename mat_traits<TV>::TSCAL>
{
private:
  ARRAY<VVector<TV>*> vec;
  
public:
  typedef typename mat_traits<TV>::TSCAL TSCAL;
  enum { VDIM = mat_traits<TV>::HEIGHT };

  T_GridFunction (const FESpace & afespace, const string & aname, const Flags & flags);
  virtual ~T_GridFunction ();

  virtual void Update ();

  virtual BaseVector & GetVector (int comp = 0);

  ///
  virtual void GetElementVector (const ARRAY<int> & dnums,
				 FlatVector<TSCAL> & elvec) const;

  ///
  virtual void SetElementVector (const ARRAY<int> & dnums,
				 const FlatVector<TSCAL> & elvec);


  ///
  virtual void GetElementVector (int comp,
				 const ARRAY<int> & dnums,
				 FlatVector<TSCAL> & elvec) const;

  ///
  virtual void SetElementVector (int comp,
				 const ARRAY<int> & dnums,
				 const FlatVector<TSCAL> & elvec);


  /*
  // only implemented double element-const gf
  template <int S, int R>
  double Evaluate (const SpecificIntegrationPoint<S,R> & ip, int comp = 1)
  {
    
  const FESpace & fes = GetFESpace();
    
  const ElementTransformation & eltrans = ip.GetElementTransformation();
    
  int elnr = eltrans.GetElementNr();
    
  const FiniteElement & el = fes.GetFE(elnr);
  int nd = el.GetNDof();
    
  ARRAY<int> dnums (nd); 
  fes.GetDofNrs (elnr, dnums);
    
  Vector<double> elemu (nd);
    
  GetVector().GetIndirect (dnums, elemu);
  fes.TransformVec (elnr, false, elemu, TRANSFORM_SOL);
    
  return  InnerProduct (el.GetShape(ip), elemu);

  return 0;
  }
  */
  friend class GridFunctionCoefficientFunction;
};


extern GridFunction * CreateGridFunction (const FESpace * space,
					  const string & name, const Flags & flags);






///
// NOTE: only for L2-gridfunctions of order 0 !!
class GridFunctionCoefficientFunction : public CoefficientFunction
{
protected:
  ///
  T_GridFunction<double> & gf;
public:
  ///
  GridFunctionCoefficientFunction (T_GridFunction<double> & agf);
  ///
  virtual ~GridFunctionCoefficientFunction ();
  ///

  template <int S, int R>
  double Evaluate (const SpecificIntegrationPoint<S,R> & ip)
  {
    int elnr = ip.GetTransformation().GetElementNr();
    
    if (elnr < 0 || elnr >= gf.vec[0]->Size())
      {
	ostringstream ost;
	ost << "GridFunctionCoefficientFunction: Element nr. "
	    << elnr << " out of range 0 - " << gf.vec[0]->Size()-1 << endl;
	throw Exception (ost.str());
      }

    return (gf.vec[0]->FV())(elnr);
  }
  
  virtual double Evaluate (const BaseSpecificIntegrationPoint & ip)
  {
    int elnr = ip.GetTransformation().GetElementNr();

    if (elnr < 0 || elnr >= gf.vec[0]->Size())
      {
	ostringstream ost;
	ost << "GridFunctionCoefficientFunction: Element nr. "
	    << elnr << " out of range 0 - " << gf.vec[0]->Size()-1 << endl;
	throw Exception (ost.str());
      }

    return (gf.vec[0]->FV())(elnr);
  }
  
};







template <class SCAL>
class VisualizeGridFunction : public netgen::SolutionData
{
  const MeshAccess & ma;
  const S_GridFunction<SCAL> * gf;
  const BilinearFormIntegrator * bfi2d;
  const BilinearFormIntegrator * bfi3d;
  bool applyd;
  //
  int cache_elnr;
  bool cache_bound;
  LocalHeap lh;
  ElementTransformation eltrans;
  ARRAY<int> dnums;
  FlatVector<SCAL> elu;

public:
  VisualizeGridFunction (const MeshAccess & ama,
			 const GridFunction * agf,
			 const BilinearFormIntegrator * abfi2d,
			 const BilinearFormIntegrator * abfi3d,
			 bool aapplyd);
  /*

  : SolutionData (agf->GetName(), -1, agf->GetFESpace().IsComplex()),
  ma(ama), gf(dynamic_cast<const S_GridFunction<SCAL>*> (agf)), 
  bfi2d(abfi2d), bfi3d(abfi3d), applyd(aapplyd), cache_elnr(-1), lh(1000000)
  { 
  if (bfi2d) components = bfi2d->DimFlux();
  if (bfi3d) components = bfi3d->DimFlux();
  if (iscomplex) components *= 2;
  multidimcomponent = 0;
  }
  */

  virtual ~VisualizeGridFunction ();
  
  virtual bool GetValue (int elnr, 
			 double lam1, double lam2, double lam3,
			 double * values) ;

  /*
    { 
    if (!bfi3d) return 0;

    LocalHeapMem<10000> lh2;

    const FESpace & fes = gf->GetFESpace();
    const FiniteElement & fel = fes.GetFE (elnr, lh2);
    int dim     = fes.GetDimension();

    // added 07.04.2004 (FB): don't draw, if not defined on this domain
    if ( !fes.DefinedOn(ma.GetElIndex(elnr)) ) return 0;

    if (cache_elnr != elnr || cache_bound)
    {
    lh.CleanUp();
    ma.GetElementTransformation (elnr, eltrans, lh);
    fes.GetDofNrs (elnr, dnums);
    
    elu.AssignMemory (dnums.Size() * dim, lh);
    gf->GetElementVector (multidimcomponent, dnums, elu);
    fes.TransformVec (elnr, 0, elu, TRANSFORM_SOL);

    cache_elnr = elnr;
    cache_bound = 0;
    }

    void * hp = lh.GetPointer();

    FlatVector<SCAL> flux;
    
    IntegrationPoint ip(lam1, lam2, lam3, 0);
    bfi3d->CalcFlux (fel, eltrans, ip, elu, flux, applyd, lh);

    (*testout) << "visualize element " << elnr 
    << "dnums = " << endl << dnums << endl
    << "elu = " << endl << elu << endl;

    (*testout) << "flux = " << flux << endl;
    
    for (int i = 0; i < components; i++)
    values[i] = ((double*)(void*)&flux(0))[i];
    
    lh.CleanUp(hp);
    return 1; 
    }
  */

  virtual bool GetSurfValue (int elnr,
			     double lam1, double lam2, 
			     double * values) ;

  /*
    { 
    if (!bfi2d) return 0;

    bool bound = (ma.GetDimension() == 3);
    const FESpace & fes = gf->GetFESpace();

    LocalHeapMem<10000> lh2;

    const FiniteElement & fel = 
    bound ? fes.GetSFE(elnr, lh2) : fes.GetFE (elnr, lh2);

    int dim     = fes.GetDimension();

    // added 07.04.2004 (FB): don't draw, if not defined on this domain
    if ( bound ? 
    !fes.DefinedOnBoundary(ma.GetSElIndex(elnr)) :
    !fes.DefinedOn(ma.GetElIndex(elnr)) ) return 0;


    if (cache_elnr != elnr || !cache_bound)
    {
    lh.CleanUp();

    if (bound)
    {
    ma.GetSurfaceElementTransformation (elnr, eltrans, lh);
    fes.GetSDofNrs (elnr, dnums);
    }
    else
    {
    ma.GetElementTransformation (elnr, eltrans, lh);
    fes.GetDofNrs (elnr, dnums);
    }

    elu.AssignMemory (dnums.Size() * dim, lh);
    gf->GetElementVector (multidimcomponent, dnums, elu);
    fes.TransformVec (elnr, bound, elu, TRANSFORM_SOL);

    cache_elnr = elnr;
    cache_bound = 1;
    }

    void * hp = lh.GetPointer();

    FlatVector<SCAL> flux;

    IntegrationPoint ip(lam1, lam2, 0, 0);
    //    bfi2d->CalcFlux (fel, eltrans, ip, elu, flux, applyd, lh);
    if (bound)
    {
    SpecificIntegrationPoint<2,3> sip (ip, eltrans, lh);
    bfi2d->CalcFlux (fel, sip, elu, flux, applyd, lh);
    }
    else
    {
    SpecificIntegrationPoint<2,2> sip (ip, eltrans, lh);
    bfi2d->CalcFlux (fel, sip, elu, flux, applyd, lh);
    }

    for (int i = 0; i < components; i++)
    values[i] = ((double*)(void*)&flux(0))[i];

    lh.CleanUp(hp);
    return 1; 
    }
  */

  void Analyze(ARRAY<double> & minima, ARRAY<double> & maxima, ARRAY<double> & averages, int component = -1);
  void Analyze(ARRAY<double> & minima, ARRAY<double> & maxima, ARRAY<double> & averages_times_volumes, ARRAY<double> & volumes, int component = -1);
    

};






#endif