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 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++) diff --git a/include/bout/coordinates.hxx b/include/bout/coordinates.hxx index 8f1420bca4..f10086d1fe 100644 --- a/include/bout/coordinates.hxx +++ b/include/bout/coordinates.hxx @@ -102,7 +102,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 diff --git a/include/bout/mesh.hxx b/include/bout/mesh.hxx index 326c2f2fa3..650de1e14c 100644 --- a/include/bout/mesh.hxx +++ b/include/bout/mesh.hxx @@ -64,6 +64,7 @@ class Mesh; #include "sys/range.hxx" // RangeIterator #include +#include #include "coordinates.hxx" // Coordinates class @@ -438,15 +439,30 @@ 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 call the method Mesh::addCoordinates(location) " + "before initializing fields staggered to 'CELL_LOC location'.", + CELL_LOC_STRING(location).c_str()); } +#endif + return coords_map.at(location).get(); + } + + /// 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); + + /// 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(); } /// Returns the non-CELL_CENTRE location @@ -790,7 +806,7 @@ class Mesh { GridDataSource *source; ///< Source for grid data - std::map > coords_map; ///< Coordinate systems at different CELL_LOCs + std::map > coords_map; ///< Coordinate systems at different CELL_LOCs Options *options; ///< Mesh options section @@ -808,9 +824,6 @@ class Mesh { private: - /// Allocates default Coordinates objects - std::shared_ptr createDefaultCoordinates(const CELL_LOC location); - //Internal region related information std::map> regionMap3D; std::map> regionMap2D; 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 diff --git a/src/field/field3d.cxx b/src/field/field3d.cxx index 65511f428a..1aeed497f4 100644 --- a/src/field/field3d.cxx +++ b/src/field/field3d.cxx @@ -209,6 +209,12 @@ void Field3D::setLocation(CELL_LOC new_location) { new_location = CELL_CENTRE; } +#if CHECK > 1 + // Check Coordinates for location have been added + // For CHECK > 0, getCoordinates will throw if location has not been added. + getMesh()->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 d6ce7a38a0..ef2aaf4915 100644 --- a/src/mesh/coordinates.cxx +++ b/src/mesh/coordinates.cxx @@ -324,7 +324,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"); @@ -462,13 +462,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 = bout::derivatives::index::DDX(1. / dx); // d/di(1/dx) } else { + // Shift d2x to our location + d2x = interp_to(d2x, location); + d1_dx = -d2x / (dx * dx); } @@ -477,9 +480,43 @@ int Coordinates::geometry() { "\tWARNING: differencing quantity 'd2y' not found. Calculating from dy\n"); d1_dy = bout::derivatives::index::DDY(1. / dy); // d/di(1/dy) } else { + // Shift d2y to our location + d2y = interp_to(d2y, location); + d1_dy = -d2y / (dy * dy); } + 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(). 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, " + "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()); + } + } + + } + return 0; } diff --git a/src/mesh/impls/bout/boutmesh.cxx b/src/mesh/impls/bout/boutmesh.cxx index 0bf59de09a..f274badfc5 100644 --- a/src/mesh/impls/bout/boutmesh.cxx +++ b/src/mesh/impls/bout/boutmesh.cxx @@ -839,7 +839,10 @@ int BoutMesh::load() { // Add boundary regions addBoundaryRegions(); - output_info.write(_("\tdone\n")); + // Create CELL_CENTRE Coordinates object + addCoordinates(CELL_CENTRE); + + output_info.write("\tdone\n"); return 0; } @@ -2208,7 +2211,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(0, xstart-1, ystart, yend, 0, LocalNz-1, LocalNy, LocalNz, maxregionblocksize)); addRegion2D("RGN_INNER_X", Region(0, xstart-1, ystart, yend, 0, 0, @@ -2225,7 +2228,7 @@ void BoutMesh::addBoundaryRegions() { } // Outer X - if(mesh->firstX() && !mesh->periodicX) { + if(firstX() && !periodicX) { addRegion3D("RGN_OUTER_X", Region(xend+1, LocalNx-1, ystart, yend, 0, LocalNz-1, LocalNy, LocalNz, maxregionblocksize)); addRegion2D("RGN_OUTER_X", Region(xend+1, LocalNx-1, ystart, yend, 0, 0, diff --git a/src/mesh/mesh.cxx b/src/mesh/mesh.cxx index 941bfc6580..9b52124d73 100644 --- a/src/mesh/mesh.cxx +++ b/src/mesh/mesh.cxx @@ -325,13 +325,44 @@ ParallelTransform& Mesh::getParallelTransform() { return *transform; } -std::shared_ptr Mesh::createDefaultCoordinates(const CELL_LOC location) { - if (location == CELL_CENTRE || location == CELL_DEFAULT) - // Initialize coordinates from input - return std::make_shared(this); - else - // Interpolate coordinates from CELL_CENTRE version - return std::make_shared(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, bout::utils::make_unique(this)); + } else if (replace_coords) { + // location does already exist 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 (!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, bout::utils::make_unique(this, location, getCoordinates(CELL_CENTRE))); + } + } } 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 41a9b04573..0078e7341a 100755 --- a/tests/MMS/upwinding3/runtest +++ b/tests/MMS/upwinding3/runtest @@ -34,6 +34,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 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-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') diff --git a/tests/integrated/test-drift-instability/2fluid.cxx b/tests/integrated/test-drift-instability/2fluid.cxx index 75bade9a81..43f082890c 100644 --- a/tests/integrated/test-drift-instability/2fluid.cxx +++ b/tests/integrated/test-drift-instability/2fluid.cxx @@ -298,6 +298,7 @@ int physics_init(bool UNUSED(restarting)) { dump.add(wci, "wci", 0); if (mesh->StaggerGrids) { + mesh->addCoordinates(CELL_YLOW); maybe_ylow = CELL_YLOW; } else { maybe_ylow = CELL_CENTRE; 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(mesh)->setCoordinates(std::make_shared( + dynamic_cast(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(mesh)->setCoordinates(std::make_shared( + dynamic_cast(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(mesh)->setCoordinates(std::make_shared( + dynamic_cast(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/mesh/test_boutmesh.cxx b/tests/unit/mesh/test_boutmesh.cxx index 3891227e48..967eda0773 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 class FakeGridDataSource : public GridDataSource { public: FakeGridDataSource(){}; ~FakeGridDataSource(){}; bool hasVar(const std::string &UNUSED(name)) { return false; }; - bool get(Mesh *UNUSED(m), int &UNUSED(ival), const std::string &UNUSED(name)) { - return true; + bool get(Mesh *UNUSED(m), int &ival, const std::string &name) { + if (intvars.count(name)>0) { + ival = intvars.at(name); + return true; + } + return false; }; bool get(Mesh *UNUSED(m), BoutReal &UNUSED(rval), const std::string &UNUSED(name)) { - return true; + return false; } bool get(Mesh *UNUSED(m), Field2D &UNUSED(var), const std::string &UNUSED(name), BoutReal UNUSED(def) = 0.0) { - return true; + return false; } bool get(Mesh *UNUSED(m), Field3D &UNUSED(var), const std::string &UNUSED(name), BoutReal UNUSED(def) = 0.0) { - return true; + return false; } bool get(Mesh *UNUSED(m), std::vector &UNUSED(var), const std::string &UNUSED(name), int UNUSED(len), int UNUSED(offset) = 0, Direction UNUSED(dir) = GridDataSource::X) { - return true; + return false; } bool get(Mesh *UNUSED(m), std::vector &UNUSED(var), const std::string &UNUSED(name), int UNUSED(len), int UNUSED(offset) = 0, Direction UNUSED(dir) = GridDataSource::X) { - return true; + return false; } + + std::unordered_map 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); } diff --git a/tests/unit/test_extras.hxx b/tests/unit/test_extras.hxx index a4593c8657..e9fdc55b6f 100644 --- a/tests/unit/test_extras.hxx +++ b/tests/unit/test_extras.hxx @@ -85,10 +85,15 @@ public: StaggerGrids = false; IncIntShear = false; maxregionblocksize = MAXREGIONBLOCKSIZE; + + coords_map.emplace(CELL_CENTRE, std::unique_ptr(nullptr)); + coords_map.emplace(CELL_XLOW, std::unique_ptr(nullptr)); + coords_map.emplace(CELL_YLOW, std::unique_ptr(nullptr)); + coords_map.emplace(CELL_ZLOW, std::unique_ptr(nullptr)); } - void setCoordinates(std::shared_ptr 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; }; diff --git a/tools/pylib/_boutcore_build/boutcore.pyx.in b/tools/pylib/_boutcore_build/boutcore.pyx.in index 4c67e1a5f9..e8af7a8c73 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 cdef extern from "boutexception_helper.hxx": @@ -630,7 +631,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): """ @@ -658,7 +658,6 @@ cdef class Mesh: self.isGlobal=False self.isNormalised=-1 self.factory = FieldFactory() - self._coords = None if create: if options: opt = (options).cobj @@ -743,14 +742,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 + """ + 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): + """ + Initialize Coordinates object at location """ - if self._coords is None: - self._coords = coordsFromObj(self.cobj.getCoordinates()) - 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 f5342511c3..727e4939d4 100755 --- a/tools/pylib/_boutcore_build/boutcpp.pxd.in +++ b/tools/pylib/_boutcore_build/boutcpp.pxd.in @@ -64,7 +64,8 @@ cdef extern from "bout/mesh.hxx": int ystart int LocalNx int LocalNy - Coordinates * getCoordinates() + Coordinates * getCoordinates(benum.CELL_LOC location) except + + void addCoordinates(benum.CELL_LOC location) except + cdef extern from "bout/coordinates.hxx": cppclass Coordinates: