pdelab.hh

Go to the documentation of this file.

// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
#ifndef DUNE_PDELAB_HH
#define DUNE_PDELAB_HH

// first of all we include a lot of dune grids and pdelab files
#include <iostream>

#include <dune/common/parallel/mpihelper.hh> // include mpi helper class
#include <dune/common/parametertreeparser.hh>
#include <dune/common/shared_ptr.hh>
#include <dune/common/classname.hh>
#include <dune/common/exceptions.hh>
#include <dune/common/fvector.hh>

#include <dune/geometry/type.hh>
#include <dune/geometry/quadraturerules.hh>

#include <dune/grid/sgrid.hh> // load sgrid definition
#include <dune/grid/onedgrid.hh>
#include <dune/grid/io/file/vtk/subsamplingvtkwriter.hh>
#include <dune/grid/yaspgrid.hh>
#if HAVE_UG
#include <dune/grid/uggrid.hh>
#endif
#if HAVE_ALBERTA
#include<dune/grid/albertagrid.hh>
#include <dune/grid/albertagrid/dgfparser.hh>
#endif
#if HAVE_UG
#include<dune/grid/uggrid.hh>
#endif
#if HAVE_ALUGRID
#include<dune/grid/alugrid.hh>
#include<dune/grid/io/file/dgfparser/dgfalu.hh>
#include<dune/grid/io/file/dgfparser/dgfparser.hh>
#endif
#include <dune/grid/utility/structuredgridfactory.hh>
#include <dune/grid/io/file/gmshreader.hh>

#include <dune/istl/bvector.hh>
#include <dune/istl/operators.hh>
#include <dune/istl/solvers.hh>
#include <dune/istl/solvercategory.hh>
#include <dune/istl/preconditioners.hh>
#include <dune/istl/io.hh>

#include <dune/istl/paamg/amg.hh>
#include <dune/pdelab/common/function.hh>
#include <dune/pdelab/common/functionutilities.hh>
#include <dune/pdelab/common/vtkexport.hh>
#include <dune/pdelab/backend/istlvectorbackend.hh>
#include <dune/pdelab/backend/istl/bcrsmatrixbackend.hh>
#include <dune/pdelab/backend/istlsolverbackend.hh>
#include <dune/pdelab/constraints/conforming.hh>
#include <dune/pdelab/constraints/hangingnode.hh>
#include <dune/pdelab/constraints/p0.hh>
#include <dune/pdelab/constraints/p0ghost.hh>
#include <dune/pdelab/constraints/common/constraints.hh>
#include <dune/pdelab/gridfunctionspace/gridfunctionspace.hh>
#include <dune/pdelab/gridfunctionspace/gridfunctionspaceutilities.hh>
#include <dune/pdelab/gridfunctionspace/interpolate.hh>
#include <dune/pdelab/gridoperator/gridoperator.hh>
#include <dune/pdelab/gridoperator/onestep.hh>
#include <dune/pdelab/stationary/linearproblem.hh>
#include <dune/pdelab/finiteelementmap/pkfem.hh>
#include <dune/pdelab/finiteelementmap/p0fem.hh>
#include <dune/pdelab/finiteelementmap/opbfem.hh>
#include <dune/pdelab/finiteelementmap/qkfem.hh>
#include <dune/pdelab/finiteelementmap/qkdg.hh>
#include <dune/pdelab/finiteelementmap/qkdggl.hh>
#include <dune/pdelab/adaptivity/adaptivity.hh>
#include <dune/pdelab/instationary/onestep.hh>
#include <dune/pdelab/newton/newton.hh>

namespace Dune {
    namespace PDELab {

        // make grids
        template<typename T>
        class StructuredGrid
        {
        public:
            // export types
            typedef T Grid;
            typedef typename T::ctype ctype;
            static const int dim = T::dimension;
            static const int dimworld = T::dimensionworld;

            // constructors
            StructuredGrid (Dune::GeometryType::BasicType meshtype, unsigned int cells)
            {
                FieldVector<ctype,dimworld> lowerLeft(0.0);
                FieldVector<ctype,dimworld> upperRight(1.0);
                array<unsigned int,dim> elements; elements.fill(cells);

                StructuredGridFactory<T> factory;

                if (meshtype==Dune::GeometryType::cube)
                    gridp = factory.createCubeGrid(lowerLeft,upperRight,elements);
                else if (meshtype==Dune::GeometryType::simplex)
                    gridp = factory.createSimplexGrid(lowerLeft,upperRight,elements);
                else
                    {
                        DUNE_THROW(GridError, className<StructuredGrid>()
                                   << "::StructuredGrid(): grid type must be simplex or cube ");
                    }
            }


            StructuredGrid (Dune::GeometryType::BasicType meshtype,
                            array<double,dimworld> lower_left, array<double,dimworld> upper_right,
                            array<unsigned int,dim> cells)
            {
                FieldVector<ctype,dimworld> lowerLeft;
                FieldVector<ctype,dimworld> upperRight;
                array<unsigned int,dim> elements;

                // copy data to correct types for StructuredGridFactory
                for (size_t i=0; i<dimworld; i++)
                    {
                        lowerLeft[i] = lower_left[i];
                        upperRight[i] = upper_right[i];
                    }
                for (size_t i=0; i<dim; i++)
                    {
                        elements[i] = cells[i];
                    }

                StructuredGridFactory<T> factory;

                if (meshtype==Dune::GeometryType::cube)
                    gridp = factory.createCubeGrid(lowerLeft,upperRight,elements);
                else if (meshtype==Dune::GeometryType::simplex)
                    gridp = factory.createSimplexGrid(lowerLeft,upperRight,elements);
                else
                    {
                        DUNE_THROW(GridError, className<StructuredGrid>()
                                   << "::StructuredGrid(): grid type must be simplex or cube ");
                    }
            }

            // return shared pointer
            Dune::shared_ptr<T> getSharedPtr ()
            {
                return gridp;
            }

            // return grid reference
            T& getGrid ()
            {
                return *gridp;
            }

            // return grid reference const version
            const T& getGrid () const
            {
                return *gridp;
            }

            T& operator*()
            {
                return *gridp;
            }

            T* operator->()
            {
                return gridp.operator->();
            }

            const T& operator*() const
            {
                return *gridp;
            }

            const T* operator->() const
            {
                return gridp.operator->();
            }


        private:
            Dune::shared_ptr<T> gridp; // hold a shared pointer to a grid
        };

        // specialization for yaspgrid; treats paralle case right
        template<int dim>
        class StructuredGrid<YaspGrid<dim> >
        {
        public:

            // export types
            typedef YaspGrid<dim> Grid;
            typedef typename Grid::ctype ctype;
            static const int dimworld = Grid::dimensionworld;

            // simple constructor for the unit cube
            StructuredGrid (Dune::GeometryType::BasicType meshtype, unsigned int cells, int overlap=1)
            {
                // check element type
                if (meshtype!=Dune::GeometryType::cube)
                    std::cout << "StructuredGrid(): element type " << meshtype << " is ignored" << std::endl;

                // copy data to correct types for YaspGrid
                Dune::FieldVector<double,dimworld> L(1.0);
                Dune::array<int,dimworld> N(Dune::fill_array<int,dimworld>(cells));
                std::bitset<dimworld> B(false);

                // instantiate the grid
#if HAVE_MPI
                gridp = shared_ptr<Grid>(new Grid(Dune::MPIHelper::getCommunicator(),L,N,B,overlap));
#else
                gridp = shared_ptr<Grid>(new Grid(L,N,B,overlap));
#endif
            }

