
#include <TNL/Solvers/Solver.h>
#include <TNL/Solvers/BuildConfigTags.h>
#include <TNL/Operators/DirichletBoundaryConditions.h>
#include <TNL/Operators/NeumannBoundaryConditions.h>
#include <TNL/Functions/Analytic/Constant.h>
#include "{problemBaseName}Problem.h"
#include "{operatorName}.h"
#include "{problemBaseName}Rhs.h"
#include "{problemBaseName}BuildConfigTag.h"

using namespace TNL;

typedef {problemBaseName}BuildConfigTag BuildConfig;
    
/****    
 * Uncomment the following (and comment the previous line) for the complete build.    
 * This will include support for all floating point precisions, all indexing types    
 * and more solvers. You may then choose between them from the command line.    
 * The compile time may, however, take tens of minutes or even several hours,    
 * especially if CUDA is enabled. Use this, if you want, only for the final build,        
 * not in the development phase.    
 */    
//typedef TNL::Solvers::DefaultConfigTag BuildConfig;    

template< typename ConfigTag >
class {problemBaseName}Config
{{
   public:
      static void configSetup( Config::ConfigDescription& config )
      {{
         config.addDelimiter( "{problemName} settings:" );
         config.addEntry< String >( "boundary-conditions-type", "Choose the boundary conditions type.", "dirichlet");
            config.addEntryEnum< String >( "dirichlet" );
            config.addEntryEnum< String >( "neumann" );
         config.addEntry< double >( "boundary-conditions-constant", "This sets a value in case of the constant boundary conditions." );

         /****
          * Add definition of your solver command line arguments.
          */

      }}
}};

template< typename Real,
          typename Device,
          typename Index,
          typename MeshType,
          typename ConfigTag,
          typename SolverStarter,
          typename Communicator >
class {problemBaseName}Setter
{{
   public:

      typedef Real RealType;
      typedef Device DeviceType;
      typedef Index IndexType;

      static bool run( const Config::ParameterContainer & parameters )
      {{
          enum {{ Dimension = MeshType::getMeshDimension() }};
          typedef {operatorName}< MeshType, Real, Index > ApproximateOperator;
          typedef {problemBaseName}Rhs< MeshType, Real > RightHandSide;    
          typedef Containers::StaticVector < MeshType::getMeshDimension(), Real > Vertex;

         /****
          * Resolve the template arguments of your solver here.
          * The following code is for the Dirichlet and the Neumann boundary conditions.
          * Both can be constant or defined as discrete functions as Functions::MeshFunction.    
          */    
          String boundaryConditionsType = parameters.getParameter< String >( "boundary-conditions-type" );
          if( parameters.checkParameter( "boundary-conditions-constant" ) )
          {{
             typedef Functions::Analytic::Constant< Dimension, Real > ConstantFunction;
             if( boundaryConditionsType == "dirichlet" )
             {{
                typedef Operators::DirichletBoundaryConditions< MeshType, ConstantFunction, MeshType::getMeshDimension(), Real, Index > BoundaryConditions;
                typedef {problemBaseName}Problem< MeshType, Communicator, BoundaryConditions, RightHandSide, ApproximateOperator > Problem;
                SolverStarter solverStarter;
                return solverStarter.template run< Problem >( parameters );
             }}
             typedef Operators::NeumannBoundaryConditions< MeshType, ConstantFunction, Real, Index > BoundaryConditions;
             typedef {problemBaseName}Problem< MeshType, Communicator, BoundaryConditions, RightHandSide, ApproximateOperator > Problem;
             SolverStarter solverStarter;
             return solverStarter.template run< Problem >( parameters );
          }}
          typedef Functions::MeshFunction< MeshType > MeshFunction;
          if( boundaryConditionsType == "dirichlet" )
          {{    
             typedef Operators::DirichletBoundaryConditions< MeshType, MeshFunction, MeshType::getMeshDimension(), Real, Index > BoundaryConditions;
             typedef {problemBaseName}Problem< MeshType, Communicator, BoundaryConditions, RightHandSide, ApproximateOperator > Problem;
             SolverStarter solverStarter;
             return solverStarter.template run< Problem >( parameters );
          }}
          typedef Operators::NeumannBoundaryConditions< MeshType, MeshFunction, Real, Index > BoundaryConditions;
          typedef {problemBaseName}Problem< MeshType, Communicator, BoundaryConditions, RightHandSide, ApproximateOperator > Problem;
          SolverStarter solverStarter;
          return solverStarter.template run< Problem >( parameters );
      }}
}};

int main( int argc, char* argv[] )
{{
   Solvers::Solver< {problemBaseName}Setter, {problemBaseName}Config, BuildConfig > solver;
   if( ! solver. run( argc, argv ) )
      return EXIT_FAILURE;
   return EXIT_SUCCESS;
}}
