From bdb56d88f8fe52fa1715db31b6db1949820b4dd0 Mon Sep 17 00:00:00 2001
From: John Omotani <john.omotani@ukaea.uk>
Date: Tue, 20 Nov 2018 14:29:13 +0000
Subject: [PATCH 01/24] Use localmesh for d2x and d2y, and interpolate them to
 location

Previously the localmesh was not passed through to the constructor of
d2x and d2y in Coordinates::geometry(), so they used the global 'mesh'.

Also, when d2x/d2y are read from the mesh, they are at CELL_CENTRE, so
we need to interpolate them to the location of the Coordinates.
---
 src/mesh/coordinates.cxx | 8 +++++++-
 1 file changed, 7 insertions(+), 1 deletion(-)

diff --git a/src/mesh/coordinates.cxx b/src/mesh/coordinates.cxx
index 81c0e921e4..8b287638df 100644
--- a/src/mesh/coordinates.cxx
+++ b/src/mesh/coordinates.cxx
@@ -442,13 +442,16 @@ int Coordinates::geometry() {
 
   OPTION(Options::getRoot(), non_uniform, true);
 
-  Field2D d2x, d2y; // d^2 x / d i^2
+  Field2D d2x(localmesh), d2y(localmesh); // d^2 x / d i^2
   // Read correction for non-uniform meshes
   if (localmesh->get(d2x, "d2x")) {
     output_warn.write(
         "\tWARNING: differencing quantity 'd2x' not found. Calculating from dx\n");
     d1_dx = localmesh->indexDDX(1. / dx); // d/di(1/dx)
   } else {
+    // Shift d2x to our location
+    d2x = interp_to(d2x, location);
+
     d1_dx = -d2x / (dx * dx);
   }
 
@@ -457,6 +460,9 @@ int Coordinates::geometry() {
         "\tWARNING: differencing quantity 'd2y' not found. Calculating from dy\n");
     d1_dy = localmesh->indexDDY(1. / dy); // d/di(1/dy)
   } else {
+    // Shift d2y to our location
+    d2y = interp_to(d2y, location);
+
     d1_dy = -d2y / (dy * dy);
   }
 

From 8795acad2b5ec34f53bf815604113348849e502a Mon Sep 17 00:00:00 2001
From: John Omotani <john.omotani@ukaea.uk>
Date: Tue, 20 Nov 2018 14:32:15 +0000
Subject: [PATCH 02/24] Remove uses of global 'mesh' in
 BoutMesh::addBoundaryRegions()

Could cause bugs, for example in unit tests of BoutMesh where
mesh==nullptr.
---
 src/mesh/impls/bout/boutmesh.cxx | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/mesh/impls/bout/boutmesh.cxx b/src/mesh/impls/bout/boutmesh.cxx
