Removed unused files

#pragma once

namespace TNL {

namespace Meshes {

/* Index1 represents the cell of the index in the original grid

 * Index2 represents the refinement of a particular cell

 * Index3 represents the refinement level of the cell */

template< typename Index1 = unsigned int, typename Index2 = unsigned int, typename Index3 = unsigned int >

class AdaptiveCell

{

public:

   Index1 meshPos;

   Index2 refinementPos;

   Index3 levelOfRef;

   AdaptiveCell( Index1 meshPos = 0, Index2 refinementPos = 0, Index3 levelOfRef = 0 )

   : meshPos( meshPos ), refinementPos( refinementPos ), levelOfRef( levelOfRef )

   {}

};

/* The datatype used to store vertex information */

template< typename Index1 = unsigned int, typename Index2 = unsigned int >

class AdaptiveVertex

{

public:

   Index1 meshPos;

   Index2 refinementPos;

};

}  // namespace Meshes

}  // namespace TNL

#pragma once

namespace TNL {

namespace Meshes {

template< typename Grid, typename Entity >

class AdaptiveEntityGetter

{};

template< typename Grid >

class AdaptiveEntityGetter< Grid, typename Grid ::Vertex >

{

   using Vertex = typename Grid ::Vertex;

   using Index = typename Grid ::GlobalIndexType;

public:

   static Vertex

   getEntity( const Grid& grid, Index vertexToPull )

   {

      return grid.arrayOfEntityVertices[ vertexToPull ];

   }

};

/* Spot for one more specilization - spec. for cells */

template< typename Grid >

class AdaptiveEntityGetter< Grid, typename Grid ::Cell >

{

   using Cell = typename Grid ::Cell;

   using Index = typename Grid ::GlobalIndexType;

public:

   static Cell

   getEntity( const Grid& grid, Index cellToPull )

   {

      return grid.arrayOfEntityCells[ cellToPull ];

   }

};

}  // namespace Meshes

}  // namespace TNL

/***************************************************************************

                          Grid2D.h  -  description

                             -------------------

    begin                : Feb 13, 2015

    copyright            : (C) 2015 by Tomas Oberhuber

    email                : tomas.oberhuber@fjfi.cvut.cz

 ***************************************************************************/

/* See Copyright Notice in tnl/Copyright */

#pragma once

#include <TNL/Containers/StaticVector.h>

#include <TNL/Matrices/SparseMatrix.h>

#include <TNL/Meshes/Writers/VTKWriter.h>

#include <TNL/Meshes/Writers/VTUWriter.h>

/* 	The vector structure is used to hold the elements of the refined grid,

   AdaptiveCell.h holds the sctructure that represents a single adaptive cell */

#include <climits>

#include <typeinfo>

#include <vector>

#include "AdaptiveCell.h"

#include "AdaptiveEntityGetter.h"

#include "AdaptiveGridEntity.h"

#include "GridRefiner.h"

// Debugger include

#include <csignal>

#include <iostream>

#include <unordered_map>

#include <utility>

namespace TNL {

namespace Meshes {

/* Enumeration used to label directions when calculating the adjacent cells */

enum baseDirections

{

   X_NEGATIVE,

   X_POSITIVE,

   Y_NEGATIVE,

   Y_POSITIVE,

   Z_NEGATIVE,

   Z_POSITIVE

};

/* all possible directions are ordered counter-clock-wise. Z is put first

 * because we view everything in a layer kinda way */

enum allDirections

{

   Z_NEG,

   Z_NEG_X_NEG,

   Z_NEG_X_NEG_Y_NEG,

   Z_NEG_Y_NEG,

   Z_NEG_X_POS_Y_NEG,

   Z_NEG_X_POS,

   Z_NEG_X_POS_Y_POS,

   Z_NEG_Y_POS,

   Z_NEG_X_NEG_Y_POS,

   X_NEG,

   X_NEG_Y_NEG,

   Y_NEG,

   X_POS_Y_NEG,

   X_POS,

   X_POS_Y_POS,

   Y_POS,

   X_NEG_Y_POS,

   Z_POS,

   Z_POS_X_NEG,

