//  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: BaseFEMDiscretization.hpp,v 1.10 2007/06/07 22:42:02 delpinux Exp $

#ifndef BASE_FEM_DISCRETIZATION_HPP
#define BASE_FEM_DISCRETIZATION_HPP

#include <FiniteElementTraits.hpp>
#include <Discretization.hpp>

#include <DiscretizedOperators.hpp>
#include <DegreeOfFreedomSet.hpp>

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

#include <VariationalProblem.hpp>

#include <VariationalOperatorFV.hpp>
#include <VariationalOperatorFdxGV.hpp>
#include <VariationalOperatorFdxV.hpp>
#include <VariationalOperatorFgradGgradV.hpp>

#include <ConformTransformation.hpp>

#include <ErrorHandler.hpp>

/**
 * @file   BaseFEMDiscretization.hpp
 * @author Stphane Del Pino
 * @date   Mon Apr 14 00:09:03 2003
 * 
 * @brief  Finite element base class
 * 
 * This class is used to define finite element method standard tools
 * [elementary matrices generation...]. It main purpose is to
 * factorise pieces of code.
 */


template <typename GivenMeshType,
	  DiscretizationType::Type TypeOfDiscretization>
class BaseFEMDiscretization
  : public Discretization
{
protected:
  /// The type of mesh used for discretization
  typedef GivenMeshType MeshType;

  /// The geometry of finite elements
  typedef typename MeshType::CellType CellType;

  /// Finite element type
  typedef typename FiniteElementTraits<CellType,
				       TypeOfDiscretization>::Type FiniteElement;

  /// Elementary matrices type
  typedef typename FiniteElement::ElementaryMatrix ElementaryMatrixType;

  /// Elementary vector type
  typedef typename FiniteElement::ElementaryVector ElementaryVectorType;

  /// type of transformation from reference element
  typedef
  typename FiniteElementTraits<CellType,
			       TypeOfDiscretization>::Transformation
  ConformTransformation;

  /// Associated jacobian
  typedef
  typename FiniteElementTraits<CellType,
			       TypeOfDiscretization>::JacobianTransformation
  JacobianTransformation;

  /// Mesh used to perform discretization
  MeshType& __mesh;

  /// Set of elementary matrices
  mutable ElementaryMatrixSet <ElementaryMatrixType> __eSet;

  /// Operators that are discretized
  mutable DiscretizedOperators<ElementaryMatrixType> __discretizedOperators;

  /// Set of degrees of freedom
  const DegreeOfFreedomSet& __degreeOfFreedomSet;

  /** 
   * Generates elementary vector
   * 
   * @param eVector the generated elementary vector
   * @param J the jacobian of the transformation
   * @param f the function to discretize
   */
  void
  generatesElementaryVector(ElementaryVectorType& eVector,
			    const JacobianTransformation& J,
			    const TinyVector<FiniteElement::numberOfQuadraturePoints,real_t>& f) const
  {
    FiniteElement::instance().integrateWj(eVector,J,f);
  }

  /** 
   * Generates elementary matrices set for a given element
   * 
   * @param eSet the set of elementary matrices
   * @param J the jacobian of the transformation
   */
  void
  generatesElementaryMatrix(ElementaryMatrixSet<ElementaryMatrixType>& eSet,
			    const JacobianTransformation& J) const
  {
    if (eSet.isMassOperator()) {
      generatesElementaryMatrix(PDEOperator::massop,
				J, eSet.massOperator());
    }
    if (eSet.isFirstOrderOperator()) {
      for (size_t i=0; i<3; ++i) {
	if (eSet.isFirstOrderUdxV(i)) {
	  generatesElementaryMatrix(PDEOperator::firstorderopTransposed,
				    J,eSet.firstOrderOperatorUdxV(i),i);
	}
	if (eSet.isFirstOrderDxUV(i)) {
	  generatesElementaryMatrix(PDEOperator::firstorderop, J,
				    eSet.firstOrderOperatorDxUV(i),i);
	}
      }
    }
    if (eSet.isSecondOrderOperator()) {
      for (size_t i=0; i<3; ++i)
	for (size_t j=0; j<3; ++j) {
	  if (eSet.isSecondOrderOperator(i,j)) {
	    generatesElementaryMatrix(PDEOperator::secondorderop, J,
				      eSet.secondOrderOperator(i,j),i,j);
	  }
	}
    }
    if (eSet.isDivMuGrad()) {
      generatesElementaryMatrix(PDEOperator::divmugrad, J, eSet.divMuGrad());
    }
  }

  /** 
   * Generates an elementary matrix for a given operator in an
   * element. The row and column number can be specified when operator
   * is not scalar: \f$ \partial_{x_i}(w_l)\partial_{x_j}(w_k)\f$ for instance.
   * 
   * @param operatorType type of the operator
   * @param J jacobian of the transformation
   * @param matelem generated elementary matrix
   * @param i row number
   * @param j column number
   */
  void
  generatesElementaryMatrix(const PDEOperator::Type operatorType,
			    const JacobianTransformation& J,
			    ElementaryMatrixType& matelem,
			    const size_t i = 0, const size_t j = 0) const
  {
    matelem = 0;
    switch(operatorType) {
    case PDEOperator::firstorderop: {
      FiniteElement::instance().integrateDWjWi(matelem,i,J);
      matelem *= J.jacobianDet();
      break;
    }
    case PDEOperator::firstorderopTransposed: {
      FiniteElement::instance().integrateWjDWi(matelem,i,J);
      matelem *= J.jacobianDet();
      break;
    }
    case PDEOperator::divmugrad: {
      FiniteElement::instance().integrateDWjDWi(matelem,0,0,J);
      FiniteElement::instance().integrateDWjDWi(matelem,1,1,J);
      FiniteElement::instance().integrateDWjDWi(matelem,2,2,J);
      matelem *= J.jacobianDet();
      break;
    }
    case PDEOperator::secondorderop: {
      FiniteElement::instance().integrateDWjDWi(matelem,i,j,J);
      matelem *= J.jacobianDet();
      break;
    }
    case PDEOperator::massop: {
      FiniteElement::instance().integrateWjWi(matelem,J);
      matelem *= J.jacobianDet();
      break;
    }
    default: {
      throw ErrorHandler(__FILE__,__LINE__,
			 "unable to built elementary matrix for this operator",
			 ErrorHandler::unexpected);
    }
    }
  }

  /** 
   * Constructor of the discretization
   * 
   * @param p the problem
   * @param m the mesh used for discretization
   * @param a matrix storing discretization
   * @param bb vector that stores second member discretization
   * @param dof degrees of freedom set
   * 
   */
  BaseFEMDiscretization(const Problem& p,
			MeshType& m,
			BaseMatrix& a,
			BaseVector& bb,
			const DegreeOfFreedomSet& dof)
    : Discretization(Discretization::FEM, p, a, bb),
      __mesh(m),
      __eSet(problem()),
      __discretizedOperators(__eSet,problem()),
      __degreeOfFreedomSet(dof)
  {
    ;
  }

  /** 
   * virtual destructor
   * 
   */
  virtual ~BaseFEMDiscretization()
  {
    ;
  }
};

#endif // BASE_FEM_DISCRETIZATION_HPP