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

Refactoring the fast sweeping method.

parent bc957149

examples/hamilton-jacobi/tnl-direct-eikonal-solver.h

+1 −0
Original line number Diff line number Diff line
@@ -32,6 +32,7 @@ class tnlDirectEikonalSolverConfig 
      static void configSetup( tnlConfigDescription& config )

      {

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

         config.addRequiredEntry< tnlString >( "input-file", "Input file." );

      };

};

examples/hamilton-jacobi/tnlDirectEikonalMethodsBase.h

0 → 100644
+85 −0
Original line number Diff line number Diff line 
/* 

 * File:   tnlDirectEikonalMethodsBase.h

 * Author: oberhuber

 *

 * Created on July 14, 2016, 3:17 PM

 */



#pragma once 



#include <mesh/tnlGrid.h>

#include <functions/tnlMeshFunction.h>



template< typename Mesh >

class tnlDirectEikonalMethodsBase

{   

};



template< typename Real,

          typename Device,

          typename Index >

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

{

   public:

      

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

      typedef Real RealType;

      typedef Device DevcieType;

      typedef Index IndexType;

      

      template< typename MeshFunction >

      void initInterface( const MeshFunction& input,

                          MeshFunction& output );

      

      template< typename MeshFunction, typename MeshEntity >

      void updateCell( MeshFunction& u,

                       const MeshEntity& cell );

      

};





template< typename Real,

          typename Device,

          typename Index >

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

{

   public:

      

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

      typedef Real RealType;

      typedef Device DevcieType;

      typedef Index IndexType;

      typedef tnlMeshFunction< MeshType > MeshFunctionType;

      

      template< typename MeshFunction >

      void initInterface( const MeshFunction& input,

                          MeshFunction& output );

      

      template< typename MeshFunction, typename MeshEntity >

      void updateCell( MeshFunction& u,

                       const MeshEntity& cell );

};



template< typename Real,

          typename Device,

          typename Index >

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

{

   public:

      

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

      typedef Real RealType;

      typedef Device DevcieType;

      typedef Index IndexType;

      typedef tnlMeshFunction< MeshType > MeshFunctionType;

      

      template< typename MeshFunction >

      void initInterface( const MeshFunction& input,

                          MeshFunction& output );

      

      template< typename MeshFunction, typename MeshEntity >

      void updateCell( MeshFunction& u,

                       const MeshEntity& cell );      

};



#include "tnlDirectEikonalMethodsBase_impl.h"

examples/hamilton-jacobi/tnlDirectEikonalMethodsBase_impl.h

0 → 100644
+192 −0
Original line number Diff line number Diff line 
/* 

 * File:   tnlDirectEikonalMethodsBase_impl.h

 * Author: oberhuber

 *

 * Created on July 14, 2016, 3:22 PM

 */



#pragma once



#include <core/tnlTypeInfo.h>



template< typename Real,

          typename Device,

          typename Index >

   template< typename MeshFunction >

void

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

initInterface( const MeshFunction& input,

               MeshFunction& output )

{

   const MeshType& mesh = input.getMesh();

   typedef typename MeshType::Cell Cell;

   Cell cell( mesh );

   for( cell.getCoordinates().x() = 1;

        cell.getCoordinates().x() < mesh.getDimensions().x() - 1;

        cell.getCoordinates().x() ++ )

   {

      cell.refresh();

      const auto& neighbours = cell.getNeighbourEntities();

      //const IndexType& c = cell.getIndex();

      const IndexType e = neighbours.template getEntityIndex<  1 >();

      const IndexType w = neighbours.template getEntityIndex< -1 >();

      const RealType& c = input( cell );

      if( c * input[ e ] <= 0 || c * input[ w ] <= 0 )

         output[ cell.getIndex() ] = c;

      else output[ cell.getIndex() ] =

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

                -tnlTypeInfo< RealType >::getMaxValue();

   }

}



template< typename Real,

          typename Device,

          typename Index >

   template< typename MeshFunction, typename MeshEntity >

void

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

updateCell( MeshFunction& u,

            const MeshEntity& cell )

