field_model.hh

Go to the documentation of this file.

/*!
 *
 * Copyright (C) 2015 Technical University of Liberec.  All rights reserved.
 *
 * This program is free software; you can redistribute it and/or modify it under
 * the terms of the GNU General Public License version 3 as published by the
 * Free Software Foundation. (http://www.gnu.org/licenses/gpl-3.0.en.html)
 *
 * This program is distributed in the hope that it will be useful, but WITHOUT
 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
 * FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
 *
 *
 * @file    field_model.hh
 * @brief
 */
 
#ifndef FIELD_MODEL_HH_
#define FIELD_MODEL_HH_
 
#include <armadillo>
#include <iostream>
#include <tuple>
#include <string>
#include <vector>
#include <utility>
#include <type_traits>
 
#include "fields/field.hh"
#include "fields/field_common.hh"
#include "fields/field_algo_base.hh"
#include "fields/field_values.hh"
#include "fields/field_value_cache.hh"
#include "fields/field.hh"
#include "fields/multi_field.hh"
#include "system/fmt/posix.h"                          // for FMT_UNUSED
 
template <int spacedim> class ElementAccessor;
 
 
/**
 * Wrapper for resolution of the overloaded functions passed as a parameter to the FieldModel.
 *
 * Example:
 @code
   double product(double x, double y) {...}
   Vec product(double x, Vec y) {...}
   Model<3, FieldValue<3>::VectorFixed>::create(wrapper_overload(product), f_scal, f_vec);
 @endcode
 *
 * Should automaticaly resolve the second model_cache_item::eval function.
 */
#define wrap_overload(func) [](auto&&... ps){ return func( std::forward<decltype(ps)>(ps)... ); }
 
namespace detail
{
    /**
     * Extract cached values on index i_cache from the input fields and call given function on these values.
     */
 
    /// base case for building up arguments for the function call
    template< typename CALLABLE, typename FIELD_TUPLE, int INDEX >
    struct model_cache_item
    {
        template< typename... Vs >
        static auto eval(int i_cache, CALLABLE f, FIELD_TUPLE fields, Vs&&... args) -> decltype(auto) {
            const auto &single_field = std::get < INDEX - 1 > (std::forward<decltype(fields)>(fields));
            return model_cache_item<CALLABLE, FIELD_TUPLE, INDEX - 1>::eval(
                    i_cache, f, std::forward<decltype(fields)>(fields),
                    single_field[i_cache], std::forward<Vs>(args)...);
 
        }
    };
 
    /// terminal case - do the actual function call
    template< typename CALLABLE, typename FIELD_TUPLE >
    struct model_cache_item< CALLABLE, FIELD_TUPLE, 0 >
    {
        template< typename... Vs >
        static auto eval(FMT_UNUSED int i_cache, CALLABLE f, FMT_UNUSED FIELD_TUPLE fields, Vs&&... args) -> decltype(auto)
        {
            return f(std::forward<Vs>(args)...);
        }
    };
 
    /**
     * Check common number of components of the input fields/multifields.
     * Return number of components.
     * Return 0 for no multifields.
     * Throw for different number of components.
     */
    template<typename FIELD_TUPLE, int INDEX >
    struct n_components {
        static uint eval(FIELD_TUPLE fields, uint n_comp) {
            const auto &single_field = std::get < INDEX - 1 > (std::forward<decltype(fields)>(fields));
            uint n_comp_new = single_field.n_comp();
            if (n_comp == 0) {
                n_comp = n_comp_new;
            } else {
                ASSERT_DBG(n_comp == n_comp_new);
            }
            return n_components<FIELD_TUPLE, INDEX - 1>::eval(std::forward<decltype(fields)>(fields), n_comp);
        };
    };
 
    template<typename FIELD_TUPLE>
    struct n_components<FIELD_TUPLE, 0> {
        static uint eval(FMT_UNUSED FIELD_TUPLE fields, uint n_comp)
        {
            return n_comp;
        };
    };
 
 
    /**
     * Return component 'i_comp' of the multifield 'f' or the field 'f'.
     */
    /*template<class FIELD>
    auto field_component(FIELD f, uint i_comp) -> decltype(auto)
    {
        if (f.is_multifield()) {
            return f[i_comp];
        } else {
            return f;
        }
    }*/
 
    /**
     * Return component 'i_comp' of the multifield 'f'.
     */
    template<int spacedim, class Value>
    auto field_component(const MultiField<spacedim, Value> &f, uint i_comp) -> decltype(auto)
    {
        ASSERT(f.is_multifield());
        return f[i_comp];
    }
 
    /**
     * Return the field 'f'. Variant to previous method.
     */
    template<int spacedim, class Value>
    auto field_component(const Field<spacedim, Value> &f, FMT_UNUSED uint i_comp) -> decltype(auto)
    {
        ASSERT(!f.is_multifield());
        return f;
    }
 
