Refactoring the fast sweeping method. (ace01c90) · Commits · TNL / tnl-dev

examples/hamilton-jacobi/tnlDirectEikonalMethodsBase_impl.h

+35 −33

Original line number	Diff line number	Diff line
		@@ -8,6 +8,7 @@
		#pragma once

		#include <core/tnlTypeInfo.h>
		#include <functions/tnlFunctions.h>

		template< typename Real,
		typename Device,
		@@ -48,38 +49,6 @@ 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 );
		}


		@@ -133,7 +102,40 @@ tnlDirectEikonalMethodsBase< tnlGrid< 2, Real, Device, Index > >::
		updateCell( MeshFunction& u,
		const MeshEntity& cell )
		{
		const auto& neighbourEntities = cell.template getNeighbourEntities< 2 >();
		const MeshType& mesh = cell.getMesh();

		const RealType& h = mesh.getSpaceSteps().x();
		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 = ArgAbsMin( 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 = ArgAbsMin( u[ neighbourEntities.template getEntityIndex< 0, -1 >() ],
		u[ neighbourEntities.template getEntityIndex< 0, 1 >() ] );
		}

		if( fabs( a - b ) >= h )
		tmp = ArgAbsMin( 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() ] = ArgAbsMin( value, tmp );
		std::cerr << ArgAbsMin( value, tmp ) << " ";
		}

examples/hamilton-jacobi/tnlFastSweepingMethod2D_impl.h

+15 −9

Original line number	Diff line number	Diff line
		@@ -62,7 +62,7 @@ solve( const MeshType& mesh,
		aux.setMesh( mesh );
		std::cout << "Initiating the interface cells ..." << std::endl;
		BaseType::initInterface( u, aux );
		aux.save( "aux.tnl" );
		aux.save( "aux-ini.tnl" );

		typename MeshType::Cell cell( mesh );

		@@ -74,8 +74,9 @@ solve( const MeshType& mesh,
		cell.getCoordinates().x() < mesh.getDimensions().x();
		cell.getCoordinates().x()++ )
		{
		std::cerr << "1 -> ";
		cell.refresh();
		this->updateValue( aux, cell );
		this->updateCell( aux, cell );
		}
		}

		@@ -87,8 +88,9 @@ solve( const MeshType& mesh,
		cell.getCoordinates().x() >= 0 ;
		cell.getCoordinates().x()-- )
		{
		std::cerr << "2 -> ";
		cell.refresh();
		this->updateValue( aux, cell );
		this->updateCell( aux, cell );
		}
		}

		@@ -100,8 +102,9 @@ solve( const MeshType& mesh,
		cell.getCoordinates().x() < mesh.getDimensions().x();
		cell.getCoordinates().x()++ )
		{
		std::cerr << "3 -> ";
		cell.refresh();
		this->updateValue( aux, cell );
		this->updateCell( aux, cell );
		}
		}

		@@ -114,9 +117,12 @@ solve( const MeshType& mesh,
		cell.getCoordinates().x() >= 0 ;
		cell.getCoordinates().x()-- )
		{
		std::cerr << "4 -> ";
		cell.refresh();
		this->updateValue( aux, cell );
		this->updateCell( aux, cell );
		}
		}

		aux.save( "aux-fin.tnl" );
		}

src/functions/tnlFunctions.h

+13 −3

Original line number	Diff line number	Diff line
		@@ -5,8 +5,7 @@
		* Created on July 11, 2016, 6:01 PM
		*/

		#ifndef TNLFUNCTIONS_H
		#define TNLFUNCTIONS_H
		#pragma once

		#include <core/tnlCuda.h>

		@@ -37,6 +36,17 @@ Real negativePart( const Real& arg)
		return arg < 0.0 ? arg : 0.0;
		}

		template< typename Real >
		__cuda_callable__
		Real ArgAbsMin( const Real& x, const Real& y )
		{
		return fabs( x ) < fabs( y ) ? x : y;
		}

		#endif /* TNLFUNCTIONS_H */
		template< typename Real >
		__cuda_callable__
		Real ArgAbsMax( const Real& x, const Real& y )
		{
		return fabs( x ) > fabs( y ) ? x : y;
		}