Removed VTKReader_libvtk.h and MeshReaderTest.h

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

                          VTKReader_libvtk.h  -  description

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

    begin                : Nov 6, 2016

    copyright            : (C) 2016 by Tomas Oberhuber et al.

    email                : tomas.oberhuber@fjfi.cvut.cz

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

/* See Copyright Notice in tnl/Copyright */

#pragma once

#include <fstream>

#include <istream>

#include <vector>

#include <map>

#include <unordered_map>

#include <type_traits>

#include <TNL/Containers/StaticVector.h>

#include <TNL/Meshes/MeshBuilder.h>

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

#include <TNL/Meshes/VTKTraits.h>

#ifdef HAVE_VTK

#include <vtkSmartPointer.h>

#include <vtkCell.h>

#include <vtkUnstructuredGrid.h>

#include <vtkUnstructuredGridReader.h>

#endif

namespace TNL {

namespace Meshes {

namespace Readers {

// types used in VTK

using VTKRealType = double;

#ifdef HAVE_VTK

using VTKIndexType = vtkIdType;

#else

using VTKIndexType = int;

#endif

// wrapper type for physical coordinates to preserve my sanity

using VTKPointType = Containers::StaticVector< 3, VTKRealType >;

template< typename Index = int >

class VTKReader_libvtk

{

public:

   using IndexType = Index;

   bool

   detectMesh( const String& fileName )

   {

      resetAll();

      if( ! loadVTKFile( fileName ) ) {

         resetAll();

         return false;

      }

      reset();

      return true;

   }

   template< typename MeshType >

   bool

   readMesh( const String& fileName, MeshType& mesh )

   {

      static_assert( std::is_same< IndexType, typename MeshType::GlobalIndexType >::value, "VTKReader_libvtk::IndexType and MeshType::GlobalIndexType must be the same type." );

      resetAll();

      if( ! loadVTKFile( fileName ) ) {

         resetAll();

         return false;

      }

      // GOTCHA: The unary "+" is a workaround due to odr-use of undefined static data member. See:

      // https://stackoverflow.com/questions/39646958/constexpr-static-member-before-after-c17

      // https://stackoverflow.com/questions/272900/undefined-reference-to-static-class-member/272996#272996

      TNL_ASSERT_EQ( this->worldDimension, + MeshType::Config::worldDimension, "world dimensions do not match" );

      TNL_ASSERT_EQ( this->meshDimension, + MeshType::Config::meshDimension, "mesh dimensions do not match" );

      const int subvertices = Topologies::Subtopology< typename MeshType::Config::CellTopology, 0 >::count;

      TNL_ASSERT_EQ( this->verticesInEntities.at( this->meshDimension ), subvertices, "numbers of cell subvertices do not match" );

      using MeshBuilder = MeshBuilder< MeshType >;

      using GlobalIndexType = typename MeshType::GlobalIndexType;

      const GlobalIndexType numberOfPoints = this->pointsData.size();

      const GlobalIndexType numberOfCells = this->entityIdMappings.at( this->meshDimension ).size();

      MeshBuilder meshBuilder;

      meshBuilder.setPointsCount( numberOfPoints );

      meshBuilder.setCellsCount( numberOfCells );

      for( GlobalIndexType i = 0; i < numberOfPoints; i++ ) {

         typename MeshType::PointType p;

         for( int j = 0; j < p.getSize(); j++ )

            p[ j ] = this->pointsData.at( i )[ j ];

         meshBuilder.setPoint( i, p );

      }

      const auto& cellIdMap = this->entityIdMappings.at( this->meshDimension );

      for( GlobalIndexType i = 0; i < numberOfCells; i++ ) {

         const VTKIndexType vtkCellIndex = cellIdMap.at( i );

         const auto& vtkCellSeeds = this->entitySeeds.at( vtkCellIndex );

         using CellSeedType = typename MeshBuilder::CellSeedType;

         TNL_ASSERT_EQ( CellSeedType::getCornersCount(), vtkCellSeeds.size(), "wrong number of subvertices" );

         CellSeedType& seed = meshBuilder.getCellSeed( i );

         for( int v = 0; v < CellSeedType::getCornersCount(); v++ ) {

            seed.setCornerId( v, vtkCellSeeds[ v ] );

         }

      }

      // drop all data from the VTK file since it's not needed anymore

      this->resetAll();

      return meshBuilder.build( mesh );

   }

