Add reader for raw files

src/geometry/structured/BinaryRawFile.cpp

+//======================================================================================================================
+//
+//  This file is part of waLBerla. waLBerla is free software: you can 
+//  redistribute it and/or modify it under the terms of the GNU General Public
+//  License as published by the Free Software Foundation, either version 3 of 
+//  the License, or (at your option) any later version.
+//  
+//  waLBerla is distributed in the hope that it will be useful, but WITHOUT 
+//  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or 
+//  FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License 
+//  for more details.
+//  
+//  You should have received a copy of the GNU General Public License along
+//  with waLBerla (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
+//
+//! \file BinaryRawFileReader.cpp
+//! \ingroup geometry
+//! \author Christian Godenschwager <christian.godenschwager@fau.de>
+//
+//======================================================================================================================
+
+#include "BinaryRawFile.h"
+
+
+namespace walberla {
+namespace geometry {
+
+BinaryRawFile::BinaryRawFile( const std::string & filename, const Vector3< uint_t > & size, const std::string & datatype )
+   : size_( size )
+{
+   init( filename, datatype );
+}
+
+BinaryRawFile::BinaryRawFile( const Config::BlockHandle & configBlock )
+{
+   size_ = configBlock.getParameter< Vector3< uint_t > >( "size" );
+   const std::string filename = configBlock.getParameter< std::string >( "filename" );
+   const std::string datatype = configBlock.getParameter< std::string >( "datatype" );
+
+   init( filename, datatype );
+}
+
+void BinaryRawFile::init( const std::string & filename, const std::string & datatype )
+{
+   if (datatype == "uint8")
+      init<uint8_t>( filename );
+   else if (datatype == "int8")
+      init<int8_t>( filename );
+   else if (datatype == "uint16")
+      init<uint16_t>( filename );
+   else if (datatype == "int16")
+      init<int16_t>( filename );
+   else if (datatype == "uint32")
+      init<uint32_t>( filename );
+   else if (datatype == "int32")
+      init<int32_t>( filename );
+   else if (datatype == "uint64")
+      init<uint64_t>( filename );
+   else if (datatype == "int64")
+      init<int64_t>( filename );
+   else if (datatype == "float")
+      init<float>( filename );
+   else if (datatype == "double")
+      init<double>( filename );
+   else
+      WALBERLA_ABORT( "Unknown datatype!" );
+}
+
+
+} // namespace geometry
+} // namespace walberla
\ No newline at end of file

src/geometry/structured/BinaryRawFile.h

