statistics_data.h

Go to the documentation of this file.

//    |  /           |
//    ' /   __| _` | __|  _ \   __|
//    . \  |   (   | |   (   |\__ `
//   _|\_\_|  \__,_|\__|\___/ ____/
//                   Multi-Physics
//
//  License:         BSD License
//                   Kratos default license: kratos/license.txt
//
//  Main authors:    Jordi Cotela
//
 
#ifndef KRATOS_STATISTICS_DATA_H_INCLUDED
#define KRATOS_STATISTICS_DATA_H_INCLUDED
 
// System includes
#include <string>
#include <iostream>
#include <functional>
 
// External includes
 
// Project includes
#include "includes/define.h"
#include "includes/element.h"
 
namespace Kratos
{
 
 
 
class StatisticsData
{
public:
 
    KRATOS_CLASS_POINTER_DEFINITION(StatisticsData);
 
    typedef std::vector<double> ValueContainerType;
 
    typedef Matrix::iterator1 IntegrationPointDataView;
    typedef Matrix::const_iterator1 IntegrationPointDataConstView;
 
 
    StatisticsData() {}
 
    StatisticsData(StatisticsData const &rOther) {
        mData = rOther.mData;
    }
 
    virtual ~StatisticsData() {}
 
 
    StatisticsData &operator=(StatisticsData const &rOther) {
        mData = rOther.mData;
        return *this;
    }
 
 
 
    void InitializeStorage(Element& rElement, std::size_t MeasurementSize)
    {
        const GeometryData::IntegrationMethod integration_method = rElement.GetIntegrationMethod();
        std::size_t number_of_integration_points = rElement.GetGeometry().IntegrationPointsNumber(integration_method);
 
        mData.resize(number_of_integration_points, MeasurementSize, false);
        mData = ZeroMatrix(number_of_integration_points, MeasurementSize);
    }
 
 
    void UpdateMeasurement(
        const Element& rElement,
        const PointerVector<StatisticsSampler>& rStatisticsSamplers,
        const PointerVector<StatisticsSampler>& rHigherOrderStatistics,
        ValueContainerType& rUpdate,
        const std::size_t NumMeasurements)
    {
        KRATOS_DEBUG_ERROR_IF(NumMeasurements == 0)
        << "Trying to update statistics, but providied number of recorded steps is zero" << std::endl;
 
        const Geometry<Node> &r_geometry = rElement.GetGeometry();
        const GeometryData::IntegrationMethod integration_method = rElement.GetIntegrationMethod();
        Matrix shape_functions;
        typename Geometry<Node>::ShapeFunctionsGradientsType shape_gradients;
        this->CalculateGeometryData(r_geometry, integration_method, shape_functions, shape_gradients);
 
        for (unsigned int g = 0; g < shape_functions.size1(); g++)
        {
            auto N = row(shape_functions, g);
            auto &rDN_DN = shape_gradients[g];
 
            auto it_update_buffer = rUpdate.begin();
            for (auto it_sampler = rStatisticsSamplers.begin(); it_sampler != rStatisticsSamplers.end(); ++it_sampler)
            {
                it_sampler->SampleDataPoint(r_geometry, N, rDN_DN, it_update_buffer);
            }
 
            if (NumMeasurements > 1) { // Second order and higher statistics start from the second iteration
                // loop on higer order statistics. They require the already written rMeasurements, (const) mData and the number of steps as input
                for (auto it_statistic = rHigherOrderStatistics.begin(); it_statistic != rHigherOrderStatistics.end(); ++it_statistic)
                {
                    it_statistic->SampleDataPoint(it_update_buffer, IntegrationPointData(g), rUpdate, NumMeasurements);
                }
            }
 
            for (unsigned int i = 0; i < rUpdate.size(); i++)
            {
                mData(g,i) += rUpdate[i];
            }
        }
    }
 
 
    void WriteToCSVOutput(
        std::ofstream& rOutputStream,
        const Element& rElement,
        const PointerVector<StatisticsSampler>& rRecordedStatistics,
        const PointerVector<StatisticsSampler>& rHigherOrderStatistics,
        const std::size_t NumberOfMeasurements,
        const std::string& rSeparator) const
    {
        const Geometry<Node> &r_geometry = rElement.GetGeometry();
        const GeometryData::IntegrationMethod integration_method = rElement.GetIntegrationMethod();
        Matrix shape_functions;
        typename Geometry<Node>::ShapeFunctionsGradientsType shape_gradients;
        this->CalculateGeometryData(r_geometry, integration_method, shape_functions, shape_gradients);
 
        for (unsigned int g = 0; g < shape_functions.size1(); g++)
        {
            // Print element ID and integration point index
            rOutputStream << rElement.Id() << rSeparator << g << rSeparator;
 
            // Print integration point coordinates
            array_1d<double,3> coordinates(3,0.0);
            for (unsigned int n = 0; n < shape_functions.size2(); n++)
                coordinates += shape_functions(g,n) * r_geometry[n].Coordinates();
            rOutputStream << coordinates[0] << rSeparator << coordinates[1] << rSeparator << coordinates[2];
 
            auto data_iterator = IntegrationPointData(g).begin();
            for (auto it_sampler = rRecordedStatistics.begin(); it_sampler != rRecordedStatistics.end(); ++it_sampler)
            {
                it_sampler->OutputResult(rOutputStream,data_iterator,NumberOfMeasurements,rSeparator);
            }
 
            for (auto it_sampler = rHigherOrderStatistics.begin(); it_sampler != rHigherOrderStatistics.end(); ++it_sampler)
            {
                it_sampler->OutputResult(rOutputStream,data_iterator,NumberOfMeasurements,rSeparator);
            }
 
            rOutputStream << "\n";
        }
    }
 
 
    std::size_t NumberOfIntegrationPoints() const
    {
        return mData.size1();
    }
 
 
    std::size_t NumberOfStatisticalQuantities() const
    {
        return mData.size2();
    }
 
 
 
