Commit e2582089 authored by Tomas Sobotik's avatar Tomas Sobotik
Browse files

Minor tweaks, CPU version for fast-sweeping-map.

parent a3ec3adb
Loading
Loading
Loading
Loading
+2 −2
Original line number Diff line number Diff line
@@ -17,9 +17,9 @@

#include "MainBuildConfig.h"
	//for HOST versions:
//#include "tnlFastSweepingMap.h"
#include "tnlFastSweepingMap.h"
	//for DEVICE versions:
#include "tnlFastSweepingMap_CUDA.h"
//#include "tnlFastSweepingMap_CUDA.h"
#include "fastSweepingMapConfig.h"
#include <solvers/tnlBuildConfigTags.h>

+3 −1
Original line number Diff line number Diff line
@@ -99,8 +99,10 @@ protected:

	bool exactInput;

	int something_changed;

	tnlMeshFunction<MeshType> dofVector, dofVector2;
	DofVectorType data;
	DofVectorType data,map;

	RealType h;

+5 −1
Original line number Diff line number Diff line
@@ -247,7 +247,7 @@ void tnlFastSweepingMap< tnlGrid< 2,MeshReal, Device, MeshIndex >, Real, Index >
	//	cudaDofVector2[Entity.getIndex()]  = fabsMin(value, tmp);
		atomicFabsMin(&(cudaDofVector2[Entity.getIndex()]), tmp);

		if(abs(value)-abs(tmp) > 0.1*h)
		if(abs(value)-abs(tmp) > 0.0)
			atomicMax(something_changed,1);
	}
	else
@@ -280,7 +280,11 @@ bool tnlFastSweepingMap< tnlGrid< 2,MeshReal, Device, MeshIndex >, Real, Index >
	cudaDofVector2[gid] = INT_MAX*Sign(cudaDofVector[gid]);

	if(abs(cudaDofVector[gid]) < 1.01*h)
	{
		cudaDofVector2[gid] = cudaDofVector[gid];
		if(map_cuda[gid] != 0.0)
			cudaDofVector2[gid] /=map_cuda[gid];
	}



+166 −267
Original line number Diff line number Diff line
@@ -16,6 +16,10 @@
#ifndef TNLFASTSWEEPING2D_IMPL_H_
#define TNLFASTSWEEPING2D_IMPL_H_


#define MAP_SOLVER_MAX_VALUE 3


#include "tnlFastSweepingMap.h"

template< typename MeshReal,
@@ -70,6 +74,10 @@ bool tnlFastSweepingMap< tnlGrid< 2,MeshReal, Device, MeshIndex >, Real, Index >
	}
	dofVector2.load(initialCondition);

	const tnlString& mapFile = parameters.getParameter <tnlString>("map");
	if(! this->map.load( mapFile ))
		cout << "Failed to load map file : " << mapFile << endl;

	h = Mesh.template getSpaceStepsProducts< 1, 0 >();
	Entity.refresh();

@@ -81,6 +89,8 @@ bool tnlFastSweepingMap< tnlGrid< 2,MeshReal, Device, MeshIndex >, Real, Index >
		exactInput=true;

	cout << "a" << endl;

	something_changed = 1;
	return initGrid();
}

