VTIWriter.h

Go to the documentation of this file.

// ==============================================================================================
// Microflow 3D, http://www.microflow.pwr.edu.pl/
// Created by Roman Szafran on 26.03.19.
// Copyright (c) 2019 Wroclaw University of Science and Technology.
// Distributed under the Apache License, Version 2.0. You may obtain a copy of the License at
// http://www.apache.org/licenses/LICENSE-2.0 or see accompanying file license.txt.
// Redistributions of source code must retain the above copyright and license notice.
// ==============================================================================================

#pragma once

#include <iostream>
#include <string>
#include <cstring>
#include <vector>
#include <cmath>

// OpenVDB
#include <openvdb/openvdb.h>

// VTK
#include <vtk-7.1/vtkSmartPointer.h>
#include <vtk-7.1/vtkXMLImageDataWriter.h>
#include <vtk-7.1/vtkXMLImageDataReader.h>
#include <vtk-7.1/vtkImageData.h>
#include <vtk-7.1/vtkPointData.h>
#include <vtk-7.1/vtkCellData.h>
#include <vtk-7.1/vtkDataArray.h>
#include <vtk-7.1/vtkDoubleArray.h>
#include <vtk-7.1/vtkFloatArray.h>
#include <vtk-7.1/vtkIntArray.h>
#include <vtk-7.1/vtkTypeUInt64Array.h>

#include "MFDatabase/ThreadArray.h"
#include "MFDatabase/Vec4.h"
#include "LBParametersConverter/ParametersConversion.h"

namespace MF::RW {


        template<class T_VDBGridType> // T_VDBGridType = <openvdb::FloatGrid::Ptr>, <openvdb::Int32Grid::Ptr>, <openvdb::DobleGrid::Ptr> etc.
        class VTIWriter {
        public:
            VTIWriter(std::shared_ptr<T_VDBGridType> VDBGrid_Ptr, const std::string &FileName, float voxelSize) {
                m_FileNameExt = FileName + ".vti";
                // The bounding box for active voxels
                VDBGrid_Ptr->treePtr()->evalActiveVoxelBoundingBox(m_bbox);
                // Smart pointer to vtkImageData object
                m_VTKimageData_Ptr = vtkSmartPointer<vtkImageData>::New();
                m_VTKimageData_Ptr->SetDimensions(m_bbox.dim().x(), m_bbox.dim().y(), m_bbox.dim().z());
                m_VTKimageData_Ptr->SetExtent(m_bbox.min().x(), m_bbox.max().x(), m_bbox.min().y(), m_bbox.max().y(), m_bbox.min().z(), m_bbox.max().z());
                m_VTKimageData_Ptr->SetSpacing(voxelSize, voxelSize, voxelSize);
                m_pointNumber = m_VTKimageData_Ptr->GetNumberOfPoints();
                m_VDBGrid_Ptr = VDBGrid_Ptr;
            };

            ~VTIWriter() = default;

            // T_VDBGridType = <openvdb::FloatGrid::Ptr>, <openvdb::Int32Grid::Ptr>, <openvdb::DoubleGrid::Ptr> etc.
            static std::shared_ptr<MF::RW::VTIWriter<T_VDBGridType>> New(std::shared_ptr<T_VDBGridType> VDBGrid_Ptr, const std::string &FileName, float voxelSize) {
                auto VTIWriter_Ptr = std::make_shared<MF::RW::VTIWriter<T_VDBGridType>>(VDBGrid_Ptr, FileName, voxelSize);
                return VTIWriter_Ptr;
            };

        private:

            vtkSmartPointer<vtkImageData> m_VTKimageData_Ptr;
            std::shared_ptr<T_VDBGridType>  m_VDBGrid_Ptr;
            std::string m_FileNameExt;
            openvdb::CoordBBox m_bbox;
            uint64_t m_pointNumber;

        public:

