Connectivity.h

Go to the documentation of this file.

/*
 * (C) Copyright 2013 ECMWF.
 *
 * This software is licensed under the terms of the Apache Licence Version 2.0
 * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
 * In applying this licence, ECMWF does not waive the privileges and immunities
 * granted to it by virtue of its status as an intergovernmental organisation
 * nor does it submit to any jurisdiction.
 */


#pragma once

#include <array>
#include <cstring>
#include <initializer_list>

#include "atlas/array.h"
#include "atlas/array/ArrayView.h"
#include "atlas/array/DataType.h"
#include "atlas/array/IndexView.h"
#include "atlas/array/SVector.h"
#include "atlas/array/Vector.h"
#include "atlas/array_fwd.h"
#include "atlas/library/config.h"
#include "atlas/util/Object.h"

namespace eckit {
class Stream;
}

namespace atlas {
namespace mesh {

constexpr size_t MAX_STRING_SIZE() {
    return 60;
}

template <typename ConnectivityImpl>
class ConnectivityInterface : public util::Object, DOXYGEN_HIDE( public ConnectivityImpl ) {
    using ConnectivityImpl::ConnectivityImpl;
    using util::Object::Object;
};

// Classes defined in this file:
class IrregularConnectivityImpl;
class BlockConnectivityImpl;
class MultiBlockConnectivityImpl;

// --------------------------------------------------------------------------

#ifndef DOXYGEN_SHOULD_SKIP_THIS
#if ATLAS_HAVE_FORTRAN
#define INDEX_DEREF Index
#define INDEX_REF Index
#define FROM_FORTRAN -1
#define TO_FORTRAN +1
#else
#define INDEX_DEREF *
#define INDEX_REF idx_t&
#define FROM_FORTRAN
#define TO_FORTRAN
#endif
#endif


namespace detail {
// FortranIndex:
// Helper class that does +1 and -1 operations on stored values

class ConnectivityIndex {
public:
    enum
    {
        BASE = 1
    };

public:
    ATLAS_HOST_DEVICE ConnectivityIndex( const ConnectivityIndex& other ) { set( other.get() ); }
    ATLAS_HOST_DEVICE ConnectivityIndex( idx_t* idx ) : idx_( idx ) {}
    ATLAS_HOST_DEVICE void set( const idx_t& value ) { *( idx_ ) = value + BASE; }
    ATLAS_HOST_DEVICE idx_t get() const { return *(idx_)-BASE; }
    ATLAS_HOST_DEVICE void operator=( const idx_t& value ) { set( value ); }
    ATLAS_HOST_DEVICE ConnectivityIndex& operator+( const idx_t& value ) {
        *( idx_ ) += value;
        return *this;
    }
    ATLAS_HOST_DEVICE ConnectivityIndex& operator-( const idx_t& value ) {
        *( idx_ ) -= value;
        return *this;
    }
    ATLAS_HOST_DEVICE ConnectivityIndex& operator--() {
        --( *( idx_ ) );
        return *this;
    }
    ATLAS_HOST_DEVICE ConnectivityIndex& operator++() {
        ++( *( idx_ ) );
        return *this;
    }

    // implicit conversion
    ATLAS_HOST_DEVICE operator idx_t() const { return get(); }

private:
    idx_t* idx_;
};
}  // namespace detail

class ConnectivityRow {
#if ATLAS_HAVE_FORTRAN
    using Index = detail::ConnectivityIndex;
#else
    using Index = idx_t&;
#endif
public:
    ATLAS_HOST_DEVICE
    ConnectivityRow( idx_t* data, idx_t size ) : data_( data ), size_( size ) {}

    template <typename Int>
    ATLAS_HOST_DEVICE idx_t operator()( Int i ) const {
        return data_[i] FROM_FORTRAN;
    }

    template <typename Int>
    ATLAS_HOST_DEVICE INDEX_REF operator()( Int i ) {
        return INDEX_DEREF( data_ + i );
    }

    ATLAS_HOST_DEVICE
    idx_t size() const { return size_; }

private:
    idx_t* data_;
    idx_t size_;
};

class IrregularConnectivityImpl {
public:
    typedef ConnectivityRow Row;

public:
    //-- Constructors

    IrregularConnectivityImpl( const std::string& name = "" );

    IrregularConnectivityImpl( idx_t values[], idx_t rows, idx_t displs[], idx_t counts[] );

#ifdef __CUDACC__
    // This ctr has to be defined in the header, since __CUDACC__ will identify
    // whether
    // it is compiled it for a GPU kernel
    IrregularConnectivityImpl( const IrregularConnectivityImpl& other ) :
        owns_( false ),
        values_( other.values_ ),
        displs_( other.displs_ ),
        counts_( other.counts_ ),
        missing_value_( other.missing_value_ ),
        rows_( other.rows_ ),
        maxcols_( other.maxcols_ ),
        mincols_( other.mincols_ ),
        ctxt_( nullptr ) {}
#endif

