Separate initialization of ParallelTransform objects from constructors

include/bout/paralleltransform.hxx

@@ -101,7 +101,13 @@ public:
 class ShiftedMetric : public ParallelTransform {
 public:
   ShiftedMetric() = delete;
-  ShiftedMetric(Mesh &mesh);
+  ShiftedMetric(Mesh &m) : mesh(m) {}
+
+  /// initialize after constructor is finished, so that
+  /// localmesh->getCoordinateSystem() can call
+  /// ParallelTransform::getCoordinateSystem() for fields belonging to the
+  /// ParallelTransform
+  void initialize();
 
   /*!
    * Calculates the yup() and ydown() fields of f
@@ -151,6 +157,10 @@ private:
   Tensor<dcomplex> yupPhs;   ///< Cache of phase shifts for calculating yup fields
   Tensor<dcomplex> ydownPhs; ///< Cache of phase shifts for calculating ydown fields
 
+#if CHECK>0
+  bool isinitialized = false;
+#endif
+
   /*!
    * Shift a 2D field in Z. 
    * Since 2D fields are constant in Z, this has no effect

src/mesh/mesh.cxx

@@ -299,15 +299,17 @@ void Mesh::setParallelTransform() {
 
   }else if(ptstr == "shifted") {
     // Shifted metric method
-  transform = bout::utils::make_unique<ShiftedMetric>(*this);
+    transform = bout::utils::make_unique<ShiftedMetric>(*this);
+    static_cast<ShiftedMetric*>(transform.get())->initialize();
 
   }else if(ptstr == "fci") {
 
     Options *fci_options = Options::getRoot()->getSection("fci");
     // Flux Coordinate Independent method
     bool fci_zperiodic;
     fci_options->get("z_periodic", fci_zperiodic, true);
-    transform = bout::utils::make_unique<FCITransform>(*this, fci_zperiodic);
+    transform = bout::utils::make_unique<FCITransform>(*this);
+    static_cast<FCITransform*>(transform.get())->initialize(fci_zperiodic);
 
   }else {
     throw BoutException(_("Unrecognised paralleltransform option.\n"

src/mesh/parallel/fci.cxx

@@ -53,8 +53,11 @@ inline BoutReal sgn(BoutReal val) { return (BoutReal(0) < val) - (val < BoutReal
 
 // Calculate all the coefficients needed for the spline interpolation
 // dir MUST be either +1 or -1
-FCIMap::FCIMap(Mesh &mesh, int dir, bool zperiodic)
-  : dir(dir), boundary_mask(mesh), corner_boundary_mask(mesh), y_prime(&mesh) {
+void FCIMap::initialize(Mesh &mesh, int dir_in, bool zperiodic) {
+  dir = dir_in;
+  boundary_mask = BoutMask(mesh);
+  corner_boundary_mask = BoutMask(mesh);
+  y_prime = Field3D(0., &mesh);
 
   interp = InterpolationFactory::getInstance()->create(&mesh);
   interp->setYOffset(dir);
@@ -238,7 +241,7 @@ FCIMap::FCIMap(Mesh &mesh, int dir, bool zperiodic)
 
 const Field3D FCIMap::integrate(Field3D &f) const {
   TRACE("FCIMap::integrate");
-  
+
   // Cell centre values
   Field3D centre = interp->interpolate(f);
 
@@ -296,6 +299,8 @@ const Field3D FCIMap::integrate(Field3D &f) const {
 void FCITransform::calcYUpDown(Field3D &f) {
   TRACE("FCITransform::calcYUpDown");
 
+  ASSERT1(isinitialized);
+
   // Ensure that yup and ydown are different fields
   f.splitYupYdown();
 
@@ -306,7 +311,9 @@ void FCITransform::calcYUpDown(Field3D &f) {
 
 void FCITransform::integrateYUpDown(Field3D &f) {
   TRACE("FCITransform::integrateYUpDown");
-
+
+  ASSERT1(isinitialized);
+
   // Ensure that yup and ydown are different fields
   f.splitYupYdown();
 

src/mesh/parallel/fci.hxx

@@ -39,12 +39,10 @@ class FCIMap {
   /// Interpolation object
   Interpolation *interp;        // Cell centre
   Interpolation *interp_corner; // Cell corner at (x+1, z+1)
-
-  /// Private constructor - must be initialised with mesh
-  FCIMap();
 public:
+  /// initialize method must be called before using FCIMap
   /// dir MUST be either +1 or -1
-  FCIMap(Mesh& mesh, int dir, bool zperiodic);
+  void initialize(Mesh& mesh, int dir_in, bool zperiodic_in);
 
   int dir;                     /**< Direction of map */
 
@@ -65,9 +63,21 @@ public:
  */
 class FCITransform : public ParallelTransform {
 public:
-  FCITransform(Mesh &mesh, bool zperiodic = true)
-      : mesh(mesh), forward_map(mesh, +1, zperiodic), backward_map(mesh, -1, zperiodic),
-        zperiodic(zperiodic) {}
+  FCITransform(Mesh& m) : mesh(m) {}
+
+  /// initialize after constructor is finished, so that
+  /// localmesh->getCoordinateSystem() can call
+  /// ParallelTransform::getCoordinateSystem() for fields belonging to the
+  /// ParallelTransform
+  void initialize(bool zperiodic_in = true) {
+    forward_map.initialize(mesh, +1, zperiodic);
+    backward_map.initialize(mesh, -1, zperiodic);
+    zperiodic = zperiodic_in;
+
+#if CHECK>0
+    isinitialized = true;
+#endif
+  }
 
   void calcYUpDown(Field3D &f) override;
 
@@ -97,6 +107,10 @@ private:
   FCIMap backward_map;          /**< FCI map for field lines in -ve y */
 
   bool zperiodic;               /**< Is the z-direction periodic? */
+
+#if CHECK>0
+  bool isinitialized = false;
+#endif
 };
 
 #endif // __FCITRANSFORM_H__

src/mesh/parallel/shiftedmetric.cxx

@@ -15,9 +15,9 @@
 
 #include <output.hxx>
 
-ShiftedMetric::ShiftedMetric(Mesh &m) : mesh(m), zShift(&m) {
+void ShiftedMetric::initialize() {
   // Read the zShift angle from the mesh
-
+  zShift = Field2D(&mesh);
   if(mesh.get(zShift, "zShift")) {
     // No zShift variable. Try qinty in BOUT grid files
     mesh.get(zShift, "qinty");
@@ -76,13 +76,17 @@ ShiftedMetric::ShiftedMetric(Mesh &m) : mesh(m), zShift(&m) {
     }
   }
 
+#if CHECK > 1
+  isinitialized = true;
+#endif
 }
 
 /*!
  * Calculate the Y up and down fields
  */
 void ShiftedMetric::calcYUpDown(Field3D &f) {
   ASSERT1(&mesh == f.getMesh());
+  ASSERT1(isinitialized);
 
   f.splitYupYdown();
 
@@ -106,7 +110,9 @@ void ShiftedMetric::calcYUpDown(Field3D &f) {
  * and Y is then field aligned.
  */
 const Field3D ShiftedMetric::toFieldAligned(const Field3D &f, const REGION region) {
+  ASSERT1(isinitialized);
   ASSERT2(f.getCoordinateSystem() == COORDINATE_SYSTEM::Orthogonal);
+
   Field3D result = shiftZ(f, toAlignedPhs, region);
   result.setCoordinateSystem(COORDINATE_SYSTEM::FieldAligned);
   return result;
@@ -117,7 +123,9 @@ const Field3D ShiftedMetric::toFieldAligned(const Field3D &f, const REGION regio
  * but Y is not field aligned.
  */
 const Field3D ShiftedMetric::fromFieldAligned(const Field3D &f, const REGION region) {
+  ASSERT1(isinitialized);
   ASSERT2(f.getCoordinateSystem() == COORDINATE_SYSTEM::FieldAligned);
+
   Field3D result = shiftZ(f, fromAlignedPhs, region);
   result.setCoordinateSystem(f.getMesh()->getCoordinateSystem());
   return result;
@@ -126,6 +134,7 @@ const Field3D ShiftedMetric::fromFieldAligned(const Field3D &f, const REGION reg
 const Field3D ShiftedMetric::shiftZ(const Field3D& f, const Tensor<dcomplex>& phs,
                                     const REGION region) {
   ASSERT1(&mesh == f.getMesh());
+  ASSERT1(isinitialized);
   if(mesh.LocalNz == 1)
     return f; // Shifting makes no difference
 
@@ -141,6 +150,8 @@ const Field3D ShiftedMetric::shiftZ(const Field3D& f, const Tensor<dcomplex>& ph
 }
 
 void ShiftedMetric::shiftZ(const BoutReal* in, const dcomplex* phs, BoutReal* out) {
+  ASSERT1(isinitialized);
+
   Array<dcomplex> cmplx(nmodes);
 
   // Take forward FFT
@@ -160,6 +171,7 @@ void ShiftedMetric::shiftZ(const BoutReal* in, const dcomplex* phs, BoutReal* ou
 
 //Old approach retained so we can still specify a general zShift
 const Field3D ShiftedMetric::shiftZ(const Field3D &f, const Field2D &zangle, const REGION region) {
+  ASSERT1(isinitialized);
   ASSERT1(&mesh == f.getMesh());
   ASSERT1(f.getLocation() == zangle.getLocation());
   if(mesh.LocalNz == 1)
@@ -182,6 +194,8 @@ const Field3D ShiftedMetric::shiftZ(const Field3D &f, const Field2D &zangle, con
 }
 
 void ShiftedMetric::shiftZ(const BoutReal *in, int len, BoutReal zangle,  BoutReal *out) {
+  ASSERT1(isinitialized);
+
   int nmodes = len/2 + 1;
 
   // Complex array used for FFTs

tests/integrated/test-yupdown/test_yupdown.cxx

@@ -39,6 +39,7 @@ int main(int argc, char** argv) {
   BoutInitialise(argc, argv);
 
   ShiftedMetric s(*mesh);
+  s.initialize();
 
   // Read variable from mesh
   Field3D var;