fs_measurement.hpp

Go to the documentation of this file.

/* 
 * SPDX-FileCopyrightText: 2025 Gesellschaft fuer Anlagen- und Reaktorsicherheit gGmbH
 * SPDX-License-Identifier: EUPL-1.2
 * SPDX-FileContributor: Tim Schoen
 * SPDX-FileContributor: Goethe Universität Frankfurt
 * SPDX-FileType: SOURCE
 * 
 * This file is part of d3f++.
 * d3f++ is an extension for UG4. Licensing information and citation requirements of UG4 are provided in LICENSES/UG4-LGPL_2.1
 */
 
#ifndef __H__UG__PLUGINS__D3F__FS_MEASUREMENT
#define __H__UG__PLUGINS__D3F__FS_MEASUREMENT
 
 
#include "level_set_pos.h"
#include <lib_disc/time_disc/time_integrator_observers/time_integrator_observer_interface.h>
#include <lib_disc/time_disc/finished_conditions.hpp>
#include <fstream>
#include <cmath>
 
#ifdef UG_PARALLEL
#include "pcl/pcl_base.h"
#include "pcl/pcl_process_communicator.h"
#endif
 
namespace ug {
 
 
template<class TDomain, class TAlgebra>
class FSFileMeasurer;
 
template<class TDomain, class TAlgebra>
class FSFileMeasurerFinishedObserver : public ITimeIntegratorObserver<TDomain, TAlgebra>
{
  typedef GridFunction<TDomain, TAlgebra> grid_function_type;
    typedef FSFileMeasurer<TDomain, TAlgebra> measurer_type;
 
    public:
        FSFileMeasurerFinishedObserver(SmartPtr<measurer_type> parent) : m_parent(parent) {}
 
        bool step_process(SmartPtr<grid_function_type> u, int step, number time, number dt)
        {
      return m_parent.get()->on_finished(u, step, time, dt);
        }
 
    private:
    SmartPtr<measurer_type> m_parent;
};
 
template<class TDomain, class TAlgebra>
class FSFileMeasurer : public ITimeIntegratorObserver<TDomain, TAlgebra>
{
    private:   
 
        static const int dim = TDomain::dim;
 
        struct t_pnt_data
        {
             MathVector<dim-1> xy; 
            number z; 
            number z_spec; 
            bool z_given; 
            
            bool valid; 
            
            t_pnt_data (MathVector<dim-1> & the_xy, number the_z_spec, bool if_z_given)
            : xy (the_xy), z (0), z_spec (the_z_spec), z_given (if_z_given), valid (false) {};
        };
 
        class IOutputAdapter
        {
            public:
                virtual void output(const std::vector<t_pnt_data>& measurementPoints, int step, number time) {  };
                virtual void finish() {};
                virtual ~IOutputAdapter() {};
 
        };
 
        class ShellOutput : public IOutputAdapter
        {
            public:
                ShellOutput(bool output_difference, std::string& title)
                : m_outputDifference(output_difference),
                  m_title(title)
                {}
 
                void output(const std::vector<t_pnt_data>& measurementPoints, int step, number time) override
                {
                    std::cout << "---------- FS Height --------" << std::endl;
                    for (auto& pnt : measurementPoints)
                    {
                        std::cout << "(";
                        for (size_t i = 0; i < dim-1; i++ )
                        {
                            std::cout << pnt.xy[i];
                            if(i != dim-2)
                            {
                                std::cout << ", ";
                            }
                        }
                        std::cout << ") \t " << time  << "\t";
                        if(pnt.valid)
                        {
                            std::cout << pnt.z;
 
                            if(this->m_outputDifference && pnt.z_given)
                            {
                                std::cout << " - " << pnt.z_spec-pnt.z;
                            }
                        }
                        else
                        {
                            std::cout << "not valid";
                        }
                        
                        if (m_title.length () != 0)
                        {
                          std::cout << " - " << m_title;
                        }
 
                        std::cout << std::endl;
                    }
                    std::cout << "-----------------------------" << std::endl;                    
                }
 
                void finish() override {}
 
 
            private:
                bool m_outputDifference;
                std::string& m_title;
        };
 
