Add output of RBC barycentres

Code/configuration/SimConfig.h

@@ -198,12 +198,18 @@ namespace hemelb::configuration
         std::size_t exponent;
     };
 
+    struct CellOutputConfig {
+        LatticeTimeStep output_period;
+        bool physical_units = false;
+    };
+
     struct RBCConfig {
         LatticeDistance boxSize;
         std::map<std::string, TemplateCellConfig> meshes;
         NodeForceConfig cell2cell;
         NodeForceConfig cell2wall;
-        LatticeTimeStep output_period;
+        std::optional<CellOutputConfig> full_output;
+        std::optional<CellOutputConfig> summary_output;
     };
 
     class SimConfig

Code/configuration/SimConfigReader.cc

@@ -958,10 +958,19 @@ namespace hemelb::configuration {
             throw Exception() << "Box-size < cell-cell interaction size: "
                                  "cell-cell interactions cannot be all accounted for.";
 
-        ans.output_period = PopDimensionalValue<LatticeTimeStep>(
-                *rbcEl.PopChildOrThrow("output")->PopChildOrThrow("period"),
-                "lattice"
-        );
+        auto outEl = rbcEl.PopChildOrThrow("output");
+        if (auto vtkEl = outEl->PopChildOrNull("vtk")) {
+            ans.full_output = CellOutputConfig{
+                    .output_period = PopDimensionalValue<LatticeTimeStep>(*vtkEl->PopChildOrThrow("period"), "lattice"),
+                    .physical_units = PopDimensionalValue<bool>(*vtkEl->PopChildOrThrow("physical"), "bool")
+            };
+        }
+        if (auto summEl = outEl->PopChildOrNull("summary")) {
+            ans.summary_output = CellOutputConfig{
+                    .output_period = PopDimensionalValue<LatticeTimeStep>(*summEl->PopChildOrThrow("period"), "lattice"),
+                    .physical_units = PopDimensionalValue<bool>(*summEl->PopChildOrThrow("physical"), "bool")
+            };
+        }
 
         return ans;
     }

Code/io/PathManager.cc

@@ -59,26 +59,24 @@ namespace hemelb::io
         return cpDir;
     }
 
-    fs::path PathManager::GetRBCOutputPathWithSubdir(std::string const& subdirectoryName) const
+    fs::path PathManager::GetRBCOutputPathWithSubdir(std::string const& subdirectoryName, bool create, bool allow_existing) const
     {
-      auto rbcSubdir = rbcDir / subdirectoryName;
-
-      if (doIo)
-      {
-        if (fs::exists(rbcSubdir))
-          throw Exception() << "Output directory '" << rbcSubdir << "' already exists.";
-
-        std::error_code ec;
-        if (!fs::create_directory(rbcSubdir, rbcDir, ec))
-        if (ec)
-        {
-          throw Exception() << "Output directory '" << rbcSubdir << "' could not be created because: "
-                            << ec.message();
+        auto rbcSubdir = rbcDir / subdirectoryName;
+        if (doIo) {
+            bool does_exist = fs::exists(rbcSubdir);
+            if (create) {
+                if (does_exist && !allow_existing)
+                    throw Exception() << "Output directory '" << rbcSubdir << "' already exists.";
+
+                std::error_code ec;
+                if (!fs::create_directory(rbcSubdir, rbcDir, ec))
+                    if (ec) {
+                        throw Exception() << "Output directory '" << rbcSubdir << "' could not be created because: "
+                                          << ec.message();
+                    }
+            }
         }
-      }
-
-      return rbcSubdir;
+        return rbcSubdir;
     }
 
 }
-

Code/io/PathManager.h

@@ -67,7 +67,7 @@ namespace hemelb::io
          * @param subdirectoryName Name of the subdirectory to be created
          * @return Path to the newly created subdirectory
          */
-        [[nodiscard]] path GetRBCOutputPathWithSubdir(std::string const& subdirectoryName) const;
+        [[nodiscard]] path GetRBCOutputPathWithSubdir(std::string const& subdirectoryName, bool create, bool allow_existing) const;
 
       private:
         path outputDir;

Code/io/formats/rbc.h

@@ -0,0 +1,32 @@
+// This file is part of HemeLB and is Copyright (C)
+// the HemeLB team and/or their institutions, as detailed in the
+// file AUTHORS. This software is provided under the terms of the
+// license in the file LICENSE.
+
+#ifndef HEMELB_IO_FORMATS_RBC_H
+#define HEMELB_IO_FORMATS_RBC_H
+
+#include <cstddef>
+#include <cstdint>
+#include "io/formats/formats.h"
+
+// Information on per-cell summary data
+namespace hemelb::io::formats::rbc
+{
+    // Magic number to identify extraction data files.
+    // ASCII for 'rbc' + EOF
+    static constexpr std::uint32_t MagicNumber = 0x72626304;
+
+    // The version number of the file format.
+    static constexpr std::uint32_t VersionNumber = 1;
+
+    // Header is 4 x uint32_t:
+    // HemeLB magic, RBC magic, version, number of rows
+    constexpr std::size_t header_size = 16U;
+
+    // Each row is:
+    //  - string representation of UUID (36 bytes)
+    //  - cell barycentre (three doubles)
+    constexpr std::size_t row_size = 36 + 3*8;
+}
+#endif

Code/io/readers/XdrReader.h

