//  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: DiscretizedOperators.hpp,v 1.6 2006/07/20 19:08:54 delpinux Exp $

#ifndef DISCRETIZED_OPERATORS_HPP
#define DISCRETIZED_OPERATORS_HPP

#include <TinyMatrix.hpp>

#include <ReferenceCounting.hpp>

#include <Problem.hpp>
#include <ScalarFunctionBase.hpp>

#include <ElementaryMatrixSet.hpp>

#include <PDESystem.hpp>

#include <MassOperator.hpp>
#include <FirstOrderOperator.hpp>
#include <DivMuGrad.hpp>
#include <SecondOrderOperator.hpp>


#include <VariationalProblem.hpp>
#include <VariationalOperatorMuGradUGradV.hpp>
#include <VariationalOperatorAlphaDxUDxV.hpp>
#include <VariationalOperatorNuUdxV.hpp>
#include <VariationalOperatorNuDxUV.hpp>
#include <VariationalOperatorAlphaUV.hpp>

#include <ErrorHandler.hpp>

#include <map>

/**
 * @file   DiscretizedOperators.hpp
 * @author Stephane Del Pino
 * @date   Wed Jul 19 15:57:01 2006
 * 
 * @brief Builds lists of discretized operators according to a problem
 * definition
 * 
 */
template <typename ElementaryMatrixType>
class DiscretizedOperators
{
public:
  /*! This sub-class is used to describe the position of the parameter
    function associated to a PDE Operator in the System.
   */
  class FunctionAndPosition
  {
  private:
    const size_t __i;		/**< line of the function in the linear system */
    const size_t __j;		/**< column of the function in the linear system */
    ConstReferenceCounting<ScalarFunctionBase>
    __function;			/**< coefficient function */
  public:
    /** 
     * Read-only access to the line of the function
     * 
     * @return __i
     */
    inline const size_t& i() const
    {
      return __i;
    }

    /** 
     * Read-only access to the column of the function
     * 
     * @return __j
     */
    inline const size_t& j() const
    {
      return __j;
    }

    /** 
     * Evaluates the function at the point @f$ x @f$
     * 
     * @return @f$ f(x) @f$
     */
    inline real_t operator()(const TinyVector<3, real_t>& x) const
    {
      return (*__function)(x);
    }

    /** 
     * Copy constructor
     * 
     * @param fap given FunctionAndPosition
     */
    FunctionAndPosition(const FunctionAndPosition& fap)
      : __i(fap.__i),
	__j(fap.__j),
	__function(fap.__function)
    {
      ;
    }

    /** 
     * Constructor
     * 
     * @param line line in the linear system
     * @param column column in the linear system
     * @param function given function
     */
    FunctionAndPosition(const size_t& line,
			const size_t& column,
			ConstReferenceCounting<ScalarFunctionBase> function)
      : __i(line),
	__j(column),
	__function(function)
    {
      ;
    }

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

  // ================================
  // Now we describe the class itself 
  // ================================

  typedef std::pair<const ElementaryMatrixType*,
		    FunctionAndPosition > ListPair;
  typedef std::multimap<const ElementaryMatrixType*,
			FunctionAndPosition > List;
  typedef typename List::iterator iterator;
  typedef typename List::const_iterator const_iterator;

  /**Begining of the list 
   * 
   * @return begin iterator
   */
  typename DiscretizedOperators::iterator
  begin()
  {
    return __list.begin();
  }

  /** 
   * Begining of the list
   * 
   * @return begin const_iterator
   */
  const typename DiscretizedOperators::const_iterator
  begin() const
  {
    return __list.begin();
  }

  /** 
   * End of the list
   * 
   * @return end iterator
   */
  typename DiscretizedOperators::iterator
  end()
  {
    return __list.end();
  }

  /** 
   * End of the list
   * 
   * @return end const_iterator
   */
  const typename DiscretizedOperators::const_iterator
  end() const
  {
    return __list.end();
  }

private:
  const ElementaryMatrixSet<ElementaryMatrixType>&
  __elementaryMatrixSet;	/**< Reference to the elementary
				   matrices set */

  List __list;			/**< list of discretized operators */