        class CSVOutput : public IOutputAdapter
        {
            public:
                CSVOutput(std::string filename, std::string delimiter, bool output_difference) : m_delim(delimiter), m_outputDifference(output_difference)
                {                   
                    m_outputfilename = filename;
 
                    std::fstream output_file;
                    output_file.open(m_outputfilename.c_str(), std::fstream::trunc | std::fstream::out);
                    
                    // print header
                    output_file << "step"<< m_delim << "time" << m_delim;
                    for (size_t i = 0; i < dim-1; i++ )
                    {
                        output_file << "dim" << i << m_delim;
                    }
                    output_file << "z";
 
                    if(this->m_outputDifference)
                        output_file << m_delim << "difference";
                    
                    output_file << std::endl;
 
                    output_file.close();
                }
 
                void output(const std::vector<t_pnt_data>& measurementPoints, int step, number time)  override
                {
                    std::fstream output_file;
                    output_file.open(m_outputfilename.c_str(), std::fstream::app | std::fstream::out);     
                    output_file.precision(15);
 
                    
                    for (auto& pnt : measurementPoints)
                    {
                        output_file << step << m_delim << time << m_delim;
                        for (size_t i = 0; i < dim-1; i++ )
                        {
                            output_file << pnt.xy[i] << m_delim;
                        }
                        if(pnt.valid)
                        {
                            output_file << pnt.z;
                            
                        }
                        else
                        {
                            output_file << "NaN";
                        }
 
                        if(this->m_outputDifference)
                        {
                            if(pnt.valid && pnt.z_given)
                                output_file << m_delim << pnt.z_spec-pnt.z;
                            else
                                output_file << m_delim << "NaN";
                        }
 
                        output_file << std::endl;
                    }    
                    
                    output_file.close();             
                }
 
                void finish() override
                {
                    std::fstream output_file;
                    output_file.open(m_outputfilename.c_str(), std::fstream::app | std::fstream::out);  
                    output_file << "FINISHED" << m_delim << m_delim;
                    for (size_t i = 0; i < dim-1; i++ )
                    {
                        output_file << m_delim;
                    }
                    output_file << m_delim;
 
                    if(this->m_outputDifference)
                        output_file << m_delim;
                    
                    output_file << std::endl;
                }
 
            private:
                std::string m_delim;
                std::string m_outputfilename;
                bool m_outputDifference;
        };
 
        class StepFileOutput : public IOutputAdapter
        {
            public:
                StepFileOutput(std::string filename, std::string delimiter, bool output_difference) : m_delim(delimiter), m_outputfilename(filename), m_outputDifference(output_difference)
                {
                }
 
                void output(const std::vector<t_pnt_data>& measurementPoints, int step, number time) override
                {
                    // Format the integer with leading zeros to a width of 8 digits
                  std::ostringstream oss;
                  oss << std::setw(8) << std::setfill('0') << step;
 
                    std::string filename = m_outputfilename + "-ts-" + oss.str() + ".dat";
 
                    std::fstream output_file;
                    output_file.open(filename.c_str(), std::fstream::trunc | std::fstream::out);     
                    output_file.precision(15);
 
                    for (auto& pnt : measurementPoints)
                    {
                        for (size_t j = 0; j < dim - 1; j++)
                            output_file << pnt.xy[j] << m_delim;
                        if (pnt.valid)
                            output_file << pnt.z;
                        else
                            output_file << "nAn";
                        if (pnt.z_given)
                        {
                            output_file << m_delim << pnt.z_spec;
                            if (this->m_outputDifference)
                            {
                                if (pnt.valid)
                                    output_file << m_delim << pnt.z - pnt.z_spec;
                                else
                                    output_file << m_delim << "nAn";
                            }
                        }
                        output_file << std::endl;
                    }     
                    
                    output_file.close();           
 
                }
 
                void finish() override
                {
                }
 
