//  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: ScalarFunctionIntegrate.cpp,v 1.5 2006/09/12 22:34:39 delpinux Exp $

#include <ScalarFunctionIntegrate.hpp>
#include <ErrorHandler.hpp>

#include <Structured3DMesh.hpp>
#include <ConnectivityBuilder.hpp>
#include <ConformTransformation.hpp>

std::ostream& 
ScalarFunctionIntegrate::
__put(std::ostream& os) const
{
  switch(__direction) {
  case ScalarFunctionIntegrate::x: {
    os << "int[x](";
    break;
  }
  case ScalarFunctionIntegrate::y: {
    os << "int[y](";
    break;
  }
  case ScalarFunctionIntegrate::z: {
    os << "int[z](";
    break;
  }
  default: {
    throw ErrorHandler(__FILE__,__LINE__,
		       "unknown integrate direction",
		       ErrorHandler::unexpected);
  }
  }
  os << *__functionToIntegrate << ')';
  return os;
}

template <typename MeshType>
real_t
ScalarFunctionIntegrate::
__evaluate(const TinyVector<3,real_t>& X) const
{
  throw ErrorHandler(__FILE__,__LINE__,
		     "mesh type not supported",
		     ErrorHandler::unexpected);
  return 0;
}

template <>
real_t
ScalarFunctionIntegrate::
__evaluate<Structured3DMesh>(const TinyVector<3,real_t>& X) const
{
  const Mesh& baseMesh = *__functionToIntegrate->baseMesh();
  const Structured3DMesh& mesh = static_cast<const Structured3DMesh&>(baseMesh);

  const real_t lowerBound = (*__lowerBound)(X);
  const real_t upperBound = (*__upperBound)(X);

  const real_t minBound = std::min(lowerBound,upperBound);
  const real_t maxBound = std::max(lowerBound,upperBound);

  TinyVector<3, real_t> a = X;
  TinyVector<3, real_t> b = X;

  a[__direction] = std::max(mesh.shape().a()[__direction], minBound);
  b[__direction] = std::min(mesh.shape().b()[__direction], maxBound);

  Structured3DMesh::const_iterator h = mesh.find(a);

  if (h.end()) {
    return 0; // we are outside the mesh
  }

  const Connectivity<Structured3DMesh>& connectivity = mesh.connectivity();
  if (not(connectivity.hasCellToCells())) {
    ConnectivityBuilder<Structured3DMesh> c(mesh);
    c.generates(Connectivity<Structured3DMesh>::CellToCells);
  }

  TinyVector<3, real_t> X0 = a;

  TinyVector<3, real_t> Xhat0;
  {
    ConformTransformationQ1CartesianHexahedron T0(*h);
    T0.invertT(X0, Xhat0);
  }

  const TinyVector<3,real_t> u = b-a;
  real_t dt0 = 1;

  const size_t faceNumberConverter[6] = {4,2,1,3,0,5};

  real_t integrale = 0;

  do {
    const CartesianHexahedron& H = (*h);
    ConformTransformationQ1CartesianHexahedron T(H);

    a=X0;

    T.value(Xhat0,X0);

    integrale += (*__functionToIntegrate)(0.5*(X0+a))*(X0[__direction]-a[__direction]);

    TinyVector<3, real_t> X1 = X0+dt0*u;

    TinyVector<3, real_t> Xhat1;
    T.invertT(b, Xhat1);

    if((   Xhat1[0] >= 0)
       and(Xhat1[1] >= 0)
       and(Xhat1[2] >= 0)
       and(Xhat1[0] <= 1)
       and(Xhat1[1] <= 1)
       and(Xhat1[2] <= 1)) {

      dt0 = 0;
      integrale += (*__functionToIntegrate)(0.5*(X0+X1))*Norm(X1-X0);
    } else {

      const TinyVector<3, real_t> deltaX = Xhat1-Xhat0;
	
      // Get the planes which are required to compute intersection
      TinyVector<3, real_t> x0 = 0;
      for (size_t i=0; i<3; ++i) {
	if (deltaX[i]>0) x0[i] = 1;
      }

      // Computes the linear abcisse required to reach faces
      TinyVector<3, real_t> coef = std::numeric_limits<real_t>::max();
      for (size_t i=0; i<3; ++i) {
	if (std::abs(deltaX[i])>0) {
	  coef[i] = std::abs((x0[i]-Xhat0[i])/deltaX[i]);
	}
      }

      size_t minimumCoefficentNumber = 0;
      for (size_t i=1; i<3; ++i) {
	minimumCoefficentNumber
	  = (coef[minimumCoefficentNumber]<coef[i])?minimumCoefficentNumber:i;
      }

      const size_t outGoingFace
	= faceNumberConverter[2*minimumCoefficentNumber + ((deltaX[minimumCoefficentNumber]>0)?1:0)];

      const real_t dt = dt0*coef[minimumCoefficentNumber];

      Xhat0 += coef[minimumCoefficentNumber] * deltaX;
      dt0-=dt;

      h = connectivity.cells(*h)[outGoingFace];
      if (h.end()) {
	throw ErrorHandler(__FILE__,__LINE__,
			   "cannot find next cell",
			   ErrorHandler::unexpected);
      } // takes the value on the border when going outside

      // X0 is not updated. We compute the value of Xhat0 in the next cell
      Xhat0[minimumCoefficentNumber] = 1-x0[minimumCoefficentNumber];
    }
  } while (dt0>0);

  if (lowerBound == minBound) {
    return integrale;
  } else {
    return -integrale;
  }
}

real_t
ScalarFunctionIntegrate::
operator()(const TinyVector<3, real_t>& X) const
{
  const Mesh& baseMesh = *__functionToIntegrate->baseMesh();
  if (baseMesh.type() != Mesh::cartesianHexahedraMesh) {
    throw ErrorHandler(__FILE__,__LINE__,
		       "Can only evaluate integrate functions on cartesian structured meshes",
		       ErrorHandler::normal);
  }

  return this->__evaluate<Structured3DMesh>(X);
}

ScalarFunctionIntegrate::
ScalarFunctionIntegrate(ConstReferenceCounting<ScalarFunctionBase> lowerBound,
			ConstReferenceCounting<ScalarFunctionBase> upperBound,
			ConstReferenceCounting<ScalarFunctionBase> functionToIntegrate,
			const ScalarFunctionIntegrate::Direction& direction)
  : ScalarFunctionBase(ScalarFunctionBase::integrate),
    __lowerBound(lowerBound),
    __upperBound(upperBound),
    __direction(direction)
{
  if (functionToIntegrate->type() != ScalarFunctionBase::femfunction) {
    std::stringstream errorMsg;

    errorMsg << "cannot compute 1D integrale of non FEM functions :-(\n";
    errorMsg << "the function " << (*functionToIntegrate) << " is not of that kind"
	     << std::ends;

    throw ErrorHandler(__FILE__,__LINE__,
		       errorMsg.str(),
		       ErrorHandler::normal);

  }
  const ScalarFunctionBase* f = functionToIntegrate;
  __functionToIntegrate = dynamic_cast<const FEMFunctionBase*>(f);
}

ScalarFunctionIntegrate::
ScalarFunctionIntegrate(const ScalarFunctionIntegrate& f)
  : ScalarFunctionBase(f),
    __lowerBound(f.__lowerBound),
    __upperBound(f.__upperBound),
    __functionToIntegrate(f.__functionToIntegrate),
    __direction(f.__direction)
{
  ;
}
    
ScalarFunctionIntegrate::
~ScalarFunctionIntegrate()
{
  ;
}