{

	const auto& neighbourEntities = cell.getNeighbourEntities< 2 >();

   const MeshType& mesh = cell.getMesh();

   

	const RealType& value = u( cell );

	Real a,b, tmp;



	if( cell.getCoordinates().x() == 0 )

		a = u[ neighbourEntities.template getEntityIndex< 1,  0 >() ];

	else if( cell.getCoordinates().x() == mesh.getDimensions().x() - 1 )

		a = u[ neighbourEntities.template getEntityIndex< -1,  0 >() ];

	else

	{

		a = fabsMin( u[ neighbourEntities.template getEntityIndex< -1,  0 >() ],

				       u[ neighbourEntities.template getEntityIndex<  1,  0 >() ] );

	}



	if( cell.getCoordinates().y() == 0 )

		b = u[ neighbourEntities.template getEntityIndex< 0,  1 >()];

	else if( cell.getCoordinates().y() == mesh.getDimensions().y() - 1 )

		b = u[ neighbourEntities.template getEntityIndex< 0,  -1 >() ];

	else

	{

		b = fabsMin( u[ neighbourEntities.template getEntityIndex< 0,  -1 >() ],

				       u[ neighbourEntities.template getEntityIndex< 0,   1 >() ] );

	}



	if( fabs( a - b ) >= h )

		tmp = fabsMin( a, b ) + Sign( value ) * h;

	else

		tmp = 0.5 * (a + b + Sign(value)*sqrt(2.0 * h * h - (a - b) * (a - b) ) );

	

   u[ cell.getIndex() ] = fabsMin( value, tmp );

}





template< typename Real,

          typename Device,

          typename Index >

      template< typename MeshFunction >

void

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

initInterface( const MeshFunction& input,

               MeshFunction& output )

{

   const MeshType& mesh = input.getMesh();

   typedef typename MeshType::Cell Cell;

   Cell cell( mesh );

   for( cell.getCoordinates().y() = 0;

        cell.getCoordinates().y() < mesh.getDimensions().y();

        cell.getCoordinates().y() ++ )

      for( cell.getCoordinates().x() = 0;

           cell.getCoordinates().x() < mesh.getDimensions().x();

           cell.getCoordinates().x() ++ )

      {

         cell.refresh();

         const RealType& c = input( cell );

         if( ! cell.isBoundaryEntity()  )

         {

            auto neighbours = cell.getNeighbourEntities();

            const IndexType e = neighbours.template getEntityIndex<  1,  0 >();

            const IndexType w = neighbours.template getEntityIndex< -1,  0 >();

            const IndexType n = neighbours.template getEntityIndex<  0,  1 >();

            const IndexType s = neighbours.template getEntityIndex<  0, -1 >();            

            if( c * input[ e ] <= 0 || c * input[ w ] <= 0 ||

                c * input[ n ] <= 0 || c * input[ s ] <= 0 )

            {

               output[ cell.getIndex() ] = c;

               continue;

            }

         }

         output[ cell.getIndex() ] =

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

                   -tnlTypeInfo< RealType >::getMaxValue();         

      }

}



template< typename Real,

          typename Device,

          typename Index >

   template< typename MeshFunction, typename MeshEntity >

void

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

updateCell( MeshFunction& u,

            const MeshEntity& cell )

{

   

}





template< typename Real,

          typename Device,

          typename Index >

   template< typename MeshFunction >

void

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

initInterface( const MeshFunction& input,

               MeshFunction& output )