   String

   getMeshType() const

   {

      // we can read only the UNSTRUCTURED_GRID dataset

      return "Meshes::Mesh";

   }

   int

   getWorldDimension() const

   {

      return this->worldDimension;

   }

   int

   getMeshDimension() const

   {

      return this->meshDimension;

   }

   VTK::EntityShape

   getCellShape() const

   {

      return this->entityTypes.at( this->meshDimension );

   }

//   int

//   getVerticesInCell() const

//   {

//      return this->verticesInEntities.at( this->getMeshDimension() );

//   }

//   VTK::EntityShape

//   getEntityType( int entityDimension ) const

//   {

//      return this->entityTypes.at( entityDimension );

//   }

   String

   getRealType() const

   {

      // TODO: how to extract it from the VTK object?

//      return "float";

      return "double";

   }

   String

   getGlobalIndexType() const

   {

      // not stored in the VTK file

      return "int";

   }

   String

   getLocalIndexType() const

   {

      // not stored in the VTK file

      return "short int";

   }

protected:

   int worldDimension = 0;

   int meshDimension = 0;

   // maps vertex indices to physical coordinates

   std::unordered_map< VTKIndexType, VTKPointType > pointsData;

   // maps entity dimension to the number of vertices in the entity

   std::map< int, int > verticesInEntities;

   // type for mapping TNL entity indices to VTK cell indices

   using TNL2VTKindexmap = std::unordered_map< IndexType, VTKIndexType >;

   // maps dimension to maps of TNL entity IDs to corresponding VTK cell IDs

   std::unordered_map< int, TNL2VTKindexmap > entityIdMappings;

   // maps VTK cell indices to entity seeds (set of indices of subvertices)

   std::unordered_map< VTKIndexType, std::vector< VTKIndexType > > entitySeeds;

   // maps dimension to VTK type of the entity with given dimension

   std::unordered_map< int, VTK::EntityShape > entityTypes;

   void reset()

   {

      pointsData.clear();

      verticesInEntities.clear();

      entityIdMappings.clear();

      entitySeeds.clear();

   }

   void resetAll()

   {

      reset();

      worldDimension = 0;

      meshDimension = 0;

      entityTypes.clear();

   }

   bool loadVTKFile( const String& fileName )