+//======================================================================================================================
+//
+//  This file is part of waLBerla. waLBerla is free software: you can 
+//  redistribute it and/or modify it under the terms of the GNU General Public
+//  License as published by the Free Software Foundation, either version 3 of 
+//  the License, or (at your option) any later version.
+//  
+//  waLBerla is distributed in the hope that it will be useful, but WITHOUT 
+//  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or 
+//  FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License 
+//  for more details.
+//  
+//  You should have received a copy of the GNU General Public License along
+//  with waLBerla (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
+//
+//! \file BinaryRawFileReader.h
+//! \ingroup geometry
+//! \author Christian Godenschwager <christian.godenschwager@fau.de>
+//
+//======================================================================================================================
+
+#pragma once
+
+#include "core/DataTypes.h"
+#include "core/mpi/Broadcast.h"
+#include "core/mpi/MPIManager.h"
+#include "core/mpi/BufferDataTypeExtensions.h"
+#include "core/math/AABB.h"
+#include "core/math/Vector3.h"
+#include "core/config/Config.h"
+
+#include <fstream>
+#include <string>
+#include <vector>
+#include <iterator>
+
+
+namespace walberla {
+namespace geometry {
+
+class BinaryRawFile
+{
+public:
+   BinaryRawFile( const std::string & filename, const Vector3< uint_t > & size, const std::string & datatype );
+
+   template< typename T >
+   static BinaryRawFile loadFile( const std::string & filename, const Vector3< uint_t > & size );
+
+   BinaryRawFile( const Config::BlockHandle & configBlock );
+
+   inline bool get( const Vector3< uint_t > & pos ) const;
+   inline bool get( const uint_t x, const uint_t y, const uint_t z ) const;
+
+   const Vector3< uint_t> & size() const { return size_; }
+
+private:
+   template< typename T >
+   BinaryRawFile( const std::string & filename, const Vector3< uint_t > & size, const T );
+
+   void init( const std::string & filename, const std::string & datatype );
+
+   template< typename T >
+   void init( const std::string & filename );
+
+   Vector3< uint_t > size_;
+   std::vector< bool > data_;
+};
+
+class BinaryRawFileInterpolator
+{
+public:
+   enum Interpolator { NEAREST_NEIGHBOR };
+
+   BinaryRawFileInterpolator( const AABB & aabb, const BinaryRawFile & binaryRawFile, Interpolator interpolator )
+      : aabb_(aabb), interpolator_(interpolator), binaryRawFile_( binaryRawFile ) {}
+
+   inline bool get( const Vector3< real_t > & pos ) const;
+   inline bool get( const real_t x, const real_t y, const real_t z ) const;
+
+   const AABB & aabb() const { return aabb_; }
+
+   Interpolator interpolator() const { return interpolator_; }
+
+private:
+   inline bool getNearestNeighbor( const real_t x, const real_t y, const real_t z ) const;
+
+   AABB aabb_;
+   Interpolator interpolator_;
+   const BinaryRawFile & binaryRawFile_;
+};
+
+
+template<typename T>
+BinaryRawFile::BinaryRawFile( const std::string & filename, const Vector3<uint_t> & size, const T )
+   : size_(size)
+{
+   init<T>( filename );
+}
+
+
+template<typename T>
+void BinaryRawFile::init( const std::string & filename )
+{
+   const uint_t numElements = size_[0] * size_[1] * size_[2];
+   data_.reserve( numElements );
+   WALBERLA_ROOT_SECTION()
+   {
+      std::ifstream ifs( filename, std::ifstream::in | std::ifstream::binary );
+      std::transform( std::istream_iterator<T>( ifs ), std::istream_iterator<T>(),
+                      std::back_inserter( data_ ), []( const T v ) { return v > T( 0 ); } );
+      WALBERLA_CHECK_EQUAL( numElements, data_.size(), "Error reading file \"" << filename << "\"!" );
+   }
+   mpi::broadcastObject( data_ );
+}
+
+template< typename T >
+static BinaryRawFile loadFile( const std::string & filename, const Vector3< uint_t > & size )
+{
+   return BinaryRawFileReader( filename, size, T() );
+}
+
+bool BinaryRawFile::get( const Vector3< uint_t > & pos ) const
+{
+   return get( pos[0], pos[1], pos[2] );
+}
+
+inline bool BinaryRawFile::get( const uint_t x, const uint_t y, const uint_t z ) const
+{
+   WALBERLA_ASSERT_LESS( x, size_[0] );
+   WALBERLA_ASSERT_LESS( y, size_[1] );
+   WALBERLA_ASSERT_LESS( z, size_[2] );
+
+   const uint_t i = z * size_[0] * size_[1] + y * size_[0] + x;
+   
+   WALBERLA_ASSERT_LESS( i, data_.size() );
+
+   return data_[i];
+}
+
+
+bool BinaryRawFileInterpolator::get( const Vector3< real_t > & pos ) const
+{
+   return get( pos[0], pos[1], pos[2] );
+}
+
+
+bool BinaryRawFileInterpolator::get( const real_t x, const real_t y, const real_t z ) const
+{
+   WALBERLA_ASSERT( aabb_.contains( x, y, z ) );
+
+   switch (interpolator_)
+   {
+   case NEAREST_NEIGHBOR:
+      return getNearestNeighbor( x, y, z );
+   default:
+      WALBERLA_ABORT( "Unknown Interpolator!" );
+   }
+}
+
+bool BinaryRawFileInterpolator::getNearestNeighbor( const real_t x, const real_t y, const real_t z ) const
+{
+   uint_t xInt = uint_c( (x - aabb_.xMin()) / aabb_.xSize() * real_t( binaryRawFile_.size()[0] ) + real_t( 0.5 ) );
+   uint_t yInt = uint_c( (y - aabb_.yMin()) / aabb_.ySize() * real_t( binaryRawFile_.size()[1] ) + real_t( 0.5 ) );
+   uint_t zInt = uint_c( (z - aabb_.zMin()) / aabb_.zSize() * real_t( binaryRawFile_.size()[2] ) + real_t( 0.5 ) );
+
+   return binaryRawFile_.get( xInt, yInt, zInt );
+}
+
+} // namespace geometry
+} // namespace walberla
\ No newline at end of file

