From d2449492ad250776d07c2931bb30a098ac054f2d Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Tue, 8 Feb 2022 17:32:40 +0000 Subject: [PATCH 1/3] ATLAS-344 Fix PolygonXY::contains, requiring tolerances at edges --- .../MatchingMeshPartitionerLonLatPolygon.cc | 20 +++++++++++++------ src/atlas/util/PolygonXY.cc | 9 +++++++-- 2 files changed, 21 insertions(+), 8 deletions(-) diff --git a/src/atlas/grid/detail/partitioner/MatchingMeshPartitionerLonLatPolygon.cc b/src/atlas/grid/detail/partitioner/MatchingMeshPartitionerLonLatPolygon.cc index dee4a9ac4..5dbc571fd 100644 --- a/src/atlas/grid/detail/partitioner/MatchingMeshPartitionerLonLatPolygon.cc +++ b/src/atlas/grid/detail/partitioner/MatchingMeshPartitionerLonLatPolygon.cc @@ -76,9 +76,13 @@ void MatchingMeshPartitionerLonLatPolygon::partition(const Grid& grid, int parti int min = compute(east - 360.); constexpr double eps = 1.e-10; - if (min < 0 && east - west > 360. + eps) { - min = compute(west - eps); - } + bool second_try = [&]() { + if (min < 0 && east - west > 360. + eps) { + min = compute(west - eps); + return true; + } + return false; + }(); if (min < 0) { size_t i = 0; size_t max_failures = grid.size(); @@ -95,10 +99,14 @@ void MatchingMeshPartitionerLonLatPolygon::partition(const Grid& grid, int parti } size_t nb_failures = failed_index.size(); std::stringstream err; + err.precision(20); err << "Could not find partition of " << nb_failures - << " target grid points (source " - "mesh does not contain all target grid points)\n" - "Failed target grid points with global index:\n"; + << " target grid points (source mesh does not contain all target grid points)\n" + << "Tried first normalizing coordinates with west=" << east - 360.; + if (second_try) { + err << "Tried second time normalizing coordinates with west=" << west - eps << "\n"; + } + err << "Failed target grid points with global index:\n"; for (size_t n = 0; n < nb_failures; ++n) { err << " - " << std::setw(10) << std::left << failed_index[n] + 1 << " {lon,lat} : " << failed_lonlat[n] << "\n"; diff --git a/src/atlas/util/PolygonXY.cc b/src/atlas/util/PolygonXY.cc index 190188904..71735becb 100644 --- a/src/atlas/util/PolygonXY.cc +++ b/src/atlas/util/PolygonXY.cc @@ -101,9 +101,10 @@ bool PolygonXY::contains(const Point2& P) const { return dx * dx + dy * dy; }; + constexpr double eps = 1.e-10; // check first bounding box - if (coordinatesMax_.y() < P.y() || P.y() < coordinatesMin_.y() || coordinatesMax_.x() < P.x() || - P.x() < coordinatesMin_.x()) { + if (coordinatesMax_.y() + eps < P.y() || P.y() < coordinatesMin_.y() - eps || coordinatesMax_.x() + eps < P.x() || + P.x() < coordinatesMin_.x() - eps) { return false; } @@ -137,6 +138,10 @@ bool PolygonXY::contains(const Point2& P) const { if (APB != BPA) { const double side = cross_product_analog(P, A, B); + bool on_edge = std::abs(side) < eps; + if (on_edge) { + return true; + } if (APB && side > 0) { ++wn; } From 09cde23dbc49e1711973c02c2abacad873511781 Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Thu, 17 Feb 2022 10:20:39 +0000 Subject: [PATCH 2/3] Add macro to check eckit version --- src/atlas/library/config.h | 4 +++- src/atlas/util/detail/KDTree.h | 2 +- src/tests/util/test_kdtree.cc | 2 +- 3 files changed, 5 insertions(+), 3 deletions(-) diff --git a/src/atlas/library/config.h b/src/atlas/library/config.h index 08fa218b8..4b43fe9c8 100644 --- a/src/atlas/library/config.h +++ b/src/atlas/library/config.h @@ -46,7 +46,9 @@ using idx_t = long; /// Integer type for unique indices typedef gidx_t uidx_t; -#if ATLAS_ECKIT_VERSION_INT >= 11900 // eckit >= 1.19 +#define ATLAS_ECKIT_VERSION_AT_LEAST(x, y, z) (ATLAS_ECKIT_VERSION_INT >= x * 10000 + y * 100 + z) + +#if ATLAS_ECKIT_VERSION_AT_LEAST(1, 19, 0) #define ATLAS_ECKIT_HAVE_ECKIT_585 1 #else #define ATLAS_ECKIT_HAVE_ECKIT_585 0 diff --git a/src/atlas/util/detail/KDTree.h b/src/atlas/util/detail/KDTree.h index 46ee97ce5..086edebd1 100644 --- a/src/atlas/util/detail/KDTree.h +++ b/src/atlas/util/detail/KDTree.h @@ -309,7 +309,7 @@ class KDTree_eckit : public KDTreeBase { } idx_t size() const override { -#if ATLAS_ECKIT_VERSION_INT >= 11302 // v1.13.2 +#if ATLAS_ECKIT_VERSION_AT_LEAST(1, 13, 2) return static_cast(tree_->size()); #else // Assume ECKIT-515 not implemented. diff --git a/src/tests/util/test_kdtree.cc b/src/tests/util/test_kdtree.cc index 37e28e50a..782a0f7b6 100644 --- a/src/tests/util/test_kdtree.cc +++ b/src/tests/util/test_kdtree.cc @@ -153,7 +153,7 @@ static const IndexKDTree& search() { } // namespace //------------------------------------------------------------------------------------------------ -static bool ECKIT_515_implemented = (ATLAS_ECKIT_VERSION_INT >= 11302); // version 1.13.2 +static bool ECKIT_515_implemented = ATLAS_ECKIT_VERSION_AT_LEAST(1, 13, 2); // --> implements eckit::KDTree::size() and eckit::KDTree::empty() CASE("test kdtree") { From e4b12c728da736acf32976c003b0fe97d37649d3 Mon Sep 17 00:00:00 2001 From: odlomax Date: Thu, 17 Feb 2022 09:58:58 +0000 Subject: [PATCH 3/3] Added methods to projection::Jacobian (#93) --- src/atlas/projection/detail/ProjectionImpl.h | 127 ++++++++----------- src/tests/projection/test_jacobian.cc | 55 ++++++++ 2 files changed, 106 insertions(+), 76 deletions(-) diff --git a/src/atlas/projection/detail/ProjectionImpl.h b/src/atlas/projection/detail/ProjectionImpl.h index 655134b11..45ccce5ef 100644 --- a/src/atlas/projection/detail/ProjectionImpl.h +++ b/src/atlas/projection/detail/ProjectionImpl.h @@ -10,9 +10,12 @@ #pragma once +#include +#include #include #include +#include "atlas/runtime/Exception.h" #include "atlas/util/Factory.h" #include "atlas/util/NormaliseLongitude.h" #include "atlas/util/Object.h" @@ -41,97 +44,68 @@ class ProjectionImpl : public util::Object { using Spec = atlas::util::Config; class Jacobian : public std::array, 2> { public: - static Jacobian identity() { - Jacobian id; - id[0] = {1.0, 0.0}; - id[1] = {0.0, 1.0}; - return id; - } + using std::array, 2>::array; - Jacobian inverse() const { - const Jacobian& jac = *this; - Jacobian inv; - double det = jac[0][0] * jac[1][1] - jac[0][1] * jac[1][0]; - inv[0][0] = +jac[1][1] / det; - inv[0][1] = -jac[0][1] / det; - inv[1][0] = -jac[1][0] / det; - inv[1][1] = +jac[0][0] / det; - return inv; - }; + static Jacobian identity() { return Jacobian{1., 0., 0., 1.}; } - Jacobian transpose() const { - Jacobian tra = *this; - std::swap(tra[0][1], tra[1][0]); - return tra; - }; + Jacobian(double j00, double j01, double j10, double j11): + array{std::array{j00, j01}, std::array{j10, j11}} {} - Jacobian operator-(const Jacobian& jac2) const { - const Jacobian& jac1 = *this; - Jacobian jac; - jac[0][0] = jac1[0][0] - jac2[0][0]; - jac[0][1] = jac1[0][1] - jac2[0][1]; - jac[1][0] = jac1[1][0] - jac2[1][0]; - jac[1][1] = jac1[1][1] - jac2[1][1]; - return jac; - } + Jacobian(std::initializer_list> list): + Jacobian{*(list.begin()->begin()), *(list.begin()->begin() + 1), *((list.begin() + 1)->begin()), + *((list.begin() + 1)->begin() + 1)} {} - Jacobian operator+(const Jacobian& jac2) const { - const Jacobian& jac1 = *this; - Jacobian jac; - jac[0][0] = jac1[0][0] + jac2[0][0]; - jac[0][1] = jac1[0][1] + jac2[0][1]; - jac[1][0] = jac1[1][0] + jac2[1][0]; - jac[1][1] = jac1[1][1] + jac2[1][1]; - return jac; + Jacobian operator-(const Jacobian& jac) const { + return Jacobian{(*this)[0][0] - jac[0][0], (*this)[0][1] - jac[0][1], (*this)[1][0] - jac[1][0], + (*this)[1][1] - jac[1][1]}; } - double norm() const { - const Jacobian& jac = *this; - return sqrt(jac[0][0] * jac[0][0] + jac[0][1] * jac[0][1] + jac[1][0] * jac[1][0] + jac[1][1] * jac[1][1]); + Jacobian operator+(const Jacobian& jac) const { + return Jacobian{(*this)[0][0] + jac[0][0], (*this)[0][1] + jac[0][1], (*this)[1][0] + jac[1][0], + (*this)[1][1] + jac[1][1]}; } - Jacobian operator*(const Jacobian& jac2) const { - const Jacobian& jac1 = *this; - Jacobian jac; - jac[0][0] = jac1[0][0] * jac2[0][0] + jac1[0][1] * jac2[1][0]; - jac[0][1] = jac1[0][0] * jac2[0][1] + jac1[0][1] * jac2[1][1]; - jac[1][0] = jac1[1][0] * jac2[0][0] + jac1[1][1] * jac2[1][0]; - jac[1][1] = jac1[1][0] * jac2[0][1] + jac1[1][1] * jac2[1][1]; - return jac; + Jacobian operator*(double a) const { + return Jacobian{(*this)[0][0] * a, (*this)[0][1] * a, (*this)[1][0] * a, (*this)[1][1] * a}; } - double dx_dlon() const { - const Jacobian& jac = *this; - return jac[JDX][JDLON]; - } - double dy_dlon() const { - const Jacobian& jac = *this; - return jac[JDY][JDLON]; - } - double dx_dlat() const { - const Jacobian& jac = *this; - return jac[JDX][JDLAT]; - } - double dy_dlat() const { - const Jacobian& jac = *this; - return jac[JDY][JDLAT]; + Jacobian operator*(const Jacobian& jac) const { + return Jacobian{(*this)[0][0] * jac[0][0] + (*this)[0][1] * jac[1][0], + (*this)[0][0] * jac[0][1] + (*this)[0][1] * jac[1][1], + (*this)[1][0] * jac[0][0] + (*this)[1][1] * jac[1][0], + (*this)[1][0] * jac[0][1] + (*this)[1][1] * jac[1][1]}; } - double dlon_dx() const { - const Jacobian& jac = *this; - return jac[JDLON][JDX]; + Point2 operator*(const Point2& x) const { + return Point2{(*this)[0][0] * x[0] + (*this)[0][1] * x[1], (*this)[1][0] * x[0] + (*this)[1][1] * x[1]}; } - double dlon_dy() const { - const Jacobian& jac = *this; - return jac[JDLON][JDY]; + + double norm() const { + return std::sqrt((*this)[0][0] * (*this)[0][0] + (*this)[0][1] * (*this)[0][1] + + (*this)[1][0] * (*this)[1][0] + (*this)[1][1] * (*this)[1][1]); } - double dlat_dx() const { - const Jacobian& jac = *this; - return jac[JDLAT][JDX]; + + double determinant() const { return (*this)[0][0] * (*this)[1][1] - (*this)[0][1] * (*this)[1][0]; } + + Jacobian inverse() const { + return Jacobian{(*this)[1][1], -(*this)[0][1], -(*this)[1][0], (*this)[0][0]} * (1. / determinant()); } - double dlat_dy() const { - const Jacobian& jac = *this; - return jac[JDLAT][JDY]; + + Jacobian transpose() const { return Jacobian{(*this)[1][1], (*this)[0][1], (*this)[1][0], (*this)[0][0]}; } + + double dx_dlon() const { return (*this)[JDX][JDLON]; } + double dy_dlon() const { return (*this)[JDY][JDLON]; } + double dx_dlat() const { return (*this)[JDX][JDLAT]; } + double dy_dlat() const { return (*this)[JDY][JDLAT]; } + + double dlon_dx() const { return (*this)[JDLON][JDX]; } + double dlon_dy() const { return (*this)[JDLON][JDY]; } + double dlat_dx() const { return (*this)[JDLAT][JDX]; } + double dlat_dy() const { return (*this)[JDLAT][JDY]; } + + friend std::ostream& operator<<(std::ostream& os, const Jacobian& jac) { + os << jac[0][0] << " " << jac[0][1] << "\n" << jac[1][0] << " " << jac[1][1]; + return os; } private: @@ -147,6 +121,7 @@ class ProjectionImpl : public util::Object { }; }; + public: static const ProjectionImpl* create(const eckit::Parametrisation& p); static const ProjectionImpl* create(const std::string& type, const eckit::Parametrisation& p); diff --git a/src/tests/projection/test_jacobian.cc b/src/tests/projection/test_jacobian.cc index 9206998ec..3abfa774e 100644 --- a/src/tests/projection/test_jacobian.cc +++ b/src/tests/projection/test_jacobian.cc @@ -163,6 +163,61 @@ CASE("test_lonlat") { doTaylorTest(StructuredGrid("N16"), 1e-9); } +CASE("test_constructors") { + Projection::Jacobian a = {}; + Projection::Jacobian b = {1., 2., 3., 4.}; + Projection::Jacobian c = {{5., 6.}, {7., 8.}}; + + Log::info() << "a =" << std::endl; + Log::info() << a << std::endl; + Log::info() << "b =" << std::endl; + Log::info() << b << std::endl; + Log::info() << "c =" << std::endl; + Log::info() << c << std::endl; + + EXPECT_EQ(a[0][0], 0.); + EXPECT_EQ(a[0][1], 0.); + EXPECT_EQ(a[1][0], 0.); + EXPECT_EQ(a[1][1], 0.); + + EXPECT_EQ(b[0][0], 1.); + EXPECT_EQ(b[0][1], 2.); + EXPECT_EQ(b[1][0], 3.); + EXPECT_EQ(b[1][1], 4.); + + EXPECT_EQ(c[0][0], 5.); + EXPECT_EQ(c[0][1], 6.); + EXPECT_EQ(c[1][0], 7.); + EXPECT_EQ(c[1][1], 8.); +} + +CASE("test_multiplication") { + const Projection::Jacobian A = {{0., -1.}, {1., 0.}}; + const Point2 x = {2., 1.}; + + const auto Ax = A * x; + const auto AinvAx = A.inverse() * Ax; + const auto Atimes2 = A * 2.; + + Log::info() << "A =" << std::endl; + Log::info() << A << std::endl; + Log::info() << "2A =" << std::endl; + Log::info() << Atimes2 << std::endl; + Log::info() << "x =" << std::endl; + Log::info() << x << std::endl; + Log::info() << "Ax =" << std::endl; + Log::info() << Ax << std::endl; + Log::info() << "A^-1 Ax =" << std::endl; + Log::info() << AinvAx << std::endl; + + EXPECT_EQ(Atimes2[0][0], 0.); + EXPECT_EQ(Atimes2[0][1], -2.); + EXPECT_EQ(Atimes2[1][0], 2.); + EXPECT_EQ(Atimes2[1][1], 0.); + EXPECT_EQ(Ax, Point2(-1., 2.)); + EXPECT_EQ(AinvAx, x); +} + //-----------------------------------------------------------------------------