            template<typename T3> // T3 = <vtkFloatArray>, <vtkDoubleArray>, <vtkIntArray> etc.
            void addDataFromGrid(const std::string &DataName, unsigned int numComponents, std::shared_ptr<T_VDBGridType> VDBGrid_Ptr){
                auto newData = vtkSmartPointer<T3>::New();
                newData->SetNumberOfComponents(numComponents);
                newData->SetNumberOfTuples(m_VTKimageData_Ptr->GetNumberOfPoints());
                newData->SetName(DataName.c_str());

                // Get an accessor for coordinate-based access to voxels.
                auto accessor = VDBGrid_Ptr->getAccessor();
                openvdb::Coord ijk;
                uint64_t nr = 0;
                double tuple[numComponents];
                for (unsigned int n=0; n<numComponents; ++n) tuple[n] = NAN;
                int32_t &i = ijk[0], &j = ijk[1], &k = ijk[2];

                for (k = m_bbox.min().z(); k <= m_bbox.max().z(); ++k) {
                    for (j = m_bbox.min().y(); j <= m_bbox.max().y(); ++j) {
                        for (i = m_bbox.min().x(); i <= m_bbox.max().x(); ++i)  {
                            for(unsigned int n=0; n<numComponents; ++n)
                                tuple[n] = accessor.getValue(ijk);
                            newData->SetTuple(nr, tuple);
                            nr ++;
                        }
                    }
                }
                m_VTKimageData_Ptr->GetPointData()->AddArray(newData);
            };

            template<typename T3> // T3 = <vtkFloatArray>, <vtkDoubleArray>, <vtkIntArray> etc.
            void addDataFromGrid(const std::string &DataName, unsigned int numComponents){
                auto newData = vtkSmartPointer<T3>::New();
                newData->SetNumberOfComponents(numComponents);
                newData->SetNumberOfTuples(m_pointNumber);
                newData->SetName(DataName.c_str());

                // Get an accessor for coordinate-based access to voxels.
                auto accessor = m_VDBGrid_Ptr->getAccessor();
                openvdb::Coord ijk;
                uint64_t nr = 0;
                double tuple[numComponents]; //?! Why only double or float allowed?
                for (unsigned int n=0; n<numComponents; ++n)
                    tuple[n] = NAN;

                int32_t &i = ijk[0], &j = ijk[1], &k = ijk[2];

                for (k = m_bbox.min().z(); k <= m_bbox.max().z(); ++k) {
                    for (j = m_bbox.min().y(); j <= m_bbox.max().y(); ++j) {
                        for (i = m_bbox.min().x(); i <= m_bbox.max().x(); ++i)  {
                            for(unsigned int n=0; n<numComponents; ++n)
                                tuple[n] = accessor.getValue(ijk);
                            newData->SetTuple(nr, tuple);
                            nr ++;
                        }
                    }
                }
                m_VTKimageData_Ptr->GetPointData()->AddArray(newData);
            };

            void addNodeIDFromGrid(){
                auto newData = vtkSmartPointer<vtkTypeUInt64Array>::New();
                newData->SetNumberOfComponents(1);
                newData->SetNumberOfTuples(m_pointNumber);
                newData->SetName("Node ID");

                // Get an accessor for coordinate-based access to voxels.
                auto accessor = m_VDBGrid_Ptr->getAccessor();
                openvdb::Coord ijk;
                uint64_t nr = 0;
                double tuple[1] = {0};
                int32_t &i = ijk[0], &j = ijk[1], &k = ijk[2];

                for (k = m_bbox.min().z(); k <= m_bbox.max().z(); ++k) {
                    for (j = m_bbox.min().y(); j <= m_bbox.max().y(); ++j) {
                        for (i = m_bbox.min().x(); i <= m_bbox.max().x(); ++i) {
                            tuple[0] = accessor.getValue(ijk);
                            newData->SetTuple(nr, tuple);
                            nr ++;
                        }
                    }
                }
                m_VTKimageData_Ptr->GetPointData()->AddArray(newData);
            };