    /**
     * For given tuple of fields/multifields and given component index
     * return  the tuple of the selected components.
     */
    template<typename FIELD_TUPLE, int INDEX>
    struct get_components {
        static auto eval(FIELD_TUPLE fields, uint i_comp) -> decltype(auto)
        {
            const auto &single_field = std::get < INDEX - 1 > (std::forward<decltype(fields)>(fields));
            return std::tuple_cat(
                    get_components<FIELD_TUPLE, INDEX - 1>::eval(
                            std::forward<decltype(fields)>(fields), i_comp),
                    std::forward_as_tuple(field_component(single_field, i_comp))
                    );
        };
    };
 
    template<typename FIELD_TUPLE>
    struct get_components<FIELD_TUPLE, 0> {
        static auto eval(FMT_UNUSED FIELD_TUPLE fields, FMT_UNUSED uint n_comp) -> decltype(auto)
        {
            return std::forward_as_tuple<>();
        };
    };
 
    /**
     * Check common number of components of the input fields/multifields.
     * Return number of components.
     * Return 0 for no multifields.
     * Throw for different number of components.
     */
    template<typename FIELD_TUPLE, int INDEX >
    struct get_dependency {
        static std::vector<const FieldCommon *> eval(FIELD_TUPLE fields) {
            const auto &single_field = std::get < INDEX - 1 > (std::forward<decltype(fields)>(fields));
            auto vec = get_dependency<FIELD_TUPLE, INDEX - 1>::eval(std::forward<decltype(fields)>(fields));
            vec.push_back(&single_field);
            return vec;
        };
    };
 
    template<typename FIELD_TUPLE>
    struct get_dependency<FIELD_TUPLE, 0> {
        static std::vector<const FieldCommon *> eval(FMT_UNUSED FIELD_TUPLE fields)
        {
            return std::vector<const FieldCommon *>();
        };
    };
 
 
 
    template< int spacedim, class Value, typename CALLABLE, typename FIELD_TUPLE, int INDEX >
    struct field_value
    {
        typedef typename FieldAlgorithmBase<spacedim, Value>::Point Point;
 
        template< typename... Vs >
        static auto eval(const Point &p, const ElementAccessor<spacedim> &elm, CALLABLE f, FIELD_TUPLE fields, Vs&&... args) -> decltype(auto) {
            const auto &single_field = std::get < INDEX - 1 > (std::forward<decltype(fields)>(fields));
            return field_value<spacedim, Value, CALLABLE, FIELD_TUPLE, INDEX - 1>::eval(
                    p, elm, f, std::forward<decltype(fields)>(fields),
                    single_field.value(p, elm), std::forward<Vs>(args)...);
        }
    };
 
 
    /// terminal case - do the actual function call
    template< int spacedim, class Value, typename CALLABLE, typename FIELD_TUPLE >
    struct field_value< spacedim, Value, CALLABLE, FIELD_TUPLE, 0 >
    {
        typedef typename FieldAlgorithmBase<spacedim, Value>::Point Point;
 
        template< typename... Vs >
        static auto eval(FMT_UNUSED const Point &p, FMT_UNUSED const ElementAccessor<spacedim> &elm, CALLABLE f, FMT_UNUSED FIELD_TUPLE fields,
                Vs&&... args) -> decltype(auto)
        {
            return f(std::forward<Vs>(args)...);
        }
    };
 
 
    template<typename Function, typename Tuple>
    auto call(Function f, Tuple t)
    {
    }
}
 
 
/**
 * Class representing field computing form results of other fields.
 *
 * Example of usage:
   @code
    // Functor class with defined operator ()
    class FnProduct {
    public:
        Vector operator() (Scalar a, Vector v) {
             return a(0) * v;
        }
    };
 
    ...
 
    // Create ElementCacheMap
    std::shared_ptr<EvalPoints> eval_points = std::make_shared<EvalPoints>();
    eval_poinzs->add_bulk<3>( quad ); // initialize EvalPoints
    ElementCacheMap elm_cache_map;
    elm_cache_map.init(eval_points);
 
    // Definition of fields
    Field<3, FieldValue<3>::Scalar > f_scal;
    Field<3, FieldValue<3>::VectorFixed > f_vec;
    Field<3, FieldValue<3>::VectorFixed > result;
    ... // fill data to fields f_scal, f_vec
 
    // create instance FieldModel class, use helper method Model::create to simply passsing of parameters
    auto f_product = Model<3, FieldValue<3>::VectorFixed>::create(FnProduct(), f_scal, f_vec);
    // set field on all regions
    result.set_mesh( *mesh );
    result.set(f_product, time);
    result.cache_reallocate(elm_cache_map);
    result.set_time(tg.step(), LimitSide::right);
 
    // cache_update
    result.cache_update(elm_cache_map);
   @endcode
 *
 */
