//  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.2 2004/01/24 13:33:39 delpinux Exp $

#ifndef DISCRETIZEDOPERATORS_HPP
#define DISCRETIZEDOPERATORS_HPP

#include <TinyMatrix.hpp>

#include <ReferenceCounting.hpp>

#include <Problem.hpp>
#include <UserFunction.hpp>

#include <map>

#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>

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;
    const size_t __j;
    ConstReferenceCounting<UserFunction> __function;
  public:
    inline const size_t i() const
    {
      return __i;
    }

    inline const size_t j() const
    {
      return __j;
    }

    inline const real_t operator()(const TinyVector<3>& X) const
    {
      return (*__function)(X);
    }

    FunctionAndPosition(const FunctionAndPosition& fap)
      : __i(fap.__i),
	__j(fap.__j),
	__function(fap.__function)
    {
      ;
    }

    FunctionAndPosition(const size_t line, const size_t column,
			ConstReferenceCounting<UserFunction> function)
      : __i(line),
	__j(column),
	__function(function)
    {
      ;
    }

    ~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;

  typename DiscretizedOperators::iterator begin()
  {
    return __list.begin();
  }

  const typename DiscretizedOperators::const_iterator begin() const
  {
    return __list.begin();
  }

  typename DiscretizedOperators::iterator end()
  {
    return __list.end();
  }

  const typename DiscretizedOperators::const_iterator end() const
  {
    return __list.end();
  }

private:
  const ElementaryMatrixSet<ElementaryMatrixType>& __elementaryMatrixSet;

  List __list;

public:
  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.nbOperators(); ++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) {
		__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: {
	      fferr(2) << '\n' << __FILE__ << ':' << __LINE__
		       << ':' << "Not implemented\n";
	      std::exit(1);
	    }
	    }
	  }
	}
      }
      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.number();
	  __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.number();
	  __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: {
	  fferr(2) << __FILE__ << ':' << __LINE__ << ": Not implemented\n";
	  std::exit(1);
	}
	}
      }
      break;
    }
    default: {
      fferr(2) << __FILE__ << ':' << __LINE__ << ": Not implemented\n";
      std::exit(1);
    }
    }
  }
};


#endif