            void addNodeCoordinateFromGrid(){
                auto newData = vtkSmartPointer<vtkIntArray>::New();
                newData->SetNumberOfComponents(3);
                newData->SetNumberOfTuples(m_pointNumber);
                newData->SetName("MF NodeCoord");

                // Get an accessor for coordinate-based access to voxels.
                auto accessor = m_VDBGrid_Ptr->getAccessor();
                openvdb::Coord ijk;
                uint64_t nr = 0;
                double tuple[3] ={NAN,NAN,NAN};
                for (auto i = 0; i < m_pointNumber; ++i) {
                    newData->SetTuple(i, tuple);
                }
                int32_t &i = ijk[0], &j = ijk[1], &k = ijk[2];

                for (k = m_bbox.min().z(); k <= m_bbox.max().z(); ++k) {
                    for (j = m_bbox.min().y(); j <= m_bbox.max().y(); ++j) {
                        for (i = m_bbox.min().x(); i <= m_bbox.max().x(); ++i) {
                            if (accessor.getValue(ijk) != 0) {
                                tuple[0] = i;
                                tuple[1] = j;
                                tuple[2] = k;
                                newData->SetTuple(nr, tuple);
                            }
                            nr++;
                        }
                    }
                }
                m_VTKimageData_Ptr->GetPointData()->AddArray(newData);
            };

            template<typename T3_VTKArrayType, typename T4> // T1_VTKArrayType = <vtkFloatArray>, <vtkDoubleArray>, <vtkIntArray> etc. T4 = <foat>, <int>, <double> etc.
            void addDataFromTable(const std::string &DataName, unsigned int numComponents, std::vector<T4> &table_Ptr) {
                vtkSmartPointer<T3_VTKArrayType> newData = vtkSmartPointer<T3_VTKArrayType>::New();
                newData->SetNumberOfComponents(numComponents);
                newData->SetNumberOfTuples(m_pointNumber);
                newData->SetName(DataName.c_str());

                double tuple[numComponents];
                for (unsigned int n=0; n<numComponents; ++n)
                    tuple[n] = NAN;

                for (auto i = 0; i < m_pointNumber; ++i) {
                    for(unsigned int n=0; n<numComponents; ++n)
                        tuple[n] = table_Ptr[i*numComponents+n];
                    newData->SetTuple(i, tuple);
                }
                m_VTKimageData_Ptr->GetPointData()->AddArray(newData);
            };

            void addNodeCoordinateFromThread(const std::shared_ptr<MF::Database::ThreadArray>& ThreadArray_Ptr) {
                vtkSmartPointer<vtkIntArray> newData = vtkSmartPointer<vtkIntArray>::New();
                newData->SetNumberOfComponents(3);
                newData->SetNumberOfTuples(m_pointNumber);
                newData->SetName("MF NodeCoord");
                const uint32_t Nx = m_bbox.dim().x();
                const uint32_t Ny = m_bbox.dim().y();

                double tuple[3] = {NAN,NAN,NAN};
                for (auto i = 0; i < m_pointNumber; ++i) {
                    newData->SetTuple(i, tuple);
                }
                for (auto & thread : *ThreadArray_Ptr->m_ThreadsTable_Ptr) {
                    for (size_t node = 0; node < thread->m_NodeArray_Ptr->size(); node++) {
                        tuple[0] = thread->m_Coordinate_Ptr->operator[](node).x;
                        tuple[1] = thread->m_Coordinate_Ptr->operator[](node).y;
                        tuple[2] = thread->m_Coordinate_Ptr->operator[](node).z;
                        const uint32_t X = thread->m_Coordinate_Ptr->operator[](node).x - m_bbox.min().x();
                        const uint32_t Y = thread->m_Coordinate_Ptr->operator[](node).y - m_bbox.min().y();
                        const uint32_t Z = thread->m_Coordinate_Ptr->operator[](node).z - m_bbox.min().z();
                        const uint64_t node_nr =  (Z * Ny + Y) * Nx + X;
                        newData->SetTuple(node_nr, tuple);
                    }
                }
                m_VTKimageData_Ptr->GetPointData()->AddArray(newData);
            };

