mesh_with_data_packages.h

Go to the documentation of this file.

/* -----------------------------------------------------------------------------*
 *                               SPHinXsys                                      *
 * -----------------------------------------------------------------------------*
 * SPHinXsys (pronunciation: s'finksis) is an acronym from Smoothed Particle    *
 * Hydrodynamics for industrial compleX systems. It provides C++ APIs for       *
 * physical accurate simulation and aims to model coupled industrial dynamic    *
 * systems including fluid, solid, multi-body dynamics and beyond with SPH      *
 * (smoothed particle hydrodynamics), a meshless computational method using     *
 * particle discretization.                                                     *
 *                                                                              *
 * SPHinXsys is partially funded by German Research Foundation                  *
 * (Deutsche Forschungsgemeinschaft) DFG HU1527/6-1, HU1527/10-1,               *
 * HU1527/12-1 and HU1527/12-4.                                                 *
 *                                                                              *
 * Portions copyright (c) 2017-2022 Technical University of Munich and          *
 * the authors' affiliations.                                                   *
 *                                                                              *
 * Licensed under the Apache License, Version 2.0 (the "License"); you may      *
 * not use this file except in compliance with the License. You may obtain a    *
 * copy of the License at http://www.apache.org/licenses/LICENSE-2.0.           *
 *                                                                              *
 * -----------------------------------------------------------------------------*/
#ifndef MESH_WITH_DATA_PACKAGES_H
#define MESH_WITH_DATA_PACKAGES_H

#include "base_mesh.h"
#include "my_memory_pool.h"

#include <fstream>
#include <algorithm>
#include <mutex>
#include <functional>
using namespace std::placeholders;

namespace SPH
{
    class Kernel;

    template <class ReturnType, class DataPackageType>
    using PackageFunctor = std::function<ReturnType(DataPackageType *, Real)>;
    template <class DataPackageType>
    void PackageIterator(ConcurrentVector<DataPackageType *> &data_pkgs,
                         PackageFunctor<void, DataPackageType> &pkg_functor, Real dt = 0.0)
    {
        for (size_t i = 0; i != data_pkgs.size(); ++i)
            pkg_functor(data_pkgs[i], dt);
    };
    template <class DataPackageType>
    void PackageIterator_parallel(ConcurrentVector<DataPackageType *> &data_pkgs,
                                  PackageFunctor<void, DataPackageType> &pkg_functor, Real dt = 0.0)
    {
        parallel_for(
            blocked_range<size_t>(0, data_pkgs.size()),
            [&](const blocked_range<size_t> &r)
            {
                for (size_t i = r.begin(); i != r.end(); ++i)
                {
                    pkg_functor(data_pkgs[i], dt);
                }
            },
            ap);
    };
    template <class ReturnType, typename ReduceOperation, class DataPackageType>
    ReturnType ReducePackageIterator(ConcurrentVector<DataPackageType *> &data_pkgs, ReturnType temp,
                                     PackageFunctor<ReturnType, DataPackageType> &reduce_pkg_functor,
                                     ReduceOperation &reduce_operation, Real dt = 0.0)
    {
        for (size_t i = 0; i < data_pkgs.size(); ++i)
        {
            temp = reduce_operation(temp, reduce_pkg_functor(data_pkgs[i], dt));
        }
        return temp;
    };
    template <class ReturnType, typename ReduceOperation, class DataPackageType>
    ReturnType ReducePackageIterator_parallel(ConcurrentVector<DataPackageType *> &data_pkgs, ReturnType temp,
                                              PackageFunctor<ReturnType, DataPackageType> &reduce_pkg_functor,
                                              ReduceOperation &reduce_operation, Real dt = 0.0)
    {
        return parallel_reduce(
            blocked_range<size_t>(0, data_pkgs.size()),
            temp, [&](const blocked_range<size_t> &r, ReturnType temp0) -> ReturnType
            {
                for (size_t i = r.begin(); i != r.end(); ++i)
                {
                    temp0 = reduce_operation(temp, reduce_pkg_functor(data_pkgs[i], dt));
                }
                return temp0; },
            [&](ReturnType x, ReturnType y) -> ReturnType
            {
                return reduce_operation(x, y);
            });
    };