@@ -97,226 +107,101 @@ bool tnlFastSweepingMap< tnlGrid< 2,MeshReal, Device, MeshIndex >, Real, Index >
	for(int i=0; i< Mesh.getDimensions().x()*Mesh.getDimensions().x();i++)
	{
		dofVector2[i]=INT_MAX*Sign(dofVector[i]);
	}

	for(int i = 0 ; i < Mesh.getDimensions().x()-1; i++)
	{
		for(int j = 0 ; j < Mesh.getDimensions().x()-1; j++)
			{
			this->Entity.setCoordinates(CoordinatesType(i,j));
			this->Entity.refresh();
			neighbourEntities.refresh(Mesh,Entity.getIndex());

				if(dofVector[this->Entity.getIndex()] > 0)
				{
					if(dofVector[neighbourEntities.template getEntityIndex< 1,  0 >()] > 0)
					{
						if(dofVector[neighbourEntities.template getEntityIndex< 0,  1 >()] > 0)
						{
							if(dofVector[neighbourEntities.template getEntityIndex< 1,  1 >()] > 0)
								setupSquare1111(i,j);
							else
								setupSquare1110(i,j);
						}
						else
						{
							if(dofVector[neighbourEntities.template getEntityIndex< 1,  1 >()] > 0)
								setupSquare1101(i,j);
							else
								setupSquare1100(i,j);
						}
					}
					else
					{
						if(dofVector[neighbourEntities.template getEntityIndex< 0,  1 >()] > 0)
						{
							if(dofVector[neighbourEntities.template getEntityIndex< 1,  1 >()] > 0)
								setupSquare1011(i,j);
							else
								setupSquare1010(i,j);
						}
						else
		if(abs(dofVector[i]) < 1.01*h)
		{
							if(dofVector[neighbourEntities.template getEntityIndex< 1,  1 >()] > 0)
								setupSquare1001(i,j);
							else
								setupSquare1000(i,j);
			dofVector2[i] = dofVector[i];
			if(map[i] != 0.0)
				dofVector2[i] /= map[i];
		}
	}
				}
				else
				{
					if(dofVector[neighbourEntities.template getEntityIndex< 1,  0 >()] > 0)
					{
						if(dofVector[neighbourEntities.template getEntityIndex< 0,  1 >()] > 0)
						{
							if(dofVector[neighbourEntities.template getEntityIndex< 1,  1 >()] > 0)
								setupSquare0111(i,j);
							else
								setupSquare0110(i,j);
						}
						else
						{
							if(dofVector[neighbourEntities.template getEntityIndex< 1,  1 >()] > 0)
								setupSquare0101(i,j);
							else
								setupSquare0100(i,j);
						}
					}
					else
					{
						if(dofVector[neighbourEntities.template getEntityIndex< 0,  1 >()] > 0)
						{
							if(dofVector[neighbourEntities.template getEntityIndex< 1,  1 >()] > 0)
								setupSquare0011(i,j);
							else
								setupSquare0010(i,j);
						}
						else
						{
							if(dofVector[neighbourEntities.template getEntityIndex< 1,  1 >()] > 0)
								setupSquare0001(i,j);
							else
								setupSquare0000(i,j);
						}
					}
				}

			}
	}
	cout << "a" << endl;