            void addNodeTypeFromThread(const std::shared_ptr<MF::Database::ThreadArray>& ThreadArray_Ptr) {
                vtkSmartPointer<vtkIntArray> newData = vtkSmartPointer<vtkIntArray>::New();
                newData->SetNumberOfComponents(1);
                newData->SetNumberOfTuples(m_pointNumber);
                newData->SetName("Node type");
                const uint32_t Nx = m_bbox.dim().x();
                const uint32_t Ny = m_bbox.dim().y();

                double tuple[1] = {0};
                for (auto i = 0; i < m_pointNumber; ++i) {
                    newData->SetTuple(i, tuple);
                }
                for (auto & thread : *ThreadArray_Ptr->m_ThreadsTable_Ptr) {
                    tuple[0] = thread.get()->m_NodeType;
                    for (size_t node = 0; node < thread->m_NodeArray_Ptr->size(); node++) {
                        const uint32_t X = thread->m_Coordinate_Ptr->operator[](node).x - m_bbox.min().x();
                        const uint32_t Y = thread->m_Coordinate_Ptr->operator[](node).y - m_bbox.min().y();
                        const uint32_t Z = thread->m_Coordinate_Ptr->operator[](node).z - m_bbox.min().z();
                        const uint64_t node_nr =  (Z * Ny + Y) * Nx + X;
                        newData->SetTuple(node_nr, tuple);
                    }
                }
                m_VTKimageData_Ptr->GetPointData()->AddArray(newData);
            };

            void addUidThreadFromThread(const std::shared_ptr<MF::Database::ThreadArray>& ThreadArray_Ptr) {
                vtkSmartPointer<vtkIntArray> newData = vtkSmartPointer<vtkIntArray>::New();
                newData->SetNumberOfComponents(1);
                newData->SetNumberOfTuples(m_pointNumber);
                newData->SetName("Node uidThreadNr");
                const uint32_t Nx = m_bbox.dim().x();
                const uint32_t Ny = m_bbox.dim().y();

                double tuple[1] = {NAN};
                for (auto i = 0; i < m_pointNumber; ++i) {
                    newData->SetTuple(i, tuple);
                }
                for (auto & thread : *ThreadArray_Ptr->m_ThreadsTable_Ptr) {
                    tuple[0] = thread.get()->m_uidThreadNr;
                    for (size_t node = 0; node < thread->m_NodeArray_Ptr->size(); node++) {
                        const uint32_t X = thread->m_Coordinate_Ptr->operator[](node).x - m_bbox.min().x();
                        const uint32_t Y = thread->m_Coordinate_Ptr->operator[](node).y - m_bbox.min().y();
                        const uint32_t Z = thread->m_Coordinate_Ptr->operator[](node).z - m_bbox.min().z();
                        const uint64_t node_nr =  (Z * Ny + Y) * Nx + X;
                        newData->SetTuple(node_nr, tuple);
                    }
                }
                m_VTKimageData_Ptr->GetPointData()->AddArray(newData);
            };

            void addNodeThreadsByNameThread(const std::shared_ptr<MF::Database::ThreadArray>& ThreadArray_Ptr) {

                auto threadVectorName = std::vector<std::string>();
                for (auto & thread : *ThreadArray_Ptr->m_ThreadsTable_Ptr) {
                    std::string threadName = thread.get()->m_ThreadName;
                    if (!threadName.empty()) {
                        unsigned int number = 0;
                        for (auto & vec_it : threadVectorName) {
                            if (vec_it == threadName)
                                number++;
                        }
                        if (number == 0)
                            threadVectorName.push_back(threadName);
                    }
                }

                for (auto & vec_it : threadVectorName) {
                    vtkSmartPointer<vtkIntArray> newData = vtkSmartPointer<vtkIntArray>::New();
                    newData->SetNumberOfComponents(1);
                    newData->SetNumberOfTuples(m_pointNumber);
                    newData->SetName(vec_it.c_str());
                    const uint32_t Nx = m_bbox.dim().x();
                    const uint32_t Ny = m_bbox.dim().y();

                    double tuple[1] = {0};
                    for (auto i = 0; i < m_pointNumber; ++i)
                        newData->SetTuple(i, tuple);
                    tuple[0] = 1;
                    for (auto & thread : *ThreadArray_Ptr->m_ThreadsTable_Ptr) {
                        if (vec_it == thread.get()->m_ThreadName) {
                            for (size_t node = 0; node < thread->m_NodeArray_Ptr->size(); node++) {
                                const uint32_t X = thread->m_Coordinate_Ptr->operator[](node).x - m_bbox.min().x();
                                const uint32_t Y = thread->m_Coordinate_Ptr->operator[](node).y - m_bbox.min().y();
                                const uint32_t Z = thread->m_Coordinate_Ptr->operator[](node).z - m_bbox.min().z();
                                const uint64_t node_nr = (Z * Ny + Y) * Nx + X;
                                newData->SetTuple(node_nr, tuple);
                            }
                        }
                    }
                    m_VTKimageData_Ptr->GetPointData()->AddArray(newData);
                }
            };