            // constructor with sizes given
            StructuredGrid (Dune::GeometryType::BasicType meshtype,
                            array<double,dimworld> lower_left, array<double,dimworld> upper_right,
                            array<unsigned int,dim> cells, int overlap=1)
            {
                // check that lower right corner is the origin
                for(int d = 0; d < dimworld; ++d)
                    if(std::abs(lower_left[d]) > std::abs(upper_right[d])*1e-10)
                        DUNE_THROW(GridError, className<StructuredGrid>()
                                   << "::createCubeGrid(): The lower coordinates "
                                   "must be at the origin for YaspGrid.");

                // check element type
                if (meshtype!=Dune::GeometryType::cube)
                    std::cout << "StructuredGrid(): element type " << meshtype << " is ignored" << std::endl;

                // copy data to correct types for YaspGrid
                Dune::FieldVector<double,dimworld> L;
                Dune::array<int,dimworld> N;
                std::bitset<dimworld> B(false);
                for (size_t i=0; i<dimworld; i++)
                    {
                        L[i] = upper_right[i];
                        N[i] = cells[i];
                    }

                // instantiate the grid
#if HAVE_MPI
                gridp = shared_ptr<Grid>(new Grid(Dune::MPIHelper::getCommunicator(),L,N,B,overlap));
#else
                gridp = shared_ptr<Grid>(new Grid(L,N,B,overlap));
#endif
            }

            // constructor with periodicity argument
            StructuredGrid (Dune::GeometryType::BasicType meshtype,
                            array<double,dimworld> lower_left, array<double,dimworld> upper_right,
                            array<unsigned int,dim> cells, array<bool,dim> periodic, int overlap=1)
            {
                // check that lower right corner is the origin
                for(int d = 0; d < dimworld; ++d)
                    if(std::abs(lower_left[d]) > std::abs(upper_right[d])*1e-10)
                        DUNE_THROW(GridError, className<StructuredGrid>()
                                   << "::createCubeGrid(): The lower coordinates "
                                   "must be at the origin for YaspGrid.");

                // check element type
                if (meshtype!=Dune::GeometryType::cube)
                    std::cout << "StructuredGrid(): element type " << meshtype << " is ignored" << std::endl;

                // copy data to correct types for YaspGrid
                Dune::FieldVector<double,dimworld> L;
                Dune::array<int,dimworld> N;
                std::bitset<dimworld> B(false);
                for (size_t i=0; i<dimworld; i++)
                    {
                        L[i] = upper_right[i];
                        N[i] = cells[i];
                        B[i] = periodic[i];
                    }

                // instantiate the grid
#if HAVE_MPI
                gridp = shared_ptr<Grid>(new Grid(Dune::MPIHelper::getCommunicator(),L,N,B,overlap));
#else
                gridp = shared_ptr<Grid>(new Grid(L,N,B,overlap));
#endif
            }

            // return shared pointer
            Dune::shared_ptr<Grid> getSharedPtr ()
            {
                return gridp;
            }

            // return grid reference
            Grid& getGrid ()
            {
                return *gridp;
            }

            // return grid reference const version
            const Grid& getGrid () const
            {
                return *gridp;
            }

            Grid& operator*()
            {
                return *gridp;
            }

            Grid* operator->()
            {
                return gridp.operator->();
            }

            const Grid& operator*() const
            {
                return *gridp;
            }

            const Grid* operator->() const
            {
                return gridp.operator->();
            }

        private:
            Dune::shared_ptr<Grid> gridp; // hold a shared pointer to a grid
        };

        // unstructured grid read from gmsh file
        template<typename T>
        class UnstructuredGrid
        {
        public:
            // export types
            typedef T Grid;
            typedef typename T::ctype ctype;
            static const int dim = T::dimension;
            static const int dimworld = T::dimensionworld;

            // constructors
            UnstructuredGrid (std::string filename, bool verbose = true, bool insert_boundary_segments=true)
            {
                Dune::GridFactory<T> factory;
                Dune::GmshReader<T>::read(factory,filename,verbose,insert_boundary_segments);
                gridp = shared_ptr<T>(factory.createGrid());
            }

            // return shared pointer
            Dune::shared_ptr<T> getSharedPtr ()
            {
                return gridp;
            }

            // return grid reference
            T& getGrid ()
            {
                return *gridp;
            }

            // return grid reference const version
            const T& getGrid () const
            {
                return *gridp;
            }

            T& operator*()
            {
                return *gridp;
            }

            T* operator->()
            {
                return gridp.operator->();
            }

            const T& operator*() const
            {
                return *gridp;
            }

            const T* operator->() const
            {
                return gridp.operator->();
            }

        private:
            Dune::shared_ptr<T> gridp; // hold a shared pointer to a grid
        };


        //============================================================================
        // Continuous Lagrange Finite Element Space
        //============================================================================

        // finite element map base template
        template<typename GV, typename C, typename R, unsigned int degree, unsigned int dim, Dune::GeometryType::BasicType gt>
        class CGFEMBase
        {};

        template<typename GV, typename C, typename R, unsigned int degree, unsigned int dim>
        class CGFEMBase<GV,C,R,degree,dim,Dune::GeometryType::simplex>
        {
        public:
            typedef PkLocalFiniteElementMap<GV,C,R,degree> FEM;

            CGFEMBase (const GV& gridview)
            {
                femp = shared_ptr<FEM>(new FEM(gridview));
            }

            FEM& getFEM() {return *femp;}
            const FEM& getFEM() const {return *femp;}

        private:
            shared_ptr<FEM> femp;
        };

        template<typename GV, typename C, typename R, unsigned int degree, unsigned int dim>
        class CGFEMBase<GV,C,R,degree,dim,Dune::GeometryType::cube>
        {
        public:
            typedef QkLocalFiniteElementMap<GV,C,R,degree> FEM;

            CGFEMBase (const GV& gridview)
            {
                femp = shared_ptr<FEM>(new FEM(gridview));
            }

            FEM& getFEM() {return *femp;}
            const FEM& getFEM() const {return *femp;}

        private:
            shared_ptr<FEM> femp;
        };

        //============================================================================

        // define enumeration type that differentiate conforming and nonconforming meshes
        enum MeshType {
            conforming,
            nonconforming
        };

        // constraints base template
        template<typename Grid, unsigned int degree, Dune::GeometryType::BasicType gt, MeshType mt, SolverCategory::Category st, typename BCType, typename GV = typename Grid::LeafGridView>
        class CGCONBase
        {};

        template<typename Grid, typename BCType, typename GV>
        class CGCONBase<Grid,1,Dune::GeometryType::simplex,MeshType::nonconforming,SolverCategory::sequential,BCType,GV>
        {
        public:
            typedef HangingNodesDirichletConstraints<Grid,HangingNodesConstraintsAssemblers::SimplexGridP1Assembler,BCType> CON;

            CGCONBase (Grid& grid, const BCType& bctype, const GV& gv)
            {
                conp = shared_ptr<CON>(new CON(grid,true,bctype));
            }

            CGCONBase (Grid& grid, const BCType& bctype)
            {
                conp = shared_ptr<CON>(new CON(grid,true,bctype));
            }

            template<typename GFS>
            void postGFSHook (const GFS& gfs) {}
            CON& getCON() {return *conp;}
            const CON& getCON() const {return *conp;}
            template<typename GFS, typename DOF>
            void make_consistent (const GFS& gfs, DOF& x) const {}
        private:
            shared_ptr<CON> conp;
        };

        template<typename Grid, typename BCType, typename GV>
        class CGCONBase<Grid,1,Dune::GeometryType::cube,MeshType::nonconforming,SolverCategory::sequential,BCType,GV>
        {
        public:
            typedef HangingNodesDirichletConstraints<Grid,HangingNodesConstraintsAssemblers::CubeGridQ1Assembler,BCType> CON;

            CGCONBase (Grid& grid, const BCType& bctype, const GV& gv)
            {
                conp = shared_ptr<CON>(new CON(grid,true,bctype));
            }

            CGCONBase (Grid& grid, const BCType& bctype)
            {
                conp = shared_ptr<CON>(new CON(grid,true,bctype));
            }

