diff --git a/advection/.gitignore b/advection/.gitignore deleted file mode 100644 index 34ddd4b..0000000 --- a/advection/.gitignore +++ /dev/null @@ -1,7 +0,0 @@ -.DS_Store -src/.DS_Store -src/.cache -src/build -.idea -src/cmake-build-debug -src/cmake-build-release diff --git a/advection/README.md b/advection/README.md deleted file mode 100644 index 0abcf14..0000000 --- a/advection/README.md +++ /dev/null @@ -1,46 +0,0 @@ -## What is new? -There is one new added component: `AdvectionKernel`s which is an "interface" (i.e an abstract class). -There are two implementations simple Euler integration called `EulerIntegrationKernel` and -Runge Kutta integration called `RK4AdvectionKernel`. - -Main function gives a good example of how to use the library. Especially the following function which prints the -position of the particle at every time step. -```Cpp -template -void advectForSomeTime(const UVGrid &uvGrid, const AdvectionKernelImpl &kernel, double latstart, double lonstart) { - - // Require at compile time that kernel derives from the abstract class AdvectionKernel - static_assert(std::is_base_of::value, NotAKernelError); - - double lat1 = latstart, lon1 = lonstart; - for(int time = 100; time <= 10000; time += AdvectionKernel::DT) { - cout << "lat = " << lat1 << " lon = " << lon1 << endl; - auto [templat, templon] = kernel.advect(time, lat1, lon1); - lat1 = templat; - lon1 = templon; - } -} -``` - - -## Location of data -The data path is hardcoded such that the following tree structure is assumed: -The current assumption is that the name of the `u`s and `v`s are flipped since this is the way the data was given to us. -``` -data/ - grid.h5 - hydrodynamic_U.h5 - hydrodynamic_V.h5 -interactive-track-and-trace/ - opening-hdf5/ - ... -``` - -## Compiling -Let the current directory be the `src` directory. Run: -```shell -mkdir build -cd build -cmake .. -make -``` \ No newline at end of file diff --git a/advection/src/CMakeLists.txt b/advection/src/CMakeLists.txt deleted file mode 100644 index b8c8758..0000000 --- a/advection/src/CMakeLists.txt +++ /dev/null @@ -1,42 +0,0 @@ -cmake_minimum_required (VERSION 3.28) -project (Advection) - -set(CMAKE_CXX_STANDARD 23) -set(CMAKE_CXX_STANDARD_REQUIRED ON) - -set(CMAKE_EXPORT_COMPILE_COMMANDS ON) - -find_package(netCDF REQUIRED) - -add_executable(Advection main.cpp - readdata.cpp - readdata.h - interpolate.cpp - interpolate.h - UVGrid.cpp - UVGrid.h - Vel.h - Vel.cpp - AdvectionKernel.h - EulerAdvectionKernel.cpp - EulerAdvectionKernel.h - RK4AdvectionKernel.cpp - RK4AdvectionKernel.h -) - -execute_process( - COMMAND nc-config --includedir - OUTPUT_VARIABLE NETCDF_INCLUDE_DIR - OUTPUT_STRIP_TRAILING_WHITESPACE -) - -execute_process( - COMMAND ncxx4-config --libdir - OUTPUT_VARIABLE NETCDFCXX_LIB_DIR - OUTPUT_STRIP_TRAILING_WHITESPACE -) - -target_include_directories(Advection PUBLIC ${netCDF_INCLUDE_DIR}) - -find_library(NETCDF_LIB NAMES netcdf-cxx4 netcdf_c++4 PATHS ${NETCDFCXX_LIB_DIR} NO_DEFAULT_PATH) -target_link_libraries(Advection ${NETCDF_LIB}) diff --git a/advection/src/EulerAdvectionKernel.cpp b/advection/src/EulerAdvectionKernel.cpp deleted file mode 100644 index 1a38944..0000000 --- a/advection/src/EulerAdvectionKernel.cpp +++ /dev/null @@ -1,13 +0,0 @@ - -#include "EulerAdvectionKernel.h" -#include "interpolate.h" - -using namespace std; - -EulerAdvectionKernel::EulerAdvectionKernel(std::shared_ptr grid): grid(grid) { } - -std::pair EulerAdvectionKernel::advect(int time, double latitude, double longitude) const { - auto [u, v] = bilinearinterpolate(*grid, time, latitude, longitude); - - return {latitude+metreToDegrees(v*DT), longitude+metreToDegrees(u*DT)}; -} diff --git a/advection/src/UVGrid.cpp b/advection/src/UVGrid.cpp deleted file mode 100644 index 1febbfe..0000000 --- a/advection/src/UVGrid.cpp +++ /dev/null @@ -1,66 +0,0 @@ -#include - -#include "UVGrid.h" -#include "readdata.h" - -#define sizeError2 "The sizes of the hydrodynamic data files are different" -#define sizeError "The sizes of the hydrodynamicU or -V files does not correspond with the sizes of the grid file" - -using namespace std; - -UVGrid::UVGrid() { - auto us = readHydrodynamicU(); - auto vs = readHydrodynamicV(); - if (us.size() != vs.size()) { - throw domain_error(sizeError2); - } - - tie(times, lats, lons) = readGrid(); - - timeSize = times.size(); - latSize = lats.size(); - lonSize = lons.size(); - - size_t gridSize = timeSize * latSize * lonSize; - if (gridSize != us.size()) { - throw domain_error(sizeError); - } - - uvData.reserve(gridSize); - - for (auto vel: views::zip(us, vs)) { - uvData.push_back(Vel(vel)); - } -} - -const Vel &UVGrid::operator[](size_t timeIndex, size_t latIndex, size_t lonIndex) const { - if(timeIndex < 0 or timeIndex >= timeSize - or latIndex < 0 or latIndex >= latSize - or lonIndex < 0 or lonIndex >= lonSize) { - throw std::out_of_range("Index out of bounds"); - } - size_t index = timeIndex * (latSize * lonSize) + latIndex * lonSize + lonIndex; - return uvData[index]; -} - -double UVGrid::lonStep() const { - return lons[1] - lons[0]; -} - -double UVGrid::latStep() const { - return lats[1] - lats[0]; -} - -int UVGrid::timeStep() const { - return times[1] - times[0]; -} - -void UVGrid::streamSlice(ostream &os, size_t t) { - for (int x = 0; x < latSize; x++) { - for (int y = 0; y < lonSize; y++) { - auto vel = (*this)[t,x,y]; - os << vel << " "; - } - os << endl; - } -} \ No newline at end of file diff --git a/advection/src/main.cpp b/advection/src/main.cpp deleted file mode 100644 index 0e0fc01..0000000 --- a/advection/src/main.cpp +++ /dev/null @@ -1,95 +0,0 @@ -#include -#include -#include - -#include "interpolate.h" -#include "Vel.h" -#include "EulerAdvectionKernel.h" -#include "RK4AdvectionKernel.h" -#include "interpolate.h" - -#define NotAKernelError "Template parameter T must derive from AdvectionKernel" - -using namespace std; - -template -void advectForSomeTime(const UVGrid &uvGrid, const AdvectionKernelImpl &kernel, double latstart, double lonstart, int i, char colour[10]) { - - // Require at compile time that kernel derives from the abstract class AdvectionKernel - static_assert(std::is_base_of::value, NotAKernelError); - - double lat1 = latstart, lon1 = lonstart; - for(int time = 0; time <= 31536000.; time += AdvectionKernel::DT) { -// cout << setprecision(8) << lat1 << "," << setprecision(8) << lon1 << ",end" << i << "," << colour << endl; - try { - auto [templat, templon] = kernel.advect(time, lat1, lon1); - lat1 = templat; - lon1 = templon; - } catch (const out_of_range& e) { - cerr << "broke out of loop!" << endl; - time = 31536001; - } - } - cout << setprecision(8) << latstart << "," << setprecision(8) << lonstart << ",begin" << i << "," << colour << endl; - cout << setprecision(8) << lat1 << "," << setprecision(8) << lon1 << ",end" << i << "," << colour << endl; -} - -void testGridIndexing(const UVGrid *uvGrid) { - int time = 20000; - cout << "=== land === (should all give 0)" << endl; - cout << bilinearinterpolate(*uvGrid, time, 53.80956379699079, -1.6496306344654406) << endl; - cout << bilinearinterpolate(*uvGrid, time, 55.31428895563707, -2.851581041325997) << endl; - cout << bilinearinterpolate(*uvGrid, time, 47.71548983067583, -1.8704054037408626) << endl; - cout << bilinearinterpolate(*uvGrid, time, 56.23521060314398, 8.505979324950573) << endl; - cout << bilinearinterpolate(*uvGrid, time, 53.135645440244375, 8.505979324950573) << endl; - cout << bilinearinterpolate(*uvGrid, time, 56.44761278775708, -4.140629303756164) << endl; - cout << bilinearinterpolate(*uvGrid, time, 52.67625153110339, 0.8978569759455872) << endl; - cout << bilinearinterpolate(*uvGrid, time, 52.07154079279377, 4.627951041411331) << endl; - - cout << "=== ocean === (should give not 0)" << endl; - cout << bilinearinterpolate(*uvGrid, time, 47.43923166616274, -4.985451481829083) << endl; - cout << bilinearinterpolate(*uvGrid, time, 50.68943556852362, -9.306162999561733) << endl; - cout << bilinearinterpolate(*uvGrid, time, 53.70606799886677, -4.518347647034465) << endl; - cout << bilinearinterpolate(*uvGrid, time, 60.57987114267971, -12.208262973672621) << endl; - cout << bilinearinterpolate(*uvGrid, time, 46.532221548197285, -13.408189172582638) << endl; - cout << bilinearinterpolate(*uvGrid, time, 50.92725094937812, 1.3975824052375256) << endl; - cout << bilinearinterpolate(*uvGrid, time, 51.4028921682209, 2.4059571950925203) << endl; - cout << bilinearinterpolate(*uvGrid, time, 53.448445236769004, 0.7996966058017515) << endl; -// cout << bilinearinterpolate(*uvGrid, time, ) << endl; -} - -int main() { - std::shared_ptr uvGrid = std::make_shared(); - - uvGrid->streamSlice(cout, 900); - - auto kernelRK4 = RK4AdvectionKernel(uvGrid); - -// You can use https://maps.co/gis/ to visualise these points - cout << "======= RK4 Integration =======" << endl; - advectForSomeTime(*uvGrid, kernelRK4, 53.53407391652826, 6.274975037862238, 0, "#ADD8E6"); - advectForSomeTime(*uvGrid, kernelRK4, 53.494053820479365, 5.673454142386921, 1, "#DC143C"); - advectForSomeTime(*uvGrid, kernelRK4, 53.49321966653616, 5.681867022043919, 2, "#50C878"); - advectForSomeTime(*uvGrid, kernelRK4, 53.581548701694324, 6.552600066543153, 3, "#FFEA00"); - advectForSomeTime(*uvGrid, kernelRK4, 53.431446729744124, 5.241592961691523, 4, "#663399"); - advectForSomeTime(*uvGrid, kernelRK4, 53.27913608324572, 4.82094897884165, 5, "#FFA500"); - advectForSomeTime(*uvGrid, kernelRK4, 53.18597595482688, 4.767667388308705, 6, "#008080"); - advectForSomeTime(*uvGrid, kernelRK4, 53.01592078792383, 4.6064205160882, 7, "#FFB6C1"); - advectForSomeTime(*uvGrid, kernelRK4, 52.72816940158886, 4.5853883152993635, 8, "#36454F"); // on land - advectForSomeTime(*uvGrid, kernelRK4, 52.56142091881038, 4.502661662924255, 9, "#1E90FF"); // Dodger Blue - advectForSomeTime(*uvGrid, kernelRK4, 52.23202593893584, 4.2825246383181845, 10, "#FFD700"); // Gold - advectForSomeTime(*uvGrid, kernelRK4, 52.08062567609582, 4.112864890830927, 11, "#6A5ACD"); // Slate Blue - advectForSomeTime(*uvGrid, kernelRK4, 51.89497719759734, 3.8114033568921686, 12, "#20B2AA"); // Light Sea Green - advectForSomeTime(*uvGrid, kernelRK4, 51.752848503723634, 3.664177951809339, 13, "#FF69B4"); // Hot Pink - advectForSomeTime(*uvGrid, kernelRK4, 51.64595756528835, 3.626319993352851, 14, "#800080"); // Purple - advectForSomeTime(*uvGrid, kernelRK4, 51.55140730645238, 3.4326152213887986, 15, "#FF4500"); // Orange Red - advectForSomeTime(*uvGrid, kernelRK4, 51.45679776223422, 3.4452813365018384, 16, "#A52A2A"); // Brown - advectForSomeTime(*uvGrid, kernelRK4, 51.41444662720727, 3.4648562416765363, 17, "#4682B4"); // Steel Blue - advectForSomeTime(*uvGrid, kernelRK4, 51.37421261203866, 3.2449264214689455, 18, "#FF6347"); // Tomato - advectForSomeTime(*uvGrid, kernelRK4, 51.29651848898365, 2.9547572241424773, 19, "#008000"); // Green - advectForSomeTime(*uvGrid, kernelRK4, 51.19705098468974, 2.7647654914530024, 20, "#B8860B"); // Dark Goldenrod - advectForSomeTime(*uvGrid, kernelRK4, 51.114719857442665, 2.577076679365129, 21, "#FFC0CB"); // Pink -// advectForSomeTime(*uvGrid, kernelRK4, ,0); - - return 0; -} diff --git a/vtk/src/CMakeLists.txt b/vtk/src/CMakeLists.txt index ef2a3aa..cf0b9d3 100644 --- a/vtk/src/CMakeLists.txt +++ b/vtk/src/CMakeLists.txt @@ -2,6 +2,9 @@ cmake_minimum_required(VERSION 3.12 FATAL_ERROR) project(VtkBase) +set(CMAKE_CXX_STANDARD 23) +set(CMAKE_CXX_STANDARD_REQUIRED ON) + set(CMAKE_EXPORT_COMPILE_COMMANDS ON) find_package(VTK COMPONENTS @@ -47,6 +50,19 @@ add_executable(VtkBase MACOSX_BUNDLE main.cpp commands/SpawnPointCallback.h commands/SpawnPointCallback.cpp CartographicTransformation.cpp + advection/AdvectionKernel.h + advection/EulerAdvectionKernel.cpp + advection/EulerAdvectionKernel.h + advection/interpolate.cpp + advection/interpolate.h + advection/readdata.cpp + advection/readdata.h + advection/RK4AdvectionKernel.cpp + advection/RK4AdvectionKernel.h + advection/UVGrid.cpp + advection/UVGrid.h + advection/Vel.cpp + advection/Vel.h ) execute_process( diff --git a/vtk/src/CartographicTransformation.cpp b/vtk/src/CartographicTransformation.cpp index 54c528e..5518d61 100644 --- a/vtk/src/CartographicTransformation.cpp +++ b/vtk/src/CartographicTransformation.cpp @@ -15,11 +15,11 @@ vtkSmartPointer createNormalisedCamera() { return camera; } -vtkSmartPointer getCartographicTransformMatrix() { - const double XMin = -15.875; - const double XMax = 12.875; - const double YMin = 46.125; - const double YMax = 62.625; +vtkSmartPointer getCartographicTransformMatrix(const std::shared_ptr uvGrid) { + const double XMin = uvGrid->lons.front(); + const double XMax = uvGrid->lons.back(); + const double YMin = uvGrid->lats.front(); + const double YMax = uvGrid->lats.back(); double eyeTransform[] = { 2/(XMax-XMin), 0, 0, -(XMax+XMin)/(XMax-XMin), @@ -34,10 +34,10 @@ vtkSmartPointer getCartographicTransformMatrix() { } // Assumes Normalised camera is used -vtkSmartPointer createCartographicTransformFilter() { +vtkSmartPointer createCartographicTransformFilter(const std::shared_ptr uvGrid) { vtkNew transform; - transform->SetMatrix(getCartographicTransformMatrix()); + transform->SetMatrix(getCartographicTransformMatrix(uvGrid)); vtkSmartPointer transformFilter = vtkSmartPointer::New(); transformFilter->SetTransform(transform); diff --git a/vtk/src/CartographicTransformation.h b/vtk/src/CartographicTransformation.h index 56ffbeb..334cfb4 100644 --- a/vtk/src/CartographicTransformation.h +++ b/vtk/src/CartographicTransformation.h @@ -1,5 +1,7 @@ +#include #include #include +#include "advection/UVGrid.h" #ifndef VTKBASE_NORMALISEDCARTOGRAPHICCAMERA_H #define VTKBASE_NORMALISEDCARTOGRAPHICCAMERA_H @@ -16,15 +18,14 @@ vtkSmartPointer createNormalisedCamera(); /** * Constructs a 4x4 projection matrix that maps homogenious (longitude, latitude, 0, 1) points * to the normalised space. - * TODO: This will soon require UVGrid as a parameter after the advection code is merged properly. * TODO: This transformation has room for improvement see: * https://github.com/MakeNEnjoy/interactive-track-and-trace/issues/12 * @return pointer to 4x4 matrix */ -vtkSmartPointer getCartographicTransformMatrix(); +vtkSmartPointer getCartographicTransformMatrix(const std::shared_ptr uvGrid); /** * Convenience function that converts the 4x4 projection matrix into a vtkTransformFilter * @return pointer to transform filter */ -vtkSmartPointer createCartographicTransformFilter(); \ No newline at end of file +vtkSmartPointer createCartographicTransformFilter(const std::shared_ptr uvGrid); diff --git a/vtk/src/Program.cpp b/vtk/src/Program.cpp index bcb3e61..5f1cf3c 100644 --- a/vtk/src/Program.cpp +++ b/vtk/src/Program.cpp @@ -33,21 +33,22 @@ void Program::setWinProperties() { interact->SetInteractorStyle(style); } -void Program::setupTimer() { +void Program::setupTimer(int dt) { auto callback = vtkSmartPointer::New(this); callback->SetClientData(this); + callback->setDt(dt); this->interact->AddObserver(vtkCommand::TimerEvent, callback); this->interact->AddObserver(vtkCommand::KeyPressEvent, callback); this->interact->CreateRepeatingTimer(17); // 60 fps == 1000 / 60 == 16.7 ms per frame } -Program::Program() { +Program::Program(int timerDT) { this->win = vtkSmartPointer::New(); this->interact = vtkSmartPointer::New(); this->win->SetNumberOfLayers(0); setWinProperties(); - setupTimer(); + setupTimer(timerDT); } diff --git a/vtk/src/Program.h b/vtk/src/Program.h index 4097e4c..fbe5023 100644 --- a/vtk/src/Program.h +++ b/vtk/src/Program.h @@ -31,14 +31,14 @@ private: /** This function sets up and connects a TimerCallbackCommand with the program. */ - void setupTimer(); + void setupTimer(int dt); void setupInteractions(); public: /** Constructor. */ - Program(); + Program(int timerDT); /** This function adds a new layer (and thus vtkRenderer) to the program. * The layer is expected to set its own position in the vtkRenderWindow layer system. diff --git a/advection/src/AdvectionKernel.h b/vtk/src/advection/AdvectionKernel.h similarity index 88% rename from advection/src/AdvectionKernel.h rename to vtk/src/advection/AdvectionKernel.h index 0170c79..3515e6c 100644 --- a/advection/src/AdvectionKernel.h +++ b/vtk/src/advection/AdvectionKernel.h @@ -5,13 +5,11 @@ #include "Vel.h" - /* * Implement this class for every integration method. */ class AdvectionKernel { public: - const static int DT = 60 * 15; // 60 sec/min * 15 mins /** * This function must take a time, latitude and longitude of a particle and must output * a new latitude and longitude after being advected once for AdvectionKernel::DT time as defined above. @@ -20,13 +18,14 @@ public: * @param longitude Longitude of particle * @return A pair of latitude and longitude of particle. */ - virtual std::pair advect(int time, double latitude, double longitude) const = 0; + virtual std::pair advect(int time, double latitude, double longitude, int dt) const = 0; // Taken from Parcels https://github.com/OceanParcels/parcels/blob/daa4b062ed8ae0b2be3d87367d6b45599d6f95db/parcels/tools/converters.py#L155 const static double metreToDegrees(double metre) { return metre / 1000. / 1.852 / 60.; } + virtual ~AdvectionKernel() = default; // Apparently I need this, idk why }; #endif //ADVECTION_ADVECTIONKERNEL_H diff --git a/vtk/src/advection/EulerAdvectionKernel.cpp b/vtk/src/advection/EulerAdvectionKernel.cpp new file mode 100644 index 0000000..f9a4b52 --- /dev/null +++ b/vtk/src/advection/EulerAdvectionKernel.cpp @@ -0,0 +1,13 @@ + +#include "EulerAdvectionKernel.h" +#include "interpolate.h" + +using namespace std; + +EulerAdvectionKernel::EulerAdvectionKernel(std::shared_ptr grid) : grid(grid) {} + +std::pair EulerAdvectionKernel::advect(int time, double latitude, double longitude, int dt) const { + auto [u, v] = bilinearinterpolate(*grid, time, latitude, longitude); + + return {latitude + metreToDegrees(v * dt), longitude + metreToDegrees(u * dt)}; +} diff --git a/advection/src/EulerAdvectionKernel.h b/vtk/src/advection/EulerAdvectionKernel.h similarity index 91% rename from advection/src/EulerAdvectionKernel.h rename to vtk/src/advection/EulerAdvectionKernel.h index d9893cb..df3d3d8 100644 --- a/advection/src/EulerAdvectionKernel.h +++ b/vtk/src/advection/EulerAdvectionKernel.h @@ -3,6 +3,7 @@ #include "AdvectionKernel.h" #include "UVGrid.h" +#include /** * Implementation of AdvectionKernel. @@ -17,7 +18,7 @@ private: std::shared_ptr grid; public: explicit EulerAdvectionKernel(std::shared_ptr grid); - std::pair advect(int time, double latitude, double longitude) const override; + std::pair advect(int time, double latitude, double longitude, int dt) const override; }; diff --git a/advection/src/RK4AdvectionKernel.cpp b/vtk/src/advection/RK4AdvectionKernel.cpp similarity index 61% rename from advection/src/RK4AdvectionKernel.cpp rename to vtk/src/advection/RK4AdvectionKernel.cpp index 1edc9f9..85104b5 100644 --- a/advection/src/RK4AdvectionKernel.cpp +++ b/vtk/src/advection/RK4AdvectionKernel.cpp @@ -5,31 +5,31 @@ using namespace std; RK4AdvectionKernel::RK4AdvectionKernel(std::shared_ptr grid): grid(grid) { } -std::pair RK4AdvectionKernel::advect(int time, double latitude, double longitude) const { +std::pair RK4AdvectionKernel::advect(int time, double latitude, double longitude, int dt) const { auto [u1, v1] = bilinearinterpolate(*grid, time, latitude, longitude); // lon1, lat1 = (particle.lon + u1*.5*particle.dt, particle.lat + v1*.5*particle.dt); - double lon1 = longitude + metreToDegrees(u1 * 0.5*DT); - double lat1 = latitude + metreToDegrees(v1 * 0.5*DT); + double lon1 = longitude + metreToDegrees(u1 * 0.5*dt); + double lat1 = latitude + metreToDegrees(v1 * 0.5*dt); // (u2, v2) = fieldset.UV[time + .5 * particle.dt, particle.depth, lat1, lon1, particle] - auto [u2, v2] = bilinearinterpolate(*grid, time + 0.5 * DT, lat1, lon1); + auto [u2, v2] = bilinearinterpolate(*grid, time + 0.5 * dt, lat1, lon1); // lon2, lat2 = (particle.lon + u2*.5*particle.dt, particle.lat + v2*.5*particle.dt) - double lon2 = longitude + metreToDegrees(u2 * 0.5 * DT); - double lat2 = latitude + metreToDegrees(v2 * 0.5 * DT); + double lon2 = longitude + metreToDegrees(u2 * 0.5 * dt); + double lat2 = latitude + metreToDegrees(v2 * 0.5 * dt); // (u3, v3) = fieldset.UV[time + .5 * particle.dt, particle.depth, lat2, lon2, particle] - auto [u3, v3] = bilinearinterpolate(*grid, time + 0.5 * DT, lat2, lon2); + auto [u3, v3] = bilinearinterpolate(*grid, time + 0.5 * dt, lat2, lon2); // lon3, lat3 = (particle.lon + u3*particle.dt, particle.lat + v3*particle.dt) - double lon3 = longitude + metreToDegrees(u3 * DT); - double lat3 = latitude + metreToDegrees(v3 * DT); + double lon3 = longitude + metreToDegrees(u3 * dt); + double lat3 = latitude + metreToDegrees(v3 * dt); // (u4, v4) = fieldset.UV[time + particle.dt, particle.depth, lat3, lon3, particle] - auto [u4, v4] = bilinearinterpolate(*grid, time + DT, lat3, lon3); + auto [u4, v4] = bilinearinterpolate(*grid, time + dt, lat3, lon3); - double lonFinal = longitude + metreToDegrees((u1 + 2 * u2 + 2 * u3 + u4) / 6.0 * DT); - double latFinal = latitude + metreToDegrees((v1 + 2 * v2 + 2 * v3 + v4) / 6.0 * DT); + double lonFinal = longitude + metreToDegrees((u1 + 2 * u2 + 2 * u3 + u4) / 6.0 * dt); + double latFinal = latitude + metreToDegrees((v1 + 2 * v2 + 2 * v3 + v4) / 6.0 * dt); return {latFinal, lonFinal}; } diff --git a/advection/src/RK4AdvectionKernel.h b/vtk/src/advection/RK4AdvectionKernel.h similarity index 91% rename from advection/src/RK4AdvectionKernel.h rename to vtk/src/advection/RK4AdvectionKernel.h index 6b6c88d..af44cdd 100644 --- a/advection/src/RK4AdvectionKernel.h +++ b/vtk/src/advection/RK4AdvectionKernel.h @@ -3,6 +3,7 @@ #include "AdvectionKernel.h" #include "UVGrid.h" +#include /** * Implementation of Advection kernel using RK4 integration @@ -14,7 +15,7 @@ private: std::shared_ptr grid; public: explicit RK4AdvectionKernel(std::shared_ptr grid); - std::pair advect(int time, double latitude, double longitude) const override; + std::pair advect(int time, double latitude, double longitude, int dt) const override; }; diff --git a/vtk/src/advection/UVGrid.cpp b/vtk/src/advection/UVGrid.cpp new file mode 100644 index 0000000..f35aa72 --- /dev/null +++ b/vtk/src/advection/UVGrid.cpp @@ -0,0 +1,66 @@ +#include + +#include "UVGrid.h" +#include "readdata.h" + +#define sizeError2 "The sizes of the hydrodynamic data files are different" +#define sizeError "The sizes of the hydrodynamicU or -V files does not correspond with the sizes of the grid file" + +using namespace std; + +UVGrid::UVGrid() { + auto us = readHydrodynamicU(); + auto vs = readHydrodynamicV(); + if (us.size() != vs.size()) { + throw domain_error(sizeError2); + } + + tie(times, lats, lons) = readGrid(); + + timeSize = times.size(); + latSize = lats.size(); + lonSize = lons.size(); + + size_t gridSize = timeSize * latSize * lonSize; + if (gridSize != us.size()) { + throw domain_error(sizeError); + } + + uvData.reserve(gridSize); + + for (auto vel: views::zip(us, vs)) { + uvData.push_back(Vel(vel)); + } +} + +const Vel &UVGrid::operator[](size_t timeIndex, size_t latIndex, size_t lonIndex) const { + if (timeIndex < 0 or timeIndex >= timeSize + or latIndex < 0 or latIndex >= latSize + or lonIndex < 0 or lonIndex >= lonSize) { + throw std::out_of_range("Index out of bounds"); + } + size_t index = timeIndex * (latSize * lonSize) + latIndex * lonSize + lonIndex; + return uvData[index]; +} + +double UVGrid::lonStep() const { + return lons[1] - lons[0]; +} + +double UVGrid::latStep() const { + return lats[1] - lats[0]; +} + +int UVGrid::timeStep() const { + return times[1] - times[0]; +} + +void UVGrid::streamSlice(ostream &os, size_t t) { + for (int x = 0; x < latSize; x++) { + for (int y = 0; y < lonSize; y++) { + auto vel = (*this)[t, x, y]; + os << vel << " "; + } + os << endl; + } +} \ No newline at end of file diff --git a/advection/src/UVGrid.h b/vtk/src/advection/UVGrid.h similarity index 100% rename from advection/src/UVGrid.h rename to vtk/src/advection/UVGrid.h diff --git a/advection/src/Vel.cpp b/vtk/src/advection/Vel.cpp similarity index 100% rename from advection/src/Vel.cpp rename to vtk/src/advection/Vel.cpp diff --git a/advection/src/Vel.h b/vtk/src/advection/Vel.h similarity index 100% rename from advection/src/Vel.h rename to vtk/src/advection/Vel.h diff --git a/advection/src/interpolate.cpp b/vtk/src/advection/interpolate.cpp similarity index 100% rename from advection/src/interpolate.cpp rename to vtk/src/advection/interpolate.cpp diff --git a/advection/src/interpolate.h b/vtk/src/advection/interpolate.h similarity index 100% rename from advection/src/interpolate.h rename to vtk/src/advection/interpolate.h diff --git a/advection/src/readdata.cpp b/vtk/src/advection/readdata.cpp similarity index 100% rename from advection/src/readdata.cpp rename to vtk/src/advection/readdata.cpp diff --git a/advection/src/readdata.h b/vtk/src/advection/readdata.h similarity index 100% rename from advection/src/readdata.h rename to vtk/src/advection/readdata.h diff --git a/vtk/src/commands/SpawnPointCallback.cpp b/vtk/src/commands/SpawnPointCallback.cpp index a0f38cd..21a0bcb 100644 --- a/vtk/src/commands/SpawnPointCallback.cpp +++ b/vtk/src/commands/SpawnPointCallback.cpp @@ -9,68 +9,74 @@ #include "../CartographicTransformation.h" -void convertDisplayToWorld(vtkRenderer* renderer, int x, int y, double *worldPos) { - double displayPos[3] = {static_cast(x), static_cast(y), 0.0}; - renderer->SetDisplayPoint(displayPos); - renderer->DisplayToWorld(); - renderer->GetWorldPoint(worldPos); +void convertDisplayToWorld(vtkRenderer *renderer, int x, int y, double *worldPos) { + double displayPos[3] = {static_cast(x), static_cast(y), 0.0}; + renderer->SetDisplayPoint(displayPos); + renderer->DisplayToWorld(); + renderer->GetWorldPoint(worldPos); } void SpawnPointCallback::Execute(vtkObject *caller, unsigned long evId, void *callData) { - // Note the use of reinterpret_cast to cast the caller to the expected type. - auto interactor = reinterpret_cast(caller); + // Note the use of reinterpret_cast to cast the caller to the expected type. + auto interactor = reinterpret_cast(caller); - if (evId == vtkCommand::LeftButtonPressEvent) { - dragging = true; - } - if (evId == vtkCommand::LeftButtonReleaseEvent) { - dragging = false; - } - if (!dragging) { - return; - } + if (evId == vtkCommand::LeftButtonPressEvent) { + dragging = true; + } + if (evId == vtkCommand::LeftButtonReleaseEvent) { + dragging = false; + } + if (!dragging) { + return; + } - int x, y; - interactor->GetEventPosition(x, y); + int x, y; + interactor->GetEventPosition(x, y); - double worldPos[4] = {2, 0 ,0, 0}; - double displayPos[3] = {static_cast(x), static_cast(y), 0.0}; - ren->SetDisplayPoint(displayPos); - ren->DisplayToWorld(); - ren->GetWorldPoint(worldPos); - inverseCartographicProjection->MultiplyPoint(worldPos, worldPos); - cout << "clicked on lon = " << worldPos[0] << " and lat = " << worldPos[1] << endl; + double worldPos[4] = {2, 0, 0, 0}; + double displayPos[3] = {static_cast(x), static_cast(y), 0.0}; + ren->SetDisplayPoint(displayPos); + ren->DisplayToWorld(); + ren->GetWorldPoint(worldPos); + inverseCartographicProjection->MultiplyPoint(worldPos, worldPos); + cout << "clicked on lon = " << worldPos[0] << " and lat = " << worldPos[1] << endl; - vtkIdType id = points->InsertNextPoint(worldPos[0], worldPos[1], 0); - data->SetPoints(points); + vtkIdType id = points->InsertNextPoint(worldPos[0], worldPos[1], 0); + data->SetPoints(points); - vtkSmartPointer vertex = vtkSmartPointer::New(); - vertex->GetPointIds()->SetId(0, id); + vtkSmartPointer vertex = vtkSmartPointer::New(); + vertex->GetPointIds()->SetId(0, id); - vtkSmartPointer vertices = vtkSmartPointer::New(); - vertices->InsertNextCell(vertex); - data->SetVerts(vertices); - ren->GetRenderWindow()->Render(); + vtkSmartPointer vertices = vtkSmartPointer::New(); + vertices->InsertNextCell(vertex); + data->SetVerts(vertices); + ren->GetRenderWindow()->Render(); } -SpawnPointCallback::SpawnPointCallback() : data(nullptr), points(nullptr), inverseCartographicProjection(nullptr) { - inverseCartographicProjection = getCartographicTransformMatrix(); - inverseCartographicProjection->Invert(); -} +SpawnPointCallback::SpawnPointCallback() : data(nullptr), + points(nullptr), + inverseCartographicProjection(nullptr), + uvGrid(nullptr) { } SpawnPointCallback *SpawnPointCallback::New() { return new SpawnPointCallback; } void SpawnPointCallback::setData(const vtkSmartPointer &data) { - this->data = data; + this->data = data; } void SpawnPointCallback::setPoints(const vtkSmartPointer &points) { - this->points = points; + this->points = points; } void SpawnPointCallback::setRen(const vtkSmartPointer &ren) { - this->ren = ren; + this->ren = ren; +} + +void SpawnPointCallback::setUVGrid(const std::shared_ptr &uvGrid) { + this->uvGrid = uvGrid; + inverseCartographicProjection = getCartographicTransformMatrix(uvGrid); + inverseCartographicProjection->Invert(); } diff --git a/vtk/src/commands/SpawnPointCallback.h b/vtk/src/commands/SpawnPointCallback.h index bef6ca4..a7dec2e 100644 --- a/vtk/src/commands/SpawnPointCallback.h +++ b/vtk/src/commands/SpawnPointCallback.h @@ -2,31 +2,39 @@ #define VTKBASE_SPAWNPOINTCALLBACK_H +#include #include #include #include #include #include +#include "../advection/UVGrid.h" class SpawnPointCallback : public vtkCallbackCommand { public: - static SpawnPointCallback *New(); - SpawnPointCallback(); + static SpawnPointCallback *New(); - void setPoints(const vtkSmartPointer &points); + SpawnPointCallback(); - void setData(const vtkSmartPointer &data); + void setPoints(const vtkSmartPointer &points); + + void setData(const vtkSmartPointer &data); + + void setRen(const vtkSmartPointer &ren); + + void setUVGrid(const std::shared_ptr &uvGrid); - void setRen(const vtkSmartPointer &ren); private: - vtkSmartPointer data; - vtkSmartPointer points; - vtkSmartPointer ren; - vtkSmartPointer inverseCartographicProjection; + vtkSmartPointer data; + vtkSmartPointer points; + vtkSmartPointer ren; + std::shared_ptr uvGrid; + vtkSmartPointer inverseCartographicProjection; - void Execute(vtkObject *caller, unsigned long evId, void *callData) override; - bool dragging = false; + void Execute(vtkObject *caller, unsigned long evId, void *callData) override; + + bool dragging = false; }; diff --git a/vtk/src/commands/TimerCallbackCommand.cpp b/vtk/src/commands/TimerCallbackCommand.cpp index 8e91b58..6fad2ca 100644 --- a/vtk/src/commands/TimerCallbackCommand.cpp +++ b/vtk/src/commands/TimerCallbackCommand.cpp @@ -38,3 +38,7 @@ void TimerCallbackCommand::setProgram(Program *program) { void TimerCallbackCommand::setPaused(const bool val) { this->paused = val; } + +void TimerCallbackCommand::setDt(int dt) { + this->dt = dt; +} diff --git a/vtk/src/commands/TimerCallbackCommand.h b/vtk/src/commands/TimerCallbackCommand.h index 2ff65ea..7123f4d 100644 --- a/vtk/src/commands/TimerCallbackCommand.h +++ b/vtk/src/commands/TimerCallbackCommand.h @@ -13,6 +13,8 @@ public: void setProgram(Program *program); void setPaused(const bool val); + void setDt(int dt); + private: int time; int dt; diff --git a/vtk/src/layers/EGlyphLayer.cpp b/vtk/src/layers/EGlyphLayer.cpp index cb91825..93db577 100644 --- a/vtk/src/layers/EGlyphLayer.cpp +++ b/vtk/src/layers/EGlyphLayer.cpp @@ -11,43 +11,20 @@ #include #include #include -#include #include #include "../CartographicTransformation.h" +#include "../advection/readdata.h" +#include "../advection/interpolate.h" -using namespace netCDF; using namespace std; -template -vector getVarVector(const NcVar &var) { - int length = 1; - for (NcDim dim : var.getDims()) { - length *= dim.getSize(); - } - - vector vec(length); - - var.getVar(vec.data()); - - return vec; -} - -tuple, vector, vector> readGrid() { - netCDF::NcFile data("../../../../data/grid.h5", netCDF::NcFile::read); - multimap< string, NcVar > vars = data.getVars(); - vector time = getVarVector(vars.find("times")->second); - vector longitude = getVarVector(vars.find("longitude")->second); - vector latitude = getVarVector(vars.find("latitude")->second); - - return {time, latitude, longitude}; -} - - -EGlyphLayer::EGlyphLayer() { +EGlyphLayer::EGlyphLayer(std::shared_ptr uvGrid) { this->ren = vtkSmartPointer::New(); this->ren->SetLayer(1); this->ren->InteractiveOff(); + this->uvGrid = uvGrid; + this->data = vtkSmartPointer::New(); this->direction = vtkSmartPointer::New(); this->direction->SetName("direction"); @@ -57,31 +34,33 @@ EGlyphLayer::EGlyphLayer() { void EGlyphLayer::readCoordinates() { vtkNew points; - auto [times, lats, lons] = readGrid(); // FIXME: import Robin's readData function and use it - this->numLats = lats.size(); - this->numLons = lons.size(); + this->numLats = uvGrid->lats.size(); + this->numLons = uvGrid->lons.size(); this->direction->SetNumberOfComponents(3); - this->direction->SetNumberOfTuples(numLats*numLons); - points->Allocate(numLats*numLons); + this->direction->SetNumberOfTuples(numLats * numLons); + points->Allocate(numLats * numLons); auto camera = createNormalisedCamera(); ren->SetActiveCamera(camera); int i = 0; - for (double lat : lats) { - for (double lon : lons) { - cout << "lon: " << lon << " lat: " << lat << endl; - direction->SetTuple3(i, 0.45, 0.90, 0); //FIXME: read this info from file + int latIndex = 0; + for (double lat: uvGrid->lats) { + int lonIndex = 0; + for (double lon: uvGrid->lons) { + auto [u, v] = (*uvGrid)[0, latIndex, lonIndex]; + direction->SetTuple3(i, 5*u, 5*v, 0); points->InsertPoint(i++, lon, lat, 0); - // see also https://vtk.org/doc/nightly/html/classvtkPolyDataMapper2D.html + lonIndex++; } + latIndex++; } this->data->SetPoints(points); this->data->GetPointData()->AddArray(this->direction); this->data->GetPointData()->SetActiveVectors("direction"); - vtkSmartPointer transformFilter = createCartographicTransformFilter(); + vtkSmartPointer transformFilter = createCartographicTransformFilter(uvGrid); transformFilter->SetInputData(data); vtkNew arrowSource; @@ -94,8 +73,8 @@ void EGlyphLayer::readCoordinates() { glyph2D->SetInputConnection(transformFilter->GetOutputPort()); glyph2D->OrientOn(); glyph2D->ClampingOn(); - glyph2D->SetScaleModeToScaleByVector(); - glyph2D->SetVectorModeToUseVector(); + glyph2D->SetScaleModeToScaleByVector(); + glyph2D->SetVectorModeToUseVector(); glyph2D->Update(); // vtkNew coordinate; @@ -109,16 +88,22 @@ void EGlyphLayer::readCoordinates() { vtkNew actor; actor->SetMapper(mapper); - actor->GetProperty()->SetColor(0,0,0); + actor->GetProperty()->SetColor(0, 0, 0); actor->GetProperty()->SetOpacity(0.2); - this->ren->AddActor(actor) ; + this->ren->AddActor(actor); } void EGlyphLayer::updateData(int t) { - for (int i=0; i < numLats*numLons; i++) { - this->direction->SetTuple3(i, std::cos(t), std::sin(t), 0); // FIXME: fetch data from file. + int i = 0; + for (int lat = 0; lat < uvGrid->latSize; lat++) { + for (int lon = 0; lon < uvGrid->lonSize; lon++) { + auto [u, v] = (*uvGrid)[t/3600, lat, lon]; + // TODO: 5*v scaling stuff should really be a filter transform + this->direction->SetTuple3(i, 5*u, 5*v, 0); + i++; + } } this->direction->Modified(); } diff --git a/vtk/src/layers/EGlyphLayer.h b/vtk/src/layers/EGlyphLayer.h index 8631f32..7704af5 100644 --- a/vtk/src/layers/EGlyphLayer.h +++ b/vtk/src/layers/EGlyphLayer.h @@ -2,8 +2,11 @@ #define EGLYPHLAYER_H #include "Layer.h" +#include #include +#include "../advection/UVGrid.h" + /** Implements the Layer class for the case of a Eulerian visualization. * Specifically, this class models the eulerian flow-fields of the simulation using the 'glyph' mark and 'direction' and 'form' channels to denote direction and strength of velocities. */ @@ -11,6 +14,7 @@ class EGlyphLayer : public Layer { private: vtkSmartPointer data; vtkSmartPointer direction; + std::shared_ptr uvGrid; int numLats; int numLons; @@ -22,7 +26,7 @@ private: public: /** Constructor. */ - EGlyphLayer(); + EGlyphLayer(std::shared_ptr uvGrid); /** updates the glyphs to reflect the given timestamp in the dataset. * @param t : the time at which to fetch the data. diff --git a/vtk/src/layers/LGlyphLayer.cpp b/vtk/src/layers/LGlyphLayer.cpp index 845d628..35acafd 100644 --- a/vtk/src/layers/LGlyphLayer.cpp +++ b/vtk/src/layers/LGlyphLayer.cpp @@ -17,13 +17,13 @@ #include "../CartographicTransformation.h" - vtkSmartPointer LGlyphLayer::createSpawnPointCallback() { - auto newPointCallBack = vtkSmartPointer::New(); - newPointCallBack->setData(data); - newPointCallBack->setPoints(points); - newPointCallBack->setRen(ren); - return newPointCallBack; + auto newPointCallBack = vtkSmartPointer::New(); + newPointCallBack->setData(data); + newPointCallBack->setPoints(points); + newPointCallBack->setRen(ren); + newPointCallBack->setUVGrid(uvGrid); + return newPointCallBack; } // Further notes; current thinking is to allow tracking a particle's age by using a scalar array in the VtkPolyData. This would be incremented for every tick/updateData function call. @@ -31,75 +31,72 @@ vtkSmartPointer LGlyphLayer::createSpawnPointCallback() { // // TODO: modelling all this in vtkClasses is workable, but ideally i would want to work with a native C++ class. See if this is doable and feasible. -LGlyphLayer::LGlyphLayer() { - this->ren = vtkSmartPointer::New(); - this->ren->SetLayer(2); +LGlyphLayer::LGlyphLayer(std::shared_ptr uvGrid, std::unique_ptr advectionKernel) { + this->ren = vtkSmartPointer::New(); + this->ren->SetLayer(2); - this->points = vtkSmartPointer::New(); - this->data = vtkSmartPointer::New(); - this->data->SetPoints(this->points); + this->points = vtkSmartPointer::New(); + this->data = vtkSmartPointer::New(); + this->data->SetPoints(this->points); - auto camera = createNormalisedCamera(); - ren->SetActiveCamera(camera); + advector = std::move(advectionKernel); + this->uvGrid = uvGrid; - auto transform = createCartographicTransformFilter(); + auto camera = createNormalisedCamera(); + ren->SetActiveCamera(camera); - vtkSmartPointer transformFilter = createCartographicTransformFilter(); - transformFilter->SetInputData(data); + vtkSmartPointer transformFilter = createCartographicTransformFilter(uvGrid); + transformFilter->SetInputData(data); - vtkNew circleSource; - circleSource->SetGlyphTypeToCircle(); - circleSource->SetScale(0.05); - circleSource->Update(); + vtkNew circleSource; + circleSource->SetGlyphTypeToCircle(); + circleSource->SetScale(0.05); + circleSource->Update(); - vtkNew glyph2D; - glyph2D->SetSourceConnection(circleSource->GetOutputPort()); - glyph2D->SetInputConnection(transformFilter->GetOutputPort()); - glyph2D->SetColorModeToColorByScalar(); - glyph2D->Update(); + vtkNew glyph2D; + glyph2D->SetSourceConnection(circleSource->GetOutputPort()); + glyph2D->SetInputConnection(transformFilter->GetOutputPort()); + glyph2D->SetColorModeToColorByScalar(); + glyph2D->Update(); - vtkNew mapper; - mapper->SetInputConnection(glyph2D->GetOutputPort()); - mapper->Update(); + vtkNew mapper; + mapper->SetInputConnection(glyph2D->GetOutputPort()); + mapper->Update(); - vtkNew actor; - actor->SetMapper(mapper); + vtkNew actor; + actor->SetMapper(mapper); - this->ren->AddActor(actor); + this->ren->AddActor(actor); } // creates a few points so we can test the updateData function void LGlyphLayer::spoofPoints() { - this->points->InsertNextPoint(-4.125, 61.375 , 0); - this->points->InsertNextPoint(6.532949683882039, 53.24308582564463, 0); // Coordinates of Zernike - this->points->InsertNextPoint(5.315307819255385, 60.40001057122271, 0); // Coordinates of Bergen - this->points->InsertNextPoint( 6.646210231365825, 46.52346296009023, 0); // Coordinates of Lausanne - this->points->InsertNextPoint(-6.553894313570932, 62.39522131195857, 0); // Coordinates of the top of the Faroe islands + this->points->InsertNextPoint(-4.125, 61.375, 0); + this->points->InsertNextPoint(6.532949683882039, 53.24308582564463, 0); // Coordinates of Zernike + this->points->InsertNextPoint(5.315307819255385, 60.40001057122271, 0); // Coordinates of Bergen + this->points->InsertNextPoint(6.646210231365825, 46.52346296009023, 0); // Coordinates of Lausanne + this->points->InsertNextPoint(-6.553894313570932, 62.39522131195857,0); // Coordinates of the top of the Faroe islands - this->points->Modified(); -} - -// returns new coords for a point; used to test the updateData function -std::pair advect(int time, double lat, double lon) { - return {lat + 0.01, lon + 0.01}; + this->points->Modified(); } void LGlyphLayer::updateData(int t) { - double point[3]; - for (vtkIdType n = 0; n < this->points->GetNumberOfPoints(); n++) { - this->points->GetPoint(n, point); - auto [xNew, yNew] = advect(n, point[0], point[1]); - this->points->SetPoint(n, xNew, yNew, 0); + const int SUPERSAMPLINGRATE = 4; + double point[3]; + for (vtkIdType n = 0; n < this->points->GetNumberOfPoints(); n++) { + this->points->GetPoint(n, point); + for (int i = 0; i < SUPERSAMPLINGRATE; i++) { + std::tie(point[1], point[0]) = advector->advect(t, point[1], point[0], (t-lastT)/SUPERSAMPLINGRATE); } - this->points->Modified(); + this->points->SetPoint(n, point[0], point[1], 0); + } + lastT = t; + this->points->Modified(); } void LGlyphLayer::addObservers(vtkSmartPointer interactor) { - auto newPointCallBack = vtkSmartPointer::New(); - newPointCallBack->setData(data); - newPointCallBack->setPoints(points); - newPointCallBack->setRen(ren); - interactor->AddObserver(vtkCommand::LeftButtonPressEvent, newPointCallBack); - interactor->AddObserver(vtkCommand::LeftButtonReleaseEvent, newPointCallBack); - interactor->AddObserver(vtkCommand::MouseMoveEvent, newPointCallBack); + auto newPointCallBack = createSpawnPointCallback(); + interactor->AddObserver(vtkCommand::LeftButtonPressEvent, newPointCallBack); + interactor->AddObserver(vtkCommand::LeftButtonReleaseEvent, newPointCallBack); + interactor->AddObserver(vtkCommand::MouseMoveEvent, newPointCallBack); } diff --git a/vtk/src/layers/LGlyphLayer.h b/vtk/src/layers/LGlyphLayer.h index 22993b4..422cfcb 100644 --- a/vtk/src/layers/LGlyphLayer.h +++ b/vtk/src/layers/LGlyphLayer.h @@ -2,6 +2,7 @@ #define LGLYPHLAYER_H #include "Layer.h" +#include "../advection/AdvectionKernel.h" #include "../commands/SpawnPointCallback.h" #include #include @@ -13,13 +14,14 @@ class LGlyphLayer : public Layer { private: vtkSmartPointer points; vtkSmartPointer data; - - + std::unique_ptr advector; + std::shared_ptr uvGrid; + int lastT = 1000; public: /** Constructor. */ - LGlyphLayer(); + LGlyphLayer(std::shared_ptr uvGrid, std::unique_ptr advectionKernel); /** This function spoofs a few points in the dataset. Mostly used for testing. */ diff --git a/vtk/src/main.cpp b/vtk/src/main.cpp index fb89ea4..49e44bb 100644 --- a/vtk/src/main.cpp +++ b/vtk/src/main.cpp @@ -7,22 +7,30 @@ #include #include #include +#include #include "layers/BackgroundImage.h" #include "layers/EGlyphLayer.h" #include "layers/LGlyphLayer.h" #include "Program.h" +#include "advection/UVGrid.h" +#include "advection/RK4AdvectionKernel.h" using namespace std; +#define DT 60 * 60 // 60 sec/min * 60 mins int main() { - auto l = new LGlyphLayer(); - l->spoofPoints(); + cout << "reading data..." << endl; + shared_ptr uvGrid = std::make_shared(); + auto kernelRK4 = make_unique(uvGrid); + cout << "Starting vtk..." << endl; - Program *program = new Program(); + auto l = new LGlyphLayer(uvGrid, std::move(kernelRK4)); + + Program *program = new Program(DT); program->addLayer(new BackgroundImage("../../../../data/map_661-661.png")); - program->addLayer(new EGlyphLayer()); + program->addLayer(new EGlyphLayer(uvGrid)); program->addLayer(l); program->render();