            void addVelocityLBFromThread(const std::shared_ptr<MF::Database::ThreadArray>& ThreadArray_Ptr) {
                vtkSmartPointer<vtkDoubleArray> newData = vtkSmartPointer<vtkDoubleArray>::New();
                newData->SetNumberOfComponents(3);
                newData->SetNumberOfTuples(m_pointNumber);
                newData->SetName("Velocity LB [-]");
                const uint32_t Nx = m_bbox.dim().x();
                const uint32_t Ny = m_bbox.dim().y();

                double tuple[3] = {NAN,NAN,NAN};
                for (auto i = 0; i < m_pointNumber; ++i) {
                    newData->SetTuple(i, tuple);
                }
                for (auto & thread : *ThreadArray_Ptr->m_ThreadsTable_Ptr) {
                    for (size_t node = 0; node < thread->m_NodeArray_Ptr->size(); node++) {
                        const MF::Database::Vec4<double> Vector4 = (*thread->m_pVRLBFunction)(&thread->m_NodeArray_Ptr->operator[](node));
                        tuple[0] = Vector4.x;
                        tuple[1] = Vector4.y;
                        tuple[2] = Vector4.z;
                        const uint32_t X = thread->m_Coordinate_Ptr->operator[](node).x - m_bbox.min().x();
                        const uint32_t Y = thread->m_Coordinate_Ptr->operator[](node).y - m_bbox.min().y();
                        const uint32_t Z = thread->m_Coordinate_Ptr->operator[](node).z - m_bbox.min().z();
                        const uint64_t node_nr =  (Z * Ny + Y) * Nx + X;
                        newData->SetTuple(node_nr, tuple);
                    }
                }
                m_VTKimageData_Ptr->GetPointData()->AddArray(newData);
            };

            void addInitialVelocityLBFromThread(const std::shared_ptr<MF::Database::ThreadArray>& ThreadArray_Ptr) {
                vtkSmartPointer<vtkDoubleArray> newData = vtkSmartPointer<vtkDoubleArray>::New();
                newData->SetNumberOfComponents(3);
                newData->SetNumberOfTuples(m_pointNumber);
                newData->SetName("Initial velocity LB [-]");
                const uint32_t Nx = m_bbox.dim().x();
                const uint32_t Ny = m_bbox.dim().y();

                double tuple[3] = {NAN,NAN,NAN};
                for (auto i = 0; i < m_pointNumber; ++i) {
                    newData->SetTuple(i, tuple);
                }
                for (auto & thread : *ThreadArray_Ptr->m_ThreadsTable_Ptr) {
                    tuple[0] = thread.get()->m_BoundaryValue_ULB.x;
                    tuple[1] = thread.get()->m_BoundaryValue_ULB.y;
                    tuple[2] = thread.get()->m_BoundaryValue_ULB.z;
                    for (size_t node = 0; node < thread->m_NodeArray_Ptr->size(); node++) {
                        const uint32_t X = thread->m_Coordinate_Ptr->operator[](node).x - m_bbox.min().x();
                        const uint32_t Y = thread->m_Coordinate_Ptr->operator[](node).y - m_bbox.min().y();
                        const uint32_t Z = thread->m_Coordinate_Ptr->operator[](node).z - m_bbox.min().z();
                        const uint64_t node_nr =  (Z * Ny + Y) * Nx + X;
                        newData->SetTuple(node_nr, tuple);
                    }
                }
                m_VTKimageData_Ptr->GetPointData()->AddArray(newData);
            };