            template<typename GFS>
            void postGFSHook (const GFS& gfs) {}
            CON& getCON() {return *conp;}
            const CON& getCON() const {return *conp;}
            template<typename GFS, typename DOF>
            void make_consistent (const GFS& gfs, DOF& x) const {}
        private:
            shared_ptr<CON> conp;
        };

        template<typename Grid, unsigned int degree, Dune::GeometryType::BasicType gt,typename BCType, typename GV>
        class CGCONBase<Grid,degree,gt,MeshType::conforming,SolverCategory::sequential,BCType,GV>
        {
        public:
            typedef ConformingDirichletConstraints CON;

            CGCONBase (Grid& grid, const BCType& bctype, const GV& gv)
            {
                conp = shared_ptr<CON>(new CON());
            }

            CGCONBase (Grid& grid, const BCType& bctype)
            {
                conp = shared_ptr<CON>(new CON());
            }

            template<typename GFS>
            void postGFSHook (const GFS& gfs) {}
            CON& getCON() {return *conp;}
            const CON& getCON() const {return *conp;}
            template<typename GFS, typename DOF>
            void make_consistent (const GFS& gfs, DOF& x) const {}
        private:
            shared_ptr<CON> conp;
        };

        template<typename Grid, unsigned int degree, Dune::GeometryType::BasicType gt,typename BCType, typename GV>
        class CGCONBase<Grid,degree,gt,MeshType::conforming,SolverCategory::overlapping,BCType,GV>
        {
        public:
            typedef OverlappingConformingDirichletConstraints CON;

            CGCONBase (Grid& grid, const BCType& bctype, const GV& gv)
            {
                conp = shared_ptr<CON>(new CON());
            }

            CGCONBase (Grid& grid, const BCType& bctype)
            {
                conp = shared_ptr<CON>(new CON());
            }

            template<typename GFS>
            void postGFSHook (const GFS& gfs) {}
            CON& getCON() {return *conp;}
            const CON& getCON() const {return *conp;}
            template<typename GFS, typename DOF>
            void make_consistent (const GFS& gfs, DOF& x) const
            {
                // make vector consistent; this is needed for all overlapping solvers
                istl::ParallelHelper<GFS> helper(gfs);
                helper.maskForeignDOFs(istl::raw(x));
                Dune::PDELab::AddDataHandle<GFS,DOF> adddh(gfs,x);
                if (gfs.gridView().comm().size()>1)
                    gfs.gridView().communicate(adddh,Dune::InteriorBorder_All_Interface,Dune::ForwardCommunication);
            }
        private:
            shared_ptr<CON> conp;
        };

        template<typename Grid, unsigned int degree, Dune::GeometryType::BasicType gt,typename BCType, typename GV>
        class CGCONBase<Grid,degree,gt,MeshType::conforming,SolverCategory::nonoverlapping,BCType,GV>
        {
        public:
            typedef NonoverlappingConformingDirichletConstraints<GV> CON;

            CGCONBase (Grid& grid, const BCType& bctype, const GV& gv)
            {
                conp = shared_ptr<CON>(new CON(gv));
            }

            CGCONBase (Grid& grid, const BCType& bctype)
            {
                conp = shared_ptr<CON>(new CON(grid.leafGridView()));
            }

            template<typename GFS>
            void postGFSHook (const GFS& gfs) { conp->compute_ghosts(gfs); }
            CON& getCON() {return *conp;}
            const CON& getCON() const {return *conp;}
            template<typename GFS, typename DOF>
            void make_consistent (const GFS& gfs, DOF& x) const {}
        private:
            shared_ptr<CON> conp;
        };


        // continuous Lagrange finite elements
        template<typename T, typename N, unsigned int degree, typename BCType,
                 Dune::GeometryType::BasicType gt, MeshType mt, SolverCategory::Category st = SolverCategory::sequential,
                 typename VBET=ISTLVectorBackend<> >
        class CGSpace {
        public:

            // export types
            typedef T Grid;
            typedef typename T::LeafGridView GV;
            typedef typename T::ctype ctype;
            static const int dim = T::dimension;
            static const int dimworld = T::dimensionworld;

            typedef CGFEMBase<GV,ctype,N,degree,dim,gt> FEMB;
            typedef CGCONBase<Grid,degree,gt,mt,st,BCType> CONB;

            typedef typename FEMB::FEM FEM;
            typedef typename CONB::CON CON;

            typedef VBET VBE;
            typedef GridFunctionSpace<GV,FEM,CON,VBE> GFS;

            typedef N NT;
            typedef typename BackendVectorSelector<GFS,N>::Type DOF;
            typedef Dune::PDELab::DiscreteGridFunction<GFS,DOF> DGF;
            typedef typename GFS::template ConstraintsContainer<N>::Type CC;
            typedef VTKGridFunctionAdapter<DGF> VTKF;

            // constructor making the grid function space an all that is needed
            CGSpace (Grid& grid, const BCType& bctype)
                : gv(grid.leafGridView()), femb(gv), conb(grid,bctype)
            {
                gfsp = shared_ptr<GFS>(new GFS(gv,femb.getFEM(),conb.getCON()));
                gfsp->name("cgspace");
                // initialize ordering
                gfsp->update();
                conb.postGFSHook(*gfsp);
                ccp = shared_ptr<CC>(new CC());
            }

            FEM& getFEM()
            {
                return femb.getFEM();
            }

            const FEM& getFEM() const
            {
                return femb.getFEM();
            }

            // return gfs reference
            GFS& getGFS ()
            {
                return *gfsp;
            }

            // return gfs reference const version
            const GFS& getGFS () const
            {
                return *gfsp;
            }

            // return gfs reference
            CC& getCC ()
            {
                return *ccp;
            }

            // return gfs reference const version
            const CC& getCC () const
            {
                return *ccp;
            }

            void assembleConstraints (const BCType& bctype)
            {
                ccp->clear();
                constraints(bctype,*gfsp,*ccp);
            }

            void clearConstraints ()
            {
                ccp->clear();
            }

            void setConstrainedDOFS (DOF& x, NT nt) const
            {
                set_constrained_dofs(*ccp,nt,x);
                conb.make_consistent(*gfsp,x);
            }

            void setNonConstrainedDOFS (DOF& x, NT nt) const
            {
                set_nonconstrained_dofs(*ccp,nt,x);
                conb.make_consistent(*gfsp,x);
            }

            void copyConstrainedDOFS (const DOF& xin, DOF& xout) const
            {
                copy_constrained_dofs(*ccp,xin,xout);
                conb.make_consistent(*gfsp,xout);
            }

            void copyNonConstrainedDOFS (const DOF& xin, DOF& xout) const
            {
                copy_nonconstrained_dofs(*ccp,xin,xout);
                conb.make_consistent(*gfsp,xout);
            }

        private:
            GV gv; // need this object here because FEM and GFS store a const reference !!
            FEMB femb;
            CONB conb;
            shared_ptr<GFS> gfsp;
            shared_ptr<CC> ccp;
        };


        //============================================================================
        // Discontinuous Finite Element Space
        //============================================================================

        // constraints base template
        template<SolverCategory::Category st>
        class DGCONBase
        {};

        template<>
        class DGCONBase<SolverCategory::sequential>
        {
        public:
            typedef NoConstraints CON;
            DGCONBase ()
            {
                conp = shared_ptr<CON>(new CON());
            }
            CON& getCON() {return *conp;}
            const CON& getCON() const {return *conp;}
            template<typename GFS, typename DOF>
            void make_consistent (const GFS& gfs, DOF& x) const {}
        private:
            shared_ptr<CON> conp;
        };

        template<>
        class DGCONBase<SolverCategory::nonoverlapping>
        {
        public:
            typedef P0ParallelGhostConstraints CON;
            DGCONBase ()
            {
                conp = shared_ptr<CON>(new CON());
            }
            CON& getCON() {return *conp;}
            const CON& getCON() const {return *conp;}
            template<typename GFS, typename DOF>
            void make_consistent (const GFS& gfs, DOF& x) const {}
        private:
            shared_ptr<CON> conp;
        };