index e4f3c487a0..e2516fd771 100644
--- a/src/mesh/impls/bout/boutmesh.cxx
+++ b/src/mesh/impls/bout/boutmesh.cxx
@@ -2208,7 +2208,7 @@ void BoutMesh::addBoundaryRegions() {
   all_boundaries.emplace_back("RGN_UPPER_Y");
   
   // Inner X
-  if(mesh->firstX() && !mesh->periodicX) {
+  if(firstX() && !periodicX) {
     addRegion3D("RGN_INNER_X", Region<Ind3D>(0, xstart-1, ystart, yend, 0, LocalNz-1,
                                              LocalNy, LocalNz, maxregionblocksize));
     addRegion2D("RGN_INNER_X", Region<Ind2D>(0, xstart-1, ystart, yend, 0, 0,
@@ -2225,7 +2225,7 @@ void BoutMesh::addBoundaryRegions() {
   }
 
   // Outer X
-  if(mesh->firstX() && !mesh->periodicX) {
+  if(firstX() && !periodicX) {
     addRegion3D("RGN_OUTER_X", Region<Ind3D>(xend+1, LocalNx-1, ystart, yend, 0, LocalNz-1,
                                              LocalNy, LocalNz, maxregionblocksize));
     addRegion2D("RGN_OUTER_X", Region<Ind2D>(xend+1, LocalNx-1, ystart, yend, 0, 0,

From 21a0fc548b3221e30869c3a3b17e170f82b5ca5e Mon Sep 17 00:00:00 2001
From: John Omotani <john.omotani@ukaea.uk>
Date: Wed, 14 Nov 2018 15:07:23 +0000
Subject: [PATCH 03/24] Allow initialization of Coordinates with
 Mesh::addCoordinates(location)

Instead of creating Coordinates objects when Mesh::getCoordinates is
called, create them during initialization. As a temporary workaround to
maintain backward compatibility, initialize all locations in
BoutMesh::load() so that they are always available.

Allow the user to explicitly request a location to be available, by
calling Mesh::addCoordinates(location) for each location required. This
is optional in v4.2 but will be required from v4.3.

When Coordinates::geometry() is called, assume that staggered location
Coordinates objects need to be recalculated because the user has changed
the CELL_CENTRE version.

As a temporary workaround, replace geometry() with
geometryNoRecalculate() and geometry() now calls geometryNoRecalculate
and then recalculates the staggered versions. This means geometry() can
be reset back to not recalculating in next (which is OK because we will
add the requirement to call addCoordinates() after calling geometry())
without changing the function signature.
---
 include/bout/coordinates.hxx     |  3 +++
 include/bout/mesh.hxx            | 25 +++++++++++--------
 src/field/field3d.cxx            |  6 +++++
 src/mesh/coordinates.cxx         | 41 +++++++++++++++++++++++++++++---
 src/mesh/impls/bout/boutmesh.cxx | 20 ++++++++++++++++
 src/mesh/mesh.cxx                | 39 ++++++++++++++++++++++++------
 6 files changed, 114 insertions(+), 20 deletions(-)

diff --git a/include/bout/coordinates.hxx b/include/bout/coordinates.hxx
index e649cdfcc7..624b15b745 100644
--- a/include/bout/coordinates.hxx
+++ b/include/bout/coordinates.hxx
@@ -143,6 +143,9 @@ public:
   const Field3D Laplace(const Field3D &f, CELL_LOC outloc=CELL_DEFAULT);
   
 private:
+  // temporary work-around method to allow 'geometry' to be called without
+  // trying to re-calculate the non-CELL_CENTRE fields
+  int geometryNoRecalculate();
   int nz; // Size of mesh in Z. This is mesh->ngz-1
   Mesh * localmesh;
   CELL_LOC location;
diff --git a/include/bout/mesh.hxx b/include/bout/mesh.hxx
index fcee0e58e3..2c76177640 100644
--- a/include/bout/mesh.hxx
+++ b/include/bout/mesh.hxx
@@ -434,21 +434,26 @@ class Mesh {
     ASSERT1(location != CELL_DEFAULT);
     ASSERT1(location != CELL_VSHIFT);
 
-    if (coords_map.count(location)) { // True branch most common, returns immediately
-      return coords_map[location].get();
-    } else {
-      // No coordinate system set. Create default
-      // Note that this can't be allocated here due to incomplete type
-      // (circular dependency between Mesh and Coordinates)
-      coords_map.emplace(location, createDefaultCoordinates(location));
-      return coords_map[location].get();
+#if CHECK > 0
+    if (!coords_map.count(location)) {
+      throw BoutException("Error: Coordinates for %s have not been added to "
+          "this Mesh. You should use REQUEST_STAGGER(location) before "
+          "initializing fields staggered to 'CELL_LOC location'.",
+          CELL_LOC_STRING(location).c_str());
     }
+#endif
+    return coords_map.at(location).get();
   }
 
   Coordinates *DEPRECATED(coordinates(const CELL_LOC location = CELL_CENTRE)) {
     return getCoordinates(location);
   }
 
+  /// Add Coordinates object at a certain location.
+  /// If replace_coords is set to true, reset the object in coords_map if it
+  /// already exists, otherwise add a new one
+  void addCoordinates(const CELL_LOC location, bool replace_coords = false);
+
   // First derivatives in index space
   // Implemented in src/mesh/index_derivs.hxx
 
@@ -712,7 +717,7 @@ class Mesh {
   
   GridDataSource *source; ///< Source for grid data
   
-  std::map<CELL_LOC, std::shared_ptr<Coordinates> > coords_map; ///< Coordinate systems at different CELL_LOCs
+  std::map<CELL_LOC, std::unique_ptr<Coordinates> > coords_map; ///< Coordinate systems at different CELL_LOCs
 
   Options *options; ///< Mesh options section
   
@@ -751,7 +756,7 @@ class Mesh {
 
 private:
   /// Allocates default Coordinates objects
-  std::shared_ptr<Coordinates> createDefaultCoordinates(const CELL_LOC location);
+  std::unique_ptr<Coordinates> createDefaultCoordinates(const CELL_LOC location);
 
   //Internal region related information
   std::map<std::string, Region<Ind3D>> regionMap3D;
diff --git a/src/field/field3d.cxx b/src/field/field3d.cxx
index 1143eb9fdc..022fd06aca 100644
--- a/src/field/field3d.cxx
+++ b/src/field/field3d.cxx
@@ -240,6 +240,12 @@ void Field3D::setLocation(CELL_LOC new_location) {
     }
     location = new_location;
 
+#if CHECK > 1
+    // Check Coordinates for location have been added
+    // For CHECK > 0, getCoordinates will throw if location has not been added.
+    fieldmesh->getCoordinates(location);
+#endif
+
     // Invalidate the coordinates pointer
     if (new_location != location)
       fieldCoordinates = nullptr;
diff --git a/src/mesh/coordinates.cxx b/src/mesh/coordinates.cxx
index 8b287638df..11075b37df 100644
--- a/src/mesh/coordinates.cxx
+++ b/src/mesh/coordinates.cxx
@@ -151,7 +151,7 @@ Coordinates::Coordinates(Mesh *mesh)
 
   //////////////////////////////////////////////////////
   /// Calculate Christoffel symbols. Needs communication
-  if (geometry()) {
+  if (geometryNoRecalculate()) {
     throw BoutException("Differential geometry failed\n");
   }
 
@@ -269,7 +269,7 @@ Coordinates::Coordinates(Mesh *mesh, const CELL_LOC loc, const Coordinates* coor
 
   //////////////////////////////////////////////////////
   /// Calculate Christoffel symbols. Needs communication
-  if (geometry()) {
+  if (geometryNoRecalculate()) {
     throw BoutException("Differential geometry failed\n");
   }
 
@@ -304,7 +304,7 @@ void Coordinates::outputVars(Datafile &file) {
   file.add(J, "J", false);
 }
 
-int Coordinates::geometry() {
+int Coordinates::geometryNoRecalculate() {
   TRACE("Coordinates::geometry");
 
   output_progress.write("Calculating differential geometry terms\n");
@@ -469,6 +469,41 @@ int Coordinates::geometry() {
   return 0;
 }
 
+int Coordinates::geometry() {
+
+  geometryNoRecalculate();
+
+  if (location != CELL_CENTRE) {
+    throw BoutException("geometry() called from a location other than "
+      "CELL_CENTRE. This is an error as the other Coordinates are calculated "
+      "from the CELL_CENTRE version, so the changes you have made to this "
+      "object would be overwritten.");
+  }
+  // Coordinates objects at staggered location were calculated from
+  // CELL_CENTRE ones. geometry() has been called on the CELL_CENTRE
+  // Coordinates, so they must have changed; we need to re-calculate the
+  // staggered location Coordinates objects.
+
+  if (localmesh->StaggerGrids) {
+    // Replace Coordinates objects at staggered locations, if there are
+    // enough grid points.
+    // This is a temporary workaround. In v4.3 we will users to call
+    // REQUEST_LOCATION(location) for each location that is needed and change
+    // this so that we don't waste memory on unneeded Coordinates.
+    if (localmesh->LocalNx >= 4) {
+      localmesh->addCoordinates(CELL_XLOW, true);
+    }
+    if (localmesh->LocalNy >= 4) {
+      localmesh->addCoordinates(CELL_YLOW, true);
+    }
+    // Can always add ZLOW Coordinates, since z-interpolation on Field2D is a
+    // null operation
+    localmesh->addCoordinates(CELL_ZLOW, true);
+  }
+
+  return 0;
+}
+
 int Coordinates::calcCovariant() {
   TRACE("Coordinates::calcCovariant");
 
diff --git a/src/mesh/impls/bout/boutmesh.cxx b/src/mesh/impls/bout/boutmesh.cxx
index e2516fd771..62b6a11979 100644
--- a/src/mesh/impls/bout/boutmesh.cxx
+++ b/src/mesh/impls/bout/boutmesh.cxx
@@ -839,6 +839,26 @@ int BoutMesh::load() {
   // Add boundary regions
   addBoundaryRegions();
 
+  // Create CELL_CENTRE Coordinates object
+  addCoordinates(CELL_CENTRE);
+
+  if (StaggerGrids) {
+    // Add Coordinates objects at staggered locations, if there are enough grid
+    // points.
+    // This is a temporary workaround. In v4.3 we will users to call
+    // REQUEST_LOCATION(location) for each location that is needed and remove
+    // this so that we don't waste memory on unneeded Coordinates.
+    if (LocalNx >= 4) {
+      addCoordinates(CELL_XLOW);
+    }
+    if (LocalNy >= 4) {
+      addCoordinates(CELL_YLOW);
+    }
+    // Can always add ZLOW Coordinates, since z-interpolation on Field2D is a
+    // null operation
+    addCoordinates(CELL_ZLOW);
+  }
+
   output_info.write("\tdone\n");
 
   return 0;
diff --git a/src/mesh/mesh.cxx b/src/mesh/mesh.cxx
index 6a857f57a7..023beabccc 100644
--- a/src/mesh/mesh.cxx
+++ b/src/mesh/mesh.cxx
@@ -325,13 +325,38 @@ ParallelTransform& Mesh::getParallelTransform() {
   return *transform;
 }
 
-std::shared_ptr<Coordinates> Mesh::createDefaultCoordinates(const CELL_LOC location) {
-  if (location == CELL_CENTRE || location == CELL_DEFAULT)
-    // Initialize coordinates from input
-    return std::make_shared<Coordinates>(this);
-  else
-    // Interpolate coordinates from CELL_CENTRE version
-    return std::make_shared<Coordinates>(this, location, getCoordinates(CELL_CENTRE));
+void Mesh::addCoordinates(const CELL_LOC location, bool replace_coords) {
+  ASSERT1(location != CELL_DEFAULT);
+
+  if (location == CELL_VSHIFT) {
+    // CELL_VSHIFT puts vector components at CELL_XLOW, CELL_YLOW and
+    // CELL_ZLOW, so require Coordinates at all three.
+    addCoordinates(CELL_XLOW, replace_coords);
+    addCoordinates(CELL_YLOW, replace_coords);
+    addCoordinates(CELL_ZLOW, replace_coords);
+  } else {
+    // No coordinate system set. Create default
+    if (location == CELL_CENTRE) {
+      // Initialize coordinates from input
+      if (!coords_map.count(location)) {
+        // location does not exist in coords_map, so create new entry
+        coords_map.emplace(location, std::unique_ptr<Coordinates>(new Coordinates(this)));
+      } else if (replace_coords) {
+        // location does already exists in coords_map, so reset it
+        coords_map.at(location).reset(new Coordinates(this));
+      }
+    } else {
+      // Interpolate coordinates from CELL_CENTRE version
+      ASSERT1(StaggerGrids); // If StaggerGrids==false, it doesn't make sense to have non-CELL_CENTRE Coordinates
+      if (!coords_map.count(location)) {
+        // location does not exist in coords_map, so create new entry
+        coords_map.emplace(location, std::unique_ptr<Coordinates>(new Coordinates(this, location, getCoordinates(CELL_CENTRE))));
+      } else if (replace_coords) {
+        // location does already exists in coords_map, so reset it
+        coords_map.at(location).reset(new Coordinates(this, location, getCoordinates(CELL_CENTRE)));
+      }
+    }
+  }
 }
 
 

From f528e6aba79be79f13f09cc213c937fa2b4e2e05 Mon Sep 17 00:00:00 2001
From: John Omotani <john.omotani@ukaea.uk>
Date: Sun, 18 Nov 2018 20:50:22 +0000
Subject: [PATCH 04/24] Add nullptr entries to coords_map in FakeMesh of unit
 tests

Allows Coordinates pointers to be retrieved where necessary in the unit
tests.
---
 tests/unit/test_extras.hxx | 5 +++++
 1 file changed, 5 insertions(+)

diff --git a/tests/unit/test_extras.hxx b/tests/unit/test_extras.hxx
index 6e6c9dc2ad..b143edef1f 100644
--- a/tests/unit/test_extras.hxx
+++ b/tests/unit/test_extras.hxx
@@ -73,6 +73,11 @@ public:
     StaggerGrids = false;
     IncIntShear = false;
     maxregionblocksize = MAXREGIONBLOCKSIZE;
+
+    coords_map.emplace(CELL_CENTRE, std::unique_ptr<Coordinates>(nullptr));
+    coords_map.emplace(CELL_XLOW, std::unique_ptr<Coordinates>(nullptr));
+    coords_map.emplace(CELL_YLOW, std::unique_ptr<Coordinates>(nullptr));
+    coords_map.emplace(CELL_ZLOW, std::unique_ptr<Coordinates>(nullptr));
   }
 
   comm_handle send(FieldGroup &UNUSED(g)) { return nullptr; };

From 85e054d78b2476e31b0cb7d6fee6c8ad371f07d5 Mon Sep 17 00:00:00 2001
From: John Omotani <john.omotani@ukaea.uk>
Date: Wed, 14 Nov 2018 15:07:23 +0000
Subject: [PATCH 05/24] Allow initialization of Coordinates with
 Mesh::addCoordinates(location)

Instead of creating Coordinates objects when Mesh::getCoordinates is
called, create them during initialization. As a temporary workaround to
maintain backward compatibility, initialize all locations in
BoutMesh::load() so that they are always available.

Allow the user to explicitly request a location to be available, by
calling Mesh::addCoordinates(location) for each location required. This
is optional in v4.2 but will be required from v4.3.

When Coordinates::geometry() is called, assume that staggered location
Coordinates objects need to be recalculated because the user has changed
the CELL_CENTRE version.

As a temporary workaround, replace geometry() with
geometryNoRecalculate() and geometry() now calls geometryNoRecalculate
and then recalculates the staggered versions. This means geometry() can
be reset back to not recalculating in next (which is OK because we will
add the requirement to call addCoordinates() after calling geometry())
without changing the function signature.
---
 include/bout/mesh.hxx | 1 +
 1 file changed, 1 insertion(+)

diff --git a/include/bout/mesh.hxx b/include/bout/mesh.hxx
index 2c76177640..d474bdf95d 100644
--- a/include/bout/mesh.hxx
+++ b/include/bout/mesh.hxx
@@ -60,6 +60,7 @@ class Mesh;
 #include "sys/range.hxx" // RangeIterator
 
 #include <bout/griddata.hxx>
+#include <bout/macro_for_each.hxx>
 
 #include "coordinates.hxx"    // Coordinates class
 

From b93a9560966e6fa4a8606880968a9adf2df21cc2 Mon Sep 17 00:00:00 2001
From: John Omotani <john.omotani@ukaea.uk>
Date: Fri, 16 Nov 2018 09:58:31 +0000
Subject: [PATCH 06/24] Require use of Mesh::addCoordinates(location)

Remove initialization of Coordinates for all staggered locations in
BoutMesh::load(). User is now required to initialize the Coordinates for
any needed locations using Mesh::addCoordinates(location).

Removes recalculation of staggered Coordinates in
Coordinates::geometry(). The user should instead call
addCoordinates(location) after the CELL_CENTRE Coordinates object is
updated and its geometry() method is called. geometry() checks that
either it has been called at CELL_CENTRE and CELL_CENTRE is the only
entry in coords_map, or that the location of the current Coordinates
object that called geometry() has not yet been added to coords_map; this
check can be overridden by setting the option
mesh:allow_geometry_without_recalculate_staggered=true so that the
staggered Coordinates objects can be explicitly updated if necessary.

If Mesh::addCoordinates(location) is called without explicitly asking
for the Coordinates at location to be replaced (using the replace_coords
argument) then throw an exception if location is already in coords_map.
---
 include/bout/coordinates.hxx     |  4 +-
 include/bout/mesh.hxx            | 13 ++++++
 src/mesh/coordinates.cxx         | 71 ++++++++++++++++----------------
 src/mesh/impls/bout/boutmesh.cxx | 17 --------
 src/mesh/mesh.cxx                | 16 +++++--
 5 files changed, 62 insertions(+), 59 deletions(-)

diff --git a/include/bout/coordinates.hxx b/include/bout/coordinates.hxx
index 624b15b745..33179db8e9 100644
--- a/include/bout/coordinates.hxx
+++ b/include/bout/coordinates.hxx
@@ -143,12 +143,10 @@ public:
   const Field3D Laplace(const Field3D &f, CELL_LOC outloc=CELL_DEFAULT);
   
 private:
-  // temporary work-around method to allow 'geometry' to be called without
-  // trying to re-calculate the non-CELL_CENTRE fields
-  int geometryNoRecalculate();
   int nz; // Size of mesh in Z. This is mesh->ngz-1
   Mesh * localmesh;
   CELL_LOC location;
+  bool allow_geometry_without_recalculate_staggered;
 };
 
 /*
diff --git a/include/bout/mesh.hxx b/include/bout/mesh.hxx
index d474bdf95d..81b103bed5 100644
--- a/include/bout/mesh.hxx
+++ b/include/bout/mesh.hxx
@@ -455,6 +455,19 @@ class Mesh {
   /// already exists, otherwise add a new one
   void addCoordinates(const CELL_LOC location, bool replace_coords = false);
 
+  /// Check if Coordinates object at location has been added
+  bool hasCoordinates(const CELL_LOC location) {
+    return coords_map.count(location);
+  }
+
+  /// Count how many Coordinates objects have been added
+  int countCoordinates() {
+    return coords_map.size();
+  }
+
+  /// switch to pass to Coordinates objects
+  bool allow_geometry_without_recalculate_staggered;
+
   // First derivatives in index space
   // Implemented in src/mesh/index_derivs.hxx
 
diff --git a/src/mesh/coordinates.cxx b/src/mesh/coordinates.cxx
index 11075b37df..9b8d1be823 100644
--- a/src/mesh/coordinates.cxx
+++ b/src/mesh/coordinates.cxx
@@ -26,7 +26,9 @@ Coordinates::Coordinates(Mesh *mesh)
       G1_23(mesh), G2_11(mesh), G2_22(mesh), G2_33(mesh), G2_12(mesh), G2_13(mesh),
       G2_23(mesh), G3_11(mesh), G3_22(mesh), G3_33(mesh), G3_12(mesh), G3_13(mesh),
       G3_23(mesh), G1(mesh), G2(mesh), G3(mesh), ShiftTorsion(mesh),
-      IntShiftTorsion(mesh), localmesh(mesh), location(CELL_CENTRE) {
+      IntShiftTorsion(mesh), localmesh(mesh), location(CELL_CENTRE),
+      allow_geometry_without_recalculate_staggered(
+          mesh->allow_geometry_without_recalculate_staggered) {
 
   if (mesh->get(dx, "dx")) {
     output_warn.write("\tWARNING: differencing quantity 'dx' not found. Set to 1.0\n");
@@ -151,7 +153,7 @@ Coordinates::Coordinates(Mesh *mesh)
 
   //////////////////////////////////////////////////////
   /// Calculate Christoffel symbols. Needs communication
-  if (geometryNoRecalculate()) {
+  if (geometry()) {
     throw BoutException("Differential geometry failed\n");
   }
 
@@ -269,7 +271,7 @@ Coordinates::Coordinates(Mesh *mesh, const CELL_LOC loc, const Coordinates* coor
 
   //////////////////////////////////////////////////////
   /// Calculate Christoffel symbols. Needs communication
-  if (geometryNoRecalculate()) {
+  if (geometry()) {
     throw BoutException("Differential geometry failed\n");
   }
 
@@ -304,7 +306,7 @@ void Coordinates::outputVars(Datafile &file) {
   file.add(J, "J", false);
 }
 
-int Coordinates::geometryNoRecalculate() {
+int Coordinates::geometry() {
   TRACE("Coordinates::geometry");
 
   output_progress.write("Calculating differential geometry terms\n");
@@ -466,39 +468,36 @@ int Coordinates::geometryNoRecalculate() {
     d1_dy = -d2y / (dy * dy);
   }
 
-  return 0;
-}
-
-int Coordinates::geometry() {
-
-  geometryNoRecalculate();
-
-  if (location != CELL_CENTRE) {
-    throw BoutException("geometry() called from a location other than "
-      "CELL_CENTRE. This is an error as the other Coordinates are calculated "
-      "from the CELL_CENTRE version, so the changes you have made to this "
-      "object would be overwritten.");
-  }
-  // Coordinates objects at staggered location were calculated from
-  // CELL_CENTRE ones. geometry() has been called on the CELL_CENTRE
-  // Coordinates, so they must have changed; we need to re-calculate the
-  // staggered location Coordinates objects.
-
-  if (localmesh->StaggerGrids) {
-    // Replace Coordinates objects at staggered locations, if there are
-    // enough grid points.
-    // This is a temporary workaround. In v4.3 we will users to call
-    // REQUEST_LOCATION(location) for each location that is needed and change
-    // this so that we don't waste memory on unneeded Coordinates.
-    if (localmesh->LocalNx >= 4) {
-      localmesh->addCoordinates(CELL_XLOW, true);
-    }
-    if (localmesh->LocalNy >= 4) {
-      localmesh->addCoordinates(CELL_YLOW, true);
+  if (!allow_geometry_without_recalculate_staggered) {
+    // Only CELL_CENTRE Coordinates should ever be changed (and therefore need
+    // geometry() calling). If CELL_CENTRE Coordinates are changed they should
+    // be changed before other Coordinates are added to localmesh->coords_map,
+    // since Coordinates at other locations are calculated from the CELL_CENTRE
+    // ones. The allow_geometry_without_recalculate_staggered option overrides
+    // this check, in case for some reason the user does need to call
+    // geometry() in other situations.
+    if ( (localmesh->hasCoordinates(location))                                            // this location already added
+         && !(localmesh->countCoordinates()==1 && localmesh->hasCoordinates(CELL_CENTRE)) // OK to be added already only if CELL_CENTRE and no other locations added
+       ) {
+      if (location == CELL_CENTRE) {
+        throw BoutException("Coordinates::geometry() called at CELL_CENTRE, but "
+          "other locations have already been added to the Mesh. These would "
+          "need recalculating. If possible, call Mesh::addCoordinates() after "
+          "this call to geometry(). To recalculate the other Coordinates objects "
+          "from the CELL_CENTRE Coordinates, pass true for the "
+          "recalculate_staggered argument, e.g.  coords->geometry(true). If you "
+          "need to recalculate multiple Coordinates objects explicitly, set "
+          "mesh:allow_geometry_without_recalculate_staggered=true in the input "
+          "file to disable this check.");
+      } else {
+        throw BoutException("Coordinates::geometry() called at location %s, but "
+          "this location has already been initialized and added to the Mesh. "
+          "If you need to recalculate multiple Coordinates objects explicitly, "
+          "set mesh:allow_geometry_without_recalculate_staggered=true in the "
+          "input file to disable this check.", CELL_LOC_STRING(location).c_str());
+      }
     }
-    // Can always add ZLOW Coordinates, since z-interpolation on Field2D is a
-    // null operation
-    localmesh->addCoordinates(CELL_ZLOW, true);
+
   }
 
   return 0;
diff --git a/src/mesh/impls/bout/boutmesh.cxx b/src/mesh/impls/bout/boutmesh.cxx
index 62b6a11979..b29239a600 100644
--- a/src/mesh/impls/bout/boutmesh.cxx
+++ b/src/mesh/impls/bout/boutmesh.cxx
@@ -842,23 +842,6 @@ int BoutMesh::load() {
   // Create CELL_CENTRE Coordinates object
   addCoordinates(CELL_CENTRE);
 
-  if (StaggerGrids) {
-    // Add Coordinates objects at staggered locations, if there are enough grid
-    // points.
-    // This is a temporary workaround. In v4.3 we will users to call
-    // REQUEST_LOCATION(location) for each location that is needed and remove
-    // this so that we don't waste memory on unneeded Coordinates.
-    if (LocalNx >= 4) {
-      addCoordinates(CELL_XLOW);
-    }
-    if (LocalNy >= 4) {
-      addCoordinates(CELL_YLOW);
-    }
-    // Can always add ZLOW Coordinates, since z-interpolation on Field2D is a
-    // null operation
-    addCoordinates(CELL_ZLOW);
-  }
-
   output_info.write("\tdone\n");
 
   return 0;
diff --git a/src/mesh/mesh.cxx b/src/mesh/mesh.cxx
index 023beabccc..604f33a3b8 100644
--- a/src/mesh/mesh.cxx
+++ b/src/mesh/mesh.cxx
@@ -31,6 +31,7 @@ Mesh::Mesh(GridDataSource *s, Options* opt) : source(s), options(opt) {
   /// Get mesh options
   OPTION(options, StaggerGrids,   false); // Stagger grids
   OPTION(options, maxregionblocksize, MAXREGIONBLOCKSIZE);
+  OPTION(options, allow_geometry_without_recalculate_staggered, false);
   // Initialise derivatives
   derivs_init(options);  // in index_derivs.cxx for now
 }
@@ -342,7 +343,7 @@ void Mesh::addCoordinates(const CELL_LOC location, bool replace_coords) {
         // location does not exist in coords_map, so create new entry
         coords_map.emplace(location, std::unique_ptr<Coordinates>(new Coordinates(this)));
       } else if (replace_coords) {
-        // location does already exists in coords_map, so reset it
+        // location does already exist in coords_map, so reset it
         coords_map.at(location).reset(new Coordinates(this));
       }
     } else {
@@ -352,8 +353,17 @@ void Mesh::addCoordinates(const CELL_LOC location, bool replace_coords) {
         // location does not exist in coords_map, so create new entry
         coords_map.emplace(location, std::unique_ptr<Coordinates>(new Coordinates(this, location, getCoordinates(CELL_CENTRE))));
       } else if (replace_coords) {
-        // location does already exists in coords_map, so reset it
-        coords_map.at(location).reset(new Coordinates(this, location, getCoordinates(CELL_CENTRE)));
+        // location does already exist in coords_map, so reset it
+
+        // first erase the existing entry to avoid throwing an exception from
+        // Coordinates::geometry().
+        // The check that would throw the exception is not needed because
+        // replacement of the Coordinates object has been explicitly requested.
+        coords_map.erase(location);
+        coords_map.emplace(location, std::unique_ptr<Coordinates>(new Coordinates(this, location, getCoordinates(CELL_CENTRE))));
+      } else {
+        throw BoutException("Coordinates at %s already added to Mesh",
+            CELL_LOC_STRING(location).c_str());
       }
     }
   }

From d5675d2aea35e89009820b4906757da405799d4a Mon Sep 17 00:00:00 2001
From: John Omotani <john.omotani@ukaea.uk>
Date: Wed, 14 Nov 2018 15:19:00 +0000
Subject: [PATCH 07/24] Add Mesh::addCoordinates(location) calls in tests

tests/integrated/test-drift-instability, tests/MMS/wave-1d and
tests/MMS/wave-1d-y use staggered grids and so need to call
Mesh::addCoordinates(location).
---
 tests/MMS/wave-1d-y/wave.cxx                       | 1 +
 tests/MMS/wave-1d/wave.cxx                         | 1 +
 tests/integrated/test-drift-instability/2fluid.cxx | 1 +
 3 files changed, 3 insertions(+)

diff --git a/tests/MMS/wave-1d-y/wave.cxx b/tests/MMS/wave-1d-y/wave.cxx
index db4c931665..1b359be200 100644
--- a/tests/MMS/wave-1d-y/wave.cxx
+++ b/tests/MMS/wave-1d-y/wave.cxx
@@ -9,6 +9,7 @@ class Wave1D : public PhysicsModel {
 protected:
   int init(bool restarting) {
 
+    mesh->addCoordinates(CELL_YLOW);
     g.setLocation(CELL_YLOW); // g staggered 
     
     // Tell BOUT++ to solve f and g
diff --git a/tests/MMS/wave-1d/wave.cxx b/tests/MMS/wave-1d/wave.cxx
index d2543a98b1..244adc18ac 100644
--- a/tests/MMS/wave-1d/wave.cxx
+++ b/tests/MMS/wave-1d/wave.cxx
@@ -83,6 +83,7 @@ class Wave1D : public PhysicsModel {
     coord->g_23 = 0.0;
     coord->geometry();
 
+    mesh->addCoordinates(CELL_XLOW);
     g.setLocation(CELL_XLOW); // g staggered to the left of f
     
     //Dirichlet everywhere except inner x-boundary Neumann
diff --git a/tests/integrated/test-drift-instability/2fluid.cxx b/tests/integrated/test-drift-instability/2fluid.cxx
index 3e55a2c6fb..ff805f3676 100644
--- a/tests/integrated/test-drift-instability/2fluid.cxx
+++ b/tests/integrated/test-drift-instability/2fluid.cxx
@@ -297,6 +297,7 @@ int physics_init(bool restarting) {
   dump.add(wci,   "wci", 0);
 
   if (mesh->StaggerGrids) {
+    mesh->addCoordinates(CELL_YLOW);
     maybe_ylow = CELL_YLOW;
   } else {
     maybe_ylow = CELL_CENTRE;

From 73480c4c66bb4809bc4c980d661ea47ae8c8891a Mon Sep 17 00:00:00 2001
From: John Omotani <john.omotani@ukaea.uk>
Date: Wed, 14 Nov 2018 16:01:43 +0000
Subject: [PATCH 08/24] Update examples/staggered_grid

Move StaggerGrids option into [mesh] section of input files.

Call mesh->addCoordinates(CELL_YLOW).

Only use CELL_YLOW if mesh->StaggerGrids=true to avoid exceptions.
---
 examples/staggered_grid/data/BOUT.inp      |  2 ++
 examples/staggered_grid/test/BOUT.inp      |  2 +-
 examples/staggered_grid/test_staggered.cxx | 11 +++++++++--
 3 files changed, 12 insertions(+), 3 deletions(-)

diff --git a/examples/staggered_grid/data/BOUT.inp b/examples/staggered_grid/data/BOUT.inp
index 82261a569f..2d1e29bed6 100644
--- a/examples/staggered_grid/data/BOUT.inp
+++ b/examples/staggered_grid/data/BOUT.inp
@@ -6,6 +6,8 @@ MZ = 1
 
 grid = "test-staggered.nc"
 
+[mesh]
+
 StaggerGrids = true
 
 [mesh:ddy]
diff --git a/examples/staggered_grid/test/BOUT.inp b/examples/staggered_grid/test/BOUT.inp
index 2e6ec630c2..ce28b7edac 100644
--- a/examples/staggered_grid/test/BOUT.inp
+++ b/examples/staggered_grid/test/BOUT.inp
@@ -4,9 +4,9 @@ timestep = 0.02
 
 MZ = 1
 
+[mesh]
 StaggerGrids = true
 
-[mesh]
 nx = 5
 ny = 16
 
diff --git a/examples/staggered_grid/test_staggered.cxx b/examples/staggered_grid/test_staggered.cxx
index d49a854deb..7eacdb2282 100644
--- a/examples/staggered_grid/test_staggered.cxx
+++ b/examples/staggered_grid/test_staggered.cxx
@@ -8,10 +8,17 @@
 #include <derivs.hxx>
 
 Field3D n, v;
+CELL_LOC maybe_ylow = CELL_CENTRE;
 
 int physics_init(bool restart) {
   
-  v.setLocation(CELL_YLOW); // Staggered relative to n
+  if (mesh->StaggerGrids) {
+    maybe_ylow = CELL_YLOW;
+
+    mesh->addCoordinates(CELL_YLOW);
+
+    v.setLocation(CELL_YLOW); // Staggered relative to n
+  }
   
   SOLVE_FOR(n, v);
 
@@ -24,7 +31,7 @@ int physics_run(BoutReal time) {
   //ddt(n) = -Div_par_flux(v, n, CELL_CENTRE);
   ddt(n) = -n*Grad_par(v, CELL_CENTRE) - Vpar_Grad_par(v, n, CELL_CENTRE);
   
-  ddt(v) = -Grad_par(n, CELL_YLOW);
+  ddt(v) = -Grad_par(n, maybe_ylow);
  
   // Have to manually apply the lower Y boundary region, using a width of 3
   for( RangeIterator rlow = mesh->iterateBndryLowerY(); !rlow.isDone(); rlow++)

From cc0c57d1a57855a528b393994516781def5fec22 Mon Sep 17 00:00:00 2001
From: John Omotani <john.omotani@ukaea.uk>
Date: Wed, 14 Nov 2018 16:07:26 +0000
Subject: [PATCH 09/24] Update manual with addCoordinates

---
 manual/sphinx/user_docs/staggered_grids.rst | 43 ++++++++++++++++++---
 1 file changed, 38 insertions(+), 5 deletions(-)

diff --git a/manual/sphinx/user_docs/staggered_grids.rst b/manual/sphinx/user_docs/staggered_grids.rst
index 275d082766..70d60a2e64 100644
--- a/manual/sphinx/user_docs/staggered_grids.rst
+++ b/manual/sphinx/user_docs/staggered_grids.rst
@@ -16,8 +16,9 @@ staggered grids, set::
 
     StaggerGrids = true
 
-in the top section of the ``BOUT.inp`` file. The **test-staggered**
-example illustrates how to use staggered grids in BOUT++.
+in the top section of the ``BOUT.inp`` file and enable the locations you will
+use with a call to mesh->addCoordinates(location) (see below).  The
+**test-staggered** example illustrates how to use staggered grids in BOUT++.
 
 There are four possible locations in a grid cell where a quantity can be
 defined in BOUT++: centre, lower X, lower Y, and lower Z. These are
@@ -29,9 +30,10 @@ illustrated in :numref:`staggergrids-location`.
 
    The four possible cell locations for defining quantities
 
-To specify the location of a variable, use the method
-`Field3D::setLocation` with one of the `CELL_LOC` locations
-`CELL_CENTRE`, `CELL_XLOW`, `CELL_YLOW`, or `CELL_ZLOW`.
+The possible locations are specified with the `CELL_LOC` type, which has the
+possible values `CELL_CENTRE`, `CELL_XLOW`, `CELL_YLOW`, or `CELL_ZLOW`.
+`CELL_CENTRE` is enabled by default, but a call to
+mesh->addCoordinates(location) is required to enable the others.
 
 The key lines in the **staggered_grid** example which specify the
 locations of the evolving variables are::
@@ -39,7 +41,11 @@ locations of the evolving variables are::
     Field3D n, v;
 
     int init(bool restart) {
+
+      mesh->addCoordinates(CELL_YLOW);
+
       v.setLocation(CELL_YLOW); // Staggered relative to n
+
       SOLVE_FOR(n, v);
       ...
 
@@ -62,6 +68,33 @@ in Y, whilst the density :math:`n` remains cell centred.
           `CELL_CENTRE` if staggered grids are off, regardless of what
           you pass it.
 
+.. note:: For advanced users:
+          If you change members of the Coordinates object manually, you should
+          change the CELL_CENTRE Coordinates and only call addCoordinates() for
+          other locations after you call Coordinates::geometry() on the
+          CELL_CENTRE Coordinates. Then the Coordinates at staggered locations
+          will be interpolated from the correct, final CELL_CENTRE version.
+
+          The example uses the global Mesh object 'mesh'. If you are using any
+          other Mesh objects, you need to initialize the Coordinates objects in
+          their coords_map members by calling the addCoordinates(CELL_LOC
+          location) method for each location you will use.
+
+          An exception will be thrown if you call Coordinates::geometry() from
+          any Coordinates object at a staggered location, since these are
+          expected to be consistent with (and calculated from) the CELL_CENTRE
+          Coordinates. If you need to change them, you can set the option
+          mesh:allow_geometry_without_recalculate_staggered=true to disable
+          this check; you must then ensure that all the Coordinates objects in
+          Mesh::coords_map are consistent with each other. Setting this option
+          also allows changes to be made and geometry() to be called on the
+          CELL_CENTRE Coordinates after other locations have been added to
+          coords_map; in this case you will need to update the other locations
+          explicitly by calling Mesh::addCoordinates(location, true) - the
+          optional second argument causes addCoordinates to overwrite any
+          existing Coordinates object at location and replace it with one
+          calculated from the current CELL_CENTRE Coordinates.
+
 
 Arithmetic operations can only be performed between variables with the same
 location. When performing a calculation at one location, to include a variable

From 9f6f9c4865ae160ddf7970a4e9da5217df9a9858 Mon Sep 17 00:00:00 2001
From: John Omotani <john.omotani@ukaea.uk>
Date: Wed, 14 Nov 2018 17:13:11 +0000
Subject: [PATCH 10/24] Update boutcore with getCoordinates and addCoordinates

---
 tools/pylib/_boutcore_build/boutcore.pyx.in | 26 ++++++++++++++++-----
 tools/pylib/_boutcore_build/boutcpp.pxd.in  |  3 ++-
 2 files changed, 22 insertions(+), 7 deletions(-)

diff --git a/tools/pylib/_boutcore_build/boutcore.pyx.in b/tools/pylib/_boutcore_build/boutcore.pyx.in
index 6e20f61ae6..21deb239c1 100755
--- a/tools/pylib/_boutcore_build/boutcore.pyx.in
+++ b/tools/pylib/_boutcore_build/boutcore.pyx.in
@@ -627,7 +627,6 @@ cdef class Mesh:
     cdef c.bool isGlobal
     cdef double isNormalised
     cdef FieldFactory factory
-    cdef Coordinates _coords
     #factory=FieldFactory()
     def __init__(self, create=True, section=None, options=None):
         """
@@ -655,7 +654,6 @@ cdef class Mesh:
         self.isGlobal=False
         self.isNormalised=-1
         self.factory = FieldFactory()
-        self._coords = None
         if create:
             if options:
                 opt = (<Options?>options).cobj
@@ -740,14 +738,30 @@ cdef class Mesh:
         del fg
         return self
 
+    def getCoordinates(self, location = "CENTRE"):
+        """
+        Get a Coordinates object from this mesh
+        """
+        # resolve the location string to a CELL_LOC first, to check it is valid
+        loc_ = benum.resolve_cell_loc(location)
+        return coordsFromObj(self.cobj.getCoordinates(loc_))
+
     @property
     def coordinates(self):
         """
-        Get the Coordinates object of this mesh
+        Deprecated version of getCoordinates
+        """
+        print("Warning Mesh.coordinates is deprecated, and does not handle "
+              "staggered grid locations. Use Mesh.getCoordinates(location) "
+              "method instead")
+        return self.getCoordinates()
+
+    def addCoordinates(self, location):
+        """
+        Initialize Coordinates object at location
         """
-        if self._coords is None:
-            self._coords = coordsFromObj(self.cobj.coordinates())
-        return self._coords
+        loc_ = benum.resolve_cell_loc(location)
+        self.cobj.addCoordinates(loc_)
 
 cdef Coordinates coordsFromObj(c.Coordinates * obj):
     coords = Coordinates()
diff --git a/tools/pylib/_boutcore_build/boutcpp.pxd.in b/tools/pylib/_boutcore_build/boutcpp.pxd.in
index 49edc12332..6eb8836d64 100644
--- a/tools/pylib/_boutcore_build/boutcpp.pxd.in
+++ b/tools/pylib/_boutcore_build/boutcpp.pxd.in
@@ -60,7 +60,8 @@ cdef extern from "bout/mesh.hxx":
         int ystart
         int LocalNx
         int LocalNy
-        Coordinates * coordinates()
+        Coordinates * getCoordinates(benum.CELL_LOC location) except +
+        void addCoordinates(benum.CELL_LOC location) except +
 
 cdef extern from "bout/coordinates.hxx":
     cppclass Coordinates:

From 5353de5d14c47d965ef5d036b69ab235ca963fac Mon Sep 17 00:00:00 2001
From: John Omotani <john.omotani@ukaea.uk>
Date: Wed, 14 Nov 2018 17:13:28 +0000
Subject: [PATCH 11/24] Fix MMS tests

Need to call Mesh.addCoordinates(location) to initialize staggered
grids in tests/MMS/derivatives3 and tests/MMS/upwinding3.
---
 tests/MMS/derivatives3/runtest | 2 ++
 tests/MMS/upwinding3/runtest   | 2 ++
 2 files changed, 4 insertions(+)

diff --git a/tests/MMS/derivatives3/runtest b/tests/MMS/derivatives3/runtest
index c3f64624cb..e19e4a1c1f 100755
--- a/tests/MMS/derivatives3/runtest
+++ b/tests/MMS/derivatives3/runtest
@@ -33,6 +33,8 @@ def runtests(functions,derivatives,directions,stag,msg):
                                ,"2*pi/(%d)"%(nz),force=True)
             dirnfac=direction+"*"+fac
             mesh=boutcore.Mesh(section="mesh"+direction)
+            for loc in locations[1:]:
+                mesh.addCoordinates(loc)
             f=boutcore.create3D(infunc.replace("%s",dirnfac),mesh
                                 ,outloc=inloc)
             sim=diff_func(f,method=diff,outloc=outloc)
diff --git a/tests/MMS/upwinding3/runtest b/tests/MMS/upwinding3/runtest
index c0fd2129a3..eae8667582 100755
--- a/tests/MMS/upwinding3/runtest
+++ b/tests/MMS/upwinding3/runtest
@@ -33,6 +33,8 @@ def runtests(functions,derivatives,directions,stag,msg):
                                ,"2*pi/(%d)"%(nz),force=True)
             dirnfac=direction+"*"+fac
             mesh=boutcore.Mesh(section="mesh"+direction)
+            for loc in locations[1:]:
+                mesh.addCoordinates(loc)
             f=boutcore.create3D(ffunc.replace("%s",dirnfac),mesh
                                 ,outloc=floc)
             v=boutcore.create3D(vfunc.replace("%s",dirnfac),mesh

From b9ae93d9d44a83c4561035b4c416082f670dccb6 Mon Sep 17 00:00:00 2001
From: John Omotani <john.omotani@ukaea.uk>
Date: Tue, 20 Nov 2018 14:42:50 +0000
Subject: [PATCH 12/24] Test addCoordinates in BoutMeshTest

Unit tests for addCoordinates() methods.
---
 tests/unit/mesh/test_boutmesh.cxx | 77 ++++++++++++++++++++++++++-----
 1 file changed, 65 insertions(+), 12 deletions(-)

diff --git a/tests/unit/mesh/test_boutmesh.cxx b/tests/unit/mesh/test_boutmesh.cxx
index 3510bffacf..5d7e2d6dfd 100644
--- a/tests/unit/mesh/test_boutmesh.cxx
+++ b/tests/unit/mesh/test_boutmesh.cxx
@@ -4,43 +4,96 @@
 #include "bout/mesh.hxx"
 #include "output.hxx"
 #include "unused.hxx"
+#include <unordered_map>
 
 class FakeGridDataSource : public GridDataSource {
 public:
   FakeGridDataSource(){};
   ~FakeGridDataSource(){};
   bool hasVar(const string &UNUSED(name)) { return false; };
-  bool get(Mesh *UNUSED(m), int &UNUSED(ival), const string &UNUSED(name)) {
-    return true;
+  bool get(Mesh *UNUSED(m), int &ival, const string &name) {
+    if (intvars.count(name)>0) {
+      ival = intvars.at(name);
+      return true;
+    }
+    return false;
   };
   bool get(Mesh *UNUSED(m), BoutReal &UNUSED(rval), const string &UNUSED(name)) {
-    return true;
+    return false;
   }
   bool get(Mesh *UNUSED(m), Field2D &UNUSED(var), const string &UNUSED(name),
            BoutReal UNUSED(def) = 0.0) {
-    return true;
+    return false;
   }
   bool get(Mesh *UNUSED(m), Field3D &UNUSED(var), const string &UNUSED(name),
            BoutReal UNUSED(def) = 0.0) {
-    return true;
+    return false;
   }
   bool get(Mesh *UNUSED(m), vector<int> &UNUSED(var), const string &UNUSED(name),
            int UNUSED(len), int UNUSED(offset) = 0,
            Direction UNUSED(dir) = GridDataSource::X) {
-    return true;
+    return false;
   }
   bool get(Mesh *UNUSED(m), vector<BoutReal> &UNUSED(var), const string &UNUSED(name),
            int UNUSED(len), int UNUSED(offset) = 0,
            Direction UNUSED(dir) = GridDataSource::X) {
-    return true;
+    return false;
   }
+
+  std::unordered_map<std::string, int> intvars { {"nx", 6}, {"ny", 7 }, {"nz", 5}};
+};
+
+/// Test fixture with a BoutMesh
+class BoutMeshTest : public ::testing::Test {
+protected:
+  static void SetUpTestCase() {
+    output_info.disable();
+    output_warn.disable();
+    output_progress.disable();
+  }
+
+  static void TearDownTestCase() {
+    output_info.enable();
+    output_warn.enable();
+    output_progress.enable();
+  }
+
+public:
+  BoutMeshTest() : source(), localmesh(Mesh::create(&source)) {
+    localmesh->StaggerGrids = true;
+    output_info.disable();
+    localmesh->load();
+  }
+  FakeGridDataSource source;
+  Mesh* localmesh;
 };
 
-TEST(BoutMeshTest, NullOptionsCheck) {
+TEST_F(BoutMeshTest, NullOptionsCheck) {
   // Temporarily turn off outputs to make test quiet
-  output_info.disable();
-  output_warn.disable();
   EXPECT_NO_THROW(BoutMesh mesh(new FakeGridDataSource, nullptr));
-  output_info.enable();
-  output_warn.enable();
+}
+
+TEST_F(BoutMeshTest, AddCoordinatesToMeshCENTRE) {
+  EXPECT_NO_THROW(localmesh->addCoordinates(CELL_CENTRE));
+  EXPECT_NO_THROW(localmesh->getCoordinates(CELL_CENTRE)->geometry());
+  EXPECT_NO_THROW(localmesh->addCoordinates(CELL_YLOW));
+  EXPECT_THROW(localmesh->getCoordinates(CELL_CENTRE)->geometry(), BoutException);
+}
+
+TEST_F(BoutMeshTest, AddCoordinatesToMeshXLOW) {
+  EXPECT_NO_THROW(localmesh->addCoordinates(CELL_XLOW));
+  EXPECT_THROW(localmesh->addCoordinates(CELL_XLOW), BoutException);
+  EXPECT_THROW(localmesh->getCoordinates(CELL_XLOW)->geometry(), BoutException);
+}
+
+TEST_F(BoutMeshTest, AddCoordinatesToMeshYLOW) {
+  EXPECT_NO_THROW(localmesh->addCoordinates(CELL_YLOW));
+  EXPECT_THROW(localmesh->addCoordinates(CELL_YLOW), BoutException);
+  EXPECT_THROW(localmesh->getCoordinates(CELL_YLOW)->geometry(), BoutException);
+}
+
+TEST_F(BoutMeshTest, AddCoordinatesToMeshZLOW) {
+  EXPECT_NO_THROW(localmesh->addCoordinates(CELL_ZLOW));
+  EXPECT_THROW(localmesh->addCoordinates(CELL_ZLOW), BoutException);
+  EXPECT_THROW(localmesh->getCoordinates(CELL_ZLOW)->geometry(), BoutException);
 }

From c2a6fcc33156725cb6b91a32a905ffb3448940cc Mon Sep 17 00:00:00 2001
From: John Omotani <john.omotani@ukaea.uk>
Date: Wed, 21 Nov 2018 11:58:07 +0000
Subject: [PATCH 13/24] Call addCoordinates in
 tests/integrated/test-boutcore/collect-staggered

Call to addCoordinates('YLOW') is now needed to allow 'YLOW' location to
be used.
---
 tests/integrated/test-boutcore/collect-staggered/runtest | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/tests/integrated/test-boutcore/collect-staggered/runtest b/tests/integrated/test-boutcore/collect-staggered/runtest
index 6767618531..28b5e20cf9 100755
--- a/tests/integrated/test-boutcore/collect-staggered/runtest
+++ b/tests/integrated/test-boutcore/collect-staggered/runtest
@@ -6,6 +6,8 @@ import boutcore as bc
 
 bc.init("-q -q -q")
 
+bc.Mesh().getGlobal().addCoordinates('YLOW')
+
 fail=0
 f=bc.create3D("sin(y)",outloc='YLOW')
 

From 1495a35b0a991eeaff13272c55f6d243f4de33d0 Mon Sep 17 00:00:00 2001
From: John Omotani <john.omotani@ukaea.uk>
Date: Wed, 21 Nov 2018 12:09:24 +0000
Subject: [PATCH 14/24] Tidy up Mesh::addCoordinates()

---
 src/mesh/mesh.cxx | 21 +++++++++------------
 1 file changed, 9 insertions(+), 12 deletions(-)

diff --git a/src/mesh/mesh.cxx b/src/mesh/mesh.cxx
index 604f33a3b8..ea34f9850a 100644
--- a/src/mesh/mesh.cxx
+++ b/src/mesh/mesh.cxx
@@ -349,22 +349,19 @@ void Mesh::addCoordinates(const CELL_LOC location, bool replace_coords) {
     } else {
       // Interpolate coordinates from CELL_CENTRE version
       ASSERT1(StaggerGrids); // If StaggerGrids==false, it doesn't make sense to have non-CELL_CENTRE Coordinates
-      if (!coords_map.count(location)) {
-        // location does not exist in coords_map, so create new entry
-        coords_map.emplace(location, std::unique_ptr<Coordinates>(new Coordinates(this, location, getCoordinates(CELL_CENTRE))));
-      } else if (replace_coords) {
-        // location does already exist in coords_map, so reset it
 
-        // first erase the existing entry to avoid throwing an exception from
-        // Coordinates::geometry().
-        // The check that would throw the exception is not needed because
-        // replacement of the Coordinates object has been explicitly requested.
-        coords_map.erase(location);
-        coords_map.emplace(location, std::unique_ptr<Coordinates>(new Coordinates(this, location, getCoordinates(CELL_CENTRE))));
-      } else {
+      if (!replace_coords and (coords_map.count(location) > 0)) {
         throw BoutException("Coordinates at %s already added to Mesh",
             CELL_LOC_STRING(location).c_str());
       }
+      if (replace_coords) {
+        // first erase the existing entry to avoid throwing an exception from
+        // Coordinates::geometry(): replacement of the Coordinates object has
+        // been explicitly requested.
+        coords_map.erase(location);
+      }
+
+      coords_map.emplace(location, std::unique_ptr<Coordinates>(new Coordinates(this, location, getCoordinates(CELL_CENTRE))));
     }
   }
 }

From 0f5ebf434d9581b028c8221789e4e09fdf1a5e1c Mon Sep 17 00:00:00 2001
From: John Omotani <john.omotani@ukaea.uk>
Date: Wed, 21 Nov 2018 12:10:50 +0000
Subject: [PATCH 15/24] Use defaultwarn("...") rather than print("Warning ...")

---
 tools/pylib/_boutcore_build/boutcore.pyx.in | 7 ++++---
 1 file changed, 4 insertions(+), 3 deletions(-)

diff --git a/tools/pylib/_boutcore_build/boutcore.pyx.in b/tools/pylib/_boutcore_build/boutcore.pyx.in
index 21deb239c1..49d17fe55f 100755
--- a/tools/pylib/_boutcore_build/boutcore.pyx.in
+++ b/tools/pylib/_boutcore_build/boutcore.pyx.in
@@ -59,6 +59,7 @@ cimport numpy as np
 #import atexit
 cimport resolve_enum as benum
 from libc.stdlib cimport malloc, free
+from boututils.boutwarnings import defaultwarn
 import copy
 EOF
 
@@ -751,9 +752,9 @@ cdef class Mesh:
         """
         Deprecated version of getCoordinates
         """
-        print("Warning Mesh.coordinates is deprecated, and does not handle "
-              "staggered grid locations. Use Mesh.getCoordinates(location) "
-              "method instead")
+        defaultwarn("Mesh.coordinates is deprecated, and does not handle "
+                    "staggered grid locations. Use "
+                    "Mesh.getCoordinates(location) method instead")
         return self.getCoordinates()
 
     def addCoordinates(self, location):

From 755ce7590f0e6668376fc804e018abe30357d783 Mon Sep 17 00:00:00 2001
From: John Omotani <john.omotani@ukaea.uk>
Date: Wed, 21 Nov 2018 12:14:15 +0000
Subject: [PATCH 16/24] Make allow_geometry_without_recalculate_staggered
 argument not option

Make the allow_geometry_without_recalculate_staggered flag an argument
to Coordinates::geometry() rather than an option that has to be read by
Mesh and passed to Coordinates.
---
 include/bout/coordinates.hxx |  3 +--
 src/mesh/coordinates.cxx     | 23 ++++++++++-------------
 src/mesh/mesh.cxx            |  1 -
 3 files changed, 11 insertions(+), 16 deletions(-)

diff --git a/include/bout/coordinates.hxx b/include/bout/coordinates.hxx
index 33179db8e9..d23365264a 100644
--- a/include/bout/coordinates.hxx
+++ b/include/bout/coordinates.hxx
@@ -93,7 +93,7 @@ public:
   Field2D IntShiftTorsion; ///< Integrated shear (I in BOUT notation)
 
   /// Calculate differential geometry quantities from the metric tensor
-  int geometry();
+  int geometry(bool allow_geometry_without_recalculate_staggered = false);
   int calcCovariant(); ///< Inverts contravatiant metric to get covariant
   int calcContravariant(); ///< Invert covariant metric to get contravariant
   int jacobian(); ///< Calculate J and Bxy
@@ -146,7 +146,6 @@ private:
   int nz; // Size of mesh in Z. This is mesh->ngz-1
   Mesh * localmesh;
   CELL_LOC location;
-  bool allow_geometry_without_recalculate_staggered;
 };
 
 /*
diff --git a/src/mesh/coordinates.cxx b/src/mesh/coordinates.cxx
index 9b8d1be823..4887c6c77c 100644
--- a/src/mesh/coordinates.cxx
+++ b/src/mesh/coordinates.cxx
@@ -26,9 +26,7 @@ Coordinates::Coordinates(Mesh *mesh)
       G1_23(mesh), G2_11(mesh), G2_22(mesh), G2_33(mesh), G2_12(mesh), G2_13(mesh),
       G2_23(mesh), G3_11(mesh), G3_22(mesh), G3_33(mesh), G3_12(mesh), G3_13(mesh),
       G3_23(mesh), G1(mesh), G2(mesh), G3(mesh), ShiftTorsion(mesh),
-      IntShiftTorsion(mesh), localmesh(mesh), location(CELL_CENTRE),
-      allow_geometry_without_recalculate_staggered(
-          mesh->allow_geometry_without_recalculate_staggered) {
+      IntShiftTorsion(mesh), localmesh(mesh), location(CELL_CENTRE) {
 
   if (mesh->get(dx, "dx")) {
     output_warn.write("\tWARNING: differencing quantity 'dx' not found. Set to 1.0\n");
@@ -306,7 +304,7 @@ void Coordinates::outputVars(Datafile &file) {
   file.add(J, "J", false);
 }
 
-int Coordinates::geometry() {
+int Coordinates::geometry(bool allow_geometry_without_recalculate_staggered) {
   TRACE("Coordinates::geometry");
 
   output_progress.write("Calculating differential geometry terms\n");
@@ -482,19 +480,18 @@ int Coordinates::geometry() {
       if (location == CELL_CENTRE) {
         throw BoutException("Coordinates::geometry() called at CELL_CENTRE, but "
           "other locations have already been added to the Mesh. These would "
-          "need recalculating. If possible, call Mesh::addCoordinates() after "
-          "this call to geometry(). To recalculate the other Coordinates objects "
-          "from the CELL_CENTRE Coordinates, pass true for the "
-          "recalculate_staggered argument, e.g.  coords->geometry(true). If you "
-          "need to recalculate multiple Coordinates objects explicitly, set "
-          "mesh:allow_geometry_without_recalculate_staggered=true in the input "
-          "file to disable this check.");
+          "need recalculating. If possible call Mesh::addCoordinates() after "
+          "this call to geometry(). If you need to recalculate multiple "
+          "Coordinates objects explicitly, call geometry(true) [i.e. setting the "
+          "argument allow_geometry_without_recalculate_staggered=true] file to "
+          "disable this check.");
       } else {
         throw BoutException("Coordinates::geometry() called at location %s, but "
           "this location has already been initialized and added to the Mesh. "
           "If you need to recalculate multiple Coordinates objects explicitly, "
-          "set mesh:allow_geometry_without_recalculate_staggered=true in the "
-          "input file to disable this check.", CELL_LOC_STRING(location).c_str());
+          "call geometry(true) [i.e. setting the argument "
+          "allow_geometry_without_recalculate_staggered=true] input file to "
+          "disable this check.", CELL_LOC_STRING(location).c_str());
       }
     }
 
diff --git a/src/mesh/mesh.cxx b/src/mesh/mesh.cxx
index ea34f9850a..16d863ca32 100644
--- a/src/mesh/mesh.cxx
+++ b/src/mesh/mesh.cxx
@@ -31,7 +31,6 @@ Mesh::Mesh(GridDataSource *s, Options* opt) : source(s), options(opt) {
   /// Get mesh options
   OPTION(options, StaggerGrids,   false); // Stagger grids
   OPTION(options, maxregionblocksize, MAXREGIONBLOCKSIZE);
-  OPTION(options, allow_geometry_without_recalculate_staggered, false);
   // Initialise derivatives
   derivs_init(options);  // in index_derivs.cxx for now
 }

From 64442061253e02426bf2a5f333087ca4026d48c3 Mon Sep 17 00:00:00 2001
From: John Omotani <john.omotani@ukaea.uk>
Date: Sun, 25 Nov 2018 13:16:29 +0000
Subject: [PATCH 17/24] Use getMesh() instead of fieldmesh in check for
 Coordinates

If Field3D was created in global scope, its fieldmesh pointer is null.
getMesh() returns the global mesh in that case, so use getMesh() instead
of fieldmesh to avoid segfaults.
---
 src/field/field3d.cxx | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/field/field3d.cxx b/src/field/field3d.cxx
index 022fd06aca..86f5e32f8c 100644
--- a/src/field/field3d.cxx
+++ b/src/field/field3d.cxx
@@ -243,7 +243,7 @@ void Field3D::setLocation(CELL_LOC new_location) {
 #if CHECK > 1
     // Check Coordinates for location have been added
     // For CHECK > 0, getCoordinates will throw if location has not been added.
-    fieldmesh->getCoordinates(location);
+    getMesh()->getCoordinates(location);
 #endif
 
     // Invalidate the coordinates pointer

From 032b4fd5b40918756d53778c18b207ba88ebebcb Mon Sep 17 00:00:00 2001
From: John Omotani <john.omotani@ukaea.uk>
Date: Tue, 27 Nov 2018 11:58:42 +0000
Subject: [PATCH 18/24] Use make_unique for Coordinates

---
 include/bout/mesh.hxx | 3 ---
 src/mesh/mesh.cxx     | 4 ++--
 2 files changed, 2 insertions(+), 5 deletions(-)

diff --git a/include/bout/mesh.hxx b/include/bout/mesh.hxx
index d957d0c18a..2f01ee7e33 100644
--- a/include/bout/mesh.hxx
+++ b/include/bout/mesh.hxx
@@ -765,9 +765,6 @@ class Mesh {
                            REGION region = RGN_NOBNDRY);
 
 private:
-  /// Allocates default Coordinates objects
-  std::unique_ptr<Coordinates> createDefaultCoordinates(const CELL_LOC location);
-
   //Internal region related information
   std::map<std::string, Region<Ind3D>> regionMap3D;
   std::map<std::string, Region<Ind2D>> regionMap2D;
diff --git a/src/mesh/mesh.cxx b/src/mesh/mesh.cxx
index e4ecc2923b..0a33d871ac 100644
--- a/src/mesh/mesh.cxx
+++ b/src/mesh/mesh.cxx
@@ -340,7 +340,7 @@ void Mesh::addCoordinates(const CELL_LOC location, bool replace_coords) {
       // Initialize coordinates from input
       if (!coords_map.count(location)) {
         // location does not exist in coords_map, so create new entry
-        coords_map.emplace(location, std::unique_ptr<Coordinates>(new Coordinates(this)));
+        coords_map.emplace(location, bout::utils::make_unique<Coordinates>(this));
       } else if (replace_coords) {
         // location does already exist in coords_map, so reset it
         coords_map.at(location).reset(new Coordinates(this));
@@ -360,7 +360,7 @@ void Mesh::addCoordinates(const CELL_LOC location, bool replace_coords) {
         coords_map.erase(location);
       }
 
-      coords_map.emplace(location, std::unique_ptr<Coordinates>(new Coordinates(this, location, getCoordinates(CELL_CENTRE))));
+      coords_map.emplace(location, bout::utils::make_unique<Coordinates>(this, location, getCoordinates(CELL_CENTRE)));
     }
   }
 }

From b5b80f24caae6b6c1cc6136c2159444015c2e2c2 Mon Sep 17 00:00:00 2001
From: John Omotani <john.omotani@ukaea.uk>
Date: Wed, 9 Jan 2019 19:40:43 +0000
Subject: [PATCH 19/24] Fix error message, replacing reference to
 'REQUEST_STAGGER'

Also similarly fix comments in Coordinates::geometry() and
BoutMesh::load().
---
 include/bout/mesh.hxx            | 4 ++--
 src/mesh/coordinates.cxx         | 6 +++---
 src/mesh/impls/bout/boutmesh.cxx | 6 +++---
 3 files changed, 8 insertions(+), 8 deletions(-)

diff --git a/include/bout/mesh.hxx b/include/bout/mesh.hxx
index 2c76177640..d8780deffc 100644
--- a/include/bout/mesh.hxx
+++ b/include/bout/mesh.hxx
@@ -437,8 +437,8 @@ class Mesh {
 #if CHECK > 0
     if (!coords_map.count(location)) {
       throw BoutException("Error: Coordinates for %s have not been added to "
-          "this Mesh. You should use REQUEST_STAGGER(location) before "
-          "initializing fields staggered to 'CELL_LOC location'.",
+          "this Mesh. You should call the method Mesh::addCoordinates(location) "
+          "before initializing fields staggered to 'CELL_LOC location'.",
           CELL_LOC_STRING(location).c_str());
     }
 #endif
diff --git a/src/mesh/coordinates.cxx b/src/mesh/coordinates.cxx
index 11075b37df..0dd32d2d33 100644
--- a/src/mesh/coordinates.cxx
+++ b/src/mesh/coordinates.cxx
@@ -487,9 +487,9 @@ int Coordinates::geometry() {
   if (localmesh->StaggerGrids) {
     // Replace Coordinates objects at staggered locations, if there are
     // enough grid points.
-    // This is a temporary workaround. In v4.3 we will users to call
-    // REQUEST_LOCATION(location) for each location that is needed and change
-    // this so that we don't waste memory on unneeded Coordinates.
+    // This is a temporary workaround. In v4.3 we will require users to call
+    // Mesh::addCoordinates(location) for each location that is needed and
+    // change this so that we don't waste memory on unneeded Coordinates.
     if (localmesh->LocalNx >= 4) {
       localmesh->addCoordinates(CELL_XLOW, true);
     }
diff --git a/src/mesh/impls/bout/boutmesh.cxx b/src/mesh/impls/bout/boutmesh.cxx
index 62b6a11979..54359833aa 100644
--- a/src/mesh/impls/bout/boutmesh.cxx
+++ b/src/mesh/impls/bout/boutmesh.cxx
@@ -845,9 +845,9 @@ int BoutMesh::load() {
   if (StaggerGrids) {
     // Add Coordinates objects at staggered locations, if there are enough grid
     // points.
-    // This is a temporary workaround. In v4.3 we will users to call
-    // REQUEST_LOCATION(location) for each location that is needed and remove
-    // this so that we don't waste memory on unneeded Coordinates.
+    // This is a temporary workaround. In v4.3 we will require users to call
+    // Mesh::addCoordinates(location) for each location that is needed and
+    // change this so that we don't waste memory on unneeded Coordinates.
     if (LocalNx >= 4) {
       addCoordinates(CELL_XLOW);
     }

From 0bf16afe372c6506ad714602c37b5929c515f545 Mon Sep 17 00:00:00 2001
From: John Omotani <john.omotani@ukaea.uk>
Date: Wed, 9 Jan 2019 19:42:27 +0000
Subject: [PATCH 20/24] Remove createDefaultCoordinates() method from header

Implementation was already removed from Mesh.
---
 include/bout/mesh.hxx | 3 ---
 1 file changed, 3 deletions(-)

diff --git a/include/bout/mesh.hxx b/include/bout/mesh.hxx
index d8780deffc..705281727f 100644
--- a/include/bout/mesh.hxx
+++ b/include/bout/mesh.hxx
@@ -755,9 +755,6 @@ class Mesh {
                            REGION region = RGN_NOBNDRY);
 
 private:
-  /// Allocates default Coordinates objects
-  std::unique_ptr<Coordinates> createDefaultCoordinates(const CELL_LOC location);
-
   //Internal region related information
   std::map<std::string, Region<Ind3D>> regionMap3D;
   std::map<std::string, Region<Ind2D>> regionMap2D;

From 75ce52256d2815fbd674afddfd4434270a045466 Mon Sep 17 00:00:00 2001
From: John Omotani <john.omotani@ukaea.uk>
Date: Wed, 9 Jan 2019 20:37:05 +0000
Subject: [PATCH 21/24] Revert change of Mesh::coords_map to use
 std::unique_ptr

---
 include/bout/mesh.hxx      | 2 +-
 src/mesh/mesh.cxx          | 4 ++--
 tests/unit/test_extras.hxx | 8 ++++----
 3 files changed, 7 insertions(+), 7 deletions(-)

diff --git a/include/bout/mesh.hxx b/include/bout/mesh.hxx
index 705281727f..f5b459d675 100644
--- a/include/bout/mesh.hxx
+++ b/include/bout/mesh.hxx
@@ -717,7 +717,7 @@ class Mesh {
   
   GridDataSource *source; ///< Source for grid data
   
-  std::map<CELL_LOC, std::unique_ptr<Coordinates> > coords_map; ///< Coordinate systems at different CELL_LOCs
+  std::map<CELL_LOC, std::shared_ptr<Coordinates> > coords_map; ///< Coordinate systems at different CELL_LOCs
 
   Options *options; ///< Mesh options section
   
diff --git a/src/mesh/mesh.cxx b/src/mesh/mesh.cxx
index 023beabccc..6941d10e33 100644
--- a/src/mesh/mesh.cxx
+++ b/src/mesh/mesh.cxx
@@ -340,7 +340,7 @@ void Mesh::addCoordinates(const CELL_LOC location, bool replace_coords) {
       // Initialize coordinates from input
       if (!coords_map.count(location)) {
         // location does not exist in coords_map, so create new entry
-        coords_map.emplace(location, std::unique_ptr<Coordinates>(new Coordinates(this)));
+        coords_map.emplace(location, std::make_shared<Coordinates>(this));
       } else if (replace_coords) {
         // location does already exists in coords_map, so reset it
         coords_map.at(location).reset(new Coordinates(this));
@@ -350,7 +350,7 @@ void Mesh::addCoordinates(const CELL_LOC location, bool replace_coords) {
       ASSERT1(StaggerGrids); // If StaggerGrids==false, it doesn't make sense to have non-CELL_CENTRE Coordinates
       if (!coords_map.count(location)) {
         // location does not exist in coords_map, so create new entry
-        coords_map.emplace(location, std::unique_ptr<Coordinates>(new Coordinates(this, location, getCoordinates(CELL_CENTRE))));
+        coords_map.emplace(location, std::make_shared<Coordinates>(this, location, getCoordinates(CELL_CENTRE)));
       } else if (replace_coords) {
         // location does already exists in coords_map, so reset it
         coords_map.at(location).reset(new Coordinates(this, location, getCoordinates(CELL_CENTRE)));
diff --git a/tests/unit/test_extras.hxx b/tests/unit/test_extras.hxx
index b143edef1f..b644a3bac8 100644
--- a/tests/unit/test_extras.hxx
+++ b/tests/unit/test_extras.hxx
@@ -74,10 +74,10 @@ public:
     IncIntShear = false;
     maxregionblocksize = MAXREGIONBLOCKSIZE;
 
-    coords_map.emplace(CELL_CENTRE, std::unique_ptr<Coordinates>(nullptr));
-    coords_map.emplace(CELL_XLOW, std::unique_ptr<Coordinates>(nullptr));
-    coords_map.emplace(CELL_YLOW, std::unique_ptr<Coordinates>(nullptr));
-    coords_map.emplace(CELL_ZLOW, std::unique_ptr<Coordinates>(nullptr));
+    coords_map.emplace(CELL_CENTRE, std::shared_ptr<Coordinates>(nullptr));
+    coords_map.emplace(CELL_XLOW, std::shared_ptr<Coordinates>(nullptr));
+    coords_map.emplace(CELL_YLOW, std::shared_ptr<Coordinates>(nullptr));
+    coords_map.emplace(CELL_ZLOW, std::shared_ptr<Coordinates>(nullptr));
   }
 
   comm_handle send(FieldGroup &UNUSED(g)) { return nullptr; };

From c8b4dd5f1e9bc5919043f0d75e2936337459dec7 Mon Sep 17 00:00:00 2001
From: John Omotani <john.omotani@ukaea.uk>
Date: Wed, 9 Jan 2019 21:00:02 +0000
Subject: [PATCH 22/24] Use std::unique_ptr for Mesh::coords_map

---
 include/bout/mesh.hxx      | 2 +-
 tests/unit/test_extras.hxx | 8 ++++----
 2 files changed, 5 insertions(+), 5 deletions(-)

diff --git a/include/bout/mesh.hxx b/include/bout/mesh.hxx
index 36d44dd2e0..674746f2c2 100644
--- a/include/bout/mesh.hxx
+++ b/include/bout/mesh.hxx
@@ -727,7 +727,7 @@ class Mesh {
   
   GridDataSource *source; ///< Source for grid data
   
-  std::map<CELL_LOC, std::shared_ptr<Coordinates> > coords_map; ///< Coordinate systems at different CELL_LOCs
+  std::map<CELL_LOC, std::unique_ptr<Coordinates> > coords_map; ///< Coordinate systems at different CELL_LOCs
 
   Options *options; ///< Mesh options section
   
diff --git a/tests/unit/test_extras.hxx b/tests/unit/test_extras.hxx
index b0745d7afc..dd19dc5421 100644
--- a/tests/unit/test_extras.hxx
+++ b/tests/unit/test_extras.hxx
@@ -74,10 +74,10 @@ public:
     IncIntShear = false;
     maxregionblocksize = MAXREGIONBLOCKSIZE;
 
-    coords_map.emplace(CELL_CENTRE, std::shared_ptr<Coordinates>(nullptr));
-    coords_map.emplace(CELL_XLOW, std::shared_ptr<Coordinates>(nullptr));
-    coords_map.emplace(CELL_YLOW, std::shared_ptr<Coordinates>(nullptr));
-    coords_map.emplace(CELL_ZLOW, std::shared_ptr<Coordinates>(nullptr));
+    coords_map.emplace(CELL_CENTRE, std::unique_ptr<Coordinates>(nullptr));
+    coords_map.emplace(CELL_XLOW, std::unique_ptr<Coordinates>(nullptr));
+    coords_map.emplace(CELL_YLOW, std::unique_ptr<Coordinates>(nullptr));
+    coords_map.emplace(CELL_ZLOW, std::unique_ptr<Coordinates>(nullptr));
   }
 
   comm_handle send(FieldGroup &UNUSED(g)) { return nullptr; };

From 2452c52ca482c2858c31d3e8aaf1253135c5dec1 Mon Sep 17 00:00:00 2001
From: John Omotani <john.omotani@ukaea.uk>
Date: Wed, 9 Jan 2019 21:11:01 +0000
Subject: [PATCH 23/24] Remove allow_geometry_without_recalculate_staggered
 from Mesh

This member variable was used for reading an option, but the interface
was changed to use a method argument passed by the user instead of an
option.
---
 include/bout/mesh.hxx | 3 ---
 1 file changed, 3 deletions(-)

diff --git a/include/bout/mesh.hxx b/include/bout/mesh.hxx
index 674746f2c2..c283f3f513 100644
--- a/include/bout/mesh.hxx
+++ b/include/bout/mesh.hxx
@@ -461,9 +461,6 @@ class Mesh {
     return coords_map.size();
   }
 
-  /// switch to pass to Coordinates objects
-  bool allow_geometry_without_recalculate_staggered;
-
   // First derivatives in index space
   // Implemented in src/mesh/index_derivs.hxx
 

From 4c58b834d09bcda83add479132f4626ca7fe601d Mon Sep 17 00:00:00 2001
From: John Omotani <john.omotani@ukaea.uk>
Date: Wed, 9 Jan 2019 22:13:29 +0000
Subject: [PATCH 24/24] Fix FakeMesh::setCoordinates() for std::unique_ptr

Using std::unique_ptr for Mesh::coords map requires a few changes to
FakeMesh.
---
 tests/unit/field/test_vector2d.cxx              | 2 +-
 tests/unit/field/test_vector3d.cxx              | 2 +-
 tests/unit/mesh/parallel/test_shiftedmetric.cxx | 2 +-
 tests/unit/test_extras.hxx                      | 4 ++--
 4 files changed, 5 insertions(+), 5 deletions(-)

diff --git a/tests/unit/field/test_vector2d.cxx b/tests/unit/field/test_vector2d.cxx
index a6ec890c34..28b8f09b5d 100644
--- a/tests/unit/field/test_vector2d.cxx
+++ b/tests/unit/field/test_vector2d.cxx
@@ -36,7 +36,7 @@ class Vector2DTest : public ::testing::Test {
     mesh->addBoundary(new BoundaryRegionYUp("upper_target", 1, nx - 2, mesh));
     mesh->addBoundary(new BoundaryRegionYDown("lower_target", 1, nx - 2, mesh));
 
-    dynamic_cast<FakeMesh*>(mesh)->setCoordinates(std::make_shared<Coordinates>(
+    dynamic_cast<FakeMesh*>(mesh)->setCoordinates(new Coordinates(
         mesh, Field2D{1.0}, Field2D{1.0}, BoutReal{1.0}, Field2D{1.0}, Field2D{0.0},
         Field2D{1.0}, Field2D{2.0}, Field2D{3.0}, Field2D{4.0}, Field2D{5.0},
         Field2D{6.0}, Field2D{1.0}, Field2D{2.0}, Field2D{3.0}, Field2D{4.0},
diff --git a/tests/unit/field/test_vector3d.cxx b/tests/unit/field/test_vector3d.cxx
index eae149755d..ed08c255e9 100644
--- a/tests/unit/field/test_vector3d.cxx
+++ b/tests/unit/field/test_vector3d.cxx
@@ -35,7 +35,7 @@ class Vector3DTest : public ::testing::Test {
     mesh->addBoundary(new BoundaryRegionYUp("upper_target", 1, nx - 2, mesh));
     mesh->addBoundary(new BoundaryRegionYDown("lower_target", 1, nx - 2, mesh));
 
-    dynamic_cast<FakeMesh*>(mesh)->setCoordinates(std::make_shared<Coordinates>(
+    dynamic_cast<FakeMesh*>(mesh)->setCoordinates(new Coordinates(
         mesh, Field2D{1.0}, Field2D{1.0}, BoutReal{1.0}, Field2D{1.0}, Field2D{0.0},
         Field2D{1.0}, Field2D{2.0}, Field2D{3.0}, Field2D{4.0}, Field2D{5.0},
         Field2D{6.0}, Field2D{1.0}, Field2D{2.0}, Field2D{3.0}, Field2D{4.0},
diff --git a/tests/unit/mesh/parallel/test_shiftedmetric.cxx b/tests/unit/mesh/parallel/test_shiftedmetric.cxx
index 4d8c2885a7..664fbfd28e 100644
--- a/tests/unit/mesh/parallel/test_shiftedmetric.cxx
+++ b/tests/unit/mesh/parallel/test_shiftedmetric.cxx
@@ -22,7 +22,7 @@ class ShiftedMetricTest : public ::testing::Test {
 
     fillField(zShift, {{1., 2., 3., 4., 5.}, {1., 2., 3., 4., 5.}, {1., 2., 3., 4., 5.}});
 
-    dynamic_cast<FakeMesh*>(mesh)->setCoordinates(std::make_shared<Coordinates>(
+    dynamic_cast<FakeMesh*>(mesh)->setCoordinates(new Coordinates(
         mesh, Field2D{1.0}, Field2D{1.0}, BoutReal{1.0}, Field2D{1.0}, Field2D{0.0},
         Field2D{1.0}, Field2D{1.0}, Field2D{1.0}, Field2D{0.0}, Field2D{0.0},
         Field2D{0.0}, Field2D{1.0}, Field2D{1.0}, Field2D{1.0}, Field2D{0.0},
diff --git a/tests/unit/test_extras.hxx b/tests/unit/test_extras.hxx
index 76af8d0638..e9fdc55b6f 100644
--- a/tests/unit/test_extras.hxx
+++ b/tests/unit/test_extras.hxx
@@ -92,8 +92,8 @@ public:
     coords_map.emplace(CELL_ZLOW, std::unique_ptr<Coordinates>(nullptr));
   }
 
-  void setCoordinates(std::shared_ptr<Coordinates> coords, CELL_LOC location = CELL_CENTRE) {
-    coords_map[location] = coords;
+  void setCoordinates(Coordinates* coords, CELL_LOC location = CELL_CENTRE) {
+    coords_map[location].reset(coords);
   }
 
   comm_handle send(FieldGroup &UNUSED(g)) { return nullptr; };