HamiltonJacobiProblem_impl.h

/***************************************************************************
                          HamiltonJacobiProblem_impl.h  -  description
                             -------------------
    begin                : Jul 8 , 2014
    copyright            : (C) 2014 by Tomas Sobotik
 ***************************************************************************/

/***************************************************************************
 *                                                                         *
 *   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 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 ***************************************************************************/

#pragma once 

#include <core/mfilename.h>
#include <matrices/tnlMatrixSetter.h>
#include <exception>

template< typename Mesh,
		  typename HamiltonJacobi,
		  typename BoundaryCondition,
		  typename RightHandSide>
String HamiltonJacobiProblem< Mesh,HamiltonJacobi,BoundaryCondition,RightHandSide > :: getType()
{
   return String( "hamiltonJacobiSolver< " ) + Mesh  :: getType() + " >";
}

template< typename Mesh,
		  typename HamiltonJacobi,
		  typename BoundaryCondition,
		  typename RightHandSide>
String HamiltonJacobiProblem< Mesh,HamiltonJacobi,BoundaryCondition,RightHandSide > :: getPrologHeader() const
{
   return String( "Hamilton-Jacobi" );
}

template< typename Mesh,
		  typename HamiltonJacobi,
		  typename BoundaryCondition,
		  typename RightHandSide>
void HamiltonJacobiProblem< Mesh,HamiltonJacobi,BoundaryCondition,RightHandSide > :: writeProlog( tnlLogger& logger,
                                                												 const Config::ParameterContainer& parameters ) const
{
   //logger. WriteParameter< typename String >( "Problem name:", "problem-name", parameters );
   //logger. WriteParameter< String >( "Used scheme:", "scheme", parameters );
}

template< typename Mesh,
		  typename HamiltonJacobi,
		  typename BoundaryCondition,
		  typename RightHandSide>
bool HamiltonJacobiProblem< Mesh,HamiltonJacobi,BoundaryCondition,RightHandSide > :: setup( const Config::ParameterContainer& parameters )
{
   this->boundaryCondition.getFunction().setConstant( 1.0 );
   //this->rightHandSide.getFunction().setConstant( 0.0 );
   //return true;
/*
   const String& problemName = parameters. GetParameter< String >( "problem-name" );

   this->schemeTest = parameters. GetParameter< int >( "scheme-test" );
   this->tested = false;

   const String& meshFile = parameters.GetParameter< String >( "mesh" );
   if( ! this->mesh.load( meshFile ) )
   {
	  std::cerr << "I am not able to load the mesh from the file " << meshFile << "." <<std::endl;
	   return false;
   }

   const IndexType& dofs = this->mesh.getDofs();
   dofVector. setSize( dofs );

   this -> u. bind( & dofVector. getData()[ 0 * dofs ], dofs );
   this -> v. bind( & dofVector. getData()[ 1 * dofs ], dofs );

*/
   differentialOperator.getAnisotropy().setConstant( 1.0 ); //setup( parameters );
   differentialOperator.setSmoothing( 1.0 );
   return true;

}


template< typename Mesh,
          typename HamiltonJacobi,
          typename BoundaryCondition,
          typename RightHandSide>
typename HamiltonJacobiProblem< Mesh,HamiltonJacobi,BoundaryCondition,RightHandSide >::IndexType
         HamiltonJacobiProblem< Mesh,HamiltonJacobi,BoundaryCondition,RightHandSide >::getDofs( const MeshType& mesh ) const
{
   /****
    * Set-up DOFs and supporting grid functions
    */
   return mesh.template getEntitiesCount< typename MeshType::Cell >();
}

template< typename Mesh,
          typename HamiltonJacobi,
          typename BoundaryCondition,
          typename RightHandSide>