    explicit IrregularConnectivityImpl( eckit::Stream& );

    virtual ~IrregularConnectivityImpl();

    //-- Accessors

    const std::string name() const { return std::string( name_ ); }

    void rename( const std::string& name ) {
        std::strncpy( name_, name.c_str(), std::max( name.size(), MAX_STRING_SIZE() ) );
    }

    ATLAS_HOST_DEVICE
    idx_t rows() const { return rows_; }

    ATLAS_HOST_DEVICE
    idx_t cols( idx_t row_idx ) const { return counts_[row_idx]; }

    ATLAS_HOST_DEVICE
    idx_t maxcols() const { return maxcols_; }

    ATLAS_HOST_DEVICE
    idx_t mincols() const { return mincols_; }

    ATLAS_HOST_DEVICE
    idx_t operator()( idx_t row_idx, idx_t col_idx ) const;

    idx_t size() const { return values_.size(); }

    ATLAS_HOST_DEVICE
    idx_t missing_value() const { return missing_value_; }

    ATLAS_HOST_DEVICE
    Row row( idx_t row_idx ) const;

    // -- Modifiers

    void set( idx_t row_idx, const idx_t column_values[] );

    void set( idx_t row_idx, idx_t col_idx, const idx_t value );

    virtual void resize( idx_t old_size, idx_t size, bool initialize, const idx_t values[], bool fortran_array );

    virtual void add( idx_t rows, idx_t cols, const idx_t values[], bool fortran_array = false );

    virtual void add( idx_t rows, idx_t cols );

    virtual void add( idx_t rows, const idx_t cols[] );

    virtual void add( const BlockConnectivityImpl& block );

    virtual void insert( idx_t position, idx_t rows, idx_t cols, const idx_t values[], bool fortran_array = false );

    virtual void insert( idx_t position, idx_t rows, idx_t cols );

    virtual void insert( idx_t position, idx_t rows, const idx_t cols[] );

    virtual void clear();

    virtual size_t footprint() const;

    idx_t displs( const idx_t row ) const { return displs_[row]; }

    void dump( std::ostream& os ) const;
    friend std::ostream& operator<<( std::ostream& os, const IrregularConnectivityImpl& p ) {
        p.dump( os );
        return os;
    }

    virtual void encode( eckit::Stream& s ) const { encode_( s ); }
    virtual void decode( eckit::Stream& s ) { decode_( s ); }

protected:
    bool owns() { return owns_; }
    const idx_t* displs() const { return displs_.data(); }
    const idx_t* counts() const { return counts_.data(); }

    void encode_( eckit::Stream& ) const;

    void decode_( eckit::Stream& );

    friend eckit::Stream& operator<<( eckit::Stream& s, const IrregularConnectivityImpl& x ) {
        x.encode( s );
        return s;
    }

    friend eckit::Stream& operator>>( eckit::Stream& s, IrregularConnectivityImpl& x ) {
        x.decode( s );
        return s;
    }

private:
    void on_delete();
    void on_update();

private:
    char name_[MAX_STRING_SIZE()];
    bool owns_;

protected:
    array::SVector<idx_t> values_;
    array::SVector<idx_t> displs_;
    array::SVector<idx_t> counts_;

    idx_t missing_value_;
    idx_t rows_;
    idx_t maxcols_;
    idx_t mincols_;

public:
    typedef void* ctxt_t;
    typedef void ( *callback_t )( ctxt_t );

private:
    friend class ConnectivityPrivateAccess;
    ctxt_t ctxt_;
    callback_t callback_update_;
    callback_t callback_delete_;
};

// ----------------------------------------------------------------------------------------------

class MultiBlockConnectivityImpl : public IrregularConnectivityImpl {
public:
    //-- Constructors

    MultiBlockConnectivityImpl( const std::string& name = "" );

    explicit MultiBlockConnectivityImpl( eckit::Stream& );

    virtual ~MultiBlockConnectivityImpl();

    //-- Accessors

    ATLAS_HOST_DEVICE
    idx_t blocks() const { return blocks_; }

    ATLAS_HOST_DEVICE
    const BlockConnectivityImpl& block( idx_t block_idx ) const { return block_[block_idx]; }
    ATLAS_HOST_DEVICE
    BlockConnectivityImpl& block( idx_t block_idx ) { return block_[block_idx]; }

    ATLAS_HOST_DEVICE
    idx_t operator()( idx_t row_idx, idx_t col_idx ) const;

    ATLAS_HOST_DEVICE
    idx_t operator()( idx_t block_idx, idx_t block_row_idx, idx_t block_col_idx ) const;

    //-- Modifiers