   {

#ifdef HAVE_VTK

      vtkSmartPointer< vtkUnstructuredGridReader > vtkReader

         = vtkSmartPointer<vtkUnstructuredGridReader>::New();

      vtkReader->SetFileName( fileName.getString() );

      vtkReader->Update();

      if( ! vtkReader->IsFileUnstructuredGrid() ) {

         std::cerr << "The file '" << fileName << "' is not a VTK Legacy file with an UNSTRUCTURED_GRID dataset type." << std::endl;

         return false;

      }

      auto& vtkMesh = *vtkReader->GetOutput();

      /* To determine the mesh type, we need:

       * - world dimension (default 1D, if some 2nd coordinate is non-zero -> 2D, some 3rd coordinate non-zero -> 3D)

       * - mesh dimension (highest dimension of the entities present in mesh)

       * - check that it is homogeneous - all entities of the same dimension should have the same number of vertices and type

       * - types of entities for each dimension

       *

       * To initialize the mesh, we need:

       * - points data (vertex indices and world coordinates)

       * - cell seeds (indices of vertices composing each cell)

       */

      const VTKIndexType numberOfPoints = vtkMesh.GetNumberOfPoints();

      const VTKIndexType numberOfCells = vtkMesh.GetNumberOfCells();

      if( numberOfPoints == 0 ) {

         std::cerr << "There are no points data in the file '" << fileName << "'." << std::endl;

         return false;

      }

      if( numberOfCells == 0 ) {

         std::cerr << "There are no cells data in the file '" << fileName << "'." << std::endl;

         return false;

      }

      for( VTKIndexType i = 0; i < numberOfPoints; i++ ) {

         VTKRealType* p = vtkMesh.GetPoint( i );

         // get world dimension

         for( int j = 0; j < 3; j++ )

            if( p[ j ] != 0.0 )

               this->worldDimension = std::max( this->worldDimension, j + 1 );

         // copy points data

         this->pointsData[ i ] = VTKPointType( p[ 0 ], p[ 1 ], p[ 2 ] );

      }

      for( VTKIndexType i = 0; i < numberOfCells; i++ ) {

         vtkCell* cell = vtkMesh.GetCell( i );

         const int dimension = cell->GetCellDimension();

         const int points = cell->GetNumberOfPoints();

         const VTK::EntityShape type = (VTK::EntityShape) cell->GetCellType();

         // number of vertices in entities

         if( this->verticesInEntities.find( dimension ) == this->verticesInEntities.cend() )

            this->verticesInEntities[ dimension ] = points;

         else if( this->verticesInEntities[ dimension ] != points ) {

            std::cerr << "Mixed unstructured meshes are not supported. There are elements of dimension " << dimension

                      << " with " << this->verticesInEntities[ dimension ] << " vertices and with " << points

                      << " vertices. The number of vertices per entity must be constant." << std::endl;

            this->reset();

            return false;

         }

         // entity types

         if( this->entityTypes.find( dimension ) == this->entityTypes.cend() )

            this->entityTypes.emplace( std::make_pair( dimension, type ) );

         else if( this->entityTypes[ dimension ] != type ) {

            std::cerr << "Mixed unstructured meshes are not supported. There are elements of dimension " << dimension

                      << " with type " << VTK::getShapeName( this->entityTypes[ dimension ] ) << " and " << VTK::getShapeName( type )

                      << ". The type of all entities with the same dimension must be the same." << std::endl;

            this->reset();

            return false;

         }

         // mapping between TNL and VTK indices

         auto& map = this->entityIdMappings[ dimension ];

         map[ map.size() ] = i;

         // copy seed

         auto& seed = this->entitySeeds[ i ];

         for( int j = 0; j < points; j++ )

            seed.push_back( cell->GetPointId( j ) );

      }

      // std::map is sorted, so the last key is the maximum

      this->meshDimension = this->verticesInEntities.rbegin()->first;

      if( this->meshDimension > this->worldDimension ) {

         std::cerr << "Invalid mesh: world dimension is " << this->worldDimension

                   << ", but mesh dimension is " << this->meshDimension << "." << std::endl;

         this->reset();

         return false;

      }

      return true;

#else

      std::cerr << "The VTKReader_libvtk needs to be compiled with the VTK library." << std::endl;

      return false;

#endif

   }

};

} // namespace Readers

} // namespace Meshes

} // namespace TNL

#TARGET_COMPILE_OPTIONS( MeshEntityTest PRIVATE ${CXX_TESTS_FLAGS} )

#TARGET_LINK_LIBRARIES( MeshEntityTest ${GTEST_BOTH_LIBRARIES} )

#ADD_TEST( MeshEntityTest ${EXECUTABLE_OUTPUT_PATH}/MeshEntityTest${CMAKE_EXECUTABLE_SUFFIX} )

##

## Tests with VTK

##

#find_package( VTK )

#if( VTK_FOUND )

#   include(${VTK_USE_FILE})

#

#   AddCompilerFlag( "-DHAVE_VTK " )

#   SET( VTK_COMMON_LIBRARIES vtkCommonCore ; vtkIOLegacy )

#endif( VTK_FOUND )

## MeshReaderTest is not a unit test so we disable it, because it takes

## a long time to compile.

##

## Mesh cannot be compiled by nvcc < 9 due to bugs in the compiler

#if( ${BUILD_CUDA} AND ${CUDA_VERSION_MAJOR} GREATER_EQUAL 9 )

#   CUDA_ADD_EXECUTABLE( MeshReaderTest MeshReaderTest.cu

#                        OPTIONS ${CXX_TESTS_FLAGS} )

#   TARGET_LINK_LIBRARIES( MeshReaderTest

#                           ${GTEST_BOTH_LIBRARIES}

#                           ${VTK_COMMON_LIBRARIES} )

#else()

#   ADD_EXECUTABLE( MeshReaderTest MeshReaderTest.cpp )

