From b983023d49061848dce125512c6a13ae686eee41 Mon Sep 17 00:00:00 2001
From: Michael Uleysky <uleysky@poi.dvo.ru>
Date: Wed, 14 Feb 2024 13:37:10 +1000
Subject: [PATCH] New source GRIDVEL for reading pair of Surfer 6 grd files for
 velocity

---
 sources/GRIDVEL.cpp | 86 +++++++++++++++++++++++++++++++++++++++++++++
 sources/GRIDVEL.h   | 49 ++++++++++++++++++++++++++
 2 files changed, 135 insertions(+)
 create mode 100644 sources/GRIDVEL.cpp
 create mode 100644 sources/GRIDVEL.h

diff --git a/sources/GRIDVEL.cpp b/sources/GRIDVEL.cpp
new file mode 100644
index 0000000..fd7ec8c
--- /dev/null
+++ b/sources/GRIDVEL.cpp
@@ -0,0 +1,86 @@
+#define MICHLIB_NOSOURCE
+#include "GRIDVEL.h"
+
+MString GRIDVELData::Info() const
+{
+ if(!isOk()) return "";
+
+ // clang-format off
+ return MString() +
+  " Region: (" + u.Xb() + " : " + u.Xe() + ") x (" + u.Yb() + " : " + u.Ye() + ")\n" +
+  " Grid: " + u.Nx() + "x" + u.Ny() + "\n" +
+  " Supported variables: u, v, U2";
+ // clang-format on
+}
+
+MString GRIDVELData::Open(const CLArgs& args)
+{
+ MString uname = args.contains("u") ? args.at("u") : "";
+ MString vname = args.contains("v") ? args.at("v") : "";
+ if(!uname.Exist()) return "File with u component not specified";
+ if(!vname.Exist()) return "File with v component not specified";
+
+ unit = args.contains("unit") ? args.at("unit") : "cm/s";
+ fill = args.contains("fill") ? args.at("fill").ToReal() : 1e10;
+
+ u.Open(uname);
+ if(!u.Opened()) return "Can't open file: " + uname;
+ v.Open(vname);
+ if(!v.Opened()) return "Can't open file: " + vname;
+ if(u != v) return "Files " + uname + " and " + vname + " have different grids";
+
+ return "";
+}
+
+bool GRIDVELData::Read(const MString& vname, std::map<MString, GRIDVELData::Data>& cache, size_t i) const
+{
+ if(cache.contains(vname)) return true;
+ if(!isOk()) return false;
+
+ // Only rectangular grids are supported
+ real xs = (u.Xe() - u.Xb()) / (u.Nx() - 1);
+ real ys = (u.Ye() - u.Yb()) / (u.Ny() - 1);
+
+ Data out(u.Nx(), u.Ny(), u.Xb(), u.Yb(), xs, ys, MString(unit));
+
+ // U and U2 from u and v
+ if(vname == "U" || vname == "U2")
+ {
+  bool square = vname == "U2";
+  if(!(Read("u", cache, i) && Read("v", cache, i))) return false;
+  cache[vname]      = cache.at("u").CopyGrid();
+  auto&       U     = cache.at(vname);
+  const auto& udata = cache.at("u");
+  const auto& vdata = cache.at("v");
+  if(udata.Unit().Exist()) U.SetUnit(square ? ("(" + udata.Unit() + ")2") : udata.Unit());
+  for(size_t ind = 0; ind < U.N(); ind++)
+  {
+   if(udata.IsFill(ind) || vdata.IsFill(ind))
+    U.V(ind) = U.Fillval();
+   else
+    U.V(ind) = square ? (udata(ind) * udata(ind) + vdata(ind) * vdata(ind)) : michlib::Hypot(udata(ind), vdata(ind));
+  }
+  return true;
+ }
+
+ if(vname == "u" || vname == "v")
+ {
+  bool isu = vname == "u";
+  for(IType ix = 0; ix < u.Nx(); ix++)
+   for(IType iy = 0; iy < u.Ny(); iy++)
+   {
+    if(isu ? (u(ix, iy) == NAN || u(ix, iy) >= 1e10) : (v(ix, iy) == NAN || v(ix, iy) >= 1e10))
+     out.V(ix, iy) = out.Fillval();
+    else
+     out.V(ix, iy) = isu ? u(ix, iy) : v(ix, iy);
+   }
+
+  if(out)
+  {
+   cache[vname] = std::move(out);
+   return true;
+  }
+ }
+
+ return false;
+}
diff --git a/sources/GRIDVEL.h b/sources/GRIDVEL.h
new file mode 100644
index 0000000..8cddcd4
--- /dev/null
+++ b/sources/GRIDVEL.h
@@ -0,0 +1,49 @@
+#include "DataAdapters/gridfile.h"
+#include "mdatetime.h"
+#include "simple2ddata.h"
+
+using michlib::MDateTime;
+
+class GRIDVELData
+{
+ michlib::GridFile u,v;
+ MString unit;
+ real fill;
+
+ using IType=decltype(u.Nx());
+
+ public:
+ static constexpr const char* name = "GRIDVEL";
+
+ static constexpr const char* disabledactions = "genintfile";
+
+ using Data = Simple2DData;
+
+ GRIDVELData() = default;
+
+ MString Info() const;
+ // TODO: RetVal
+ MString Open(const CLArgs& args);
+
+ bool isOk() const { return u.Opened() &&v.Opened(); }
+
+ size_t NTimes() const { return 1; }
+
+ MDateTime Time(size_t i) const
+ {
+  return MDateTime();
+ }
+
+ //time_t Timestep() const { return isOk() ? (times[1] - times[0]) : 0; }
+
+ explicit operator bool() const { return u.Opened() &&v.Opened(); }
+
+ VarPresence CheckVar(const MString& vname) const
+ {
+  if(vname == "u" || vname == "v") return VarPresence::INTERNAL;
+  if(vname == "U" || vname == "U2") return VarPresence::DERIVED;
+  return VarPresence::NONE;
+ }
+
+ bool Read(const MString& vname, std::map<MString, Data>& cache, size_t i) const;
+};