            void addVelocityPhysFromThread(const std::shared_ptr<MF::Database::ThreadArray>& ThreadArray_Ptr) {
                vtkSmartPointer<vtkDoubleArray> newData = vtkSmartPointer<vtkDoubleArray>::New();
                newData->SetNumberOfComponents(3);
                newData->SetNumberOfTuples(m_pointNumber);
                newData->SetName("Velocity [m/s]");
                const uint32_t Nx = m_bbox.dim().x();
                const uint32_t Ny = m_bbox.dim().y();

                double tuple[3] = {NAN,NAN,NAN};
                for (auto i = 0; i < m_pointNumber; ++i) {
                    newData->SetTuple(i, tuple);
                }
                for (auto & thread : *ThreadArray_Ptr->m_ThreadsTable_Ptr) {
                    for (size_t node = 0; node < thread->m_NodeArray_Ptr->size(); node++) {
                        const  MF::Database::Vec4<double> Vector4 = (*thread.get()->m_pVRLBFunction)(&thread->m_NodeArray_Ptr->operator[](node));
                        tuple[0] = MF::LBPC::ParametersConversion::ULB_toPhys(Vector4.x);
                        tuple[1] = MF::LBPC::ParametersConversion::ULB_toPhys(Vector4.y);
                        tuple[2] = MF::LBPC::ParametersConversion::ULB_toPhys(Vector4.z);
                        const uint32_t X = thread->m_Coordinate_Ptr->operator[](node).x - m_bbox.min().x();
                        const uint32_t Y = thread->m_Coordinate_Ptr->operator[](node).y - m_bbox.min().y();
                        const uint32_t Z = thread->m_Coordinate_Ptr->operator[](node).z - m_bbox.min().z();
                        const uint64_t node_nr =  (Z * Ny + Y) * Nx + X;
                        newData->SetTuple(node_nr, tuple);
                    }
                }
                m_VTKimageData_Ptr->GetPointData()->AddArray(newData);
            };

            void addRhoLBFromThread(const std::shared_ptr<MF::Database::ThreadArray>& ThreadArray_Ptr) {
                vtkSmartPointer<vtkDoubleArray> newData = vtkSmartPointer<vtkDoubleArray>::New();
                newData->SetNumberOfComponents(1);
                newData->SetNumberOfTuples(m_pointNumber);
                newData->SetName("Rho LB [-]");
                const uint32_t Nx = m_bbox.dim().x();
                const uint32_t Ny = m_bbox.dim().y();

                double tuple[1] = {NAN};
                for (size_t i = 0; i < m_pointNumber; ++i) {
                    newData->SetTuple(i, tuple);
                }
                for (auto & thread : *ThreadArray_Ptr->m_ThreadsTable_Ptr) {
                    for (size_t node = 0; node < thread->m_NodeArray_Ptr->size(); node++) {
                        const  MF::Database::Vec4<double> Vector4 = (*thread.get()->m_pVRLBFunction)(&thread->m_NodeArray_Ptr->operator[](node));
                        tuple[0] = Vector4.rho;
                        const uint32_t X = thread->m_Coordinate_Ptr->operator[](node).x - m_bbox.min().x();
                        const uint32_t Y = thread->m_Coordinate_Ptr->operator[](node).y - m_bbox.min().y();
                        const uint32_t Z = thread->m_Coordinate_Ptr->operator[](node).z - m_bbox.min().z();
                        const uint64_t node_nr =  (Z * Ny + Y) * Nx + X;
                        newData->SetTuple(node_nr, tuple);
                    }
                }
                m_VTKimageData_Ptr->GetPointData()->AddArray(newData);
            };

            void addInitialRhoLBFromThread(const std::shared_ptr<MF::Database::ThreadArray>& ThreadArray_Ptr) {
                vtkSmartPointer<vtkDoubleArray> newData = vtkSmartPointer<vtkDoubleArray>::New();
                newData->SetNumberOfComponents(1);
                newData->SetNumberOfTuples(m_pointNumber);
                newData->SetName("Initial rho LB [-]");
                const uint32_t Nx = m_bbox.dim().x();
                const uint32_t Ny = m_bbox.dim().y();

                double tuple[1] = {NAN};
                for (size_t i = 0; i < m_pointNumber; ++i) {
                    newData->SetTuple(i, tuple);
                }
                for (auto & thread : *ThreadArray_Ptr->m_ThreadsTable_Ptr) {
                    tuple[0] = thread.get()->m_BoundaryValue_RhoLB;
                    for (size_t node = 0; node < thread->m_NodeArray_Ptr->size(); node++) {
                        const uint32_t X = thread->m_Coordinate_Ptr->operator[](node).x - m_bbox.min().x();
                        const uint32_t Y = thread->m_Coordinate_Ptr->operator[](node).y - m_bbox.min().y();
                        const uint32_t Z = thread->m_Coordinate_Ptr->operator[](node).z - m_bbox.min().z();
                        const uint64_t node_nr =  (Z * Ny + Y) * Nx + X;
                        newData->SetTuple(node_nr, tuple);
                    }
                }
                m_VTKimageData_Ptr->GetPointData()->AddArray(newData);
            };