#   TARGET_COMPILE_OPTIONS( MeshReaderTest PRIVATE ${CXX_TESTS_FLAGS} )

#   TARGET_LINK_LIBRARIES( MeshReaderTest

#                           ${GTEST_BOTH_LIBRARIES}

#                           ${VTK_COMMON_LIBRARIES} )

#endif()

#include "MeshReaderTest.h"

#include "MeshReaderTest.h"

#include <TNL/Meshes/Mesh.h>

#include <TNL/Meshes/DefaultConfig.h>

#include <TNL/Meshes/BuildConfigTags.h>

#include <TNL/Meshes/TypeResolver/TypeResolver.h>

#include <TNL/Communicators/NoDistrCommunicator.h>

#ifdef HAVE_VTK

#include <TNL/Meshes/Readers/VTKReader_libvtk.h>

#endif

#include <TNL/Debugging/MemoryUsage.h>

#include "MeshTest.h"

// TODO: remove this after refactoring with clang-rename

#include "MeshEntityTest.h"

#include "BoundaryTagsTest.h"

using namespace TNL;

using namespace TNL::Meshes;

template< typename Cell,

          int WorldDimension = Cell::dimension,

          typename Real = double,

          typename GlobalIndex = int,

          typename LocalIndex = GlobalIndex,

          typename Id = void >

struct MyMeshConfig

   : public DefaultConfig< Cell, WorldDimension, Real, GlobalIndex, LocalIndex, Id >

{

   static constexpr bool entityStorage( int dimension )

   {

      return true;

//      return dimension == 0 || dimension == Cell::dimension;

   }

   template< typename EntityTopology >

   static constexpr bool subentityStorage( EntityTopology, int SubentityDimension )

   {

//      return entityStorage( EntityTopology::dimension );

      return entityStorage( EntityTopology::dimension ) &&

             SubentityDimension == 0;

   }

   template< typename EntityTopology >

   static constexpr bool subentityOrientationStorage( EntityTopology, int SubentityDimension ) 

   {

      return false;

   }

   template< typename EntityTopology >

   static constexpr bool superentityStorage( EntityTopology, int SuperentityDimension )

   {

//      return entityStorage( EntityTopology::dimension );

      return entityStorage( EntityTopology::dimension ) &&

             SuperentityDimension == Cell::dimension;

   }

   template< typename EntityTopology >

   static constexpr bool boundaryTagsStorage( EntityTopology )

   {

      using BaseType = DefaultConfig< Cell, WorldDimension, Real, GlobalIndex, LocalIndex, Id >;

      using FaceTopology = typename Topologies::Subtopology< Cell, BaseType::meshDimension - 1 >::Topology;

      return entityStorage( BaseType::meshDimension - 1 ) &&

             superentityStorage( FaceTopology(), BaseType::meshDimension ) &&

             ( EntityTopology::dimension >= BaseType::meshDimension - 1 || subentityStorage( FaceTopology(), EntityTopology::dimension ) );

      //return false;

   }

};

// specialization of BuildConfigTags

struct MyBuildConfigTag {};

namespace TNL {

namespace Meshes {

namespace BuildConfigTags {

// disable grids

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

struct GridTag< MyBuildConfigTag, Grid< Dimension, Real, Device, Index > >

{ enum { enabled = false }; };

// enable all cell topologies

template<> struct MeshCellTopologyTag< MyBuildConfigTag, Topologies::Edge > { enum { enabled = true }; };

template<> struct MeshCellTopologyTag< MyBuildConfigTag, Topologies::Triangle > { enum { enabled = true }; };

template<> struct MeshCellTopologyTag< MyBuildConfigTag, Topologies::Quadrilateral > { enum { enabled = true }; };

template<> struct MeshCellTopologyTag< MyBuildConfigTag, Topologies::Tetrahedron > { enum { enabled = true }; };

template<> struct MeshCellTopologyTag< MyBuildConfigTag, Topologies::Hexahedron > { enum { enabled = true }; };

template< typename CellTopology, int WorldDimension >

struct MeshWorldDimensionTag< MyBuildConfigTag, CellTopology, WorldDimension >

{ enum { enabled = ( WorldDimension == CellTopology::dimension ) }; };

template< typename Real > struct MeshRealTag< MyBuildConfigTag, Real > { enum { enabled = false }; };

template<> struct MeshRealTag< MyBuildConfigTag, float > { enum { enabled = true }; };

template<> struct MeshRealTag< MyBuildConfigTag, double > { enum { enabled = true }; };

template< typename GlobalIndex > struct MeshGlobalIndexTag< MyBuildConfigTag, GlobalIndex > { enum { enabled = false }; };

template<> struct MeshGlobalIndexTag< MyBuildConfigTag, int > { enum { enabled = true }; };

template< typename LocalIndex > struct MeshLocalIndexTag< MyBuildConfigTag, LocalIndex > { enum { enabled = false }; };

template<> struct MeshLocalIndexTag< MyBuildConfigTag, short int > { enum { enabled = true }; };

template<>

struct MeshConfigTemplateTag< MyBuildConfigTag >

{

