Cylindrical to Cartesian

include/bout/tokamak_coordinates.hxx

@@ -11,18 +11,31 @@ using FieldMetric = MetricTensor::FieldMetric;
 
 namespace bout {
 
+    struct Coordinates3D {
+
+        Field3D x;
+        Field3D y;
+        Field3D z;
+
+        Coordinates3D(Field3D& x, Field3D& y, Field3D& z) : x(x), y(y), z(z) {}
+    };
+
     struct TokamakOptions {
         Field2D Rxy;
+        Field2D Zxy;
         Field2D Bpxy;
         Field2D Btxy;
         Field2D Bxy;
         Field2D hthe;
         FieldMetric I;
         FieldMetric dx;
+        std::vector<double> toroidal_angles;
 
-        TokamakOptions(Mesh &mesh);
+        TokamakOptions(Mesh& mesh);
 
         void normalise(BoutReal Lbar, BoutReal Bbar, BoutReal ShearFactor);
+
+        Coordinates3D CylindricalCoordinatesToCartesian();
     };
 
     BoutReal get_sign_of_bp(const Field2D &Bpxy);

src/mesh/tokamak_coordinates.cxx

@@ -4,6 +4,7 @@
 #include "bout/bout_types.hxx"
 #include "bout/field2d.hxx"
 #include "bout/utils.hxx"
+#include "bout/constants.hxx"
 
 
 namespace bout {
@@ -17,7 +18,7 @@ namespace bout {
 
     TokamakOptions::TokamakOptions(Mesh &mesh) {
         mesh.get(Rxy, "Rxy");
-        //    mesh->get(Zxy, "Zxy");
+        mesh.get(Zxy, "Zxy");
         mesh.get(Bpxy, "Bpxy");
         mesh.get(Btxy, "Btxy");
         mesh.get(Bxy, "Bxy");
@@ -26,6 +27,34 @@ namespace bout {
         if (mesh.get(dx, "dpsi")) {
             dx = mesh.getCoordinates()->dx();
         }
+//        mesh.get(toroidal_angle, "z");
+        const auto d_phi = TWOPI / mesh.LocalNz;
+        auto current_phi = 0.0;
+        for (int k = 0; k < mesh.LocalNz; k++) {
+            toroidal_angles.push_back(current_phi);
+            current_phi += d_phi;
+        }
+    }
+
+    Coordinates3D TokamakOptions::CylindricalCoordinatesToCartesian() {
+
+        auto* mesh = Rxy.getMesh();
+        Field3D x = Field3D(0.0, mesh);
+        Field3D y = Field3D(0.0, mesh);
+        Field3D z = Field3D(0.0, mesh);
+
+        for (int i = 0; i < Rxy.getNx(); i++) {
+            for (int j = 0; j < Rxy.getNy(); j++) {
+                int k = 0;
+                for (auto angle : toroidal_angles) {
+                    x(i, j, k) = Rxy(i, j) * std::cos(angle);
+                    y(i, j, k) = Rxy(i, j) * std::sin(angle);
+                    z(i, j, k) = Zxy(i, j);
+                    k++;
+                }
+            }
+        }
+        return Coordinates3D(x, y, z);
     }
 
     void TokamakOptions::normalise(BoutReal Lbar, BoutReal Bbar, BoutReal ShearFactor) {
@@ -47,7 +76,7 @@ namespace bout {
 
         const BoutReal sign_of_bp = get_sign_of_bp(tokamak_options.Bpxy);
 
-        auto *coord = mesh.getCoordinates();
+        auto* coord = mesh.getCoordinates();
 
         const auto g11 = SQ(tokamak_options.Rxy * tokamak_options.Bpxy);
         const auto g22 = 1.0 / SQ(tokamak_options.hthe);

tests/unit/CMakeLists.txt

@@ -67,6 +67,7 @@ set(serial_tests_source
   ./mesh/test_boundary_factory.cxx
   ./mesh/test_boutmesh.cxx
   ./mesh/test_coordinates.cxx
+  ./mesh/test_change_coordinate_system.cxx
   ./mesh/test_coordinates_accessor.cxx
   ./mesh/test_interpolation.cxx
   ./mesh/test_invert3x3.cxx

tests/unit/fake_mesh_fixture.hxx

@@ -25,15 +25,19 @@
 /// alias to make a new test:
 ///
 ///     using MyTest = FakeMeshFixture;
-class FakeMeshFixture : public ::testing::Test {
+///
+///     Type alias FakeMeshFixture = FakeMeshFixture_tmpl<3, 5, 7>;
+///     is used as a shim to allow FakeMeshFixture to be used with default values for nx, ny, nz
+template<int NX, int NY, int NZ>
+class FakeMeshFixture_tmpl : public ::testing::Test {
 public:
-    FakeMeshFixture() {
+    FakeMeshFixture_tmpl() {
         WithQuietOutput quiet_info{output_info};
         WithQuietOutput quiet_warn{output_warn};
 
         delete bout::globals::mesh;
         bout::globals::mpi = new MpiWrapper();
-        bout::globals::mesh = new FakeMesh(nx, ny, nz);
+        bout::globals::mesh = new FakeMesh(NX, NY, NZ);
         bout::globals::mesh->createDefaultRegions();
         static_cast<FakeMesh*>(bout::globals::mesh)->setCoordinates(nullptr);
         test_coords = std::make_shared<Coordinates>(
@@ -71,7 +75,7 @@ public:
         dynamic_cast<FakeMesh*>(bout::globals::mesh)->createBoundaryRegions();
 
         delete mesh_staggered;
-        mesh_staggered = new FakeMesh(nx, ny, nz);
+        mesh_staggered = new FakeMesh(NX, NY, NZ);
         mesh_staggered->StaggerGrids = true;
         dynamic_cast<FakeMesh*>(mesh_staggered)->setCoordinates(nullptr);
         dynamic_cast<FakeMesh*>(mesh_staggered)->setCoordinates(nullptr, CELL_XLOW);
@@ -123,7 +127,7 @@ public:
                 ->setCoordinates(test_coords_staggered, CELL_ZLOW);
     }
 
-    ~FakeMeshFixture() override {
+    ~FakeMeshFixture_tmpl() override {
         delete bout::globals::mesh;
         bout::globals::mesh = nullptr;
         delete mesh_staggered;
@@ -134,12 +138,14 @@ public:
         Options::cleanup();
     }
 
-    static constexpr int nx = 3;
-    static constexpr int ny = 5;
-    static constexpr int nz = 7;
+    static constexpr int nx = NX;
+    static constexpr int ny = NY;
+    static constexpr int nz = NZ;
 
     Mesh* mesh_staggered = nullptr;
 
     std::shared_ptr<Coordinates> test_coords{nullptr};
     std::shared_ptr<Coordinates> test_coords_staggered{nullptr};
-};
+};
+
+using FakeMeshFixture = FakeMeshFixture_tmpl<3, 5, 7>;

tests/unit/mesh/test_change_coordinate_system.cxx

@@ -0,0 +1,77 @@
+#include "gtest/gtest.h"
+#include "bout/coordinates.hxx"
+#include "bout/mesh.hxx"
+#include "fake_mesh_fixture.hxx"
+#include "bout/constants.hxx"
+#include <bout/tokamak_coordinates.hxx>
+
+static constexpr int NX = 3;
+static constexpr int NY = 5;
+static constexpr int NZ = 8;
+
+using bout::globals::mesh;
+
+class CoordinateTransformTest : public FakeMeshFixture_tmpl<NX, NY, NZ> {
+public:
+    using FieldMetric = Coordinates::FieldMetric;
+
+    CoordinateTransformTest() : FakeMeshFixture_tmpl() {}
+};
+
+TEST_F(CoordinateTransformTest, CylindricalToCartesian) {
+
+    // arrange
+
+    // Set up test values
+    // Calculate cylindrical coordinates (Rxy, Zxy)
+    // from (2D) orthogonal poloidal coordinates (r, theta)
+
+    const double R0 = 2.0;  // major radius
+    const std::array<double, NX> r_values = {0.1, 0.2, 0.3};  // minor radius
+    const std::array<double, NY> theta_values = {  // poloidal angle
+            0.0,
+            PI / 2,
+            PI,
+            3 * PI / 2,
+            2 * PI};
+
+    auto tokamak_options = bout::TokamakOptions(*mesh);
+
+    int i = 0;
+    for (auto r: r_values) {
+        int j = 0;
+        for (auto theta: theta_values) {
+            tokamak_options.Rxy(i, j) = R0 + r * std::cos(theta);
+            tokamak_options.Zxy(i, j) = r * std::sin(theta);
+            j++;
+        }
+        i++;
+    }
+
+    // act
+    bout::Coordinates3D cartesian_coords = tokamak_options.CylindricalCoordinatesToCartesian();
+
+    // assert
+    const auto max_r = *std::max_element(begin(r_values), end(r_values));
+    const auto expected_max_x = R0 + max_r;
+    const auto expected_max_y = (R0 + max_r);
+    const auto expected_max_z = max_r;
+
+    const auto expected_min_x = -1 * (R0 + max_r);
+    const auto expected_min_y = -1 * (R0 + max_r);
+    const auto expected_min_z = -1 * expected_max_z;
+
+    const auto actual_max_x = max(cartesian_coords.x, false, "RGN_ALL");
+    const auto actual_max_y = max(cartesian_coords.y, false, "RGN_ALL");
+    const auto actual_max_z = max(cartesian_coords.z, false, "RGN_ALL");
+    const auto actual_min_x = min(cartesian_coords.x, false, "RGN_ALL");
+    const auto actual_min_y = min(cartesian_coords.y, false, "RGN_ALL");
+    const auto actual_min_z = min(cartesian_coords.z, false, "RGN_ALL");
+
+    EXPECT_EQ(expected_max_x, actual_max_x);
+    EXPECT_EQ(expected_max_y, actual_max_y);
+    EXPECT_EQ(expected_max_z, actual_max_z);
+    EXPECT_EQ(expected_min_x, actual_min_x);
+    EXPECT_EQ(expected_min_y, actual_min_y);
+    EXPECT_EQ(expected_min_z, actual_min_z);
+}