template<int spacedim, class Value, typename Fn, class ... InputFields>
class FieldModel : public FieldAlgorithmBase<spacedim, Value>
{
private:
    Fn fn;
    typedef std::tuple<InputFields...> FieldsTuple;
    FieldsTuple input_fields;
 
public:
    typedef typename FieldAlgorithmBase<spacedim, Value>::Point Point;
 
    FieldModel(Fn func, InputFields... args)
    : fn(func), input_fields( std::forward_as_tuple((args)...) )
    {}
 
    /// Implements FieldAlgoBase::set_dependency
    std::vector<const FieldCommon *> set_dependency(FMT_UNUSED FieldSet &field_set) {
        return detail::get_dependency<
                            decltype(input_fields),
                            std::tuple_size<FieldsTuple>::value
                        >::eval(input_fields);
    }
 
    /// Implements FieldAlgoBase::cache_update
    void cache_update(FieldValueCache<typename Value::element_type> &data_cache,
                ElementCacheMap &cache_map, unsigned int region_patch_idx) override {
        unsigned int reg_chunk_begin = cache_map.region_chunk_begin(region_patch_idx);
        unsigned int reg_chunk_end = cache_map.region_chunk_end(region_patch_idx);
        for(unsigned int i_cache=reg_chunk_begin; i_cache<reg_chunk_end; ++i_cache) {
            data_cache.set(i_cache) =
                detail::model_cache_item<
                    Fn,
                    decltype(input_fields),
                    std::tuple_size<FieldsTuple>::value
                >::eval(i_cache, fn, input_fields);
        }
    }
 
    /// Implementation of virtual method
    typename Value::return_type const &value(const Point &p, const ElementAccessor<spacedim> &elm) override {
        this->r_value_ = detail::field_value<
                spacedim, Value, Fn, decltype(input_fields), std::tuple_size<FieldsTuple>::value
                >::eval(p, elm, fn, input_fields);
        return this->r_value_;
    }
 
    /// Implementation of virtual method
    void value_list(const Armor::array &point_list, const ElementAccessor<spacedim> &elm,
                std::vector<typename Value::return_type> &value_list) override {
        ASSERT_EQ_DBG(point_list.size(), value_list.size()).error("Different size of point list and value list.\n");
 
        for (uint i=0; i<point_list.size(); ++i)
            value_list[i] = this->value(point_list.template mat<Value::NRows_, Value::NCols_>(i), elm);
    }
 
};
 
 
/**
 * Auxiliary class to avoid explicit specification of constructor template parameters.
 */
template<int spacedim, class Value>
class Model {
public:
    typedef FieldAlgorithmBase<spacedim, Value> FieldBaseType;
    typedef std::shared_ptr< FieldBaseType > FieldBasePtr;
 
    /**
     * Fn is a functor class and fn its instance.
     */
    template<typename Fn, class ... InputFields>
    //static auto create(Fn *fn, InputFields&&... inputs) -> decltype(auto)
    static auto create(Fn fn,  InputFields&&... inputs) -> decltype(auto)
    {
        return std::make_shared<FieldModel<spacedim, Value, Fn, InputFields...>>(fn, std::forward<InputFields>(inputs)...);
    }
 
 
 
    template<typename Function, typename Tuple, size_t ... I>
    static auto call_create(Function f, Tuple t, std::index_sequence<I ...>)
    {
        return create(f, std::get<I>(t) ...);
    }
 
    /**
     * Fn is a functor class and fn its instance.
     */
    template<typename Fn, class ... InputFields>
    //static auto create_multi(Fn *fn, InputFields&&... inputs) -> decltype(auto)
    static auto create_multi(Fn fn,  InputFields&&... inputs) -> decltype(auto)
    {
        typedef std::tuple<InputFields...> FieldTuple;
        FieldTuple field_tuple = std::forward_as_tuple((inputs)...);
        constexpr uint n_inputs = sizeof...(InputFields);
        uint n_comp = detail::n_components< FieldTuple, n_inputs>::eval(field_tuple, 0);
        ASSERT_DBG(n_comp > 0);
        std::vector<FieldBasePtr> result_components;
        for(uint i=0; i<n_comp; i++) {
            const auto & component_of_inputs = detail::get_components< FieldTuple, n_inputs>::eval(field_tuple, i);
            //const auto & all_args = std::tuple_cat(std::make_tuple(fn), component_of_inputs);
            //FieldBasePtr component_field = detail::call(create<Fn, InputFields&&...>, all_args);
 
            FieldBasePtr component_field = call_create(fn, component_of_inputs, std::make_index_sequence<n_inputs>{});
            result_components.push_back(component_field);
        }
 
        return result_components;
    }
};
 
 
 
 
 
 
 
 
 
 
#endif /* FIELD_MODEL_HH_ */