        template<>
        class DGCONBase<SolverCategory::overlapping>
        {
        public:
            typedef P0ParallelConstraints CON;
            DGCONBase ()
            {
                conp = shared_ptr<CON>(new CON());
            }
            CON& getCON() {return *conp;}
            const CON& getCON() const {return *conp;}
            template<typename GFS, typename DOF>
            void make_consistent (const GFS& gfs, DOF& x) const
            {
                // make vector consistent; this is needed for all overlapping solvers
                istl::ParallelHelper<GFS> helper(gfs);
                helper.maskForeignDOFs(istl::raw(x));
                Dune::PDELab::AddDataHandle<GFS,DOF> adddh(gfs,x);
                if (gfs.gridView().comm().size()>1)
                    gfs.gridView().communicate(adddh,Dune::InteriorBorder_All_Interface,Dune::ForwardCommunication);
            }
        private:
            shared_ptr<CON> conp;
        };

        // Discontinuous space
        // default implementation, use only specializations below
        template<typename T, typename N, unsigned int degree,
                 Dune::GeometryType::BasicType gt, SolverCategory::Category st = SolverCategory::sequential,
                 typename VBET=ISTLVectorBackend<ISTLParameters::static_blocking,Dune::PB::PkSize<degree,T::dimension>::value> >
        class DGPkSpace
        {
        public:

            // export types
            typedef T Grid;
            typedef typename T::LeafGridView GV;
            typedef typename T::ctype ctype;
            static const int dim = T::dimension;
            static const int dimworld = T::dimensionworld;
            typedef N NT;
#if HAVE_GMP
            typedef OPBLocalFiniteElementMap<ctype,NT,degree,dim,gt,Dune::GMPField<512>,Dune::PB::BasisType::Pk> FEM;
#else
            typedef OPBLocalFiniteElementMap<ctype,NT,degree,dim,gt> FEM;
#endif
            typedef DGCONBase<st> CONB;
            typedef typename CONB::CON CON;
            typedef VBET VBE;
            typedef GridFunctionSpace<GV,FEM,CON,VBE> GFS;
            typedef typename BackendVectorSelector<GFS,N>::Type DOF;
            typedef Dune::PDELab::DiscreteGridFunction<GFS,DOF> DGF;
            typedef typename GFS::template ConstraintsContainer<N>::Type CC;
            typedef VTKGridFunctionAdapter<DGF> VTKF;

            // constructor making the grid function space an all that is needed
            DGPkSpace (const GV& gridview) : gv(gridview), conb()
            {
                femp = shared_ptr<FEM>(new FEM());
                gfsp = shared_ptr<GFS>(new GFS(gv,*femp));
                // initialize ordering
                gfsp->update();
                ccp = shared_ptr<CC>(new CC());
            }

            FEM& getFEM() { return *femp; }
            const FEM& getFEM() const { return *femp; }

            // return gfs reference
            GFS& getGFS () { return *gfsp; }

            // return gfs reference const version
            const GFS& getGFS () const {return *gfsp;}

            // return gfs reference
            CC& getCC () { return *ccp;}

            // return gfs reference const version
            const CC& getCC () const { return *ccp;}

            template<class BCTYPE>
            void assembleConstraints (const BCTYPE& bctype)
            {
                ccp->clear();
                constraints(bctype,*gfsp,*ccp);
            }

            void clearConstraints ()
            {
                ccp->clear();
            }

            void setConstrainedDOFS (DOF& x, NT nt) const
            {
                set_constrained_dofs(*ccp,nt,x);
                conb.make_consistent(*gfsp,x);
            }

            void setNonConstrainedDOFS (DOF& x, NT nt) const
            {
                set_nonconstrained_dofs(*ccp,nt,x);
                conb.make_consistent(*gfsp,x);
            }

            void copyConstrainedDOFS (const DOF& xin, DOF& xout) const
            {
                copy_constrained_dofs(*ccp,xin,xout);
                conb.make_consistent(*gfsp,xout);
            }

            void copyNonConstrainedDOFS (const DOF& xin, DOF& xout) const
            {
                copy_nonconstrained_dofs(*ccp,xin,xout);
                conb.make_consistent(*gfsp,xout);
            }

        private:
            GV gv; // need this object here because FEM and GFS store a const reference !!
            CONB conb;
            shared_ptr<FEM> femp;
            shared_ptr<GFS> gfsp;
            shared_ptr<CC> ccp;
        };

        // Discontinuous space
        // default implementation, use only specializations below
        template<typename T, typename N, unsigned int degree,
                 Dune::GeometryType::BasicType gt, SolverCategory::Category st = SolverCategory::sequential,
                 //typename VBET=ISTLVectorBackend<ISTLParameters::static_blocking,Dune::PB::PkSize<degree,T::dimension>::value> >
                 typename VBET=ISTLVectorBackend<> >
        class DGQkOPBSpace
        {
        public:

            // export types
            typedef T Grid;
            typedef typename T::LeafGridView GV;
            typedef typename T::ctype ctype;
            static const int dim = T::dimension;
            static const int dimworld = T::dimensionworld;
            typedef N NT;
#if HAVE_GMP
            typedef OPBLocalFiniteElementMap<ctype,NT,degree,dim,gt,Dune::GMPField<512>,Dune::PB::BasisType::Qk> FEM;
#else
            typedef OPBLocalFiniteElementMap<ctype,NT,degree,dim,gt,N,Dune::PB::BasisType::Qk> FEM;
#endif
            typedef DGCONBase<st> CONB;
            typedef typename CONB::CON CON;
            typedef VBET VBE;
            typedef GridFunctionSpace<GV,FEM,CON,VBE> GFS;
            typedef typename BackendVectorSelector<GFS,N>::Type DOF;
            typedef Dune::PDELab::DiscreteGridFunction<GFS,DOF> DGF;
            typedef typename GFS::template ConstraintsContainer<N>::Type CC;
            typedef VTKGridFunctionAdapter<DGF> VTKF;

            // constructor making the grid function space an all that is needed
            DGQkOPBSpace (const GV& gridview) : gv(gridview), conb()
            {
                femp = shared_ptr<FEM>(new FEM());
                gfsp = shared_ptr<GFS>(new GFS(gv,*femp));
                // initialize ordering
                gfsp->update();
                ccp = shared_ptr<CC>(new CC());
            }

            FEM& getFEM() { return *femp; }
            const FEM& getFEM() const { return *femp; }

            // return gfs reference
            GFS& getGFS () { return *gfsp; }

            // return gfs reference const version
            const GFS& getGFS () const {return *gfsp;}

            // return gfs reference
            CC& getCC () { return *ccp;}

            // return gfs reference const version
            const CC& getCC () const { return *ccp;}

            template<class BCTYPE>
            void assembleConstraints (const BCTYPE& bctype)
            {
                ccp->clear();
                constraints(bctype,*gfsp,*ccp);
            }

            void clearConstraints ()
            {
                ccp->clear();
            }

            void setConstrainedDOFS (DOF& x, NT nt) const
            {
                set_constrained_dofs(*ccp,nt,x);
                conb.make_consistent(*gfsp,x);
            }

            void setNonConstrainedDOFS (DOF& x, NT nt) const
            {
                set_nonconstrained_dofs(*ccp,nt,x);
                conb.make_consistent(*gfsp,x);
            }

            void copyConstrainedDOFS (const DOF& xin, DOF& xout) const
            {
                copy_constrained_dofs(*ccp,xin,xout);
                conb.make_consistent(*gfsp,xout);
            }

            void copyNonConstrainedDOFS (const DOF& xin, DOF& xout) const
            {
                copy_nonconstrained_dofs(*ccp,xin,xout);
                conb.make_consistent(*gfsp,xout);
            }