   Z_POS_X_NEG_Y_NEG,

   Z_POS_Y_NEG,

   Z_POS_X_POS_Y_NEG,

   Z_POS_X_POS,

   Z_POS_X_POS_Y_POS,

   Z_POS_Y_POS,

   Z_POS_X_NEG_Y_POS

};

template< typename Real, typename Device, typename Index >

class AdaptiveGrid< 3, Real, Device, Index >

{

public:

   typedef Real RealType;  // For VTK writer

   typedef Index GlobalIndexTyppe;

   typedef Containers::StaticVector< 3, Real > PointType;

   typedef Containers::StaticVector< 3, Index > CoordinatesType;

   /* Structure that hold cells of the adaptive grid */

   /* Note that for now the adaptive cell features all template parameters of

    * type Index. In the future a new template parameter should be created so

    * that the last index does not have to have the same size as the others */

   typedef std::vector< AdaptiveCell< Index, Index, Index > > CellContainer;

   /* This structure is used to hold information (indecies) about adjacent cells

    */

   typedef std::vector< std::vector< short int > > DirectionsArray;

   using GlobalIndexType = Index;  // For VTK writer

   /* Adjacency matrix that holds the information about the adjacency of cells

    */

   using AdjacencyMatrix = TNL::Matrices::SparseMatrix< short int, Device, Index >;

   /* The array that holds information on the vertex */

   typedef std::vector< AdaptiveVertex< Index, Index > > VertexListType;

   /* memory map - links every cell in the "big grid" to it's "children"*/

   typedef std::vector< std::vector< Index > > MemoryMap;

   static constexpr int

   getMeshDimension()

   {

      return 3;

   };

   template< int EntityDimension >

   using EntityType = AdaptiveGridEntity< AdaptiveGrid, EntityDimension >;

   typedef EntityType< 3 > Cell;

   typedef EntityType< 0 > Vertex;

   /* A way to move this into the private interface should be devised */

   /* VTK export related array*/

   std::vector< Cell > arrayOfEntityCells;       // Dependency: VTK export

   std::vector< Vertex > arrayOfEntityVertices;  // Dependency: VTK export

   AdaptiveGrid();

   void

   setDimensions( const Index xSize, const Index ySize, const Index zSize );

   void

   setDimensions( const CoordinatesType& dimensions );

   __cuda_callable__

   const CoordinatesType&

   getDimensions() const;

   void

   setDomain( const PointType& origin, const PointType& proportions );

   __cuda_callable__

   inline const PointType&

   getOrigin() const;

   __cuda_callable__

   inline const PointType&

   getProportions() const;

   template< int EntityDimension >

   __cuda_callable__

   Index

   getEntitiesCount() const;

   template< typename Entity >

   __cuda_callable__

   inline Index

   getEntitiesCount() const;

   template< typename Entity >

   __cuda_callable__

   inline Entity

   getEntity( const Index& entityIndex ) const;

   template< typename Entity >

   __cuda_callable__

   inline Index

   getEntityIndex( const Entity& entity ) const;

   /* Adaptive grid methods */

   /* TODO: is it useful? */

   void

   appendElement( const AdaptiveCell< Index, Index, Index >& cellToAppend );

   AdaptiveCell< Index, Index, Index >

   getCell( Index ind ) const;

   void

   generateInitialGrid();

   void

   refineGrid( GridRefiner< Real, Device, Index >& gridRefiner, int numberOfLevels );

   /* This method returns the current number of cells - usefull for potential

    * loops */

   Index

   getNumberOfCells() const;

   Index

   getNumberOfVertices() const;

   /* Methods used for calculating the adjacent cells - these should be moved

    * into private later on */

   /* New helper method for 3D */

   /* TODO: move this out of public space after the debug getCoordinatesOnBasicGridphase */

   void

   getAdjacentCellsUlt( Index cell,

                        std::vector< Index >& adjacentCellIndicesAllDirections,

                        std::vector< Index >& directionsOfAdjacentCellsFlat );

   /* helper methods for component adjacency matrices */

   /* Method that adds a vertex artificially to the grid */

   void

   addVertex( Index meshPos, Index refinementPos );

