// ===========================================================================
//
//       Filename:  FEM4EllipticPDE.h
//
//        Created:  08/16/12 18:19:41
//       Compiler:  Tested with g++, icpc, and MSVC 2010
//
//         Author:  Trevor Irons (ti)
//
//   Organisation:  Colorado School of Mines (CSM)
//                  United States Geological Survey (USGS)
//
//          Email:  tirons@mines.edu, tirons@usgs.gov
//
//  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 3 of the License, 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, see <http://www.gnu.org/licenses/>.
//
// ===========================================================================

/**
  @file
  @author   Trevor Irons
  @date     08/16/12
  @version   $Rev$
 **/

#ifndef  FEM4ELLIPTICPDE_INC
#define  FEM4ELLIPTICPDE_INC

#include "LemmaObject.h"
#include "vtkImplicitFunction.h"
#include "vtkDataSet.h"
#include "vtkXMLDataSetWriter.h"
#include "vtkSmartPointer.h"
#include "vtkIdList.h"
#include "vtkCell.h"
#include "vtkDoubleArray.h"
#include "vtkPointData.h"
#include "vtkDataSetSurfaceFilter.h"
#include "vtkCellArray.h"
#include "vtkCellData.h"
#include <vtkXMLUnstructuredGridReader.h>
#include <vtkUnstructuredGrid.h>

#include <Eigen/IterativeLinearSolvers>
#include <Eigen/SparseLU>
#include "DCSurvey.h"

namespace Lemma {

    /** @defgroup FEM4EllipticPDE
     @brief Finite element solver for elliptic PDE's
     @details Uses Galerkin method.
    */

    /** \addtogroup FEM4EllipticPDE
     @{
     */

    // ===================================================================
    //  Class:  FEM4EllipticPDE
    /**
      @class FEM4EllipticPDE
      @brief Galerkin FEM solver for Elliptic PDE problems
      @details Solves problems of the type
        \f[ - \nabla  \cdot ( \sigma(\mathbf{r}) \nabla u(\mathbf{r}) ) = g(\mathbf{r})
        \f]
     */
    class FEM4EllipticPDE : public LemmaObject {

        friend std::ostream &operator<<(std::ostream &stream,
			const FEM4EllipticPDE  &ob);

        public:

            // ====================  LIFECYCLE     =======================

            /** Returns a pointer to a new object of type FEM4EllipticPDE.
             * It allocates all necessary memory.
             */
            static std::shared_ptr<FEM4EllipticPDE> NewSP();

            /**
             * Default locked constructor.
             */
            explicit FEM4EllipticPDE (const ctor_key& key);

            /** Protected DeDerializing constructor, use factory DeSerialize  method*/
            FEM4EllipticPDE (const YAML::Node& node, const ctor_key& key);

            /** Default protected constructor. */
            ~FEM4EllipticPDE ();

            /**
             *  Uses YAML to serialize this object.
             *  @return a YAML::Node
             */
            YAML::Node Serialize() const;

            /**
             *   Constructs an object from a YAML::Node.
             */
            static std::shared_ptr<FEM4EllipticPDE> DeSerialize(const YAML::Node& node);

            // ====================  OPERATORS     =======================

            // ====================  OPERATIONS    =======================

            /** Sets the vtkImplictFunction that is used as sigma in
            \f$ - \nabla  \cdot ( \sigma(\mathbf{r}) \nabla u(\mathbf{r}) ) = g(\mathbf{r})
            \f$. If this is not set, sigma evaluates as 1.
             */
            void SetSigmaFunction(vtkImplicitFunction* sigma);

            /** Sets the step to be used in applying boundary conditions.
             */
            void SetBoundaryStep(const Real& h);

            /** Sets the vtkImplictFunction that is used
             */
            void SetGFunction(vtkImplicitFunction* G);

            /** Uses a function pointer instead of a vtkImplicitFunction. For functions that can
             *  be written in closed form, this can be much faster and more accurate then interpolating
             *  off a grid.
             */
            void SetGFunction( Real (*pFcn3)(const Real&, const Real&, const Real&) );

            /* Sets the vtkDataSet that defines the calculation grid.
             */
            //void SetGrid(vtkDataSet* Grid);

            /** Sets the vtkDataSet that defines the calculation grid.
             */
            void SetGrid(vtkUnstructuredGrid* Grid);

            /**
             *  Read grid in from VTK file
             */
            void SetVTUGridFile( const std::string& vtkGridFile );

            /**
             *  Read grid in from VTK file
             */
            void SetVTUGridFile( const char* vtkGridFile );

            /** Sets up a DC problem with a survey
             *  @param[in] ij is the injection index
             */
            void SetupDC(DCSurvey* Survey, const int& ij);

            /** Sets up the potential problem from the VTK file
             */
            void SetupPotential();

            /** Performs solution */
            void Solve(const std::string& fname);

            /* Performs solution */
            //void SolveOLD2(const std::string& fname);

            // ====================  ACCESS        =======================

            // ====================  INQUIRY       =======================

            /** Returns the name of the underlying class, similiar to Python's type */
            virtual std::string GetName() const {
                return this->CName;
            }

        protected:

