initial rough working

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

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/advection/AdvectionKernel.h

@@ -11,7 +11,7 @@
  */
 class AdvectionKernel {
 public:
-    const static int DT = 60 * 15; // 60 sec/min * 15 mins
+    const static int DT = 15 * 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.

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/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>
-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() {
+EGlyphLayer::EGlyphLayer(std::shared_ptr<UVGrid> uvGrid) {
   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,25 +34,28 @@ EGlyphLayer::EGlyphLayer() {
 
 void EGlyphLayer::readCoordinates() {
   vtkNew<vtkPoints> 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, 2*v, 2*u, 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);
@@ -94,8 +74,8 @@ void EGlyphLayer::readCoordinates() {
   glyph2D->SetInputConnection(transformFilter->GetOutputPort());
   glyph2D->OrientOn();
   glyph2D->ClampingOn();
-  glyph2D->SetScaleModeToScaleByVector(); 
-  glyph2D->SetVectorModeToUseVector(); 
+  glyph2D->SetScaleModeToScaleByVector();
+  glyph2D->SetVectorModeToUseVector();
   glyph2D->Update();
 
 //  vtkNew<vtkCoordinate> coordinate;
@@ -109,16 +89,21 @@ 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++) {
-    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];
+      this->direction->SetTuple3(i, 5*v, 5*u, 0); // FIXME: fetch data from file.
+      i++;
+    }
   }
   this->direction->Modified();
 }

vtk/src/layers/EGlyphLayer.h

@@ -4,6 +4,8 @@
 #include "Layer.h"
 #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 +13,7 @@ class EGlyphLayer : public Layer {
 private:
   vtkSmartPointer<vtkPolyData> data;
   vtkSmartPointer<vtkDoubleArray> direction;
+  std::shared_ptr<UVGrid> uvGrid;
   int numLats;
   int numLons;
 
@@ -22,7 +25,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

@@ -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::unique_ptr<AdvectionKernel> advectionKernel) {
     this->ren = vtkSmartPointer<vtkRenderer>::New();
     this->ren->SetLayer(2);
 
@@ -39,6 +39,8 @@ LGlyphLayer::LGlyphLayer() {
     this->data = vtkSmartPointer<vtkPolyData>::New();
     this->data->SetPoints(this->points);
 
+    advector = std::move(advectionKernel);
+
     auto camera = createNormalisedCamera();
     ren->SetActiveCamera(camera);
 
@@ -88,7 +90,7 @@ 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]);
+        auto [yNew, xNew] = advector->advect(t, point[1], point[0]);
         this->points->SetPoint(n, xNew, yNew, 0);
     }
     this->points->Modified();

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;
 
 
 
 public:
   /** Constructor.
     */
-  LGlyphLayer();
+  LGlyphLayer(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,28 @@
 #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;
 
 
 int main() {
-  auto l = new LGlyphLayer();
-  l->spoofPoints();
+  shared_ptr<UVGrid> uvGrid = std::make_shared<UVGrid>();
+  auto kernelRK4 = make_unique<RK4AdvectionKernel>(uvGrid);
+
+  auto l = new LGlyphLayer(move(kernelRK4));
+//  l->spoofPoints();
 
   Program *program = new Program();
   program->addLayer(new BackgroundImage("../../../../data/map_661-661.png"));
-  program->addLayer(new EGlyphLayer());
+  program->addLayer(new EGlyphLayer(uvGrid));
   program->addLayer(l);
 
   program->render();