tests/geometry/BinaryRawFileTest.cpp

+//======================================================================================================================
+//
+//  This file is part of waLBerla. waLBerla is free software: you can 
+//  redistribute it and/or modify it under the terms of the GNU General Public
+//  License as published by the Free Software Foundation, either version 3 of 
+//  the License, or (at your option) any later version.
+//  
+//  waLBerla is distributed in the hope that it will be useful, but WITHOUT 
+//  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or 
+//  FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License 
+//  for more details.
+//  
+//  You should have received a copy of the GNU General Public License along
+//  with waLBerla (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
+//
+//! \file BinaryRawFileTest.cpp
+//! \ingroup geometry
+//! \author Christian Godenschwager <christian.godenschwager@fau.de>
+//
+//======================================================================================================================
+
+#include "blockforest/Initialization.h"
+
+#include "core/DataTypes.h"
+#include "core/Abort.h"
+#include "core/logging/Logging.h"
+#include "core/mpi/Environment.h"
+
+#include "field/AddToStorage.h"
+#include "field/vtk/VTKWriter.h"
+
+#include "geometry/structured/BinaryRawFile.h"
+
+#include "vtk/Initialization.h"
+
+#include <string>
+#include <vector>
+
+namespace walberla {
+namespace geometry {
+
+void test( const std::string & filename, const Vector3<uint_t> & size, const std::string & datatype, const bool vtkOut )
+{
+   WALBERLA_LOG_INFO( "Testing unscaled" );
+   BinaryRawFile brf( filename, size, datatype );
+
+   auto blocks = blockforest::createUniformBlockGrid( uint_t(1), uint_t( 1 ), uint_t( 1 ), 
+      size[0], size[1], size[2], 
+      real_t( 1 ), 
+      uint_t( 1 ), uint_t( 1 ), uint_t( 1 ) );
+
+   typedef GhostLayerField< uint8_t, uint_t( 1 ) > ScalarField;
+
+   BlockDataID scalarFieldID = field::addToStorage<ScalarField>( blocks, "BinaryRawFile" );
+
+   for ( auto & block : *blocks)
+   {
+      auto field = block.getData<ScalarField>( scalarFieldID );
+      
+      CellInterval ci( 0, 0, 0, cell_idx_c( size[0] ) - 1, cell_idx_c( size[1] ) - 1, cell_idx_c( size[2] ) - 1 );
+
+      for (const Cell & c : ci)
+      {
+         field->get( c ) = brf.get( uint_c( c[0] ), uint_c( c[1] ), uint_c( c[2] ) );
+      }
+   }
+
+   if (vtkOut)
+   {
+      WALBERLA_LOG_INFO( "Writing unscaled" );
+      auto vtkOutput = vtk::createVTKOutput_BlockData( blocks, "BinaryRawFile" );
+      vtkOutput->addCellDataWriter( make_shared< field::VTKWriter< ScalarField, uint8_t > >( scalarFieldID, "BinaryRawFile" ) );
+      writeFiles( vtkOutput, true )();
+   }
+}
+
+void testScaled( const std::string & filename, const Vector3<uint_t> & size, const std::string & datatype, const bool vtkOut )
+{
+   WALBERLA_LOG_INFO( "Testing scaled" );
+   BinaryRawFile brf( filename, size, datatype );
+
+   Vector3<uint_t> scaledSize( std::max( uint_t( 1 ), size[0] / uint_t( 2 ) ),
+                               std::max( uint_t( 1 ), size[1] / uint_t( 3 ) ),
+                               std::max( uint_t( 1 ), size[2] / uint_t( 5 ) ) );
+
+   auto blocks = blockforest::createUniformBlockGrid( uint_t( 1 ), uint_t( 1 ), uint_t( 1 ),
+      scaledSize[0], scaledSize[1], scaledSize[2],
+      real_t( 1 ),
+      uint_t( 1 ), uint_t( 1 ), uint_t( 1 ) );
+
+   BinaryRawFileInterpolator brfi( blocks->getDomain(), brf, BinaryRawFileInterpolator::NEAREST_NEIGHBOR );
+
+   typedef GhostLayerField< uint8_t, uint_t( 1 ) > ScalarField;
+
+   BlockDataID scalarFieldID = field::addToStorage<ScalarField>( blocks, "BinaryRawFile" );
+
+   for (auto & block : *blocks)
+   {
+      auto field = block.getData<ScalarField>( scalarFieldID );
+
+      CellInterval ci( 0, 0, 0, cell_idx_c( scaledSize[0] ) - 1, cell_idx_c( scaledSize[1] ) - 1, cell_idx_c( scaledSize[2] ) - 1 );
+
+      for (const Cell & c : ci)
+      {
+         auto pos = blocks->getBlockLocalCellCenter( block, c );
+         field->get( c ) = brfi.get( pos );
+      }
+   }
+
+   if (vtkOut)
+   {
+      WALBERLA_LOG_INFO( "Writing scaled" );
+      auto vtkOutput = vtk::createVTKOutput_BlockData( blocks, "BinaryRawFileScaled" );
+      vtkOutput->addCellDataWriter( make_shared< field::VTKWriter< ScalarField, uint8_t > >( scalarFieldID, "BinaryRawFile" ) );
+      writeFiles( vtkOutput, true )();
+   }
+}
+
+} // namespace geometry
+} // namespace walberla
+
+int main( int argc, char * argv[] )
+{
+   walberla::mpi::Environment env( argc, argv );
+
+   std::vector< std::string > args( argv, argv + argc );
+
+   auto it = std::find( args.begin(), args.end(), std::string( "--vtk" ) );
+   const bool vtk = it != args.end();
+   if (it != args.end())
+      args.erase( it );
+
+   if (args.size() != 6u)
+      WALBERLA_ABORT_NO_DEBUG_INFO( "USAGE: " << args.front() << " [--vtk] FILENAME XSIZE YSIZE ZSIZE DATATYPE" );
+
+   const std::string filename = args[1];
+
+   walberla::Vector3< walberla::uint_t > size;
+   try
+   {
+      size.set( std::stoull( args[2] ), 
+                std::stoull( args[3] ),
+                std::stoull( args[4] ) );
+   }
+   catch (const std::invalid_argument &)
+   {
+      WALBERLA_ABORT_NO_DEBUG_INFO( "USAGE: " << args.front() << " FILENAME XSIZE YSIZE ZSIZE DATATYPE" );
+   }
+
+   const std::string datatype = args[5];
+
+   walberla::geometry::test( filename, size, datatype, vtk );
+   walberla::geometry::testScaled( filename, size, datatype, vtk );
+}
+

tests/geometry/CMakeLists.txt

@@ -28,3 +28,5 @@ waLBerla_execute_test( NAME ScalarFieldFromGrayScaleImageTest )
 
 waLBerla_compile_test( FILES ImageTest.cpp )
 waLBerla_execute_test( NAME ImageTest )
+
+waLBerla_compile_test( FILES BinaryRawFileTest.cpp )