            // ====================  LIFECYCLE     =======================

            // ====================  OPERATIONS    =======================

            /** Used internally to construct stiffness matrix, A
             */
            void ConstructAMatrix();

            /* Used internally to construct the load vector, g
             */
            //void ConstructLoadVector();

            // This function is copyright (C) 2012 Joseph Capriotti
            /**
                Returns a vtkIdList of points in member vtkGrid connected to point ido.
                @param[in] ido is the point of interest
                @param[in] grid is a pointer to a vtkDataSet
            */
            vtkSmartPointer<vtkIdList> GetConnectedPoints(const int& id0);

            /** Uses Simpson's rule to numerically integrate
             * (Shamelessly adapted from http://en.wikipedia.org/wiki/Simpson's_rule
             */
            Real CompositeSimpsons(vtkImplicitFunction* f, Real r0[3], Real r1[3], int numInt);

            /** Uses Simpson's rule to numerically integrate
             * (Shamelessly adapted from http://en.wikipedia.org/wiki/Simpson's_rule
             *  this method is for a static f
             */
            Real CompositeSimpsons(const Real& f, Real r0[3], Real r1[3], int numInt);

            /**
             *  Uses Simpson's rule to numerically integrate two functions in 1 dimension over the points
             *   r0 and r1
             */
            Real CompositeSimpsons2(vtkImplicitFunction* f, Real r0[3], Real r1[3], int numInt);

            /**
             *  Uses Simpson's rule to numerically integrate two functions in 1 dimension over the points
             *   r0 and r1, uses the Hat function and the function pointer fPtr
             */
            Real CompositeSimpsons2(  Real (*fPtr)(const Real&, const Real&, const Real&), Real r0[3], Real r1[3], int numInt);

            /**
             *  Uses Simpson's rule to numerically integrate two functions in 1 dimension over the points
             *   r0 and r1, uses the Hat function. Assumes a constand function value f
             */
            Real CompositeSimpsons2( const  Real& f, Real r0[3], Real r1[3], int numInt);

            /**
             *  Uses Simpson's rule to numerically integrate two functions in 1 dimension over the points
             *   r0 and r1, uses the Hat function. Assumes a linear function from f0 to f1.
             */
            Real CompositeSimpsons2( const  Real& f0, const Real& f1, Real r0[3], Real r1[3], int numInt);

            /**
             *  Uses Simpson's rule to numerically integrate three functions in 1 dimension over the points
             *   r0 and r1, uses two Hat functions. Assumes a linear function from f0 to f1.
             */
            Real CompositeSimpsons3( const  Real& f0, const Real& f1, Real r0[3], Real r1[3], int numInt);

            /** Hat function for use in Galerkin FEM method, actually a 1/2 Hat.
             */
            //Real hat(Real r[3], Real ri[3], Real rm1[3] );
            Real Hat(const Vector3r &r, const Vector3r& ri, const Vector3r& rm1 );

            /** Computes distance between r0 and r1
             */
            Real dist(Real r0[3], Real r1[3]);

            /** Computes distance between r0 and r1
             */
            Real dist(const Vector3r& r0, const Vector3r& r1);

            /** Sets up the potential problem from the VTK file
             */
            void SetupLineSourcePotential();

            /** Sets up the potential problem from the VTK file
             */
            void SetupSurfaceSourcePotential();

            /** Sets up the potential problem from the VTK file
             */
            void SetupPointSourcePotential();

            /** Sets up the potential problem from the VTK file
             */
            void SetupVolumeSourcePotential();

            // ====================  DATA MEMBERS  =========================

            /** The effective step at the boundary condition. Defaults to 1
             */
            Real BndryH;

            /** The sigma value at the boundary condition. Defaults to 1
             */
            Real BndrySigma;

            /** Implicit function defining the \f$\sigma\f$ term in
                \f$ - \nabla  \cdot ( \sigma(\mathbf{r}) \nabla u(\mathbf{r}) ) = g(\mathbf{r})
                \f$.
                If set to 1, the Poisson equation is solved. Any type of vtkImplicitFunction may
                be used, including vtkImplicitDataSet.
            */
            vtkSmartPointer<vtkImplicitFunction>    vtkSigma;

            /** Implicit function defining the \f$g\f$ term in
                \f$ - \nabla  \cdot ( \sigma(\mathbf{r}) \nabla u(\mathbf{r}) ) = g(\mathbf{r})
                \f$.
                Typically this is used to define the source term.
            */
            vtkSmartPointer<vtkImplicitFunction>    vtkG;

            /** Any type of vtkDataSet may be used as a calculation Grid */
            vtkSmartPointer<vtkUnstructuredGrid>             vtkGrid;

            /** Function pointer to function describing g */
            Real (*gFcn3)(const Real&, const Real&, const Real&);

            Eigen::SparseMatrix<Real> A;

            VectorXr                  g;

        private:

            static constexpr auto CName = "FEM4EllipticPDE";

    }; // -----  end of class  FEM4EllipticPDE  -----

    /*! @} */

}		// -----  end of Lemma  name  -----

#endif   // ----- #ifndef FEM4ELLIPTICPDE_INC  -----