   /* Method that must be called before the grid may be exported into VTK */

   void

   prepareForExport();

   void

   exportGrid( const std::string& fileName ) const;

   /* Just print all inforamtion about grid to console */

   void

   displayDebugInfo() const;

   void

   refineSpecificCell( Index cellIndex );

private:

   /* Comment on this better */

   /* ?Why cant I use CellContainer */

   std::vector< AdaptiveCell< Index, Index, Index > >

   refineCell( AdaptiveCell< Index, Index, Index >& adaptiveCell );

   /* Ult variant */

   AdaptiveCell< Index, Index, Index >

   getVirtualAdjacentCell( const AdaptiveCell< Index, Index, Index >& baseCell, short int direction ) const;

   /* Methods used in the interior of the getAdjacentCells(...) method */

   AdaptiveCell< Index, Index, Index >

   getAdjacentSameLevel( const AdaptiveCell< Index, Index, Index >& baseCell, int direction ) const;

   AdaptiveCell< Index, Index, Index >

   adjacentCellSameLevelTwoDirections( const AdaptiveCell< Index, Index, Index >& baseCell,

                                       int direction1,

                                       int direction2 ) const;

   AdaptiveCell< Index, Index, Index >

   adjacentCellSameLevelThreeDirections( const AdaptiveCell< Index, Index, Index >& baseCell,

                                         int direction1,

                                         int direction2,

                                         int direction3 ) const;

   /* NOTE TO SELF: Use this method to pimp up (refactor) some of the code */

   std::vector< AdaptiveCell< Index, Index, Index > >

   getParentCells( const AdaptiveCell< Index, Index, Index >& baseCell ) const;

   /* NOTE: rewrite the following functions using typedefs, I often use

    * std::vector<Index> instead of IndexArray etc. this is good for now, as it

    * improves readabiliy. */

   /* Delete the i-th column of the cell-vertex incidence matrix */

   void

   deleteBatchOfColumnsOfIncidenceMat( std::vector< Index >& columnsToDelete );

   std::vector< AdaptiveVertex< Index, Index > >

   calculateNewVertices( const std::vector< AdaptiveVertex< Index, Index > >& fourBaseVertices, int refinementLevel ) const;

   /* This method returns the indices for the base vertices as the first 1 in a

    * particular order, all the other vertices follow */

   std::vector< Index >

   returnArrayOfBaseVertices( const std::vector< Index >& arrayOfInputIndices, Index levelOfRefinement ) const;

   Index

   getIndexOfVertex( std::vector< Index > vectorsToChoseFrom, const AdaptiveVertex< Index, Index >& adaptiveVertex ) const;

   Index

   getLowerBotLeftIndex( const std::vector< AdaptiveVertex< Index, Index > >& inputVectors ) const;

   void

   recalculateIncidenceMatrixFundamental( const std::vector< Index >& indexArray );

   void

   reclaculateIncidenceMatrixAdditional( const std::vector< Index >& indexArray,

                                         int refinementLevel,

                                         std::vector< Index >& adjacentCellIndicesAllDirections,

                                         std::vector< Index >& directionsOfAdjacentCellsFlat );

   /* Helper methods for working with diagonal adjacency */

   bool

   cellIsOnBoundary( const AdaptiveCell< Index, Index, Index >& baseCell, int direction ) const;

   PointType

   returnCoordinatesOfVertex( Index vertexNumber ) const;

   /* More helper methods */

   void

   addVerticesAdjacentToCell( Index adjacentCell,

                              std::vector< int >& localIndicesOfVectors,

                              std::unordered_map< int, Index >& localToGlobalVertexMap,

                              std::vector< int >& localVectorFlags );

protected:

   void

   setCapacity( Index capacity );

   // Maybe do: clean up the set dimensions methods< Del stuff?

   // Next up: Implement vector that will reserve enough memory (use a method

   // for this) Implement a refine a particular cell method = method that will

   // generate elements for the refinement array > Refine Element Implement a

   // refine methrowLengthVectorod that does one step of refinement with a

   // particular threshhold

   // adjacency cell-to-cell (value: direction, column index: cell)