@@ -7,64 +7,86 @@
 #define HEMELB_IO_READERS_XDRREADER_H
 
 #include <cstdint>
+#include <span>
 #include <string>
 
 #include "Exception.h"
 #include "io/XdrSerialisation.h"
+#include "util/Vector3D.h"
 
 namespace hemelb::io
 {
-  // Base class to read XDR data. This does the XDR related
-  // stuff. Derived classes must implement GetPosition and get_bytes
-  // member functions to allow it to work.
-  class XdrReader
-  {
-  public:
-    // Virtual destructor.
-    virtual ~XdrReader() = default;
+    // Base class to read XDR data. This does the XDR related
+    // stuff. Derived classes must implement GetPosition and get_bytes
+    // member functions to allow it to work.
+    class XdrReader
+    {
+    public:
+        // Virtual destructor.
+        virtual ~XdrReader() = default;
 
-    // Main function for reading the next bit of the stream.
-    //
-    // Stores result of deserialising in the argument supplied and
-    // indicates success with return value.
-    template<class T>
-    bool read(T& val) {
-      constexpr auto n = xdr::xdr_serialised_size<T>();
-      auto buf = get_bytes(n);
-      xdr::xdr_deserialise(val, buf);
-      return true;
-    }
+        // Main function for reading the next bit of the stream.
+        //
+        // Stores result of deserialising in the argument supplied and
+        // indicates success with return value.
+        template<class T>
+        bool read(T& val) {
+            constexpr auto n = xdr::xdr_serialised_size<T>();
+            auto buf = get_bytes(n);
+            xdr::xdr_deserialise(val, buf);
+            return true;
+        }
 
-    // Return result of deserialising - will throw on error.
-    template <class T>
-    T read() {
-      T ans;
-      if (!read<T>(ans)) {
-	throw Exception() << "Error reading type from XDR";
-      }
-      return ans;
-    }
+        // Like xdr_opaque
+        template<std::size_t N>
+        bool read(std::span<std::byte, N>& bytes);
+
+        template <typename T>
+        bool read(util::Vector3D<T>& vec) {
+            for (int i = 0; i < 3; ++i)
+                read(vec[i]);
+            return true;
+        }
 
-    // Get the position in the stream.
-    virtual unsigned GetPosition() = 0;
+        // Return result of deserialising - will throw on error.
+        template <class T>
+        T read() {
+            T ans;
+            if (!read<T>(ans)) {
+                throw Exception() << "Error reading type from XDR";
+            }
+            return ans;
+        }
 
-  protected:
-    // Get some bytes from the underlying storage
-    virtual const std::byte* get_bytes(size_t n) = 0;
-  };
+        // Get the position in the stream.
+        virtual unsigned GetPosition() = 0;
 
-  // Specialisation for strings
-  template<>
-  inline bool XdrReader::read(std::string& val) {
-    uint32_t len = 0;
-    read(len);
-    // p == padded
-    uint32_t plen = 4 * ((len - 1)/4 + 1);
-    auto pstr = reinterpret_cast<const char*>(get_bytes(plen));
-    val.assign(pstr, len);
-    return true;
-  }
+    protected:
+        // Get some bytes from the underlying storage
+        virtual const std::byte* get_bytes(size_t n) = 0;
+    };
 
+    // Specialisation for strings
+    template<>
+    inline bool XdrReader::read(std::string& val) {
+        uint32_t len = 0;
+        read(len);
+        // p == padded
+        uint32_t plen = 4 * ((len - 1)/4 + 1);
+        auto pstr = reinterpret_cast<const char*>(get_bytes(plen));
+        val.assign(pstr, len);
+        return true;
+    }
+    // Like xdr opaque
+    template<std::size_t N>
+    inline bool XdrReader::read(std::span<std::byte, N>& bytes) {
+        const auto len = bytes.size();
+        const auto req_nwords = (len - 1)/4 + 1;
+        const auto space = req_nwords * 4;
+        auto data = get_bytes(space);
+        std::copy(data, data + len, bytes.begin());
+        return true;
+    }
 }
 
 #endif // HEMELB_IO_READERS_XDRREADER_H

Code/io/xml.cc

@@ -295,7 +295,7 @@ namespace hemelb::io::xml
 
     std::ostringstream Element::MakeAttributeStream() const {
         auto doc = static_cast<Document const*>(el->GetDocument()->GetUserData());
-        return doc->attr_stream_factory();
+        return doc->attr_stream_factory() << std::boolalpha;
     }
 
     Element::operator bool() const {

Code/io/xml.h

@@ -759,7 +759,8 @@ namespace hemelb::io::xml
         std::istringstream attrStream(s, std::ios_base::in);
 
         // Don't skip whitespace as that could indicate a malformed value
-        attrStream >> std::noskipws;
+        // Read bools as true/false
+        attrStream >> std::noskipws >> std::boolalpha;
 
         attrStream >> out;
         if (attrStream.fail())

Code/net/MpiCommunicator.cc

@@ -81,6 +81,12 @@ namespace hemelb::net
       HEMELB_MPI_CALL(MPI_Barrier, (*commPtr));
     }
 
+    MpiRequest MpiCommunicator::Ibarrier() const {
+        MpiRequest ans;
+        HEMELB_MPI_CALL(MPI_Ibarrier, (*commPtr, &ans.req));
+        return ans;
+    }
+
     MpiGroup MpiCommunicator::Group() const
     {
       MPI_Group grp;