        private:
            GV gv; // need this object here because FEM and GFS store a const reference !!
            CONB conb;
            shared_ptr<FEM> femp;
            shared_ptr<GFS> gfsp;
            shared_ptr<CC> ccp;
        };

        // Discontinuous space
        // default implementation, use only specializations below
        template<typename T, typename N, unsigned int degree,
                 Dune::GeometryType::BasicType gt, SolverCategory::Category st = SolverCategory::sequential,
                 typename VBET=ISTLVectorBackend<ISTLParameters::static_blocking,Dune::QkStuff::QkSize<degree,T::dimension>::value> >
        class DGQkSpace
        {
        public:

            // export types
            typedef T Grid;
            typedef typename T::LeafGridView GV;
            typedef typename T::ctype ctype;
            static const int dim = T::dimension;
            static const int dimworld = T::dimensionworld;
            typedef N NT;
            typedef QkDGLocalFiniteElementMap<ctype,NT,degree,dim> FEM;
            typedef DGCONBase<st> CONB;
            typedef typename CONB::CON CON;
            typedef VBET VBE;
            typedef GridFunctionSpace<GV,FEM,CON,VBE> GFS;
            typedef typename BackendVectorSelector<GFS,N>::Type DOF;
            typedef Dune::PDELab::DiscreteGridFunction<GFS,DOF> DGF;
            typedef typename GFS::template ConstraintsContainer<N>::Type CC;
            typedef VTKGridFunctionAdapter<DGF> VTKF;

            // constructor making the grid function space an all that is needed
            DGQkSpace (const GV& gridview) : gv(gridview), conb()
            {
                femp = shared_ptr<FEM>(new FEM());
                gfsp = shared_ptr<GFS>(new GFS(gv,*femp));
                // initialize ordering
                gfsp->update();
                ccp = shared_ptr<CC>(new CC());
            }

            FEM& getFEM() { return *femp; }
            const FEM& getFEM() const { return *femp; }

            // return gfs reference
            GFS& getGFS () { return *gfsp; }

            // return gfs reference const version
            const GFS& getGFS () const {return *gfsp;}

            // return gfs reference
            CC& getCC () { return *ccp;}

            // return gfs reference const version
            const CC& getCC () const { return *ccp;}

            template<class BCTYPE>
            void assembleConstraints (const BCTYPE& bctype)
            {
                ccp->clear();
                constraints(bctype,*gfsp,*ccp);
            }

            void clearConstraints ()
            {
                ccp->clear();
            }

            void setConstrainedDOFS (DOF& x, NT nt) const
            {
                set_constrained_dofs(*ccp,nt,x);
                conb.make_consistent(*gfsp,x);
            }

            void setNonConstrainedDOFS (DOF& x, NT nt) const
            {
                set_nonconstrained_dofs(*ccp,nt,x);
                conb.make_consistent(*gfsp,x);
            }

            void copyConstrainedDOFS (const DOF& xin, DOF& xout) const
            {
                copy_constrained_dofs(*ccp,xin,xout);
                conb.make_consistent(*gfsp,xout);
            }

            void copyNonConstrainedDOFS (const DOF& xin, DOF& xout) const
            {
                copy_nonconstrained_dofs(*ccp,xin,xout);
                conb.make_consistent(*gfsp,xout);
            }

        private:
            GV gv; // need this object here because FEM and GFS store a const reference !!
            CONB conb;
            shared_ptr<FEM> femp;
            shared_ptr<GFS> gfsp;
            shared_ptr<CC> ccp;
        };


        // Discontinuous space using QK with Gauss Lobatto points (use only for cube elements)
        template<typename T, typename N, unsigned int degree,
                 Dune::GeometryType::BasicType gt, SolverCategory::Category st = SolverCategory::sequential,
                 //typename VBET=ISTLVectorBackend<ISTLParameters::static_blocking,Dune::QkStuff::QkSize<degree,T::dimension>::value> >
                 typename VBET=ISTLVectorBackend<> >
        class DGQkGLSpace
        {
        public:

            // export types
            typedef T Grid;
            typedef typename T::LeafGridView GV;
            typedef typename T::ctype ctype;
            static const int dim = T::dimension;
            static const int dimworld = T::dimensionworld;
            typedef N NT;
            typedef QkDGGLLocalFiniteElementMap<ctype,NT,degree,dim> FEM;
            typedef DGCONBase<st> CONB;
            typedef typename CONB::CON CON;
            typedef VBET VBE;
            typedef GridFunctionSpace<GV,FEM,CON,VBE> GFS;
            typedef typename BackendVectorSelector<GFS,N>::Type DOF;
            typedef Dune::PDELab::DiscreteGridFunction<GFS,DOF> DGF;
            typedef typename GFS::template ConstraintsContainer<N>::Type CC;
            typedef VTKGridFunctionAdapter<DGF> VTKF;

            // constructor making the grid function space an all that is needed
            DGQkGLSpace (const GV& gridview) : gv(gridview), conb()
            {
                femp = shared_ptr<FEM>(new FEM());
                gfsp = shared_ptr<GFS>(new GFS(gv,*femp));
                // initialize ordering
                gfsp->update();
                ccp = shared_ptr<CC>(new CC());
            }

            FEM& getFEM() { return *femp; }
            const FEM& getFEM() const { return *femp; }

            // return gfs reference
            GFS& getGFS () { return *gfsp; }

            // return gfs reference const version
            const GFS& getGFS () const {return *gfsp;}

            // return gfs reference
            CC& getCC () { return *ccp;}

            // return gfs reference const version
            const CC& getCC () const { return *ccp;}

            template<class BCTYPE>
            void assembleConstraints (const BCTYPE& bctype)
            {
                ccp->clear();
                constraints(bctype,*gfsp,*ccp);
            }

            void clearConstraints ()
            {
                ccp->clear();
            }

            void setConstrainedDOFS (DOF& x, NT nt) const
            {
                set_constrained_dofs(*ccp,nt,x);
                conb.make_consistent(*gfsp,x);
            }

            void setNonConstrainedDOFS (DOF& x, NT nt) const
            {
                set_nonconstrained_dofs(*ccp,nt,x);
                conb.make_consistent(*gfsp,x);
            }

            void copyConstrainedDOFS (const DOF& xin, DOF& xout) const
            {
                copy_constrained_dofs(*ccp,xin,xout);
                conb.make_consistent(*gfsp,xout);
            }

            void copyNonConstrainedDOFS (const DOF& xin, DOF& xout) const
            {
                copy_nonconstrained_dofs(*ccp,xin,xout);
                conb.make_consistent(*gfsp,xout);
            }

        private:
            GV gv; // need this object here because FEM and GFS store a const reference !!
            CONB conb;
            shared_ptr<FEM> femp;
            shared_ptr<GFS> gfsp;
            shared_ptr<CC> ccp;
        };



        // Discontinuous P0 space
        template<typename T, typename N,
                 Dune::GeometryType::BasicType gt, SolverCategory::Category st = SolverCategory::sequential,
                 typename VBET=ISTLVectorBackend<> >
        class P0Space
        {
        public:

            // export types
            typedef T Grid;
            typedef typename T::LeafGridView GV;
            typedef typename T::ctype ctype;
            static const int dim = T::dimension;
            static const int dimworld = T::dimensionworld;
            typedef N NT;
            typedef Dune::PDELab::P0LocalFiniteElementMap<ctype,NT,dim> FEM;
            typedef DGCONBase<st> CONB;
            typedef typename CONB::CON CON;
            typedef VBET VBE;
            typedef GridFunctionSpace<GV,FEM,CON,VBE> GFS;
            typedef typename BackendVectorSelector<GFS,N>::Type DOF;
            typedef Dune::PDELab::DiscreteGridFunction<GFS,DOF> DGF;
            typedef typename GFS::template ConstraintsContainer<N>::Type CC;
            typedef VTKGridFunctionAdapter<DGF> VTKF;