            private:
                std::string m_delim;
                std::string m_outputfilename;
                bool m_outputDifference;
        };
 
 
        class TableOutput : public IOutputAdapter
        {
            public:
                TableOutput(const std::vector<t_pnt_data>& measurementPoints, std::string filename, std::string delimiter, bool output_difference, std::string& title)
                : m_delim(delimiter),
                  m_outputDifference(output_difference)
                {
                    if (filename.find(".csv") == std::string::npos)
                        m_outputfilename = filename + ".csv";
                    else
                        m_outputfilename = filename;
 
                    std::fstream output_file;
                    output_file.open(m_outputfilename.c_str(), std::fstream::trunc | std::fstream::out);
                    
                    if (title.length () != 0)
                      output_file << "# " << title << std::endl;
                    
                    for (size_t k = 0; k < dim - 1; k++)
                    {
                        output_file << "# ----" + m_delim + "----";
                        for (size_t i = 0; i < measurementPoints.size (); i++)
                        {
                            const t_pnt_data & pnt = measurementPoints[i];
                            output_file << m_delim << pnt.xy[k];
                        }
                        output_file << std::endl;
                    }
                    output_file << "# step" + m_delim + "time";
                    for (size_t i = 0; i < measurementPoints.size (); i++)
                    {
                        const t_pnt_data & pnt = measurementPoints[i];
                        if (pnt.z_given)
                            output_file << m_delim << pnt.z_spec;
                        else
                            output_file << m_delim << "----";
                    }
                    output_file << std::endl;
 
                    output_file.close();
                }
 
                void output(const std::vector<t_pnt_data>& measurementPoints, int step, number time) override
                {
                    std::fstream output_file;
                    output_file.open(m_outputfilename.c_str(), std::fstream::app | std::fstream::out);     
                    output_file.precision(15);
                    
                    for (size_t i = 0; i < measurementPoints.size (); i++)
                    {
                        output_file << step << m_delim << time << m_delim;
                        const t_pnt_data & pnt = measurementPoints[i];
                        if (pnt.valid)
                            output_file << pnt.z;
                        else
                            output_file << "nAn";
                        output_file << std::endl;
                    }  
 
                    output_file.close();         
 
                }
 
                void finish() override
                {
                }
 
            private:
                std::string m_delim;
                std::string m_outputfilename;
                bool m_outputDifference;
        };
 
        class BinaryOutput : public IOutputAdapter
        {
            public:
                BinaryOutput(std::string filename, bool output_difference) : m_outputDifference(output_difference)
                {
                    if (filename.find(".bin") == std::string::npos)
                        m_outputfilename = filename + ".bin";
                    else
                        m_outputfilename = filename;
 
                    std::fstream bin_output_file;
                    bin_output_file.open(m_outputfilename.c_str(), std::fstream::trunc | std::fstream::out | std::fstream::binary); 
                    int toWrite =  dim;
                    bin_output_file.write((char *)&toWrite, sizeof(int));
                    bin_output_file.close();                    
                }
 
                void output(const std::vector<t_pnt_data>& measurementPoints, int step, number time) override
                {
                    std::fstream bin_output_file;
                    bin_output_file.open(m_outputfilename.c_str(), std::fstream::app | std::fstream::out | std::fstream::binary);
 
                    for (auto& pnt : measurementPoints)
                    {
                        char toWrite = 0x01;
                        bin_output_file.write(&toWrite, sizeof(char));
                        bin_output_file.write((char*)&step, sizeof(int));
                        bin_output_file.write((char*)&time, sizeof(number));
 
                        for (size_t i = 0; i < dim-1; i++ )
                        {
                            bin_output_file.write((char*)&pnt.xy[i], sizeof(number));
                        }
 
                        if(pnt.valid)
                            bin_output_file.write((char*)&pnt.z, sizeof(number));
                        else
                        {
                            number NaN = nan("");
                            bin_output_file.write((char*)&NaN, sizeof(number));
                        }
                             
 
                        if(this->m_outputDifference)
                        {
                            if(pnt.valid && pnt.z_given)                            
                            {
                                number diff = pnt.z_spec-pnt.z;
                                bin_output_file.write((char*)&diff, sizeof(number));
 
                            }
                            else
                            {
                                number NaN = nan("");
                                bin_output_file.write((char*)&NaN, sizeof(number));
                            }
                        } 
                    }  
 
                    bin_output_file.close();
                }
 
                void finish() override
                {
                    std::fstream bin_output_file;
                    bin_output_file.open(m_outputfilename.c_str(), std::fstream::app | std::fstream::out | std::fstream::binary);
                    char toWrite = 0x02;
                    bin_output_file.write(&toWrite, sizeof(char));
                    bin_output_file.close();
                }
 
            private:
                std::string m_outputfilename;
                bool m_outputDifference;
        };
 
    public:
  
        typedef ug::GridFunction<TDomain, TAlgebra> grid_function_type;
 
        FSFileMeasurer (SmartPtr<grid_function_type> lsf)
        : m_measurement_interval (1),
          m_height_measurer(lsf)
        {}
 
