//  This file is part of ff3d - http://www.freefem.org/ff3d
//  Copyright (C) 2001-2005 Stphane Del Pino

//  This program is free software; you can redistribute it and/or modify
//  it under the terms of the GNU General Public License as published by
//  the Free Software Foundation; either version 2, or (at your option)
//  any later version.

//  This program is distributed in the hope that it will be useful,
//  but WITHOUT ANY WARRANTY; without even the implied warranty of
//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
//  GNU General Public License for more details.

//  You should have received a copy of the GNU General Public License
//  along with this program; if not, write to the Free Software Foundation,
//  Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.  

//  $Id: FEMFunction.hpp,v 1.7 2007/05/20 23:02:47 delpinux Exp $

#ifndef FEM_FUNCTION_HPP
#define FEM_FUNCTION_HPP

#include <FiniteElementTraits.hpp>
#include <FEMFunctionBase.hpp>

/**
 * @file   FEMFunction.hpp
 * @author Stephane Del Pino
 * @date   Wed Jul 19 16:15:18 2006
 * 
 * @brief This class manages fem functions for various finite element
 * types
 * 
 */
template <typename MeshType,
	  typename FiniteElementTraits>
class FEMFunction
  : public FEMFunctionBase
{
private:
  typedef typename MeshType::CellType 
  CellType;			/**< @typedef type of cells of the
				   mesh */
  typedef

  typename FiniteElementTraits::Transformation
  Transformation;		/**< @typedef type of the conform
				   transformation */

  typedef
  typename FiniteElementTraits::Type
  FiniteElementType;		/**< @typedef type of the finite
				   element */

  ConstReferenceCounting<MeshType>
  __mesh;			/**< the mesh where the function
				   leaves */

  /** 
   * Copy constructor
   * 
   * @param f given function
   * @note forbidden to avoid wild copies
   */
  FEMFunction(const FEMFunction<MeshType,FiniteElementTraits>& f);
public:
  /** 
   * Evaluates the FEMFunction at point \a x.
   * 
   * @param x the position of evaluation
   * 
   * @return \f$ f(x) \f$
   */
  real_t operator()(const TinyVector<3,real_t>& x) const
  {
    typename MeshType::const_iterator icell = __mesh->find(x);
    if (icell.end()) {
      return __outsideValue;
    }

    const CellType& K = *icell;
    Transformation T(K);

    TinyVector<3, real_t> xhat;
    T.invertT(x, xhat);

    real_t value = 0;
    for (size_t l=0; l<FiniteElementType::numberOfDegreesOfFreedom; ++l) {
      value += __values[__dofPositionsSet(icell.number(),l)]*FiniteElementType::instance().W(l,xhat);
    }

    return value;
  }

  /** 
   * Affects a function to a FEMFunction
   * 
   * @param f original function
   */
  void operator=(const ScalarFunctionBase& f)
  {
    // during an affectation, outside values are set to 0 (the fem
    // function is "pure")
    __outsideValue = 0;

    if (f.type() == this->type()) {
      const FEMFunctionBase& femBase = dynamic_cast<const FEMFunctionBase&>(f);
      if ((femBase.discretizationType() == this->discretizationType()) and
	  (femBase.baseMesh() == this->baseMesh())) {
	// if the function is of the same kind: just copy values
	const FEMFunction<MeshType,FiniteElementTraits>& fem
	  = dynamic_cast<const FEMFunction<MeshType,FiniteElementTraits>&>(femBase);
	__values = fem.__values;
	return;
      }
    }
    for (size_t i=0; i<__values.size(); i++) {
      const TinyVector<3>& x = __dofPositionsSet.vertex(i);
      __values[i] = f(x);
    }
  }

  /** 
   * Evaluates the function's gradient at position @a x
   * 
   * @param x position of evaluation
   * 
   * @return @f$ \nabla f (x) \nabla @f$
   */
  TinyVector<3,real_t>
  gradient(const TinyVector<3,real_t>& x) const
  {
    typename MeshType::const_iterator icell = __mesh->find(x);
    if (icell.end()) {
      return 0;
    }

    const CellType& K = *icell;
    Transformation T(K);

    TinyVector<3, real_t> xhat;
    T.invertT(x, xhat);

    TinyVector<3, real_t> referenceGradient = 0;
    for (size_t l=0; l<FiniteElementType::numberOfDegreesOfFreedom; ++l) {
      const real_t value = __values[__dofPositionsSet(icell.number(),l)];
      referenceGradient[0] += value*FiniteElementType::instance().dxW(l,xhat);
      referenceGradient[1] += value*FiniteElementType::instance().dyW(l,xhat);
      referenceGradient[2] += value*FiniteElementType::instance().dzW(l,xhat);
    }

    TinyMatrix<3,3, real_t> J;
    {
      TinyVector<3, real_t> temp;

      T.dx(x,temp);
      for(size_t i=0; i<3; ++i) {
	J(0,i) = temp[i];
      }

      T.dy(x,temp);
      for(size_t i=0; i<3; ++i) {
	J(1,i) = temp[i];
      }

      T.dz(x,temp);
      for(size_t i=0; i<3; ++i) {
	J(2,i) = temp[i];
      }
    }
    // now we use 
    TinyVector<3, real_t> result;

    gaussPivot(J, referenceGradient, result);
    return result;
  }

  /** 
   * Evaluates first derivative of the function
   * 
   * @param x position of evaluation
   * 
   * @return @f$ \partial_x f at position x @f$
   */
  real_t dx(const TinyVector<3>& x) const
  {
    return gradient(x)[0];
  }

  /** 
   * Evaluates second derivative of the function
   * 
   * @param x position of evaluation
   * 
   * @return @f$ \partial_y f at position x @f$
   */
  real_t dy(const TinyVector<3>& x) const
  {
    return gradient(x)[1];
  }

  /** 
   * Evaluates third derivative of the function
   * 
   * @param x position of evaluation
   * 
   * @return @f$ \partial_z f at position x @f$
   */
  real_t dz(const TinyVector<3>& x) const
  {
    return gradient(x)[2];
  }

  /** 
   * Constructor
   * 
   * @param mesh mesh supporting the function
   */
  FEMFunction(ConstReferenceCounting<MeshType> mesh)
    : FEMFunctionBase(mesh, DiscretizationType::Type(FiniteElementTraits::discretizationType)),
      __mesh(mesh)
  {
    ;
  }

  /** 
   * Constructor
   * 
   * @param mesh mesh supporting the function
   * @param f function of initialization
   */
  FEMFunction(ConstReferenceCounting<MeshType> mesh,
	      const ScalarFunctionBase& f)
    : FEMFunctionBase(mesh, DiscretizationType::Type(FiniteElementTraits::discretizationType)),
      __mesh(mesh)
  {
    (*this) = f;
  }

  /** 
   * Constructor
   * 
   * @param mesh mesh supporting the function
   * @param d value of initialization
   */
  FEMFunction(ConstReferenceCounting<MeshType> mesh,
	      const real_t& d)
    : FEMFunctionBase(mesh, DiscretizationType::Type(FiniteElementTraits::discretizationType)),
      __mesh(mesh)
  {
    __values = d;
  }

  /** 
   * Constructor
   * 
   * @param mesh given mesh
   * @param values given values
   */
  FEMFunction(ConstReferenceCounting<MeshType> mesh,
	      const Vector<real_t>& values)
    : FEMFunctionBase(mesh, DiscretizationType::Type(FiniteElementTraits::discretizationType)),
      __mesh(mesh)
  {
    ASSERT(__values.size() == values.size());
    __values = values;
  }

  /** 
   * Destructor
   * 
   */
  ~FEMFunction()
  {
    ;
  }
};

#endif // FEM_FUNCTION_HPP