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

examples/hamilton-jacobi/tnlDirectEikonalMethodsBase.h

+21 −17

Original line number	Diff line number	Diff line
		@@ -26,13 +26,15 @@ class tnlDirectEikonalMethodsBase< tnlGrid< 1, Real, Device, Index > >
		typedef Real RealType;
		typedef Device DevcieType;
		typedef Index IndexType;
		typedef tnlMeshFunction< MeshType > MeshFunctionType;
		typedef tnlMeshFunction< MeshType, 1, bool > InterfaceMapType;

		template< typename MeshFunction >
		void initInterface( const MeshFunction& input,
		MeshFunction& output );
		void initInterface( const MeshFunctionType& input,
		MeshFunctionType& output,
		InterfaceMapType& interfaceMap );

		template< typename MeshFunction, typename MeshEntity >
		void updateCell( MeshFunction& u,
		template< typename MeshEntity >
		void updateCell( MeshFunctionType& u,
		const MeshEntity& cell );

		};
		@@ -50,13 +52,14 @@ class tnlDirectEikonalMethodsBase< tnlGrid< 2, Real, Device, Index > >
		typedef Device DevcieType;
		typedef Index IndexType;
		typedef tnlMeshFunction< MeshType > MeshFunctionType;
		typedef tnlMeshFunction< MeshType, 2, bool > InterfaceMapType;

		template< typename MeshFunction >
		void initInterface( const MeshFunction& input,
		MeshFunction& output );
		void initInterface( const MeshFunctionType& input,
		MeshFunctionType& output,
		InterfaceMapType& interfaceMap );

		template< typename MeshFunction, typename MeshEntity >
		void updateCell( MeshFunction& u,
		template< typename MeshEntity >
		void updateCell( MeshFunctionType& u,
		const MeshEntity& cell );
		};

		@@ -72,13 +75,14 @@ class tnlDirectEikonalMethodsBase< tnlGrid< 3, Real, Device, Index > >
		typedef Device DevcieType;
		typedef Index IndexType;
		typedef tnlMeshFunction< MeshType > MeshFunctionType;
		typedef tnlMeshFunction< MeshType, 3, bool > InterfaceMapType;

		template< typename MeshFunction >
		void initInterface( const MeshFunction& input,
		MeshFunction& output );
		void initInterface( const MeshFunctionType& input,
		MeshFunctionType& output,
		InterfaceMapType& interfaceMap );

		template< typename MeshFunction, typename MeshEntity >
		void updateCell( MeshFunction& u,
		template< typename MeshEntity >
		void updateCell( MeshFunctionType& u,
		const MeshEntity& cell );
		};

examples/hamilton-jacobi/tnlDirectEikonalMethodsBase_impl.h

+63 −40