   template< typename Cell, int WorldDimension, typename Real, typename GlobalIndex, typename LocalIndex, typename Id >

   using MeshConfig = MyMeshConfig< Cell, WorldDimension, Real, GlobalIndex, LocalIndex, Id >;

};

} // namespace BuildConfigTags

} // namespace Meshes

} // namespace TNL

template< typename MeshType >

class MeshTester

{

public:

   static bool run( const String& fileName )

   {

      std::cout << "pre-init\t";

      Debugging::printMemoryUsage();

#ifdef HAVE_VTK

      MeshType mesh_libvtk;

      Readers::VTKReader_libvtk<> reader;

      if( ! reader.readMesh( fileName, mesh_libvtk ) )

         return false;

      std::cout << "libvtk vertices: " << mesh_libvtk.template getEntitiesCount< 0 >() << std::endl;

      std::cout << "libvtk faces: " << mesh_libvtk.template getEntitiesCount< MeshType::getMeshDimension() - 1 >() << std::endl;

      std::cout << "libvtk cells: " << mesh_libvtk.template getEntitiesCount< MeshType::getMeshDimension() >() << std::endl;

#endif

      Timer timer;

      timer.start();

      MeshType mesh;

      Meshes::DistributedMeshes::DistributedMesh<MeshType> distributedMesh;

      if( ! loadMesh<Communicators::NoDistrCommunicator>( fileName, mesh, distributedMesh ) )

         return false;

      timer.stop();

      std::cout << "Loading took " << timer.getRealTime() << " seconds." << std::endl;

      std::cout << "vertices: " << mesh.template getEntitiesCount< 0 >() << std::endl;

      std::cout << "faces: " << mesh.template getEntitiesCount< MeshType::getMeshDimension() - 1 >() << std::endl;

      std::cout << "cells: " << mesh.template getEntitiesCount< MeshType::getMeshDimension() >() << std::endl;

      std::cout << "post-init\t";

      Debugging::printMemoryUsage();

#ifdef HAVE_VTK

      std::cout << "mesh_libvtk == mesh: " << std::boolalpha << (mesh_libvtk == mesh) << std::endl;

      if( mesh_libvtk != mesh ) {

         std::cerr << "mesh_libvtk:\n" << mesh_libvtk << "\n\nmesh:\n" << mesh << std::endl;

         return false;

      }

#endif

#ifdef HAVE_GTEST 

      std::cout << "Running basic I/O tests..." << std::endl;

      MeshTest::testFinishedMesh( mesh );

#endif

//      mesh.save( "mesh-test.tnl" );

      return true;

   }

};

int main( int argc, char* argv[] )

{

   if( argc < 2 ) {

      std::cerr << "Usage: " << argv[ 0 ] << " filename.[tnl|ng|vtk] ..." << std::endl;

      return EXIT_FAILURE;

   }

   bool result = true;

   for( int i = 1; i < argc; i++ ) {

      String fileName = argv[ i ];

      result &= resolveMeshType< MyBuildConfigTag, Devices::Host, MeshTester >

                  ( fileName,

                    fileName  // passed to MeshTester::run

                  );

   }

   std::cout << "final\t";

   Debugging::printMemoryUsage();

   return ! result;

}