    class BaseDataPackage
    {
    public:
        Vecu pkg_index_;    
        bool is_inner_pkg_; 
        BaseDataPackage() : pkg_index_(0), is_inner_pkg_(false){};
        virtual ~BaseDataPackage(){};
    };

    template <int PKG_SIZE, int ADDRS_SIZE>
    class GridDataPackage : public BaseDataPackage, public BaseMesh
    {
    public:
        template <typename DataType>
        using PackageData = PackageDataMatrix<DataType, PKG_SIZE>;
        template <typename DataType>
        using PackageDataAddress = PackageDataMatrix<DataType *, ADDRS_SIZE>;
        GeneralDataAssemble<PackageData> all_pkg_data_;
        GeneralDataAssemble<PackageDataAddress> all_pkg_data_addrs_;

        template <typename DataType>
        using PackageTemporaryData = PackageDataMatrix<DataType, ADDRS_SIZE>;

        GridDataPackage() : BaseDataPackage(), BaseMesh(Vecu(ADDRS_SIZE)){};
        virtual ~GridDataPackage(){};

        constexpr int PackageSize() { return PKG_SIZE; };
        constexpr int AddressSize() { return ADDRS_SIZE; };
        constexpr int AddressBufferWidth() { return (ADDRS_SIZE - PKG_SIZE) / 2; };
        constexpr int OperationUpperBound() { return PKG_SIZE + AddressBufferWidth(); };
        Vecd DataLowerBound() { return mesh_lower_bound_ + Vecd(grid_spacing_) * (Real)AddressBufferWidth(); };
        void initializePackageGeometry(const Vecd &pkg_lower_bound, Real data_spacing)
        {
            mesh_lower_bound_ = pkg_lower_bound - Vecd(data_spacing) * ((Real)AddressBufferWidth() - 0.5);
            grid_spacing_ = data_spacing;
        };
        template <typename DataType>
        DataType probeDataPackage(PackageDataAddress<DataType> &pkg_data_addrs, const Vecd &position);
        template <typename InDataType, typename OutDataType>
        void computeGradient(PackageDataAddress<InDataType> &in_pkg_data_addrs,
                             PackageDataAddress<OutDataType> out_pkg_data_addrs, Real dt = 0.0);
        template <typename InDataType, typename OutDataType>
        void computeNormalizedGradient(PackageDataAddress<InDataType> &in_pkg_data_addrs,
                                       PackageDataAddress<OutDataType> out_pkg_data_addrs, Real dt = 0.0);

    protected:
        template <typename DataType>
        void registerPackageData(PackageData<DataType> &pkg_data,
                                 PackageDataAddress<DataType> &pkg_data_addrs)
        {
            constexpr int type_index = DataTypeIndex<DataType>::value;
            std::get<type_index>(all_pkg_data_).push_back(&pkg_data);
            std::get<type_index>(all_pkg_data_addrs_).push_back(&pkg_data_addrs);
        };
        template <typename DataType>
        struct initializePackageDataAddress
        {
            void operator()(GeneralDataAssemble<PackageData> &all_pkg_data,
                            GeneralDataAssemble<PackageDataAddress> &all_pkg_data_addrs);
        };
        DataAssembleOperation<initializePackageDataAddress> initialize_pkg_data_addrs_;
        template <typename DataType>
        struct assignPackageDataAddress
        {
            void operator()(GeneralDataAssemble<PackageDataAddress> &all_pkg_data_addrs,
                            const Vecu &addrs_index,
                            GeneralDataAssemble<PackageData> &all_pkg_data,
                            const Vecu &data_index);
        };
        DataAssembleOperation<assignPackageDataAddress> assign_pkg_data_addrs_;