            // constructor making the grid function space an all that is needed
            P0Space (const GV& gridview) : gv(gridview), conb()
            {
                femp = shared_ptr<FEM>(new FEM(Dune::GeometryType(gt,dim)));
                gfsp = shared_ptr<GFS>(new GFS(gv,*femp));
                // initialize ordering
                gfsp->update();
                ccp = shared_ptr<CC>(new CC());
            }

            FEM& getFEM() { return *femp; }
            const FEM& getFEM() const { return *femp; }

            // return gfs reference
            GFS& getGFS () { return *gfsp; }

            // return gfs reference const version
            const GFS& getGFS () const {return *gfsp;}

            // return gfs reference
            CC& getCC () { return *ccp;}

            // return gfs reference const version
            const CC& getCC () const { return *ccp;}

            template<class BCTYPE>
            void assembleConstraints (const BCTYPE& bctype)
            {
                ccp->clear();
                constraints(bctype,*gfsp,*ccp);
            }

            void clearConstraints ()
            {
                ccp->clear();
            }

            void setConstrainedDOFS (DOF& x, NT nt) const
            {
                set_constrained_dofs(*ccp,nt,x);
                conb.make_consistent(*gfsp,x);
            }

            void setNonConstrainedDOFS (DOF& x, NT nt) const
            {
                set_nonconstrained_dofs(*ccp,nt,x);
                conb.make_consistent(*gfsp,x);
            }

            void copyConstrainedDOFS (const DOF& xin, DOF& xout) const
            {
                copy_constrained_dofs(*ccp,xin,xout);
                conb.make_consistent(*gfsp,xout);
            }

            void copyNonConstrainedDOFS (const DOF& xin, DOF& xout) const
            {
                copy_nonconstrained_dofs(*ccp,xin,xout);
                conb.make_consistent(*gfsp,xout);
            }

        private:
            GV gv; // need this object here because FEM and GFS store a const reference !!
            CONB conb;
            shared_ptr<FEM> femp;
            shared_ptr<GFS> gfsp;
            shared_ptr<CC> ccp;
        };


        // how can we most easily specify a grid function
        // pass a function space as parameter
        template<typename FS, typename Functor>
        class UserFunction
            : public GridFunctionBase<GridFunctionTraits<typename FS::GV, typename FS::NT,
                                                         1,FieldVector<typename FS::NT,1> >
                                      ,UserFunction<FS,Functor> >
        {
        public:
            typedef GridFunctionTraits<typename FS::GV, typename FS::NT,
                                       1,FieldVector<typename FS::NT,1> > Traits;

            UserFunction (const FS& fs_, const Functor& f_)
                : fs(fs_), f(f_)
            {}

            inline void evaluate (const typename Traits::ElementType& e,
                                  const typename Traits::DomainType& x,
                                  typename Traits::RangeType& y) const
            {
                typename Traits::DomainType x_ = e.geometry().global(x);
                std::vector<double> x__(x.size());
                for (size_t i=0; i<x.size(); ++i) x__[i]=x_[i];
                y = f(x__);
            }

            inline const typename FS::GV& getGridView () const
            {
                return fs.getGFS().gridView();
            }

        private:
            const FS fs; // store a copy of the function space
            const Functor f;
        };


        template<typename FS, typename LOP, SolverCategory::Category st = SolverCategory::sequential>
        class GalerkinGlobalAssembler
        {
        public:
            // export types
            typedef ISTLMatrixBackend MBE;
            typedef Dune::PDELab::GridOperator<typename FS::GFS,typename FS::GFS,LOP,MBE,
                                               typename FS::NT,typename FS::NT,typename FS::NT,
                                               typename FS::CC,typename FS::CC> GO;
            typedef typename GO::Jacobian MAT;

            GalerkinGlobalAssembler (const FS& fs, LOP& lop)
            {
                gop = shared_ptr<GO>(new GO(fs.getGFS(),fs.getCC(),fs.getGFS(),fs.getCC(),lop));
            }

            // return grid reference
            GO& getGO ()
            {
                return *gop;
            }

            // return grid reference const version
            const GO& getGO () const
            {
                return *gop;
            }

            GO& operator*()
            {
                return *gop;
            }

            GO* operator->()
            {
                return gop.operator->();
            }

            const GO& operator*() const
            {
                return *gop;
            }

            const GO* operator->() const
            {
                return gop.operator->();
            }

        private:
            shared_ptr<GO> gop;
        };

        // nonoverlapping variant
        template<typename FS, typename LOP>
        class GalerkinGlobalAssembler<FS,LOP,SolverCategory::nonoverlapping>
        {
        public:
            // export types
            typedef ISTLMatrixBackend MBE;
            typedef Dune::PDELab::GridOperator<typename FS::GFS,typename FS::GFS,LOP,MBE,
                                               typename FS::NT,typename FS::NT,typename FS::NT,
                                               typename FS::CC,typename FS::CC,true> GO;
            typedef typename GO::Jacobian MAT;

            GalerkinGlobalAssembler (const FS& fs, LOP& lop)
            {
                gop = shared_ptr<GO>(new GO(fs.getGFS(),fs.getCC(),fs.getGFS(),fs.getCC(),lop));
            }

            // return grid reference
            GO& getGO ()
            {
                return *gop;
            }

            // return grid reference const version
            const GO& getGO () const
            {
                return *gop;
            }

            GO& operator*()
            {
                return *gop;
            }

            GO* operator->()
            {
                return gop.operator->();
            }

            const GO& operator*() const
            {
                return *gop;
            }

            const GO* operator->() const
            {
                return gop.operator->();
            }

        private:
            shared_ptr<GO> gop;
        };


        template<typename FS, typename LOP, SolverCategory::Category st = SolverCategory::sequential>
        class GalerkinGlobalAssemblerNewBackend
        {
        public:
            // export types
            typedef Dune::PDELab::istl::BCRSMatrixBackend<> MBE;
            typedef Dune::PDELab::GridOperator<typename FS::GFS,typename FS::GFS,LOP,MBE,
                                               typename FS::NT,typename FS::NT,typename FS::NT,
                                               typename FS::CC,typename FS::CC> GO;
            typedef typename GO::Jacobian MAT;

            GalerkinGlobalAssemblerNewBackend (const FS& fs, LOP& lop, const MBE& mbe)
            {
                gop = shared_ptr<GO>(new GO(fs.getGFS(),fs.getCC(),fs.getGFS(),fs.getCC(),lop,mbe));
            }

            // return grid reference
            GO& getGO ()
            {
                return *gop;
            }

            // return grid reference const version
            const GO& getGO () const
            {
                return *gop;
            }

            GO& operator*()
            {
                return *gop;
            }

            GO* operator->()
            {
                return gop.operator->();
            }

            const GO& operator*() const
            {
                return *gop;
            }

            const GO* operator->() const
            {
                return gop.operator->();
            }

        private:
            shared_ptr<GO> gop;
        };

        // nonoverlapping variant
        template<typename FS, typename LOP>
        class GalerkinGlobalAssemblerNewBackend<FS,LOP,SolverCategory::nonoverlapping>
        {
        public:
            // export types
            typedef Dune::PDELab::istl::BCRSMatrixBackend<> MBE;
            typedef Dune::PDELab::GridOperator<typename FS::GFS,typename FS::GFS,LOP,MBE,
                                               typename FS::NT,typename FS::NT,typename FS::NT,
                                               typename FS::CC,typename FS::CC,true> GO;
            typedef typename GO::Jacobian MAT;

            GalerkinGlobalAssemblerNewBackend (const FS& fs, LOP& lop, const MBE& mbe)
            {
                gop = shared_ptr<GO>(new GO(fs.getGFS(),fs.getCC(),fs.getGFS(),fs.getCC(),lop,mbe));
            }

            // return grid reference
            GO& getGO ()
            {
                return *gop;
            }

