Removed unused methods and atributes (65c2a48b) · Commits · TNL / 3D adaptive grid

Grid3DAdaptive.h

+5 −28

Original line number	Diff line number	Diff line
		@@ -81,13 +81,11 @@ template< typename Real, typename Device, typename Index >
		class GridA< 3, Real, Device, Index >
		{
		public:
		typedef Real RealType;
		typedef Real RealType; // For VTK writer
		typedef Device DeviceType;
		typedef Index GlobalIndexTyppe;
		typedef Containers::StaticVector< 3, Real > PointType;
		typedef Containers::StaticVector< 3, Index > CoordinatesType;
		typedef GridA< 3, Real, Devices::Host, Index > HostType;
		typedef GridA< 3, Real, Devices::Cuda, Index > CudaType;

		/* Structure that hold cells of the adaptive grid */
		/* Note that for now the adaptive cell features all template parameters of
		@@ -97,7 +95,6 @@ public:

		/* This structure is used to hold information (indecies) about adjacent cells
		*/
		typedef std::vector< Index > IndexArray;
		typedef std::vector< std::vector< short int > > DirectionsArray;

		using GlobalIndexType = Index;
		@@ -109,17 +106,9 @@ public:
		/* The array that holds information on the vertex */
		typedef std::vector< AdaptiveVertex< Index, Index > > VertexListType;

		/* Sparse matrix format */
		typedef std::pair< Index, Index > PairOfIndices;
		typedef std::map< short int, short int > UpTo4ShortInts;

		typedef std::map< PairOfIndices, UpTo4ShortInts > SparseMatrixOfShortInts;

		/* memory map - links every cell in the "big grid" to it's "children"*/
		typedef std::vector< std::vector< Index > > MemoryMap;

		typedef std::vector< Index > InactiveColumnsMap; // TODO: IndexArray

		static constexpr int
		getMeshDimension()
		{
		@@ -132,14 +121,10 @@ public:
		typedef EntityType< 3 > Cell;
		typedef EntityType< 0 > Vertex;

		/* VTK export related types */
		typedef std::vector< Cell > ArrayOfEntityCellsType;
		typedef std::vector< Vertex > ArrayOfEntityVerticesType;

		/* A way to move this into the private interface should be devised */
		/* VTK export related array*/
		ArrayOfEntityCellsType arrayOfEntityCells;
		ArrayOfEntityVerticesType arrayOfEntityVertices;
		std::vector< Cell > arrayOfEntityCells; // Dependency: VTK export
		std::vector< Vertex > arrayOfEntityVertices; // Dependency: VTK export

		GridA();

		@@ -209,6 +194,7 @@ public:
		getNumberOfVertices() const;

		/* Method used to sort all entries in the grid */
		/* TODO: is it useful? */
		void
		sortGrid();

		@@ -330,12 +316,6 @@ private:

		PointType
		returnCoordinatesOfVertex( Index vertexNumber ) const;
		CoordinatesType
		returnDiscreteCoordinatesOfVertex( Index vertexNumber ) const;

		Index
		calculateRefinementLevelBetween2Vertices( const AdaptiveVertex< Index, Index > vertex1,
		const AdaptiveVertex< Index, Index > vertex2 );

		/* Methods used to generate the dual mesh (virtual cells) */
		bool
		@@ -344,9 +324,6 @@ private:
		coarserCellHasVirtualCounterparts( Index cellIndex );

		/* More helper methods */
		Index
		returnIncrementationOfCellIndex( Index baseCellIndex, int x, int y, int z ) const;

		void
		addVerticesAdjacentToCell( Index adjacentCell,
		std::vector< int >& localIndicesOfVectors,
		@@ -393,7 +370,7 @@ protected:

		/* Current adjacent cells list */
		// TODO: avoid global variables, pass as a parameter when necessary
		IndexArray adjacentCellIndicesAllDirections;
		std::vector< Index > adjacentCellIndicesAllDirections;

		// NOTE: this contains directions of the adjacent cells in adjacentCellsIndices (adjacentCellIndicesAllDirections)
		std::vector< short int > directionsOfAdjacentCellsFlat;

Grid3DAdaptive_impl.h

+1 −223

Original line number	Diff line number	Diff line
		@@ -365,7 +365,7 @@ GridA< 3, Real, Device, Index >::refineGrid( GridRefiner< Real, Device, Index >&

		/* On each refinement level, we iterate through all the cells in the grid
		*/
		InactiveColumnsMap inactiveColumnsMap;
		std::vector< Index > inactiveColumnsMap;
		for( j = 0; j < numberOfCellsInGridNow; j++ ) {
		/* Use the passed instance of grid refiner to evaluate wether we refine
		* a cell or not */
		@@ -896,25 +896,6 @@ GridA< 3, Real, Device, Index >::calculateAdjacencyMatrix( int mode )
		/* !!!DEBUG!!! - NOTE: now the counter is behind one!*/
		addCentreCellToAdjacencyList( k );
		remapDirections();

		/* New code - START HERE! check this!*/

		SparseMatrixOfShortInts miniMatrixToPush;

		Index numberOfAdjacentCells = directionsOfAdjacentCellsFlat.size();

		for( Index w = 0; w < numberOfAdjacentCells; w++ ) {
		/* all the data should be ready in the array - just load it up */
		auto itr = miniMatrixToPush.find( std::make_pair( k, adjacentCellIndicesAllDirections[ w ] ) );
		if( itr == miniMatrixToPush.end() ) {
		UpTo4ShortInts intsToPutInMatrix;
		intsToPutInMatrix[ 0 ] = directionsOfAdjacentCellsFlat[ w ] + 1;
		miniMatrixToPush[ std::make_pair( k, adjacentCellIndicesAllDirections[ w ] ) ] = intsToPutInMatrix;
		}
		else {
		itr->second[ itr->second.size() ] = directionsOfAdjacentCellsFlat[ w ] + 1;
		}
		}
		}
		}

		@@ -2062,80 +2043,6 @@ GridA< 3, Real, Device, Index >::cellIsOnBoundary( const AdaptiveCell< Index, In
		return cellOnBoundaryInGrid;
		}

		template< typename Real, typename Device, typename Index >
		Index
		GridA< 3, Real, Device, Index >::calculateRefinementLevelBetween2Vertices( const AdaptiveVertex< Index, Index > vertex1,
		const AdaptiveVertex< Index, Index > vertex2 )
		{
		/* Load the two vectors into an array */
		std::vector< AdaptiveVertex< Index, Index > > arrayOfTwoVectors;
		arrayOfTwoVectors.push_back( vertex1 );
		arrayOfTwoVectors.push_back( vertex2 );

		/* first, we compute the maximal non-zero refinement index for our set of
		* vectors */
		Index firstNonZeroBit = -1;
		Index upperBound = sizeof( Index ) * 8;

		for( Index j = 0; j < 2; j++ ) {
		for( Index i = 0; i < ( upperBound / 2 ); i++ ) {
		/* if a non-zero bit is found and it is larger than the current largest
		* bit, it is recorded */
		if( ( getBit( arrayOfTwoVectors[ j ].refinementPos, 2 * i ) == 1 ) && ( 2 * i > firstNonZeroBit ) ) {
		firstNonZeroBit = 2 * i;
		}
		}
		}

		/* keep the case where fisrtNonZeroBit = -1 in mind! */

		/* The resolution gives information on how many steps of refinement it took
		* to get the finest vertex in the array*/
		Index resolution = -1;
		if( firstNonZeroBit != -1 ) {
		resolution = (Index) ( firstNonZeroBit / 2 );
		}

		std::vector< Index > normsOfIndices( arrayOfTwoVectors.size(), 0 );

		Index counter = resolution;
		for( Index j = 0; j < 2; j++ ) {
		counter = resolution;
		while( counter >= 0 ) {
		/* Only the x position is taken into account */
		if( getBit( arrayOfTwoVectors[ j ].refinementPos, counter * 2 ) ) {
		normsOfIndices[ j ] = normsOfIndices[ j ] + pow( 2, ( resolution - counter ) );
		}
		counter = counter - 1;
		}
		}

		/* In any case, the counter is -1 at this point, now is the right time to
		* increment based on meshPos */
		for( Index j = 0; j < 2; j++ ) {
		/* First the x and y coordinates are calculated for each vector */
		Index yCoordinate = (Index) ( (Real) arrayOfTwoVectors[ j ].meshPos / (Real) ( dimensions.x() + 1 ) );
		Index xCoordinate = (Index) ( (Real) arrayOfTwoVectors[ j ].meshPos - (Real) ( yCoordinate * ( dimensions.x() + 1 ) ) );

		normsOfIndices[ j ] = normsOfIndices[ j ] + pow( 2, resolution + 1 ) * ( xCoordinate );
		}

		/* Calculate the refinement level between the vertices and return it */
		Index differenceInNorms = abs( normsOfIndices[ 1 ] - normsOfIndices[ 0 ] );
		Index refinementLevel = -1;
		if( resolution == -1 ) {
		refinementLevel = 0;
		}
		else {
		for( Index k = 0; k <= ( resolution + 1 ); k++ ) {
		if( differenceInNorms == pow( 2, k ) ) {
		refinementLevel = resolution + 1 - k;
		}
		}
		}
		return refinementLevel;
		}

		template< typename Real, typename Device, typename Index >
		void
		GridA< 3, Real, Device, Index >::prepareForVtkExport()
		@@ -2231,73 +2138,6 @@ GridA< 3, Real, Device, Index >::returnCoordinatesOfVertex( Index vertexNumber )
		return coordinatesToReturn;
		}

		template< typename Real, typename Device, typename Index >
		typename GridA< 3, Real, Device, Index >::CoordinatesType
		GridA< 3, Real, Device, Index >::returnDiscreteCoordinatesOfVertex( Index vertexNumber ) const
		{
		/* First we pull the vertex from the list */
		AdaptiveVertex< Index, Index > adaptiveVertex = listOfVertices[ vertexNumber ];

		/* The max discrete size of the grid is scaled so that the refinement level
		* is accommodated */
		Index sizeOfGridInX = dimensions.x() * pow( 2, maximumRefinementLevel );
		Index sizeOfGridInY = dimensions.y() * pow( 2, maximumRefinementLevel );
		Index sizeOfGridInZ = dimensions.z() * pow( 2, maximumRefinementLevel );

		/* This could very well be just written as pow(2,maximumRefinementLevel) but
		* we keep this to maintain structural similarities with the previous method
		* -
		* TODO: merge them */
		Index sizeOfStepInX = sizeOfGridInX / dimensions.x();
		Index sizeOfStepInY = sizeOfGridInY / dimensions.y();
		Index sizeOfStepInZ = sizeOfGridInZ / dimensions.z();

		/* We do not need this line since we have maximum refinement level */
		// Index upperBoundOfLoop = sizeof(Index) * (CHAR_BIT/2);

		Index numberOfVerticesPerLayer = ( dimensions.x() + 1 ) * ( dimensions.y() + 1 );

		Index zCoordinateCoarse = (Index) ( adaptiveVertex.meshPos / numberOfVerticesPerLayer );

		Index projectionToLowestLayer = adaptiveVertex.meshPos % numberOfVerticesPerLayer;
		Index yCoordinateCoarse = (Index) ( projectionToLowestLayer / ( dimensions.x() + 1 ) );
		Index xCoordinateCoarse = projectionToLowestLayer - yCoordinateCoarse * ( dimensions.x() + 1 );

		/* Add the base value of the coarse grid */
		Index xCoordinateToReturn = sizeOfStepInX * xCoordinateCoarse;
		Index yCoordinateToReturn = sizeOfStepInY * yCoordinateCoarse;
		Index zCoordinateToReturn = sizeOfStepInZ * zCoordinateCoarse;

		Index currentStepSizeX = ( sizeOfStepInX / 2 );
		Index currentStepSizeY = ( sizeOfStepInY / 2 );
		Index currentStepSizeZ = ( sizeOfStepInZ / 2 );

		for( Index i = 0; i < maximumRefinementLevel; i++ ) {
		/* The coordinates are incremented according using the bit string
		* adaptiveVertex.refinementPos */
		if( getBit( adaptiveVertex.refinementPos, 3 * i ) ) {
		xCoordinateToReturn += currentStepSizeX;
		}

		if( getBit( adaptiveVertex.refinementPos, 3 * i + 1 ) ) {
		yCoordinateToReturn += currentStepSizeY;
		}

		if( getBit( adaptiveVertex.refinementPos, 3 * i + 2 ) ) {
		zCoordinateToReturn += currentStepSizeZ;
		}
		currentStepSizeX = ( currentStepSizeX / 2 );
		currentStepSizeY = ( currentStepSizeY / 2 );
		currentStepSizeZ = ( currentStepSizeZ / 2 );
		}

		PointType coordinatesToReturn;
		coordinatesToReturn.x() = xCoordinateToReturn;
		coordinatesToReturn.y() = yCoordinateToReturn;
		coordinatesToReturn.z() = zCoordinateToReturn;
		return coordinatesToReturn;
		}

		/* a more complete version of the previous method */
		template< typename Real, typename Device, typename Index >
		AdaptiveCell< Index, Index, Index >
		@@ -2801,22 +2641,6 @@ GridA< 3, Real, Device, Index >::calculateComponentAdjacencyMatrix()
		* extra data */
		addCentreCellToAdjacencyList( k );
		remapDirections();

		SparseMatrixOfShortInts miniMatrixToPush;

		Index numberOfAdjacentCells = directionsOfAdjacentCellsFlat.size();
		for( Index w = 0; w < numberOfAdjacentCells; w++ ) {
		/* all the data should be ready in the array - just load it up */
		auto itr = miniMatrixToPush.find( std::make_pair( k, adjacentCellIndicesAllDirections[ w ] ) );
		if( itr == miniMatrixToPush.end() ) {
		UpTo4ShortInts intsToPutInMatrix;
		intsToPutInMatrix[ 0 ] = directionsOfAdjacentCellsFlat[ w ] + 1;
		miniMatrixToPush[ std::make_pair( k, adjacentCellIndicesAllDirections[ w ] ) ] = intsToPutInMatrix;
		}
		else {
		itr->second[ itr->second.size() ] = directionsOfAdjacentCellsFlat[ w ] + 1;
		}
		}
		}
		}

		@@ -2928,52 +2752,6 @@ GridA< 3, Real, Device, Index >::recalculateMemoryMapsByLayer()
		}
		}

		/* more helper methods */
		template< typename Real, typename Device, typename Index >
		Index
		GridA< 3, Real, Device, Index >::returnIncrementationOfCellIndex( Index baseCellIndex, int x, int y, int z ) const
		{
		/* with regards to the use of this method we return -1 if the incrementation
		* of the direction is out of bounds */
		Index valueToReturn = baseCellIndex;

		Index coordinateValueOfZ = (Index) ( valueToReturn / ( dimensions.y() * dimensions.x() ) );
		Index coordinateValueOfY =
		(Index) ( ( valueToReturn - coordinateValueOfZ * dimensions.y() * dimensions.x() ) / dimensions.x() );
		Index coordinateValueOfX =
		(Index) ( valueToReturn - coordinateValueOfY * dimensions.x() - coordinateValueOfZ * dimensions.y() * dimensions.x() );
		// std::cout << "coord vals: " << coordinateValueOfX << ", " <<
		// coordinateValueOfY << ", " << coordinateValueOfZ <<std::endl;

		/* increment if it's still in bounds */
		if( coordinateValueOfX + x >= 0 && coordinateValueOfX + x < dimensions.x() ) {
		valueToReturn = valueToReturn + x;
		}
		else {
		// std::cout << " X out of bounds! " << std::endl;
		return -1;
		}

		if( coordinateValueOfY + y >= 0 && coordinateValueOfY + y < dimensions.y() ) {
		valueToReturn = valueToReturn + y * dimensions.x();
		}
		else {
		// std::cout << " Y out of bounds! " << std::endl;
		return -1;
		}

		if( coordinateValueOfZ + z >= 0 && coordinateValueOfZ + z < dimensions.z() ) {
		valueToReturn = valueToReturn + z * dimensions.x() * dimensions.y();
		}
		else {
		// std::cout << " Z out of bounds! " << std::endl;
		return -1;
		}

		// std::cout << " returning value: " << valueToReturn << std::endl;
		return valueToReturn;
		}

		/* remaps directions for adjacency matrix export - adds directions for export */
		template< typename Real, typename Device, typename Index >
		void