   CoordinatesType dimensions;

   PointType origin, proportions;

   /* The structure that holds the cells of the adaptive grid. The structure

      gridCellsTemp holds the chronological order in which the cells were added.

      The structure gridCells holds the sorted version of the data */

   // the main container - contains the cells

   CellContainer gridCellsTemp;

   /* List that holds all the vectors of the grid */

   VertexListType listOfVertices;

   /* Variable that tracks maximum refinement level */

   int maximumRefinementLevel;

   // TODO: rename MemoryMap to something else

   // MemoryMap data structures -  for optimal lookup to adjacent cells - performance critical

   // (each "row" contains all smaller cells of the given cell on coarse level)

   /* Memory map - holds information about the hierarchy */

   MemoryMap memoryMap;        // this is the memory map for cells

   MemoryMap vertexMemoryMap;  // this mapps from cells to vertices

   template< typename, int >

   friend class GridEntityGetter;

};

}  // namespace Meshes

}  // namespace TNL

#include "AdaptiveGrid_impl.h"

#pragma once

#include <TNL/Containers/StaticVector.h>

#include <TNL/Meshes/Topologies/Hexahedron.h>

namespace TNL {

namespace Meshes {

// forward declaration

template< int Dimension, typename Real, typename Device, typename Index >

class AdaptiveGrid;

template< typename Grid, int EntityDimension >

class AdaptiveGridEntity

{};

template< int Dimension, typename Real, typename Device, typename Index, int EntityDimension >

class AdaptiveGridEntity< Meshes::AdaptiveGrid< Dimension, Real, Device, Index >, EntityDimension >

{};

/* Specialization for cells */

template< int Dimension, typename Real, typename Device, typename Index >

class AdaptiveGridEntity< Meshes::AdaptiveGrid< Dimension, Real, Device, Index >, Dimension >

{

   std::vector< Index > ListOfVertices;

public:

   using EntityTopology = Meshes::Topologies::Hexahedron;

   AdaptiveGridEntity();

   void

   setVertices( const std::vector< Index >& inputVertices );

   void

   printContent() const;

   template< int DimensionOfEntity >

   Index

   getSubentityIndex( const int& vertexNumber ) const;

   static constexpr int

   getVerticesCount()

   {

      return 8;

   }

   static constexpr int

   getEntityDimension()

   {

      return 3;

   }

   template< int Dimension2 >

   static constexpr int

   getSubentitiesCount()

   {

      static_assert( Dimension2 == 0, "not implemented for Dimension > 0" );

      return 8;

   }

};

/****

 * Specialization for vertices

 */

template< int Dimension, typename Real, typename Device, typename Index >

class AdaptiveGridEntity< Meshes::AdaptiveGrid< Dimension, Real, Device, Index >, 0 >

{

   typedef Containers::StaticVector< 3, Real > PointType;

   PointType coordinatesOfVertex;

public:

   AdaptiveGridEntity();

   void

   setCoordinates( const Containers::StaticVector< 3, Real >& inputCoordinates );

   void

   printContent() const;

   PointType

   getPoint() const;

   static constexpr int

   getEntityDimension()