            // return grid reference const version
            const GO& getGO () const
            {
                return *gop;
            }

            GO& operator*()
            {
                return *gop;
            }

            GO* operator->()
            {
                return gop.operator->();
            }

            const GO& operator*() const
            {
                return *gop;
            }

            const GO* operator->() const
            {
                return gop.operator->();
            }

        private:
            shared_ptr<GO> gop;
        };




        // variant with two different function spaces
        template<typename FSU, typename FSV, typename LOP, SolverCategory::Category st>
        class GlobalAssembler
        {
        public:
            // export types
            typedef ISTLMatrixBackend MBE;
            typedef Dune::PDELab::GridOperator<typename FSU::GFS,typename FSV::GFS,LOP,MBE,
                                               typename FSU::NT,typename FSU::NT,typename FSU::NT,
                                               typename FSU::CC,typename FSV::CC> GO;
            typedef typename GO::Jacobian MAT;

            GlobalAssembler (const FSU& fsu, const FSV& fsv, LOP& lop)
            {
                gop = shared_ptr<GO>(new GO(fsu.getGFS(),fsu.getCC(),fsv.getGFS(),fsv.getCC(),lop));
            }

            // return grid reference
            GO& getGO ()
            {
                return *gop;
            }

            // return grid reference const version
            const GO& getGO () const
            {
                return *gop;
            }

            GO& operator*()
            {
                return *gop;
            }

            GO* operator->()
            {
                return gop.operator->();
            }

            const GO& operator*() const
            {
                return *gop;
            }

            const GO* operator->() const
            {
                return gop.operator->();
            }

        private:
            shared_ptr<GO> gop;
        };

        // nonoverlapping variant
        template<typename FSU, typename FSV, typename LOP>
        class GlobalAssembler<FSU,FSV,LOP,SolverCategory::nonoverlapping>
        {
        public:
            // export types
            typedef ISTLMatrixBackend MBE;
            typedef Dune::PDELab::GridOperator<typename FSU::GFS,typename FSV::GFS,LOP,MBE,
                                               typename FSU::NT,typename FSU::NT,typename FSU::NT,
                                               typename FSU::CC,typename FSV::CC,true> GO;
            typedef typename GO::Jacobian MAT;

            GlobalAssembler (const FSU& fsu, const FSV& fsv, LOP& lop)
            {
                gop = shared_ptr<GO>(new GO(fsu.getGFS(),fsu.getCC(),fsv.getGFS(),fsv.getCC(),lop));
            }

            // return grid reference
            GO& getGO ()
            {
                return *gop;
            }

            // return grid reference const version
            const GO& getGO () const
            {
                return *gop;
            }

            GO& operator*()
            {
                return *gop;
            }

            GO* operator->()
            {
                return gop.operator->();
            }

            const GO& operator*() const
            {
                return *gop;
            }

            const GO* operator->() const
            {
                return gop.operator->();
            }

        private:
            shared_ptr<GO> gop;
        };



        template<typename GO1, typename GO2, bool implicit = true>
        class OneStepGlobalAssembler
        {
        public:
            // export types
            typedef ISTLMatrixBackend MBE;
            typedef Dune::PDELab::OneStepGridOperator<typename GO1::GO,typename GO2::GO,implicit> GO;
            typedef typename GO::Jacobian MAT;

            OneStepGlobalAssembler (GO1& go1, GO2& go2)
            {
                gop = shared_ptr<GO>(new GO(*go1,*go2));
            }

            // return grid reference
            GO& getGO ()
            {
                return *gop;
            }

            // return grid reference const version
            const GO& getGO () const
            {
                return *gop;
            }

            GO& operator*()
            {
                return *gop;
            }

            GO* operator->()
            {
                return gop.operator->();
            }

            const GO& operator*() const
            {
                return *gop;
            }

            const GO* operator->() const
            {
                return gop.operator->();
            }

        private:
            shared_ptr<GO> gop;
        };


        // packaging of the CG_AMG_SSOR solver: default version is sequential
        template<typename FS, typename ASS, SolverCategory::Category st = SolverCategory::sequential>
        class ISTLSolverBackend_CG_AMG_SSOR
        {
        public:
            // types exported
            typedef Dune::PDELab::ISTLBackend_SEQ_CG_AMG_SSOR<typename ASS::GO> LS;

            ISTLSolverBackend_CG_AMG_SSOR (const FS& fs, const ASS& ass, unsigned maxiter_=5000,
                                           int verbose_=1, bool reuse_=false, bool usesuperlu_=true)
            {
                lsp = shared_ptr<LS>(new LS(maxiter_,verbose_,reuse_,usesuperlu_));
            }

            LS& getLS () {return *lsp;}
            const LS& getLS () const { return *lsp;}
            LS& operator*(){return *lsp;}
            LS* operator->() { return lsp.operator->(); }
            const LS& operator*() const{return *lsp;}
            const LS* operator->() const{ return lsp.operator->();}

       private:
            shared_ptr<LS> lsp;
        };

        // packaging of the CG_AMG_SSOR solver: nonoverlapping version
        template<typename FS, typename ASS>
        class ISTLSolverBackend_CG_AMG_SSOR<FS,ASS, SolverCategory::nonoverlapping>
        {
        public:
            // types exported
            typedef Dune::PDELab::ISTLBackend_NOVLP_CG_AMG_SSOR<typename ASS::GO> LS;

            ISTLSolverBackend_CG_AMG_SSOR (const FS& fs, const ASS& ass, unsigned maxiter_=5000,
                                           int verbose_=1, bool reuse_=false, bool usesuperlu_=true)
            {
                lsp = shared_ptr<LS>(new LS(fs.getGFS(),maxiter_,verbose_,reuse_,usesuperlu_));
            }

            LS& getLS () {return *lsp;}
            const LS& getLS () const { return *lsp;}
            LS& operator*(){return *lsp;}
            LS* operator->() { return lsp.operator->(); }
            const LS& operator*() const{return *lsp;}
            const LS* operator->() const{ return lsp.operator->();}

       private:
            shared_ptr<LS> lsp;
        };

        // packaging of the CG_AMG_SSOR solver: overlapping version
        template<typename FS, typename ASS>
        class ISTLSolverBackend_CG_AMG_SSOR<FS,ASS, SolverCategory::overlapping>
        {
        public:
            // types exported
            typedef Dune::PDELab::ISTLBackend_CG_AMG_SSOR<typename ASS::GO> LS;

            ISTLSolverBackend_CG_AMG_SSOR (const FS& fs, const ASS& ass, unsigned maxiter_=5000,
                                           int verbose_=1, bool reuse_=false, bool usesuperlu_=true)
            {
                lsp = shared_ptr<LS>(new LS(fs.getGFS(),maxiter_,verbose_,reuse_,usesuperlu_));
            }

            LS& getLS () {return *lsp;}
            const LS& getLS () const { return *lsp;}
            LS& operator*(){return *lsp;}
            LS* operator->() { return lsp.operator->(); }
            const LS& operator*() const{return *lsp;}
            const LS* operator->() const{ return lsp.operator->();}

        private:
            shared_ptr<LS> lsp;
        };

        // packaging of the CG_SSOR solver: default version is sequential
        template<typename FS, typename ASS, SolverCategory::Category st = SolverCategory::sequential>
        class ISTLSolverBackend_CG_SSOR
        {
        public:
            // types exported
            typedef ISTLBackend_SEQ_CG_SSOR LS;

            ISTLSolverBackend_CG_SSOR (const FS& fs, const ASS& ass, unsigned maxiter_=5000,
                                       int steps_=5, int verbose_=1)
            {
                lsp = shared_ptr<LS>(new LS(maxiter_,verbose_));
            }

            LS& getLS () {return *lsp;}
            const LS& getLS () const { return *lsp;}
            LS& operator*(){return *lsp;}
            LS* operator->() { return lsp.operator->(); }
            const LS& operator*() const{return *lsp;}
            const LS* operator->() const{ return lsp.operator->();}