        FSFileMeasurer (SmartPtr<grid_function_type> lsf, const char* cmp, number zero)
        : m_measurement_interval (1),
          m_height_measurer(lsf, cmp, zero)
        {}
        
        void set_title_text (const char* text)
        {
          m_titleText = text;
        }
 
        void enable_difference_output(bool differenceOutput)
        {
            m_outputDifference = differenceOutput;                      
        }
 
        void enable_csv_output(const char * output_filename)
        {
            #ifdef UG_PARALLEL
            if(pcl::ProcRank() != 0)
                return;
            #endif    
 
            std::string filename = std::string(output_filename);
 
            if (filename.find(".csv") == std::string::npos)
                filename = filename + ".csv";
 
            IOutputAdapter* adapter = new CSVOutput(filename, ",", m_outputDifference);
            m_output_adapters.push_back(make_sp(adapter));
        }
 
        void enable_tsv_output(const char * output_filename)
        {
            #ifdef UG_PARALLEL
            if(pcl::ProcRank() != 0)
                return;
            #endif    
 
 
            std::string filename = std::string(output_filename);
 
 
            if (filename.find(".tsv") == std::string::npos)
                filename = filename + ".tsv";
 
            IOutputAdapter* adapter = new CSVOutput(filename, "\t", m_outputDifference);
            m_output_adapters.push_back(make_sp(adapter));
        }
 
        void enable_table_output(const char * output_filename)
        {
            #ifdef UG_PARALLEL
            if(pcl::ProcRank() != 0)
                return;
            #endif    
 
            IOutputAdapter* adapter = new TableOutput(m_measurement_points, std::string(output_filename), "\t", m_outputDifference, m_titleText);
            m_output_adapters.push_back(make_sp(adapter));
        }
 
        void enable_step_file_output(const char * output_filename)
        {
            #ifdef UG_PARALLEL
            if(pcl::ProcRank() != 0)
                return;
            #endif    
 
            IOutputAdapter* adapter = new StepFileOutput(std::string(output_filename), "\t", m_outputDifference);
            m_output_adapters.push_back(make_sp(adapter));
        }
 
        void enable_binary_output(const char * output_filename)
        {
            #ifdef UG_PARALLEL
            if(pcl::ProcRank() != 0)
                return;
            #endif    
 
            IOutputAdapter* adapter = new BinaryOutput(std::string(output_filename), m_outputDifference);
            m_output_adapters.push_back(make_sp(adapter));
        }
 
        void enable_print_output()
        {
            #ifdef UG_PARALLEL
            if(pcl::ProcRank() != 0)
                return;
            #endif
 
            IOutputAdapter* adapter = new ShellOutput(m_outputDifference, m_titleText);
            m_output_adapters.push_back(make_sp(adapter));
        }
 
        void set_measurement_interval(int measurement_interval)
        {
            if (measurement_interval < 1)
            {
                std::cout << "FSFileMeasurer: Illegal measurement interval " << measurement_interval << std::endl;
            }
            m_measurement_interval = measurement_interval;
        }
 
