From e811513b095a2e558afeed5e2f9917a70e62b61e Mon Sep 17 00:00:00 2001 From: Ivan Pribec Date: Thu, 11 May 2023 22:56:19 +0200 Subject: [PATCH 01/16] Add Makefile variables --- Makefile | 20 ++++++++++++++++---- 1 file changed, 16 insertions(+), 4 deletions(-) diff --git a/Makefile b/Makefile index 22ff3f5..183c0bb 100644 --- a/Makefile +++ b/Makefile @@ -1,4 +1,16 @@ -hemeXtract: - g++ -Wall -O3 hemeXtract.cc -o hemeXtract -debug: - g++ -Wall -g hemeXtract.cc -o hemeXtract + +CXX=g++ -Wall +CXXFLAGS=-O3 +CXXFLAGS_DEBUG=-g + +SRC=hemeXtract.cc + +hemeXtract: $(SRC) + $(CXX) $(CXXFLAGS) -o $@ $< + +debug: $(SRC) + $(CXX) $(CXXFLAGS_DEBUG) -o $@ $< + +.phony: clean +clean: + rm -f hemeXtract \ No newline at end of file From 0db28b38ab429a8f50c7e184f7861e497dd3d043 Mon Sep 17 00:00:00 2001 From: Ivan Pribec Date: Tue, 23 May 2023 14:14:51 +0200 Subject: [PATCH 02/16] Add missing header from rpc for MacOS --- hemeXtract.cc | 17 ++++++++++------- 1 file changed, 10 insertions(+), 7 deletions(-) diff --git a/hemeXtract.cc b/hemeXtract.cc index 1fb2999..56c621d 100644 --- a/hemeXtract.cc +++ b/hemeXtract.cc @@ -1,12 +1,15 @@ -#include -#include -#include -#include -#include -#include -#include +#include +#include +#include +#include #include #include +#include +#include + +#include +#include +#include #include "Mapping.h" #include "Compare.h" From 5164c176e8e9870a0227257c3b4378413c94c73d Mon Sep 17 00:00:00 2001 From: Ivan Pribec Date: Tue, 23 May 2023 17:20:56 +0200 Subject: [PATCH 03/16] Add basic CMake build --- CMakeLists.txt | 16 +++++++++ cmake/Findargp.cmake | 78 ++++++++++++++++++++++++++++++++++++++++++++ hemeXtract.cc | 4 +-- 3 files changed, 96 insertions(+), 2 deletions(-) create mode 100644 CMakeLists.txt create mode 100644 cmake/Findargp.cmake diff --git a/CMakeLists.txt b/CMakeLists.txt new file mode 100644 index 0000000..8a03730 --- /dev/null +++ b/CMakeLists.txt @@ -0,0 +1,16 @@ +cmake_minimum_required(VERSION 3.15 FATAL_ERROR) + +project(hemeXtract + VERSION 0.1.0 + LANGUAGES CXX) + +set(CMAKE_CXX_STANDARD 11) + +set(CMAKE_MODULE_PATH "${CMAKE_SOURCE_DIR}/cmake") + +find_package(argp REQUIRED) + +add_executable(hemeXtract hemeXtract.cc) + +target_include_directories(hemeXtract PRIVATE ${argp_INCLUDE_DIRS}) +target_link_libraries(hemeXtract PRIVATE ${argp_LIBRARIES}) diff --git a/cmake/Findargp.cmake b/cmake/Findargp.cmake new file mode 100644 index 0000000..fae190c --- /dev/null +++ b/cmake/Findargp.cmake @@ -0,0 +1,78 @@ +#[=======================================================================[.rst: +Findargp +------- + +Finds the argp library. + +Result Variables +^^^^^^^^^^^^^^^^ + +This will define the following variables: + +``Foo_FOUND`` + True if the system has the Foo library. +``Foo_INCLUDE_DIRS`` + Include directories needed to use Foo. +``Foo_LIBRARIES`` + Libraries needed to link to Foo. + +Cache Variables +^^^^^^^^^^^^^^^ + +The following cache variables may also be set: + +``Foo_INCLUDE_DIR`` + The directory containing ``foo.h``. +``Foo_LIBRARY`` + The path to the Foo library. + +#]=======================================================================] + +include(FindPackageHandleStandardArgs) +include(CheckCXXSymbolExists) + +set(CMAKE_REQUIRED_QUITE TRUE) + +# Check if header file exists +find_path(argp_INCLUDE_DIR + NAMES argp.h + HINTS /usr/include /usr/local/include) + +set(CMAKE_REQUIRED_INCLUDES ${argp_INCLUDE_DIR}) + +if(argp_INCLUDE_DIR) + # Check for function argp_parse in default libraries (libc) + check_cxx_symbol_exists(argp_parse "argp.h" ARGP_LIBC) + if(ARGP_LIBC) + set(argp_LIBRARY c CACHE STRING + "argp library file") + else() + find_library(argp_LIBRARY "argp" + HINTS /usr/local/lib + DOC "argp library file") + if (argp_LIBRARY) + set(CMAKE_REQUIRED_LIBRARIES ${argp_LIBRARY}) + + # Check for function argp_parse + check_cxx_symbol_exists(argp_parse "argp.h" ARGP_EXTERNAL) + if(NOT ARGP_EXTERNAL) + message(FATAL_ERROR + "The argp library was found in ${argp_LIBRARY}," + "but it doesn't contain a symbol named argp_parse.") + endif() + endif() + endif() +endif() + +find_package_handle_standard_args(argp + DEFAULT_MSG argp_LIBRARY argp_INCLUDE_DIR) + +if(argp_FOUND) + set(argp_LIBRARIES ${argp_LIBRARY}) + set(argp_INCLUDE_DIRS ${argp_INCLUDE_DIR}) +endif() + +mark_as_advanced( + argp_LIBRARY + argp_INCLUDE_DIR +) \ No newline at end of file diff --git a/hemeXtract.cc b/hemeXtract.cc index 56c621d..0283c8f 100644 --- a/hemeXtract.cc +++ b/hemeXtract.cc @@ -1,11 +1,11 @@ #include #include #include +#include #include +#include #include #include -#include -#include #include #include From 1df1af6eb1bf9b2fb874397f6a3cc3e1bea23503 Mon Sep 17 00:00:00 2001 From: Ivan Pribec Date: Tue, 23 May 2023 17:22:27 +0200 Subject: [PATCH 04/16] Fix argp module docstring --- cmake/Findargp.cmake | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/cmake/Findargp.cmake b/cmake/Findargp.cmake index fae190c..f9a3cf6 100644 --- a/cmake/Findargp.cmake +++ b/cmake/Findargp.cmake @@ -9,22 +9,22 @@ Result Variables This will define the following variables: -``Foo_FOUND`` - True if the system has the Foo library. -``Foo_INCLUDE_DIRS`` - Include directories needed to use Foo. -``Foo_LIBRARIES`` - Libraries needed to link to Foo. +``argp_FOUND`` + True if the system has the argp library. +``argp_INCLUDE_DIRS`` + Include directories needed to use argp. +``argp_LIBRARIES`` + Libraries needed to link to argp. Cache Variables ^^^^^^^^^^^^^^^ The following cache variables may also be set: -``Foo_INCLUDE_DIR`` - The directory containing ``foo.h``. -``Foo_LIBRARY`` - The path to the Foo library. +``argp_INCLUDE_DIR`` + The directory containing ``argp.h``. +``argp_LIBRARY`` + The path to the argp library. #]=======================================================================] From 8069cda67c7fce2de4ae65e210c34303008d161f Mon Sep 17 00:00:00 2001 From: Ivan Pribec Date: Tue, 23 May 2023 18:36:12 +0200 Subject: [PATCH 05/16] Add build workflow --- .github/workflows/build.yml | 28 ++++++++++++++++++++++++++++ 1 file changed, 28 insertions(+) create mode 100644 .github/workflows/build.yml diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml new file mode 100644 index 0000000..f0aaba4 --- /dev/null +++ b/.github/workflows/build.yml @@ -0,0 +1,28 @@ +name: build +on: [push, pull_request] + +jobs: + build_debug: + runs_on: ${{ matrix.os }} + strategy: + fail-fast: false + matrix: + os: [macos-latest] + toolchain: + - {compiler: gcc, version: 12} + - {compiler: clang, version: 14} + env: + CXX: ${{ matrix.toolchain.compiler }} + steps: + - name: Checkout source code + uses: actions/checkout@v3 + + - name: Install argp-standalone + if: contains(matrix.os, 'macos') + run: | + brew install argp-standalone + + - name: Build hemeXtract + run: | + mkdir build && cd build + cmake .. \ No newline at end of file From bbb5b80b5935830b709a38421c344cf69ac75b7a Mon Sep 17 00:00:00 2001 From: Ivan Pribec Date: Tue, 23 May 2023 18:38:49 +0200 Subject: [PATCH 06/16] Fix workflow typo --- .github/workflows/build.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml index f0aaba4..a0a26a6 100644 --- a/.github/workflows/build.yml +++ b/.github/workflows/build.yml @@ -3,7 +3,7 @@ on: [push, pull_request] jobs: build_debug: - runs_on: ${{ matrix.os }} + runs-on: ${{ matrix.os }} strategy: fail-fast: false matrix: From 48ddf3b6eedd6d34b0727c4fa2a5093203aee338 Mon Sep 17 00:00:00 2001 From: Ivan Pribec Date: Tue, 23 May 2023 19:07:31 +0200 Subject: [PATCH 07/16] Adjust cxx variable --- .github/workflows/build.yml | 15 ++++++++++----- 1 file changed, 10 insertions(+), 5 deletions(-) diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml index a0a26a6..f822b0c 100644 --- a/.github/workflows/build.yml +++ b/.github/workflows/build.yml @@ -9,14 +9,18 @@ jobs: matrix: os: [macos-latest] toolchain: - - {compiler: gcc, version: 12} - - {compiler: clang, version: 14} + - {compiler: gcc, cxx: g++-12} + - {compiler: clang, cxx: clang++} env: - CXX: ${{ matrix.toolchain.compiler }} + CXX: ${{ matrix.toolchain.cxx }} + steps: + - name: Modify gcc symlink + run: ln -s /usr/local/bin/g + - name: Checkout source code uses: actions/checkout@v3 - + - name: Install argp-standalone if: contains(matrix.os, 'macos') run: | @@ -25,4 +29,5 @@ jobs: - name: Build hemeXtract run: | mkdir build && cd build - cmake .. \ No newline at end of file + cmake .. -DCMAKE_BUILD_TYPE=Debug + make hemeXtract \ No newline at end of file From 3fc4d9a715642e765e3e10ff10bf1c3b83f5ee36 Mon Sep 17 00:00:00 2001 From: Ivan Pribec Date: Tue, 23 May 2023 19:10:52 +0200 Subject: [PATCH 08/16] Add ubuntu-latest runner --- .github/workflows/build.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml index f822b0c..f71644e 100644 --- a/.github/workflows/build.yml +++ b/.github/workflows/build.yml @@ -7,7 +7,7 @@ jobs: strategy: fail-fast: false matrix: - os: [macos-latest] + os: [ubuntu-latest,macos-latest] toolchain: - {compiler: gcc, cxx: g++-12} - {compiler: clang, cxx: clang++} From 80a5f496d7f5dfba0f42464da6663553dd6ad493 Mon Sep 17 00:00:00 2001 From: Ivan Pribec Date: Tue, 23 May 2023 19:29:06 +0200 Subject: [PATCH 09/16] Add TIRPC build support --- .github/workflows/build.yml | 9 +++++---- CMakeLists.txt | 10 ++++++++++ cmake/FindTIRPC.cmake | 24 ++++++++++++++++++++++++ 3 files changed, 39 insertions(+), 4 deletions(-) create mode 100644 cmake/FindTIRPC.cmake diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml index f71644e..e3f75ad 100644 --- a/.github/workflows/build.yml +++ b/.github/workflows/build.yml @@ -15,17 +15,18 @@ jobs: CXX: ${{ matrix.toolchain.cxx }} steps: - - name: Modify gcc symlink - run: ln -s /usr/local/bin/g - - name: Checkout source code uses: actions/checkout@v3 - - name: Install argp-standalone + - name: Install argp-standalone (MacOS) if: contains(matrix.os, 'macos') run: | brew install argp-standalone + - name: Install libtirpc (Ubuntu) + if: contains(matrix.os,'ubuntu') + run: sudo apt-get install libtirpc + - name: Build hemeXtract run: | mkdir build && cd build diff --git a/CMakeLists.txt b/CMakeLists.txt index 8a03730..721d621 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -10,7 +10,17 @@ set(CMAKE_MODULE_PATH "${CMAKE_SOURCE_DIR}/cmake") find_package(argp REQUIRED) +if(CMAKE_SYSTEM_NAME STREQUAL "Linux") + find_package(RPC REQUIRED) +endif() + + add_executable(hemeXtract hemeXtract.cc) target_include_directories(hemeXtract PRIVATE ${argp_INCLUDE_DIRS}) target_link_libraries(hemeXtract PRIVATE ${argp_LIBRARIES}) + +if(TIRPC_FOUND) + target_include_directories(hemeXtract PRIVATE ${TIRPC_INCLUDE_DIRS}) + target_link_libraries(hemeXtract PRIVATE ${TIRPC_LIBRARY}) +endif() diff --git a/cmake/FindTIRPC.cmake b/cmake/FindTIRPC.cmake new file mode 100644 index 0000000..c5b0037 --- /dev/null +++ b/cmake/FindTIRPC.cmake @@ -0,0 +1,24 @@ +find_package(PkgConfig QUIET) +pkg_check_modules(PC_TIRPC libtirpc) + +# check for header file +find_path(TIRPC_INCLUDE_DIRS + NAMES netconfig.h + HINTS ${TIRPC_DIR}/include $ENV{TIRPC_DIR}/include ${PC_TIRPC_INCLUDE_DIRS} + PATH_SUFFIXES tirpc + DOC "directory where the TIRPC header is located") + +# check for tirpc +find_library(TIRPC_LIBRARY + NAMES tirpc + HINTS ${TIRPC_DIR}/lib $ENV{TIRPC_DIR}/lib ${PC_TIRPC_LIBRARY_DIRS} + NO_DEFAULT_PATH + DOC "the TIRPC library") +find_library(TIRPC_LIBRARY + NAMES tirpc + DOC "the TIRPC library") + +include(FindPackageHandleStandardArgs) +find_package_handle_standard_args(TIRPC + REQUIRED_VARS TIRPC_INCLUDE_DIRS TIRPC_LIBRARY + VERSION_VAR TIRPC_VERSION) \ No newline at end of file From 92c5f443c09500f4e6033ef01efd6784593e6135 Mon Sep 17 00:00:00 2001 From: Ivan Pribec Date: Tue, 23 May 2023 19:38:35 +0200 Subject: [PATCH 10/16] Fix typo --- CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 721d621..0cdfec6 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -11,7 +11,7 @@ set(CMAKE_MODULE_PATH "${CMAKE_SOURCE_DIR}/cmake") find_package(argp REQUIRED) if(CMAKE_SYSTEM_NAME STREQUAL "Linux") - find_package(RPC REQUIRED) + find_package(TIRPC REQUIRED) endif() From 27c2881cd4c3a8d581c9d0539f347bda0321badb Mon Sep 17 00:00:00 2001 From: Ivan Pribec Date: Tue, 23 May 2023 19:40:46 +0200 Subject: [PATCH 11/16] Fix libtirpc package name --- .github/workflows/build.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml index e3f75ad..1d22b18 100644 --- a/.github/workflows/build.yml +++ b/.github/workflows/build.yml @@ -25,7 +25,7 @@ jobs: - name: Install libtirpc (Ubuntu) if: contains(matrix.os,'ubuntu') - run: sudo apt-get install libtirpc + run: sudo apt-get install libtirpc-dev - name: Build hemeXtract run: | From d20c3149e55b80143607b89b27ab5a9e9e8a5396 Mon Sep 17 00:00:00 2001 From: Ivan Pribec Date: Tue, 23 May 2023 20:08:30 +0200 Subject: [PATCH 12/16] Export TIRPC_DIR environment variable --- .github/workflows/build.yml | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml index 1d22b18..d3bc4d4 100644 --- a/.github/workflows/build.yml +++ b/.github/workflows/build.yml @@ -25,7 +25,9 @@ jobs: - name: Install libtirpc (Ubuntu) if: contains(matrix.os,'ubuntu') - run: sudo apt-get install libtirpc-dev + run: | + sudo apt-get install libtirpc-dev + export TIRPC_DIR=/usr - name: Build hemeXtract run: | From dce4c5c23bf428d18f24754b6de92b3a195c15b2 Mon Sep 17 00:00:00 2001 From: Ivan Pribec Date: Tue, 23 May 2023 20:17:10 +0200 Subject: [PATCH 13/16] Remove manual makefile --- Makefile | 16 ---------------- 1 file changed, 16 deletions(-) delete mode 100644 Makefile diff --git a/Makefile b/Makefile deleted file mode 100644 index 183c0bb..0000000 --- a/Makefile +++ /dev/null @@ -1,16 +0,0 @@ - -CXX=g++ -Wall -CXXFLAGS=-O3 -CXXFLAGS_DEBUG=-g - -SRC=hemeXtract.cc - -hemeXtract: $(SRC) - $(CXX) $(CXXFLAGS) -o $@ $< - -debug: $(SRC) - $(CXX) $(CXXFLAGS_DEBUG) -o $@ $< - -.phony: clean -clean: - rm -f hemeXtract \ No newline at end of file From 9b36028b185ccbd275a525742b657bc91d718b94 Mon Sep 17 00:00:00 2001 From: Ivan Pribec Date: Tue, 23 May 2023 20:26:19 +0200 Subject: [PATCH 14/16] Remove xdr header --- hemeXtract.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/hemeXtract.cc b/hemeXtract.cc index 0283c8f..cccbe21 100644 --- a/hemeXtract.cc +++ b/hemeXtract.cc @@ -8,7 +8,7 @@ #include #include -#include +//#include #include #include "Mapping.h" From 43b02bf5200a743af0fafc40e25a8ef66c7a397c Mon Sep 17 00:00:00 2001 From: Ivan Pribec Date: Tue, 23 May 2023 23:39:36 +0200 Subject: [PATCH 15/16] Decouple and reorganize headers --- CMakeLists.txt | 9 +- HemeLBExtractionFile.h | 606 ------------------ Snapshot.h | 238 ------- include/Compare.h | 28 + include/HemeLBExtractionFile.h | 150 +++++ .../HemeLBExtractionFileTypes.h | 16 +- include/Mapping.h | 12 + Site.h => include/Site.h | 33 +- include/Snapshot.h | 52 ++ Vector3.h => include/Vector3.h | 36 +- Compare.h => src/Compare.cc | 13 +- src/HemeLBExtractionFile.cc | 550 ++++++++++++++++ Mapping.h => src/Mapping.cc | 41 +- src/Snapshot.cc | 210 ++++++ hemeXtract.cc => src/hemeXtract.cc | 20 +- 15 files changed, 1064 insertions(+), 950 deletions(-) delete mode 100644 HemeLBExtractionFile.h delete mode 100644 Snapshot.h create mode 100644 include/Compare.h create mode 100644 include/HemeLBExtractionFile.h rename HemeLBExtractionFileTypes.h => include/HemeLBExtractionFileTypes.h (78%) create mode 100644 include/Mapping.h rename Site.h => include/Site.h (59%) create mode 100644 include/Snapshot.h rename Vector3.h => include/Vector3.h (69%) rename Compare.h => src/Compare.cc (97%) create mode 100644 src/HemeLBExtractionFile.cc rename Mapping.h => src/Mapping.cc (61%) create mode 100644 src/Snapshot.cc rename hemeXtract.cc => src/hemeXtract.cc (97%) diff --git a/CMakeLists.txt b/CMakeLists.txt index 0cdfec6..af85fdc 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -14,8 +14,13 @@ if(CMAKE_SYSTEM_NAME STREQUAL "Linux") find_package(TIRPC REQUIRED) endif() - -add_executable(hemeXtract hemeXtract.cc) +add_executable(hemeXtract + src/hemeXtract.cc + src/Mapping.cc + src/Snapshot.cc + src/Compare.cc + src/HemeLBExtractionFile.cc) +target_include_directories(hemeXtract PRIVATE include/) target_include_directories(hemeXtract PRIVATE ${argp_INCLUDE_DIRS}) target_link_libraries(hemeXtract PRIVATE ${argp_LIBRARIES}) diff --git a/HemeLBExtractionFile.h b/HemeLBExtractionFile.h deleted file mode 100644 index 5131eb0..0000000 --- a/HemeLBExtractionFile.h +++ /dev/null @@ -1,606 +0,0 @@ -#ifndef INCLUDED_HEMELBEXTRACTIONFILE_H -#define INCLUDED_HEMELBEXTRACTIONFILE_H - -#include -#include -#include -#include -#include -#include -#include -#include -//#include -#include -#include -#include - -#include "HemeLBExtractionFileTypes.h" -#include "Snapshot.h" -#include "Vector3.h" - -class HemeLBExtractionFile -{ - public: - HemeLBExtractionFile(char *fname, double step_length, double scaling, bool verbose, Vector3 *translate) - { - this->step_length = step_length; - this->scaling = scaling; - this->bool_verbose = verbose; - bool_correctly_initialised = false; - snapshot = NULL; - in = fopen(fname, "rb"); - if (in == NULL) { - fprintf(stderr, "Could not open file '%s'. Does it exist?\n", fname); - return; - } - xdrstdio_create(&xdrs, in, XDR_DECODE); - - has_velocity = false; - has_pressure = false; - has_shearstress = false; - - // Read the file header, populating the 'header' class variable. Determines whether file is "normal" or colloids - int r; - r = read_header(); - if(r != 0) { - fprintf(stderr, "Problem reading header for file: '%s'.\n", fname); - return; - } - - // If file is a colloids file, return now (no more info to be read in) - if(header->is_colloids_file == true) { - bool_correctly_initialised = true; - return; - } - - /// Otherwise, if file is "normal", read the field header - r = read_field_header(); - if(r != 0) { - fprintf(stderr, "Problem reading field header for file: '%s'.\n", fname); - return; - } - - // If specified, convert lattice translation vector (real units) into grid units (via rounding based on voxel size) - if(translate != NULL) { - this->translate_grid_x = (int32_t)(translate->get_x()/header->voxelsz); - this->translate_grid_y = (int32_t)(translate->get_y()/header->voxelsz); - this->translate_grid_z = (int32_t)(translate->get_z()/header->voxelsz); - if(bool_verbose == true) { - fprintf(stderr, "Lattice translation vector (%e,%e,%e) in real units, converted to (%d,%d,%d) in grid units (voxel size = %e).\n", translate->get_x(), translate->get_y(), translate->get_z(), translate_grid_x, translate_grid_y, translate_grid_z, header->voxelsz); - } - } else { - this->translate_grid_x = 0; - this->translate_grid_y = 0; - this->translate_grid_z = 0; - } - - // Allocate a snapshot of the right dimensions - snapshot = new Snapshot(header, field_header); - - // Get the time of the next snapshot - bool next = read_time_next(); - if(next == false) { - bool_no_more_snapshots = true; - } else { - bool_no_more_snapshots = false; - load_next_snapshot(); // Load in the first snapshot - bool_correctly_initialised = true; - } - } - - ~HemeLBExtractionFile() - { - xdr_destroy(&xdrs); - fclose(in); - } - - bool read_time_next() - { - // Check if we have reached the end of the file - if (feof(in)) { - if(bool_verbose == true) { - fprintf(stderr, "# Reached end of file.\n"); - } - return false; - } - - // Work around for problem with xdr_u_long() not working - uint32_t dt1, dt2; - uint64_t step; - xdr_u_int(&xdrs, &dt1); - xdr_u_int(&xdrs, &dt2); - step = ((uint64_t)dt1)<<32 | dt2; - - time_next = step * step_length; - return true; - } - - int load_next_snapshot() - { - if(bool_no_more_snapshots == true) { - if(bool_verbose == true) { - fprintf(stderr, "Warning: No more snapshots left in file.\n"); - } - return 1; - } - - uint32_t gridposx, gridposy, gridposz; - float value; - - snapshot->set_timestep(time_next); - for(unsigned int s = 0; s < header->num_sites; s++) { - // Read in grid position of site - xdr_u_int(&xdrs, &gridposx); - xdr_u_int(&xdrs, &gridposy); - xdr_u_int(&xdrs, &gridposz); - - // Apply any requested translation of this xtr file directly to the site positions - gridposx += translate_grid_x; - gridposy += translate_grid_y; - gridposz += translate_grid_z; - - snapshot->site_set(s, gridposx, gridposy, gridposz); - for(unsigned int i = 0; i < header->num_columns; i++) { - // Read value and add it to the appropriate data column (scaling it by the specified scaling) - xdr_float(&xdrs, &value); - snapshot->column_set_plus_offset(i, s, value * this->scaling); - } - } - bool next = read_time_next(); - if(next == false) { - bool_no_more_snapshots = true; - if(bool_verbose == true) {fprintf(stderr, "No more snapshots. (Note that final snapshot may have been incomplete)\n");} - } - return 0; - } - - // Workaround for XDR's broken long reading ability - void xdr_long(XDR *xdrs, uint64_t *ret) - { - uint32_t half1, half2; - xdr_u_int(xdrs, &half1); - xdr_u_int(xdrs, &half2); - *ret = ((uint64_t)half1)<<32 | half2; // Assume big-endianness - } - - - void read_and_print_colloids(FILE *outfile) - { - uint32_t headerLen, recordLen; - uint64_t dsetLen, timeStep; - uint64_t id, rank; - double X, Y, Z, VX, VY, VZ; - - while(!feof(in)) { - - // Reinitialise reading vars to catch partially written header - headerLen = 0; - recordLen = 0; - dsetLen = 0; - timeStep = 0; - - // Read timestep header - xdr_u_int(&xdrs, &headerLen); - xdr_u_int(&xdrs, &recordLen); - xdr_long(&xdrs, &dsetLen); - xdr_long(&xdrs, &timeStep); - - // Check for end of file here (to avoid spurious EOF output) - if(feof(in)) { - break; - } - - if(recordLen <= 0) { - break; - } - - if(bool_verbose == true) { - fprintf(outfile, "# headerLen: %u recordLen %u dsetLen %ld timeStep: %ld\n", headerLen, recordLen, dsetLen, timeStep); - } - - // Calc. num particles from the record length data - uint32_t num_particles = dsetLen/recordLen; - - for (uint32_t i = 0; i < num_particles; i++) { - // Read particle data - xdr_long(&xdrs, &id); - xdr_long(&xdrs, &rank); - xdr_double(&xdrs, &X); - xdr_double(&xdrs, &Y); - xdr_double(&xdrs, &Z); - xdr_double(&xdrs, &VX); - xdr_double(&xdrs, &VY); - xdr_double(&xdrs, &VZ); - - // Rescale the positions by the scaling factor provided to hemeXtract - X *= this->scaling; - Y *= this->scaling; - Z *= this->scaling; - -// fprintf(outfile, "TIME: %ld ID: %ld RANK: %ld A0: %e Ah: %e X: %e Y: %e Z: %e\n", timeStep, id, rank, A0, Ah, X, Y, Z); - fprintf(outfile, "%ld %ld %ld %.13e %.13e %.13e %.13e %.13e %.13e\n", timeStep, id, rank, X, Y, Z, VX, VY, VZ); - } - } - if(bool_verbose == true) { - fprintf(stderr, "# Reached end of file.\n"); - } - } - - bool correctly_initialised() - { - return bool_correctly_initialised; - } - - bool is_colloids_file() - { - return header->is_colloids_file; - } - - bool no_more_snapshots() { - return bool_no_more_snapshots; - } - - uint64_t get_num_sites() - { - return header->num_sites; - } - - double get_voxelsz() - { - return header->voxelsz; - } - - bool hasVelocity() - { - return has_velocity; - } - - bool hasShearStress() - { - return has_shearstress; - } - - bool hasPressure() - { - return has_pressure; - } - - double get_scaling() - { - return scaling; - } - - double get_scalar_quantity(uint32_t column_index, uint64_t site_index) - { - return snapshot->get(column_index, site_index); - } - - double get_interpolated_scalar_quantity(uint32_t column_index, lattice_map *map) - { - double average_existing = 0; - uint32_t num_existing = 0; - for(uint32_t i = 0; i < 8; i++) { - if(map->index[i].exists == true) { - average_existing += get_scalar_quantity(column_index, map->index[i].index); - num_existing++; - } - } - if(num_existing == 0) { - return 0; - } - average_existing /= num_existing; - - // Get the scalars at the 8 lattice sites forming the cube - double s[8]; - for(uint32_t i = 0; i < 8; i++) { - // If there was no corresponding site_index, the site was not measured, so set it to the average - // value of the existing sites in the map - if(map->index[i].exists == false) { - s[i] = average_existing; - } else { - s[i] = get_scalar_quantity(column_index, map->index[i].index); - } - } - - // Trilinear interpolation within the cube - return s[0] * (1 - map->a) * (1 - map->b) * (1 - map->c) - + s[1] * (map->a * (1 - map->b) * (1 - map->c)) - + s[2] * ((1 - map->a) * map->b * (1 - map->c)) - + s[3] * ((1 - map->a) * (1 - map->b) * map->c) - + s[4] * (map->a * (1 - map->b) * map->c) - + s[5] * ((1 - map->a) * map->b * map->c) - + s[6] * (map->a * map->b * (1 - map->c)) - + s[7] * (map->a * map->b * map->c); - } - void get_vector_quantity(uint32_t column_index, uint64_t site_index, Vector3 *returned_val) - { - double vx, vy, vz; - vx = get_scalar_quantity(column_index, site_index); - vy = get_scalar_quantity(column_index + 1, site_index); - vz = get_scalar_quantity(column_index + 2, site_index); - returned_val->set(vx, vy, vz); - } - - void get_interpolated_vector_quantity(uint32_t column_index, lattice_map *map, Vector3 *returned_val) - { - double vx, vy, vz; - vx = get_interpolated_scalar_quantity(column_index, map); - vy = get_interpolated_scalar_quantity(column_index + 1, map); - vz = get_interpolated_scalar_quantity(column_index + 2, map); - returned_val->set(vx, vy, vz); - } - - void get_pressure(uint64_t site_index, double *returned_val) - { - *returned_val = get_scalar_quantity(column_pressure, site_index); - } - - void get_velocity(uint64_t site_index, Vector3 *returned_val) - { - get_vector_quantity(column_velocity, site_index, returned_val); - } - - void get_shearstress(uint64_t site_index, double *returned_val) - { - *returned_val = get_scalar_quantity(column_shearstress, site_index); - } - - void get_interpolated_pressure(lattice_map *map, double *returned_val) - { - *returned_val = get_interpolated_scalar_quantity(column_pressure, map); - } - - void get_interpolated_velocity(lattice_map *map, Vector3 *returned_val) - { - get_interpolated_vector_quantity(column_velocity, map, returned_val); - } - - void get_interpolated_shearstress(lattice_map *map, double *returned_val) - { - *returned_val = get_interpolated_scalar_quantity(column_shearstress, map); - } - - Site * get_sites() - { - return snapshot->get_sites(); - } - - SiteIndex * get_site_indices(Site *list, uint64_t list_size) - { - return snapshot->get_site_indices(list, list_size); - } - - SiteIndex * get_site_indices_hashed_lookup(Site *list, uint64_t list_size) - { - return snapshot->get_site_indices_hashed_lookup(list, list_size, bool_verbose); - } - - double get_time() - { - return snapshot->get_timestep(); - } - - double get_time_next() - { - return time_next; - } - - void print_header(FILE *outfile) - { - fprintf(outfile, "\n# HEADER:\n# -------\n"); - fprintf(outfile, "# version=%u\n", header->version); - fprintf(outfile, "# voxelsz=%f\n", header->voxelsz); - fprintf(outfile, "# originx=%f\n", header->originx); - fprintf(outfile, "# originy=%f\n", header->originy); - fprintf(outfile, "# originz=%f\n", header->originz); - fprintf(outfile, "# num_sites=%lu\n", header->num_sites); - fprintf(outfile, "# field_count=%u\n", header->field_count); - fprintf(outfile, "# field_header_length=%u\n", header->field_header_length); - } - - void print_field_header(FILE *outfile) - { - fprintf(outfile, "\n# FIELD HEADERS:\n# -------\n"); - for(unsigned int i = 0; i < header->field_count; i++) { - fprintf(outfile, "# fieldname=%s\n# num_floats=%u\n# offset=%lf\n# ---\n", field_header[i].name, field_header[i].num_floats, field_header[i].offset); - } - fprintf(outfile, "# Number of 'columns' per snapshot = %d\n", header->num_columns); - } - - void print_column_headings(FILE *outfile) - { - fprintf(outfile, "# step | grid_x | grid_y | grid_z"); - for(unsigned int i = 0; i < header->field_count; i++) { - if(field_header[i].num_floats == 1) { - fprintf(outfile, " | %s", field_header[i].name); - } else { - for(unsigned int j = 0; j < field_header[i].num_floats; j++) { - fprintf(outfile, " | %s(%d)", field_header[i].name, j); - } - } - } - fprintf(outfile, "\n"); - } - void print_stats_column_headings(FILE *outfile) - { - fprintf(outfile, "# step"); - for(unsigned int i = 0; i < header->field_count; i++) { - if(field_header[i].num_floats == 1) { - fprintf(outfile, " | %s [average | stdev | max | min]", field_header[i].name); - } else { - for(unsigned int j = 0; j < field_header[i].num_floats; j++) { - fprintf(outfile, " | %s(%d) [average | stdev | max | min]", field_header[i].name, j); - } - } - } - fprintf(outfile, "\n"); - } - - void print_all(FILE *outfile) - { - snapshot->print(outfile); - } - - void print_stats(FILE *outfile) - { - snapshot->print_stats(outfile); - } - - private: - /* True if file has been opened and the headers read without problem */ - bool bool_correctly_initialised; - - /* True if user wants verbose output. If false, means quiet. */ - bool bool_verbose; - - /* All info contained in the file header */ - HEADER *header; - FIELD_HEADER *field_header; - - /** The step length for this file (has to be provided by the user since Heme extraction files don't store this value...) */ - double step_length; - - /** The velocity scaling factor */ - double scaling; - - /** The currently loaded snapshot */ - Snapshot *snapshot; - - /** The open extraction file's stream */ - FILE *in; - - /** The XDR reader for this file */ - XDR xdrs; - - /** The time of the next snapshot (next to be loaded) */ - double time_next; - - /** True when no more snapshots can be loaded */ - bool bool_no_more_snapshots; - - /** True if the file contains these field types */ - bool has_velocity, has_shearstress, has_pressure; - - /** The column indices in which to find each field type */ - uint32_t column_velocity; - uint32_t column_pressure; - uint32_t column_shearstress; - - /** The translation vector (in grid units) to be applied to this lattice (if requested by --translate option) */ - int32_t translate_grid_x; - int32_t translate_grid_y; - int32_t translate_grid_z; - - /** Magic numbers designating a heme file, a heme extraction file and a colloid file */ - static const uint32_t heme_magic = 0x686C6221; // hlb! - static const uint32_t extract_magic = 0x78747204; // xtr - static const uint32_t colloid_magic = 0x636F6C04; // colloid - - /** Reads a HemeLB extraction file header. Checks that the magic numbers are correct. */ - int read_header() - { - uint32_t magic1 = 0, magic2 = 0; - uint32_t num_sites1, num_sites2; - - // Read heme magic number - xdr_u_int(&xdrs, &magic1); - if(magic1 != heme_magic) { - fprintf(stderr, "First uint32_t does not match heme magic number.\n"); - return 1; - } - - header = new HEADER(); - - // Read extraction magic number - xdr_u_int(&xdrs, &magic2); - if(magic2 == extract_magic) { - if(bool_verbose == true) { - fprintf(stderr, "Reading a normal extraction file.\n"); - } - - header->is_colloids_file = false; - - xdr_u_int(&xdrs, &(header->version)); - xdr_double(&xdrs, &(header->voxelsz)); - xdr_double(&xdrs, &(header->originx)); - xdr_double(&xdrs, &(header->originy)); - xdr_double(&xdrs, &(header->originz)); - - // There is something wrong with either Heme or XDR which means that xdr_u_long is not working (possibly even reading more/less than 8 bytes...) - xdr_u_int(&xdrs, &num_sites1); - xdr_u_int(&xdrs, &num_sites2); - - // Combine the two ints into a long (assuming bigendian-ness) - header->num_sites = ((uint64_t)num_sites1)<<32 | num_sites2; - - xdr_u_int(&xdrs, &(header->field_count)); - xdr_u_int(&xdrs, &(header->field_header_length)); - - } else if (magic2 == colloid_magic) { - if(bool_verbose == true) { - fprintf(stderr, "Reading a colloids extraction file.\n"); - } - header->is_colloids_file = true; - - // Check version - uint32_t version=0; - xdr_u_int(&xdrs, &version); - if(bool_verbose == true) { - printf("Colloids version %u\n", version); - } - } else { - fprintf(stderr, "Unknown format: Second uint32_t does not match extraction file or colloids file magic number.\n"); - return 1; - } - - // Everything went well - return 0; - } - - /** Reads the field headers and remembers which columns correspond to velocity, pressure and shearstress. */ - int read_field_header() - { - if(bool_verbose == true) { - printf("Reading field header...\n"); - } - int const maxlength = 100; // surely it can't be bigger than this...? - field_header = new FIELD_HEADER[header->field_count]; - for(unsigned int i = 0; i < header->field_count; i++) { - field_header[i].name = new char[maxlength]; - xdr_string(&xdrs, &(field_header[i].name), maxlength); - xdr_u_int(&xdrs, &(field_header[i].num_floats)); - xdr_double(&xdrs, &(field_header[i].offset)); - } - if(bool_verbose == true) { - printf("Finished reading field header...\n"); - } - - // Calc. number of 'columns' needed to represent a snapshot - uint32_t num_columns = 0; - for(unsigned int i = 0; i < header->field_count; i++) { - // If file contains velocity, pressure or shearstress, remember the corresponding column index for that field type. - if(strcmp(field_header[i].name, "velocity") == 0) { - has_velocity = true; - column_velocity = num_columns; - } else if(strcmp(field_header[i].name, "pressure") == 0) { - has_pressure = true; - column_pressure = num_columns; - } else if(strcmp(field_header[i].name, "shearstress") == 0) { - has_shearstress = true; - column_shearstress = num_columns; - } else if(strcmp(field_header[i].name, "d_shearstress") == 0) { // For "legacy" reasons only. Some old MCA runs have this field. - has_shearstress = true; - column_shearstress = num_columns; - } - - // This field has 'num_floats' columns to it - num_columns += field_header[i].num_floats; - } - header->num_columns = num_columns; - - // Everything went well - return 0; - } -}; - -#endif diff --git a/Snapshot.h b/Snapshot.h deleted file mode 100644 index 13e14c5..0000000 --- a/Snapshot.h +++ /dev/null @@ -1,238 +0,0 @@ -#ifndef INCLUDED_SNAPSHOT_H -#define INCLUDED_SNAPSHOT_H - -#include -#include -#include -#include -#include -#include -#include -#include -//#include -#include -#include -#include - -#include "HemeLBExtractionFileTypes.h" -#include "Site.h" -#include "Mapping.h" - -class Column { - public: - Column(unsigned int num_records, double offset) - { - this->num_records = num_records; - this->offset = offset; - stats_calcd = false; - records = new double[num_records]; - } - ~Column() - { - delete[] records; - } - - double get(uint64_t index) - { - if(index >= num_records) { - fprintf(stderr, "Error: index (%lu) exceeds number of records in column (%lu)\n", index, num_records); - exit(1); - } - return records[index]; - } - void set_plus_offset(uint64_t index, float value) - { - if(index >= num_records) { - fprintf(stderr, "Error: index (%lu) exceeds number of records in column (%lu)\n", index, num_records); - exit(1); - } - records[index] = value + offset; - stats_calcd = false; - } - void print(FILE *outfile, uint64_t index) { - if(index >= num_records) { - fprintf(stderr, "Error: index (%lu) exceeds number of records in column (%lu)\n", index, num_records); - exit(1); - } - fprintf(outfile, "%e", records[index]); - } - void print_stats(FILE *outfile) { - calc_stats(); - fprintf(outfile, "%e %e %e %e", average, stdev, max, min); - } - - private: - uint64_t num_records; - double *records; - double offset; - bool stats_calcd; - double average; - double stdev; - double max, min; - void calc_stats() - { - if(!stats_calcd) { - double sum = 0, sum2 = 0; - max = -INFINITY; - min = INFINITY; - for(uint64_t i = 0; i < num_records; i++) { - sum += records[i]; - sum2 += records[i] * records[i]; - if(max < records[i]) { - max = records[i]; - } - if(min > records[i]) { - min = records[i]; - } - } - average = sum / num_records; - stdev = sqrt(sum2/num_records - average*average); - stats_calcd = true; - } - } -}; - -class Snapshot { - public: - Snapshot(HEADER *header, FIELD_HEADER *field_header) - { - this->num_sites = header->num_sites; - this->num_columns = header->num_columns; - this->voxelsz = header->voxelsz; - sites = new Site[num_sites]; - columns = new Column*[num_columns]; - unsigned int c = 0; - for(unsigned int i = 0; i < header->field_count; i++) { - for(unsigned int bro=0; bro < field_header[i].num_floats; bro++) { - columns[c] = new Column(num_sites, field_header[i].offset); - c++; - } - } - } - - void set_timestep(double timestep) - { - this->timestep = timestep; - } - double get_timestep() - { - return timestep; - } - - void site_set(unsigned int site_index, uint32_t x, uint32_t y, uint32_t z) - { - sites[site_index].set(x, y, z); - } - - void column_set_plus_offset(uint32_t column_index, uint64_t site_index, double value) - { - if(column_index >= num_columns) { - fprintf(stderr, "Error: column index (%u) exceeds number of columns (%u)\n", column_index, num_columns); - exit(1); - } - columns[column_index]->set_plus_offset(site_index, value); - } - - double get(uint32_t column_index, uint64_t site_index) - { - return columns[column_index]->get(site_index); - } - - bool get_site_index(uint32_t a, uint32_t b, uint32_t c, uint64_t *site_index) - { - for(unsigned int i = 0; i < num_sites; i++) { - if(sites[i].equals(a, b, c)) { - *site_index = i; - return true; - } - } - - // If no such site can be found, return false - return false; - } - - int get_site_index(Site *s, uint64_t *site_index) - { - return get_site_index(s->x, s->y, s->z, site_index); - } - - Site * get_sites() - { - return sites; - } - - /* Very inefficient calculation of site_index corresponding to the given site. Only suitable for small data sets. */ - SiteIndex * get_site_indices(Site *list, uint64_t list_size) - { - SiteIndex *indices = new SiteIndex[list_size]; - for(uint64_t i = 0; i < list_size; i++) { - indices[i].exists = get_site_index(&list[i], &(indices[i].index)); - } - return indices; - } - - /** Builds hashtable to get indices corresponding to the given site. Necessary for very large data sets. */ - SiteIndex * get_site_indices_hashed_lookup(Site *list, uint64_t list_size, bool bool_verbose) - { - // Build the hash table - if(bool_verbose == true) { fprintf(stderr, "# Building hashtable...\n"); } - //std::tr1::unordered_map hashtable; //JM for C++11 capability - std::unordered_map hashtable; - for(uint64_t i = 0; i < num_sites; i++) { - std::ostringstream oss; - oss << sites[i].x << "," << sites[i].y << "," << sites[i].z; - std::string hashstr = oss.str(); - hashtable[hashstr] = i; - } - if(bool_verbose == true) { fprintf(stderr, "# ...done building hashtable\n"); } - - if(bool_verbose == true) { fprintf(stderr, "# Calculating mapping...\n"); } - SiteIndex *indices = new SiteIndex[list_size]; - for(uint64_t i = 0; i < list_size; i++) { - std::ostringstream oss; - oss << list[i].x << "," << list[i].y << "," << list[i].z; - std::string hashstr = oss.str(); - if(hashtable.find(hashstr) != hashtable.end()) { - indices[i].exists = true; - indices[i].index = hashtable[hashstr]; - } else { - indices[i].exists = false; - } - } - if(bool_verbose == true) { fprintf(stderr, "# ... done calculating mapping\n"); } - - return indices; - } - - void print(FILE *outfile) - { - for(uint64_t s = 0; s < num_sites; s++) { - fprintf(outfile, "%f ", timestep); - sites[s].print(outfile, this->voxelsz); - for(uint32_t c = 0; c < num_columns; c++) { - fprintf(outfile, " "); - columns[c]->print(outfile, s); - } - fprintf(outfile, "\n"); - } - fprintf(outfile, "\n"); - } - void print_stats(FILE *outfile) - { - fprintf(outfile, "%f", timestep); - for(uint32_t c = 0; c < num_columns; c++) { - fprintf(outfile, " "); - columns[c]->print_stats(outfile); - } - fprintf(outfile, "\n"); - } - private: - double timestep; - uint64_t num_sites; - uint32_t num_columns; - double voxelsz; - Column **columns; - Site *sites; -}; - -#endif diff --git a/include/Compare.h b/include/Compare.h new file mode 100644 index 0000000..e8c0458 --- /dev/null +++ b/include/Compare.h @@ -0,0 +1,28 @@ +#ifndef INCLUDED_COMPARE_H +#define INCLUDED_COMPARE_H + +#include + +#include "Vector3.h" +#include "HemeLBExtractionFileTypes.h" // lattice_map +#include "HemeLBExtractionFile.h" // HemeLBExtractionFile + +/** Compares two lattices, A and B, with mapping mapA_to_B, by calculating their crosscorrelation, and their L2 norm difference. */ +void compare(FILE *outfile, + lattice_map *mapA_to_B, + HemeLBExtractionFile *A, + HemeLBExtractionFile *B, + int minexistent, + bool normalize_correl); + +/** Calculates the difference between each site in A, and the corresponding (trilinearly interpolated) point in B */ +void diff(FILE *outfile, + lattice_map *mapA_to_B, + HemeLBExtractionFile *A, + HemeLBExtractionFile *B, + int minexistent, + Vector3 *project, + bool relativeErr, + bool verbose); + +#endif diff --git a/include/HemeLBExtractionFile.h b/include/HemeLBExtractionFile.h new file mode 100644 index 0000000..4f38ecd --- /dev/null +++ b/include/HemeLBExtractionFile.h @@ -0,0 +1,150 @@ +#ifndef INCLUDED_HEMELBEXTRACTIONFILE_H +#define INCLUDED_HEMELBEXTRACTIONFILE_H + +#include +#include + +#include // must precede xdr.h +#include // XDR + +#include "Vector3.h" // Vector3 +#include "Site.h" // Site +#include "HemeLBExtractionFileTypes.h" // lattice_map, SiteIndex +#include "Snapshot.h" // Snapshot + +class HemeLBExtractionFile +{ + public: + // Constructor + HemeLBExtractionFile(char *fname, + double step_length, + double scaling, + bool verbose, + Vector3 *translate); + // Destructor + ~HemeLBExtractionFile(); + + bool read_time_next(); + + int load_next_snapshot(); + + // Workaround for XDR's broken long reading ability + void xdr_long(XDR *xdrs, uint64_t *ret); + + void read_and_print_colloids(FILE *outfile); + + bool correctly_initialised(); + + bool is_colloids_file(); + + bool no_more_snapshots(); + + uint64_t get_num_sites(); + + double get_voxelsz(); + + bool hasVelocity(); + + bool hasShearStress(); + + bool hasPressure(); + + double get_scaling(); + + double get_scalar_quantity(uint32_t column_index, uint64_t site_index); + + double get_interpolated_scalar_quantity( + uint32_t column_index, + lattice_map *map); + + void get_vector_quantity( + uint32_t column_index, + uint64_t site_index, + Vector3 *returned_val); + + void get_interpolated_vector_quantity( + uint32_t column_index, + lattice_map *map, + Vector3 *returned_val); + + void get_pressure(uint64_t site_index, double *returned_val); + + void get_velocity(uint64_t site_index, Vector3 *returned_val); + + void get_shearstress(uint64_t site_index, double *returned_val); + + void get_interpolated_pressure(lattice_map *map, double *returned_val); + + void get_interpolated_velocity(lattice_map *map, Vector3 *returned_val); + + void get_interpolated_shearstress(lattice_map *map, double *returned_val); + + Site * get_sites(); + + SiteIndex * get_site_indices(Site *list, uint64_t list_size); + + SiteIndex * get_site_indices_hashed_lookup(Site *list, uint64_t list_size); + + double get_time(); + double get_time_next(); + + void print_header(FILE *outfile); + void print_field_header(FILE *outfile); + void print_column_headings(FILE *outfile); + void print_stats_column_headings(FILE *outfile); + void print_all(FILE *outfile); + void print_stats(FILE *outfile); + + private: + /* True if file has been opened and the headers read without problem */ + bool bool_correctly_initialised; + + /* True if user wants verbose output. If false, means quiet. */ + bool bool_verbose; + + /* All info contained in the file header */ + HEADER *header; + FIELD_HEADER *field_header; + + /** The step length for this file (has to be provided by the user since Heme extraction files don't store this value...) */ + double step_length; + + /** The velocity scaling factor */ + double scaling; + + /** The currently loaded snapshot */ + Snapshot *snapshot; + + /** The open extraction file's stream */ + FILE *in; + + /** The XDR reader for this file */ + XDR xdrs; + + /** The time of the next snapshot (next to be loaded) */ + double time_next; + + /** True when no more snapshots can be loaded */ + bool bool_no_more_snapshots; + + /** True if the file contains these field types */ + bool has_velocity, has_shearstress, has_pressure; + + /** The column indices in which to find each field type */ + uint32_t column_velocity; + uint32_t column_pressure; + uint32_t column_shearstress; + + /** The translation vector (in grid units) to be applied to this lattice (if requested by --translate option) */ + int32_t translate_grid_x; + int32_t translate_grid_y; + int32_t translate_grid_z; + + /** Reads a HemeLB extraction file header. Checks that the magic numbers are correct. */ + int read_header(); + + /** Reads the field headers and remembers which columns correspond to velocity, pressure and shearstress. */ + int read_field_header(); +}; + +#endif diff --git a/HemeLBExtractionFileTypes.h b/include/HemeLBExtractionFileTypes.h similarity index 78% rename from HemeLBExtractionFileTypes.h rename to include/HemeLBExtractionFileTypes.h index d9ab48c..11a980c 100644 --- a/HemeLBExtractionFileTypes.h +++ b/include/HemeLBExtractionFileTypes.h @@ -1,6 +1,8 @@ #ifndef INCLUDED_HEMELBEXTRACTIONFILETYPES_H #define INCLUDED_HEMELBEXTRACTIONFILETYPES_H +#include + typedef struct { uint32_t version; double voxelsz; @@ -25,19 +27,15 @@ typedef struct { class SiteIndex { public: - SiteIndex() - { - index = 0; - exists = false; - } - uint64_t index; - bool exists; + SiteIndex() {} + uint64_t index{0}; + bool exists{false}; }; typedef struct { - SiteIndex index[8]; - double a, b, c; + SiteIndex index[8]; + double a, b, c; } lattice_map; #endif diff --git a/include/Mapping.h b/include/Mapping.h new file mode 100644 index 0000000..4f7c41b --- /dev/null +++ b/include/Mapping.h @@ -0,0 +1,12 @@ +#ifndef INCLUDED_MAPPING_H +#define INCLUDED_MAPPING_H + +#include "HemeLBExtractionFileTypes.h" // lattice_map +#include "HemeLBExtractionFile.h" // HemeLBExtractionFile + +/* Returns a mapping for each site in A to a corresponding + * trilinearly interpolated position in B. + */ +lattice_map *get_mapping(HemeLBExtractionFile *A, HemeLBExtractionFile *B); + +#endif diff --git a/Site.h b/include/Site.h similarity index 59% rename from Site.h rename to include/Site.h index 9afd5c9..49deda5 100644 --- a/Site.h +++ b/include/Site.h @@ -1,27 +1,14 @@ #ifndef INCLUDED_SITE_H #define INCLUDED_SITE_H -#include -#include -#include -#include -#include -#include -#include -#include -//#include -#include -#include -#include - -class Site { +#include +#include + +class Site +{ public: - Site() - { - this->x = 0; - this->y = 0; - this->z = 0; - } + + Site() {} void set(uint32_t x, uint32_t y, uint32_t z) { this->x = x; @@ -29,8 +16,6 @@ class Site { this->z = z; } - uint32_t x, y, z; - void print(FILE *outfile) { fprintf(outfile, "%u %u %u", x, y, z); } @@ -39,7 +24,7 @@ class Site { fprintf(outfile, "%f %f %f", x * voxelsz, y * voxelsz, z * voxelsz); } - bool equals(uint32_t a, uint32_t b, uint32_t c) + inline bool equals(uint32_t a, uint32_t b, uint32_t c) { if(a != x) { return false; @@ -57,6 +42,8 @@ class Site { { return equals(s->x, s->y, s->z); } + + uint32_t x{0}, y{0}, z{0}; }; #endif diff --git a/include/Snapshot.h b/include/Snapshot.h new file mode 100644 index 0000000..fab81d8 --- /dev/null +++ b/include/Snapshot.h @@ -0,0 +1,52 @@ +#ifndef INCLUDED_SNAPSHOT_H +#define INCLUDED_SNAPSHOT_H + +#include +#include + +#include "Site.h" // Site +#include "HemeLBExtractionFileTypes.h" // HEADER, FIELD_HEADER, SiteIndex + +class Column; // Forward declaration (could be private potentially) + +class Snapshot { + public: + Snapshot(HEADER *header, FIELD_HEADER *field_header); + + void set_timestep(double timestep); + double get_timestep(); + + void site_set(unsigned int site_index, uint32_t x, uint32_t y, uint32_t z); + + void column_set_plus_offset(uint32_t column_index, uint64_t site_index, double value); + + double get(uint32_t column_index, uint64_t site_index); + + bool get_site_index(uint32_t a, uint32_t b, uint32_t c, uint64_t *site_index); + + int get_site_index(Site *s, uint64_t *site_index); + + Site *get_sites(); + + /* Very inefficient calculation of site_index corresponding to the + * given site. Only suitable for small data sets. */ + SiteIndex *get_site_indices(Site *list, uint64_t list_size); + + /* Builds hashtable to get indices corresponding to the given + * site. Necessary for very large data sets. */ + SiteIndex *get_site_indices_hashed_lookup(Site *list, + uint64_t list_size, bool bool_verbose); + + void print(FILE *outfile); + void print_stats(FILE *outfile); + + private: + double timestep; + uint64_t num_sites; + uint32_t num_columns; + double voxelsz; + Column **columns; + Site *sites; +}; + +#endif diff --git a/Vector3.h b/include/Vector3.h similarity index 69% rename from Vector3.h rename to include/Vector3.h index 5ad8cb3..0455a72 100644 --- a/Vector3.h +++ b/include/Vector3.h @@ -1,35 +1,18 @@ #ifndef INCLUDED_VECTOR3_H #define INCLUDED_VECTOR3_H -#include -#include -#include -#include -#include -#include -#include -#include -//#include -#include -#include -#include +#include +#include /** Reinventing the wheel... */ class Vector3 { public: - Vector3() - { - x = 0; - y = 0; - z = 0; - } - Vector3(double x, double y, double z) - { - this->x = x; - this->y = y; - this->z = z; - } + + Vector3() {} + + Vector3(double x_, double y_, double z_) + : x(x_), y(y_), z(z_) {} void set(double x, double y, double z) { @@ -107,9 +90,8 @@ class Vector3 } private: - //static const double tol = 1e-10; - static constexpr double tol = 1e-10; //JM for C++ 11 - double x, y, z; + static constexpr double tol{1e-10}; //JM for C++ 11 + double x{0}, y{0}, z{0}; }; #endif diff --git a/Compare.h b/src/Compare.cc similarity index 97% rename from Compare.h rename to src/Compare.cc index 05896d6..7cfec33 100644 --- a/Compare.h +++ b/src/Compare.cc @@ -1,13 +1,6 @@ -#ifndef INCLUDED_COMPARE_H -#define INCLUDED_COMPARE_H +#include "Compare.h" -#include -#include -#include -#include -#include -#include -#include +#include /** Compares two lattices, A and B, with mapping mapA_to_B, by calculating their crosscorrelation, and their L2 norm difference. */ void compare(FILE *outfile, lattice_map *mapA_to_B, HemeLBExtractionFile *A, HemeLBExtractionFile *B, int minexistent, bool normalize_correl) @@ -225,5 +218,3 @@ void diff(FILE *outfile, lattice_map *mapA_to_B, HemeLBExtractionFile *A, HemeLB } } } - -#endif diff --git a/src/HemeLBExtractionFile.cc b/src/HemeLBExtractionFile.cc new file mode 100644 index 0000000..b3e5efd --- /dev/null +++ b/src/HemeLBExtractionFile.cc @@ -0,0 +1,550 @@ +#include "HemeLBExtractionFile.h" + +#include + + +/* Magic numbers designating: + * - a heme file, + * - a heme extraction file, + * - and a colloid file + */ +static const uint32_t heme_magic = 0x686C6221; // hlb! +static const uint32_t extract_magic = 0x78747204; // xtr +static const uint32_t colloid_magic = 0x636F6C04; // colloid + +HemeLBExtractionFile::HemeLBExtractionFile( + char *fname, + double step_length, + double scaling, + bool verbose, + Vector3 *translate) +{ + this->step_length = step_length; + this->scaling = scaling; + this->bool_verbose = verbose; + bool_correctly_initialised = false; + snapshot = NULL; + in = fopen(fname, "rb"); + if (in == NULL) { + fprintf(stderr, "Could not open file '%s'. Does it exist?\n", fname); + return; + } + xdrstdio_create(&xdrs, in, XDR_DECODE); + + has_velocity = false; + has_pressure = false; + has_shearstress = false; + + // Read the file header, populating the 'header' class variable. Determines whether file is "normal" or colloids + int r; + r = read_header(); + if(r != 0) { + fprintf(stderr, "Problem reading header for file: '%s'.\n", fname); + return; + } + + // If file is a colloids file, return now (no more info to be read in) + if(header->is_colloids_file == true) { + bool_correctly_initialised = true; + return; + } + + /// Otherwise, if file is "normal", read the field header + r = read_field_header(); + if(r != 0) { + fprintf(stderr, "Problem reading field header for file: '%s'.\n", fname); + return; + } + + // If specified, convert lattice translation vector (real units) into grid units (via rounding based on voxel size) + if(translate != NULL) { + this->translate_grid_x = (int32_t)(translate->get_x()/header->voxelsz); + this->translate_grid_y = (int32_t)(translate->get_y()/header->voxelsz); + this->translate_grid_z = (int32_t)(translate->get_z()/header->voxelsz); + if(bool_verbose == true) { + fprintf(stderr, "Lattice translation vector (%e,%e,%e) in real units, converted to (%d,%d,%d) in grid units (voxel size = %e).\n", translate->get_x(), translate->get_y(), translate->get_z(), translate_grid_x, translate_grid_y, translate_grid_z, header->voxelsz); + } + } else { + this->translate_grid_x = 0; + this->translate_grid_y = 0; + this->translate_grid_z = 0; + } + + // Allocate a snapshot of the right dimensions + snapshot = new Snapshot(header, field_header); + + // Get the time of the next snapshot + bool next = read_time_next(); + if(next == false) { + bool_no_more_snapshots = true; + } else { + bool_no_more_snapshots = false; + load_next_snapshot(); // Load in the first snapshot + bool_correctly_initialised = true; + } +} + +HemeLBExtractionFile::~HemeLBExtractionFile() +{ + xdr_destroy(&xdrs); + fclose(in); +} + +bool HemeLBExtractionFile::read_time_next() +{ + // Check if we have reached the end of the file + if (feof(in)) { + if(bool_verbose == true) { + fprintf(stderr, "# Reached end of file.\n"); + } + return false; + } + + // Work around for problem with xdr_u_long() not working + uint32_t dt1, dt2; + uint64_t step; + xdr_u_int(&xdrs, &dt1); + xdr_u_int(&xdrs, &dt2); + step = ((uint64_t)dt1)<<32 | dt2; + + time_next = step * step_length; + return true; +} + +int HemeLBExtractionFile::load_next_snapshot() +{ + if(bool_no_more_snapshots == true) { + if(bool_verbose == true) { + fprintf(stderr, "Warning: No more snapshots left in file.\n"); + } + return 1; + } + + uint32_t gridposx, gridposy, gridposz; + float value; + + snapshot->set_timestep(time_next); + for(unsigned int s = 0; s < header->num_sites; s++) { + // Read in grid position of site + xdr_u_int(&xdrs, &gridposx); + xdr_u_int(&xdrs, &gridposy); + xdr_u_int(&xdrs, &gridposz); + + // Apply any requested translation of this xtr file directly to the site positions + gridposx += translate_grid_x; + gridposy += translate_grid_y; + gridposz += translate_grid_z; + + snapshot->site_set(s, gridposx, gridposy, gridposz); + for(unsigned int i = 0; i < header->num_columns; i++) { + // Read value and add it to the appropriate data column (scaling it by the specified scaling) + xdr_float(&xdrs, &value); + snapshot->column_set_plus_offset(i, s, value * this->scaling); + } + } + bool next = read_time_next(); + if(next == false) { + bool_no_more_snapshots = true; + if(bool_verbose == true) {fprintf(stderr, "No more snapshots. (Note that final snapshot may have been incomplete)\n");} + } + return 0; +} + +// Workaround for XDR's broken long reading ability +void HemeLBExtractionFile::xdr_long(XDR *xdrs, uint64_t *ret) +{ + uint32_t half1, half2; + xdr_u_int(xdrs, &half1); + xdr_u_int(xdrs, &half2); + *ret = ((uint64_t)half1)<<32 | half2; // Assume big-endianness +} + + +void HemeLBExtractionFile::read_and_print_colloids(FILE *outfile) +{ + uint32_t headerLen, recordLen; + uint64_t dsetLen, timeStep; + uint64_t id, rank; + double X, Y, Z, VX, VY, VZ; + + while(!feof(in)) { + + // Reinitialise reading vars to catch partially written header + headerLen = 0; + recordLen = 0; + dsetLen = 0; + timeStep = 0; + + // Read timestep header + xdr_u_int(&xdrs, &headerLen); + xdr_u_int(&xdrs, &recordLen); + xdr_long(&xdrs, &dsetLen); + xdr_long(&xdrs, &timeStep); + + // Check for end of file here (to avoid spurious EOF output) + if(feof(in)) { + break; + } + + if(recordLen <= 0) { + break; + } + + if(bool_verbose == true) { + fprintf(outfile, "# headerLen: %u recordLen %u dsetLen %ld timeStep: %ld\n", headerLen, recordLen, dsetLen, timeStep); + } + + // Calc. num particles from the record length data + uint32_t num_particles = dsetLen/recordLen; + + for (uint32_t i = 0; i < num_particles; i++) { + // Read particle data + xdr_long(&xdrs, &id); + xdr_long(&xdrs, &rank); + xdr_double(&xdrs, &X); + xdr_double(&xdrs, &Y); + xdr_double(&xdrs, &Z); + xdr_double(&xdrs, &VX); + xdr_double(&xdrs, &VY); + xdr_double(&xdrs, &VZ); + + // Rescale the positions by the scaling factor provided to hemeXtract + X *= this->scaling; + Y *= this->scaling; + Z *= this->scaling; + +// fprintf(outfile, "TIME: %ld ID: %ld RANK: %ld A0: %e Ah: %e X: %e Y: %e Z: %e\n", timeStep, id, rank, A0, Ah, X, Y, Z); + fprintf(outfile, "%ld %ld %ld %.13e %.13e %.13e %.13e %.13e %.13e\n", timeStep, id, rank, X, Y, Z, VX, VY, VZ); + } + } + if(bool_verbose == true) { + fprintf(stderr, "# Reached end of file.\n"); + } +} + +bool HemeLBExtractionFile::correctly_initialised() +{ + return bool_correctly_initialised; +} + +bool HemeLBExtractionFile::is_colloids_file() +{ + return header->is_colloids_file; +} + +bool HemeLBExtractionFile::no_more_snapshots() { + return bool_no_more_snapshots; +} + +uint64_t HemeLBExtractionFile::get_num_sites() +{ + return header->num_sites; +} + +double HemeLBExtractionFile::get_voxelsz() +{ + return header->voxelsz; +} + +bool HemeLBExtractionFile::hasVelocity() +{ + return has_velocity; +} + +bool HemeLBExtractionFile::hasShearStress() +{ + return has_shearstress; +} + +bool HemeLBExtractionFile::hasPressure() +{ + return has_pressure; +} + +double HemeLBExtractionFile::get_scaling() +{ + return scaling; +} + +double HemeLBExtractionFile::get_scalar_quantity(uint32_t column_index, uint64_t site_index) +{ + return snapshot->get(column_index, site_index); +} + +double HemeLBExtractionFile::get_interpolated_scalar_quantity(uint32_t column_index, lattice_map *map) +{ + double average_existing = 0; + uint32_t num_existing = 0; + for(uint32_t i = 0; i < 8; i++) { + if(map->index[i].exists == true) { + average_existing += get_scalar_quantity(column_index, map->index[i].index); + num_existing++; + } + } + if(num_existing == 0) { + return 0; + } + average_existing /= num_existing; + + // Get the scalars at the 8 lattice sites forming the cube + double s[8]; + for(uint32_t i = 0; i < 8; i++) { + // If there was no corresponding site_index, the site was not measured, so set it to the average + // value of the existing sites in the map + if(map->index[i].exists == false) { + s[i] = average_existing; + } else { + s[i] = get_scalar_quantity(column_index, map->index[i].index); + } + } + + // Trilinear interpolation within the cube + return s[0] * (1 - map->a) * (1 - map->b) * (1 - map->c) + + s[1] * (map->a * (1 - map->b) * (1 - map->c)) + + s[2] * ((1 - map->a) * map->b * (1 - map->c)) + + s[3] * ((1 - map->a) * (1 - map->b) * map->c) + + s[4] * (map->a * (1 - map->b) * map->c) + + s[5] * ((1 - map->a) * map->b * map->c) + + s[6] * (map->a * map->b * (1 - map->c)) + + s[7] * (map->a * map->b * map->c); +} +void HemeLBExtractionFile::get_vector_quantity(uint32_t column_index, uint64_t site_index, Vector3 *returned_val) +{ + double vx, vy, vz; + vx = get_scalar_quantity(column_index, site_index); + vy = get_scalar_quantity(column_index + 1, site_index); + vz = get_scalar_quantity(column_index + 2, site_index); + returned_val->set(vx, vy, vz); +} + +void HemeLBExtractionFile::get_interpolated_vector_quantity(uint32_t column_index, lattice_map *map, Vector3 *returned_val) +{ + double vx, vy, vz; + vx = get_interpolated_scalar_quantity(column_index, map); + vy = get_interpolated_scalar_quantity(column_index + 1, map); + vz = get_interpolated_scalar_quantity(column_index + 2, map); + returned_val->set(vx, vy, vz); +} + +void HemeLBExtractionFile::get_pressure(uint64_t site_index, double *returned_val) +{ + *returned_val = get_scalar_quantity(column_pressure, site_index); +} + +void HemeLBExtractionFile::get_velocity(uint64_t site_index, Vector3 *returned_val) +{ + get_vector_quantity(column_velocity, site_index, returned_val); +} + +void HemeLBExtractionFile::get_shearstress(uint64_t site_index, double *returned_val) +{ + *returned_val = get_scalar_quantity(column_shearstress, site_index); +} + +void HemeLBExtractionFile::get_interpolated_pressure(lattice_map *map, double *returned_val) +{ + *returned_val = get_interpolated_scalar_quantity(column_pressure, map); +} + +void HemeLBExtractionFile::get_interpolated_velocity(lattice_map *map, Vector3 *returned_val) +{ + get_interpolated_vector_quantity(column_velocity, map, returned_val); +} + +void HemeLBExtractionFile::get_interpolated_shearstress(lattice_map *map, double *returned_val) +{ + *returned_val = get_interpolated_scalar_quantity(column_shearstress, map); +} + +Site * HemeLBExtractionFile::get_sites() +{ + return snapshot->get_sites(); +} + +SiteIndex * HemeLBExtractionFile::get_site_indices(Site *list, uint64_t list_size) +{ + return snapshot->get_site_indices(list, list_size); +} + +SiteIndex * HemeLBExtractionFile::get_site_indices_hashed_lookup(Site *list, uint64_t list_size) +{ + return snapshot->get_site_indices_hashed_lookup(list, list_size, bool_verbose); +} + +double HemeLBExtractionFile::get_time() +{ + return snapshot->get_timestep(); +} + +double HemeLBExtractionFile::get_time_next() +{ + return time_next; +} + +void HemeLBExtractionFile::print_header(FILE *outfile) +{ + fprintf(outfile, "\n# HEADER:\n# -------\n"); + fprintf(outfile, "# version=%u\n", header->version); + fprintf(outfile, "# voxelsz=%f\n", header->voxelsz); + fprintf(outfile, "# originx=%f\n", header->originx); + fprintf(outfile, "# originy=%f\n", header->originy); + fprintf(outfile, "# originz=%f\n", header->originz); + fprintf(outfile, "# num_sites=%lu\n", header->num_sites); + fprintf(outfile, "# field_count=%u\n", header->field_count); + fprintf(outfile, "# field_header_length=%u\n", header->field_header_length); +} + +void HemeLBExtractionFile::print_field_header(FILE *outfile) +{ + fprintf(outfile, "\n# FIELD HEADERS:\n# -------\n"); + for(unsigned int i = 0; i < header->field_count; i++) { + fprintf(outfile, "# fieldname=%s\n# num_floats=%u\n# offset=%lf\n# ---\n", field_header[i].name, field_header[i].num_floats, field_header[i].offset); + } + fprintf(outfile, "# Number of 'columns' per snapshot = %d\n", header->num_columns); +} + +void HemeLBExtractionFile::print_column_headings(FILE *outfile) +{ + fprintf(outfile, "# step | grid_x | grid_y | grid_z"); + for(unsigned int i = 0; i < header->field_count; i++) { + if(field_header[i].num_floats == 1) { + fprintf(outfile, " | %s", field_header[i].name); + } else { + for(unsigned int j = 0; j < field_header[i].num_floats; j++) { + fprintf(outfile, " | %s(%d)", field_header[i].name, j); + } + } + } + fprintf(outfile, "\n"); +} +void HemeLBExtractionFile::print_stats_column_headings(FILE *outfile) +{ + fprintf(outfile, "# step"); + for(unsigned int i = 0; i < header->field_count; i++) { + if(field_header[i].num_floats == 1) { + fprintf(outfile, " | %s [average | stdev | max | min]", field_header[i].name); + } else { + for(unsigned int j = 0; j < field_header[i].num_floats; j++) { + fprintf(outfile, " | %s(%d) [average | stdev | max | min]", field_header[i].name, j); + } + } + } + fprintf(outfile, "\n"); +} + +void HemeLBExtractionFile::print_all(FILE *outfile) +{ + snapshot->print(outfile); +} + +void HemeLBExtractionFile::print_stats(FILE *outfile) +{ + snapshot->print_stats(outfile); +} + +/** Reads a HemeLB extraction file header. Checks that the magic numbers are correct. */ +int HemeLBExtractionFile::read_header() +{ + uint32_t magic1 = 0, magic2 = 0; + uint32_t num_sites1, num_sites2; + + // Read heme magic number + xdr_u_int(&xdrs, &magic1); + if(magic1 != heme_magic) { + fprintf(stderr, "First uint32_t does not match heme magic number.\n"); + return 1; + } + + // TODO: Fix memory leak + header = new HEADER(); + + // Read extraction magic number + xdr_u_int(&xdrs, &magic2); + if(magic2 == extract_magic) { + if(bool_verbose == true) { + fprintf(stderr, "Reading a normal extraction file.\n"); + } + + header->is_colloids_file = false; + + xdr_u_int(&xdrs, &(header->version)); + xdr_double(&xdrs, &(header->voxelsz)); + xdr_double(&xdrs, &(header->originx)); + xdr_double(&xdrs, &(header->originy)); + xdr_double(&xdrs, &(header->originz)); + + // There is something wrong with either Heme or XDR which means that xdr_u_long is not working (possibly even reading more/less than 8 bytes...) + xdr_u_int(&xdrs, &num_sites1); + xdr_u_int(&xdrs, &num_sites2); + + // Combine the two ints into a long (assuming bigendian-ness) + header->num_sites = ((uint64_t)num_sites1)<<32 | num_sites2; + + xdr_u_int(&xdrs, &(header->field_count)); + xdr_u_int(&xdrs, &(header->field_header_length)); + + } else if (magic2 == colloid_magic) { + if(bool_verbose == true) { + fprintf(stderr, "Reading a colloids extraction file.\n"); + } + header->is_colloids_file = true; + + // Check version + uint32_t version=0; + xdr_u_int(&xdrs, &version); + if(bool_verbose == true) { + printf("Colloids version %u\n", version); + } + } else { + fprintf(stderr, "Unknown format: Second uint32_t does not match extraction file or colloids file magic number.\n"); + return 1; + } + + // Everything went well + return 0; +} + +/** Reads the field headers and remembers which columns correspond to velocity, pressure and shearstress. */ +int HemeLBExtractionFile::read_field_header() +{ + if(bool_verbose == true) { + printf("Reading field header...\n"); + } + int const maxlength = 100; // surely it can't be bigger than this...? + field_header = new FIELD_HEADER[header->field_count]; + for(unsigned int i = 0; i < header->field_count; i++) { + field_header[i].name = new char[maxlength]; + xdr_string(&xdrs, &(field_header[i].name), maxlength); + xdr_u_int(&xdrs, &(field_header[i].num_floats)); + xdr_double(&xdrs, &(field_header[i].offset)); + } + if(bool_verbose == true) { + printf("Finished reading field header...\n"); + } + + // Calculate number of 'columns' needed to represent a snapshot + uint32_t num_columns = 0; + for(unsigned int i = 0; i < header->field_count; i++) { + // If file contains velocity, pressure or shearstress, remember the corresponding column index for that field type. + if(strcmp(field_header[i].name, "velocity") == 0) { + has_velocity = true; + column_velocity = num_columns; + } else if(strcmp(field_header[i].name, "pressure") == 0) { + has_pressure = true; + column_pressure = num_columns; + } else if(strcmp(field_header[i].name, "shearstress") == 0) { + has_shearstress = true; + column_shearstress = num_columns; + } else if(strcmp(field_header[i].name, "d_shearstress") == 0) { // For "legacy" reasons only. Some old MCA runs have this field. + has_shearstress = true; + column_shearstress = num_columns; + } + + // This field has 'num_floats' columns to it + num_columns += field_header[i].num_floats; + } + header->num_columns = num_columns; + + // Everything went well + return 0; +} diff --git a/Mapping.h b/src/Mapping.cc similarity index 61% rename from Mapping.h rename to src/Mapping.cc index 6a951fd..7f863fe 100644 --- a/Mapping.h +++ b/src/Mapping.cc @@ -1,25 +1,21 @@ -#ifndef INCLUDED_MAPPING_H -#define INCLUDED_MAPPING_H +#include "Mapping.h" -#include -#include -#include -#include -#include -#include -#include +#include +#include +#include -#include "HemeLBExtractionFileTypes.h" -#include "HemeLBExtractionFile.h" - -void map(double old_voxelsz, double new_voxelsz, uint32_t in_coord, uint32_t *out_root, double *out_interpolate) +inline void map( + double old_voxelsz, + double new_voxelsz, + uint32_t in_coord, + uint32_t *out_root, + double *out_interpolate) { double real_pos = in_coord * old_voxelsz; - *out_root = floor(real_pos/new_voxelsz); + *out_root = std::floor(real_pos/new_voxelsz); *out_interpolate = (real_pos - *out_root * new_voxelsz) / new_voxelsz; } -/** Returns a mapping for each site in A to a corresponding trilinearly interpolated position in B. */ lattice_map * get_mapping(HemeLBExtractionFile *A, HemeLBExtractionFile *B) { // Work out mapping between the two lattices @@ -46,13 +42,13 @@ lattice_map * get_mapping(HemeLBExtractionFile *A, HemeLBExtractionFile *B) map(voxelA, voxelB, sitesA[i].z, &mapped_z, &mapping_A_B[i].c); // Add the 8 sites forming the cube this point lies within - mapped_sites[i * 8].set(mapped_x, mapped_y, mapped_z); // root (000) - mapped_sites[i * 8 + 1].set(mapped_x + 1, mapped_y, mapped_z); // (100) - mapped_sites[i * 8 + 2].set(mapped_x, mapped_y + 1, mapped_z); // (010) - mapped_sites[i * 8 + 3].set(mapped_x, mapped_y, mapped_z + 1); // (001) - mapped_sites[i * 8 + 4].set(mapped_x + 1, mapped_y, mapped_z + 1); // (101) - mapped_sites[i * 8 + 5].set(mapped_x, mapped_y + 1, mapped_z + 1); // (011) - mapped_sites[i * 8 + 6].set(mapped_x + 1, mapped_y + 1, mapped_z); // (110) + mapped_sites[i * 8 ].set(mapped_x , mapped_y , mapped_z ); // (000) root + mapped_sites[i * 8 + 1].set(mapped_x + 1, mapped_y , mapped_z ); // (100) + mapped_sites[i * 8 + 2].set(mapped_x , mapped_y + 1, mapped_z ); // (010) + mapped_sites[i * 8 + 3].set(mapped_x , mapped_y , mapped_z + 1); // (001) + mapped_sites[i * 8 + 4].set(mapped_x + 1, mapped_y , mapped_z + 1); // (101) + mapped_sites[i * 8 + 5].set(mapped_x , mapped_y + 1, mapped_z + 1); // (011) + mapped_sites[i * 8 + 6].set(mapped_x + 1, mapped_y + 1, mapped_z ); // (110) mapped_sites[i * 8 + 7].set(mapped_x + 1, mapped_y + 1, mapped_z + 1); // (111) } @@ -66,4 +62,3 @@ lattice_map * get_mapping(HemeLBExtractionFile *A, HemeLBExtractionFile *B) return mapping_A_B; } -#endif diff --git a/src/Snapshot.cc b/src/Snapshot.cc new file mode 100644 index 0000000..1e43e55 --- /dev/null +++ b/src/Snapshot.cc @@ -0,0 +1,210 @@ +#include "Snapshot.h" + +#include +#include +#include +#include +#include +#include + +// Calculate statistics (helper function) +inline void calc_stats( + const std::vector &records, + double &average, double &stdev, double &max, double &min) +{ + double sum{0}, sum2{0}; + max = -INFINITY; + min = INFINITY; + for(const double& rec : records) { + sum += rec; + sum2 += rec * rec; + if(max < rec) max = rec; + if(min > rec) min = rec; + } + average = sum/records.size(); + stdev = std::sqrt(sum2/records.size() - average*average); +} + +// +// Column (helper class) +// +class Column +{ + public: + Column(unsigned int num_records_, double offset_) + : records(num_records_), offset(offset_) {} + + double get(uint64_t index) + { + if(index >= records.size()) { + fprintf(stderr, "Error: index (%lu) exceeds number of records in column (%lu)\n", index, records.size()); + exit(1); + } + return records[index]; + } + void set_plus_offset(uint64_t index, float value) + { + if(index >= records.size()) { + fprintf(stderr, "Error: index (%lu) exceeds number of records in column (%lu)\n", index, records.size()); + exit(1); + } + records[index] = value + offset; + stats_calcd = false; + } + void print(FILE *outfile, uint64_t index) { + if(index >= records.size()) { + fprintf(stderr, "Error: index (%lu) exceeds number of records in column (%lu)\n", index, records.size()); + exit(1); + } + fprintf(outfile, "%e", records[index]); + } + void print_stats(FILE *outfile) { + if (!stats_calcd) { + calc_stats(records,average,stdev,max,min); + stats_calcd = true; + } + fprintf(outfile, "%e %e %e %e", average, stdev, max, min); + } + + private: + std::vector records; + double offset; + bool stats_calcd{false}; + double average; + double stdev; + double max, min; +}; + +Snapshot::Snapshot(HEADER *header, FIELD_HEADER *field_header) + : num_sites(header->num_sites), + num_columns(header->num_columns), + voxelsz(header->voxelsz) +{ + sites = new Site[num_sites]; + columns = new Column*[num_columns]; + unsigned int c = 0; + for(unsigned int i = 0; i < header->field_count; i++) { + for(unsigned int bro=0; bro < field_header[i].num_floats; bro++) { + columns[c] = new Column(num_sites, field_header[i].offset); + c++; + } + } +} + +void Snapshot::set_timestep(double timestep) +{ + this->timestep = timestep; +} +double Snapshot::get_timestep() +{ + return timestep; +} + +void Snapshot::site_set(unsigned int site_index, uint32_t x, uint32_t y, uint32_t z) +{ + sites[site_index].set(x, y, z); +} + +void Snapshot::column_set_plus_offset(uint32_t column_index, uint64_t site_index, double value) +{ + if(column_index >= num_columns) { + fprintf(stderr, "Error: column index (%u) exceeds number of columns (%u)\n", column_index, num_columns); + exit(1); + } + columns[column_index]->set_plus_offset(site_index, value); +} + +double Snapshot::get(uint32_t column_index, uint64_t site_index) +{ + return columns[column_index]->get(site_index); +} + +bool Snapshot::get_site_index(uint32_t a, uint32_t b, uint32_t c, uint64_t *site_index) +{ + for(unsigned int i = 0; i < num_sites; i++) { + if(sites[i].equals(a, b, c)) { + *site_index = i; + return true; + } + } + + // If no such site can be found, return false + return false; +} + +int Snapshot::get_site_index(Site *s, uint64_t *site_index) +{ + return get_site_index(s->x, s->y, s->z, site_index); +} + +Site * Snapshot::get_sites() +{ + return sites; +} + +/* Very inefficient calculation of site_index corresponding to the given site. Only suitable for small data sets. */ +SiteIndex * Snapshot::get_site_indices(Site *list, uint64_t list_size) +{ + SiteIndex *indices = new SiteIndex[list_size]; + for(uint64_t i = 0; i < list_size; i++) { + indices[i].exists = get_site_index(&list[i], &(indices[i].index)); + } + return indices; +} + +/** Builds hashtable to get indices corresponding to the given site. Necessary for very large data sets. */ +SiteIndex * Snapshot::get_site_indices_hashed_lookup(Site *list, uint64_t list_size, bool bool_verbose) +{ + // Build the hash table + if(bool_verbose == true) { fprintf(stderr, "# Building hashtable...\n"); } + std::unordered_map hashtable; + for(uint64_t i = 0; i < num_sites; i++) { + std::ostringstream oss; + oss << sites[i].x << "," << sites[i].y << "," << sites[i].z; + std::string hashstr = oss.str(); + hashtable[hashstr] = i; + } + if(bool_verbose == true) { fprintf(stderr, "# ...done building hashtable\n"); } + + if(bool_verbose == true) { fprintf(stderr, "# Calculating mapping...\n"); } + SiteIndex *indices = new SiteIndex[list_size]; + for(uint64_t i = 0; i < list_size; i++) { + std::ostringstream oss; + oss << list[i].x << "," << list[i].y << "," << list[i].z; + std::string hashstr = oss.str(); + if(hashtable.find(hashstr) != hashtable.end()) { + indices[i].exists = true; + indices[i].index = hashtable[hashstr]; + } else { + indices[i].exists = false; + } + } + if(bool_verbose == true) { fprintf(stderr, "# ... done calculating mapping\n"); } + + return indices; +} + +void Snapshot::print(FILE *outfile) +{ + for(uint64_t s = 0; s < num_sites; s++) { + fprintf(outfile, "%f ", timestep); + sites[s].print(outfile, this->voxelsz); + for(uint32_t c = 0; c < num_columns; c++) { + fprintf(outfile, " "); + columns[c]->print(outfile, s); + } + fprintf(outfile, "\n"); + } + fprintf(outfile, "\n"); +} + +void Snapshot::print_stats(FILE *outfile) +{ + fprintf(outfile, "%f", timestep); + for(uint32_t c = 0; c < num_columns; c++) { + fprintf(outfile, " "); + columns[c]->print_stats(outfile); + } + fprintf(outfile, "\n"); +} + diff --git a/hemeXtract.cc b/src/hemeXtract.cc similarity index 97% rename from hemeXtract.cc rename to src/hemeXtract.cc index cccbe21..2e29617 100644 --- a/hemeXtract.cc +++ b/src/hemeXtract.cc @@ -1,19 +1,14 @@ -#include -#include -#include #include -#include #include -#include #include +#include -#include -//#include #include -#include "Mapping.h" -#include "Compare.h" -#include "HemeLBExtractionFile.h" +#include "Vector3.h" +#include "HemeLBExtractionFile.h" // HemeLBExtractionFile +#include "Compare.h" // compare, diff +#include "Mapping.h" // get_mapping /* Command line parsing stuff */ const char *argp_program_version = "hemeXtract 0.0.1"; @@ -57,6 +52,7 @@ struct arguments { bool relativeErr; bool normalize_correl; }; + static Vector3 * parse_vector3(char* arg) { std::stringstream ss(arg); std::vector result; @@ -71,6 +67,7 @@ static Vector3 * parse_vector3(char* arg) { } return new Vector3(atof(result[0].c_str()), atof(result[1].c_str()), atof(result[2].c_str())); } + static error_t parse_opt(int key, char *arg, struct argp_state *state) { struct arguments *arggs = reinterpret_cast(state->input); @@ -136,10 +133,11 @@ static error_t parse_opt(int key, char *arg, struct argp_state *state) { return 0; } + static struct argp argp = { options, parse_opt, args_doc, doc, 0, 0, 0 }; /** Simple pointer swapper */ -void swap(HemeLBExtractionFile **i, HemeLBExtractionFile **j) { +inline void swap(HemeLBExtractionFile **i, HemeLBExtractionFile **j) { HemeLBExtractionFile *temp = *i; *i = *j; *j = temp; From 8aa4f4d6ec6a366491ea1ee8ab29ab47beecfad1 Mon Sep 17 00:00:00 2001 From: Ivan Pribec Date: Wed, 24 May 2023 03:49:09 +0200 Subject: [PATCH 16/16] Refactoring --- CMakeLists.txt | 2 +- include/Compare.h | 12 +-- include/HemeLBExtractionFile.h | 83 +++++++-------- include/Site.h | 12 ++- include/Snapshot.h | 21 ++-- include/Vector3.h | 14 +-- src/Compare.cc | 45 +++++--- src/HemeLBExtractionFile.cc | 187 ++++++++++++++++----------------- src/Mapping.cc | 32 +++--- src/Snapshot.cc | 57 +++++----- 10 files changed, 235 insertions(+), 230 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index af85fdc..db41814 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -4,7 +4,7 @@ project(hemeXtract VERSION 0.1.0 LANGUAGES CXX) -set(CMAKE_CXX_STANDARD 11) +set(CMAKE_CXX_STANDARD 14) set(CMAKE_MODULE_PATH "${CMAKE_SOURCE_DIR}/cmake") diff --git a/include/Compare.h b/include/Compare.h index e8c0458..c587863 100644 --- a/include/Compare.h +++ b/include/Compare.h @@ -9,17 +9,17 @@ /** Compares two lattices, A and B, with mapping mapA_to_B, by calculating their crosscorrelation, and their L2 norm difference. */ void compare(FILE *outfile, - lattice_map *mapA_to_B, - HemeLBExtractionFile *A, - HemeLBExtractionFile *B, + const lattice_map *mapA_to_B, + const HemeLBExtractionFile *A, + const HemeLBExtractionFile *B, int minexistent, bool normalize_correl); /** Calculates the difference between each site in A, and the corresponding (trilinearly interpolated) point in B */ void diff(FILE *outfile, - lattice_map *mapA_to_B, - HemeLBExtractionFile *A, - HemeLBExtractionFile *B, + const lattice_map *mapA_to_B, + const HemeLBExtractionFile *A, + const HemeLBExtractionFile *B, int minexistent, Vector3 *project, bool relativeErr, diff --git a/include/HemeLBExtractionFile.h b/include/HemeLBExtractionFile.h index 4f38ecd..5abe5ef 100644 --- a/include/HemeLBExtractionFile.h +++ b/include/HemeLBExtractionFile.h @@ -3,6 +3,7 @@ #include #include +#include // unique_ptr #include // must precede xdr.h #include // XDR @@ -28,83 +29,77 @@ class HemeLBExtractionFile int load_next_snapshot(); - // Workaround for XDR's broken long reading ability - void xdr_long(XDR *xdrs, uint64_t *ret); - void read_and_print_colloids(FILE *outfile); - bool correctly_initialised(); - - bool is_colloids_file(); + bool correctly_initialised() const; - bool no_more_snapshots(); + bool is_colloids_file() const; - uint64_t get_num_sites(); + bool no_more_snapshots() const; - double get_voxelsz(); + uint64_t get_num_sites() const; - bool hasVelocity(); + double get_voxelsz() const; - bool hasShearStress(); + bool hasVelocity() const; + bool hasShearStress() const ; + bool hasPressure() const; - bool hasPressure(); + double get_scaling() const; - double get_scaling(); - - double get_scalar_quantity(uint32_t column_index, uint64_t site_index); + double get_scalar_quantity(uint32_t column_index, uint64_t site_index) const; double get_interpolated_scalar_quantity( uint32_t column_index, - lattice_map *map); + const lattice_map &map) const; void get_vector_quantity( uint32_t column_index, uint64_t site_index, - Vector3 *returned_val); + Vector3 *returned_val) const; void get_interpolated_vector_quantity( uint32_t column_index, - lattice_map *map, - Vector3 *returned_val); + const lattice_map &map, + Vector3 *returned_val) const; - void get_pressure(uint64_t site_index, double *returned_val); + void get_pressure(uint64_t site_index, double *returned_val) const; - void get_velocity(uint64_t site_index, Vector3 *returned_val); + void get_velocity(uint64_t site_index, Vector3 *returned_val) const; - void get_shearstress(uint64_t site_index, double *returned_val); + void get_shearstress(uint64_t site_index, double *returned_val) const; - void get_interpolated_pressure(lattice_map *map, double *returned_val); + void get_interpolated_pressure(const lattice_map &map, double *returned_val) const; - void get_interpolated_velocity(lattice_map *map, Vector3 *returned_val); + void get_interpolated_velocity(const lattice_map &map, Vector3 *returned_val) const; - void get_interpolated_shearstress(lattice_map *map, double *returned_val); + void get_interpolated_shearstress(const lattice_map &map, double *returned_val) const; - Site * get_sites(); + const Site * get_sites() const; SiteIndex * get_site_indices(Site *list, uint64_t list_size); - SiteIndex * get_site_indices_hashed_lookup(Site *list, uint64_t list_size); - double get_time(); - double get_time_next(); + double get_time() const; + double get_time_next() const; - void print_header(FILE *outfile); - void print_field_header(FILE *outfile); - void print_column_headings(FILE *outfile); - void print_stats_column_headings(FILE *outfile); - void print_all(FILE *outfile); - void print_stats(FILE *outfile); + void print_header(FILE *outfile) const; + void print_field_header(FILE *outfile) const; + void print_column_headings(FILE *outfile) const; + void print_stats_column_headings(FILE *outfile) const; + void print_all(FILE *outfile) const; + void print_stats(FILE *outfile) const; private: /* True if file has been opened and the headers read without problem */ - bool bool_correctly_initialised; + bool bool_correctly_initialised{false}; /* True if user wants verbose output. If false, means quiet. */ bool bool_verbose; /* All info contained in the file header */ - HEADER *header; - FIELD_HEADER *field_header; + std::unique_ptr
header{nullptr}; + std::vector field_header; /** The step length for this file (has to be provided by the user since Heme extraction files don't store this value...) */ double step_length; @@ -113,7 +108,7 @@ class HemeLBExtractionFile double scaling; /** The currently loaded snapshot */ - Snapshot *snapshot; + std::unique_ptr snapshot{nullptr}; /** The open extraction file's stream */ FILE *in; @@ -128,7 +123,9 @@ class HemeLBExtractionFile bool bool_no_more_snapshots; /** True if the file contains these field types */ - bool has_velocity, has_shearstress, has_pressure; + bool has_velocity{false}; + bool has_shearstress{false}; + bool has_pressure{false}; /** The column indices in which to find each field type */ uint32_t column_velocity; @@ -136,9 +133,9 @@ class HemeLBExtractionFile uint32_t column_shearstress; /** The translation vector (in grid units) to be applied to this lattice (if requested by --translate option) */ - int32_t translate_grid_x; - int32_t translate_grid_y; - int32_t translate_grid_z; + int32_t translate_grid_x{0}; + int32_t translate_grid_y{0}; + int32_t translate_grid_z{0}; /** Reads a HemeLB extraction file header. Checks that the magic numbers are correct. */ int read_header(); diff --git a/include/Site.h b/include/Site.h index 49deda5..a3c8e49 100644 --- a/include/Site.h +++ b/include/Site.h @@ -9,6 +9,8 @@ class Site public: Site() {} + Site(uint32_t x_, uint32_t y_, uint32_t z_) + : x(x_), y(y_), z(z_) {} void set(uint32_t x, uint32_t y, uint32_t z) { this->x = x; @@ -16,15 +18,15 @@ class Site this->z = z; } - void print(FILE *outfile) { + void print(FILE *outfile) const { fprintf(outfile, "%u %u %u", x, y, z); } - void print(FILE *outfile, double voxelsz) { + void print(FILE *outfile, double voxelsz) const { fprintf(outfile, "%f %f %f", x * voxelsz, y * voxelsz, z * voxelsz); } - inline bool equals(uint32_t a, uint32_t b, uint32_t c) + inline bool equals(uint32_t a, uint32_t b, uint32_t c) const { if(a != x) { return false; @@ -38,9 +40,9 @@ class Site return true; } - bool equals(Site *s) + bool equals(const Site &s) const { - return equals(s->x, s->y, s->z); + return equals(s.x, s.y, s.z); } uint32_t x{0}, y{0}, z{0}; diff --git a/include/Snapshot.h b/include/Snapshot.h index fab81d8..a72882d 100644 --- a/include/Snapshot.h +++ b/include/Snapshot.h @@ -3,6 +3,7 @@ #include #include +#include #include "Site.h" // Site #include "HemeLBExtractionFileTypes.h" // HEADER, FIELD_HEADER, SiteIndex @@ -11,30 +12,28 @@ class Column; // Forward declaration (could be private potentially) class Snapshot { public: - Snapshot(HEADER *header, FIELD_HEADER *field_header); + Snapshot(const HEADER *header, const FIELD_HEADER *field_header); void set_timestep(double timestep); - double get_timestep(); + double get_timestep() const; void site_set(unsigned int site_index, uint32_t x, uint32_t y, uint32_t z); void column_set_plus_offset(uint32_t column_index, uint64_t site_index, double value); - double get(uint32_t column_index, uint64_t site_index); + double get(uint32_t column_index, uint64_t site_index) const; - bool get_site_index(uint32_t a, uint32_t b, uint32_t c, uint64_t *site_index); + bool get_site_index(const Site &s, uint64_t &site_index) const; - int get_site_index(Site *s, uint64_t *site_index); - - Site *get_sites(); + const Site *get_sites() const; /* Very inefficient calculation of site_index corresponding to the * given site. Only suitable for small data sets. */ - SiteIndex *get_site_indices(Site *list, uint64_t list_size); + SiteIndex *get_site_indices(const Site *list, uint64_t list_size); /* Builds hashtable to get indices corresponding to the given * site. Necessary for very large data sets. */ - SiteIndex *get_site_indices_hashed_lookup(Site *list, + SiteIndex *get_site_indices_hashed_lookup(const Site *list, uint64_t list_size, bool bool_verbose); void print(FILE *outfile); @@ -45,8 +44,8 @@ class Snapshot { uint64_t num_sites; uint32_t num_columns; double voxelsz; - Column **columns; - Site *sites; + Column **columns; // Memory Leak!! + std::vector sites; }; #endif diff --git a/include/Vector3.h b/include/Vector3.h index 0455a72..400560a 100644 --- a/include/Vector3.h +++ b/include/Vector3.h @@ -28,12 +28,12 @@ class Vector3 this->z *= s; } - double dot(Vector3 *v) + double dot(Vector3 *v) const { return this->x*v->x + this->y*v->y + this->z*v->z; } - double length() + double length() const { return sqrt(x*x + y*y + z*z); } @@ -54,7 +54,7 @@ class Vector3 z /= l; } - double abs_diff(Vector3 *v) + double abs_diff(Vector3 *v) const { double dx = x - v->x; double dy = y - v->y; @@ -69,22 +69,22 @@ class Vector3 this->z += v->z * scale; } - double get_x() + double get_x() const { return this->x; } - double get_y() + double get_y() const { return this->y; } - double get_z() + double get_z() const { return this->z; } - void print(FILE *outfile) + void print(FILE *outfile) const { fprintf(outfile, "(%f %f %f)\n", x, y, z); } diff --git a/src/Compare.cc b/src/Compare.cc index 7cfec33..7a7f9b1 100644 --- a/src/Compare.cc +++ b/src/Compare.cc @@ -3,7 +3,13 @@ #include /** Compares two lattices, A and B, with mapping mapA_to_B, by calculating their crosscorrelation, and their L2 norm difference. */ -void compare(FILE *outfile, lattice_map *mapA_to_B, HemeLBExtractionFile *A, HemeLBExtractionFile *B, int minexistent, bool normalize_correl) +void compare( + FILE *outfile, + const lattice_map *mapA_to_B, + const HemeLBExtractionFile *A, + const HemeLBExtractionFile *B, + int minexistent, + bool normalize_correl) { uint64_t num_sites_A = A->get_num_sites(); double timeA = A->get_time(); @@ -41,7 +47,7 @@ void compare(FILE *outfile, lattice_map *mapA_to_B, HemeLBExtractionFile *A, Hem // Velocity A->get_velocity(i, &velA); - B->get_interpolated_velocity(&mapA_to_B[i], &velB); + B->get_interpolated_velocity(mapA_to_B[i], &velB); // Max amplitude if(velA.length() > max_vel_A) { @@ -77,7 +83,7 @@ void compare(FILE *outfile, lattice_map *mapA_to_B, HemeLBExtractionFile *A, Hem for(uint64_t i = 0; i < num_sites_A; i++) { // Cross correlation shearstress A->get_shearstress(i, &shearstressA); - B->get_interpolated_shearstress(&mapA_to_B[i], &shearstressB); + B->get_interpolated_shearstress(mapA_to_B[i], &shearstressB); crosscorrelAB += shearstressA * shearstressB; // Autocorrelations @@ -87,8 +93,8 @@ void compare(FILE *outfile, lattice_map *mapA_to_B, HemeLBExtractionFile *A, Hem // L2 shearstressL2 += (shearstressA - shearstressB) * (shearstressA - shearstressB); } - correl_shearstress = crosscorrelAB / (sqrt(autocorrelA) * sqrt(autocorrelB)); - shearstressL2 = sqrt(shearstressL2); + correl_shearstress = crosscorrelAB / (std::sqrt(autocorrelA) * std::sqrt(autocorrelB)); + shearstressL2 = std::sqrt(shearstressL2); } if(do_pressure) { @@ -98,22 +104,30 @@ void compare(FILE *outfile, lattice_map *mapA_to_B, HemeLBExtractionFile *A, Hem for(uint64_t i = 0; i < num_sites_A; i++) { // Pressure A->get_pressure(i, &presA); - B->get_interpolated_pressure(&mapA_to_B[i], &presB); + B->get_interpolated_pressure(mapA_to_B[i], &presB); autocorrelA += presA * presA; autocorrelB += presB * presB; crosscorrelAB += presA * presB; presL2 += (presA - presB) * (presA - presB); } - correl_pres = crosscorrelAB / (sqrt(autocorrelA) * sqrt(autocorrelB)); - presL2 = sqrt(presL2); + correl_pres = crosscorrelAB / (std::sqrt(autocorrelA) * std::sqrt(autocorrelB)); + presL2 = std::sqrt(presL2); } fprintf(outfile, "%f %f %f %f %f %f %f %f %f\n", timeA, timeB, correl_vel, max_vel_A, max_vel_B, correl_shearstress, shearstressL2, correl_pres, presL2); } /** Calculates the difference between each site in A, and the corresponding (trilinearly interpolated) point in B */ -void diff(FILE *outfile, lattice_map *mapA_to_B, HemeLBExtractionFile *A, HemeLBExtractionFile *B, int minexistent, Vector3 *project, bool relativeErr, bool verbose) +void diff( + FILE *outfile, + const lattice_map *mapA_to_B, + const HemeLBExtractionFile *A, + const HemeLBExtractionFile *B, + int minexistent, + Vector3 *project, + bool relativeErr, + bool verbose) { uint64_t num_sites_A = A->get_num_sites(); double timeA = A->get_time(); @@ -125,7 +139,7 @@ void diff(FILE *outfile, lattice_map *mapA_to_B, HemeLBExtractionFile *A, HemeLB if(do_velocity) { - if(verbose == true) { + if(verbose) { fprintf(stderr, "Velocity difference calc\n"); fprintf(outfile, "# timeA=%f timeB=%f\n", timeA, timeB); } @@ -147,7 +161,7 @@ void diff(FILE *outfile, lattice_map *mapA_to_B, HemeLBExtractionFile *A, HemeLB // Get velocity from A and corresponding interpolated velocity from B A->get_velocity(i, &velA); - B->get_interpolated_velocity(&mapA_to_B[i], &velB); + B->get_interpolated_velocity(mapA_to_B[i], &velB); // fprintf(stderr, "velA "); // velA.print(stderr); @@ -180,7 +194,7 @@ void diff(FILE *outfile, lattice_map *mapA_to_B, HemeLBExtractionFile *A, HemeLB } if(do_shearstress) { - if(verbose == true) { + if(verbose) { fprintf(stderr, "WSS difference calc\n"); fprintf(outfile, "# timeA=%f timeB=%f\n", timeA, timeB); } @@ -202,15 +216,16 @@ void diff(FILE *outfile, lattice_map *mapA_to_B, HemeLBExtractionFile *A, HemeLB // Get shearstress from A and corresponding interpolated shearstress from B A->get_shearstress(i, &shearA); - B->get_interpolated_shearstress(&mapA_to_B[i], &shearB); + B->get_interpolated_shearstress(mapA_to_B[i], &shearB); // Print out site coords - A->get_sites()[i].print(outfile, voxelA); + auto sites = A->get_sites(); + sites[i].print(outfile, voxelA); double sheardiff = shearA - shearB; // If requested, give difference relative to local WSS from input file A - if(relativeErr == true) { + if (relativeErr) { sheardiff /= shearA; } diff --git a/src/HemeLBExtractionFile.cc b/src/HemeLBExtractionFile.cc index b3e5efd..51b571a 100644 --- a/src/HemeLBExtractionFile.cc +++ b/src/HemeLBExtractionFile.cc @@ -2,7 +2,6 @@ #include - /* Magic numbers designating: * - a heme file, * - a heme extraction file, @@ -12,18 +11,24 @@ static const uint32_t heme_magic = 0x686C6221; // hlb! static const uint32_t extract_magic = 0x78747204; // xtr static const uint32_t colloid_magic = 0x636F6C04; // colloid + +// Workaround for XDR's broken long reading ability +void xdr_long(XDR *xdrs, uint64_t *ret) +{ + uint32_t half1, half2; + xdr_u_int(xdrs, &half1); + xdr_u_int(xdrs, &half2); + *ret = ((uint64_t)half1)<<32 | half2; // Assume big-endianness +} + HemeLBExtractionFile::HemeLBExtractionFile( char *fname, - double step_length, - double scaling, + double step_length_, + double scaling_, bool verbose, Vector3 *translate) + : step_length(step_length_), scaling(scaling_), bool_verbose(verbose) { - this->step_length = step_length; - this->scaling = scaling; - this->bool_verbose = verbose; - bool_correctly_initialised = false; - snapshot = NULL; in = fopen(fname, "rb"); if (in == NULL) { fprintf(stderr, "Could not open file '%s'. Does it exist?\n", fname); @@ -31,10 +36,6 @@ HemeLBExtractionFile::HemeLBExtractionFile( } xdrstdio_create(&xdrs, in, XDR_DECODE); - has_velocity = false; - has_pressure = false; - has_shearstress = false; - // Read the file header, populating the 'header' class variable. Determines whether file is "normal" or colloids int r; r = read_header(); @@ -44,7 +45,7 @@ HemeLBExtractionFile::HemeLBExtractionFile( } // If file is a colloids file, return now (no more info to be read in) - if(header->is_colloids_file == true) { + if (header->is_colloids_file) { bool_correctly_initialised = true; return; } @@ -57,21 +58,17 @@ HemeLBExtractionFile::HemeLBExtractionFile( } // If specified, convert lattice translation vector (real units) into grid units (via rounding based on voxel size) - if(translate != NULL) { - this->translate_grid_x = (int32_t)(translate->get_x()/header->voxelsz); - this->translate_grid_y = (int32_t)(translate->get_y()/header->voxelsz); - this->translate_grid_z = (int32_t)(translate->get_z()/header->voxelsz); + if(translate != nullptr) { + translate_grid_x = (int32_t)(translate->get_x()/header->voxelsz); + translate_grid_y = (int32_t)(translate->get_y()/header->voxelsz); + translate_grid_z = (int32_t)(translate->get_z()/header->voxelsz); if(bool_verbose == true) { fprintf(stderr, "Lattice translation vector (%e,%e,%e) in real units, converted to (%d,%d,%d) in grid units (voxel size = %e).\n", translate->get_x(), translate->get_y(), translate->get_z(), translate_grid_x, translate_grid_y, translate_grid_z, header->voxelsz); } - } else { - this->translate_grid_x = 0; - this->translate_grid_y = 0; - this->translate_grid_z = 0; } // Allocate a snapshot of the right dimensions - snapshot = new Snapshot(header, field_header); + snapshot = std::make_unique(header.get(), field_header.data()); // Get the time of the next snapshot bool next = read_time_next(); @@ -94,18 +91,14 @@ bool HemeLBExtractionFile::read_time_next() { // Check if we have reached the end of the file if (feof(in)) { - if(bool_verbose == true) { + if(bool_verbose) { fprintf(stderr, "# Reached end of file.\n"); } return false; } - // Work around for problem with xdr_u_long() not working - uint32_t dt1, dt2; uint64_t step; - xdr_u_int(&xdrs, &dt1); - xdr_u_int(&xdrs, &dt2); - step = ((uint64_t)dt1)<<32 | dt2; + xdr_long(&xdrs,&step); time_next = step * step_length; return true; @@ -113,8 +106,8 @@ bool HemeLBExtractionFile::read_time_next() int HemeLBExtractionFile::load_next_snapshot() { - if(bool_no_more_snapshots == true) { - if(bool_verbose == true) { + if(bool_no_more_snapshots) { + if(bool_verbose) { fprintf(stderr, "Warning: No more snapshots left in file.\n"); } return 1; @@ -142,22 +135,17 @@ int HemeLBExtractionFile::load_next_snapshot() snapshot->column_set_plus_offset(i, s, value * this->scaling); } } - bool next = read_time_next(); - if(next == false) { + + if(!read_time_next()) { bool_no_more_snapshots = true; - if(bool_verbose == true) {fprintf(stderr, "No more snapshots. (Note that final snapshot may have been incomplete)\n");} + if(bool_verbose) { + fprintf(stderr, "No more snapshots. (Note that final snapshot may have been incomplete)\n"); + } } return 0; } -// Workaround for XDR's broken long reading ability -void HemeLBExtractionFile::xdr_long(XDR *xdrs, uint64_t *ret) -{ - uint32_t half1, half2; - xdr_u_int(xdrs, &half1); - xdr_u_int(xdrs, &half2); - *ret = ((uint64_t)half1)<<32 | half2; // Assume big-endianness -} + void HemeLBExtractionFile::read_and_print_colloids(FILE *outfile) @@ -182,15 +170,15 @@ void HemeLBExtractionFile::read_and_print_colloids(FILE *outfile) xdr_long(&xdrs, &timeStep); // Check for end of file here (to avoid spurious EOF output) - if(feof(in)) { + if (feof(in)) { break; } - if(recordLen <= 0) { + if (recordLen <= 0) { break; } - if(bool_verbose == true) { + if (bool_verbose) { fprintf(outfile, "# headerLen: %u recordLen %u dsetLen %ld timeStep: %ld\n", headerLen, recordLen, dsetLen, timeStep); } @@ -217,98 +205,100 @@ void HemeLBExtractionFile::read_and_print_colloids(FILE *outfile) fprintf(outfile, "%ld %ld %ld %.13e %.13e %.13e %.13e %.13e %.13e\n", timeStep, id, rank, X, Y, Z, VX, VY, VZ); } } - if(bool_verbose == true) { + if (bool_verbose) { fprintf(stderr, "# Reached end of file.\n"); } } -bool HemeLBExtractionFile::correctly_initialised() +bool HemeLBExtractionFile::correctly_initialised() const { return bool_correctly_initialised; } -bool HemeLBExtractionFile::is_colloids_file() +bool HemeLBExtractionFile::is_colloids_file() const { return header->is_colloids_file; } -bool HemeLBExtractionFile::no_more_snapshots() { +bool HemeLBExtractionFile::no_more_snapshots() const +{ return bool_no_more_snapshots; } -uint64_t HemeLBExtractionFile::get_num_sites() +uint64_t HemeLBExtractionFile::get_num_sites() const { return header->num_sites; } -double HemeLBExtractionFile::get_voxelsz() +double HemeLBExtractionFile::get_voxelsz() const { return header->voxelsz; } -bool HemeLBExtractionFile::hasVelocity() +bool HemeLBExtractionFile::hasVelocity() const { return has_velocity; } -bool HemeLBExtractionFile::hasShearStress() +bool HemeLBExtractionFile::hasShearStress() const { return has_shearstress; } -bool HemeLBExtractionFile::hasPressure() +bool HemeLBExtractionFile::hasPressure() const { return has_pressure; } -double HemeLBExtractionFile::get_scaling() +double HemeLBExtractionFile::get_scaling() const { return scaling; } -double HemeLBExtractionFile::get_scalar_quantity(uint32_t column_index, uint64_t site_index) +double HemeLBExtractionFile::get_scalar_quantity(uint32_t column_index, uint64_t site_index) const { return snapshot->get(column_index, site_index); } -double HemeLBExtractionFile::get_interpolated_scalar_quantity(uint32_t column_index, lattice_map *map) +double HemeLBExtractionFile::get_interpolated_scalar_quantity(uint32_t column_index, const lattice_map &map) const { double average_existing = 0; - uint32_t num_existing = 0; - for(uint32_t i = 0; i < 8; i++) { - if(map->index[i].exists == true) { - average_existing += get_scalar_quantity(column_index, map->index[i].index); + int num_existing = 0; + for (int i = 0; i < 8; i++) { + if(map.index[i].exists) { + average_existing += get_scalar_quantity(column_index, map.index[i].index); num_existing++; } } - if(num_existing == 0) { + if (num_existing == 0) { return 0; } average_existing /= num_existing; // Get the scalars at the 8 lattice sites forming the cube - double s[8]; - for(uint32_t i = 0; i < 8; i++) { + double s[8]; + for (int i = 0; i < 8; i++) { // If there was no corresponding site_index, the site was not measured, so set it to the average // value of the existing sites in the map - if(map->index[i].exists == false) { + if(map.index[i].exists == false) { s[i] = average_existing; } else { - s[i] = get_scalar_quantity(column_index, map->index[i].index); + s[i] = get_scalar_quantity(column_index, map.index[i].index); } } // Trilinear interpolation within the cube - return s[0] * (1 - map->a) * (1 - map->b) * (1 - map->c) - + s[1] * (map->a * (1 - map->b) * (1 - map->c)) - + s[2] * ((1 - map->a) * map->b * (1 - map->c)) - + s[3] * ((1 - map->a) * (1 - map->b) * map->c) - + s[4] * (map->a * (1 - map->b) * map->c) - + s[5] * ((1 - map->a) * map->b * map->c) - + s[6] * (map->a * map->b * (1 - map->c)) - + s[7] * (map->a * map->b * map->c); + return s[0] * (1 - map.a) * (1 - map.b) * (1 - map.c) + + s[1] * (map.a * (1 - map.b) * (1 - map.c)) + + s[2] * ((1 - map.a) * map.b * (1 - map.c)) + + s[3] * ((1 - map.a) * (1 - map.b) * map.c) + + s[4] * (map.a * (1 - map.b) * map.c) + + s[5] * ((1 - map.a) * map.b * map.c) + + s[6] * (map.a * map.b * (1 - map.c)) + + s[7] * (map.a * map.b * map.c); } -void HemeLBExtractionFile::get_vector_quantity(uint32_t column_index, uint64_t site_index, Vector3 *returned_val) + +void HemeLBExtractionFile::get_vector_quantity(uint32_t column_index, uint64_t site_index, Vector3 *returned_val) const { double vx, vy, vz; vx = get_scalar_quantity(column_index, site_index); @@ -317,7 +307,7 @@ void HemeLBExtractionFile::get_vector_quantity(uint32_t column_index, uint64_t s returned_val->set(vx, vy, vz); } -void HemeLBExtractionFile::get_interpolated_vector_quantity(uint32_t column_index, lattice_map *map, Vector3 *returned_val) +void HemeLBExtractionFile::get_interpolated_vector_quantity(uint32_t column_index, const lattice_map &map, Vector3 *returned_val) const { double vx, vy, vz; vx = get_interpolated_scalar_quantity(column_index, map); @@ -326,37 +316,37 @@ void HemeLBExtractionFile::get_interpolated_vector_quantity(uint32_t column_inde returned_val->set(vx, vy, vz); } -void HemeLBExtractionFile::get_pressure(uint64_t site_index, double *returned_val) +void HemeLBExtractionFile::get_pressure(uint64_t site_index, double *returned_val) const { *returned_val = get_scalar_quantity(column_pressure, site_index); } -void HemeLBExtractionFile::get_velocity(uint64_t site_index, Vector3 *returned_val) +void HemeLBExtractionFile::get_velocity(uint64_t site_index, Vector3 *returned_val) const { get_vector_quantity(column_velocity, site_index, returned_val); } -void HemeLBExtractionFile::get_shearstress(uint64_t site_index, double *returned_val) +void HemeLBExtractionFile::get_shearstress(uint64_t site_index, double *returned_val) const { *returned_val = get_scalar_quantity(column_shearstress, site_index); } -void HemeLBExtractionFile::get_interpolated_pressure(lattice_map *map, double *returned_val) +void HemeLBExtractionFile::get_interpolated_pressure(const lattice_map &map, double *returned_val) const { *returned_val = get_interpolated_scalar_quantity(column_pressure, map); } -void HemeLBExtractionFile::get_interpolated_velocity(lattice_map *map, Vector3 *returned_val) +void HemeLBExtractionFile::get_interpolated_velocity(const lattice_map &map, Vector3 *returned_val) const { get_interpolated_vector_quantity(column_velocity, map, returned_val); } -void HemeLBExtractionFile::get_interpolated_shearstress(lattice_map *map, double *returned_val) +void HemeLBExtractionFile::get_interpolated_shearstress(const lattice_map &map, double *returned_val) const { *returned_val = get_interpolated_scalar_quantity(column_shearstress, map); } -Site * HemeLBExtractionFile::get_sites() +const Site * HemeLBExtractionFile::get_sites() const { return snapshot->get_sites(); } @@ -371,17 +361,17 @@ SiteIndex * HemeLBExtractionFile::get_site_indices_hashed_lookup(Site *list, uin return snapshot->get_site_indices_hashed_lookup(list, list_size, bool_verbose); } -double HemeLBExtractionFile::get_time() +double HemeLBExtractionFile::get_time() const { return snapshot->get_timestep(); } -double HemeLBExtractionFile::get_time_next() +double HemeLBExtractionFile::get_time_next() const { return time_next; } -void HemeLBExtractionFile::print_header(FILE *outfile) +void HemeLBExtractionFile::print_header(FILE *outfile) const { fprintf(outfile, "\n# HEADER:\n# -------\n"); fprintf(outfile, "# version=%u\n", header->version); @@ -394,7 +384,7 @@ void HemeLBExtractionFile::print_header(FILE *outfile) fprintf(outfile, "# field_header_length=%u\n", header->field_header_length); } -void HemeLBExtractionFile::print_field_header(FILE *outfile) +void HemeLBExtractionFile::print_field_header(FILE *outfile) const { fprintf(outfile, "\n# FIELD HEADERS:\n# -------\n"); for(unsigned int i = 0; i < header->field_count; i++) { @@ -403,7 +393,7 @@ void HemeLBExtractionFile::print_field_header(FILE *outfile) fprintf(outfile, "# Number of 'columns' per snapshot = %d\n", header->num_columns); } -void HemeLBExtractionFile::print_column_headings(FILE *outfile) +void HemeLBExtractionFile::print_column_headings(FILE *outfile) const { fprintf(outfile, "# step | grid_x | grid_y | grid_z"); for(unsigned int i = 0; i < header->field_count; i++) { @@ -417,7 +407,7 @@ void HemeLBExtractionFile::print_column_headings(FILE *outfile) } fprintf(outfile, "\n"); } -void HemeLBExtractionFile::print_stats_column_headings(FILE *outfile) +void HemeLBExtractionFile::print_stats_column_headings(FILE *outfile) const { fprintf(outfile, "# step"); for(unsigned int i = 0; i < header->field_count; i++) { @@ -432,12 +422,12 @@ void HemeLBExtractionFile::print_stats_column_headings(FILE *outfile) fprintf(outfile, "\n"); } -void HemeLBExtractionFile::print_all(FILE *outfile) +void HemeLBExtractionFile::print_all(FILE *outfile) const { snapshot->print(outfile); } -void HemeLBExtractionFile::print_stats(FILE *outfile) +void HemeLBExtractionFile::print_stats(FILE *outfile) const { snapshot->print_stats(outfile); } @@ -455,13 +445,12 @@ int HemeLBExtractionFile::read_header() return 1; } - // TODO: Fix memory leak - header = new HEADER(); + header = std::make_unique
(); // Read extraction magic number xdr_u_int(&xdrs, &magic2); if(magic2 == extract_magic) { - if(bool_verbose == true) { + if(bool_verbose) { fprintf(stderr, "Reading a normal extraction file.\n"); } @@ -484,7 +473,7 @@ int HemeLBExtractionFile::read_header() xdr_u_int(&xdrs, &(header->field_header_length)); } else if (magic2 == colloid_magic) { - if(bool_verbose == true) { + if(bool_verbose) { fprintf(stderr, "Reading a colloids extraction file.\n"); } header->is_colloids_file = true; @@ -492,7 +481,7 @@ int HemeLBExtractionFile::read_header() // Check version uint32_t version=0; xdr_u_int(&xdrs, &version); - if(bool_verbose == true) { + if(bool_verbose) { printf("Colloids version %u\n", version); } } else { @@ -507,18 +496,20 @@ int HemeLBExtractionFile::read_header() /** Reads the field headers and remembers which columns correspond to velocity, pressure and shearstress. */ int HemeLBExtractionFile::read_field_header() { - if(bool_verbose == true) { + if(bool_verbose) { printf("Reading field header...\n"); } int const maxlength = 100; // surely it can't be bigger than this...? - field_header = new FIELD_HEADER[header->field_count]; + + field_header.resize(header->field_count); + for(unsigned int i = 0; i < header->field_count; i++) { field_header[i].name = new char[maxlength]; xdr_string(&xdrs, &(field_header[i].name), maxlength); xdr_u_int(&xdrs, &(field_header[i].num_floats)); xdr_double(&xdrs, &(field_header[i].offset)); } - if(bool_verbose == true) { + if(bool_verbose) { printf("Finished reading field header...\n"); } diff --git a/src/Mapping.cc b/src/Mapping.cc index 7f863fe..78c5fc0 100644 --- a/src/Mapping.cc +++ b/src/Mapping.cc @@ -8,19 +8,19 @@ inline void map( double old_voxelsz, double new_voxelsz, uint32_t in_coord, - uint32_t *out_root, - double *out_interpolate) + uint32_t &out_root, + double &out_interpolate) { double real_pos = in_coord * old_voxelsz; - *out_root = std::floor(real_pos/new_voxelsz); - *out_interpolate = (real_pos - *out_root * new_voxelsz) / new_voxelsz; + out_root = std::floor(real_pos/new_voxelsz); + out_interpolate = (real_pos - out_root * new_voxelsz) / new_voxelsz; } lattice_map * get_mapping(HemeLBExtractionFile *A, HemeLBExtractionFile *B) { // Work out mapping between the two lattices uint64_t num_sites_A = A->get_num_sites(); - Site *sitesA = A->get_sites(); + const Site *sitesA = A->get_sites(); double voxelA = A->get_voxelsz(); double voxelB = B->get_voxelsz(); uint32_t mapped_x, mapped_y, mapped_z; @@ -29,17 +29,18 @@ lattice_map * get_mapping(HemeLBExtractionFile *A, HemeLBExtractionFile *B) fprintf(stderr, "Could not allocate memory for lattice_map\n"); return NULL; } - Site *mapped_sites = new Site[num_sites_A * 8]; - if(mapped_sites == NULL) { - fprintf(stderr, "Could not allocate memory for mapped_sites\n"); - return NULL; - } + //Site *mapped_sites = new Site[num_sites_A * 8]; + //if(mapped_sites == NULL) { + // fprintf(stderr, "Could not allocate memory for mapped_sites\n"); + // return NULL; + //} + std::vector mapped_sites(num_sites_A * 8); for(uint64_t i = 0; i < num_sites_A; i++) { // Get the interpolation constants for this site in terms of the other lattice - map(voxelA, voxelB, sitesA[i].x, &mapped_x, &mapping_A_B[i].a); - map(voxelA, voxelB, sitesA[i].y, &mapped_y, &mapping_A_B[i].b); - map(voxelA, voxelB, sitesA[i].z, &mapped_z, &mapping_A_B[i].c); + map(voxelA, voxelB, sitesA[i].x, mapped_x, mapping_A_B[i].a); + map(voxelA, voxelB, sitesA[i].y, mapped_y, mapping_A_B[i].b); + map(voxelA, voxelB, sitesA[i].z, mapped_z, mapping_A_B[i].c); // Add the 8 sites forming the cube this point lies within mapped_sites[i * 8 ].set(mapped_x , mapped_y , mapped_z ); // (000) root @@ -52,10 +53,11 @@ lattice_map * get_mapping(HemeLBExtractionFile *A, HemeLBExtractionFile *B) mapped_sites[i * 8 + 7].set(mapped_x + 1, mapped_y + 1, mapped_z + 1); // (111) } - SiteIndex *site_indices = B->get_site_indices_hashed_lookup(mapped_sites, num_sites_A * 8); + SiteIndex *site_indices = B->get_site_indices_hashed_lookup( + mapped_sites.data(), num_sites_A * 8); for(uint64_t i = 0; i < num_sites_A; i++) { - for(uint32_t j = 0; j < 8; j++) { + for(int j = 0; j < 8; j++) { mapping_A_B[i].index[j] = site_indices[i * 8 + j]; } } diff --git a/src/Snapshot.cc b/src/Snapshot.cc index 1e43e55..521794c 100644 --- a/src/Snapshot.cc +++ b/src/Snapshot.cc @@ -6,6 +6,8 @@ #include #include #include +#include +#include // Calculate statistics (helper function) inline void calc_stats( @@ -34,7 +36,7 @@ class Column Column(unsigned int num_records_, double offset_) : records(num_records_), offset(offset_) {} - double get(uint64_t index) + double get(uint64_t index) const { if(index >= records.size()) { fprintf(stderr, "Error: index (%lu) exceeds number of records in column (%lu)\n", index, records.size()); @@ -75,12 +77,12 @@ class Column double max, min; }; -Snapshot::Snapshot(HEADER *header, FIELD_HEADER *field_header) +Snapshot::Snapshot(const HEADER *header, const FIELD_HEADER *field_header) : num_sites(header->num_sites), num_columns(header->num_columns), - voxelsz(header->voxelsz) + voxelsz(header->voxelsz), + sites(num_sites) { - sites = new Site[num_sites]; columns = new Column*[num_columns]; unsigned int c = 0; for(unsigned int i = 0; i < header->field_count; i++) { @@ -95,7 +97,7 @@ void Snapshot::set_timestep(double timestep) { this->timestep = timestep; } -double Snapshot::get_timestep() +double Snapshot::get_timestep() const { return timestep; } @@ -114,52 +116,50 @@ void Snapshot::column_set_plus_offset(uint32_t column_index, uint64_t site_index columns[column_index]->set_plus_offset(site_index, value); } -double Snapshot::get(uint32_t column_index, uint64_t site_index) +double Snapshot::get(uint32_t column_index, uint64_t site_index) const { return columns[column_index]->get(site_index); } -bool Snapshot::get_site_index(uint32_t a, uint32_t b, uint32_t c, uint64_t *site_index) +bool Snapshot::get_site_index(const Site &s, uint64_t &site_index) const { - for(unsigned int i = 0; i < num_sites; i++) { - if(sites[i].equals(a, b, c)) { - *site_index = i; - return true; - } + auto res = std::find_if(sites.begin(),sites.end(), + [=](Site st){ return st.equals(s); }); + + if (res != sites.end()) { + site_index = std::distance(sites.begin(),res); + return true; + } else { + // If no such site can be found, return false + return false; } - - // If no such site can be found, return false - return false; -} - -int Snapshot::get_site_index(Site *s, uint64_t *site_index) -{ - return get_site_index(s->x, s->y, s->z, site_index); } -Site * Snapshot::get_sites() +const Site *Snapshot::get_sites() const { - return sites; + return sites.data(); } -/* Very inefficient calculation of site_index corresponding to the given site. Only suitable for small data sets. */ -SiteIndex * Snapshot::get_site_indices(Site *list, uint64_t list_size) +/* Very inefficient calculation of site_index corresponding to the given site. + * Only suitable for small data sets. */ +SiteIndex * Snapshot::get_site_indices(const Site *list, uint64_t list_size) { SiteIndex *indices = new SiteIndex[list_size]; for(uint64_t i = 0; i < list_size; i++) { - indices[i].exists = get_site_index(&list[i], &(indices[i].index)); + indices[i].exists = get_site_index(list[i], indices[i].index); } return indices; } -/** Builds hashtable to get indices corresponding to the given site. Necessary for very large data sets. */ -SiteIndex * Snapshot::get_site_indices_hashed_lookup(Site *list, uint64_t list_size, bool bool_verbose) +/* Builds hashtable to get indices corresponding to the given site. + * Necessary for very large data sets. */ +SiteIndex * Snapshot::get_site_indices_hashed_lookup(const Site *list, uint64_t list_size, bool bool_verbose) { // Build the hash table if(bool_verbose == true) { fprintf(stderr, "# Building hashtable...\n"); } std::unordered_map hashtable; + std::ostringstream oss; for(uint64_t i = 0; i < num_sites; i++) { - std::ostringstream oss; oss << sites[i].x << "," << sites[i].y << "," << sites[i].z; std::string hashstr = oss.str(); hashtable[hashstr] = i; @@ -169,7 +169,6 @@ SiteIndex * Snapshot::get_site_indices_hashed_lookup(Site *list, uint64_t list_s if(bool_verbose == true) { fprintf(stderr, "# Calculating mapping...\n"); } SiteIndex *indices = new SiteIndex[list_size]; for(uint64_t i = 0; i < list_size; i++) { - std::ostringstream oss; oss << list[i].x << "," << list[i].y << "," << list[i].z; std::string hashstr = oss.str(); if(hashtable.find(hashstr) != hashtable.end()) {