Merge pull request #22 from MakeNEnjoy/robin-combine-advection-vtk

Robin combine advection vtk

advection/.gitignore

@@ -1,7 +0,0 @@
-.DS_Store
-src/.DS_Store
-src/.cache
-src/build
-.idea
-src/cmake-build-debug
-src/cmake-build-release

advection/README.md

@@ -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 <typename AdvectionKernelImpl>
-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<AdvectionKernel, AdvectionKernelImpl>::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
-```

advection/src/CMakeLists.txt

@@ -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})

advection/src/EulerAdvectionKernel.cpp

@@ -1,13 +0,0 @@
-
-#include "EulerAdvectionKernel.h"
-#include "interpolate.h"
-
-using namespace std;
-
-EulerAdvectionKernel::EulerAdvectionKernel(std::shared_ptr<UVGrid> grid): grid(grid) { }
-
-std::pair<double, double> 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)};
-}

advection/src/UVGrid.cpp

@@ -1,66 +0,0 @@
-#include <ranges>
-
-#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;
-    }
-}

advection/src/main.cpp

@@ -1,95 +0,0 @@
-#include <ranges>
-#include <iomanip>
-#include <stdexcept>
-
-#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 <typename AdvectionKernelImpl>
-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<AdvectionKernel, AdvectionKernelImpl>::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> uvGrid = std::make_shared<UVGrid>();
-
-    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;
-}

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(

vtk/src/CartographicTransformation.cpp

@@ -15,11 +15,11 @@ vtkSmartPointer<vtkCamera> createNormalisedCamera() {
     return camera;
 }
 
-vtkSmartPointer<vtkMatrix4x4> getCartographicTransformMatrix() {
+vtkSmartPointer<vtkMatrix4x4> getCartographicTransformMatrix(const std::shared_ptr<UVGrid> uvGrid) {
-    const double XMin = -15.875;
+    const double XMin = uvGrid->lons.front();
-    const double XMax = 12.875;
+    const double XMax = uvGrid->lons.back();
-    const double YMin = 46.125;
+    const double YMin = uvGrid->lats.front();
-    const double YMax = 62.625;
+    const double YMax = uvGrid->lats.back();
 
     double eyeTransform[] = {
             2/(XMax-XMin), 0, 0, -(XMax+XMin)/(XMax-XMin),
@@ -34,10 +34,10 @@ vtkSmartPointer<vtkMatrix4x4> getCartographicTransformMatrix() {
 }
 
 // Assumes Normalised camera is used
-vtkSmartPointer<vtkTransformFilter> createCartographicTransformFilter() {
+vtkSmartPointer<vtkTransformFilter> createCartographicTransformFilter(const std::shared_ptr<UVGrid> uvGrid) {
     vtkNew<vtkTransform> transform;
 
-    transform->SetMatrix(getCartographicTransformMatrix());
+    transform->SetMatrix(getCartographicTransformMatrix(uvGrid));
 
     vtkSmartPointer<vtkTransformFilter> transformFilter = vtkSmartPointer<vtkTransformFilter>::New();
     transformFilter->SetTransform(transform);

vtk/src/CartographicTransformation.h

@@ -1,5 +1,7 @@
+#include <memory>
 #include <vtkCamera.h>
 #include <vtkTransformFilter.h>
+#include "advection/UVGrid.h"
 
 #ifndef VTKBASE_NORMALISEDCARTOGRAPHICCAMERA_H
 #define VTKBASE_NORMALISEDCARTOGRAPHICCAMERA_H
@@ -16,15 +18,14 @@ vtkSmartPointer<vtkCamera> 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<vtkMatrix4x4> getCartographicTransformMatrix();
+vtkSmartPointer<vtkMatrix4x4> getCartographicTransformMatrix(const std::shared_ptr<UVGrid> uvGrid);
 
 /**
  * Convenience function that converts the 4x4 projection matrix into a vtkTransformFilter
  * @return pointer to transform filter
  */
-vtkSmartPointer<vtkTransformFilter> createCartographicTransformFilter();
+vtkSmartPointer<vtkTransformFilter> createCartographicTransformFilter(const std::shared_ptr<UVGrid> uvGrid);

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<TimerCallbackCommand>::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<vtkRenderWindow>::New();
   this->interact = vtkSmartPointer<vtkRenderWindowInteractor>::New();
 
   this->win->SetNumberOfLayers(0);
   setWinProperties();
-  setupTimer();
+  setupTimer(timerDT);
 }
 
 

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.

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<double, double> advect(int time, double latitude, double longitude) const = 0;
+    virtual std::pair<double, double> 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

vtk/src/advection/EulerAdvectionKernel.cpp

@@ -0,0 +1,13 @@
+
+#include "EulerAdvectionKernel.h"
+#include "interpolate.h"
+
+using namespace std;
+
+EulerAdvectionKernel::EulerAdvectionKernel(std::shared_ptr<UVGrid> grid) : grid(grid) {}
+
+std::pair<double, double> 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)};
+}

vtk/src/advection/EulerAdvectionKernel.h

@@ -3,6 +3,7 @@
 
 #include "AdvectionKernel.h"
 #include "UVGrid.h"
+#include <memory>
 
 /**
  * Implementation of AdvectionKernel.
@@ -17,7 +18,7 @@ private:
     std::shared_ptr<UVGrid> grid;
 public:
     explicit EulerAdvectionKernel(std::shared_ptr<UVGrid> grid);
-    std::pair<double, double> advect(int time, double latitude, double longitude) const override;
+    std::pair<double, double> advect(int time, double latitude, double longitude, int dt) const override;
 
 };
 

vtk/src/advection/RK4AdvectionKernel.cpp

@@ -5,31 +5,31 @@ using namespace std;
 
 RK4AdvectionKernel::RK4AdvectionKernel(std::shared_ptr<UVGrid> grid): grid(grid) { }
 
-std::pair<double, double> RK4AdvectionKernel::advect(int time, double latitude, double longitude) const {
+std::pair<double, double> 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 lon1 = longitude + metreToDegrees(u1 * 0.5*dt);
-    double lat1 = latitude + metreToDegrees(v1 * 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 lon2 = longitude + metreToDegrees(u2 * 0.5 * dt);
-    double lat2 = latitude + metreToDegrees(v2 * 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 lon3 = longitude + metreToDegrees(u3 * dt);
-    double lat3 = latitude + metreToDegrees(v3 * 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 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 latFinal = latitude + metreToDegrees((v1 + 2 * v2 + 2 * v3 + v4) / 6.0 * dt);
 
     return {latFinal, lonFinal};
 }

vtk/src/advection/RK4AdvectionKernel.h

@@ -3,6 +3,7 @@
 
 #include "AdvectionKernel.h"
 #include "UVGrid.h"
+#include <memory>
 
 /**
  * Implementation of Advection kernel using RK4 integration
@@ -14,7 +15,7 @@ private:
     std::shared_ptr<UVGrid> grid;
 public:
     explicit RK4AdvectionKernel(std::shared_ptr<UVGrid> grid);
-    std::pair<double, double> advect(int time, double latitude, double longitude) const override;
+    std::pair<double, double> advect(int time, double latitude, double longitude, int dt) const override;
 
 };
 

vtk/src/advection/UVGrid.cpp

@@ -0,0 +1,66 @@
+#include <ranges>
+
+#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;
+  }
+}

vtk/src/commands/SpawnPointCallback.cpp

@@ -9,7 +9,7 @@
 
 #include "../CartographicTransformation.h"
 
-void convertDisplayToWorld(vtkRenderer* renderer, int x, int y, double *worldPos) {
+void convertDisplayToWorld(vtkRenderer *renderer, int x, int y, double *worldPos) {
   double displayPos[3] = {static_cast<double>(x), static_cast<double>(y), 0.0};
   renderer->SetDisplayPoint(displayPos);
   renderer->DisplayToWorld();
@@ -33,7 +33,7 @@ void SpawnPointCallback::Execute(vtkObject *caller, unsigned long evId, void *ca
   int x, y;
   interactor->GetEventPosition(x, y);
 
-    double worldPos[4] = {2, 0 ,0, 0};
+  double worldPos[4] = {2, 0, 0, 0};
   double displayPos[3] = {static_cast<double>(x), static_cast<double>(y), 0.0};
   ren->SetDisplayPoint(displayPos);
   ren->DisplayToWorld();
@@ -54,10 +54,10 @@ void SpawnPointCallback::Execute(vtkObject *caller, unsigned long evId, void *ca
 }
 
 
-SpawnPointCallback::SpawnPointCallback() : data(nullptr), points(nullptr), inverseCartographicProjection(nullptr) {
+SpawnPointCallback::SpawnPointCallback() : data(nullptr),
-    inverseCartographicProjection = getCartographicTransformMatrix();
+                                           points(nullptr),
-    inverseCartographicProjection->Invert();
+                                           inverseCartographicProjection(nullptr),
-}
+                                           uvGrid(nullptr) { }
 
 SpawnPointCallback *SpawnPointCallback::New() {
     return new SpawnPointCallback;
@@ -74,3 +74,9 @@ void SpawnPointCallback::setPoints(const vtkSmartPointer<vtkPoints> &points) {
 void SpawnPointCallback::setRen(const vtkSmartPointer<vtkRenderer> &ren) {
   this->ren = ren;
 }
+
+void SpawnPointCallback::setUVGrid(const std::shared_ptr<UVGrid> &uvGrid) {
+  this->uvGrid = uvGrid;
+  inverseCartographicProjection = getCartographicTransformMatrix(uvGrid);
+  inverseCartographicProjection->Invert();
+}

vtk/src/commands/SpawnPointCallback.h

@@ -2,16 +2,19 @@
 #define VTKBASE_SPAWNPOINTCALLBACK_H
 
 
+#include <memory>
 #include <vtkCallbackCommand.h>
 #include <vtkRenderWindowInteractor.h>
 #include <vtkPoints.h>
 #include <vtkPolyData.h>
 #include <vtkMatrix4x4.h>
+#include "../advection/UVGrid.h"
 
 class SpawnPointCallback : public vtkCallbackCommand {
 
 public:
   static SpawnPointCallback *New();
+
   SpawnPointCallback();
 
   void setPoints(const vtkSmartPointer<vtkPoints> &points);
@@ -19,13 +22,18 @@ public:
   void setData(const vtkSmartPointer<vtkPolyData> &data);
 
   void setRen(const vtkSmartPointer<vtkRenderer> &ren);
+
+  void setUVGrid(const std::shared_ptr<UVGrid> &uvGrid);
+
 private:
   vtkSmartPointer<vtkPolyData> data;
   vtkSmartPointer<vtkPoints> points;
   vtkSmartPointer<vtkRenderer> ren;
+  std::shared_ptr<UVGrid> uvGrid;
   vtkSmartPointer<vtkMatrix4x4> inverseCartographicProjection;
 
   void Execute(vtkObject *caller, unsigned long evId, void *callData) override;
+
   bool dragging = false;
 };
 

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;
+}

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;

vtk/src/layers/EGlyphLayer.cpp

@@ -11,43 +11,20 @@
 #include <vtkProperty.h>
 #include <vtkProperty2D.h>
 #include <vtkVertexGlyphFilter.h>
-#include <netcdf>
 #include <vtkArrowSource.h>
 #include "../CartographicTransformation.h"
+#include "../advection/readdata.h"
+#include "../advection/interpolate.h"
 
-using namespace netCDF;
 using namespace std;
 
-template <typename T>
+EGlyphLayer::EGlyphLayer(std::shared_ptr<UVGrid> uvGrid) {
-vector<T> getVarVector(const NcVar &var) {
-    int length = 1;
-    for (NcDim dim : var.getDims()) {
-        length *= dim.getSize();
-    }
-
-    vector<T> vec(length);
-
-    var.getVar(vec.data());
-
-    return vec;
-}
-
-tuple<vector<int>, vector<double>, vector<double>> readGrid() {
-    netCDF::NcFile data("../../../../data/grid.h5", netCDF::NcFile::read);
-    multimap< string, NcVar > vars = data.getVars();
-    vector<int> time = getVarVector<int>(vars.find("times")->second);
-    vector<double> longitude = getVarVector<double>(vars.find("longitude")->second);
-    vector<double> latitude = getVarVector<double>(vars.find("latitude")->second);
-
-    return {time, latitude, longitude};
-}
-
-
-EGlyphLayer::EGlyphLayer() {
   this->ren = vtkSmartPointer<vtkRenderer>::New();
   this->ren->SetLayer(1);
   this->ren->InteractiveOff();
 
+  this->uvGrid = uvGrid;
+
   this->data = vtkSmartPointer<vtkPolyData>::New();
   this->direction = vtkSmartPointer<vtkDoubleArray>::New();
   this->direction->SetName("direction");
@@ -57,31 +34,33 @@ EGlyphLayer::EGlyphLayer() {
 
 void EGlyphLayer::readCoordinates() {
   vtkNew<vtkPoints> points;
-  auto [times, lats, lons] = readGrid(); // FIXME: import Robin's readData function and use it
+  this->numLats = uvGrid->lats.size();
-  this->numLats = lats.size();
+  this->numLons = uvGrid->lons.size();
-  this->numLons = lons.size();
 
   this->direction->SetNumberOfComponents(3);
-  this->direction->SetNumberOfTuples(numLats*numLons); 
+  this->direction->SetNumberOfTuples(numLats * numLons);
-  points->Allocate(numLats*numLons);
+  points->Allocate(numLats * numLons);
 
   auto camera = createNormalisedCamera();
   ren->SetActiveCamera(camera);
 
   int i = 0;
-  for (double lat : lats) {
+  int latIndex = 0;
-    for (double lon : lons) {
+  for (double lat: uvGrid->lats) {
-        cout << "lon: " << lon << " lat: " << lat << endl;
+    int lonIndex = 0;
-      direction->SetTuple3(i, 0.45, 0.90, 0); //FIXME: read this info from file
+    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<vtkTransformFilter> transformFilter = createCartographicTransformFilter();
+  vtkSmartPointer<vtkTransformFilter> transformFilter = createCartographicTransformFilter(uvGrid);
   transformFilter->SetInputData(data);
 
   vtkNew<vtkGlyphSource2D> arrowSource;
@@ -109,16 +88,22 @@ void EGlyphLayer::readCoordinates() {
   vtkNew<vtkActor> 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++) {
+  int i = 0;
-    this->direction->SetTuple3(i, std::cos(t), std::sin(t), 0); // FIXME: fetch data from file.
+  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();
 }

vtk/src/layers/EGlyphLayer.h

@@ -2,8 +2,11 @@
 #define EGLYPHLAYER_H
 
 #include "Layer.h"
+#include <memory>
 #include <vtkPolyData.h>
 
+#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<vtkPolyData> data;
   vtkSmartPointer<vtkDoubleArray> direction;
+  std::shared_ptr<UVGrid> uvGrid;
   int numLats;
   int numLons;
 
@@ -22,7 +26,7 @@ private:
 public:
   /** Constructor.
     */
-  EGlyphLayer();
+  EGlyphLayer(std::shared_ptr<UVGrid> uvGrid);
 
   /** updates the glyphs to reflect the given timestamp in the dataset.
     * @param t : the time at which to fetch the data.

vtk/src/layers/LGlyphLayer.cpp

@@ -17,12 +17,12 @@
 
 #include "../CartographicTransformation.h"
 
-
 vtkSmartPointer<SpawnPointCallback> LGlyphLayer::createSpawnPointCallback() {
   auto newPointCallBack = vtkSmartPointer<SpawnPointCallback>::New();
   newPointCallBack->setData(data);
   newPointCallBack->setPoints(points);
   newPointCallBack->setRen(ren);
+  newPointCallBack->setUVGrid(uvGrid);
   return newPointCallBack;
 }
 
@@ -31,7 +31,7 @@ vtkSmartPointer<SpawnPointCallback> 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() {
+LGlyphLayer::LGlyphLayer(std::shared_ptr<UVGrid> uvGrid, std::unique_ptr<AdvectionKernel> advectionKernel) {
   this->ren = vtkSmartPointer<vtkRenderer>::New();
   this->ren->SetLayer(2);
 
@@ -39,12 +39,13 @@ LGlyphLayer::LGlyphLayer() {
   this->data = vtkSmartPointer<vtkPolyData>::New();
   this->data->SetPoints(this->points);
 
+  advector = std::move(advectionKernel);
+  this->uvGrid = uvGrid;
+
   auto camera = createNormalisedCamera();
   ren->SetActiveCamera(camera);
 
-    auto transform = createCartographicTransformFilter();
+  vtkSmartPointer<vtkTransformFilter> transformFilter = createCartographicTransformFilter(uvGrid);
-
-    vtkSmartPointer<vtkTransformFilter> transformFilter = createCartographicTransformFilter();
   transformFilter->SetInputData(data);
 
   vtkNew<vtkGlyphSource2D> circleSource;
@@ -70,35 +71,31 @@ LGlyphLayer::LGlyphLayer() {
 
 // 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(-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.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(-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<double, double> advect(int time, double lat, double lon) {
-    return {lat + 0.01, lon + 0.01};
-}
-
 void LGlyphLayer::updateData(int t) {
+  const int SUPERSAMPLINGRATE = 4;
   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]);
+    for (int i = 0; i < SUPERSAMPLINGRATE; i++) {
-        this->points->SetPoint(n, xNew, yNew, 0);
+      std::tie(point[1], point[0]) = advector->advect(t, point[1], point[0], (t-lastT)/SUPERSAMPLINGRATE);
     }
+    this->points->SetPoint(n, point[0], point[1], 0);
+  }
+  lastT = t;
   this->points->Modified();
 }
 
 void LGlyphLayer::addObservers(vtkSmartPointer<vtkRenderWindowInteractor> interactor) {
-    auto newPointCallBack = vtkSmartPointer<SpawnPointCallback>::New();
+  auto newPointCallBack = createSpawnPointCallback();
-    newPointCallBack->setData(data);
-    newPointCallBack->setPoints(points);
-    newPointCallBack->setRen(ren);
   interactor->AddObserver(vtkCommand::LeftButtonPressEvent, newPointCallBack);
   interactor->AddObserver(vtkCommand::LeftButtonReleaseEvent, newPointCallBack);
   interactor->AddObserver(vtkCommand::MouseMoveEvent, newPointCallBack);

vtk/src/layers/LGlyphLayer.h

@@ -2,6 +2,7 @@
 #define LGLYPHLAYER_H
 
 #include "Layer.h"
+#include "../advection/AdvectionKernel.h"
 #include "../commands/SpawnPointCallback.h"
 #include <vtkPolyData.h>
 #include <vtkInteractorStyle.h>
@@ -13,13 +14,14 @@ class LGlyphLayer : public Layer {
 private:
   vtkSmartPointer<vtkPoints> points;
   vtkSmartPointer<vtkPolyData> data;
-
+  std::unique_ptr<AdvectionKernel> advector;
-
+  std::shared_ptr<UVGrid> uvGrid;
+  int lastT = 1000;
 
 public:
   /** Constructor.
     */
-  LGlyphLayer();
+  LGlyphLayer(std::shared_ptr<UVGrid> uvGrid, std::unique_ptr<AdvectionKernel> advectionKernel);
 
   /** This function spoofs a few points in the dataset. Mostly used for testing.
     */

vtk/src/main.cpp

@@ -7,22 +7,30 @@
 #include <vtkProperty2D.h>
 #include <vtkRenderer.h>
 #include <vtkVertexGlyphFilter.h>
+#include <memory>
 
 #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();
+  cout << "reading data..." << endl;
-  l->spoofPoints();
+  shared_ptr<UVGrid> uvGrid = std::make_shared<UVGrid>();
+  auto kernelRK4 = make_unique<RK4AdvectionKernel>(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();