       private:
            shared_ptr<LS> lsp;
        };

        // packaging of the CG_SSOR solver: nonoverlapping version
        template<typename FS, typename ASS>
        class ISTLSolverBackend_CG_SSOR<FS,ASS,SolverCategory::nonoverlapping>
        {
        public:
            // types exported
            typedef ISTLBackend_NOVLP_CG_SSORk<typename ASS::GO> LS;

            ISTLSolverBackend_CG_SSOR (const FS& fs, const ASS& ass, unsigned maxiter_=5000,
                                       int steps_=5, int verbose_=1)
            {
                lsp = shared_ptr<LS>(new LS(fs.getGFS(),maxiter_,steps_,verbose_));
            }

            LS& getLS () {return *lsp;}
            const LS& getLS () const { return *lsp;}
            LS& operator*(){return *lsp;}
            LS* operator->() { return lsp.operator->(); }
            const LS& operator*() const{return *lsp;}
            const LS* operator->() const{ return lsp.operator->();}

       private:
            shared_ptr<LS> lsp;
        };

        // packaging of the CG_SSOR solver: overlapping version
        template<typename FS, typename ASS>
        class ISTLSolverBackend_CG_SSOR<FS,ASS,SolverCategory::overlapping>
        {
        public:
            // types exported
            typedef ISTLBackend_OVLP_CG_SSORk<typename FS::GFS, typename FS::CC> LS;

            ISTLSolverBackend_CG_SSOR (const FS& fs, const ASS& ass, unsigned maxiter_=5000,
                                       int steps_=5, int verbose_=1)
            {
                lsp = shared_ptr<LS>(new LS(fs.getGFS(),fs.getCC(),maxiter_,steps_,verbose_));
            }

            LS& getLS () {return *lsp;}
            const LS& getLS () const { return *lsp;}
            LS& operator*(){return *lsp;}
            LS* operator->() { return lsp.operator->(); }
            const LS& operator*() const{return *lsp;}
            const LS* operator->() const{ return lsp.operator->();}

       private:
            shared_ptr<LS> lsp;
        };


        // packaging of a default solver that should always work
        // in the sequential case : BCGS SSOR
        template<typename FS, typename ASS, SolverCategory::Category st = SolverCategory::sequential>
        class ISTLSolverBackend_IterativeDefault
        {
        public:
            // types exported
            typedef ISTLBackend_SEQ_BCGS_SSOR LS;

            ISTLSolverBackend_IterativeDefault (const FS& fs, const ASS& ass, unsigned maxiter_=5000, int verbose_=1)
            {
                lsp = shared_ptr<LS>(new LS(maxiter_,verbose_));
            }

            LS& getLS () {return *lsp;}
            const LS& getLS () const { return *lsp;}
            LS& operator*(){return *lsp;}
            LS* operator->() { return lsp.operator->(); }
            const LS& operator*() const{return *lsp;}
            const LS* operator->() const{ return lsp.operator->();}

       private:
            shared_ptr<LS> lsp;
        };

        // in the nonoverlapping case : BCGS Jacobi
        template<typename FS, typename ASS>
        class ISTLSolverBackend_IterativeDefault<FS,ASS,SolverCategory::nonoverlapping>
        {
        public:
            // types exported
            typedef ISTLBackend_NOVLP_BCGS_Jacobi<typename FS::GFS> LS;

            ISTLSolverBackend_IterativeDefault (const FS& fs, const ASS& ass, unsigned maxiter_=5000, int verbose_=1)
            {
                lsp = shared_ptr<LS>(new LS(fs.getGFS(),maxiter_,verbose_));
            }

            LS& getLS () {return *lsp;}
            const LS& getLS () const { return *lsp;}
            LS& operator*(){return *lsp;}
            LS* operator->() { return lsp.operator->(); }
            const LS& operator*() const{return *lsp;}
            const LS* operator->() const{ return lsp.operator->();}

       private:
            shared_ptr<LS> lsp;
        };

        // in the overlapping case : BCGS SSORk
        template<typename FS, typename ASS>
        class ISTLSolverBackend_IterativeDefault<FS,ASS,SolverCategory::overlapping>
        {
        public:
            // types exported
            typedef ISTLBackend_OVLP_BCGS_SSORk<typename FS::GFS, typename FS::CC> LS;

            ISTLSolverBackend_IterativeDefault (const FS& fs, const ASS& ass, unsigned maxiter_=5000, int verbose_=1)
            {
                lsp = shared_ptr<LS>(new LS(fs.getGFS(),fs.getCC(),maxiter_,3,verbose_));
            }

            LS& getLS () {return *lsp;}
            const LS& getLS () const { return *lsp;}
            LS& operator*(){return *lsp;}
            LS* operator->() { return lsp.operator->(); }
            const LS& operator*() const{return *lsp;}
            const LS* operator->() const{ return lsp.operator->();}

       private:
            shared_ptr<LS> lsp;
        };

        // packaging of a default solver that should always work
        // in the sequential case : BCGS SSOR
        template<typename FS, typename ASS, SolverCategory::Category st = SolverCategory::sequential>
        class ISTLSolverBackend_ExplicitDiagonal
        {
        public:
            // types exported
            typedef Dune::PDELab::ISTLBackend_SEQ_ExplicitDiagonal LS;

            ISTLSolverBackend_ExplicitDiagonal (const FS& fs, const ASS& ass, unsigned maxiter_=5000, int verbose_=1)
            {
                lsp = shared_ptr<LS>(new LS());
            }

            LS& getLS () {return *lsp;}
            const LS& getLS () const { return *lsp;}
            LS& operator*(){return *lsp;}
            LS* operator->() { return lsp.operator->(); }
            const LS& operator*() const{return *lsp;}
            const LS* operator->() const{ return lsp.operator->();}

       private:
            shared_ptr<LS> lsp;
        };

        // packaging of a default solver that should always work
        // in the sequential case : BCGS SSOR
        template<typename FS, typename ASS>
        class ISTLSolverBackend_ExplicitDiagonal<FS,ASS,SolverCategory::overlapping>
        {
        public:
            // types exported
            typedef Dune::PDELab::ISTLBackend_OVLP_ExplicitDiagonal<typename FS::GFS> LS;

            ISTLSolverBackend_ExplicitDiagonal (const FS& fs, const ASS& ass, unsigned maxiter_=5000, int verbose_=1)
            {
                lsp = shared_ptr<LS>(new LS(fs.getGFS()));
            }

            LS& getLS () {return *lsp;}
            const LS& getLS () const { return *lsp;}
            LS& operator*(){return *lsp;}
            LS* operator->() { return lsp.operator->(); }
            const LS& operator*() const{return *lsp;}
            const LS* operator->() const{ return lsp.operator->();}

       private:
            shared_ptr<LS> lsp;
        };

        // packaging of a default solver that should always work
        // in the sequential case : BCGS SSOR
        template<typename FS, typename ASS>
        class ISTLSolverBackend_ExplicitDiagonal<FS,ASS,SolverCategory::nonoverlapping>
        {
        public:
            // types exported
            typedef Dune::PDELab::ISTLBackend_NOVLP_ExplicitDiagonal<typename FS::GFS> LS;

            ISTLSolverBackend_ExplicitDiagonal (const FS& fs, const ASS& ass, unsigned maxiter_=5000, int verbose_=1)
            {
                lsp = shared_ptr<LS>(new LS(fs.getGFS()));
            }

            LS& getLS () {return *lsp;}
            const LS& getLS () const { return *lsp;}
            LS& operator*(){return *lsp;}
            LS* operator->() { return lsp.operator->(); }
            const LS& operator*() const{return *lsp;}
            const LS* operator->() const{ return lsp.operator->();}

       private:
            shared_ptr<LS> lsp;
        };


} // end namespace PDELab
    } // end namespace Dune

#endif