    virtual void add( idx_t rows, idx_t cols, const idx_t values[], bool fortran_array = false );

    virtual void add( const BlockConnectivityImpl& );

    virtual void add( idx_t rows, idx_t cols );

    virtual void add( idx_t rows, const idx_t cols[] );

    virtual void insert( idx_t position, idx_t rows, idx_t cols, const idx_t values[], bool fortran_array = false );

    virtual void insert( idx_t position, idx_t rows, idx_t cols );

    virtual void insert( idx_t position, idx_t rows, const idx_t cols[] );

    virtual void clear();

    virtual size_t footprint() const;

    virtual void encode( eckit::Stream& s ) const {
        IrregularConnectivityImpl::encode_( s );
        encode_( s );
    }
    virtual void decode( eckit::Stream& s ) {
        IrregularConnectivityImpl::decode_( s );
        decode_( s );
    }

protected:
    void encode_( eckit::Stream& ) const;

    void decode_( eckit::Stream& );

    friend eckit::Stream& operator<<( eckit::Stream& s, const MultiBlockConnectivityImpl& x ) {
        x.encode( s );
        return s;
    }

    friend eckit::Stream& operator>>( eckit::Stream& s, MultiBlockConnectivityImpl& x ) {
        x.decode( s );
        return s;
    }

private:
    void rebuild_block_connectivity();

private:
    idx_t blocks_;
    array::SVector<idx_t> block_displs_;
    array::SVector<idx_t> block_cols_;

    array::SVector<BlockConnectivityImpl> block_;
};

// -----------------------------------------------------------------------------------------------------

class BlockConnectivityImpl {
private:
    friend class IrregularConnectivityImpl;
    friend class MultiBlockConnectivityImpl;
    BlockConnectivityImpl( idx_t rows, idx_t cols, idx_t values[], bool dummy );

public:
    //-- Constructors

    BlockConnectivityImpl();
    BlockConnectivityImpl( idx_t rows, idx_t cols, const std::initializer_list<idx_t>& );

    BlockConnectivityImpl( idx_t rows, idx_t cols, idx_t values[] );

#ifdef __CUDACC__
    // This ctr has to be defined in the header, since __CUDACC__ will identify
    // whether it is compiled it for a GPU kernel
    BlockConnectivityImpl( const BlockConnectivityImpl& other ) :
        owns_( false ),
        values_( other.values_ ),
        rows_( other.rows_ ),
        cols_( other.cols_ ),
        missing_value_( other.missing_value_ ) {}
#endif

    BlockConnectivityImpl( BlockConnectivityImpl&& other );

    BlockConnectivityImpl& operator=( const BlockConnectivityImpl& other );
    BlockConnectivityImpl& operator=( BlockConnectivityImpl&& other );

    explicit BlockConnectivityImpl( eckit::Stream& );

    ~BlockConnectivityImpl();

    ATLAS_HOST_DEVICE
    idx_t index( idx_t i, idx_t j ) const;

    //-- Accessors

    ATLAS_HOST_DEVICE
    idx_t operator()( idx_t row_idx, idx_t col_idx ) const;

    ATLAS_HOST_DEVICE
    idx_t rows() const { return rows_; }

    ATLAS_HOST_DEVICE
    idx_t cols() const { return cols_; }

    ATLAS_HOST_DEVICE
    const idx_t* data() const { return values_.data(); }
    ATLAS_HOST_DEVICE
    idx_t* data() { return values_.data(); }

    ATLAS_HOST_DEVICE
    idx_t missing_value() const { return missing_value_; }

    size_t footprint() const;

    //-- Modifiers

    ATLAS_HOST_DEVICE
    void set( idx_t row_idx, const idx_t column_values[] );

    ATLAS_HOST_DEVICE
    void set( idx_t row_idx, idx_t col_idx, const idx_t value );

    void add( idx_t rows, idx_t cols, const idx_t values[], bool fortran_array = false );

    bool owns() const { return owns_; }


    friend std::ostream& operator<<( std::ostream& out, const BlockConnectivityImpl& x ) {
        x.print( out );
        return out;
    }


protected:
    void rebuild( idx_t rows, idx_t cols, idx_t values[] );


    void encode( eckit::Stream& ) const;

    void decode( eckit::Stream& );

    friend eckit::Stream& operator<<( eckit::Stream& s, const BlockConnectivityImpl& x ) {
        x.encode( s );
        return s;
    }

    friend eckit::Stream& operator>>( eckit::Stream& s, BlockConnectivityImpl& x ) {
        x.decode( s );
        return s;
    }

    void print( std::ostream& ) const;

private:
    bool owns_;
    array::SVector<idx_t> values_;