  /** 
   * Copy constructor
   * @note this constructor is forbidden
   * 
   * @param c given disretized operator
   */
  DiscretizedOperators(const DiscretizedOperators& c);

public:
  /** 
   * Constructor
   * 
   * @param e elementary matrices set
   * @param problem problem to discretize
   */
  DiscretizedOperators(const ElementaryMatrixSet<ElementaryMatrixType>& e,
		       const Problem& problem)
    : __elementaryMatrixSet(e)
  {
    switch (problem.type()) {
    case Problem::pdeProblem: {
      const PDESystem& pdeSystem
	= dynamic_cast<const PDESystem&>(problem);

      for (size_t i=0; i<problem.numberOfUnknown(); ++i) {
	const PDE& pde = pdeSystem[i].pde();
	for (size_t j=0; j<pdeSystem.numberOfEquations(); ++j) {
	  const PDEOperatorSum& pdeOpSum = *pde[j];
	  for (size_t k=0; k<pdeOpSum.numberOfOperators(); ++k) {
	    const PDEOperator& pdeOperator = *pdeOpSum[k];
	    switch (pdeOperator.type()) {
	    case PDEOperator::firstorderop: {
	      const FirstOrderOperator& firstOrderOperator
		= dynamic_cast<const FirstOrderOperator&>(pdeOperator);
	      for (size_t m=0; m<3; ++m) {
		if (firstOrderOperator.isSet(m)) {
		  __list.insert(ListPair(&__elementaryMatrixSet.firstOrderOperatorDxUV(m),
					 FunctionAndPosition(i,j,firstOrderOperator.nu(m))));
		}
	      }
	      break;
	    }
	    case PDEOperator::divmugrad: {
	      const DivMuGrad& divMuGrad
		= dynamic_cast<const DivMuGrad&>(pdeOperator);
	      __list.insert(ListPair(&__elementaryMatrixSet.divMuGrad(),
				     FunctionAndPosition(i,j,divMuGrad.mu())));
	      break;
	    }
	    case PDEOperator::secondorderop: {
	      const SecondOrderOperator& secondOrderOperator
		= dynamic_cast<const SecondOrderOperator&>(pdeOperator);
	      for (size_t m=0; m<3; ++m)
		for (size_t n=0; n<3; ++n) {
		  if (secondOrderOperator.isSet(m,n)) {
		    __list.insert(ListPair(&__elementaryMatrixSet.secondOrderOperator(m,n),
					   FunctionAndPosition(i,j,secondOrderOperator.A(m,n))));
		  }
		}
	      break;
	    }
	    case PDEOperator::massop: {
	      const MassOperator& massOperator
		= dynamic_cast<const MassOperator&>(pdeOperator);
	      __list.insert(ListPair(&__elementaryMatrixSet.massOperator(),
				     FunctionAndPosition(i,j,massOperator.alpha())));
	      break;
	    }
	    default: {
	      throw ErrorHandler(__FILE__,__LINE__,
				 "unknown operator",
				 ErrorHandler::unexpected);
	    }
	    }
	  }
	}
      }
      break;
    }
    case Problem::variationalProblem: {
      const VariationalProblem& P
	= dynamic_cast<const VariationalProblem&>(problem);
      for (VariationalProblem::bilinearOperatorConst_iterator
	     i = P.beginBilinearOperator();
	   i != P.endBilinearOperator(); ++i) {
	switch ((*(*i)).type()) {
	case VariationalBilinearOperator::muGradUGradV: {
	  const VariationalMuGradUGradVOperator& O
	    = dynamic_cast<const VariationalMuGradUGradVOperator&>(*(*i));

	  __list.insert(ListPair(&__elementaryMatrixSet.divMuGrad(),
				 FunctionAndPosition(O.unknownNumber(),
						     O.testFunctionNumber(),
						     O.mu())));

	  break;
	}
	case VariationalBilinearOperator::alphaDxUDxV: {
	  const VariationalAlphaDxUDxVOperator& O
	    = dynamic_cast<const VariationalAlphaDxUDxVOperator&>(*(*i));

	  const size_t i = O.i();
	  const size_t j = O.j();
	
	  __list.insert(ListPair(&__elementaryMatrixSet.secondOrderOperator(i,j),
				 FunctionAndPosition(O.unknownNumber(),
						     O.testFunctionNumber(),
						     O.alpha())));
	  break;
	}
	case VariationalBilinearOperator::nuUdxV: {
	  const VariationalNuUdxVOperator& O
	    = dynamic_cast<const VariationalNuUdxVOperator&>(*(*i));

	  const size_t n = O.i();
	  __list.insert(ListPair(&__elementaryMatrixSet.firstOrderOperatorUdxV(n),
				 FunctionAndPosition(O.unknownNumber(),
						     O.testFunctionNumber(),
						     O.nu())));
	  break;
	}
	case VariationalBilinearOperator::nuDxUV: {
	  const VariationalNuDxUVOperator& O
	    = dynamic_cast<const VariationalNuDxUVOperator&>(*(*i));

	  const size_t n = O.i();
	  __list.insert(ListPair(&__elementaryMatrixSet.firstOrderOperatorDxUV(n),
				 FunctionAndPosition(O.unknownNumber(),
						     O.testFunctionNumber(),
						     O.nu())));
	  break;
	}
	case VariationalBilinearOperator::alphaUV: {
	  const VariationalAlphaUVOperator& O
	    = dynamic_cast<const VariationalAlphaUVOperator&>(*(*i));

	  __list.insert(ListPair(&__elementaryMatrixSet.massOperator(),
				 FunctionAndPosition(O.unknownNumber(),
						     O.testFunctionNumber(),
						     O.alpha())));

	  break;
	}
	default: {
	  throw ErrorHandler(__FILE__,__LINE__,
			     "unknown variational form",
			     ErrorHandler::unexpected);
	}
	}
      }
      break;
    }
    default: {
      throw ErrorHandler(__FILE__,__LINE__,
			 "Unknown problem type",
			 ErrorHandler::unexpected);
    }
    }
  }
};

#endif // DISCRETIZED_OPERATORS_HPP