Commit 207994b1 authored by Tomáš Oberhuber's avatar Tomáš Oberhuber
Browse files

Refactoring the direct eikonal solver.

parent 488e2376

src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/CMakeLists.txt

+3 −3
Original line number Diff line number Diff line
@@ -7,14 +7,14 @@ set( tnl_hamilton_jacobi_SOURCES 
     hamiltonJacobiProblemConfig.h

     tnl-direct-eikonal-solver.h )

               

ADD_EXECUTABLE(tnl-hamilton-jacobi${debugExt} main.cpp)

target_link_libraries (tnl-hamilton-jacobi${debugExt} tnl${debugExt}-${tnlVersion} )

#ADD_EXECUTABLE(tnl-hamilton-jacobi${debugExt} main.cpp)

#target_link_libraries (tnl-hamilton-jacobi${debugExt} tnl${debugExt}-${tnlVersion} )



ADD_EXECUTABLE(tnl-direct-eikonal-solver${debugExt} tnl-direct-eikonal-solver.cpp )

target_link_libraries (tnl-direct-eikonal-solver${debugExt} tnl${debugExt}-${tnlVersion} )





INSTALL( TARGETS tnl-hamilton-jacobi${debugExt} 

INSTALL( TARGETS #tnl-hamilton-jacobi${debugExt} 

                 tnl-direct-eikonal-solver${debugExt}

         RUNTIME DESTINATION bin

         PERMISSIONS OWNER_READ OWNER_WRITE OWNER_EXECUTE GROUP_READ GROUP_EXECUTE WORLD_READ WORLD_EXECUTE )

src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnl-direct-eikonal-solver.h

+14 −13
Original line number Diff line number Diff line
@@ -13,23 +13,24 @@ 


#pragma once



#include <solvers/tnlSolver.h>

#include <solvers/tnlFastBuildConfigTag.h>

#include <solvers/tnlBuildConfigTags.h>

#include <functions/tnlConstantFunction.h>

#include <functions/tnlMeshFunction.h>

//#include <problems/tnlHeatEquationProblem.h>

#include <mesh/tnlGrid.h>

#include <TNL/Solvers/Solver.h>

#include <TNL/Solvers/FastBuildConfigTag.h>

#include <TNL/Solvers/BuildConfigTags.h>

#include <TNL/Functions/Analytic/Constant.h>

#include <TNL/Functions/MeshFunction.h>

#include <TNL/Meshes/Grid.h>

#include "tnlDirectEikonalProblem.h"



using namespace TNL;



//typedef tnlDefaultBuildMeshConfig BuildConfig;

typedef tnlFastBuildConfig BuildConfig;

typedef Solvers::FastBuildConfig BuildConfig;



template< typename MeshConfig >

class tnlDirectEikonalSolverConfig

{

   public:

      static void configSetup( tnlConfigDescription& config )

      static void configSetup( Config::ConfigDescription& config )

      {

         config.addDelimiter( "Direct eikonal equation solver settings:" );

         config.addRequiredEntry< String >( "input-file", "Input file." );
@@ -50,12 +51,12 @@ class tnlDirectEikonalSolverSetter 
   typedef Device DeviceType;

   typedef Index IndexType;



   typedef Containers::StaticVector< MeshType::meshDimensions, Real > Point;

   typedef Containers::StaticVector< MeshType::getMeshDimension(), Real > Point;



   static bool run( const Config::ParameterContainer& parameters )

   {

      enum { Dimensions = MeshType::meshDimensions };

      typedef tnlConstantFunction< Dimensions, Real > Anisotropy;

      static const int Dimension = MeshType::getMeshDimension();

      typedef Functions::Analytic::Constant< Dimension, Real > Anisotropy;

      typedef tnlDirectEikonalProblem< MeshType, Anisotropy > Problem;

      SolverStarter solverStarter;

      return solverStarter.template run< Problem >( parameters );
@@ -64,7 +65,7 @@ class tnlDirectEikonalSolverSetter 


int main( int argc, char* argv[] )

{

   if( ! tnlSolver< tnlDirectEikonalSolverSetter, tnlDirectEikonalSolverConfig, BuildConfig >::run( argc, argv ) )

   if( ! Solvers::Solver< tnlDirectEikonalSolverSetter, tnlDirectEikonalSolverConfig, BuildConfig >::run( argc, argv ) )

      return EXIT_FAILURE;

   return EXIT_SUCCESS;

}

src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlDirectEikonalMethodsBase.h

+16 −14
Original line number Diff line number Diff line
@@ -7,8 +7,10 @@ 


#pragma once 



#include <mesh/tnlGrid.h>

#include <functions/tnlMeshFunction.h>

#include <TNL/Meshes/Grid.h>

#include <TNL/Functions/MeshFunction.h>



using namespace TNL;



template< typename Mesh >

class tnlDirectEikonalMethodsBase
@@ -18,16 +20,16 @@ class tnlDirectEikonalMethodsBase 
template< typename Real,

          typename Device,

          typename Index >

class tnlDirectEikonalMethodsBase< tnlGrid< 1, Real, Device, Index > >

class tnlDirectEikonalMethodsBase< Meshes::Grid< 1, Real, Device, Index > >

{

   public:

      

      typedef tnlGrid< 1, Real, Device, Index > MeshType;

      typedef Meshes::Grid< 1, Real, Device, Index > MeshType;

      typedef Real RealType;

      typedef Device DevcieType;

      typedef Index IndexType;

      typedef tnlMeshFunction< MeshType > MeshFunctionType;

      typedef tnlMeshFunction< MeshType, 1, bool > InterfaceMapType;

      typedef Functions::MeshFunction< MeshType > MeshFunctionType;

      typedef Functions::MeshFunction< MeshType, 1, bool > InterfaceMapType;

      

      void initInterface( const MeshFunctionType& input,

                          MeshFunctionType& output,
@@ -43,16 +45,16 @@ class tnlDirectEikonalMethodsBase< tnlGrid< 1, Real, Device, Index > > 
template< typename Real,

          typename Device,

          typename Index >

class tnlDirectEikonalMethodsBase< tnlGrid< 2, Real, Device, Index > >

class tnlDirectEikonalMethodsBase< Meshes::Grid< 2, Real, Device, Index > >

{

   public:

      

      typedef tnlGrid< 2, Real, Device, Index > MeshType;

      typedef Meshes::Grid< 2, Real, Device, Index > MeshType;

      typedef Real RealType;

      typedef Device DevcieType;

      typedef Index IndexType;

      typedef tnlMeshFunction< MeshType > MeshFunctionType;

      typedef tnlMeshFunction< MeshType, 2, bool > InterfaceMapType;

      typedef Functions::MeshFunction< MeshType > MeshFunctionType;

      typedef Functions::MeshFunction< MeshType, 2, bool > InterfaceMapType;



      void initInterface( const MeshFunctionType& input,

                          MeshFunctionType& output,
@@ -66,16 +68,16 @@ class tnlDirectEikonalMethodsBase< tnlGrid< 2, Real, Device, Index > > 
template< typename Real,

          typename Device,

          typename Index >

class tnlDirectEikonalMethodsBase< tnlGrid< 3, Real, Device, Index > >

class tnlDirectEikonalMethodsBase< Meshes::Grid< 3, Real, Device, Index > >

{

   public:

      

      typedef tnlGrid< 3, Real, Device, Index > MeshType;

      typedef Meshes::Grid< 3, Real, Device, Index > MeshType;

      typedef Real RealType;

      typedef Device DevcieType;

      typedef Index IndexType;

      typedef tnlMeshFunction< MeshType > MeshFunctionType;

      typedef tnlMeshFunction< MeshType, 3, bool > InterfaceMapType;

      typedef Functions::MeshFunction< MeshType > MeshFunctionType;

      typedef Functions::MeshFunction< MeshType, 3, bool > InterfaceMapType;



      void initInterface( const MeshFunctionType& input,

                          MeshFunctionType& output,

src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlDirectEikonalMethodsBase_impl.h

+17 −18
Original line number Diff line number Diff line
@@ -7,14 +7,13 @@ 


#pragma once



#include <core/tnlTypeInfo.h>

#include <functions/tnlFunctions.h>

#include <TNL/TypeInfo.h>



template< typename Real,

          typename Device,

          typename Index >

void

tnlDirectEikonalMethodsBase< tnlGrid< 1, Real, Device, Index > >::

tnlDirectEikonalMethodsBase< Meshes::Grid< 1, Real, Device, Index > >::

initInterface( const MeshFunctionType& input,

               MeshFunctionType& output,

               InterfaceMapType& interfaceMap  )
@@ -43,8 +42,8 @@ initInterface( const MeshFunctionType& input, 
         }

      }

      output[ cell.getIndex() ] =

      c > 0 ? tnlTypeInfo< RealType >::getMaxValue() :

             -tnlTypeInfo< RealType >::getMaxValue();

      c > 0 ? TypeInfo< RealType >::getMaxValue() :

             -TypeInfo< RealType >::getMaxValue();

      interfaceMap[ cell.getIndex() ] = false;

   }

}
@@ -54,7 +53,7 @@ template< typename Real, 
          typename Index >

   template< typename MeshEntity >

void

tnlDirectEikonalMethodsBase< tnlGrid< 1, Real, Device, Index > >::

tnlDirectEikonalMethodsBase< Meshes::Grid< 1, Real, Device, Index > >::

updateCell( MeshFunctionType& u,

            const MeshEntity& cell )

{
@@ -65,7 +64,7 @@ template< typename Real, 
          typename Device,

          typename Index >

void

tnlDirectEikonalMethodsBase< tnlGrid< 2, Real, Device, Index > >::

tnlDirectEikonalMethodsBase< Meshes::Grid< 2, Real, Device, Index > >::

initInterface( const MeshFunctionType& input,

               MeshFunctionType& output,

               InterfaceMapType& interfaceMap  )
@@ -98,8 +97,8 @@ initInterface( const MeshFunctionType& input, 
            }

         }

         output[ cell.getIndex() ] =

            c > 0 ? tnlTypeInfo< RealType >::getMaxValue() :

                   -tnlTypeInfo< RealType >::getMaxValue();  

            c > 0 ? TypeInfo< RealType >::getMaxValue() :

                   -TypeInfo< RealType >::getMaxValue();  

         interfaceMap[ cell.getIndex() ] = false;

      }

}
@@ -109,7 +108,7 @@ template< typename Real, 
          typename Index >

   template< typename MeshEntity >

void

tnlDirectEikonalMethodsBase< tnlGrid< 2, Real, Device, Index > >::

tnlDirectEikonalMethodsBase< Meshes::Grid< 2, Real, Device, Index > >::

updateCell( MeshFunctionType& u,

            const MeshEntity& cell )

{
@@ -140,11 +139,11 @@ updateCell( MeshFunctionType& u, 
                     u[ neighborEntities.template getEntityIndex< 0,   1 >() ] );

   }



   if( fabs( a ) == tnlTypeInfo< Real >::getMaxValue() && 

       fabs( b ) == tnlTypeInfo< Real >::getMaxValue() )

   if( fabs( a ) == TypeInfo< Real >::getMaxValue() && 

       fabs( b ) == TypeInfo< Real >::getMaxValue() )

      return;

   if( fabs( a ) == tnlTypeInfo< Real >::getMaxValue() ||

       fabs( b ) == tnlTypeInfo< Real >::getMaxValue() ||

   if( fabs( a ) == TypeInfo< Real >::getMaxValue() ||

       fabs( b ) == TypeInfo< Real >::getMaxValue() ||

       fabs( a - b ) >= h )

   {

      tmp = ArgAbsMin( a, b ) + sign( value ) * h;
@@ -173,7 +172,7 @@ template< typename Real, 
          typename Device,

          typename Index >

void

tnlDirectEikonalMethodsBase< tnlGrid< 3, Real, Device, Index > >::

tnlDirectEikonalMethodsBase< Meshes::Grid< 3, Real, Device, Index > >::

initInterface( const MeshFunctionType& input,

               MeshFunctionType& output,

               InterfaceMapType& interfaceMap  )
@@ -214,8 +213,8 @@ initInterface( const MeshFunctionType& input, 
               }

            }

            output[ cell.getIndex() ] =

               c > 0 ? tnlTypeInfo< RealType >::getMaxValue() :

                      -tnlTypeInfo< RealType >::getMaxValue();

               c > 0 ? TypeInfo< RealType >::getMaxValue() :

                      -TypeInfo< RealType >::getMaxValue();

            interfaceMap[ cell.getIndex() ] = false;

         }

}
@@ -225,7 +224,7 @@ template< typename Real, 
          typename Index >

   template< typename MeshEntity >

void

tnlDirectEikonalMethodsBase< tnlGrid< 3, Real, Device, Index > >::

tnlDirectEikonalMethodsBase< Meshes::Grid< 3, Real, Device, Index > >::

updateCell( MeshFunctionType& u,

            const MeshEntity& cell )

{

src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlDirectEikonalProblem.h

+27 −22
Original line number Diff line number Diff line
@@ -13,8 +13,9 @@ 


#pragma once



#include <problems/tnlPDEProblem.h>

#include <functions/tnlMeshFunction.h>

#include <TNL/Problems/PDEProblem.h>

#include <TNL/Functions/MeshFunction.h>

#include <TNL/SharedPointer.h>

#include "tnlFastSweepingMethod.h"



template< typename Mesh,
@@ -22,8 +23,7 @@ template< typename Mesh, 
          typename Real = typename Mesh::RealType,

          typename Index = typename Mesh::IndexType >

class tnlDirectEikonalProblem

   : public tnlPDEProblem< Mesh,

                           TimeIndependentProblem,

   : public Problems::PDEProblem< Mesh,

                                  Real,

                                  typename Mesh::DeviceType,

                                  Index  >
@@ -33,38 +33,43 @@ class tnlDirectEikonalProblem 
      typedef Real RealType;

      typedef typename Mesh::DeviceType DeviceType;

      typedef Index IndexType;

      typedef tnlMeshFunction< Mesh > MeshFunctionType;

      typedef tnlPDEProblem< Mesh, TimeIndependentProblem, RealType, DeviceType, IndexType > BaseType;

      typedef Anisotropy AnisotropyType;

      typedef Functions::MeshFunction< Mesh > MeshFunctionType;

      typedef Problems::PDEProblem< Mesh, RealType, DeviceType, IndexType > BaseType;

      using AnisotropyType = Anisotropy;



      using typename BaseType::MeshType;

      using typename BaseType::DofVectorType;

      using typename BaseType::MeshDependentDataType;

      using MeshPointer = SharedPointer< MeshType >;

      using DofVectorPointer = SharedPointer< DofVectorType >;

      using MeshDependentDataPointer = SharedPointer< MeshDependentDataType >;



      static String getTypeStatic();



      String getPrologHeader() const;



      void writeProlog( tnlLogger& logger,

      void writeProlog( Logger& logger,

                        const Config::ParameterContainer& parameters ) const;

      

      bool writeEpilog( tnlLogger& logger );

      bool writeEpilog( Logger& logger );





      bool setup( const Config::ParameterContainer& parameters );

      bool setup( const MeshPointer& mesh,

                  const Config::ParameterContainer& parameters,

                  const String& prefix );



      IndexType getDofs( const MeshType& mesh ) const;

      IndexType getDofs( const MeshPointer& mesh ) const;



      void bindDofs( const MeshType& mesh,

                     const DofVectorType& dofs );

      void bindDofs( const MeshPointer& mesh,

                     const DofVectorPointer& dofs );

      

      bool setInitialData( const Config::ParameterContainer& parameters,

                           const MeshType& mesh,

                           DofVectorType& dofs,

                           MeshDependentDataType& meshdependentData );

      bool setInitialCondition( const Config::ParameterContainer& parameters,

                                const MeshPointer& mesh,

                                DofVectorPointer& dofs,

                                MeshDependentDataPointer& meshdependentData );



      bool solve( const MeshType& mesh,

                  DofVectorType& dosf );

      bool solve( const MeshPointer& mesh,

                  DofVectorPointer& dosf );





      protected: