Homec4science
Diffusion tamaas (features/materials) Browse

tamaas/src/grid.hh
b445d515369bfeatures/materials

Last Change
Raw File

grid.hh
View Options

/**
*
* @author Lucas Frérot <lucas.frerot@epfl.ch>
*
* @section LICENSE
*
* Copyright (©) 2016 EPFL (Ecole Polytechnique Fédérale de
* Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des
* Solides)
*
* Tamaas is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Tamaas 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 Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Tamaas. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#ifndef __GRID_HH__
#define __GRID_HH__
/* -------------------------------------------------------------------------- */
#include "tamaas.hh"
#include "array.hh"
#include <cstddef>
/* -------------------------------------------------------------------------- */
__BEGIN_TAMAAS__
/**
* @brief Multi-dimensional & multi-component array class
*
* This class is a container for multi-component data stored on a multi-
* dimensional grid.
*
* The access function is the parenthesis operator. For a grid of dimension d,
* the operator takes d+1 arguments: the first d arguments are the position on
* the grid and the last one is the component of the stored data.
*
* It is also possible to use the access operator with only one argument, it is
* then considering the grid as a flat array, accessing the given cell of the
* array.
*/
template <typename T, UInt dim>
class Grid {
/* -------------------------------------------------------------------------- */
/* Constructors */
/* -------------------------------------------------------------------------- */
public:
/// Constructor by default (empty array)
Grid();
/// Constructor
Grid(const UInt n[dim], const Real L[dim], UInt nb_components);
/// Destructor
virtual ~Grid();
/* -------------------------------------------------------------------------- */
/* Iterator class */
/* -------------------------------------------------------------------------- */
struct iterator {
/// constructor
iterator(T * start, ptrdiff_t offset):ptr(start), offset(offset) {}
/// destructor
~iterator() {}
/// pre-increment
iterator & operator++() { ptr += offset; return *this; }
/// comparison
bool operator<(const iterator & a) { return ptr < a.ptr; }
/// increment with given offset
iterator & operator+=(ptrdiff_t a) { ptr += a*offset; return *this;}
/// dereference iterator
T * operator*() { return ptr; }
/// needed for OpenMP range calculations
Int operator-(const iterator & a) const { return (ptr - a.ptr)/offset; }
private:
T * ptr;
ptrdiff_t offset;
};
public:
/// Begin iterator with n components
iterator begin(UInt n = 1) { return iterator(data.getPtr(), n); }
/// End iterator with n components
iterator end(UInt n = 1) { return iterator(data.getPtr()+computeSize(), n); }
public:
/* -------------------------------------------------------------------------- */
/* Common operations */
/* -------------------------------------------------------------------------- */
/// Resize array
void resize(const UInt n[dim]);
/// Compute size
inline UInt computeSize() const;
/// Print
void printself(std::ostream & str) const;
/// Get sizes
const UInt * sizes() const {return n;}
/// Get total size
UInt dataSize() const {return data.size(); }
/// Get number of points
UInt getNbPoints() const {return dataSize() / getNbComponents(); }
/// Get lengths
const Real * lengths() const {return L;}
/// Get number of components
UInt getNbComponents() const {return nb_components;}
/// Set number of components
void setNbComponents(UInt n) {nb_components = n;}
/// Get internal data pointer (const)
const T * getInternalData() const {return data.getPtr();}
/// Get internal data pointer (non-const)
T * getInternalData() {return data.getPtr();}
/// Set components
void uniformSetComponents(const Grid<T, 1> & vec);
/* -------------------------------------------------------------------------- */
/* Access operators */
/* -------------------------------------------------------------------------- */
/// Variadic access operator (non-const)
template <typename... T1>
inline T & operator()(T1... args);
/// Variadic access operator
template <typename... T1>
inline const T & operator()(T1... args) const;
/// Tuple index access operator
inline T & operator()(UInt tuple[dim+1]);
inline const T & operator()(UInt tuple[dim+1]) const;
private:
/// Unpacking the arguments of variadic () operator when nargs == dim+1
template <typename... T1>
inline UInt unpackOffset(UInt offset, UInt index, T1... rest) const;
/// End case for recursion
template <typename... T1>
inline UInt unpackOffset(UInt offset, UInt index) const;
template <typename... T1>
/// Unpacking the arguments of variadic () operator when nargs == dim (for nb_components == 1)
inline UInt unpackOffsetReduced(UInt offset, UInt index, T1... rest) const;
template <typename... T1>
/// End case for recursion
inline UInt unpackOffsetReduced(UInt offset, UInt index) const;
/// Computing offset for a tuple index
inline UInt computeOffset(UInt tuple[dim+1]) const;
/* -------------------------------------------------------------------------- */
/* Arithmetic operators */
/* -------------------------------------------------------------------------- */
public:
#define GRID_VEC_OPERATOR(op) \
template <typename T1> \
void operator op(const Grid<T1, dim> & other) \
GRID_VEC_OPERATOR(+=);
GRID_VEC_OPERATOR(*=);
GRID_VEC_OPERATOR(-=);
GRID_VEC_OPERATOR(/=);
#undef GRID_VEC_OPERATOR
#define GRID_SCALAR_OPERATOR(op) \
void operator op(const T & e) \
GRID_SCALAR_OPERATOR(+=);
GRID_SCALAR_OPERATOR(*=);
GRID_SCALAR_OPERATOR(-=);
GRID_SCALAR_OPERATOR(/=);
GRID_SCALAR_OPERATOR(=);
#undef GRID_SCALAR_OPERATOR
// = operator
void operator=(const Grid<T, dim> & other);
// = operator (not const input, otherwise gcc is confused)
void operator=(Grid<T, dim> & other);
// Copy data from another grid
template <typename T1>
void copy(const Grid<T1, dim> & other);
/* -------------------------------------------------------------------------- */
/* Member variables */
/* -------------------------------------------------------------------------- */
protected:
UInt n[dim];
Real L[dim];
Array<T> data;
UInt nb_components;
};
/* -------------------------------------------------------------------------- */
/* Inline function definitions */
/* -------------------------------------------------------------------------- */
template <typename T, UInt dim>
inline UInt Grid<T, dim>::computeSize() const {
UInt size = 1;
for (UInt i = 0 ; i < dim ; i++) size *= n[i];
size *= nb_components;
return size;
}
/* -------------------------------------------------------------------------- */
template <typename T, UInt dim>
template <typename... T1>
inline UInt Grid<T, dim>::unpackOffset(UInt offset, UInt index, T1... rest) const {
UInt size = sizeof...(T1);
offset += index;
offset *= (size == 1) ? this->nb_components : this->n[dim-size+1];
return unpackOffset(offset, rest...); // tail-rec bb
}
template <typename T, UInt dim>
template <typename... T1>
inline UInt Grid<T, dim>::unpackOffset(UInt offset, UInt index) const {
return offset + index;
}
template <typename T, UInt dim>
template <typename... T1>
inline UInt Grid<T, dim>::unpackOffsetReduced(UInt offset, UInt index, T1... rest) const {
UInt size = sizeof...(T1);
offset += index;
offset *= this->n[dim-size-1];
return unpackOffsetReduced(offset, rest...);
}
template <typename T, UInt dim>
template <typename... T1>
inline UInt Grid<T, dim>::unpackOffsetReduced(UInt offset, UInt index) const {
return unpackOffset(offset, index);
}
template <typename T, UInt dim>
template <typename... T1>
inline T & Grid<T, dim>::operator()(T1... args) {
TAMAAS_ASSERT(sizeof...(T1) == dim + 1 || sizeof...(T1) == 1 ||
(sizeof...(T1) == dim && this->nb_components == 1),
"operator() does not match dimension");
UInt offset = (sizeof...(T1) == dim) ?
unpackOffsetReduced(0, args...) : unpackOffset(0, args...);
return this->data[offset];
}
template <typename T, UInt dim>
template <typename... T1>
inline const T & Grid<T, dim>::operator()(T1... args) const {
TAMAAS_ASSERT(sizeof...(T1) == dim + 1 || sizeof...(T1) == 1 ||
(sizeof...(T1) == dim && this->nb_components == 1),
"operator() does not match dimension");
UInt offset = (sizeof...(T1) == dim) ?
unpackOffsetReduced(0, args...) : unpackOffset(0, args...);
return this->data[offset];
}
template <typename T, UInt dim>
inline UInt Grid<T, dim>::computeOffset(UInt tuple[dim+1]) const {
UInt offset = 0;
for (UInt i = 0 ; i < dim-1 ; i++) {
offset += tuple[i];
offset *= n[i+1];
}
offset += tuple[dim-1];
offset *= nb_components;
return offset + tuple[dim];
}
template <typename T, UInt dim>
inline T & Grid<T, dim>::operator()(UInt tuple[dim+1]) {
UInt offset = computeOffset(tuple);
return this->data[offset];
}
template <typename T, UInt dim>
inline const T & Grid<T, dim>::operator()(UInt tuple[dim+1]) const {
UInt offset = computeOffset(tuple);
return this->data[offset];
}
/* -------------------------------------------------------------------------- */
#define GRID_VEC_OPERATOR_IMPL(op) \
template <typename T, UInt dim> \
template <typename T1> \
inline void Grid<T, dim>::operator op(const Grid<T1, dim> & other) { \
TAMAAS_ASSERT(other.computeSize() == data.size(), "surface size does not match");\
_Pragma("omp parallel for") \
for (UInt i = 0 ; i < this->data.size() ; i++) { \
data[i] op other(i); \
} \
} \
GRID_VEC_OPERATOR_IMPL(+=);
GRID_VEC_OPERATOR_IMPL(-=);
GRID_VEC_OPERATOR_IMPL(*=);
GRID_VEC_OPERATOR_IMPL(/=);
#undef GRID_VEC_OPERATOR
/* -------------------------------------------------------------------------- */
template <typename T, UInt dim>
void Grid<T, dim>::operator=(const Grid<T, dim> & other) {
this->copy(other);
}
template <typename T, UInt dim>
void Grid<T, dim>::operator=(Grid<T, dim> & other) {
this->copy(other);
}
/* -------------------------------------------------------------------------- */
template <typename T, UInt dim>
template <typename T1>
void Grid<T, dim>::copy(const Grid<T1, dim> & other) {
resize(other.sizes());
#pragma omp parallel for
for (UInt i = 0 ; i < data.size() ; i++)
data[i] = other(i);
}
/* -------------------------------------------------------------------------- */
/* Stream operator */
/* -------------------------------------------------------------------------- */
template <typename T, UInt dim>
inline std::ostream & operator<<(std::ostream & stream,
const Grid<T, dim> & _this) {
_this.printself(stream);
return stream;
}
/* -------------------------------------------------------------------------- */
__END_TAMAAS__
#endif // __GRID_HH__
c4science · Help