void HamiltonJacobiProblem< Mesh,HamiltonJacobi,BoundaryCondition,RightHandSide >::
bindDofs( const MeshType& mesh,
          DofVectorType& dofVector )
{
   //const IndexType dofs = mesh.template getEntitiesCount< typename MeshType::Cell >();
   this->solution.bind( mesh, dofVector );
}

template< typename Mesh,
		  typename HamiltonJacobi,
		  typename BoundaryCondition,
		  typename RightHandSide>
bool
HamiltonJacobiProblem< Mesh,HamiltonJacobi,BoundaryCondition,RightHandSide >::
setInitialCondition( const Config::ParameterContainer& parameters,
                     const MeshType& mesh,
                     DofVectorType& dofs,
                     MeshDependentDataType& meshDependentData  )
{
  this->bindDofs( mesh, dofs );
  const String& initialConditionFile = parameters.getParameter< String >( "initial-condition" );
  if(CommunicatorType::isDistributed())
  {
    std::cout<<"Nodes Distribution: " << uPointer->getMesh().getDistributedMesh()->printProcessDistr() << std::endl;
    if(distributedIOType==Meshes::DistributedMeshes::MpiIO)
      Meshes::DistributedMeshes::DistributedGridIO<MeshFunctionType,Meshes::DistributedMeshes::MpiIO> ::load(initialConditionFile, *uPointer );
    if(distributedIOType==Meshes::DistributedMeshes::LocalCopy)
      Meshes::DistributedMeshes::DistributedGridIO<MeshFunctionType,Meshes::DistributedMeshes::LocalCopy> ::load(initialConditionFile, *uPointer );
    uPointer->template synchronize<CommunicatorType>();
  }
  else
  {
    if( ! this->solution.boundLoad( initialConditionFile ) )
    {
      std::cerr << "I am not able to load the initial condition from the file " << initialConditionFile << "." <<std::endl;
      return false;
    }
  }
   return true;
}

template< typename Mesh,
		  typename HamiltonJacobi,
		  typename BoundaryCondition,
		  typename RightHandSide>
bool 
HamiltonJacobiProblem< Mesh,HamiltonJacobi,BoundaryCondition,RightHandSide >::
makeSnapshot( const RealType& time,
              const IndexType& step,
              const MeshType& mesh,
              DofVectorType& dofs,
              MeshDependentDataType& meshDependentData  )
{
  std::cout <<std::endl << "Writing output at time " << time << " step " << step << "." <<std::endl;

   String fileName;
   FileNameBaseNumberEnding( "u-", step, 5, ".tnl", fileName );
   if( ! this->solution.save( fileName ) )
	   return false;

   return true;
}


template< typename Mesh,
		  typename HamiltonJacobi,
		  typename BoundaryCondition,
		  typename RightHandSide>
void
HamiltonJacobiProblem< Mesh,HamiltonJacobi,BoundaryCondition,RightHandSide >::
getExplicitUpdate(  const RealType& time,
                 const RealType& tau,
                 const MeshType& mesh,
                 DofVectorType& _u,
                 DofVectorType& _fu,
                 MeshDependentDataType& meshDependentData  )
{

	/*
	if(!(this->schemeTest))
		scheme.getExplicitUpdate(time, tau, _u, _fu);
	else if(!(this->tested))
	{
		this->tested = true;
		DofVectorType tmp;
		if(tmp.setLike(_u))
			tmp = _u;
		scheme.getExplicitUpdate(time, tau, tmp, _u);

	}
	*/

	//this->bindDofs( mesh, _u );
   MeshFunctionType u, fu;
   u.bind( mesh, _u );
   fu.bind( mesh, _fu );
   explicitUpdater.template update< typename MeshType::Cell >( time,
	                                                            mesh,
	                                                            this->differentialOperator,
	                                                            this->boundaryCondition,
	                                                            this->rightHandSide,
	                                                            u,
	                                                            fu );
   //fu.save( "fu.tnl" );
   //std::cerr << "Enter." << std::endl;
   //getchar();
}