//  This file is part of ff3d - http://www.freefem.org/ff3d
//  Copyright (C) 2003 Pascal Hav
//
//  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: FEMFunctionBase.hpp,v 1.10 2007/05/20 23:02:47 delpinux Exp $

#ifndef FEM_FUNCTION_BASE_HPP
#define FEM_FUNCTION_BASE_HPP

#include <ScalarFunctionBase.hpp>

#include <DegreeOfFreedomPositionsSet.hpp>
#include <DegreeOfFreedomSetManager.hpp>

#include <DiscretizationType.hpp>

#include <Vector.hpp>

class Mesh;

/**
 * @file   FEMFunctionBase.hpp
 * @author Stephane Del Pino
 * @date   Wed Jul 19 17:34:19 2006
 * 
 * @brief  base class for finite element functions
 */
class FEMFunctionBase
  : public ScalarFunctionBase
{
protected:
  ConstReferenceCounting<Mesh>
  __baseMesh;			/**< reference to the mesh through the
				   Mesh (basis) class */

  const DegreeOfFreedomPositionsSet&
  __dofPositionsSet;		/**< degree of freedom position set */

  Vector<real_t> __values;	/**< vector of values */

  const DiscretizationType::Type
  __discretizationType;		/**< type of discretization*/

  /** 
   * output of the FEMFunction
   * 
   * @param os output stream
   * 
   * @return os
   */
  std::ostream& __put(std::ostream& os) const
  {
    os << "{fem-" << DiscretizationType::name(__discretizationType) << '}';
    return os;
  }

  /**
   * this value is returned when function is evaluated outside the
   * mesh.
   * 
   * @note It is required to allow expression like @f$ f+a @f$ (where
   * @f$ f @f$ is a fem-function and @f$ a @f$ a constant) to be
   * simplifed (FEM functions are considered to be 0 outside the mesh
   * by convention).
   */
  real_t __outsideValue;

public:
  /** 
   * set outside value
   * 
   * @param value outside value
   */
  void setOutsideValue(const real_t& value)
  {
    __outsideValue = value;
  }

  /** 
   * Access to the outside value
   * 
   * @return __outsidValue
   */
  const real_t& outsideValue() const
  {
    return __outsideValue;
  }

  /** 
   * Checks if the function can be simplified
   * 
   * @return true
   */
  bool canBeSimplified() const
  {
    return true;
  }

  /** 
   * Check if a finite element function is of the same type
   * 
   * @param f given function
   * 
   * @return true if discretization type and mesh are the same
   */
  bool hasSameType(const FEMFunctionBase& f) const
  {
    return ((__discretizationType == f.__discretizationType)
	    and (__baseMesh == f.__baseMesh));
  }

  /** 
   * Access to the mesh through its base class
   * 
   * @return __baseMesh
   */
  ConstReferenceCounting<Mesh> baseMesh() const
  {
    return __baseMesh;
  }

  /** 
   * The type of the finite elememt function
   * 
   * @return __discretizationType
   */
  const DiscretizationType::Type& discretizationType() const
  {
    return __discretizationType;
  }

  /** 
   * Read-only access to the set of values
   * 
   * @return __values
   */
  const Vector<real_t>& values() const
  {
    return __values;
  }

  /** 
   * access to the value at the ith degree of freedom
   * 
   * @param i number of the degree of freedom
   * 
   * @return __values[i]
   */
  inline real_t& operator[](const size_t& i)
  {
    return __values[i];
  }

  /** 
   * Read only Access to the value at the ith degree of freedom
   * 
   * @param i number of the degree of freedom
   * 
   * @return __values[i]
   */
  inline const real_t& operator[](const size_t& i) const
  {
    return __values[i];
  }

  /** 
   * Affects the vector values
   * 
   * @param values new values
   */
  void operator=(const Vector<real_t>& values)
  {
    ASSERT (values.size() == __values.size());
    __values = values;
  }

  /** 
   * Affects a function to a FEMFunction
   * 
   * @param f original function
   */
  virtual void operator=(const ScalarFunctionBase& f) = 0;

  /** 
   * Returns the gradient of a FEMFunction at position @a x
   * 
   * @param x position of evaluation
   * 
   * @return @f$ \nabla f (x) @f$
   */
  virtual TinyVector<3,real_t>
  gradient(const TinyVector<3,real_t>& x) const = 0;

  /** 
   * Constructor
   * 
   * @param mesh the mesh supporting the function
   * @param discretizationType the type of discretization
   */
  FEMFunctionBase(const Mesh* mesh,
		  DiscretizationType::Type discretizationType)
    : ScalarFunctionBase(ScalarFunctionBase::femfunction),
      __baseMesh(mesh),
      __dofPositionsSet(DegreeOfFreedomSetManager::instance().getDOFPositionsSet(*mesh,discretizationType)),
      __values(__dofPositionsSet.number()),
      __discretizationType(discretizationType),
      __outsideValue(0)
  {
    ;
  }

  /** 
   * Destructor
   * 
   */
  virtual ~FEMFunctionBase()
  {
    DegreeOfFreedomSetManager::instance().unsubscribe(__dofPositionsSet);
  }
};

#endif // FEM_FUNCTION_BASE_HPP