Original line number	Diff line number	Diff line
		@@ -13,11 +13,11 @@
		template< typename Real,
		typename Device,
		typename Index >
		template< typename MeshFunction >
		void
		tnlDirectEikonalMethodsBase< tnlGrid< 1, Real, Device, Index > >::
		initInterface( const MeshFunction& input,
		MeshFunction& output )
		initInterface( const MeshFunctionType& input,
		MeshFunctionType& output,
		InterfaceMapType& interfaceMap )
		{
		const MeshType& mesh = input.getMesh();
		typedef typename MeshType::Cell Cell;
		@@ -27,26 +27,35 @@ initInterface( const MeshFunction& input,
		cell.getCoordinates().x() ++ )
		{
		cell.refresh();
		const RealType& c = input( cell );
		if( ! cell.isBoundaryEntity() )
		{
		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() ] =
		interfaceMap[ cell.getIndex() ] = true;
		continue;
		}
		}
		output[ cell.getIndex() ] =
		c > 0 ? tnlTypeInfo< RealType >::getMaxValue() :
		-tnlTypeInfo< RealType >::getMaxValue();
		interfaceMap[ cell.getIndex() ] = false;
		}
		}

		template< typename Real,
		typename Device,
		typename Index >
		template< typename MeshFunction, typename MeshEntity >
		template< typename MeshEntity >
		void
		tnlDirectEikonalMethodsBase< tnlGrid< 1, Real, Device, Index > >::
		updateCell( MeshFunction& u,
		updateCell( MeshFunctionType& u,
		const MeshEntity& cell )
		{
		}
		@@ -55,11 +64,11 @@ updateCell( MeshFunction& u,
		template< typename Real,
		typename Device,
		typename Index >
		template< typename MeshFunction >
		void
		tnlDirectEikonalMethodsBase< tnlGrid< 2, Real, Device, Index > >::
		initInterface( const MeshFunction& input,
		MeshFunction& output )
		initInterface( const MeshFunctionType& input,
		MeshFunctionType& output,
		InterfaceMapType& interfaceMap )
		{
		const MeshType& mesh = input.getMesh();
		typedef typename MeshType::Cell Cell;
		@@ -84,22 +93,24 @@ initInterface( const MeshFunction& input,
		c * input[ n ] <= 0 \|\| c * input[ s ] <= 0 )
		{
		output[ cell.getIndex() ] = c;
		interfaceMap[ cell.getIndex() ] = true;
		continue;
		}
		}
		output[ cell.getIndex() ] =
		c > 0 ? tnlTypeInfo< RealType >::getMaxValue() :
		-tnlTypeInfo< RealType >::getMaxValue();
		interfaceMap[ cell.getIndex() ] = false;
		}
		}

		template< typename Real,
		typename Device,
		typename Index >
		template< typename MeshFunction, typename MeshEntity >
		template< typename MeshEntity >
		void
		tnlDirectEikonalMethodsBase< tnlGrid< 2, Real, Device, Index > >::
		updateCell( MeshFunction& u,
		updateCell( MeshFunctionType& u,
		const MeshEntity& cell )
		{
		const auto& neighbourEntities = cell.template getNeighbourEntities< 2 >();
		@@ -129,24 +140,27 @@ updateCell( MeshFunction& u,
		u[ neighbourEntities.template getEntityIndex< 0, 1 >() ] );
		}

		if( fabs( a ) == tnlTypeInfo< Real >::getMaxValue() \|\|
		fabs( a ) == tnlTypeInfo< Real >::getMaxValue() )
		return;
		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 ) << " ";
		//std::cerr << ArgAbsMin( value, tmp ) << " ";
		}


		template< typename Real,
		typename Device,
		typename Index >
		template< typename MeshFunction >
		void
		tnlDirectEikonalMethodsBase< tnlGrid< 3, Real, Device, Index > >::
		initInterface( const MeshFunction& input,
		MeshFunction& output )
		initInterface( const MeshFunctionType& input,
		MeshFunctionType& output,
		InterfaceMapType& interfaceMap )
		{
		const MeshType& mesh = input.getMesh();
		typedef typename MeshType::Cell Cell;
		@@ -162,6 +176,9 @@ initInterface( const MeshFunction& input,
		cell.getCoordinates().x() ++ )
		{
		cell.refresh();
		const RealType& c = input( cell );
		if( ! cell.isBoundaryEntity() )
		{
		auto neighbours = cell.getNeighbourEntities();
		//const IndexType& c = cell.getIndex();
		const IndexType e = neighbours.template getEntityIndex< 1, 0, 0 >();
		@@ -170,24 +187,30 @@ initInterface( const MeshFunction& input,
		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() ] =
		interfaceMap[ cell.getIndex() ] = true;
		continue;
		}
		}
		output[ cell.getIndex() ] =
		c > 0 ? tnlTypeInfo< RealType >::getMaxValue() :
		-tnlTypeInfo< RealType >::getMaxValue();
		interfaceMap[ cell.getIndex() ] = false;
		}
		}

		template< typename Real,
		typename Device,
		typename Index >
		template< typename MeshFunction, typename MeshEntity >
		template< typename MeshEntity >
		void
		tnlDirectEikonalMethodsBase< tnlGrid< 3, Real, Device, Index > >::
		updateCell( MeshFunction& u,
		updateCell( MeshFunctionType& u,
		const MeshEntity& cell )
		{

examples/hamilton-jacobi/tnlFastSweepingMethod.h

+12 −6

Original line number	Diff line number	Diff line
		@@ -39,9 +39,11 @@ class tnlFastSweepingMethod< tnlGrid< 1, Real, Device, Index >, Anisotropy >
		typedef tnlHost DeviceType;
		typedef Index IndexType;
		typedef Anisotropy AnisotropyType;
		typedef tnlMeshFunction< MeshType > MeshFunctionType;
		typedef tnlDirectEikonalMethodsBase< tnlGrid< 1, Real, Device, Index > > BaseType;

		using typename BaseType::InterfaceMapType;
		using typename BaseType::MeshFunctionType;

		tnlFastSweepingMethod();

		const IndexType& getMaxIterations() const;
		@@ -74,9 +76,11 @@ class tnlFastSweepingMethod< tnlGrid< 2, Real, Device, Index >, Anisotropy >
		typedef tnlHost DeviceType;
		typedef Index IndexType;
		typedef Anisotropy AnisotropyType;
		typedef tnlMeshFunction< MeshType > MeshFunctionType;
		typedef tnlDirectEikonalMethodsBase< tnlGrid< 2, Real, Device, Index > > BaseType;

		using typename BaseType::InterfaceMapType;
		using typename BaseType::MeshFunctionType;

		tnlFastSweepingMethod();

		const IndexType& getMaxIterations() const;
		@@ -109,9 +113,11 @@ class tnlFastSweepingMethod< tnlGrid< 3, Real, Device, Index >, Anisotropy >
		typedef tnlHost DeviceType;
		typedef Index IndexType;
		typedef Anisotropy AnisotropyType;
		typedef tnlMeshFunction< MeshType > MeshFunctionType;
		typedef tnlDirectEikonalMethodsBase< tnlGrid< 3, Real, Device, Index > > BaseType;

		using typename BaseType::InterfaceMapType;
		using typename BaseType::MeshFunctionType;

		tnlFastSweepingMethod();

		const IndexType& getMaxIterations() const;

examples/hamilton-jacobi/tnlFastSweepingMethod1D_impl.h

+5 −2

Original line number	Diff line number	Diff line
		@@ -59,9 +59,12 @@ solve( const MeshType& mesh,
		MeshFunctionType& u )
		{
		MeshFunctionType aux;
		InterfaceMapType interfaceMap;
		aux.setMesh( mesh );
		interfaceMap.setMesh( mesh );
		std::cout << "Initiating the interface cells ..." << std::endl;
		BaseType::initInterface( u, aux );
		aux.save( "aux.tnl" );
		BaseType::initInterface( u, aux, interfaceMap );
		aux.save( "aux-ini.tnl" );

		}

examples/hamilton-jacobi/tnlFastSweepingMethod2D_impl.h

+20 −12

Original line number	Diff line number	Diff line
		@@ -59,9 +59,11 @@ solve( const MeshType& mesh,
		MeshFunctionType& u )
		{
		MeshFunctionType aux;
		InterfaceMapType interfaceMap;
		aux.setMesh( mesh );
		interfaceMap.setMesh( mesh );
		std::cout << "Initiating the interface cells ..." << std::endl;
		BaseType::initInterface( u, aux );
		BaseType::initInterface( u, aux, interfaceMap );
		aux.save( "aux-ini.tnl" );

		typename MeshType::Cell cell( mesh );
		@@ -74,11 +76,13 @@ solve( const MeshType& mesh,
		cell.getCoordinates().x() < mesh.getDimensions().x();
		cell.getCoordinates().x()++ )
		{
		std::cerr << "1 -> ";
		//std::cerr << "1 -> ";
		cell.refresh();
		if( ! interfaceMap( cell ) )
		this->updateCell( aux, cell );
		}
		}
		aux.save( "aux-1.tnl" );

		for( cell.getCoordinates().y() = 0;
		cell.getCoordinates().y() < mesh.getDimensions().y();
		@@ -88,11 +92,13 @@ solve( const MeshType& mesh,
		cell.getCoordinates().x() >= 0 ;
		cell.getCoordinates().x()-- )
		{
		std::cerr << "2 -> ";
		//std::cerr << "2 -> ";
		cell.refresh();
		if( ! interfaceMap( cell ) )
		this->updateCell( aux, cell );
		}
		}
		aux.save( "aux-2.tnl" );

		for( cell.getCoordinates().y() = mesh.getDimensions().y() - 1;
		cell.getCoordinates().y() >= 0 ;
		@@ -102,11 +108,13 @@ solve( const MeshType& mesh,
		cell.getCoordinates().x() < mesh.getDimensions().x();
		cell.getCoordinates().x()++ )
		{
		std::cerr << "3 -> ";
		//std::cerr << "3 -> ";
		cell.refresh();
		if( ! interfaceMap( cell ) )
		this->updateCell( aux, cell );
		}
		}
		aux.save( "aux-3.tnl" );


		for( cell.getCoordinates().y() = mesh.getDimensions().y() - 1;
		@@ -117,12 +125,12 @@ solve( const MeshType& mesh,
		cell.getCoordinates().x() >= 0 ;
		cell.getCoordinates().x()-- )
		{
		std::cerr << "4 -> ";
		//std::cerr << "4 -> ";
		cell.refresh();
		if( ! interfaceMap( cell ) )
		this->updateCell( aux, cell );
		}
		}

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