   {

      return 0;

   };

};

/* The following "implementation" content, this should be latter separated into

 * an "impl.h" file */

template< int Dimension, typename Real, typename Device, typename Index >

AdaptiveGridEntity< Meshes::AdaptiveGrid< Dimension, Real, Device, Index >, Dimension >::AdaptiveGridEntity()

{}

template< int Dimension, typename Real, typename Device, typename Index >

void

AdaptiveGridEntity< Meshes::AdaptiveGrid< Dimension, Real, Device, Index >, Dimension >::setVertices(

   const std::vector< Index >& inputVertices )

{

   /* We assume that the vertices are ordered in a particular way, only the

    * first 4 define the corners of the cell */

   /* They need to be reordered since vtk exports work with a fixed

    * (counterclockwise) orientation */

   this->ListOfVertices.push_back( inputVertices[ 0 ] );

   this->ListOfVertices.push_back( inputVertices[ 1 ] );

   this->ListOfVertices.push_back( inputVertices[ 3 ] );

   this->ListOfVertices.push_back( inputVertices[ 2 ] );

   this->ListOfVertices.push_back( inputVertices[ 4 ] );

   this->ListOfVertices.push_back( inputVertices[ 5 ] );

   this->ListOfVertices.push_back( inputVertices[ 7 ] );

   this->ListOfVertices.push_back( inputVertices[ 6 ] );

}

template< int Dimension, typename Real, typename Device, typename Index >

void

AdaptiveGridEntity< Meshes::AdaptiveGrid< Dimension, Real, Device, Index >, Dimension >::printContent() const

{

   std::cout << "The fundamental indices of the cell are (in counterclockwise "

                "order): ";

   for( int i = 0; i < 4; i++ ) {

      std::cout << this->ListOfVertices[ i ] << ", ";

   }

   std::cout << std::endl;

}

template< int Dimension, typename Real, typename Device, typename Index >

AdaptiveGridEntity< Meshes::AdaptiveGrid< Dimension, Real, Device, Index >, 0 >::AdaptiveGridEntity()

{}

template< int Dimension, typename Real, typename Device, typename Index >

void

AdaptiveGridEntity< Meshes::AdaptiveGrid< Dimension, Real, Device, Index >, 0 >::printContent() const

{

   std::cout << "The coordinates of the vertex are: " << this->coordinatesOfVertex.x() << ", " << this->coordinatesOfVertex.y()

             << std::endl;

}

template< int Dimension, typename Real, typename Device, typename Index >

void

AdaptiveGridEntity< Meshes::AdaptiveGrid< Dimension, Real, Device, Index >, 0 >::setCoordinates(

   const Containers::StaticVector< 3, Real >& inputCoordinates )

{

   this->coordinatesOfVertex.x() = inputCoordinates.x();

   this->coordinatesOfVertex.y() = inputCoordinates.y();

   this->coordinatesOfVertex.z() = inputCoordinates.z();

}

template< int Dimension, typename Real, typename Device, typename Index >

template< int DimensionOfEntity >

Index

AdaptiveGridEntity< Meshes::AdaptiveGrid< Dimension, Real, Device, Index >, Dimension >::getSubentityIndex(

   const int& vertexNumber ) const

{

   return this->ListOfVertices[ vertexNumber ];

}

template< int Dimension, typename Real, typename Device, typename Index >

typename AdaptiveGridEntity< Meshes::AdaptiveGrid< Dimension, Real, Device, Index >, 0 >::PointType

AdaptiveGridEntity< Meshes::AdaptiveGrid< Dimension, Real, Device, Index >, 0 >::getPoint() const

{

   return coordinatesOfVertex;

}

}  // namespace Meshes

}  // namespace TNL

#include <vector>

#include <iostream>

#ifdef HAVE_GTEST

   #include <gtest/gtest.h>

#endif

#include <TNL/Containers/StaticVector.h>

#include "AdaptiveGrid.h"

#include "AdaptiveCell.h"

#include "GridRefiner.h"

using namespace std;

#ifdef HAVE_GTEST

TEST( AdaptiveGrid, CreatingSimpleGrid )

{

   TNL::Meshes::AdaptiveGrid< 3, float, TNL::Devices::Host, long int > grid;

   TNL::Containers::StaticVector< 3, float > origin{0, 0, 0};

   TNL::Containers::StaticVector< 3, float > dimensions{2, 1, 1};

   TNL::Containers::StaticVector< 3, float > proportions = dimensions;

   grid.setDimensions( dimensions.x(), dimensions.y(), dimensions.z() );

   grid.setDomain( origin, proportions );

   grid.generateInitialGrid();

   ASSERT_EQ( grid.getNumberOfCells(), 2 );

   ASSERT_EQ( grid.getDimensions().x(), dimensions.x() );

   ASSERT_EQ( grid.getDimensions().y(), dimensions.y() );

   ASSERT_EQ( grid.getDimensions().z(), dimensions.z() );

   ASSERT_EQ( grid.getNumberOfVertices(), 12 );

   for( int i = 0; i < grid.getNumberOfCells(); ++i ) {

      ASSERT_EQ( grid.getCell( i ).levelOfRef, 0 );

      ASSERT_EQ( grid.getCell( i ).meshPos, i );

      ASSERT_EQ( grid.getCell( i ).refinementPos, 0 );

   }

   vector<int> results(grid.getNumberOfCells(), 0);

   for (int i = 0; i < grid.getNumberOfCells(); ++i) {

      std::vector< long int > adjacentCellIndicesAllDirections;

      std::vector< long int > directionsOfAdjacentCellsFlat; // Unused

      auto cell = grid.getCell(i);

      grid.getAdjacentCellsUlt(i, adjacentCellIndicesAllDirections, directionsOfAdjacentCellsFlat);

      for (auto cellIdx : adjacentCellIndicesAllDirections) {

         ++results.at(i);

      }

   }

   ASSERT_EQ(results.at(0), 1);

   ASSERT_EQ(results.at(1), 1);

   ASSERT_EQ(results.size(), 2);

   /* grid export*/

   grid.prepareForExport();

   grid.exportGrid( "testSimpleGrid.vtu" );

}