            void addRhoPhysFromThread(const std::shared_ptr<MF::Database::ThreadArray>& ThreadArray_Ptr) {
                vtkSmartPointer<vtkDoubleArray> newData = vtkSmartPointer<vtkDoubleArray>::New();
                newData->SetNumberOfComponents(1);
                newData->SetNumberOfTuples(m_pointNumber);
                newData->SetName("Rho [kg/m3]");
                const uint32_t Nx = m_bbox.dim().x();
                const uint32_t Ny = m_bbox.dim().y();

                double tuple[1] = {NAN};
                for (auto i = 0; i < m_pointNumber; ++i) {
                    newData->SetTuple(i, tuple);
                }
                for (auto & thread : *ThreadArray_Ptr->m_ThreadsTable_Ptr) {
                    for (size_t node = 0; node < thread->m_NodeArray_Ptr->size(); node++) {
                        const MF::Database::Vec4<double> Vector4 = (*thread.get()->m_pVRLBFunction)(&thread->m_NodeArray_Ptr->operator[](node));
                        tuple[0] = MF::LBPC::ParametersConversion::RhoLB_toRhoPhys(Vector4.rho);
                        const uint32_t X = thread->m_Coordinate_Ptr->operator[](node).x - m_bbox.min().x();
                        const uint32_t Y = thread->m_Coordinate_Ptr->operator[](node).y - m_bbox.min().y();
                        const uint32_t Z = thread->m_Coordinate_Ptr->operator[](node).z - m_bbox.min().z();
                        const uint64_t node_nr =  (Z * Ny + Y) * Nx + X;
                        newData->SetTuple(node_nr, tuple);
                    }
                }
                m_VTKimageData_Ptr->GetPointData()->AddArray(newData);
            };

            void addPressPhysFromThread(const std::shared_ptr<MF::Database::ThreadArray>& ThreadArray_Ptr) {
                vtkSmartPointer<vtkDoubleArray> newData = vtkSmartPointer<vtkDoubleArray>::New();
                newData->SetNumberOfComponents(1);
                newData->SetNumberOfTuples(m_pointNumber);
                newData->SetName("Gauge pressure [Pa]");
                const uint32_t Nx = m_bbox.dim().x();
                const uint32_t Ny = m_bbox.dim().y();

                double tuple[1] = {NAN};
                for (auto i = 0; i < m_pointNumber; ++i) {
                    newData->SetTuple(i, tuple);
                }
                for (auto & thread : *ThreadArray_Ptr->m_ThreadsTable_Ptr) {
                    for (size_t node = 0; node < thread->m_NodeArray_Ptr->size(); node++) {
                        const MF::Database::Vec4<double> Vector4 = (*thread.get()->m_pVRLBFunction)(&thread->m_NodeArray_Ptr->operator[](node));
                        tuple[0] = MF::LBPC::ParametersConversion::RhoLB_toPressurePhys(Vector4.rho);
                        const uint32_t X = thread->m_Coordinate_Ptr->operator[](node).x - m_bbox.min().x();
                        const uint32_t Y = thread->m_Coordinate_Ptr->operator[](node).y - m_bbox.min().y();
                        const uint32_t Z = thread->m_Coordinate_Ptr->operator[](node).z - m_bbox.min().z();
                        const uint64_t node_nr =  (Z * Ny + Y) * Nx + X;
                        newData->SetTuple(node_nr, tuple);
                    }
                }
                m_VTKimageData_Ptr->GetPointData()->AddArray(newData);
            };

