#ifndef DUNE_PLANAR_ROD_ASSEMBLER_HH
#define DUNE_PLANAR_ROD_ASSEMBLER_HH

#include <dune/common/fmatrix.hh>

#include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/matrixindexset.hh>
#include <dune/istl/matrix.hh>

namespace Dune 
{

    /** \brief The FEM operator for an extensible, shearable rod
     */
    template <class GridType, int polOrd>
    class PlanarRodAssembler {
        
        typedef typename GridType::template Codim<0>::Entity EntityType;
        typedef typename GridType::template Codim<0>::LevelIterator ElementIterator;

        //! Dimension of the grid.  This needs to be one!
        enum { gridDim = GridType::dimension };

        enum { elementOrder = 1};

        //! Each block is x, y, theta
        enum { blocksize = 3 };
        
        //!
        typedef FieldMatrix<double, blocksize, blocksize> MatrixBlock;
        
        const GridType* grid_; 
        
        /** \brief Material constants */
        double B;
        double A1;
        double A3;

    public:
        
        //! ???
        PlanarRodAssembler(const GridType &grid) : 
            grid_(&grid)
        { 
            B = 1;
            A1 = 1;
            A3 = 1;
        }

        ~PlanarRodAssembler() {}

        void setParameters(double b, double a1, double a3) {
            B  = b;
            A1 = a1;
            A3 = a3;
        }

        /** \brief Assemble the tangent stiffness matrix and the right hand side
         */
        void assembleMatrix(const BlockVector<FieldVector<double, blocksize> >& sol,
                            BCRSMatrix<MatrixBlock>& matrix);
        
        void assembleGradient(const BlockVector<FieldVector<double, blocksize> >& sol,
                              BlockVector<FieldVector<double, blocksize> >& grad) const;

        /** \brief Compute the energy of a deformation state */
        double computeEnergy(const BlockVector<FieldVector<double, blocksize> >& sol) const;

        void getNeighborsPerVertex(MatrixIndexSet& nb) const;
        
    protected:

        /** \brief Compute the element tangent stiffness matrix  */
        template <class MatrixType>
        void getLocalMatrix( EntityType &entity, 
                             const BlockVector<FieldVector<double, blocksize> >& localSolution, 
                             const int matSize, MatrixType& mat) const;

        
        
        
        
        
        
    }; // end class
    
} // end namespace 

#include "planarrodassembler.cc"

#endif