//  This file is part of ff3d - http://www.freefem.org/ff3d
//  Copyright (C) 2001, 2002, 2003 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: RealExpression.cpp,v 1.16 2006/02/05 20:58:21 delpinux Exp $

#include <RealExpression.hpp>
#include <BooleanExpression.hpp>

#include <FunctionExpression.hpp>
#include <Vector3Expression.hpp>

#include <MeshExpression.hpp>

#include <Variable.hpp>

#include <Structured3DMesh.hpp>
#include <MeshOfHexahedra.hpp>
#include <MeshOfTetrahedra.hpp>
#include <SurfaceMeshOfTriangles.hpp>
#include <SurfaceMeshOfQuadrangles.hpp>

#include <Information.hpp>

#include <Hexahedron.hpp>
#include <Tetrahedron.hpp>
#include <Triangle.hpp>
#include <Quadrangle.hpp>

#include <ConformTransformation.hpp>
#include <FiniteElementTraits.hpp>

#include <limits>

ReferenceCounting<RealExpression> RealExpressionVariable::value()
{
  return (*__expression).value();
}

void RealExpressionVariable::execute()
{
  __expression = (*__realVariable).expression();
  __realValue = (*__expression).realValue();
}

RealExpressionVariable::RealExpressionVariable(ReferenceCounting<RealVariable> realVariable)
  : __realVariable(realVariable),
    __expression(0)
{
  ;
}

RealExpressionVariable::RealExpressionVariable(const RealExpressionVariable& e)
  : __realVariable(e.__realVariable),
    __expression(e.__expression)
{
  ;
}

RealExpressionVariable::~RealExpressionVariable()
{
  ;
}

/* Function Evaluate */

ReferenceCounting<RealExpression>
RealExpressionFunctionEvaluate::value()
{
  return new RealExpressionValue(__realValue);
}

void RealExpressionFunctionEvaluate::execute()
{
  (*__realFunction).execute();
  if ((*__realFunction).hasBoundaryExpression()) {
    throw ErrorHandler(__FILE__,__LINE__,
		       "Cannot evaluate \""
		       +stringify(*__realFunction)+
		       "\": missing boundary",
		       ErrorHandler::normal);
  }
  TinyVector<3,real_t> X;
  if (__v == 0) {
    (*__x).execute();
    (*__y).execute();
    (*__z).execute();

    X[0] = (*__x).realValue();
    X[1] = (*__y).realValue();
    X[2] = (*__z).realValue();
  } else {
    (*__v).execute();
    for (size_t i=0; i<3; ++i) {
      X[i] = (*__v).value(i);
    }
  }
  __realValue = (*__realFunction).value(X[0],X[1],X[2]);
}

RealExpressionFunctionEvaluate::
RealExpressionFunctionEvaluate (ReferenceCounting<FunctionExpression> f,
				ReferenceCounting<RealExpression> x,
				ReferenceCounting<RealExpression> y,
				ReferenceCounting<RealExpression> z)
  : __realFunction(f),
    __v(0),
    __x(x),
    __y(y),
    __z(z)
{
  ;
}


RealExpressionFunctionEvaluate::
RealExpressionFunctionEvaluate(ReferenceCounting<FunctionExpression> f,
			       ReferenceCounting<Vector3Expression> v)
  : __realFunction(f),
    __v(v),
    __x(0),
    __y(0),
    __z(0)
{
  ;
}

RealExpressionFunctionEvaluate::
RealExpressionFunctionEvaluate(const RealExpressionFunctionEvaluate& e)
  : __realFunction(e.__realFunction),
    __v(e.__v),
    __x(e.__x),
    __y(e.__y),
    __z(e.__z)
{
  ;
}

RealExpressionFunctionEvaluate::~RealExpressionFunctionEvaluate()
{
  ;
}

/* Integrate */

ReferenceCounting<RealExpression> RealExpressionIntegrate::value()
{
  return new RealExpressionValue(__realValue);
}

template <typename MeshType, typename QuadratureType>
real_t RealExpressionIntegrate::__integrate(const MeshType& M,
					    const QuadratureType& Q,
					    FunctionExpression& f)
{
  real_t integral=0;
  typedef typename MeshType::CellType CellType;
  for (typename MeshType::const_iterator iCell(M);
       not(iCell.end()); ++iCell) {
    const CellType& C = *iCell;
    typename FiniteElementTraits<CellType,
                                 DiscretizationType::LagrangianFEM1>::Transformation T(C);
    typename FiniteElementTraits<CellType,
                                 DiscretizationType::LagrangianFEM1>::JacobianTransformation J(T);

    TinyVector<3, real_t> X;
    for (size_t i=0; i<QuadratureType::numberOfQuadraturePoints; ++i) {
      T.value(Q[i],X);
      integral += J.jacobianDet()*Q.weight(i)*f.value(X[0],X[1],X[2]);
    }
  }
  return integral;
}