    idx_t rows_;
    idx_t cols_;
    idx_t missing_value_;
};

using IrregularConnectivity  = ConnectivityInterface<IrregularConnectivityImpl>;
using MultiBlockConnectivity = ConnectivityInterface<MultiBlockConnectivityImpl>;
using BlockConnectivity      = BlockConnectivityImpl;
using Connectivity           = IrregularConnectivity;

// -----------------------------------------------------------------------------------------------------

inline idx_t IrregularConnectivityImpl::operator()( idx_t row_idx, idx_t col_idx ) const {
    assert( counts_[row_idx] > ( col_idx ) );
    return values_[displs_[row_idx] + col_idx] FROM_FORTRAN;
}

inline void IrregularConnectivityImpl::set( idx_t row_idx, const idx_t column_values[] ) {
    const idx_t N = counts_[row_idx];
    for ( idx_t n = 0; n < N; ++n ) {
        values_[displs_[row_idx] + n] = column_values[n] TO_FORTRAN;
    }
}

inline void IrregularConnectivityImpl::set( idx_t row_idx, idx_t col_idx, const idx_t value ) {
    assert( col_idx < counts_[row_idx] );
    values_[displs_[row_idx] + col_idx] = value TO_FORTRAN;
}

inline IrregularConnectivityImpl::Row IrregularConnectivityImpl::row( idx_t row_idx ) const {
    return IrregularConnectivityImpl::Row( const_cast<idx_t*>( values_.data() ) + displs_( row_idx ),
                                           counts_( row_idx ) );
}

// -----------------------------------------------------------------------------------------------------

inline idx_t MultiBlockConnectivityImpl::operator()( idx_t row_idx, idx_t col_idx ) const {
    return IrregularConnectivityImpl::operator()( row_idx, col_idx );
}

inline idx_t MultiBlockConnectivityImpl::operator()( idx_t block_idx, idx_t block_row_idx, idx_t block_col_idx ) const {
    return block( block_idx )( block_row_idx, block_col_idx );
}

// -----------------------------------------------------------------------------------------------------

inline idx_t BlockConnectivityImpl::operator()( idx_t row_idx, idx_t col_idx ) const {
    return values_[index( row_idx, col_idx )] FROM_FORTRAN;
}

inline void BlockConnectivityImpl::set( idx_t row_idx, const idx_t column_values[] ) {
    for ( idx_t n = 0; n < cols_; ++n ) {
        values_[index( row_idx, n )] = column_values[n] TO_FORTRAN;
    }
}

inline void BlockConnectivityImpl::set( idx_t row_idx, idx_t col_idx, const idx_t value ) {
    values_[index( row_idx, col_idx )] = value TO_FORTRAN;
}

inline idx_t BlockConnectivityImpl::index( idx_t i, idx_t j ) const {
    return i * cols_ + j;
}

// ------------------------------------------------------------------------------------------------------

extern "C" {
Connectivity* atlas__Connectivity__create();
MultiBlockConnectivity* atlas__MultiBlockConnectivity__create();
const char* atlas__Connectivity__name( Connectivity* This );
void atlas__Connectivity__rename( Connectivity* This, const char* name );
void atlas__Connectivity__delete( Connectivity* This );
void atlas__Connectivity__displs( Connectivity* This, idx_t*& displs, idx_t& size );
void atlas__Connectivity__counts( Connectivity* This, idx_t*& counts, idx_t& size );
void atlas__Connectivity__values( Connectivity* This, idx_t*& values, idx_t& size );
idx_t atlas__Connectivity__rows( const Connectivity* This );
void atlas__Connectivity__add_values( Connectivity* This, idx_t rows, idx_t cols, idx_t values[] );
void atlas__Connectivity__add_missing( Connectivity* This, idx_t rows, idx_t cols );
idx_t atlas__Connectivity__missing_value( const Connectivity* This );

idx_t atlas__MultiBlockConnectivity__blocks( const MultiBlockConnectivity* This );
BlockConnectivityImpl* atlas__MultiBlockConnectivity__block( MultiBlockConnectivity* This, idx_t block_idx );

idx_t atlas__BlockConnectivity__rows( const BlockConnectivityImpl* This );
idx_t atlas__BlockConnectivity__cols( const BlockConnectivityImpl* This );
idx_t atlas__BlockConnectivity__missing_value( const BlockConnectivityImpl* This );
void atlas__BlockConnectivity__data( BlockConnectivityImpl* This, idx_t*& data, idx_t& rows, idx_t& cols );
void atlas__BlockConnectivity__delete( BlockConnectivityImpl* This );
}

#undef FROM_FORTRAN
#undef TO_FORTRAN
#undef INDEX_DEREF

//------------------------------------------------------------------------------------------------------

}  // namespace mesh

namespace array {

// TODO HACK
// template<> inline DataType::kind_t
// DataType::kind<mesh::BlockConnectivityImpl*>()   { return KIND_INT32;    }
}

}  // namespace atlas