TEST( AdaptiveGrid, Refinement1Level )

{

   TNL::Meshes::AdaptiveGrid< 3, float, TNL::Devices::Host, long int > grid;

   TNL::Containers::StaticVector< 3, float > origin{0, 0, 0};

   TNL::Containers::StaticVector< 3, float > dimensions{2, 2, 2};

   TNL::Containers::StaticVector< 3, float > proportions = dimensions;

   grid.setDimensions( dimensions.x(), dimensions.y(), dimensions.z() );

   grid.setDomain( origin, proportions );

   grid.generateInitialGrid();

   TNL::Meshes::GridRefiner< float, TNL::Devices::Host, long int > refiner;

   refiner.setDomain( origin, proportions );

   refiner.setDimensions( dimensions );

   /* Without refinement */

   ASSERT_EQ( grid.getNumberOfCells(), 8 );

   ASSERT_EQ( grid.getDimensions().x(), dimensions.x() );

   ASSERT_EQ( grid.getDimensions().y(), dimensions.y() );

   ASSERT_EQ( grid.getDimensions().z(), dimensions.z() );

   ASSERT_EQ( grid.getNumberOfVertices(), 27 );

   for( int i = 0; i < grid.getNumberOfCells(); ++i ) {

      ASSERT_EQ( grid.getCell( i ).levelOfRef, 0 );

      ASSERT_EQ( grid.getCell( i ).meshPos, i );

      ASSERT_EQ( grid.getCell( i ).refinementPos, 0 );

   }

   refiner.setGridRefPattern( 2 );

   grid.refineGrid( refiner, 1 );

   /* After refinement */

   ASSERT_EQ( grid.getNumberOfCells(), 36 );  // 4 upper cubes + 4 cubes divided by 8 smaller cubes

   ASSERT_EQ( grid.getDimensions().x(), dimensions.x() );

   ASSERT_EQ( grid.getDimensions().y(), dimensions.y() );

   ASSERT_EQ( grid.getDimensions().z(), dimensions.z() );

   ASSERT_EQ( grid.getNumberOfVertices(), 84 );

   /* First 4 cells are without refinement */

   for( size_t i = 0; i < 4; ++i ) {

      std::vector< long int > adjacentCellIndicesAllDirections;

      std::vector< long int > directionsOfAdjacentCellsFlat; // Unused

      ASSERT_EQ( grid.getCell( i ).levelOfRef, 0 );

      ASSERT_EQ( grid.getCell( i ).refinementPos, 0 );

      grid.getAdjacentCellsUlt(i, adjacentCellIndicesAllDirections, directionsOfAdjacentCellsFlat);

      ASSERT_EQ( adjacentCellIndicesAllDirections.size(), 12 );

   }

   /* Other cells are refined, every cell is divided by 8 smaller cells */

   for( size_t i = 4; i < grid.getNumberOfCells(); i += 8 ) {

      for( size_t j = 0; j < 8; ++j ) {

         ASSERT_EQ( grid.getCell( i + j ).levelOfRef, 1 );

         ASSERT_EQ( grid.getCell( i + j ).refinementPos, j );