    IntegrationPointDataView IntegrationPointData(std::size_t IntegrationPointIndex)
    {
        KRATOS_DEBUG_ERROR_IF(IntegrationPointIndex >= mData.size1())
            << "Asking for integration point number " << IntegrationPointIndex
            << " but only " << mData.size1() << " points are recorded." << std::endl;
        return (mData.begin1() + IntegrationPointIndex);
    }
 
 
    IntegrationPointDataConstView IntegrationPointData(std::size_t IntegrationPointIndex) const
    {
        KRATOS_DEBUG_ERROR_IF(IntegrationPointIndex >= mData.size1())
            << "Asking for integration point number " << IntegrationPointIndex
            << " but only " << mData.size1() << " points are recorded." << std::endl;
        return (mData.begin1() + IntegrationPointIndex);
    }
 
 
    virtual std::string Info() const
    {
        std::stringstream buffer;
        buffer << "StatisticsData";
        return buffer.str();
    }
 
    virtual void PrintInfo(std::ostream &rOStream) const { rOStream << "StatisticsData"; }
 
    virtual void PrintData(std::ostream &rOStream) const {}
 
 
protected:
 
 
    void CalculateGeometryData(
        const Geometry<Node>& rGeometry,
        const GeometryData::IntegrationMethod IntegrationMethod,
        Matrix &rN,
        typename Geometry<Node>::ShapeFunctionsGradientsType &rDN_DX) const
    {
        const unsigned int number_of_nodes = rGeometry.PointsNumber();
        const unsigned int number_of_gauss_points = rGeometry.IntegrationPointsNumber(IntegrationMethod);
 
        Vector det_J;
        rGeometry.ShapeFunctionsIntegrationPointsGradients(rDN_DX,det_J,IntegrationMethod);
 
        rN.resize(number_of_gauss_points,number_of_nodes,false);
        rN = rGeometry.ShapeFunctionsValues(IntegrationMethod);
    }
 
 
private:
 
 
    Matrix mData;
 
 
    friend class Serializer;
 
    void save(Serializer& rSerializer) const {}
 
    void load(Serializer& rSerializer) {}
 
 
}; // Class StatisticsData
 
 
inline std::istream &operator>>(std::istream &rIStream,
                                StatisticsData/*< std::vector<double> >*/ &rThis)
{
    return rIStream;
}
 
inline std::ostream &operator<<(std::ostream &rOStream,
                                const StatisticsData/*< std::vector<double> >*/ &rThis)
{
    rThis.PrintInfo(rOStream);
    rOStream << std::endl;
    rThis.PrintData(rOStream);
 
    return rOStream;
}
 
 
} // namespace Kratos.
 
#endif // KRATOS_STATISTICS_DATA_H_INCLUDED  defined