diff --git a/MeshLib/IO/XDMF/HdfData.h b/MeshLib/IO/XDMF/HdfData.h
index 275f97b97edd79ede5407343396cbfb2cbf88f6c..50f6eb17a5be138b362bf03d8b5ddac378b0c0a5 100644
--- a/MeshLib/IO/XDMF/HdfData.h
+++ b/MeshLib/IO/XDMF/HdfData.h
@@ -21,7 +21,7 @@ namespace MeshLib::IO
 {
 using Hdf5DimType = unsigned long long;
 
-struct HdfData
+struct HdfData final
 {
     HdfData(void const* data_start, std::size_t size_partitioned_dim,
             std::size_t size_tuple, std::string const& name,
@@ -31,7 +31,7 @@ struct HdfData
     std::vector<Hdf5DimType> offsets;
     std::vector<Hdf5DimType> file_space;
     std::vector<Hdf5DimType> chunk_space;
-    std::string const name;
+    std::string name;
     int64_t data_type;
 };
 
diff --git a/MeshLib/IO/XDMF/HdfWriter.cpp b/MeshLib/IO/XDMF/HdfWriter.cpp
index fcdd27f465740088d11b9429dec2ccee3f61099b..ce28e48bace9a0d4df1a38dd888919bdbd769db4 100644
--- a/MeshLib/IO/XDMF/HdfWriter.cpp
+++ b/MeshLib/IO/XDMF/HdfWriter.cpp
@@ -71,23 +71,23 @@ static hid_t createDataSet(
 {
     int const time_dim_local_size = data_dims.size() + 1;
 
-    std::vector<Hdf5DimType> time_max_dims =
+    std::vector<Hdf5DimType> const time_max_dims =
         prependDimension(H5S_UNLIMITED, max_dims);
-    std::vector<Hdf5DimType> time_data_global_dims =
+    std::vector<Hdf5DimType> const time_data_global_dims =
         prependDimension(1, max_dims);
 
-    std::vector<Hdf5DimType> time_data_chunk_dims =
+    std::vector<Hdf5DimType> const time_data_chunk_dims =
         prependDimension(1, chunk_dims);
 
-    hid_t fspace =
+    hid_t const fspace =
         H5Screate_simple(time_dim_local_size, time_data_global_dims.data(),
                          time_max_dims.data());
     assert(fspace >= 0);
 
-    hid_t dcpl = H5Pcreate(H5P_DATASET_CREATE);
+    hid_t const dcpl = H5Pcreate(H5P_DATASET_CREATE);
     assert(dcpl >= 0);
 
-    hid_t status =
+    hid_t const status =
         H5Pset_chunk(dcpl, chunk_dims.size() + 1, time_data_chunk_dims.data());
     if (status < 0)
     {
@@ -99,8 +99,8 @@ static hid_t createDataSet(
         H5Pset_deflate(dcpl, default_compression_factor);
     }
 
-    hid_t dataset = H5Dcreate2(section, dataset_name.c_str(), data_type, fspace,
-                               H5P_DEFAULT, dcpl, H5P_DEFAULT);
+    hid_t const dataset = H5Dcreate2(section, dataset_name.c_str(), data_type,
+                                     fspace, H5P_DEFAULT, dcpl, H5P_DEFAULT);
 
     assert(dataset >= 0);
     H5Pclose(dcpl);
@@ -121,17 +121,18 @@ static void writeDataSet(
     std::vector<Hdf5DimType> const& offset_dims,
     std::vector<Hdf5DimType> const& max_dims,
     [[maybe_unused]] std::vector<Hdf5DimType> const& chunk_dims,
-    std::string const& dataset_name, int step, hid_t dataset)  // where
+    std::string const& dataset_name, Hdf5DimType const step,
+    hid_t const dataset)  // where
 {
-    Hdf5DimType hdf_step = step;
-    Hdf5DimType time_steps = hdf_step + 1;
+    Hdf5DimType const time_steps = step + 1;
 
-    std::vector<Hdf5DimType> time_data_local_dims = data_dims;
-    std::vector<Hdf5DimType> time_max_dims =
+    std::vector<Hdf5DimType> const time_data_local_dims = data_dims;
+    std::vector<Hdf5DimType> const time_max_dims =
         prependDimension(time_steps, max_dims);
-    std::vector<Hdf5DimType> time_offsets =
-        prependDimension(hdf_step, offset_dims);
-    std::vector<hsize_t> count = prependDimension(1, time_data_local_dims);
+    std::vector<Hdf5DimType> const time_offsets =
+        prependDimension(step, offset_dims);
+    std::vector<hsize_t> const count =
+        prependDimension(1, time_data_local_dims);
 
     hid_t const io_transfer_property = createHDF5TransferPolicy();
 
@@ -144,7 +145,7 @@ static void writeDataSet(
     {
         OGS_FATAL("H5D set extent failed dataset '{:s}'.", dataset_name);
     }
-    hid_t fspace = H5Dget_space(dataset);
+    hid_t const fspace = H5Dget_space(dataset);
 
     H5Sselect_hyperslab(fspace, H5S_SELECT_SET, time_offsets.data(), NULL,
                         count.data(), NULL);
@@ -161,76 +162,142 @@ static void writeDataSet(
 
     return;
 }
+
+/**
+ * \brief Write vector with time values into open hdf file
+ * \details In contrast to all other hdf write methods writing is only performed
+ * by one process (is_file_manager_true). file handle is to an already opened
+ * file
+ */
+static void writeTimeSeries(hid_t const file,
+                            std::vector<double> const& step_times,
+                            bool const is_file_manager)
+{
+    hsize_t const size = step_times.size();
+    hid_t const memspace = H5Screate_simple(1, &size, NULL);
+    hid_t const filespace = H5Screate_simple(1, &size, NULL);
+
+    if (is_file_manager)
+    {
+        H5Sselect_all(memspace);
+        H5Sselect_all(filespace);
+    }
+    else
+    {
+        H5Sselect_none(memspace);
+        H5Sselect_none(filespace);
+    }
+
+    hid_t const dataset =
+        H5Dcreate2(file, "/times", H5T_NATIVE_DOUBLE, filespace, H5P_DEFAULT,
+                   H5P_DEFAULT, H5P_DEFAULT);
+
+    H5Dwrite(dataset, H5T_NATIVE_DOUBLE, memspace, filespace, H5P_DEFAULT,
+             step_times.data());
+
+    H5Dclose(dataset);
+    H5Sclose(memspace);
+    H5Sclose(filespace);
+};
 namespace MeshLib::IO
 {
-HdfWriter::HdfWriter(std::vector<HdfData> constant_attributes,
-                     std::vector<HdfData>
-                         variable_attributes,
-                     int const initial_step,
+struct HdfWriter::HdfMesh final
+{
+    hid_t const group;
+    std::string const name;
+    std::map<std::string, hid_t> const datasets;
+    std::vector<HdfData> const variable_attributes;
+};
+
+HdfWriter::HdfWriter(std::vector<MeshHdfData> meshes,
+                     unsigned long long const initial_step,
                      std::filesystem::path const& filepath,
-                     bool const use_compression)
-    : _variable_attributes(std::move(variable_attributes)),
-      _hdf5_filepath(filepath),
-      _use_compression(checkCompression() && use_compression),
+                     bool const use_compression,
+                     bool const is_file_manager)
+    : _hdf5_filepath(filepath),
       _file(createFile(filepath)),
-      _output_step(initial_step)
+      _meshes_group(
+          H5Gcreate2(_file, "/meshes", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT)),
+      _step_times{0},  // ToDo need to be initial time
+      _use_compression(checkCompression() && use_compression),
+      _is_file_manager(is_file_manager)
 {
-    _group = H5Gcreate2(_file, "data", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
-
-    auto createAndWriteDataSet = [&](auto const& attribute) -> hid_t {
-        hid_t dataset = createDataSet(
-            attribute.data_type, attribute.data_space, attribute.file_space,
-            attribute.chunk_space, _use_compression, _group, attribute.name);
-
-        checkHdfStatus(dataset, "Creating HDF5 Dataset: {:s} failed.",
-                       attribute.name);
-        writeDataSet(attribute.data_start, attribute.data_type,
-                     attribute.data_space, attribute.offsets,
-                     attribute.file_space, attribute.chunk_space,
-                     attribute.name, _output_step, dataset);
-        return dataset;
-    };
-
-    for (auto const& attribute : constant_attributes)
+    for (auto const& mesh : meshes)
     {
-        hid_t dataset = createAndWriteDataSet(attribute);
-        H5Dclose(dataset);
-    }
+        hid_t const group = H5Gcreate2(_meshes_group, mesh.name.c_str(),
+                                       H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
 
-    for (auto const& attribute : _variable_attributes)
-    {
-        hid_t dataset = createAndWriteDataSet(attribute);
-        // datasets are kept open
-        _datasets.insert({attribute.name, dataset});
+        auto const createAndWriteDataSet = [&](auto const& attribute) -> hid_t
+        {
+            hid_t const dataset = createDataSet(
+                attribute.data_type, attribute.data_space, attribute.file_space,
+                attribute.chunk_space, _use_compression, group, attribute.name);
+
+            checkHdfStatus(dataset, "Creating HDF5 Dataset: {:s} failed.",
+                           attribute.name);
+
+            writeDataSet(attribute.data_start, attribute.data_type,
+                         attribute.data_space, attribute.offsets,
+                         attribute.file_space, attribute.chunk_space,
+                         attribute.name, initial_step, dataset);
+            return dataset;
+        };
+
+        for (auto const& attribute : mesh.constant_attributes)
+        {
+            hid_t const dataset = createAndWriteDataSet(attribute);
+            H5Dclose(dataset);
+        }
+
+        std::map<std::string, hid_t> datasets;
+        for (auto const& attribute : mesh.variable_attributes)
+        {
+            hid_t const dataset = createAndWriteDataSet(attribute);
+            // datasets are kept open
+            datasets.insert({attribute.name, dataset});
+        }
+
+        _hdf_meshes.push_back(std::make_unique<HdfMesh>(
+            HdfMesh{group, mesh.name, datasets, mesh.variable_attributes}));
     }
-    _output_step++;
 }
 
 HdfWriter::~HdfWriter()
 {
-    for (auto& dataset : _datasets)
+    writeTimeSeries(_file, _step_times, _is_file_manager);
+
+    for (auto const& mesh : _hdf_meshes)
     {
-        H5Dclose(dataset.second);
+        for (auto const& dataset : mesh->datasets)
+        {
+            H5Dclose(dataset.second);
+        }
+        H5Gclose(mesh->group);
     }
-    H5Gclose(_group);
+    H5Gclose(_meshes_group);
     H5Fclose(_file);
 }
 
-void HdfWriter::writeStep()
+void HdfWriter::writeStep(double const time)
 {
-    for (auto const& attribute : _variable_attributes)
+    _step_times.push_back(time);
+    auto const output_step = _step_times.size();
+
+    for (auto const& mesh : _hdf_meshes)
     {
-        auto const& dataset_hid = _datasets.find(attribute.name);
-        if (dataset_hid == _datasets.end())
+        for (auto const& attribute : mesh->variable_attributes)
         {
-            OGS_FATAL("Writing HDF5 Dataset: {:s} failed.", attribute.name);
+            auto const& dataset_hid = mesh->datasets.find(attribute.name);
+            if (dataset_hid == mesh->datasets.end())
+            {
+                OGS_FATAL("Writing HDF5 Dataset: {:s} failed.", attribute.name);
+            }
+
+            writeDataSet(
+                attribute.data_start, attribute.data_type, attribute.data_space,
+                attribute.offsets, attribute.file_space, attribute.chunk_space,
+                attribute.name, output_step, mesh->datasets.at(attribute.name));
         }
-
-        writeDataSet(
-            attribute.data_start, attribute.data_type, attribute.data_space,
-            attribute.offsets, attribute.file_space, attribute.chunk_space,
-            attribute.name, _output_step, _datasets.at(attribute.name));
     }
-    _output_step++;
 }
 }  // namespace MeshLib::IO
diff --git a/MeshLib/IO/XDMF/HdfWriter.h b/MeshLib/IO/XDMF/HdfWriter.h
index 23edc9a50ccfa8a7a8d3f9c59add8a5bf0c60fe7..97c8707850a552134f3efd4092502089933e88ee 100644
--- a/MeshLib/IO/XDMF/HdfWriter.h
+++ b/MeshLib/IO/XDMF/HdfWriter.h
@@ -15,12 +15,21 @@
 #include <hdf5.h>
 
 #include <map>
+#include <memory>
 #include <vector>
 
 #include "HdfData.h"
 
 namespace MeshLib::IO
 {
+using HDFAttributes = std::vector<HdfData>;
+struct MeshHdfData
+{
+    HDFAttributes constant_attributes;
+    HDFAttributes variable_attributes;
+    std::string name;
+};
+
 class HdfWriter final
 {
 public:
@@ -28,40 +37,39 @@ public:
      * \brief Write file with geometry and topology data. The data
      * itself is held by a structure outside of this class. The writer assumes
      * the data holder to not change during writing
-     * @param constant_attributes vector of constant attributes (each attribute
-     * is a OGS mesh property), geometry and topology are considered as constant
-     * attributes
-     * @param variable_attributes vector of variable attributes (each attribute
-     * is a OGS mesh property
+     * @param meshes meta data of meshes to be written
      * @param initial_step number of the step (temporal collection), usually 0,
      * greater 0 with continuation of simulation
      * @param filepath absolute or relative filepath to the hdf5 file
      * @param use_compression if true gzip compression is enabled
+     * @param is_file_manager True if process (in parallel execution) is
+     * File_Manager
      */
-    HdfWriter(std::vector<HdfData> constant_attributes,
-              std::vector<HdfData>
-                  variable_attributes,
-              int const initial_step,
+    HdfWriter(std::vector<MeshHdfData> meshes,
+              unsigned long long initial_step,
               std::filesystem::path const& filepath,
-              bool const use_compression);
-
+              bool use_compression,
+              bool is_file_manager);
     /**
      * \brief Writes attributes. The data
      * itself is hold by a structure outside of this class. The writer assumes
      * the data holder to not change during writing and HdfData given to
      * constructor to be still valid
-     * @param step number of the step (temporal collection)
+     * @param time time_value of step to be written to temporal collection
      */
-    void writeStep();
+    void writeStep(double time);
     ~HdfWriter();
 
 private:
-    std::vector<HdfData> const _variable_attributes;
+    // internal data holder
+    struct HdfMesh;
+
     std::filesystem::path const _hdf5_filepath;
-    bool const _use_compression;
     hid_t const _file;
-    hid_t _group;
-    std::map<std::string, hid_t> _datasets;
-    int _output_step;
+    hid_t const _meshes_group;
+    std::vector<std::unique_ptr<HdfMesh>> _hdf_meshes;
+    std::vector<double> _step_times;
+    bool const _use_compression;
+    bool const _is_file_manager;
 };
 }  // namespace MeshLib::IO
\ No newline at end of file
diff --git a/MeshLib/IO/XDMF/XdmfData.h b/MeshLib/IO/XDMF/XdmfData.h
index b6a826eeed8c5cd00ab4514a8463fca5440817eb..5d256f248120b5ec1c44d2f5839b0e824126164e 100644
--- a/MeshLib/IO/XDMF/XdmfData.h
+++ b/MeshLib/IO/XDMF/XdmfData.h
@@ -55,12 +55,12 @@ struct XdmfData final
              unsigned int const index);
     // a hyperslab is defined by starts and strides see
     // https://www.xdmf.org/index.php/XDMF_Model_and_Format#HyperSlab
-    std::vector<XdmfDimType> const starts;
-    std::vector<XdmfDimType> const strides;
+    std::vector<XdmfDimType> starts;
+    std::vector<XdmfDimType> strides;
     std::vector<XdmfDimType> global_block_dims;
-    MeshPropertyDataType const data_type;
-    std::string const name;
-    std::optional<MeshLib::MeshItemType> const attribute_center;
+    MeshPropertyDataType data_type;
+    std::string name;
+    std::optional<MeshLib::MeshItemType> attribute_center;
     unsigned int index;
 };
 
diff --git a/MeshLib/IO/XDMF/XdmfHdfData.h b/MeshLib/IO/XDMF/XdmfHdfData.h
index b56560fc7d39770a483f03d6af3f00ae8c864c9a..992900bd7f73383067bdb2fd81a34fd623d8c60a 100644
--- a/MeshLib/IO/XDMF/XdmfHdfData.h
+++ b/MeshLib/IO/XDMF/XdmfHdfData.h
@@ -18,25 +18,9 @@
 #include "HdfData.h"
 #include "XdmfData.h"
 
-// TODO (tm) This structs are the result of transformData.cpp. This structures
-// can be eliminated.
 namespace MeshLib::IO
 {
-struct Geometry final
-{
-    std::vector<double> flattened_values;
-    HdfData hdf;
-    XdmfData xdmf;
-};
-
-struct Topology final
-{
-    std::vector<int> flattened_values;
-    HdfData hdf;
-    XdmfData xdmf;
-};
-
-struct AttributeMeta final
+struct XdmfHdfData final
 {
     HdfData hdf;
     XdmfData xdmf;
diff --git a/MeshLib/IO/XDMF/XdmfHdfWriter.cpp b/MeshLib/IO/XDMF/XdmfHdfWriter.cpp
index 9230a80ad6e6cba08ba25581804e442adaf81591..563a80c1144a714b3c9f3279a308874283b2335d 100644
--- a/MeshLib/IO/XDMF/XdmfHdfWriter.cpp
+++ b/MeshLib/IO/XDMF/XdmfHdfWriter.cpp
@@ -19,34 +19,43 @@
 
 namespace MeshLib::IO
 {
-XdmfHdfWriter::XdmfHdfWriter(MeshLib::Mesh const& mesh,
-                             std::filesystem::path const& filepath,
-                             int const time_step, double const initial_time,
-                             std::set<std::string> const& variable_output_names,
-                             bool const use_compression)
+struct TransformedMeshData final
 {
-    // transform Data into contiguous data and collect meta data
-    Geometry const& geometry = transformGeometry(mesh);
-    Topology const& topology = transformTopology(mesh, geometry.hdf.offsets[0]);
-
-    std::vector<AttributeMeta> attributes = transformAttributes(mesh);
-    std::vector<HdfData> hdf_data_attributes = {geometry.hdf, topology.hdf};
-    std::transform(attributes.begin(), attributes.end(),
-                   std::back_inserter(hdf_data_attributes),
-                   [](AttributeMeta att) -> HdfData { return att.hdf; });
-
-    // if no output name is specified, all data will be assumened to be variable
-    // over the timesteps. The xdmfhdfwriter is an alternative to other writers,
-    // that do not consider the constantness of data Callers of xdmfwriter (e.g.
-    // ogs tools) do not provide these information yet and indicate with empty
-    // list
-    std::function<bool(HdfData)> is_variable_hdf_attribute =
+    std::vector<double> flattened_geometry_values;
+    std::vector<int> flattened_topology_values;
+};
+struct XdmfHdfMesh final
+{
+    XdmfHdfData geometry;
+    XdmfHdfData topology;
+    std::vector<XdmfHdfData> attributes;
+    std::string name;
+    // TransformedMeshData may be large, ensure it is never copied
+    std::unique_ptr<TransformedMeshData> transformed_data;
+};
+
+XdmfHdfWriter::XdmfHdfWriter(
+    std::vector<std::reference_wrapper<const MeshLib::Mesh>> meshes,
+    std::filesystem::path const& filepath, unsigned long long const time_step,
+    double const initial_time,
+    std::set<std::string> const& variable_output_names,
+    bool const use_compression)
+{
+    // ogs meshes to vector of Xdmf/HDF meshes (we keep Xdmf and HDF together
+    // because XDMF depends on HDF) to meta
+
+    // if no output name is specified, all data will be assumened to be
+    // variable over the timesteps. The xdmfhdfwriter is an alternative to
+    // other writers, that do not consider the constantness of data Callers
+    // of xdmfwriter (e.g. ogs tools) do not provide these information yet
+    // and indicate with empty list
+    std::function<bool(HdfData)> const is_variable_hdf_attribute =
         [&variable_output_names](
             bool outputnames) -> std::function<bool(HdfData)>
     {
         if (outputnames)
         {
-            return [&variable_output_names](HdfData data) -> bool
+            return [&variable_output_names](HdfData const& data) -> bool
             {
                 return std::find(variable_output_names.begin(),
                                  variable_output_names.end(),
@@ -55,89 +64,167 @@ XdmfHdfWriter::XdmfHdfWriter(MeshLib::Mesh const& mesh,
         }
         else
         {
-            return [](HdfData) -> bool { return true; };
+            return [](HdfData const&) -> bool { return true; };
         }
     }(!variable_output_names.empty());
 
-    std::vector<HdfData> hdf_variable_attributes;
-    std::copy_if(hdf_data_attributes.begin(), hdf_data_attributes.end(),
-                 back_inserter(hdf_variable_attributes),
-                 is_variable_hdf_attribute);
-    std::vector<HdfData> hdf_constant_attributes;
-    std::copy_if(hdf_data_attributes.begin(), hdf_data_attributes.end(),
-                 back_inserter(hdf_constant_attributes),
-                 std::not_fn(is_variable_hdf_attribute));
+    auto is_variable_xdmf_attribute =
+        [&variable_output_names](XdmfData const& data) -> bool
+    {
+        return std::find(variable_output_names.begin(),
+                         variable_output_names.end(),
+                         data.name) != variable_output_names.end();
+    };
+
+    // Transform the data to be written into a format conforming with the rules
+    // of xdmf topology and geometry
+    auto const transform_ogs_mesh_data_to_xdmf_conforming_data =
+        [](auto const& mesh)
+    {
+        auto flattened_geometry_values = transformToXDMFGeometry(mesh);
+        // actually this line is only needed to calculate the offset
+        XdmfHdfData const& geometry =
+            transformGeometry(mesh, flattened_geometry_values.data());
+        auto const flattened_topology_values =
+            transformToXDMFTopology(mesh, geometry.hdf.offsets[0]);
+        return std::make_unique<TransformedMeshData>(
+            TransformedMeshData(std::move(flattened_geometry_values),
+                                std::move(flattened_topology_values)));
+    };
+
+    // create metadata for transformed data and original ogs mesh data
+    auto const transform_to_meta_data =
+        [&transform_ogs_mesh_data_to_xdmf_conforming_data](auto const& mesh)
+    {
+        // important: transformed data must survive and be unique, raw pointer
+        // to its memory!
+        auto xdmf_conforming_data =
+            transform_ogs_mesh_data_to_xdmf_conforming_data(mesh);
+        auto const geometry = transformGeometry(
+            mesh, xdmf_conforming_data->flattened_geometry_values.data());
+        auto const topology =
+            transformTopology(xdmf_conforming_data->flattened_topology_values);
+        auto const attributes = transformAttributes(mesh);
+        return XdmfHdfMesh{std::move(geometry), std::move(topology),
+                           std::move(attributes), mesh.get().getName(),
+                           std::move(xdmf_conforming_data)};
+    };
+
+    // extract meta data relevant for HDFWriter
+    auto const transform_metamesh_to_hdf =
+        [&is_variable_hdf_attribute](auto const& metamesh)
+    {
+        // topology and geometry can be treated as any other attribute
+        std::vector<HdfData> hdf_data_attributes = {metamesh.geometry.hdf,
+                                                    metamesh.topology.hdf};
+
+        hdf_data_attributes.reserve(hdf_data_attributes.size() +
+                                    metamesh.attributes.size());
+        std::transform(metamesh.attributes.begin(), metamesh.attributes.end(),
+                       std::back_inserter(hdf_data_attributes),
+                       [](XdmfHdfData att) -> HdfData { return att.hdf; });
+
+        HDFAttributes constant_attributes;
+        std::copy_if(hdf_data_attributes.begin(), hdf_data_attributes.end(),
+                     back_inserter(constant_attributes),
+                     std::not_fn(is_variable_hdf_attribute));
+        HDFAttributes variable_attributes;
+        std::copy_if(hdf_data_attributes.begin(), hdf_data_attributes.end(),
+                     back_inserter(variable_attributes),
+                     is_variable_hdf_attribute);
+
+        return MeshHdfData{
+            .constant_attributes = std::move(constant_attributes),
+            .variable_attributes = std::move(variable_attributes),
+            .name = std::move(metamesh.name)};
+    };
+
+    // --------------- XDMF + HDF ---------------------
+    std::vector<XdmfHdfMesh> xdmf_hdf_meshes;
+    xdmf_hdf_meshes.reserve(meshes.size());
+    std::transform(meshes.begin(), meshes.end(),
+                   std::back_inserter(xdmf_hdf_meshes), transform_to_meta_data);
 
-    // HDF5
+    std::vector<MeshHdfData> hdf_meshes;
+    hdf_meshes.reserve(xdmf_hdf_meshes.size());
+    std::transform(xdmf_hdf_meshes.begin(), xdmf_hdf_meshes.end(),
+                   std::back_inserter(hdf_meshes), transform_metamesh_to_hdf);
+
+    // --------------- HDF ---------------------
     std::filesystem::path const hdf_filepath =
         filepath.parent_path() / (filepath.stem().string() + ".h5");
 
-    // geometry hdf data refers to the local data geometry and topology vector.
-    // The hdf writer can write when data is out of scope.
-    _hdf_writer = std::make_unique<HdfWriter>(
-        std::move(hdf_constant_attributes), std::move(hdf_variable_attributes),
-        time_step, hdf_filepath, use_compression);
-    // XDMF
+    auto const is_file_manager = isFileManager();
+    _hdf_writer = std::make_unique<HdfWriter>(std::move(hdf_meshes), time_step,
+                                              hdf_filepath, use_compression,
+                                              is_file_manager);
+
+    // --------------- XDMF ---------------------
     // The light data is only written by just one process
-    if (!isFileManager())
+    if (!is_file_manager)
     {
         return;
     }
 
-    std::filesystem::path const xdmf_filepath =
-        filepath.parent_path() / (filepath.stem().string() + ".xdmf");
+    // xdmf section
+    // extract meta data relevant for XDMFWriter
+    auto const transform_metamesh_to_xdmf =
+        [&is_variable_xdmf_attribute, &filepath, &hdf_filepath,
+         &initial_time](XdmfHdfMesh& metamesh)
+    {
+        std::string const xdmf_name = metamesh.name;
+        std::filesystem::path const xdmf_filepath =
+            filepath.parent_path() / (xdmf_name + ".xdmf");
 
-    std::vector<XdmfData> xdmf_attributes;
-    std::transform(attributes.begin(), attributes.end(),
-                   std::back_inserter(xdmf_attributes),
-                   [](AttributeMeta const& att) -> XdmfData
-                   { return att.xdmf; });
+        std::vector<XdmfData> xdmf_attributes;
+        std::transform(metamesh.attributes.begin(), metamesh.attributes.end(),
+                       std::back_inserter(xdmf_attributes),
+                       [](XdmfHdfData const& att) -> XdmfData
+                       { return att.xdmf; });
 
-    for (size_t i = 0; i < attributes.size(); ++i)
-    {
-        // index 1 time,  index 2 geo, index 3 topology, attributes start at
-        // index 4
-        xdmf_attributes[i].index = i + 4;
-    }
+        for (std::size_t i = 0; i < metamesh.attributes.size(); ++i)
+        {
+            // index 1 time,  index 2 geo, index 3 topology, attributes start at
+            // index 4
+            xdmf_attributes[i].index = i + 4;
+        }
 
-    std::function<bool(XdmfData)> is_variable_xdmf_attribute =
-        [&variable_output_names](XdmfData data) -> bool
-    {
-        return std::find(variable_output_names.begin(),
-                         variable_output_names.end(),
-                         data.name) != variable_output_names.end();
+        std::vector<XdmfData> xdmf_variable_attributes;
+        std::copy_if(xdmf_attributes.begin(), xdmf_attributes.end(),
+                     back_inserter(xdmf_variable_attributes),
+                     is_variable_xdmf_attribute);
+        std::vector<XdmfData> xdmf_constant_attributes;
+        std::copy_if(xdmf_attributes.begin(), xdmf_attributes.end(),
+                     back_inserter(xdmf_constant_attributes),
+                     std::not_fn(is_variable_xdmf_attribute));
+
+        auto const xdmf_writer_fn =
+            write_xdmf(metamesh.geometry.xdmf, metamesh.topology.xdmf,
+                       xdmf_constant_attributes, xdmf_variable_attributes,
+                       hdf_filepath.filename().string(),
+                       GitInfoLib::GitInfo::ogs_version, xdmf_name);
+        auto xdmf_writer = std::make_unique<XdmfWriter>(xdmf_filepath.string(),
+                                                        xdmf_writer_fn);
+        xdmf_writer->addTimeStep(initial_time);
+        return xdmf_writer;
     };
 
-    std::vector<XdmfData> xdmf_variable_attributes;
-    std::copy_if(xdmf_attributes.begin(), xdmf_attributes.end(),
-                 back_inserter(xdmf_variable_attributes),
-                 is_variable_xdmf_attribute);
-    std::vector<XdmfData> xdmf_constant_attributes;
-    std::copy_if(xdmf_attributes.begin(), xdmf_attributes.end(),
-                 back_inserter(xdmf_constant_attributes),
-                 std::not_fn(is_variable_xdmf_attribute));
-
-    if (isFileManager())
-    {
-        auto xdmf_writer_fn = write_xdmf(
-            geometry.xdmf, topology.xdmf, xdmf_constant_attributes,
-            xdmf_variable_attributes, hdf_filepath.filename().string(),
-            GitInfoLib::GitInfo::ogs_version);
-        _xdmf_writer = std::make_unique<XdmfWriter>(xdmf_filepath.string(),
-                                                    xdmf_writer_fn);
-        _xdmf_writer->addTimeStep(initial_time);
-    }
+    std::transform(xdmf_hdf_meshes.begin(), xdmf_hdf_meshes.end(),
+                   std::back_inserter(_xdmf_writer),
+                   transform_metamesh_to_xdmf);
 }
 
-void XdmfHdfWriter::writeStep([[maybe_unused]] int const time_step,
-                              double const time)
+void XdmfHdfWriter::writeStep(double const time)
 {
     // ToDo (tm) time_step will be used for simulation continuation (restart)
-    _hdf_writer->writeStep();
+    _hdf_writer->writeStep(time);
     // The light data is only written by just one process
     if (isFileManager())
     {
-        _xdmf_writer->addTimeStep(time);
+        for (auto const& xdmf_writer : _xdmf_writer)
+        {
+            xdmf_writer->addTimeStep(time);
+        }
     }
 }
 
diff --git a/MeshLib/IO/XDMF/XdmfHdfWriter.h b/MeshLib/IO/XDMF/XdmfHdfWriter.h
index bf843dd501434fbbb83926ccc600ab264e7df236..69d8b8e233fbaff0327339a7ef4bf85300876506 100644
--- a/MeshLib/IO/XDMF/XdmfHdfWriter.h
+++ b/MeshLib/IO/XDMF/XdmfHdfWriter.h
@@ -28,7 +28,7 @@ class XdmfHdfWriter final
 public:
     /**
      * \brief Write xdmf and h5 file with geometry and topology data.
-     * @param mesh Mesh or NodePartitionedMesh to be written to file(s)
+     * @param meshes Meshes or NodePartitionedMeshes to be written to file(s)
      * @param filepath absolute or relative filepath to the hdf5 file
      * @param time_step number of the step (temporal collection)
      * @param initial_time time in seconds of the first time step
@@ -37,22 +37,21 @@ public:
      * @param use_compression if true, zlib compression in HDFWriter component
      * is used
      */
-    XdmfHdfWriter(MeshLib::Mesh const& mesh,
-                  std::filesystem::path const& filepath, int time_step,
-                  double initial_time,
-                  std::set<std::string> const& variable_output_names,
-                  bool use_compression);
+    XdmfHdfWriter(
+        std::vector<std::reference_wrapper<const MeshLib::Mesh>> meshes,
+        std::filesystem::path const& filepath, unsigned long long time_step,
+        double initial_time, std::set<std::string> const& variable_output_names,
+        bool use_compression);
 
     /**
      * \brief Adds data for either lazy (xdmf) or eager (hdf) writing algorithm
-     * @param time_step number of the step (temporal collection)
      * @param time time value of the current time_step
      */
-    void writeStep(int time_step, double time);
+    void writeStep(double time);
 
 private:
     // hdf_writer must be destructed before xdmf_writer
     std::unique_ptr<HdfWriter> _hdf_writer;
-    std::unique_ptr<XdmfWriter> _xdmf_writer;
+    std::vector<std::unique_ptr<XdmfWriter>> _xdmf_writer;
 };
 }  // namespace MeshLib::IO
\ No newline at end of file
diff --git a/MeshLib/IO/XDMF/XdmfWriter.cpp b/MeshLib/IO/XDMF/XdmfWriter.cpp
index 67ffb7cc66959e28de57f63bf0e850154b678799..e2657c5c2fbab3d0462d0b92767fb24c6fa7820e 100644
--- a/MeshLib/IO/XDMF/XdmfWriter.cpp
+++ b/MeshLib/IO/XDMF/XdmfWriter.cpp
@@ -20,10 +20,12 @@ XdmfWriter::~XdmfWriter()
     BaseLib::RunTime time_output;
     time_output.start();
     std::ofstream fout;
-    fout.open(this->filename);
+    fout.open(filename);
     fout << xdmf_writer(times);
 
-    INFO("[time] Output of XDMF took {:g} s.", time_output.elapsed());
+    INFO("[time] Output of XDMF to {:s} took {:g} s.",
+         filename,
+         time_output.elapsed());
 }
 
 void XdmfWriter::addTimeStep(double const& time_step)
diff --git a/MeshLib/IO/XDMF/transformData.cpp b/MeshLib/IO/XDMF/transformData.cpp
index 734b41c79a17eee24773f0836441e72983476aac..cd07b709a27ac7058db884c36d67e3d544e1ccf5 100644
--- a/MeshLib/IO/XDMF/transformData.cpp
+++ b/MeshLib/IO/XDMF/transformData.cpp
@@ -69,7 +69,7 @@ constexpr auto cellTypeOGS2XDMF(MeshLib::CellType const& cell_type)
     return elem_type_ogs2xdmf[to_underlying(cell_type)];
 }
 
-std::optional<AttributeMeta> transformAttribute(
+std::optional<XdmfHdfData> transformAttribute(
     std::pair<std::string, PropertyVectorBase*> const& property_pair)
 {
     // 3 data that will be captured and written by lambda f below
@@ -80,8 +80,8 @@ std::optional<AttributeMeta> transformAttribute(
     // lambda f : Collects properties from the propertyVectorBase. It captures
     // (and overwrites) data that can only be collected via the typed property.
     // It has boolean return type to allow kind of pipe using || operator.
-    auto f = [&data_type, &num_of_tuples, &data_ptr,
-              &property_pair](auto basic_type) -> bool
+    auto f = [&data_type, &num_of_tuples, &data_ptr, &property_pair](
+                 auto basic_type) -> bool
     {
         auto const property_base = property_pair.second;
         auto const typed_property =
@@ -190,30 +190,31 @@ std::optional<AttributeMeta> transformAttribute(
 
     std::string const& name = property_base->getPropertyName();
 
-    HdfData hdf =
-        HdfData(data_ptr, num_of_tuples, ui_global_components, name, data_type);
+    HdfData hdf = {data_ptr, num_of_tuples, ui_global_components, name,
+                   data_type};
 
-    XdmfData xdmf = XdmfData(num_of_tuples, ui_global_components, data_type,
-                             name, mesh_item_type, 0);
+    XdmfData xdmf = {num_of_tuples, ui_global_components, data_type,
+                     name,          mesh_item_type,       0};
 
-    return AttributeMeta{std::move(hdf), std::move(xdmf)};
+    return XdmfHdfData{std::move(hdf), std::move(xdmf)};
 }
 
-std::vector<AttributeMeta> transformAttributes(MeshLib::Mesh const& mesh)
+std::vector<XdmfHdfData> transformAttributes(MeshLib::Mesh const& mesh)
 {
     MeshLib::Properties const& properties = mesh.getProperties();
 
     // \TODO (tm) use c++20 ranges
     // a = p | filter (first!=OGS_VERSION) | filter null_opt | transformAttr |
-    std::vector<AttributeMeta> attributes;
-    for (auto [name, property_base] : properties)
+    std::vector<XdmfHdfData> attributes;
+    for (auto const& [name, property_base] : properties)
     {
         if (name == GitInfoLib::GitInfo::OGS_VERSION)
         {
             continue;
         }
 
-        if (auto attribute = transformAttribute(std::pair(name, property_base)))
+        if (auto const attribute =
+                transformAttribute(std::pair(name, property_base)))
         {
             attributes.push_back(attribute.value());
         }
@@ -225,9 +226,8 @@ std::vector<AttributeMeta> transformAttributes(MeshLib::Mesh const& mesh)
     return attributes;
 }
 
-Geometry transformGeometry(MeshLib::Mesh const& mesh)
+std::vector<double> transformToXDMFGeometry(MeshLib::Mesh const& mesh)
 {
-    std::string const name = "geometry";
     std::vector<MeshLib::Node*> const& nodes = mesh.getNodes();
 
     int const point_size = 3;
@@ -239,23 +239,29 @@ Geometry transformGeometry(MeshLib::Mesh const& mesh)
         values.insert(values.cend(), x, x + point_size);
     }
 
+    return values;
+}
+
+XdmfHdfData transformGeometry(MeshLib::Mesh const& mesh, double const* data_ptr)
+{
+    std::string const name = "geometry";
+    std::vector<MeshLib::Node*> const& nodes = mesh.getNodes();
+
+    int const point_size = 3;
     auto const& partition_dim = nodes.size();
 
-    HdfData hdf = HdfData(values.data(),
-                          partition_dim,
-                          point_size,
-                          name,
-                          MeshPropertyDataType::float64);
-    XdmfData xdmf =
-        XdmfData(partition_dim, point_size, MeshPropertyDataType::float64, name,
-                 std::nullopt, 2);
+    HdfData const hdf = {data_ptr, partition_dim, point_size, name,
+                         MeshPropertyDataType::float64};
+    XdmfData const xdmf = {
+        partition_dim, point_size,   MeshPropertyDataType::float64,
+        name,          std::nullopt, 2};
 
-    return Geometry{std::move(values), std::move(hdf), std::move(xdmf)};
+    return XdmfHdfData{std::move(hdf), std::move(xdmf)};
 }
 
-Topology transformTopology(MeshLib::Mesh const& mesh, std::size_t const offset)
+std::vector<int> transformToXDMFTopology(MeshLib::Mesh const& mesh,
+                                         std::size_t const offset)
 {
-    std::string const name = "topology";
     std::vector<MeshLib::Element*> const& elements = mesh.getElements();
     std::vector<int> values;
     values.reserve(elements.size());
@@ -277,12 +283,17 @@ Topology transformTopology(MeshLib::Mesh const& mesh, std::size_t const offset)
             values.push_back(node->getID() + offset);
         }
     }
+    return values;
+}
 
-    HdfData hdf = HdfData(values.data(), values.size(), 1, name,
-                          MeshPropertyDataType::int32);
-    XdmfData xdmf = XdmfData(values.size(), 1, MeshPropertyDataType::int32,
-                             name, std::nullopt, 3);
+XdmfHdfData transformTopology(std::vector<int> const& values)
+{
+    std::string const name = "topology";
+    HdfData const hdf = {values.data(), values.size(), 1, name,
+                         MeshPropertyDataType::int32};
+    XdmfData const xdmf = {values.size(), 1, MeshPropertyDataType::int32, name,
+                           std::nullopt,  3};
 
-    return Topology{std::move(values), std::move(hdf), std::move(xdmf)};
+    return XdmfHdfData{std::move(hdf), std::move(xdmf)};
 }
 }  // namespace MeshLib::IO
diff --git a/MeshLib/IO/XDMF/transformData.h b/MeshLib/IO/XDMF/transformData.h
index 367615f7a715999ec320bd7a21a41f7f51d10d70..16284421061aae9e1222594528f55b1344a47e0f 100644
--- a/MeshLib/IO/XDMF/transformData.h
+++ b/MeshLib/IO/XDMF/transformData.h
@@ -28,14 +28,29 @@ namespace MeshLib::IO
  * @param mesh OGS mesh can be mesh or partitionedMesh
  * @return vector of meta data
  */
-std::vector<AttributeMeta> transformAttributes(MeshLib::Mesh const& mesh);
+std::vector<XdmfHdfData> transformAttributes(MeshLib::Mesh const& mesh);
+/**
+ * \brief Create meta data for geometry used for hdf5 and xdmf
+ * @param mesh OGS mesh can be mesh or partitionedMesh
+ * @param data_ptr Memory location of geometry values.
+ * @return Geometry meta data
+ */
+XdmfHdfData transformGeometry(MeshLib::Mesh const& mesh,
+                              double const* data_ptr);
+/**
+ * \brief  Create meta data for topology used for HDF5 and XDMF
+ * @param values actual topology values to get size and memory location
+ * @return Topology meta data
+ */
+XdmfHdfData transformTopology(std::vector<int> const& values);
 /**
  * \brief Copies all node points into a new vector. Contiguous data used for
- * writing.
+ * writing. Conform with XDMF standard in
+ * https://www.xdmf.org/index.php/XDMF_Model_and_Format
  * @param mesh OGS mesh can be mesh or partitionedMesh
- * @return Geometry containing a copy of the data and the geometry meta data
+ * @return vector containing a copy of the data
  */
-Geometry transformGeometry(MeshLib::Mesh const& mesh);
+std::vector<double> transformToXDMFGeometry(MeshLib::Mesh const& mesh);
 /**
  * \brief Copies all cells into a new vector. Contiguous data used for writing.
  * The topology is specific to xdmf because it contains the xdmf cell types!!
@@ -43,8 +58,8 @@ Geometry transformGeometry(MeshLib::Mesh const& mesh);
  * @param mesh OGS mesh can be mesh or partitionedMesh
  * @param offset Local offset to transform local to global cell ID. Offset must
  * be zero in serial and must be defined for each process in parallel execution.
-
- * @return Topology containing a copy of the data and the topology meta data
+ * @return vector containing a copy of the data
  */
-Topology transformTopology(MeshLib::Mesh const& mesh, std::size_t offset);
+std::vector<int> transformToXDMFTopology(MeshLib::Mesh const& mesh,
+                                         std::size_t const offset);
 }  // namespace MeshLib::IO
diff --git a/MeshLib/IO/XDMF/writeXdmf.cpp b/MeshLib/IO/XDMF/writeXdmf.cpp
index 5829ba8483c143c98211629e95afb8f672087cc0..50742b7124dd07b07e395a455acf13eaf368dfac 100644
--- a/MeshLib/IO/XDMF/writeXdmf.cpp
+++ b/MeshLib/IO/XDMF/writeXdmf.cpp
@@ -86,13 +86,13 @@ static constexpr auto meshPropertyDatatypeSize()
     return property_sizes;
 }
 
-inline auto ogs_to_xdmf_type = meshPropertyDatatypeSize();
+inline auto ogs_to_xdmf_type_fn = meshPropertyDatatypeSize();
 
 // Transform MeshPropertyDataType into type_sizes (in bytes)
 static std::string getPropertyDataTypeSize(
     MeshPropertyDataType const& ogs_data_type)
 {
-    return fmt::format("{}", ogs_to_xdmf_type[to_underlying(ogs_data_type)]);
+    return fmt::format("{}", ogs_to_xdmf_type_fn[to_underlying(ogs_data_type)]);
 }
 
 // Collects all known data and return a function that takes all unknown data
@@ -101,22 +101,26 @@ std::function<std::string(std::vector<double>)> write_xdmf(
     XdmfData const& geometry, XdmfData const& topology,
     std::vector<XdmfData> const& constant_attributes,
     std::vector<XdmfData> const& variable_attributes,
-    std::string const& h5filename, std::string const& ogs_version)
+    std::string const& h5filename, std::string const& ogs_version,
+    std::string const& mesh_name)
 {
     // Generates function that writes <DataItem>. Late binding needed because
     // maximum number of steps unknown. Time step and h5filename are captured to
     // return a unary function
     // _a suffix is a user defined literal for fmt named arguments
-    auto const time_dataitem_genfn = [](int const time_step, int const max_step,
-                                        std::string const h5filename)
+    auto const time_dataitem_genfn =
+        [](unsigned long long const time_step, int const max_step,
+           std::string const& h5filename, std::string const& mesh_name)
     {
-        return [time_step, max_step, h5filename](auto const& xdmfdata)
+        return
+            [time_step, max_step, h5filename, mesh_name](auto const& xdmfdata)
         {
             return fmt::format(
                 "\n\t\t<DataItem DataType=\"{datatype}\" "
                 "Dimensions=\"{local_dimensions}\" "
                 "Format=\"HDF\" "
-                "Precision=\"{precision}\">{filename}:data/{datasetname}|"
+                "Precision=\"{precision}\">"
+                "{filename}:meshes/{meshname}/{datasetname}|"
                 "{time_step} {starts}:1 {strides}:1 "
                 "{local_dimensions}:{max_step} "
                 "{global_dimensions}</"
@@ -126,6 +130,7 @@ std::function<std::string(std::vector<double>)> write_xdmf(
                     fmt::join(xdmfdata.global_block_dims, " "),
                 "precision"_a = getPropertyDataTypeSize(xdmfdata.data_type),
                 "filename"_a = h5filename,
+                "meshname"_a = mesh_name,
                 "datasetname"_a = xdmfdata.name,
                 "starts"_a = fmt::join(xdmfdata.starts, " "),
                 "strides"_a = fmt::join(xdmfdata.strides, " "),
@@ -224,13 +229,15 @@ std::function<std::string(std::vector<double>)> write_xdmf(
     // Generates a function that either writes the data or points to existing
     // data
     auto const time_step_fn = [time_dataitem_genfn, pointer_transfrom,
-                               h5filename](auto const& component_transform,
-                                           int const time_step,
-                                           int const max_step,
-                                           time_attribute const variable)
+                               h5filename,
+                               mesh_name](auto const& component_transform,
+                                          unsigned long long const time_step,
+                                          int const max_step,
+                                          time_attribute const variable)
     {
         return [component_transform, time_dataitem_genfn, pointer_transfrom,
-                variable, time_step, max_step, h5filename](XdmfData const& attr)
+                variable, time_step, max_step, h5filename,
+                mesh_name](XdmfData const& attr)
         {
             // new data arrived
             bool const changed =
@@ -238,8 +245,8 @@ std::function<std::string(std::vector<double>)> write_xdmf(
                  time_step == 0);
             if (changed)
             {
-                auto dataitem =
-                    time_dataitem_genfn(time_step, max_step, h5filename);
+                auto dataitem = time_dataitem_genfn(time_step, max_step,
+                                                    h5filename, mesh_name);
                 return component_transform(attr, dataitem);
             }
             else
diff --git a/MeshLib/IO/XDMF/writeXdmf.h b/MeshLib/IO/XDMF/writeXdmf.h
index d713a4fb6605379b1386ae3598f63daab0d76935..843313b8a35f2558e57a9caaf04ff121edb2b48f 100644
--- a/MeshLib/IO/XDMF/writeXdmf.h
+++ b/MeshLib/IO/XDMF/writeXdmf.h
@@ -33,11 +33,13 @@ namespace MeshLib::IO
  * @param constant_attributes Meta data for attributes NOT changing over time
  * @param h5filename Name of the file where the actual data was written
  * @param ogs_version OGS Version to be added to XdmfInformation tag
+ * @param mesh_name Name of the output mesh
  * @return unary function with vector of time step values, returning XDMF string
  */
 std::function<std::string(std::vector<double>)> write_xdmf(
     XdmfData const& geometry, XdmfData const& topology,
     std::vector<XdmfData> const& constant_attributes,
     std::vector<XdmfData> const& variable_attributes,
-    std::string const& h5filename, std::string const& ogs_version);
+    std::string const& h5filename, std::string const& ogs_version,
+    std::string const& mesh_name);
 }  // namespace MeshLib::IO
\ No newline at end of file
diff --git a/MeshLib/IO/writeMeshToFile.cpp b/MeshLib/IO/writeMeshToFile.cpp
index 3574308774f1f3323ead21ca56a98ebcc5200ee7..36b0ef47815aee160fc0458641ff204eecfb9fa9 100644
--- a/MeshLib/IO/writeMeshToFile.cpp
+++ b/MeshLib/IO/writeMeshToFile.cpp
@@ -9,6 +9,8 @@
 
 #include "writeMeshToFile.h"
 
+#include <vector>
+
 #include "BaseLib/FileTools.h"
 #include "BaseLib/Logging.h"
 #include "BaseLib/StringTools.h"
@@ -45,7 +47,10 @@ int writeMeshToFile(const MeshLib::Mesh& mesh,
     }
     if (file_path.extension().string() == ".xdmf")
     {
-        MeshLib::IO::XdmfHdfWriter(mesh, file_path, 0, 0.0,
+        std::vector<std::reference_wrapper<const MeshLib::Mesh>> meshes;
+        const std::reference_wrapper<const MeshLib::Mesh> mr = mesh;
+        meshes.push_back(mr);
+        MeshLib::IO::XdmfHdfWriter(std::move(meshes), file_path, 0, 0.0,
                                    variable_output_names, true);
         return 0;
     }
diff --git a/ProcessLib/LiquidFlow/Tests.cmake b/ProcessLib/LiquidFlow/Tests.cmake
index a302aca454319a6fdada998be8a693b33b576a9d..401fc4ae88a0f54188f46c0f8fde178e99fd8a45 100644
--- a/ProcessLib/LiquidFlow/Tests.cmake
+++ b/ProcessLib/LiquidFlow/Tests.cmake
@@ -481,15 +481,43 @@ AddTest(
     EXECUTABLE_ARGS FunctionParameterTest_XDMF.prj
     WRAPPER time
     TESTER xdmfdiff
-    # See https://gitlab.opengeosys.org/ogs/ogs/-/merge_requests/3184#note_85104
-    REQUIREMENTS NOT OGS_USE_MPI AND NOT COMPILER_IS_APPLE_CLANG
+    REQUIREMENTS NOT OGS_USE_MPI
+    DIFF_DATA
+    square_5x5_tris_32.xdmf square_5x5_tris_32.xdmf pressure pressure 1e-7 1e-13
+    square_5x5_tris_32.xdmf square_5x5_tris_32.xdmf HydraulicFlow HydraulicFlow 1e-7 1e-13
+    square_5x5_tris_32.xdmf square_5x5_tris_32.xdmf MaterialIDs MaterialIDs 1e-7 1e-13
+    square_5x5_tris_32.xdmf square_5x5_tris_32.xdmf v v 1e-7 1e-13
+    square_5x5_tris_32_right_boundary.xdmf square_5x5_tris_32_right_boundary.xdmf pressure pressure 1e-7 1e-13
+    square_5x5_tris_32_right_boundary.xdmf square_5x5_tris_32_right_boundary.xdmf bulk_element_ids bulk_element_ids 1e-7 1e-13
+    square_5x5_tris_32_right_boundary.xdmf square_5x5_tris_32_right_boundary.xdmf bulk_node_ids bulk_node_ids 1e-7 1e-13
+    square_5x5_tris_32_left_boundary.xdmf square_5x5_tris_32_left_boundary.xdmf pressure pressure 1e-7 1e-13
+    square_5x5_tris_32_left_boundary.xdmf square_5x5_tris_32_left_boundary.xdmf bulk_element_ids bulk_element_ids 1e-7 1e-13
+    square_5x5_tris_32_left_boundary.xdmf square_5x5_tris_32_left_boundary.xdmf bulk_node_ids bulk_node_ids 1e-7 1e-13
+)
+
+AddTest(
+    NAME SimpleSynthetics_XDMF_compression_off
+    PATH Parabolic/LiquidFlow/SimpleSynthetics/XDMF_compression_off
+    EXECUTABLE ogs
+    EXECUTABLE_ARGS FunctionParameterTest_XDMF.prj
+    WRAPPER time
+    TESTER xdmfdiff
+    REQUIREMENTS NOT OGS_USE_MPI
     DIFF_DATA
     square_5x5_tris_32.xdmf square_5x5_tris_32.xdmf pressure pressure 1e-7 1e-13
     square_5x5_tris_32.xdmf square_5x5_tris_32.xdmf HydraulicFlow HydraulicFlow 1e-7 1e-13
     square_5x5_tris_32.xdmf square_5x5_tris_32.xdmf MaterialIDs MaterialIDs 1e-7 1e-13
     square_5x5_tris_32.xdmf square_5x5_tris_32.xdmf v v 1e-7 1e-13
+    square_5x5_tris_32_right_boundary.xdmf square_5x5_tris_32_right_boundary.xdmf pressure pressure 1e-7 1e-13
+    square_5x5_tris_32_right_boundary.xdmf square_5x5_tris_32_right_boundary.xdmf bulk_element_ids bulk_element_ids 1e-7 1e-13
+    square_5x5_tris_32_right_boundary.xdmf square_5x5_tris_32_right_boundary.xdmf bulk_node_ids bulk_node_ids 1e-7 1e-13
+    square_5x5_tris_32_left_boundary.xdmf square_5x5_tris_32_left_boundary.xdmf pressure pressure 1e-7 1e-13
+    square_5x5_tris_32_left_boundary.xdmf square_5x5_tris_32_left_boundary.xdmf bulk_element_ids bulk_element_ids 1e-7 1e-13
+    square_5x5_tris_32_left_boundary.xdmf square_5x5_tris_32_left_boundary.xdmf bulk_node_ids bulk_node_ids 1e-7 1e-13
 )
 
+
+
 #AddTest(
 #    NAME LiquidFlow_SimpleSynthetics_constraint_dirichlet_bc
 #    PATH Parabolic/LiquidFlow/SimpleSynthetics
diff --git a/ProcessLib/Output/Output.cpp b/ProcessLib/Output/Output.cpp
index d1340797473209176bd537debc3f0c1b3dcca218..c9376decfe85ea6dda33c9462b003300376653d6 100644
--- a/ProcessLib/Output/Output.cpp
+++ b/ProcessLib/Output/Output.cpp
@@ -233,10 +233,12 @@ struct Output::OutputFile
     }
 };
 
-void Output::outputMeshXdmf(OutputFile const& output_file,
-                            MeshLib::Mesh const& mesh,
-                            int const timestep,
-                            double const t)
+void Output::outputMeshXdmf(
+    OutputFile const& output_file,
+    std::vector<std::reference_wrapper<const MeshLib::Mesh>>
+        meshes,
+    int const timestep,
+    double const t)
 {
     // \TODO (tm) Refactor to a dedicated VTKOutput and XdmfOutput
     // The XdmfOutput will create on construction the XdmfHdfWriter
@@ -244,14 +246,13 @@ void Output::outputMeshXdmf(OutputFile const& output_file,
     {
         std::filesystem::path path(output_file.path);
         _mesh_xdmf_hdf_writer = std::make_unique<MeshLib::IO::XdmfHdfWriter>(
-            MeshLib::IO::XdmfHdfWriter(
-                mesh, path, timestep, t,
-                _output_data_specification.output_variables,
-                output_file.compression));
+            std::move(meshes), path, timestep, t,
+            _output_data_specification.output_variables,
+            output_file.compression);
     }
     else
     {
-        _mesh_xdmf_hdf_writer->writeStep(timestep, t);
+        _mesh_xdmf_hdf_writer->writeStep(t);
     };
 }
 
@@ -308,38 +309,45 @@ void Output::doOutputAlways(Process const& process,
         return;
     }
 
-    auto output_bulk_mesh = [&](MeshLib::Mesh const& mesh)
+    auto output_bulk_mesh =
+        [&](std::vector<std::reference_wrapper<const MeshLib::Mesh>> meshes)
     {
         MeshLib::IO::PVDFile* pvd_file = nullptr;
         if (_output_file_type == ProcessLib::OutputType::vtk)
         {
-            OutputFile const file(
-                _output_directory, _output_file_type, _output_file_prefix,
-                _output_file_suffix, mesh.getName(), timestep, t, iteration,
-                _output_file_data_mode, _output_file_compression,
-                _output_data_specification.output_variables);
-
-            pvd_file = findPVDFile(process, process_id, mesh.getName());
-            outputMesh(file, pvd_file, mesh, t);
+            for (auto const& mesh : meshes)
+            {
+                OutputFile const file(
+                    _output_directory, _output_file_type, _output_file_prefix,
+                    _output_file_suffix, mesh.get().getName(), timestep, t,
+                    iteration, _output_file_data_mode, _output_file_compression,
+                    _output_data_specification.output_variables);
+
+                pvd_file =
+                    findPVDFile(process, process_id, mesh.get().getName());
+                outputMesh(file, pvd_file, mesh, t);
+            }
         }
         else if (_output_file_type == ProcessLib::OutputType::xdmf)
         {
-            OutputFile const file(
-                _output_directory, _output_file_type, _output_file_prefix, "",
-                mesh.getName(), timestep, t, iteration, _output_file_data_mode,
-                _output_file_compression,
-                _output_data_specification.output_variables);
-
-            outputMeshXdmf(file, mesh, timestep, t);
+            std::string name = meshes[0].get().getName();
+            OutputFile const file(_output_directory, _output_file_type,
+                                  _output_file_prefix, "", name, timestep, t,
+                                  iteration, _output_file_data_mode,
+                                  _output_file_compression,
+                                  _output_data_specification.output_variables);
+
+            outputMeshXdmf(std::move(file), std::move(meshes), timestep, t);
         }
     };
 
+    std::vector<std::reference_wrapper<const MeshLib::Mesh>> output_meshes;
     for (auto const& mesh_output_name : _mesh_names_for_output)
     {
         // process related output
         if (process.getMesh().getName() == mesh_output_name)
         {
-            output_bulk_mesh(process.getMesh());
+            output_meshes.push_back(process.getMesh());
             continue;
         }
         // mesh related output
@@ -378,8 +386,10 @@ void Output::doOutputAlways(Process const& process,
             process.getIntegrationPointWriter(non_bulk_mesh),
             _output_data_specification);
 
-        output_bulk_mesh(non_bulk_mesh);
+        output_meshes.push_back(non_bulk_mesh);
     }
+
+    output_bulk_mesh(std::move(output_meshes));
     INFO("[time] Output of timestep {:d} took {:g} s.", timestep,
          time_output.elapsed());
 }
diff --git a/ProcessLib/Output/Output.h b/ProcessLib/Output/Output.h
index 78c0bbee5e8c6d6255fe2f411223054623a38cbd..6e530c7f220dbe3e33034f5f747ea905bbb24eb6 100644
--- a/ProcessLib/Output/Output.h
+++ b/ProcessLib/Output/Output.h
@@ -88,12 +88,13 @@ private:
     struct OutputFile;
 
     static void outputMesh(OutputFile const& output_file,
-                    MeshLib::IO::PVDFile* const pvd_file,
-                    MeshLib::Mesh const& mesh,
-                    double const t);
+                           MeshLib::IO::PVDFile* const pvd_file,
+                           MeshLib::Mesh const& mesh,
+                           double const t);
 
     void outputMeshXdmf(OutputFile const& output_file,
-                        MeshLib::Mesh const& mesh,
+                        std::vector<std::reference_wrapper<const MeshLib::Mesh>>
+                            meshes,
                         int const timestep,
                         double const t);
 
diff --git a/ProcessLib/SteadyStateDiffusion/Tests.cmake b/ProcessLib/SteadyStateDiffusion/Tests.cmake
index 9dd07c7fcee15b08bd0024f83744ee67eeba8672..688763e93af237bb9a713c0c1f4e1da8d851e9ee 100644
--- a/ProcessLib/SteadyStateDiffusion/Tests.cmake
+++ b/ProcessLib/SteadyStateDiffusion/Tests.cmake
@@ -425,8 +425,8 @@ AddTest(
     TESTER xdmfdiff
     REQUIREMENTS OGS_USE_MPI
     DIFF_DATA
-    cube_1e3_np3.xdmf cube_1e3_np3.xdmf pressure pressure 1e-3 1e-3
-    cube_1e3_np3.xdmf cube_1e3_np3.xdmf v v 1e-3 1e-3
+    cube_1e3_np3.xdmf cube_1x1x1_hex_1e3.xdmf pressure pressure 1e-3 1e-3
+    cube_1e3_np3.xdmf cube_1x1x1_hex_1e3.xdmf v v 1e-3 1e-3
 )
 
 AddTest(
@@ -439,8 +439,8 @@ AddTest(
     TESTER xdmfdiff
     REQUIREMENTS OGS_USE_MPI
     DIFF_DATA
-    cube_1e3_np2.xdmf cube_1e3_np2.xdmf pressure pressure 1e-3 1e-3
-    cube_1e3_np2.xdmf cube_1e3_np2.xdmf v v 1e-3 1e-3
+    cube_1e3_np2.xdmf cube_1x1x1_hex_1e3.xdmf pressure pressure 1e-3 1e-3
+    cube_1e3_np2.xdmf cube_1x1x1_hex_1e3.xdmf v v 1e-3 1e-3
 )
 
 AddTest(