//	Real tmp = 0.0;
//	Real ax=0.5/sqrt(2.0);
//	for(int i = 0 ; i < Mesh.getDimensions().x()-1; i++)
//	{
//		for(int j = 0 ; j < Mesh.getDimensions().x()-1; j++)
//			{
//			this->Entity.setCoordinates(CoordinatesType(i,j));
//			this->Entity.refresh();
//			neighbourEntities.refresh(Mesh,Entity.getIndex());
//
//	if(!exactInput)
//				if(dofVector[this->Entity.getIndex()] > 0)
//				{
//					if(dofVector[neighbourEntities.template getEntityIndex< 1,  0 >()] > 0)
//					{
//		for(Index i = 0; i < Mesh.getDimensions().x()*Mesh.getDimensions().y(); i++)
//				dofVector[i]=0.5*h*Sign(dofVector[i]);
//						if(dofVector[neighbourEntities.template getEntityIndex< 0,  1 >()] > 0)
//						{
//							if(dofVector[neighbourEntities.template getEntityIndex< 1,  1 >()] > 0)
//								setupSquare1111(i,j);
//							else
//								setupSquare1110(i,j);
//						}
//
//
//	for(Index i = 1; i < Mesh.getDimensions().x()-1; i++)
//						else
//						{
//							if(dofVector[neighbourEntities.template getEntityIndex< 1,  1 >()] > 0)
//								setupSquare1101(i,j);
//							else
//								setupSquare1100(i,j);
//						}
//					}
//					else
//					{
//		for(Index j = 1; j < Mesh.getDimensions().y()-1; j++)
//						if(dofVector[neighbourEntities.template getEntityIndex< 0,  1 >()] > 0)
//						{
//			 tmp = Sign(dofVector[Mesh.getCellIndex(CoordinatesType(i,j))]);
//
//			if(tmp == 0.0)
//			{}
//			else if(dofVector[Mesh.getCellIndex(CoordinatesType(i+1,j))]*tmp < 0.0 ||
//					dofVector[Mesh.getCellIndex(CoordinatesType(i-1,j))]*tmp < 0.0 ||
//					dofVector[Mesh.getCellIndex(CoordinatesType(i,j+1))]*tmp < 0.0 ||
//					dofVector[Mesh.getCellIndex(CoordinatesType(i,j-1))]*tmp < 0.0 )
//			{}
//							if(dofVector[neighbourEntities.template getEntityIndex< 1,  1 >()] > 0)
//								setupSquare1011(i,j);
//							else
//				dofVector[Mesh.getCellIndex(CoordinatesType(i,j))] = tmp*INT_MAX;
//								setupSquare1010(i,j);
//						}
//						else
//						{
//							if(dofVector[neighbourEntities.template getEntityIndex< 1,  1 >()] > 0)
//								setupSquare1001(i,j);
//							else
//								setupSquare1000(i,j);
//						}
//
//
//
//	for(int i = 1; i < Mesh.getDimensions().x()-1; i++)
//					}
//				}
//				else
//				{
//		Index j = 0;
//		tmp = Sign(dofVector[Mesh.getCellIndex(CoordinatesType(i,j))]);
//
//
//		if(tmp == 0.0)
//		{}
//		else if(dofVector[Mesh.getCellIndex(CoordinatesType(i+1,j))]*tmp < 0.0 ||
//				dofVector[Mesh.getCellIndex(CoordinatesType(i-1,j))]*tmp < 0.0 ||
//				dofVector[Mesh.getCellIndex(CoordinatesType(i,j+1))]*tmp < 0.0 )
//		{}
//					if(dofVector[neighbourEntities.template getEntityIndex< 1,  0 >()] > 0)
//					{
//						if(dofVector[neighbourEntities.template getEntityIndex< 0,  1 >()] > 0)
//						{
//							if(dofVector[neighbourEntities.template getEntityIndex< 1,  1 >()] > 0)
//								setupSquare0111(i,j);
//							else
//			dofVector[Mesh.getCellIndex(CoordinatesType(i,j))] = tmp*INT_MAX;
//								setupSquare0110(i,j);
//						}
//
//	for(int i = 1; i < Mesh.getDimensions().x()-1; i++)
//						else
//						{
//		Index j = Mesh.getDimensions().y() - 1;
//		tmp = Sign(dofVector[Mesh.getCellIndex(CoordinatesType(i,j))]);
//
//
//		if(tmp == 0.0)
//		{}
//		else if(dofVector[Mesh.getCellIndex(CoordinatesType(i+1,j))]*tmp < 0.0 ||
//				dofVector[Mesh.getCellIndex(CoordinatesType(i-1,j))]*tmp < 0.0 ||
//				dofVector[Mesh.getCellIndex(CoordinatesType(i,j-1))]*tmp < 0.0 )
//		{}
//							if(dofVector[neighbourEntities.template getEntityIndex< 1,  1 >()] > 0)
//								setupSquare0101(i,j);
//							else
//			dofVector[Mesh.getCellIndex(CoordinatesType(i,j))] = tmp*INT_MAX;
//								setupSquare0100(i,j);
//						}
//
//	for(int j = 1; j < Mesh.getDimensions().y()-1; j++)
//					}
//					else
//					{
//		Index i = 0;
//		tmp = Sign(dofVector[Mesh.getCellIndex(CoordinatesType(i,j))]);
//
//
//		if(tmp == 0.0)
//		{}
//		else if(dofVector[Mesh.getCellIndex(CoordinatesType(i+1,j))]*tmp < 0.0 ||
//				dofVector[Mesh.getCellIndex(CoordinatesType(i,j+1))]*tmp < 0.0 ||
//				dofVector[Mesh.getCellIndex(CoordinatesType(i,j-1))]*tmp < 0.0 )
//		{}
//						if(dofVector[neighbourEntities.template getEntityIndex< 0,  1 >()] > 0)
//						{
//							if(dofVector[neighbourEntities.template getEntityIndex< 1,  1 >()] > 0)
//								setupSquare0011(i,j);
//							else
//			dofVector[Mesh.getCellIndex(CoordinatesType(i,j))] = tmp*INT_MAX;
//								setupSquare0010(i,j);
//						}
//
//	for(int j = 1; j < Mesh.getDimensions().y()-1; j++)
//						else
//						{
//		Index i = Mesh.getDimensions().x() - 1;
//		tmp = Sign(dofVector[Mesh.getCellIndex(CoordinatesType(i,j))]);
//
//
//		if(tmp == 0.0)
//		{}
//		else if(dofVector[Mesh.getCellIndex(CoordinatesType(i-1,j))]*tmp < 0.0 ||
//				dofVector[Mesh.getCellIndex(CoordinatesType(i,j+1))]*tmp < 0.0 ||
//				dofVector[Mesh.getCellIndex(CoordinatesType(i,j-1))]*tmp < 0.0 )
//		{}
//							if(dofVector[neighbourEntities.template getEntityIndex< 1,  1 >()] > 0)
//								setupSquare0001(i,j);
//							else
//			dofVector[Mesh.getCellIndex(CoordinatesType(i,j))] = tmp*INT_MAX;
//								setupSquare0000(i,j);
//						}
//					}
//				}
//
//
//	Index i = Mesh.getDimensions().x() - 1;
//	Index j = Mesh.getDimensions().y() - 1;
//
//	tmp = Sign(dofVector[Mesh.getCellIndex(CoordinatesType(i,j))]);
//	if(dofVector[Mesh.getCellIndex(CoordinatesType(i-1,j))]*tmp > 0.0 &&
//			dofVector[Mesh.getCellIndex(CoordinatesType(i,j-1))]*tmp > 0.0)
//
//		dofVector[Mesh.getCellIndex(CoordinatesType(i,j))] = tmp*INT_MAX;
//
//
//
//	j = 0;
//	tmp = Sign(dofVector[Mesh.getCellIndex(CoordinatesType(i,j))]);
//	if(dofVector[Mesh.getCellIndex(CoordinatesType(i-1,j))]*tmp > 0.0 &&
//			dofVector[Mesh.getCellIndex(CoordinatesType(i,j+1))]*tmp > 0.0)
//
//		dofVector[Mesh.getCellIndex(CoordinatesType(i,j))] = tmp*INT_MAX;
//
//
//
//	i = 0;
//	j = Mesh.getDimensions().y() -1;
//	tmp = Sign(dofVector[Mesh.getCellIndex(CoordinatesType(i,j))]);
//	if(dofVector[Mesh.getCellIndex(CoordinatesType(i+1,j))]*tmp > 0.0 &&
//			dofVector[Mesh.getCellIndex(CoordinatesType(i,j-1))]*tmp > 0.0)
//
//		dofVector[Mesh.getCellIndex(CoordinatesType(i,j))] = tmp*INT_MAX;
//
//
//
//	j = 0;
//	tmp = Sign(dofVector[Mesh.getCellIndex(CoordinatesType(i,j))]);
//	if(dofVector[Mesh.getCellIndex(CoordinatesType(i+1,j))]*tmp > 0.0 &&
//			dofVector[Mesh.getCellIndex(CoordinatesType(i,j+1))]*tmp > 0.0)
//
//		dofVector[Mesh.getCellIndex(CoordinatesType(i,j))] = tmp*INT_MAX;
//			}
//	}
	cout << "a" << endl;

	//data.setLike(dofVector2.getData());
	//data=dofVector2.getData();
