macarray3.h

Go to the documentation of this file.

/*
**  macarray3.h
**
**  This is part of Shiokaze, a research-oriented fluid solver for computer graphics.
**  Created by Ryoichi Ando <rand@nii.ac.jp> on Feb 25, 2018.
**
**  Permission is hereby granted, free of charge, to any person obtaining a copy of
**  this software and associated documentation files (the "Software"), to deal in
**  the Software without restriction, including without limitation the rights to use,
**  copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the
**  Software, and to permit persons to whom the Software is furnished to do so,
**  subject to the following conditions:
**
**  The above copyright notice and this permission notice shall be included in all copies
**  or substantial portions of the Software.
**
**  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED,
**  INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A
**  PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
**  HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF
**  CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE
**  OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
*/
//
#ifndef SHKZ_MACARRAY3_H
#define SHKZ_MACARRAY3_H
//
#include <shiokaze/array/array3.h>
#include <algorithm>
#include <array>
//
SHKZ_BEGIN_NAMESPACE
//
template<class T> class macarray3 : public recursive_configurable, public messageable {
public:
    macarray3 ( recursive_configurable *parent, const shape3 &shape, vec3<T> value=vec3<T>(), std::string core_name="") : m_shape(shape), m_array_0(this,shape.face(0),value[0],core_name), m_array_1(this,shape.face(1),value[1],core_name), m_array_2(this,shape.face(2),value[2],core_name) {
        if( parent ) parent->add_child(this);
        else setup_now();
    }
    macarray3( recursive_configurable *parent, std::string core_name="" ) : macarray3(parent,shape3(0,0,0),vec3<T>(),core_name) {}
    macarray3 ( std::string core_name="") : macarray3(nullptr,shape3(0,0,0),vec3<T>(),core_name) {}
    macarray3( const shape3 &shape, vec3<T> value=vec3<T>(), std::string core_name="") : macarray3(nullptr,shape,value,core_name) {}
    macarray3 ( const macarray3& v ) : m_array_0(this), m_array_1(this), m_array_2(this){
        copy(v);
    }
    virtual bool send_message( std::string message, void *ptr ) override {
        bool handled (false);
        if( m_array_0.send_message(message,ptr)) handled = true;
        if( m_array_1.send_message(message,ptr)) handled = true;
        if( m_array_2.send_message(message,ptr)) handled = true;
        return handled;
    }
    virtual bool const_send_message( std::string message, void *ptr ) const override {
        bool handled (false);
        if( m_array_0.const_send_message(message,ptr)) handled = true;
        if( m_array_1.const_send_message(message,ptr)) handled = true;
        if( m_array_2.const_send_message(message,ptr)) handled = true;
        return handled;
    }
    macarray3& operator=(const macarray3 &array) {
        copy(array);
        return *this;
    }
    void copy( const macarray3 &array ) {
        if( this != &array ) {
            set_type(array.type());
            for( int dim : DIMS3 ) (*this)[dim].copy(array[dim]);
        }
    }
    void initialize ( const shape3 &shape, vec3<T> value=vec3<T>() ) {
        m_shape = shape;
        for( int dim : DIMS3 ) (*this)[dim].initialize(shape.face(dim),value[dim]);
    }
    void set_touch_only_actives( bool touch_only_actives ) {
        for( int dim : DIMS3 )(*this)[dim].set_touch_only_actives(touch_only_actives);
    }
    size_t count () const {
        size_t sum (0);
        for( int dim : DIMS3 ) sum += (*this)[dim].count();
        return sum;
    }
    std::array<std::vector<vec3i>,DIM3> actives() const {
        std::array<std::vector<vec3i>,DIM3> result;
        m_parallel.for_each( DIM3, [&]( size_t dim ) {
            result[dim] = (*this)[dim].actives();
        });
        return result;
    }
    void activate( const std::array<std::vector<vec3i>,DIM3> &active_entries, const std::array<vec3i,DIM3> &offsets={vec3i(),vec3i(),vec3i()}  ) {
        m_parallel.for_each( DIM3, [&]( size_t dim ) {
            (*this)[dim].activate(active_entries[dim],offsets[dim]);
        });
    }
    template <class Y> void activate_as( const macarray3<Y> &array, const std::array<vec3i,DIM3> &offsets={vec3i(),vec3i(),vec3i()} ) {
        m_parallel.for_each( DIM3, [&]( size_t dim ) {
            (*this)[dim].activate_as(array[dim],offsets[dim]);
        });
    }
    template <class Y> void activate_as_bit( const Y &array, const std::array<vec3i,DIM3> &offsets={vec3i(),vec3i()} ) {
        m_parallel.for_each( DIM3, [&]( size_t dim ) {
            (*this)[dim].activate_as_bit(array[dim],offsets[dim]);
        });
    }
    template <class Y> void activate_inside_as( const macarray3<Y> &array, const std::array<vec3i,DIM3> &offsets={vec3i(),vec3i(),vec3i()} ) {
        m_parallel.for_each( DIM3, [&]( size_t dim ) {
            (*this)[dim].activate_inside_as(array[dim],offsets[dim]);
        });
    }
    void activate_all() {
        m_parallel.for_each( DIM3, [&]( size_t dim ) {
            (*this)[dim].activate_all();
        });
    }
    void activate_inside() {
        m_parallel.for_each( DIM3, [&]( size_t dim ) {
            (*this)[dim].activate_inside();
        });
    }
    template <class Y> void copy_active_as( const macarray3<Y> &array, const vec3i &offset=vec3i() ) {
        m_parallel.for_each( DIM3, [&]( size_t dim ) {
            (*this)[dim].copy_active_as(array[dim],offset);
        });
    }
    vec3<T> get_background_value () const {
        return vec3<T>(
            m_array_0.get_background_value(),
            m_array_1.get_background_value(),
            m_array_2.get_background_value()
        );
    }
    void set_background_value( const vec3<T>& value ) {
        m_array_0.set_background_value(value[0]);
        m_array_1.set_background_value(value[1]);
        m_array_2.set_background_value(value[2]);
    }
    void convert_to_full( array3<vec3<T> > &cell_array ) const {
        std::vector<vec3i> active_cells;
        for( int dim : DIMS3 ) {
            std::vector<vec3i> active_faces = (*this)[dim].actives();
            for( const auto& fi : active_faces ) {
                vec3i pi (fi-vec3i(dim==0,dim==1,dim==2));
                if( ! m_shape.out_of_bounds(fi)) active_cells.push_back(fi);
                if( ! m_shape.out_of_bounds(pi)) active_cells.push_back(pi);
            }
        }
        std::sort( active_cells.begin(), active_cells.end(),[&](const vec3i &a, const vec3i &b) -> bool {
            return a[0]+m_shape.w*a[1]+(m_shape.w*m_shape.h)*a[2] < b[0]+m_shape.w*b[1]+(m_shape.w*m_shape.h)*b[2];
        });
        active_cells.erase( std::unique( active_cells.begin(), active_cells.end() ), active_cells.end() );
        cell_array.clear(get_background_value());
        cell_array.activate(active_cells);
        cell_array.parallel_actives([&](int i, int j, int k, auto &it, int tn) {
            vec3d v;
            char valid_count (0);
            for( unsigned dim : DIMS3 ) {
                char wsum (0);
                double value (0.0);
                if( (*this)[dim].active(i,j,k)) {
                    value += (*this)[dim](i,j,k);
                    wsum ++;
                }
                if( (*this)[dim].active(i+(dim==0),j+(dim==1),k+(dim==2))) {
                    value += (*this)[dim](i+(dim==0),j+(dim==1),k+(dim==2));
                    wsum ++;
                }
                if( wsum == 2 ) {
                    v[dim] = 0.5 * value;
                    valid_count ++;
                }
            }
            if( valid_count == 3 ) it.set(v);
            else it.set_off();
        });
    }
    void convert_to_full( macarray3<vec3<T> > &face_array ) const {
        //
        face_array.clear();
        face_array.activate_as(*this);
        face_array.parallel_actives([&](int dim, int i, int j, int k, auto &it, int tn) {
            //
            vec3d u;
            for( int u_dim : DIMS3 ) {
                if( u_dim == dim ) u[u_dim] = (*this)[u_dim](i,j,k);
                else {
                    const shape3 s = (*this)[u_dim].shape();
                    vec3i pi = vec3i(i,j,k)-vec3i(dim==0,dim==1,dim==2);
                    vec3i ivec = vec3i(dim==0,dim==1,dim==2);
                    vec3i jvec = vec3i(u_dim==0,u_dim==1,u_dim==2);
                    for( int ii=0; ii<2; ++ii ) for( int jj=0; jj<2; ++jj ) {
                        u[u_dim] += (*this)[u_dim](s.clamp(pi+ii*ivec+jj*jvec));
                    }
                    u[u_dim] /= 4.0;
                }
            }
            it.set(u);
        });
    }
    shape3 shape() const { return m_shape; }
    shape3 shape(int dim) const { return (*this)[dim].shape(); }
    void clear() {
        for( int dim : DIMS3 ) (*this)[dim].clear();
    }
    void clear(vec3<T> v) {
        for( int dim : DIMS3 ) (*this)[dim].clear(v[dim]);
    }
    bool operator!=( const macarray3<T> &v ) const {
        return ! (*this == v);
    }
    bool operator==(const macarray3<T> &v) const {
        for( int dim : DIMS3 ) {
            if ((*this)[dim] != v[dim]) return false;
        }
        return true;
    }
    void operator=(T v) {
        for( int dim : DIMS3 ) (*this)[dim] = v;
    }
    void operator+=(const macarray3<T> &v) {
        for( int dim : DIMS3 ) (*this)[dim] += v[dim];
    }
    void operator-=(const macarray3<T> &v) {
        for( int dim : DIMS3 ) (*this)[dim] -= v[dim];
    }
    void operator+=(const vec3<T> &v) {
        for( int dim : DIMS3 ) (*this)[dim] += v[dim];
    }
    void operator+=(T v) {
        for( int dim : DIMS3 ) (*this)[dim] += v;
    }
    void operator-=(const vec3<T> &v) {
        for( int dim : DIMS3 ) (*this)[dim] -= v[dim];
    }
    void operator-=(T v) {
        for( int dim : DIMS3 ) (*this)[dim] -= v;
    }
    void operator*=(T v) {
        for( int dim : DIMS3 ) (*this)[dim] *= v;
    }
    void operator/=(const T &v) {
        for( int dim : DIMS3 ) (*this)[dim] /= v;
    }
    const array3<T>& operator[](int dim) const {
        return dim==0 ? m_array_0 : (dim == 1 ? m_array_1 : m_array_2 );
    }
    array3<T>& operator[](int dim) {
        return dim==0 ? m_array_0 : (dim == 1 ? m_array_1 : m_array_2 );
    }
    void set_thread_num( int number ) {
        for( int dim : DIMS3 ) (*this)[dim].set_thread_num(number);
    }
    int get_thread_num() const {
        return m_array_0.get_thread_num();
    }
    //
    enum { ACTIVES = true, ALL = false };
    void parallel_actives( std::function<void(typename array3<T>::iterator& it)> func ) { parallel_op(func,ACTIVES); }
    void parallel_all( std::function<void(typename array3<T>::iterator& it)> func ) { parallel_op(func,ALL); }
    void parallel_op( std::function<void(typename array3<T>::iterator& it)> func, bool type=ALL ) {
        parallel_op([func](int dim, int i, int j, int k, typename array3<T>::iterator& it, int thread_index){
            func(it);
        },type);
    }
    void parallel_actives( std::function<void(int dim, int i, int j, int k, typename array3<T>::iterator& it)> func ) { parallel_op(func,ACTIVES); }
    void parallel_all( std::function<void(int dim, int i, int j, int k, typename array3<T>::iterator& it)> func ) { parallel_op(func,ALL); }
    void parallel_op( std::function<void(int dim, int i, int j, int k, typename array3<T>::iterator& it)> func, bool type=ALL ) {
        parallel_op([func](int dim, int i, int j, int k, typename array3<T>::iterator& it, int thread_index){
            func(dim,i,j,k,it);
        },type);
    }
    void parallel_actives( std::function<void(int dim, int i, int j, int k, typename array3<T>::iterator& it, int thread_index)> func ) { parallel_op(func,ACTIVES); }
    void parallel_all( std::function<void(int dim, int i, int j, int k, typename array3<T>::iterator& it, int thread_index)> func ) { parallel_op(func,ALL); }
    void parallel_op( std::function<void(int dim, int i, int j, int k, typename array3<T>::iterator& it, int thread_index)> func, bool type=ALL ) {
        for( int dim : DIMS3 ) {
            (*this)[dim].parallel_op([&](int i, int j, int k, typename array3<T>::iterator& it, int thread_index) {
                func(dim,i,j,k,it,thread_index);
            },type);
        };
    }
    void const_parallel_actives( std::function<void(const typename array3<T>::const_iterator& it)> func ) const { const_parallel_op(func,ACTIVES); }
    void const_parallel_all( std::function<void(const typename array3<T>::const_iterator& it)> func ) const { const_parallel_op(func,ALL); }
    void const_parallel_op( std::function<void(const typename array3<T>::const_iterator& it)> func, bool type=ALL ) const {
        const_parallel_op([func](int dim, int i, int j, int k, const typename array3<T>::const_iterator& it, int thread_index){
            func(it);
        },type);
    }
    void const_parallel_actives( std::function<void(int dim, int i, int j, int k, const typename array3<T>::const_iterator& it)> func ) const { const_parallel_op(func,ACTIVES); }
    void const_parallel_all( std::function<void(int dim, int i, int j, int k, const typename array3<T>::const_iterator& it)> func ) const { const_parallel_op(func,ALL); }
    void const_parallel_op( std::function<void(int dim, int i, int j, int k, const typename array3<T>::const_iterator& it)> func, bool type=ALL ) const {
        const_parallel_op([func](int dim, int i, int j, int k, const typename array3<T>::const_iterator& it, int thread_index){
            func(dim,i,j,k,it);
        },type);
    }
    void const_parallel_actives( std::function<void(int dim, int i, int j, int k, const typename array3<T>::const_iterator& it, int thread_index)> func ) const { const_parallel_op(func,ACTIVES); }
    void const_parallel_all( std::function<void(int dim, int i, int j, int k, const typename array3<T>::const_iterator& it, int thread_index)> func ) const { const_parallel_op(func,ALL); }
    void const_parallel_op( std::function<void(int dim, int i, int j, int k, const typename array3<T>::const_iterator& it, int thread_index)> func, bool type=ALL ) const {
        for( int dim : DIMS3 ) {
            (*this)[dim].const_parallel_op([&](int i, int j, int k, const typename array3<T>::const_iterator& it, int thread_index) {
                func(dim,i,j,k,it,thread_index);
            },type);
        };
    }
    void serial_actives( std::function<void(typename array3<T>::iterator& it)> func ) { serial_op(func,ACTIVES); }
    void serial_all( std::function<void(typename array3<T>::iterator& it)> func ) { serial_op(func,ALL); }
    void serial_op( std::function<void(typename array3<T>::iterator& it)> func, bool type=ALL ) {
        serial_op([func](int dim, int i, int j, int k, typename array3<T>::iterator& it) {
            func(it);
        },type);
    }
    void serial_actives( std::function<void(int dim, int i, int j, int k, typename array3<T>::iterator& it)> func ) { serial_op(func,ACTIVES); }
    void serial_all( std::function<void(int dim, int i, int j, int k, typename array3<T>::iterator& it)> func ) { serial_op(func,ALL); }
    void serial_op( std::function<void(int dim, int i, int j, int k, typename array3<T>::iterator& it)> func, bool type=ALL ) {
        for( int dim : DIMS3 ) {
            (*this)[dim].serial_op([&](int i, int j, int k, typename array3<T>::iterator& it) {
                func(dim,i,j,k,it);
            },type );
        }
    }
    void const_serial_actives( std::function<void(const typename array3<T>::const_iterator& it)> func ) const { const_serial_op(func,ACTIVES); }
    void const_serial_all( std::function<void(const typename array3<T>::const_iterator& it)> func ) const { const_serial_op(func,ALL); }
    void const_serial_op( std::function<void(const typename array3<T>::const_iterator& it)> func, bool type=ALL ) const {
        const_serial_op([func](int dim, int i, int j, int k, const typename array3<T>::const_iterator& it) {
            func(it);
        },type);
    }
    void const_serial_actives( std::function<void(int dim, int i, int j, int k, const typename array3<T>::const_iterator& it)> func ) const { const_serial_op(func,ACTIVES); }
    void const_serial_all( std::function<void(int dim, int i, int j, int k, const typename array3<T>::const_iterator& it)> func ) const { const_serial_op(func,ALL); }
    void const_serial_op( std::function<void(int dim, int i, int j, int k, const typename array3<T>::const_iterator& it)> func, bool type=ALL ) const {
        for( int dim : DIMS3 ) {
            (*this)[dim].const_serial_op([&](int i, int j, int k, const typename array3<T>::const_iterator& it) {
                func(dim,i,j,k,it);
            },type);
        }
    }
    void interruptible_serial_actives( std::function<bool(typename array3<T>::iterator& it)> func ) { serial_op(func,ACTIVES); }
    void interruptible_serial_all( std::function<bool(typename array3<T>::iterator& it)> func ) { serial_op(func,ALL); }
    void interruptible_serial_op( std::function<bool(typename array3<T>::iterator& it)> func, bool type=ALL ) {
        serial_op([func](int dim, int i, int j, int k, typename array3<T>::iterator& it) {
            func(it);
        },type);
    }
    void interruptible_serial_actives( std::function<bool(int dim, int i, int j, int k, typename array3<T>::iterator& it)> func ) { serial_op(func,ACTIVES); }
    void interruptible_serial_all( std::function<bool(int dim, int i, int j, int k, typename array3<T>::iterator& it)> func ) { serial_op(func,ALL); }
    void interruptible_serial_op( std::function<bool(int dim, int i, int j, int k, typename array3<T>::iterator& it)> func, bool type=ALL ) {
        for( int dim : DIMS3 ) {
            (*this)[dim].serial_op([&](int i, int j, int k, typename array3<T>::iterator& it) {
                func(dim,i,j,k,it);
            },type);
        }
    }
    void interruptible_const_serial_actives( std::function<bool(const typename array3<T>::const_iterator& it)> func ) const { const_serial_op(func,ACTIVES); }
    void interruptible_const_serial_all( std::function<bool(const typename array3<T>::const_iterator& it)> func ) const { const_serial_op(func,ALL); }
    void interruptible_const_serial_op( std::function<bool(const typename array3<T>::const_iterator& it)> func, bool type=ALL ) const {
        const_serial_op([func](int dim, int i, int j, int k, const typename array3<T>::const_iterator& it) {
            func(it);
        },type);
    }
    void interruptible_const_serial_actives( std::function<bool(int dim, int i, int j, int k, const typename array3<T>::const_iterator& it)> func ) const { const_serial_op(func,ACTIVES); }
    void interruptible_const_serial_all( std::function<bool(int dim, int i, int j, int k, const typename array3<T>::const_iterator& it)> func ) const { const_serial_op(func,ALL); }
    void interruptible_const_serial_op( std::function<bool(int dim, int i, int j, int k, const typename array3<T>::const_iterator& it)> func, bool type=ALL ) const {
        for( int dim : DIMS3 ) {
            (*this)[dim].const_serial_op([&](int i, int j, int k, const typename array3<T>::const_iterator& it) {
                func(dim,i,j,k,it);
            },type);
        }
    }
    void dilate( std::function<void(int dim, int i, int j, int k, typename array3<T>::iterator& it, int thread_index)> func, int count=1 ) {
        while( count -- ) {
            m_parallel.for_each(DIM3,[&]( size_t dim ) {
                operator[](dim).dilate([&](int i, int j, int k, typename array3<T>::iterator& it, int thread_index) {
                    func(dim,i,j,k,it,thread_index);
                });
            });
        }
    }
    void dilate( std::function<void(int dim, int i, int j, int k, typename array3<T>::iterator& it)> func, int count=1 ) {
        dilate([&](int dim, int i, int j, int k, typename array3<T>::iterator& it, int thread_index) {
            func(dim,i,j,k,it);
        },count);
    }
    void dilate( int count=1 ) {
        dilate([&](int dim, int i, int j, int k, typename array3<T>::iterator& it){ it.set(it()); },count);
    }
    void erode( std::function<bool(int dim, int i, int j, int k, const typename array3<T>::const_iterator& it, int thread_index)> func, int count=1 ) {
        while( count -- ) {
            m_parallel.for_each(DIM3,[&]( size_t dim ) {
                operator[](dim).erode([&](int i, int j, int k, const typename array3<T>::const_iterator& it, int thread_index) {
                    func(dim,i,j,it,thread_index);
                });
            });
        }
    }
    void erode( std::function<void(int dim, int i, int j, int k, const typename array3<T>::const_iterator& it)> func, int count=1 ) {
        erode([&](int dim, int i, int j, int k, const typename array3<T>::const_iterator& it, int thread_index) {
            func(dim,i,j,it);
        },count);
    }
    void erode( int count=1 ) {
        return erode([&](int dim, int i, int j, int k, const typename array3<T>::const_iterator& it){ it.set(it()); },count);
    }
    void set_core_name( std::string core_name ) {
        m_array_0.set_core_name(core_name);
        m_array_1.set_core_name(core_name);
        m_array_2.set_core_name(core_name);
    }
    std::string get_core_name() const {
        return m_array_0.get_core_name();
    }
    struct type3 {
        std::string core_name;
        shape3 shape;
        typename array3<T>::type3 type0;
        typename array3<T>::type3 type1;
        typename array3<T>::type3 type2;
        bool operator==( const type3 &type ) const {
            return
                core_name == type.core_name && shape == type.shape &&
                type0 == type.type0 && type1 == type.type1 && type2 == type.type2;
        }
    };
    type3 type() const {
        return { get_core_name(), m_shape, m_array_0.type(), m_array_1.type(), m_array_2.type() };
    }
    void set_type( const type3 &type ) {
        m_shape = type.shape;
        m_array_0.set_type(type.type0);
        m_array_1.set_type(type.type1);
        m_array_2.set_type(type.type2);
    }
private:
    parallel_driver m_parallel{this};
    array3<T> m_array_0;
    array3<T> m_array_1;
    array3<T> m_array_2;
    shape3 m_shape;
    //
};
//
template <class T> static inline macarray3<T> operator*(double s, const macarray3<T> &v) {
    return v*s;
}
//
SHKZ_END_NAMESPACE
//
#endif