{

   const MeshType& mesh = input.getMesh();

   typedef typename MeshType::Cell Cell;

   Cell cell( mesh );

   for( cell.getCoordinates().z() = 1;

        cell.getCoordinates().z() < mesh.getDimensions().z() - 1;

        cell.getCoordinates().z() ++ )   

      for( cell.getCoordinates().y() = 1;

           cell.getCoordinates().y() < mesh.getDimensions().y() - 1;

           cell.getCoordinates().y() ++ )

         for( cell.getCoordinates().x() = 1;

              cell.getCoordinates().x() < mesh.getDimensions().x() - 1;

              cell.getCoordinates().x() ++ )

         {

            cell.refresh();

            auto neighbours = cell.getNeighbourEntities();

            //const IndexType& c = cell.getIndex();

            const IndexType e = neighbours.template getEntityIndex<  1,  0,  0 >();

            const IndexType w = neighbours.template getEntityIndex< -1,  0,  0 >();

            const IndexType n = neighbours.template getEntityIndex<  0,  1,  0 >();

            const IndexType s = neighbours.template getEntityIndex<  0, -1,  0 >();

            const IndexType t = neighbours.template getEntityIndex<  0,  0,  1 >();

            const IndexType b = neighbours.template getEntityIndex<  0,  0, -1 >();

            const RealType& c = input( cell );

            if( c * input[ e ] <= 0 || c * input[ w ] <= 0 ||

                c * input[ n ] <= 0 || c * input[ s ] <= 0 ||

                c * input[ t ] <= 0 || c * input[ b ] <= 0 )

               output[ cell.getIndex() ] = c;

            else output[ cell.getIndex() ] =

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

                      -tnlTypeInfo< RealType >::getMaxValue();

         }

}



template< typename Real,

          typename Device,

          typename Index >

   template< typename MeshFunction, typename MeshEntity >

void

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

updateCell( MeshFunction& u,

            const MeshEntity& cell )

{

   

}

examples/hamilton-jacobi/tnlDirectEikonalProblem.h

+3 −0
Original line number Diff line number Diff line
@@ -35,6 +35,7 @@ class tnlDirectEikonalProblem 
      typedef Index IndexType;

      typedef tnlMeshFunction< Mesh > MeshFunctionType;

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

      typedef Anisotropy AnisotropyType;



      using typename BaseType::MeshType;

      using typename BaseType::DofVectorType;
@@ -72,6 +73,8 @@ class tnlDirectEikonalProblem 
         

         MeshFunctionType initialData;

         

         AnisotropyType anisotropy;



};



#include "tnlDirectEikonalProblem_impl.h"

examples/hamilton-jacobi/tnlDirectEikonalProblem_impl.h

+14 −8
Original line number Diff line number Diff line
@@ -32,7 +32,7 @@ tnlString 
tnlDirectEikonalProblem< Mesh, Anisotropy, Real, Index >::

getPrologHeader() const

{

   

   return tnlString( "Direct eikonal solver" );

}



template< typename Mesh,
@@ -55,7 +55,7 @@ bool 
tnlDirectEikonalProblem< Mesh, Anisotropy, Real, Index >::

writeEpilog( tnlLogger& logger )

{

   

   return true;

}



template< typename Mesh,
@@ -66,7 +66,7 @@ bool 
tnlDirectEikonalProblem< Mesh, Anisotropy, Real, Index >::

setup( const tnlParameterContainer& parameters )

{

   

   return true;

}



template< typename Mesh,
@@ -77,7 +77,7 @@ Index 
tnlDirectEikonalProblem< Mesh, Anisotropy, Real, Index >::

getDofs( const MeshType& mesh ) const

{

   return mesh.getEntitiesCount();

   return mesh.template getEntitiesCount< typename MeshType::Cell >();

}



template< typename Mesh,
@@ -92,7 +92,12 @@ bindDofs( const MeshType& mesh, 
   this->u.bind( mesh, dofs );

}



template< typename Mesh,

          typename Anisotropy,

          typename Real,

          typename Index >

bool

tnlDirectEikonalProblem< Mesh, Anisotropy, Real, Index >::

setInitialData( const tnlParameterContainer& parameters,

                const MeshType& mesh,

                DofVectorType& dofs,
@@ -102,7 +107,7 @@ setInitialData( const tnlParameterContainer& parameters, 
   this->initialData.setMesh( mesh );

   if( !this->initialData.boundLoad( inputFile ) )

      return false;

   

   return true;

}
@@ -113,7 +118,8 @@ template< typename Mesh, 
bool

tnlDirectEikonalProblem< Mesh, Anisotropy, Real, Index >::

solve( const MeshType& mesh,

       DofVectorType& dosf )

       DofVectorType& dofs )

{

   

   tnlFastSweepingMethod< MeshType, AnisotropyType > fsm;

   fsm.solve( mesh, anisotropy, initialData );

}