@@ -336,7 +221,10 @@ template< typename MeshReal,
          typename Index >
bool tnlFastSweepingMap< tnlGrid< 2,MeshReal, Device, MeshIndex >, Real, Index > :: run()
{

	int cntr = 0;
	while(something_changed != 0)
	{
		something_changed = 0;
		for(Index i = 0; i < Mesh.getDimensions().x(); i++)
		{
			for(Index j = 0; j < Mesh.getDimensions().y(); j++)
@@ -375,6 +263,9 @@ bool tnlFastSweepingMap< tnlGrid< 2,MeshReal, Device, MeshIndex >, Real, Index >
		}

	/*---------------------------------------------------------------------------------------------------------------------------*/
		cntr++;
		cout << "Finished set of sweeps #" << cntr << "           " << something_changed << endl;
	}


//	data.setLike(dofVector2.getData());
@@ -397,9 +288,12 @@ void tnlFastSweepingMap< tnlGrid< 2,MeshReal, Device, MeshIndex >, Real, Index >

	this->Entity.setCoordinates(CoordinatesType(i,j));
	this->Entity.refresh();
	if(map[Entity.getIndex()] != 0.0)
	{
		tnlNeighbourGridEntityGetter<tnlGridEntity< MeshType, 2, tnlGridEntityNoStencilStorage >,2> neighbourEntities(Entity);

		Real value = dofVector2[Entity.getIndex()];
		Real im = abs(1.0/map[Entity.getIndex()]);
		Real a,b, tmp;

		if( i == 0 )
@@ -423,16 +317,21 @@ void tnlFastSweepingMap< tnlGrid< 2,MeshReal, Device, MeshIndex >, Real, Index >
		}


	if(fabs(a-b) >= h)
		tmp = fabsMin(a,b) + Sign(value)*h;
		if(fabs(a-b) >= im*h)
			tmp = fabsMin(a,b) + Sign(value)*im*h;
		else
		tmp = 0.5 * (a + b + Sign(value)*sqrt(2.0 * h * h - (a - b) * (a - b) ) );
			tmp = 0.5 * (a + b + Sign(value)*sqrt(2.0 * im * h * im * h - (a - b) * (a - b) ) );

		if(abs(value)-abs(tmp) > 0.0)
			something_changed = 1;

		dofVector2[Entity.getIndex()] = fabsMin(value, tmp);

//	if(dofVector2[Entity.getIndex()] > 1.0)
//		cout << value << "    " << tmp << " " << dofVector2[Entity.getIndex()] << endl;
	}
	else
	{
		dofVector2[Entity.getIndex()] = MAP_SOLVER_MAX_VALUE;
	}
}


+57 −393

File changed.

Preview size limit exceeded, changes collapsed.

Loading