        void read_measurement_points(const char * filename, bool includesExpectedZ)
        {
 
            std::string full_file_name = FindFileInStandardPaths (filename);
            UG_COND_THROW (full_file_name.empty (), "FSFileMeasurer: Coulnd't locate file '" << filename << "'.");
            
            std::ifstream input (full_file_name.c_str (), std::ifstream::in);
            UG_COND_THROW (input.fail (), "FSFileMeasurer: Cannot open data file '" << full_file_name << "' for input!");
 
            std::string input_line;
            std::istringstream line_stream;
            
            line_stream.exceptions (std::istream::failbit | std::istream::badbit);
            
            while (! input.eof ())
            {
                std::vector<number> xy(dim-1);
                number z_spec = 0;
                
                if (input.fail ())
                    UG_THROW ("FSFileMeasurer: Could not load the points from the file!");
                std::getline (input, input_line);
                if (input_line.length () == 0)
                    continue;
                try
                {
                    line_stream.str (input_line);
                    line_stream.clear ();
                    for (size_t i = 0; i < dim-1; i++)
                        line_stream >> xy[i];
                    if (includesExpectedZ)
                        line_stream >> z_spec;
                }
                catch (std::istream::failure & e)
                {
                    UG_THROW ("FSFileMeasurer: Failed to parse line '" << input_line << "' in the input file.");
                };
                
                if(includesExpectedZ)
                    add_measurement_point (xy, z_spec);
                else
                    add_measurement_point (xy);
            }
        }
 
 
        void add_measurement_point(std::vector<number> point)
        {        
            if(point.size() != TDomain::dim-1 && point.size() != TDomain::dim)
            {
                std::cout << "FSFileMeasurer: free surface measurement point does not have correct dimension " << TDomain::dim-1 << " or " << TDomain::dim << ", if z is given, ignoring!" << std::endl;
                return;
            }
 
            MathVector<dim-1> p;
 
            for(int i = 0; i < TDomain::dim-1; i++)
            {
                p[i] = point[i];
            }
 
            if(point.size() == TDomain::dim)
            { 
                m_measurement_points.push_back(t_pnt_data(p, point[TDomain::dim-1], true));
            }
            else
            {                
                m_measurement_points.push_back(t_pnt_data(p, 0, false));
            }
        }
 
 
        void add_measurement_point(std::vector<number> point, number specified_z)
        {        
            if(point.size() != TDomain::dim-1)
            {
                std::cout << "FSFileMeasurer: free surface measurement point does not have correct dimension " << TDomain::dim-1 << ", ignoring!" << std::endl;
                return;
            }
 
            MathVector<dim-1> p;
 
            for(int i = 0; i < TDomain::dim-1; i++)
            {
                p[i] = point[i];
            }
 
            m_measurement_points.push_back(t_pnt_data(p, specified_z, true));
        }
 
        SmartPtr< FSFileMeasurerFinishedObserver<TDomain,TAlgebra> >  get_finished_observer()
        {
            FSFileMeasurerFinishedObserver<TDomain,TAlgebra> * pointer = new FSFileMeasurerFinishedObserver<TDomain,TAlgebra>(make_sp(this));
 
            return make_sp(pointer);
        }
        
      bool on_finished(SmartPtr<grid_function_type> u, int step, number time, number dt) 
        {   
            #ifdef UG_PARALLEL
            if(pcl::ProcRank() != 0)
                return true;
            #endif
 
            for (auto& outputAdapter : m_output_adapters)
            {
                outputAdapter->finish();
            }
 
            return true;
        }
       
      virtual bool step_process(SmartPtr<grid_function_type> u, int step, number time, number dt)
        {
 
            if (step % m_measurement_interval != 0)
                return true;
 
            m_height_measurer.reinit();
 
            #ifndef UG_PARALLEL
 
            for (size_t i = 0; i < m_measurement_points.size (); i++)
            {
                t_pnt_data & pnt = m_measurement_points[i];
                pnt.valid = m_height_measurer.get_height_at (pnt.xy, pnt.z);
            }
  
            #else // i.e. ifdef UG_PARALLEL
                
            {
                pcl::ProcessCommunicator proc_comm;
                std::vector<number> loc_z (m_measurement_points.size ()), red_z (m_measurement_points.size ());
                
                for (size_t i = 0; i < m_measurement_points.size (); i++)
                    if (! m_height_measurer.get_height_at (m_measurement_points[i].xy, loc_z[i]))
                        loc_z[i] = - std::numeric_limits<number>::max (); // the position is not covered by the free surface in this process
                
                proc_comm.allreduce (loc_z, red_z, PCL_RO_MAX);
                for (size_t i = 0; i < m_measurement_points.size (); i++)
                {
                    t_pnt_data & pnt = m_measurement_points[i];
                    number z_val = red_z[i];
                    if (z_val != - std::numeric_limits<number>::max ()) // i.e. is not set in any process
                    {
                        pnt.z = z_val;
                        pnt.valid = true;
                    }
                    else
                        pnt.valid = false;
                }
            }
 
            if(pcl::ProcRank() != 0)
                return true;
            #endif
 
            for (auto& outputAdapter : m_output_adapters)
            {
                outputAdapter->output(m_measurement_points, step, time);
            }
 
            return true;
        }
 
    private:   
 
        std::vector<t_pnt_data> m_measurement_points;
        std::vector< SmartPtr<IOutputAdapter> > m_output_adapters;
        int m_measurement_interval;
        ug::d3f::LSPositionZ<grid_function_type> m_height_measurer;   
        bool m_outputDifference;
    std::string m_titleText;
};
}
#endif