template <typename MeshType>
real_t RealExpressionIntegrate::__integrate(const MeshType& M,
					    FunctionExpression& f)
{
  switch(this->__discretizationType) {
  case DiscretizationType::LagrangianFEM1: {
    typedef typename FiniteElementTraits<typename MeshType::CellType,
                                         DiscretizationType::LagrangianFEM1>::Type
      FiniteElementType;
    typedef typename FiniteElementType::QuadratureType QuadratureType;
    const QuadratureType& Q  = QuadratureType::instance();
    return __integrate(M,Q,f);
  }
  case DiscretizationType::LagrangianFEM2: {
    typedef typename FiniteElementTraits<typename MeshType::CellType,
                                         DiscretizationType::LagrangianFEM2>::Type
      FiniteElementType;
    typedef typename FiniteElementType::QuadratureType QuadratureType;
    const QuadratureType& Q  = QuadratureType::instance();
    return __integrate(M,Q,f);
  }
  default: {
    throw ErrorHandler(__FILE__,__LINE__,
		       "unexpected quadrature type",
		       ErrorHandler::unexpected);    
    return 0;
  }
  }
}

//! \todo compute the integral on the domain.
void RealExpressionIntegrate::execute()
{
  (*__mesh).execute();
  Information::instance().setMesh(__mesh);


  const Mesh& M = (*(*__mesh).mesh());

  (*__realFunction).execute();
  FunctionExpression& f = (*__realFunction);

  if (M.family() ==Mesh::volume) {
    if (f.hasBoundaryExpression()) {
      throw ErrorHandler(__FILE__,__LINE__,
			 "Cannot evaluate \""
			 +stringify(f)+
			 "\": missing boundary",
			 ErrorHandler::normal);
    }
  }

  real_t& integral = __realValue;
  integral = 0;

  switch (M.type()) {
  case Mesh::cartesianHexahedraMesh: {
    const Structured3DMesh& m = static_cast<const Structured3DMesh&>(M);
    integral = __integrate(m, f);
    break;
  }
  case Mesh::hexahedraMesh: {
    const MeshOfHexahedra& m = static_cast<const MeshOfHexahedra&>(M);
    integral = __integrate(m, f);
    break;
  }
  case Mesh::tetrahedraMesh: {
    const MeshOfTetrahedra& m = static_cast<const MeshOfTetrahedra&>(M);
    integral = __integrate(m, f);
    break;
  }
  case Mesh::surfaceMeshTriangles: {
    const SurfaceMeshOfTriangles& m
      = static_cast<const SurfaceMeshOfTriangles&>(M);
    integral = __integrate(m, f);
    break;
  }
  case Mesh::surfaceMeshQuadrangles: {
    const SurfaceMeshOfQuadrangles& m
      = static_cast<const SurfaceMeshOfQuadrangles&>(M);
    integral = __integrate(m, f);
    break;
  }
  default: {
    throw ErrorHandler(__FILE__,__LINE__,
		       "unexpected mesh type",
		       ErrorHandler::unexpected);
  }
  }

  Information::instance().unsetMesh();
}

RealExpressionIntegrate::RealExpressionIntegrate(ReferenceCounting<FunctionExpression> f,
						 ReferenceCounting<MeshExpression> m,
						 const DiscretizationType::Type& discretizationType)
  : __realFunction(f),
    __mesh(m),
    __discretizationType(discretizationType)
{
  ;
}

RealExpressionIntegrate::RealExpressionIntegrate(const RealExpressionIntegrate& e)
  : __realFunction(e.__realFunction),
    __mesh(e.__mesh),
    __discretizationType(e.__discretizationType)
{
}

RealExpressionIntegrate::~RealExpressionIntegrate()
{
  ;
}


/* Max */

ReferenceCounting<RealExpression> RealExpressionMax::value()
{
  return new RealExpressionValue(__realValue);
}

void RealExpressionMax::execute()
{
  (*__mesh).execute();
  Information::instance().setMesh(__mesh);

  const Mesh& M = (*(*__mesh).mesh());

  (*__realFunction).execute();
  FunctionExpression& f = (*__realFunction);

  real_t& max = __realValue;
  max = std::numeric_limits<real_t>::min();

  for (size_t i=0; i<M.numberOfVertices(); ++i) {
    const Vertex& X = M.vertex(i);
    const real_t value = f.value(X[0], X[1], X[2]);
    max = (max>value) ? max : value;
  }

  Information::instance().unsetMesh();
}

RealExpressionMax::RealExpressionMax(ReferenceCounting<FunctionExpression> f,
				     ReferenceCounting<MeshExpression> m)
  : __realFunction(f),
    __mesh(m)
{
  ;
}

RealExpressionMax::RealExpressionMax(const RealExpressionMax& e)
  : __realFunction(e.__realFunction),
    __mesh(e.__mesh)
{
}

RealExpressionMax::~RealExpressionMax()
{
  ;
}


/* Boolean */

ReferenceCounting<RealExpression> RealExpressionBoolean::value()
{
  if((*__booleanExpression).boolValue()) {
    return new RealExpressionValue(1.);
  } else {
    return new RealExpressionValue(0.);
  }
}

void RealExpressionBoolean::execute()
{
  (*__booleanExpression).execute();
}

RealExpressionBoolean::
RealExpressionBoolean(ReferenceCounting<BooleanExpression> be)
  : __booleanExpression(be)
{
  ;
}

RealExpressionBoolean::RealExpressionBoolean(const RealExpressionBoolean& re)
  : RealExpression(re),
    __booleanExpression(re.__booleanExpression)
{
  ;
}

RealExpressionBoolean::~RealExpressionBoolean()
{
  ;
}