Add regional interpolation

src/atlas/CMakeLists.txt

@@ -533,6 +533,7 @@ interpolation/Vector2D.cc
 interpolation/Vector2D.h
 interpolation/Vector3D.cc
 interpolation/Vector3D.h
+interpolation/element/InterpElement.h
 interpolation/element/Quad2D.h
 interpolation/element/Quad2D.cc
 interpolation/element/Quad3D.cc
@@ -655,6 +656,8 @@ interpolation/method/structured/QuasiCubic2D.cc
 interpolation/method/structured/QuasiCubic2D.h
 interpolation/method/structured/QuasiCubic3D.cc
 interpolation/method/structured/QuasiCubic3D.h
+interpolation/method/structured/RegionalLinear2D.cc
+interpolation/method/structured/RegionalLinear2D.h
 interpolation/method/structured/StructuredInterpolation2D.h
 interpolation/method/structured/StructuredInterpolation2D.tcc
 interpolation/method/structured/StructuredInterpolation3D.h

src/atlas/functionspace/detail/StructuredColumns_setup.cc

@@ -588,13 +588,20 @@ void StructuredColumns::setup(const grid::Distribution& distribution, const ecki
 
         atlas_omp_parallel_for(idx_t n = 0; n < gridpoints.size(); ++n) {
             const GridPoint& gp = gridpoints[n];
-            if (gp.j >= 0 && gp.j < grid_->ny()) {
-                xy(gp.r, XX) = grid_->x(gp.i, gp.j);
-                xy(gp.r, YY) = grid_->y(gp.j);
-            }
-            else {
-                xy(gp.r, XX) = compute_x(gp.i, gp.j);
-                xy(gp.r, YY) = compute_y(gp.j);
+            if (regional) {
+              std::vector<double> lonlatVec(2);
+              grid_->lonlat(gp.i, gp.j, lonlatVec.data());
+              xy(gp.r, XX) = lonlatVec[0];
+              xy(gp.r, YY) = lonlatVec[1];
+            } else {
+              if (gp.j >= 0 && gp.j < grid_->ny()) {
+                  xy(gp.r, XX) = grid_->x(gp.i, gp.j);
+                  xy(gp.r, YY) = grid_->y(gp.j);
+              }
+              else {
+                  xy(gp.r, XX) = compute_x(gp.i, gp.j);
+                  xy(gp.r, YY) = compute_y(gp.j);
+              }
             }
 
             bool in_domain(false);

src/atlas/interpolation/element/InterpElement.h

@@ -0,0 +1,39 @@
+/*
+ * (C) Copyright 2024 Meteorlogisk Institutt
+ *
+ * This software is licensed under the terms of the Apache Licence Version 2.0
+ * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
+ */
+
+#pragma once
+
+#include <string>
+#include <utility>
+#include <vector>
+
+namespace atlas {
+namespace interpolation {
+namespace element {
+
+//----------------------------------------------------------------------------------------------------------------------
+
+class InterpElement {
+ public:
+  // Constructor
+  explicit InterpElement(const std::vector<std::pair<size_t, double>> & operations)
+  : operations_(operations) {}
+
+  // Accessors
+  const std::vector<std::pair<size_t, double>> & operations() const
+    {return operations_;}
+
+ private:
+  // Operations
+  std::vector<std::pair<size_t, double>> operations_;
+};
+
+//----------------------------------------------------------------------------------------------------------------------
+
+}  // namespace element
+}  // namespace interpolation
+}  // namespace atlas

src/atlas/interpolation/method/MethodFactory.cc

@@ -24,6 +24,7 @@
 #include "structured/Linear3D.h"
 #include "structured/QuasiCubic2D.h"
 #include "structured/QuasiCubic3D.h"
+#include "structured/RegionalLinear2D.h"
 #include "unstructured/FiniteElement.h"
 #include "unstructured/UnstructuredBilinearLonLat.h"
 
@@ -46,6 +47,7 @@ void force_link() {
             MethodBuilder<method::Cubic3D>();
             MethodBuilder<method::QuasiCubic2D>();
             MethodBuilder<method::QuasiCubic3D>();
+            MethodBuilder<method::RegionalLinear2D>();
             MethodBuilder<method::GridBoxAverage>();
             MethodBuilder<method::GridBoxMaximum>();
             MethodBuilder<method::CubedSphereBilinear>();