        template <typename DataType>
        DataType CornerAverage(PackageDataAddress<DataType> &pkg_data_addrs, Veci addrs_index, Veci corner_direction);

    public:
        void initializeSingularDataAddress()
        {
            initialize_pkg_data_addrs_(all_pkg_data_, all_pkg_data_addrs_);
        };

        void assignAllPackageDataAddress(const Vecu &addrs_index, GridDataPackage *src_pkg, const Vecu &data_index)
        {
            assign_pkg_data_addrs_(all_pkg_data_addrs_, addrs_index, src_pkg->all_pkg_data_, data_index);
        };
    };

    template <class MeshFieldType, class GridDataPackageType>
    class MeshWithGridDataPackages : public MeshFieldType, public Mesh
    {
    public:
        MyMemoryPool<GridDataPackageType> data_pkg_pool_;         
        MeshDataMatrix<GridDataPackageType *> data_pkg_addrs_;    
        ConcurrentVector<GridDataPackageType *> inner_data_pkgs_; 
        template <typename... Args>
        explicit MeshWithGridDataPackages(BoundingBox tentative_bounds, Real data_spacing, size_t buffer_size, Args &&...args)
            : MeshFieldType(std::forward<Args>(args)...),
              Mesh(tentative_bounds, GridDataPackageType().PackageSize() * data_spacing, buffer_size),
              data_spacing_(data_spacing),
              global_mesh_(this->mesh_lower_bound_ + Vecd(data_spacing) * 0.5, data_spacing, this->number_of_cells_ * pkg_size_)
        {
            allocateMeshDataMatrix();
        };
        virtual ~MeshWithGridDataPackages() { deleteMeshDataMatrix(); };

        void allocateMeshDataMatrix(); 
        void deleteMeshDataMatrix();   
        template <class DataType, typename PackageDataAddressType, PackageDataAddressType GridDataPackageType::*MemPtr>
        DataType probeMesh(const Vecd &position);
        virtual Real DataSpacing() override { return data_spacing_; };

    protected:
        Real data_spacing_;                                                       
        const int pkg_size_ = GridDataPackageType().PackageSize();                
        const int pkg_addrs_buffer_ = GridDataPackageType().AddressBufferWidth(); 
        const int pkg_operations_ = pkg_size_ + pkg_addrs_buffer_;                
        const int pkg_addrs_size_ = pkg_size_ + 2 * pkg_addrs_buffer_;            
        std::mutex mutex_my_pool;                                                 
        BaseMesh global_mesh_;                                                    
        StdVec<GridDataPackageType *> singular_data_pkgs_addrs_;

        template <typename... ConstructorArgs>
        void initializeASingularDataPackage(ConstructorArgs &&...args)
        {
            GridDataPackageType *new_data_pkg = data_pkg_pool_.malloc();
            new_data_pkg->registerAllVariables();
            new_data_pkg->initializeSingularData(std::forward<ConstructorArgs>(args)...);
            new_data_pkg->initializeSingularDataAddress();
            singular_data_pkgs_addrs_.push_back(new_data_pkg);
        };

        void assignDataPackageAddress(const Vecu &cell_index, GridDataPackageType *data_pkg);
        GridDataPackageType *DataPackageFromCellIndex(const Vecu &cell_index);
        template <typename DataType, typename PackageDataType, PackageDataType GridDataPackageType::*MemPtr>
        DataType DataValueFromGlobalIndex(const Vecu &global_grid_index);
        void initializePackageAddressesInACell(const Vecu &cell_index);
        std::pair<int, int> CellShiftAndDataIndex(int data_addrs_index_component)
        {
            std::pair<int, int> shift_and_index;
            int signed_date_index = data_addrs_index_component - pkg_addrs_buffer_;
            shift_and_index.first = (signed_date_index + pkg_size_) / pkg_size_ - 1;
            shift_and_index.second = signed_date_index - shift_and_index.first * pkg_size_;
            return shift_and_index;
        }
    };
}
#endif // MESH_WITH_DATA_PACKAGES_H