Merge branch 'develop' into feature/interpolation_binning_01

doc/example-grids/classic_gaussian_2.yml

@@ -10,7 +10,7 @@ projection :
 
 check : 
   size : 1688
-  lonlat(first) : [3,44.8796]
-  lonlat(last)  : [-172.453,-54.9736]
-  uid : 64d609c6ee4b036b209047aef97a10eb
+  lonlat(first) : [3,49.1204]
+  lonlat(last)  : [179.6485,-38.8936]
+  uid : 02ef7ed866af557e04882342a90fddfe
   bounding_box(n,w,s,e) : [90,0,-90,360]

doc/example-grids/custom_structured_5.yml

@@ -10,5 +10,5 @@ check :
   size : 144
   lonlat(first) : [3,40]
   lonlat(last)  : [-177,-40]
-  uid : d55211f27c4a8c6ce94ea67f0c706e22
+  uid : 23f0a23085bbea10277b457f482927c0
   bounding_box(n,w,s,e) : [90,0,-90,360]

doc/example-grids/octahedral_gaussian_2.yml

@@ -7,7 +7,7 @@ projection : { type: "rotated_schmidt", stretching_factor : 4.0,  north_pole : [
 
 check : 
   size : 1600
-  lonlat(first) : [3,45.9397]
-  lonlat(last)  : [-165.776,-62.6128]
-  uid : 20308244515b7c78e22709ebaa842246
+  lonlat(first) : [3,48.0603]
+  lonlat(last)  : [177.0166,-30.8]
+  uid : 1297f9ab088b831da3139eb5f44e20a9
   bounding_box(n,w,s,e) : [90,0,-90,360]

doc/example-grids/octahedral_gaussian_4.yml

@@ -6,7 +6,7 @@ projection : { type: "rotated_schmidt", stretching_factor: 4.0,  north_pole: [3.
 
 check : 
   size : 1600
-  lonlat(first) : [3,45.9397]
-  lonlat(last)  : [-165.776,-62.6128]
-  uid : 20308244515b7c78e22709ebaa842246
+  lonlat(first) : [3,48.0603]
+  lonlat(last)  : [177.0166,-30.8]
+  uid : 1297f9ab088b831da3139eb5f44e20a9
   bounding_box(n,w,s,e) : [90,0,-90,360]

doc/example-grids/regional_rotated_mercator_1.yml

@@ -14,5 +14,5 @@ check :
   size : 2000
   lonlat(first) : [-10.319,40.2109]
   lonlat(last)  : [24.4206,57.2263]
-  uid : 9d4b8c38db08c6f3fe3221c8770d8d57
+  uid : 88d856331d46d30703b38e3ed2e8b2de
   bounding_box(n,w,s,e) : [58.734,-16.4206,40.2109,24.4206]

doc/example-grids/regular_gaussian_2.yml

@@ -10,7 +10,7 @@ projection :
 
 check : 
   size : 128
-  lonlat(first) : [3,38.8589]
-  lonlat(last)  : [-135.011,-72.4668]
-  uid : d179fd0b56811242a1a0a8e83dade6a1
+  lonlat(first) : [3,55.1411]
+  lonlat(last)  : [170.8438,-16.8513]
+  uid : 2efc83bfae617ec2fe1a7cebe523cf41
   bounding_box(n,w,s,e) : [90,0,-90,360]

src/atlas/field/FieldSet.cc

@@ -22,22 +22,40 @@ namespace field {
 //------------------------------------------------------------------------------------------------------
 
 void FieldSetImpl::FieldObserver::onFieldRename(FieldImpl& field) {
-    std::string name = field.name();
-    for (auto& kv: fieldset_.index_) {
-        const auto old_name = kv.first;
-        const auto idx      = kv.second;
-        if (&field == fieldset_.fields_[idx].get()) {
+
+    for (idx_t idx=0; idx<fieldset_.size(); ++idx) {
+        if (fieldset_[idx].get() == &field) {
+            std::string old_name = fieldset_.field_names_[idx];
+            std::string name = field.name();
             if (name.empty()) {
                 std::stringstream ss;
                 ss << fieldset_.name_ << "[" << idx << "]";
                 name = ss.str();
             }
             fieldset_.index_.erase(old_name);
             fieldset_.index_[name] = idx;
+            fieldset_.field_names_[idx] = name;
+
+            auto duplicate_exists = [&](const std::string& _name) {
+                return fieldset_.duplicates_.find(_name) != fieldset_.duplicates_.end();
+            };
+
+            if (duplicate_exists(old_name)) {
+                --fieldset_.duplicates_[old_name];
+                if (fieldset_.duplicates_[old_name] == 1) {
+                    std::size_t restored_index = std::find(fieldset_.field_names_.begin(), fieldset_.field_names_.end(), old_name) - fieldset_.field_names_.begin();
+                    fieldset_.index_[old_name] = restored_index;
+                }
+            }
+            if (duplicate_exists(name)) {
+                ++fieldset_.duplicates_[name];
+            }
+            else {
+                fieldset_.duplicates_[name] = 1;
+            }
             return;
         }
     }
-    throw_AssertionFailed("Should not be here",Here());
 }
 
 
@@ -52,19 +70,34 @@ void FieldSetImpl::clear() {
     }
     index_.clear();
     fields_.clear();
+    field_names_.clear();
+    duplicates_.clear();
 }
 
 Field FieldSetImpl::add(const Field& field) {
+
+    auto update_duplicates = [&](const std::string& name) {
+        if (duplicates_.find(name) != duplicates_.end()) {
+            ++duplicates_[name];
+        }
+        else {
+            duplicates_[name] = 1;
+        }
+    };
+
+    std::string name;
     if (field.name().size()) {
-        index_[field.name()] = size();
+        name = field.name();
     }
     else {
-        std::stringstream name;
-        name << name_ << "[" << size() << "]";
-        index_[name.str()] = size();
+        std::stringstream name_ss;
+        name_ss << name_ << "[" << size() << "]";
+        name = name_ss.str();
     }
+    index_[name] = size();
     fields_.push_back(field);
-
+    field_names_.push_back(name);
+    update_duplicates(name);
     field.get()->attachObserver(field_observer_);
     return field;
 }
@@ -76,9 +109,16 @@ bool FieldSetImpl::has(const std::string& name) const {
 
 Field& FieldSetImpl::field(const std::string& name) const {
     if (!has(name)) {
-        const std::string msg("FieldSet" + (name_.length() ? " \"" + name_ + "\"" : "") + ": cannot find field \"" +
-                              name + "\"");
-        throw_Exception(msg, Here());
+        std::stringstream msg;
+        msg << "FieldSet" << (name_.length() ? " \"" + name_ + "\"" : "") << ": cannot find field \"" << name << "\"";
+        throw_Exception(msg.str(), Here());
+    }
+    if (duplicates_.at(name) > 1) {
+        std::stringstream msg;
+        msg << "FieldSet" << (name_.length() ? " \"" + name_ + "\"" : "") << ": cannot get field with ambiguous name \n" << name << "\". "
+            << duplicates_.at(name) << " fields are registered with same name. Access field by index or iterator instead.";
+        throw_Exception(msg.str(), Here());
+
     }
     return const_cast<Field&>(fields_[index_.at(name)]);
 }
@@ -101,14 +141,8 @@ void FieldSetImpl::set_dirty(bool value) const {
     }
 }
 
-std::vector<std::string> FieldSetImpl::field_names() const {
-    std::vector<std::string> ret;
-
-    for (const_iterator field = cbegin(); field != cend(); ++field) {
-        ret.push_back(field->name());
-    }
-
-    return ret;
+const std::vector<std::string>& FieldSetImpl::field_names() const {
+    return field_names_;
 }
 
 //-----------------------------------------------------------------------------

src/atlas/field/FieldSet.h

@@ -114,7 +114,7 @@ class FieldSetImpl : public util::Object {
         return fields_[i];
     }
 
-    std::vector<std::string> field_names() const;
+    const std::vector<std::string>& field_names() const;
 
     Field add(const Field&);
 
@@ -141,6 +141,8 @@ class FieldSetImpl : public util::Object {
     std::string name_;                    ///< internal name
     util::Metadata metadata_;             ///< metadata associated with the FieldSet
     std::map<std::string, idx_t> index_;  ///< name-to-index map, to refer fields by name
+    std::vector<std::string>     field_names_;  ///< field names
+    std::map<std::string, idx_t> duplicates_;  ///< name-to-duplicates map, to refer fields by name
 
     friend class FieldObserver;
     FieldObserver field_observer_;
@@ -234,7 +236,7 @@ class FieldSet : DOXYGEN_HIDE(public util::ObjectHandle<field::FieldSetImpl>) {
         return get()->field(i);
     }
 
-    std::vector<std::string> field_names() const { return get()->field_names(); }
+    const std::vector<std::string>& field_names() const { return get()->field_names(); }
 
     Field add(const Field& field) { return get()->add(field); }
 

src/atlas/interpolation/method/cubedsphere/CubedSphereBilinear.cc

@@ -33,6 +33,10 @@ void CubedSphereBilinear::do_setup(const FunctionSpace& source, const FunctionSp
     ATLAS_ASSERT(ncSource);
     ATLAS_ASSERT(target_);
 
+    // Enable or disable halo exchange.
+    this->allow_halo_exchange_ = halo_exchange_;
+
+
     // return early if no output points on this partition reserve is called on
     // the triplets but also during the sparseMatrix constructor. This won't
     // work for empty matrices
@@ -46,9 +50,6 @@ void CubedSphereBilinear::do_setup(const FunctionSpace& source, const FunctionSp
     const auto N         = CubedSphereGrid(ncSource.mesh().grid()).N();
     const auto tolerance = 2. * std::numeric_limits<double>::epsilon() * N;
 
-    // Enable or disable halo exchange.
-    this->allow_halo_exchange_ = halo_exchange_;
-
     // Loop over target at calculate interpolation weights.
     auto weights          = std::vector<Triplet>{};
     const auto ghostView  = array::make_view<int, 1>(target_.ghost());

src/atlas/util/Rotation.cc

@@ -126,9 +126,9 @@ void Rotation::precompute() {
 
 Rotation::Rotation(const PointLonLat& south_pole, double rotation_angle) {
     spole_ = south_pole;
-    npole_ = PointLonLat(spole_.lon() - 180., spole_.lat() + 180.);
-    if (npole_.lat() > 90) {
-        npole_.lon() += 180.;
+    npole_ = PointLonLat(spole_.lon() - 180., -spole_.lat());
+    if (npole_.lon() < 0.) {
+        npole_.lon() += 360.;
     }
     angle_ = wrap_angle(rotation_angle);
 
@@ -144,16 +144,16 @@ Rotation::Rotation(const eckit::Parametrisation& p) {
     std::vector<double> pole(2);
     if (p.get("north_pole", pole)) {
         npole_ = PointLonLat(pole.data());
-        spole_ = PointLonLat(npole_.lon() + 180., npole_.lat() - 180.);
-        if (spole_.lat() < -90) {
-            spole_.lon() -= 180.;
+        spole_ = PointLonLat(npole_.lon() - 180., -npole_.lat());
+        if (spole_.lon() < 0.) {
+            spole_.lon() += 360.;
         }
     }
     else if (p.get("south_pole", pole)) {
         spole_ = PointLonLat(pole.data());
-        npole_ = PointLonLat(spole_.lon() - 180., spole_.lat() + 180.);
-        if (npole_.lat() > 90) {
-            npole_.lon() += 180.;
+        npole_ = PointLonLat(spole_.lon() - 180., -spole_.lat());
+        if (npole_.lon() < 0.) {
+            npole_.lon() += 360.;
         }
     }
 

src/atlas_f/CMakeLists.txt

@@ -96,6 +96,7 @@ generate_fortran_bindings(FORTRAN_BINDINGS ../atlas/mesh/actions/BuildPeriodicBo
 generate_fortran_bindings(FORTRAN_BINDINGS ../atlas/mesh/actions/BuildHalo.h)
 generate_fortran_bindings(FORTRAN_BINDINGS ../atlas/mesh/actions/BuildEdges.h)
 generate_fortran_bindings(FORTRAN_BINDINGS ../atlas/mesh/actions/BuildDualMesh.h)
+generate_fortran_bindings(FORTRAN_BINDINGS ../atlas/mesh/actions/BuildNode2CellConnectivity.h)
 generate_fortran_bindings(FORTRAN_BINDINGS ../atlas/mesh/actions/WriteLoadBalanceReport.h)
 generate_fortran_bindings(FORTRAN_BINDINGS ../atlas/meshgenerator/detail/MeshGeneratorInterface.h
     MODULE atlas_meshgenerator_c_binding

src/atlas_f/atlas_module.F90

@@ -141,6 +141,7 @@ module atlas_module
     & atlas_build_halo, &
     & atlas_build_edges, &
     & atlas_build_pole_edges, &
+    & atlas_build_node_to_cell_connectivity, &
     & atlas_build_node_to_edge_connectivity, &
     & atlas_build_median_dual_mesh, &
     & atlas_write_load_balance_report, &

src/atlas_f/mesh/atlas_mesh_actions_module.F90

@@ -19,6 +19,7 @@ module atlas_mesh_actions_module
 public :: atlas_build_halo
 public :: atlas_build_edges
 public :: atlas_build_pole_edges
+public :: atlas_build_node_to_cell_connectivity
 public :: atlas_build_node_to_edge_connectivity
 public :: atlas_build_median_dual_mesh
 public :: atlas_write_load_balance_report
@@ -84,6 +85,13 @@ subroutine atlas_build_pole_edges(mesh)
   call atlas__build_pole_edges(mesh%CPTR_PGIBUG_A)
 end subroutine atlas_build_pole_edges
 
+subroutine atlas_build_node_to_cell_connectivity(mesh)
+  use atlas_BuildNode2CellConnectivity_c_binding
+  use atlas_Mesh_module, only: atlas_Mesh
+  type(atlas_Mesh), intent(inout) :: mesh
+  call atlas__build_node_to_cell_connectivity(mesh%CPTR_PGIBUG_A)
+end subroutine atlas_build_node_to_cell_connectivity
+
 subroutine atlas_build_node_to_edge_connectivity(mesh)
   use atlas_BuildEdges_c_binding
   use atlas_Mesh_module, only: atlas_Mesh

src/tests/field/test_fieldset.cc

@@ -59,6 +59,33 @@ CASE("test_rename") {
 
 }
 
+CASE("test_duplicate_name_throws") {
+    FieldSet fieldset;
+    auto field_0 = fieldset.add(Field("0", make_datatype<double>(), array::make_shape(10,4)));
+    auto field_1 = fieldset.add(Field("0", make_datatype<double>(), array::make_shape(10,5))); // same name as field_0, uh-oh !
+    auto field_2 = fieldset.add(Field("2", make_datatype<double>(), array::make_shape(10,6)));
+
+    Field f;
+    EXPECT_NO_THROW(f = fieldset["2"]); // OK
+    EXPECT_EQ(f.shape(1), 6);
+
+    EXPECT_THROWS(f = fieldset["0"]);   // ambiguous because field_0 and field_1 have same name, should throw
+    field_1.rename("1");                // fix ambiguity between field_0 and field_1
+    EXPECT_NO_THROW(f = fieldset["0"]); // no longer ambiguous
+    EXPECT_EQ(f.shape(1), 4);           // to be sure that we got the right field
+    EXPECT_NO_THROW(f = fieldset["1"]); // no longer ambiguous
+    EXPECT_EQ(f.shape(1), 5);           // to be sure that we got the right field
+
+    field_2.rename("0");                // Introduce new ambiguity between field_0 and field_2
+    EXPECT_THROWS(f = fieldset["0"]);   // ambiguous because field_0 and field_2 have same name, should throw
+    field_2.rename("2");                // fix ambiguity
+    EXPECT_NO_THROW(f = fieldset["0"]); // no longer ambiguous
+    EXPECT_EQ(f.shape(1), 4);           // to be sure we got the right field
+    EXPECT_NO_THROW(f = fieldset["2"]); // no longer ambiguous
+    EXPECT_EQ(f.shape(1), 6);           // to be sure we got the right field
+}
+
+
 //-----------------------------------------------------------------------------
 
 }  // namespace test