            void addFQ19FromThread(const std::shared_ptr<MF::Database::ThreadArray>& ThreadArray_Ptr) {
                vtkSmartPointer<vtkDoubleArray> newData = vtkSmartPointer<vtkDoubleArray>::New();
                newData->SetNumberOfComponents(MFQ19);
                newData->SetNumberOfTuples(m_pointNumber);
                newData->SetName("FQ19");
                const uint32_t Nx = m_bbox.dim().x();
                const uint32_t Ny = m_bbox.dim().y();

                double tuple[MFQ19] = {NAN, NAN, NAN, NAN, NAN, NAN, NAN, NAN, NAN, NAN, NAN, NAN, NAN, NAN, NAN, NAN, NAN, NAN, NAN};
                for (auto i = 0; i < m_pointNumber; ++i) {
                    newData->SetTuple(i, tuple);
                }
                for (auto & thread : *ThreadArray_Ptr->m_ThreadsTable_Ptr) {
                    for (size_t node = 0; node < thread->m_NodeArray_Ptr->size(); node++) {
                        for (unsigned int i = 0; i < MFQ19; i++)
                            tuple[i] = thread->m_NodeArray_Ptr->operator[](node).FQ19[i];
                        const uint32_t X = thread->m_Coordinate_Ptr->operator[](node).x - m_bbox.min().x();
                        const uint32_t Y = thread->m_Coordinate_Ptr->operator[](node).y - m_bbox.min().y();
                        const uint32_t Z = thread->m_Coordinate_Ptr->operator[](node).z - m_bbox.min().z();
                        const uint64_t node_nr =  (Z * Ny + Y) * Nx + X;
                        newData->SetTuple(node_nr, tuple);
                    }
                }
                m_VTKimageData_Ptr->GetPointData()->AddArray(newData);
            };

            void addPropagationDirFromThread(const std::shared_ptr<MF::Database::ThreadArray>& ThreadArray_Ptr) {
                for (uint8_t k = 0; k < MFQ19; k++) {
                    vtkSmartPointer<vtkIntArray> newData = vtkSmartPointer<vtkIntArray>::New();
                    newData->SetNumberOfComponents(3);
                    newData->SetNumberOfTuples(m_pointNumber);
                    const std::string dataName = "PD" + std::to_string(k) + " [x,y,z]";
                    newData->SetName(dataName.c_str());
                    const uint32_t Nx = m_bbox.dim().x();
                    const uint32_t Ny = m_bbox.dim().y();
                    double tuple[3] = {NAN, NAN, NAN};
                    for (auto i = 0; i < m_pointNumber; ++i)
                        newData->SetTuple(i, tuple);

                    for (auto & thread : *ThreadArray_Ptr->m_ThreadsTable_Ptr) {
                        for (size_t node = 0; node < thread->m_NodeArray_Ptr->size(); node++) {
                            if (thread->getLinkedNodePointer(node, k) != nullptr) {
                                auto neighbNodeNr = thread->getLinkedNodePointer(node, k)-> pMyThread_Ptr->operator*().getNodeNr(thread->getLinkedNodePointer(node, k));
                                tuple[0] = thread->getLinkedNodePointer(node, k)-> pMyThread_Ptr->operator*().getNodeCoordX(neighbNodeNr);
                                tuple[1] = thread->getLinkedNodePointer(node, k)-> pMyThread_Ptr->operator*().getNodeCoordY(neighbNodeNr);
                                tuple[2] = thread->getLinkedNodePointer(node, k)-> pMyThread_Ptr->operator*().getNodeCoordZ(neighbNodeNr);
                                const uint32_t X = thread->getNodeCoordX(node) - m_bbox.min().x();
                                const uint32_t Y = thread->getNodeCoordY(node) - m_bbox.min().y();
                                const uint32_t Z = thread->getNodeCoordZ(node) - m_bbox.min().z();
                                const uint64_t node_nr = (Z * Ny + Y) * Nx + X;
                                newData->SetTuple(node_nr, tuple);
                            }
                        }
                    }
                    m_VTKimageData_Ptr->GetPointData()->AddArray(newData);
                }
            };

            void writeFile(){
                vtkSmartPointer<vtkXMLImageDataWriter> writer = vtkSmartPointer<vtkXMLImageDataWriter>::New();
                writer->SetFileName(m_FileNameExt.c_str());
                writer->SetInputData(m_VTKimageData_Ptr);
                writer->Write();
                std::cout << "VTI file was saved successfully: ----> " << m_FileNameExt << std::endl;
            }
        };

} /* namespase MF::RW */