diff --git a/python/cast.hh b/python/cast.hh
index b861add..803e119 100644
--- a/python/cast.hh
+++ b/python/cast.hh
@@ -1,182 +1,182 @@
/**
* @file
* @section LICENSE
*
* Copyright (©) 2016-2021 EPFL (École Polytechnique Fédérale de Lausanne),
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU Affero General Public License as published
* by the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* 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 Affero General Public License for more details.
*
* You should have received a copy of the GNU Affero General Public License
* along with this program. If not, see .
*
*/
/* -------------------------------------------------------------------------- */
#ifndef CAST_HH
#define CAST_HH
/* -------------------------------------------------------------------------- */
#include "grid_base.hh"
#include "grid.hh"
#include "numpy.hh"
/* -------------------------------------------------------------------------- */
#include
#include
/* -------------------------------------------------------------------------- */
namespace pybind11 {
// Format descriptor necessary for correct wrap of tamaas complex type
template
struct format_descriptor<
tamaas::complex, detail::enable_if_t::value>> {
static constexpr const char c = format_descriptor::c;
static constexpr const char value[3] = {'Z', c, '\0'};
static std::string format() { return std::string(value); }
};
#ifndef PYBIND11_CPP17
template
constexpr const char format_descriptor<
tamaas::complex,
detail::enable_if_t::value>>::value[3];
#endif
namespace detail {
// declare tamaas complex as a complex type for pybind11
template
struct is_complex> : std::true_type {};
template
struct is_fmt_numeric,
detail::enable_if_t::value>>
: std::true_type {
static constexpr int index = is_fmt_numeric::index + 3;
};
static inline handle policy_switch(return_value_policy policy, handle parent) {
switch (policy) {
case return_value_policy::copy:
case return_value_policy::move:
return handle();
case return_value_policy::automatic_reference: // happens in python-derived
// classes
case return_value_policy::reference:
return none();
case return_value_policy::reference_internal:
return parent;
default:
TAMAAS_EXCEPTION("Policy is not handled");
}
}
template
handle grid_to_python(const tamaas::Grid& grid,
return_value_policy policy, handle parent) {
parent = policy_switch(policy, parent); // reusing variable
std::vector sizes(dim);
std::copy(grid.sizes().begin(), grid.sizes().end(), sizes.begin());
if (grid.getNbComponents() != 1)
sizes.push_back(grid.getNbComponents());
return array(sizes, grid.getInternalData(), parent).release();
}
/**
* Type caster for grid classes
* inspired by https://tinyurl.com/y8m47qh3 from T. De Geus
* and pybind11/eigen.h
*/
template class G, typename T,
tamaas::UInt dim>
struct type_caster> {
using type = G;
using array_type =
array_t;
public:
PYBIND11_TYPE_CASTER(type, _("GridWrap"));
/**
* Conversion part 1 (Python->C++): convert a PyObject into a grid
* instance or return false upon failure. The second argument
* indicates whether implicit conversions should be applied.
*/
bool load(handle src, bool convert) {
if (!array_type::check_(src) or !convert)
return false;
auto buf = array_type::ensure(src);
if (!buf) return false;
value.move(tamaas::wrap::GridNumpy>(buf));
return true;
}
/**
* Conversion part 2 (C++ -> Python): convert a grid instance into
* a Python object. The second and third arguments are used to
* indicate the return value policy and parent object (for
* ``return_value_policy::reference_internal``) and are generally
* ignored by implicit casters.
*/
static handle cast(const type& src, return_value_policy policy,
handle parent) {
return grid_to_python(src, policy, parent);
}
};
/**
* Type caster for GridBase classes
*/
template
struct type_caster> {
using type = tamaas::GridBase;
using array_type =
array_t;
public:
PYBIND11_TYPE_CASTER(type, _("GridBaseWrap"));
bool load(handle src, bool convert) {
if (!array_type::check_(src) or !convert)
return false;
auto buf = array_type::ensure(src);
if (!buf) return false;
value.move(tamaas::wrap::GridBaseNumpy(buf));
return true;
}
static handle cast(const type& src, return_value_policy policy,
handle parent) {
#define GRID_BASE_CASE(unused1, unused2, dim) \
case dim: { \
const auto& conv = dynamic_cast&>(src); \
return grid_to_python(conv, policy, parent); \
}
switch (src.getDimension()) {
BOOST_PP_SEQ_FOR_EACH(GRID_BASE_CASE, ~, (1)(2)(3));
default:
return nullptr;
}
#undef GRID_BASE_CASE
}
};
} // namespace detail
} // namespace pybind11
-#endif // __CAST_HH__
+#endif // CAST_HH
diff --git a/python/numpy.hh b/python/numpy.hh
index 80a840f..4ae8e31 100644
--- a/python/numpy.hh
+++ b/python/numpy.hh
@@ -1,93 +1,93 @@
/**
* @file
* @section LICENSE
*
* Copyright (©) 2016-2021 EPFL (École Polytechnique Fédérale de Lausanne),
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU Affero General Public License as published
* by the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* 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 Affero General Public License for more details.
*
* You should have received a copy of the GNU Affero General Public License
* along with this program. If not, see .
*
*/
/* -------------------------------------------------------------------------- */
-#ifndef __NUMPY_HH__
-#define __NUMPY_HH__
+#ifndef NUMPY_HH
+#define NUMPY_HH
/* -------------------------------------------------------------------------- */
#include "grid_base.hh"
#include "tamaas.hh"
/* -------------------------------------------------------------------------- */
#include
#include
namespace tamaas {
namespace wrap {
namespace py = pybind11;
/// Type alias for usual numpy array type
template
using numpy = py::array_t;
/// GridBase helper class wrapping numpy array
template
class GridBaseNumpy : public GridBase {
public:
GridBaseNumpy(numpy& buffer) : GridBase() {
this->nb_components = 1;
this->data.wrapMemory(buffer.mutable_data(), buffer.size());
}
};
/// Grid helper class wrapping numpy array
template
class GridNumpy : public Parent {
public:
GridNumpy(numpy& buffer) : Parent() {
auto* array_shape = buffer.shape();
const UInt ndim = buffer.ndim();
if (ndim > Parent::dimension + 1 || ndim < Parent::dimension)
TAMAAS_EXCEPTION(
"Numpy array dimension do not match expected grid dimensions");
if (ndim == Parent::dimension + 1)
this->nb_components = array_shape[ndim - 1];
std::copy_n(array_shape, Parent::dimension, this->n.begin());
this->computeStrides();
this->data.wrapMemory(buffer.mutable_data(), this->computeSize());
}
};
/// Surface helper class wrapping numpy array
template
class SurfaceNumpy : public Parent {
public:
SurfaceNumpy(numpy& buffer) : Parent() {
auto* array_shape = buffer.shape();
if (buffer.ndim() != 2)
TAMAAS_EXCEPTION(
"Numpy array dimension do not match expected surface dimensions");
std::copy_n(array_shape, Parent::dimension, this->n.begin());
this->computeStrides();
this->data.wrapMemory(buffer.mutable_data(), this->computeSize());
}
};
} // namespace wrap
} // namespace tamaas
-#endif // __NUMPY_HH__
+#endif // NUMPY_HH
diff --git a/python/wrap.hh b/python/wrap.hh
index 6c54053..9d182a2 100644
--- a/python/wrap.hh
+++ b/python/wrap.hh
@@ -1,77 +1,77 @@
/**
* @file
* @section LICENSE
*
* Copyright (©) 2016-2021 EPFL (École Polytechnique Fédérale de Lausanne),
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU Affero General Public License as published
* by the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* 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 Affero General Public License for more details.
*
* You should have received a copy of the GNU Affero General Public License
* along with this program. If not, see .
*
*/
/* -------------------------------------------------------------------------- */
-#ifndef __WRAP_HH__
-#define __WRAP_HH__
+#ifndef WRAP_HH
+#define WRAP_HH
/* -------------------------------------------------------------------------- */
#include "cast.hh"
#include "numpy.hh"
#include "tamaas.hh"
/* -------------------------------------------------------------------------- */
#include
/* -------------------------------------------------------------------------- */
namespace tamaas {
namespace wrap {
namespace py = pybind11;
/// For naming classes templated with dimension
inline std::string makeDimensionName(const std::string& name, UInt dim) {
std::stringstream str;
str << name << dim << "D";
return str.str();
}
/* -------------------------------------------------------------------------- */
/* Prototypes for wrapping of tamaas components */
/* -------------------------------------------------------------------------- */
void wrapCore(py::module& mod);
void wrapPercolation(py::module& mod);
void wrapSurface(py::module& mod);
void wrapModel(py::module& mod);
void wrapSolvers(py::module& mod);
void wrapTestFeatures(py::module& mod);
void wrapCompute(py::module& mod);
void wrapMPI(py::module& mod);
/* -------------------------------------------------------------------------- */
/* Deprecation macro */
/* -------------------------------------------------------------------------- */
#define TAMAAS_DEPRECATE(olds, news) \
do { \
PyErr_WarnEx(PyExc_DeprecationWarning, \
olds " is deprecated, use " news " instead.", 1); \
} while (0)
#define TAMAAS_DEPRECATE_ACCESSOR(acc, type, property) \
#acc, [](const type& m) -> decltype(m.acc()) { \
TAMAAS_DEPRECATE(#acc "()", "the " property " property"); \
return m.acc(); \
}
} // namespace wrap
} // namespace tamaas
-#endif // __WRAP_HH__
+#endif // WRAP_HH
diff --git a/src/core/array.hh b/src/core/array.hh
index d5924ba..70f4122 100644
--- a/src/core/array.hh
+++ b/src/core/array.hh
@@ -1,176 +1,176 @@
/**
* @file
* @section LICENSE
*
* Copyright (©) 2016-2021 EPFL (École Polytechnique Fédérale de Lausanne),
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU Affero General Public License as published
* by the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* 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 Affero General Public License for more details.
*
* You should have received a copy of the GNU Affero General Public License
* along with this program. If not, see .
*
*/
/* -------------------------------------------------------------------------- */
-#ifndef __ARRAY__HH__
-#define __ARRAY__HH__
+#ifndef ARRAY_HH
+#define ARRAY_HH
/* -------------------------------------------------------------------------- */
#include "logger.hh"
#include "tamaas.hh"
#include
#include
#include
#include
/* -------------------------------------------------------------------------- */
namespace tamaas {
/**
* @brief Generic storage class with wrapping capacities
*/
template
class Array final {
public:
/// Default
Array() = default;
/// Empty array of given size
Array(UInt size) : Array() { this->resize(size); }
/// Copy constructor (deep)
Array(const Array& v) : Array() {
this->resize(v._size);
thrust::copy_n(v._data, v._size, this->_data);
}
/// Move constructor (transfers data ownership)
Array(Array&& v) noexcept : Array() {
_data = std::exchange(v._data, nullptr);
_size = std::exchange(v._size, 0);
_reserved = std::exchange(v._reserved, 0);
wrapped = std::exchange(v.wrapped, false);
}
/// Wrap array on data
Array(T* data, UInt size) : Array() { this->wrapMemory(data, size); }
/// Destructor
~Array() {
if (not wrapped)
_alloc.deallocate(_data, _reserved);
}
public:
/// Copy operator
Array& operator=(const Array& v) {
this->wrapped = false;
if (v.size() != this->size())
this->resize(v.size());
thrust::copy_n(v._data, _size, _data);
return *this;
}
/// Move operator
Array& operator=(Array&& v) noexcept {
if (this != &v) {
if (not wrapped)
_alloc.deallocate(_data, _reserved);
_data = std::exchange(v._data, nullptr);
_size = std::exchange(v._size, 0);
_reserved = std::exchange(v._reserved, 0);
wrapped = std::exchange(v.wrapped, false);
}
return *this;
}
void wrap(Array& other) { this->wrapMemory(other._data, other._size); }
void wrap(const Array& other) { this->wrapMemory(other._data, other._size); }
/// Wrap a memory pointer
void wrapMemory(T* data, UInt size) {
_data = data;
_size = size;
wrapped = true;
_reserved = 0;
}
/// Assign a sigle value to the array
void assign(UInt size, const T& value) {
this->resize(size);
thrust::fill(_data, _data + size, value);
}
void setWrapped(bool w) { wrapped = w; }
/// Data pointer access (const)
const T* data() const { return _data; }
/// Data pointer access (non-const)
T* data() { return _data; }
/// Data pointer setter
void setData(T* new_ptr) { _data = new_ptr; }
/// Resize array
void resize(UInt size) {
if (wrapped)
TAMAAS_EXCEPTION("cannot resize wrapped array");
// Erase array
if (size == 0) {
_alloc.deallocate(_data, _reserved);
_data = nullptr;
_size = 0;
_reserved = 0;
return;
}
// Do nothing
if (size == this->_size)
return;
// Allocate new data
_alloc.deallocate(_data, _reserved);
_data = _alloc.allocate(size);
_size = size;
_reserved = size;
}
/// Reserve storage space
void reserve(UInt size) {
if (_reserved >= size)
return;
auto new_data = _alloc.allocate(size);
if (new_data != _data) {
thrust::copy_n(_data, _size, new_data);
_alloc.deallocate(_data, _reserved);
_data = new_data;
_reserved = size;
}
}
/// Access operator
inline T& operator[](UInt i) {
TAMAAS_ASSERT(i < _size, "memory overflow");
return _data[i];
}
/// Access operator (const)
inline const T& operator[](UInt i) const {
TAMAAS_ASSERT(i < _size, "memory overflow");
return _data[i];
}
/// Get size of array
inline UInt size() const { return _size; }
private:
T* _data = nullptr;
UInt _size = 0;
std::size_t _reserved = 0;
bool wrapped = false;
Allocator _alloc;
};
} // namespace tamaas
/* -------------------------------------------------------------------------- */
-#endif /* __ARRAY__HH__ */
+#endif /* ARRAY_HH */
diff --git a/src/core/fftw_allocator.hh b/src/core/fftw_allocator.hh
index 150bc61..fa9dcef 100644
--- a/src/core/fftw_allocator.hh
+++ b/src/core/fftw_allocator.hh
@@ -1,52 +1,52 @@
/**
* @file
* @section LICENSE
*
* Copyright (©) 2016-2021 EPFL (École Polytechnique Fédérale de Lausanne),
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU Affero General Public License as published
* by the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* 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 Affero General Public License for more details.
*
* You should have received a copy of the GNU Affero General Public License
* along with this program. If not, see .
*
*/
/* -------------------------------------------------------------------------- */
-#ifndef __FFTW_ALLOCATOR_HH__
-#define __FFTW_ALLOCATOR_HH__
+#ifndef FFTW_ALLOCATOR_HH
+#define FFTW_ALLOCATOR_HH
/* -------------------------------------------------------------------------- */
#include
#include
/* -------------------------------------------------------------------------- */
namespace tamaas {
/// Class allocating
/// [SIMD](http://www.fftw.org/fftw3_doc/SIMD-alignment-and-fftw_005fmalloc.html#SIMD-alignment-and-fftw_005fmalloc)
/// aligned memory
template
struct FFTWAllocator {
/// Allocate memory
static T* allocate(std::size_t n) noexcept {
T* p = nullptr;
p = (T*)fftw_malloc(sizeof(T) * n);
return p;
}
/// Free memory
static void deallocate(T* p, __attribute__((unused)) std::size_t n) noexcept {
fftw_free(p);
}
};
} // namespace tamaas
-#endif // __FFTW_ALLOCATOR_HH__
+#endif // FFTW_ALLOCATOR_HH
diff --git a/src/core/grid.hh b/src/core/grid.hh
index faf9fbb..a48adcb 100644
--- a/src/core/grid.hh
+++ b/src/core/grid.hh
@@ -1,222 +1,222 @@
/**
* @file
* @section LICENSE
*
* Copyright (©) 2016-2021 EPFL (École Polytechnique Fédérale de Lausanne),
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU Affero General Public License as published
* by the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* 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 Affero General Public License for more details.
*
* You should have received a copy of the GNU Affero General Public License
* along with this program. If not, see .
*
*/
/* -------------------------------------------------------------------------- */
-#ifndef __GRID_HH__
-#define __GRID_HH__
+#ifndef GRID_HH
+#define GRID_HH
/* -------------------------------------------------------------------------- */
#include "grid_base.hh"
#include "tamaas.hh"
#include
#include
#include
#include
/* -------------------------------------------------------------------------- */
namespace 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
class Grid : public GridBase {
template
using is_valid_index = fold_trait;
public:
/* ------------------------------------------------------------------------ */
/* Types */
/* ------------------------------------------------------------------------ */
using value_type = T;
using reference = value_type&;
static constexpr UInt dimension = dim;
/* ------------------------------------------------------------------------ */
/* Constructors */
/* ------------------------------------------------------------------------ */
public:
/// Constructor by default (empty array)
Grid();
/// Constructor from iterators
template
Grid(RandomAccessIterator begin, RandomAccessIterator end,
UInt nb_components);
/// Constructor
Grid(const std::array& n, UInt nb_components)
: Grid(std::begin(n), std::end(n), nb_components) {}
/// Wrap on given data
Grid(const std::array& n, UInt nb_components, value_type* data);
/// Constructor with vectors
Grid(const std::vector& n, UInt nb_components)
: Grid(std::begin(n), std::end(n), nb_components) {}
/// Constructor with initializer list
Grid(const std::initializer_list& n, UInt nb_components)
: Grid(std::begin(n), std::end(n), nb_components) {}
/// Copy constructor
Grid(const Grid& o);
/// Move constructor (transfers data ownership)
Grid(Grid&& o) noexcept;
/// Destructor
~Grid() override = default;
private:
/// Init from standard container
template
void init(const Container& n, UInt nb_components);
public:
/* ------------------------------------------------------------------------ */
/* Common operations */
/* ------------------------------------------------------------------------ */
/// Resize array
void resize(const std::array& n);
/// Resize array (from std::vector)
void resize(const std::vector& n);
/// Resize array (from initializer list)
void resize(std::initializer_list n);
/// Resize array (from iterators)
template
void resize(ForwardIt begin, ForwardIt end);
/// Compute size
inline UInt computeSize() const;
/// Get grid dimension
inline UInt getDimension() const override { return dim; }
/// Compute strides
virtual void computeStrides();
/// Print
virtual void printself(std::ostream& str) const;
/// Get sizes
const std::array& sizes() const { return n; }
/// Get strides
const std::array& getStrides() const { return this->strides; }
/* ------------------------------------------------------------------------ */
/* Access operators (these cannot be moved outside the class) */
/* ------------------------------------------------------------------------ */
/// Variadic access operator (non-const)
template
inline std::enable_if_t::value, T&>
operator()(T1... args) {
/// Checking access integrity
constexpr UInt nargs = sizeof...(T1);
static_assert(nargs == dim + 1 || nargs == 1 || nargs == dim,
"number of arguments in operator() does not match dimension");
constexpr UInt start = (nargs == 1) ? dim : 0;
UInt offset = unpackOffset(0, start, args...);
return this->data[offset];
}
/// Variadic access operator
template
inline std::enable_if_t::value, const T&>
operator()(T1... args) const {
/// Checking access integrity
constexpr UInt nargs = sizeof...(T1);
static_assert(nargs == dim + 1 || nargs == 1 || nargs == dim,
"number of arguments in operator() does not match dimension");
constexpr UInt start = (nargs == 1) ? dim : 0;
UInt offset = unpackOffset(0, start, args...);
return this->data[offset];
}
/// Tuple index access operator
template
inline T& operator()(std::array tuple);
template
inline const T& operator()(std::array tuple) const;
private:
/// Unpacking the arguments of variadic ()
template
inline UInt unpackOffset(UInt offset, UInt index_pos, UInt index,
T1... rest) const;
/// End case for recursion
template
inline UInt unpackOffset(UInt offset, UInt index_pos, UInt index) const;
/// Computing offset for a tuple index
template
inline UInt computeOffset(std::array tuple) const;
/* ------------------------------------------------------------------------ */
/* Move/Copy operators */
/* ------------------------------------------------------------------------ */
public:
using GridBase::operator=;
// = operator
Grid& operator=(const Grid& other);
// = operator (not const input, otherwise gcc is confused)
Grid& operator=(Grid& other);
// = operator (move)
Grid& operator=(Grid&& other) noexcept;
// Copy data from another grid
template
void copy(const Grid& other);
// Move data from another grid
template
void move(Grid&& other) noexcept;
// Wrap memory (non-const)
void wrap(Grid& other) {
GridBase::wrap(other);
std::copy(other.n.begin(), other.n.end(), n.begin());
this->computeStrides();
}
// Wrap memory
void wrap(const Grid& other) {
GridBase::wrap(other);
std::copy(other.n.begin(), other.n.end(), n.begin());
this->computeStrides();
}
/* ------------------------------------------------------------------------ */
/* Member variables */
/* ------------------------------------------------------------------------ */
protected:
std::array n; ///< shape of grid: size per dimension
std::array strides; ///< strides for access
};
} // namespace tamaas
/* -------------------------------------------------------------------------- */
/* Inline/template function definitions */
/* -------------------------------------------------------------------------- */
#include "grid_tmpl.hh"
/* -------------------------------------------------------------------------- */
-#endif // __GRID_HH__
+#endif // GRID_HH
diff --git a/src/core/grid_base.hh b/src/core/grid_base.hh
index 00c8410..f2bd6d4 100644
--- a/src/core/grid_base.hh
+++ b/src/core/grid_base.hh
@@ -1,345 +1,345 @@
/**
* @file
* @section LICENSE
*
* Copyright (©) 2016-2021 EPFL (École Polytechnique Fédérale de Lausanne),
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU Affero General Public License as published
* by the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* 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 Affero General Public License for more details.
*
* You should have received a copy of the GNU Affero General Public License
* along with this program. If not, see .
*
*/
/* -------------------------------------------------------------------------- */
#ifndef GRID_BASE_HH
#define GRID_BASE_HH
/* -------------------------------------------------------------------------- */
#include "array.hh"
#include "iterator.hh"
#include "loop.hh"
#include "mpi_interface.hh"
#include "static_types.hh"
#include "tamaas.hh"
#include
#include
#include
#include
/* -------------------------------------------------------------------------- */
namespace tamaas {
/**
* @brief Dimension agnostic grid with components stored per points
*/
template
class GridBase {
public:
/// Constructor by default
GridBase() = default;
/// Copy constructor
GridBase(const GridBase& o) { this->copy(o); }
/// Move constructor (transfers data ownership)
GridBase(GridBase&& o) noexcept
: data(std::move(o.data)),
nb_components(std::exchange(o.nb_components, 1)) {}
/// Initialize with size
explicit GridBase(UInt nb_points, UInt nb_components = 1) {
setNbComponents(nb_components);
this->resize(nb_points * nb_components);
}
/// Destructor
virtual ~GridBase() = default;
/* ------------------------------------------------------------------------ */
/* Iterator class */
/* ------------------------------------------------------------------------ */
public:
using iterator = iterator_::iterator;
using const_iterator = iterator_::iterator;
using value_type = T;
using reference = value_type&;
public:
/// Get grid dimension
virtual UInt getDimension() const { return 1; }
/// Get internal data pointer (const)
const T* getInternalData() const { return this->data.data(); }
/// Get internal data pointer (non-const)
T* getInternalData() { return this->data.data(); }
/// Get number of components
UInt getNbComponents() const { return nb_components; }
/// Set number of components
void setNbComponents(UInt n) { nb_components = n; }
/// Get total size
virtual UInt dataSize() const { return this->data.size(); }
/// Get global data size
UInt globalDataSize() const {
return mpi::allreduce(dataSize());
}
/// Get number of points
UInt getNbPoints() const {
return this->dataSize() / this->getNbComponents();
}
/// Get global number of points
UInt getGlobalNbPoints() const {
return this->globalDataSize() / this->getNbComponents();
}
/// Set components
void uniformSetComponents(const GridBase& vec);
/// Resize
void resize(UInt size) { this->data.assign(size, T(0.)); }
/// Reserve storage w/o changing data logic
void reserve(UInt size) { this->data.reserve(size); }
/* ------------------------------------------------------------------------ */
/* Iterator methods */
/* ------------------------------------------------------------------------ */
virtual iterator begin(UInt n = 1) {
return iterator(this->getInternalData(), n);
}
virtual iterator end(UInt n = 1) {
return iterator(this->getInternalData() + this->dataSize(), n);
}
virtual const_iterator begin(UInt n = 1) const {
return const_iterator(this->getInternalData(), n);
}
virtual const_iterator end(UInt n = 1) const {
return const_iterator(this->getInternalData() + this->dataSize(), n);
}
/* ------------------------------------------------------------------------ */
/* Operators */
/* ------------------------------------------------------------------------ */
inline T& operator()(UInt i) { return this->data[i]; }
inline const T& operator()(UInt i) const { return this->data[i]; }
#define VEC_OPERATOR(op) \
template \
GridBase& operator op(const GridBase& other)
VEC_OPERATOR(+=);
VEC_OPERATOR(*=);
VEC_OPERATOR(-=);
VEC_OPERATOR(/=);
#undef VEC_OPERATOR
template
T dot(const GridBase& other) const;
#define SCALAR_OPERATOR(op) GridBase& operator op(const T& e)
SCALAR_OPERATOR(+=);
SCALAR_OPERATOR(*=);
SCALAR_OPERATOR(-=);
SCALAR_OPERATOR(/=);
SCALAR_OPERATOR(=);
#undef SCALAR_OPERATOR
/// Broadcast operators
#define BROADCAST_OPERATOR(op) \
template \
GridBase& operator op(const StaticArray& b)
BROADCAST_OPERATOR(+=);
BROADCAST_OPERATOR(-=);
BROADCAST_OPERATOR(*=);
BROADCAST_OPERATOR(/=);
BROADCAST_OPERATOR(=);
#undef BROADCAST_OPERATOR
/* ------------------------------------------------------------------------ */
/* Min/max */
/* ------------------------------------------------------------------------ */
T min() const;
T max() const;
T sum() const;
T mean() const { return sum() / static_cast(globalDataSize()); }
T var() const;
T l2norm() const { return std::sqrt(dot(*this)); }
/* ------------------------------------------------------------------------ */
/* Move/Copy */
/* ------------------------------------------------------------------------ */
public:
GridBase& operator=(const GridBase& o) {
this->copy(o);
return *this;
}
GridBase& operator=(GridBase& o) {
this->copy(o);
return *this;
}
GridBase& operator=(GridBase&& o) noexcept {
this->move(std::move(o));
return *this;
}
template
void copy(const GridBase& other) {
if (other.dataSize() != this->dataSize())
this->resize(other.dataSize());
thrust::copy(other.begin(), other.end(), this->begin());
nb_components = other.nb_components;
}
template
void move(GridBase&& other) noexcept {
data = std::move(other.data);
nb_components = std::exchange(other.nb_components, 1);
}
GridBase& wrap(GridBase& o) {
data.wrap(o.data);
nb_components = o.nb_components;
return *this;
}
GridBase& wrap(const GridBase& o) {
data.wrap(o.data);
nb_components = o.nb_components;
return *this;
}
protected:
Array data;
UInt nb_components = 1;
};
/* -------------------------------------------------------------------------- */
template
void GridBase::uniformSetComponents(const GridBase& vec) {
if (vec.dataSize() != this->nb_components)
TAMAAS_EXCEPTION("Cannot set grid field with values of vector");
auto begin_it(begin());
auto end_it(end()); // silencing nvcc warnings
for (auto it = begin(); it < end_it; ++it) {
UInt i = it - begin_it;
*it = vec(i % this->nb_components);
}
}
/* --------------------------------------------------------------------------
*/
#define SCALAR_OPERATOR_IMPL(op) \
template \
inline GridBase& GridBase::operator op(const T& e) { \
Loop::loop([e] CUDA_LAMBDA(T& val) { val op e; }, *this); \
return *this; \
}
SCALAR_OPERATOR_IMPL(+=)
SCALAR_OPERATOR_IMPL(*=)
SCALAR_OPERATOR_IMPL(-=)
SCALAR_OPERATOR_IMPL(/=)
SCALAR_OPERATOR_IMPL(=)
#undef SCALAR_OPERATOR_IMPL
/* -------------------------------------------------------------------------- */
#define VEC_OPERATOR_IMPL(op) \
template \
template \
inline GridBase& GridBase::operator op(const GridBase& other) { \
TAMAAS_ASSERT(other.dataSize() == this->dataSize(), \
"surface size does not match"); \
Loop::loop([] CUDA_LAMBDA(T& x, const T1& y) { x op y; }, *this, other); \
return *this; \
}
VEC_OPERATOR_IMPL(+=)
VEC_OPERATOR_IMPL(-=)
VEC_OPERATOR_IMPL(*=)
VEC_OPERATOR_IMPL(/=)
#undef VEC_OPERATOR
/* -------------------------------------------------------------------------- */
#define BROADCAST_OPERATOR_IMPL(op) \
template \
template \
GridBase& GridBase::operator op(const StaticArray& b) { \
TAMAAS_ASSERT(this->dataSize() % size == 0, \
"Broadcast operator cannot broadcast" \
<< size << " on array of size " << this->dataSize()); \
Tensor::type, size> a(b); \
Loop::loop([&a] CUDA_LAMBDA(VectorProxy g) { g op a; }, \
range>(*this)); \
return *this; \
}
BROADCAST_OPERATOR_IMPL(+=)
BROADCAST_OPERATOR_IMPL(-=)
BROADCAST_OPERATOR_IMPL(*=)
BROADCAST_OPERATOR_IMPL(/=)
BROADCAST_OPERATOR_IMPL(=)
#undef BROADCAST_OPERATOR_IMPL
/* -------------------------------------------------------------------------- */
template
inline T GridBase::min() const {
// const auto id = [] CUDA_LAMBDA (const T&x){return x;};
return Loop::reduce([] CUDA_LAMBDA(const T& x) { return x; },
*this);
}
/* -------------------------------------------------------------------------- */
template
inline T GridBase::max() const {
return Loop::reduce([] CUDA_LAMBDA(const T& x) { return x; },
*this);
}
/* -------------------------------------------------------------------------- */
template
inline T GridBase::sum() const {
return Loop::reduce([] CUDA_LAMBDA(const T& x) { return x; },
*this);
}
/* -------------------------------------------------------------------------- */
template
inline T GridBase::var() const {
const T mu = mean();
const T var = Loop::reduce(
[mu] CUDA_LAMBDA(const T& x) { return (x - mu) * (x - mu); }, *this);
return var / (globalDataSize() - 1);
}
/* -------------------------------------------------------------------------- */
template
template
T GridBase::dot(const GridBase& other) const {
return Loop::reduce(
[] CUDA_LAMBDA(const T& x, const T1& y) { return x * y; }, *this, other);
}
} // namespace tamaas
-#endif // __GRID_BASE_HH__
+#endif // GRID_BASE_HH
diff --git a/src/core/grid_hermitian.hh b/src/core/grid_hermitian.hh
index 91ca2b9..93af3df 100644
--- a/src/core/grid_hermitian.hh
+++ b/src/core/grid_hermitian.hh
@@ -1,162 +1,162 @@
/**
* @file
* @section LICENSE
*
* Copyright (©) 2016-2021 EPFL (École Polytechnique Fédérale de Lausanne),
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU Affero General Public License as published
* by the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* 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 Affero General Public License for more details.
*
* You should have received a copy of the GNU Affero General Public License
* along with this program. If not, see .
*
*/
/* -------------------------------------------------------------------------- */
-#ifndef __GRID_HERMITIAN_HH__
-#define __GRID_HERMITIAN_HH__
+#ifndef GRID_HERMITIAN_HH
+#define GRID_HERMITIAN_HH
/* -------------------------------------------------------------------------- */
#include "grid.hh"
#include "tamaas.hh"
#include
#include
#include
/* -------------------------------------------------------------------------- */
namespace tamaas {
static complex dummy;
/**
* @brief Multi-dimensional, multi-component herimitian array
*
* This class represents an array of hermitian data, meaning it has dimensions
* of: n1 * n2 * n3 * ... * (nx / 2 + 1)
*
* However, it acts as a fully dimensioned array, returning a dummy reference
* for data outside its real dimension, allowing one to write full (and
* inefficient) loops without really worrying about the reduced dimension.
*
* It would however be better to just use the true dimensions of the surface
* for all intents and purposes, as it saves computation time.
*/
template
class GridHermitian : public Grid, dim> {
public:
using value_type = complex;
/* ------------------------------------------------------------------------ */
/* Constructors */
/* ------------------------------------------------------------------------ */
public:
GridHermitian() = default;
GridHermitian(const GridHermitian& o) = default;
GridHermitian(GridHermitian&& o) noexcept = default;
using Grid::Grid;
using Grid, dim>::operator=;
/* ------------------------------------------------------------------------ */
/* Access operators */
/* ------------------------------------------------------------------------ */
template
inline complex& operator()(T1... args);
template
inline const complex& operator()(T1... args) const;
static std::array hermitianDimensions(std::array n) {
n[dim - 1] /= 2;
n[dim - 1] += 1;
return n;
}
template ::value>>
static std::vector hermitianDimensions(std::vector n) {
n.back() /= 2;
n.back() += 1;
return n;
}
private:
template
inline void packTuple(UInt* t, UInt i, T1... args) const;
template
inline void packTuple(UInt* t, UInt i) const;
};
/* -------------------------------------------------------------------------- */
/* Inline function definitions */
/* -------------------------------------------------------------------------- */
template
template
inline void GridHermitian::packTuple(UInt* t, UInt i,
T1... args) const {
*t = i;
packTuple(t + 1, args...);
}
template
template
inline void GridHermitian::packTuple(UInt* t, UInt i) const {
*t = i;
}
template
template
inline complex& GridHermitian::operator()(T1... args) {
constexpr UInt nargs = sizeof...(T1);
static_assert(nargs <= dim + 1, "Too many arguments");
UInt last_index = std::get(std::forward_as_tuple(args...));
if (nargs != 1 && last_index >= this->n[dim - 1]) {
std::array tuple = {{0}};
packTuple(tuple.data(), args...);
for (UInt i = 0; i < dim; i++) {
if (tuple[i] && i != dim - 1)
tuple[i] = this->n[i] - tuple[i];
else if (tuple[i] && i == dim - 1)
tuple[i] = 2 * (this->n[i] - 1) - tuple[i];
}
dummy = conj(this->Grid, dim>::operator()(tuple));
return dummy;
}
else
return Grid, dim>::operator()(args...);
}
template
template
inline const complex& GridHermitian::operator()(T1... args) const {
constexpr UInt nargs = sizeof...(T1);
static_assert(nargs <= dim + 1, "Too many arguments");
UInt last_index = std::get(std::forward_as_tuple(args...));
if (nargs != 1 && last_index >= this->n[dim - 1]) {
std::array tuple = {{0}};
packTuple(tuple.data(), args...);
for (UInt i = 0; i < dim; i++) {
if (tuple[i] && i != dim - 1)
tuple[i] = this->n[i] - tuple[i];
else if (tuple[i] && i == dim - 1)
tuple[i] = 2 * (this->n[i] - 1) - tuple[i];
}
dummy = conj(this->Grid, dim>::operator()(tuple));
return dummy;
}
else
return Grid, dim>::operator()(args...);
}
} // namespace tamaas
-#endif // __GRID_HERMITIAN_HH__
+#endif // GRID_HERMITIAN_HH
diff --git a/src/core/grid_tmpl.hh b/src/core/grid_tmpl.hh
index 3333011..3cc57eb 100644
--- a/src/core/grid_tmpl.hh
+++ b/src/core/grid_tmpl.hh
@@ -1,159 +1,159 @@
/**
* @file
* @section LICENSE
*
* Copyright (©) 2016-2021 EPFL (École Polytechnique Fédérale de Lausanne),
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU Affero General Public License as published
* by the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* 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 Affero General Public License for more details.
*
* You should have received a copy of the GNU Affero General Public License
* along with this program. If not, see .
*
*/
/* -------------------------------------------------------------------------- */
-#ifndef __GRID_TMPL_HH__
-#define __GRID_TMPL_HH__
+#ifndef GRID_TMPL_HH
+#define GRID_TMPL_HH
/* -------------------------------------------------------------------------- */
#include "grid.hh"
#include "tamaas.hh"
namespace tamaas {
template
template
Grid::Grid(RandomAccessIterator begin, RandomAccessIterator end,
UInt nb_components)
: GridBase() {
if (end - begin != dim) {
TAMAAS_EXCEPTION("Provided sizes (" << end - begin
<< ") for grid do not match dimension ("
<< dim << ")");
}
this->nb_components = nb_components;
this->resize(begin, end);
}
/* -------------------------------------------------------------------------- */
template
inline UInt Grid::computeSize() const {
UInt size = 1;
for (UInt i = 0; i < dim; i++)
size *= n[i];
size *= this->nb_components;
return size;
}
/* -------------------------------------------------------------------------- */
template
template
void Grid::resize(ForwardIt begin, ForwardIt end) {
std::copy(begin, end, this->n.begin());
UInt size = this->computeSize();
GridBase::resize(size);
this->computeStrides();
}
/* -------------------------------------------------------------------------- */
template
template
inline UInt Grid::unpackOffset(UInt offset, UInt index_pos, UInt index,
T1... rest) const {
offset += index * strides[index_pos];
return unpackOffset(offset, index_pos + 1, rest...); // tail-rec bb
}
template
template
inline UInt Grid::unpackOffset(UInt offset, UInt index_pos,
UInt index) const {
return offset + index * strides[index_pos];
}
/* -------------------------------------------------------------------------- */
template
template
inline UInt Grid::computeOffset(std::array tuple) const {
static_assert(tdim == dim or tdim == dim + 1, "Tuple dimension is invalid");
return std::inner_product(tuple.begin(), tuple.end(), strides.begin(), 0);
}
template
template
inline T& Grid::operator()(std::array tuple) {
UInt offset = computeOffset(tuple);
return this->data[offset];
}
template
template
inline const T& Grid::operator()(std::array tuple) const {
UInt offset = computeOffset(tuple);
return this->data[offset];
}
/* -------------------------------------------------------------------------- */
template
Grid& Grid::operator=(const Grid& other) {
this->copy(other);
return *this;
}
template
Grid& Grid::operator=(Grid& other) {
this->copy(other);
return *this;
}
template
Grid& Grid::operator=(Grid&& other) noexcept {
this->move(std::move(other));
return *this;
}
/* -------------------------------------------------------------------------- */
template
template
void Grid::copy(const Grid& other) {
GridBase::copy(other);
this->n = other.n;
this->strides = other.strides;
}
template
template
void Grid::move(Grid&& other) noexcept {
GridBase::move(std::move(other));
this->n = std::move(other.n);
this->strides = std::move(other.strides);
}
/* -------------------------------------------------------------------------- */
/* Stream operator */
/* -------------------------------------------------------------------------- */
template
inline std::ostream& operator<<(std::ostream& stream,
const Grid& _this) {
_this.printself(stream);
return stream;
}
} // namespace tamaas
/* -------------------------------------------------------------------------- */
-#endif // __GRID_TMPL_HH__
+#endif // GRID_TMPL_HH
diff --git a/src/core/grid_view.hh b/src/core/grid_view.hh
index dec1488..c29afb6 100644
--- a/src/core/grid_view.hh
+++ b/src/core/grid_view.hh
@@ -1,164 +1,164 @@
/**
* @file
* @section LICENSE
*
* Copyright (©) 2016-2021 EPFL (École Polytechnique Fédérale de Lausanne),
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU Affero General Public License as published
* by the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* 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 Affero General Public License for more details.
*
* You should have received a copy of the GNU Affero General Public License
* along with this program. If not, see .
*
*/
/* -------------------------------------------------------------------------- */
-#ifndef __GRID_VIEW_HH__
-#define __GRID_VIEW_HH__
+#ifndef GRID_VIEW_HH
+#define GRID_VIEW_HH
/* -------------------------------------------------------------------------- */
#include "grid.hh"
#include "tamaas.hh"
#include
/* -------------------------------------------------------------------------- */
namespace tamaas {
/* -------------------------------------------------------------------------- */
/**
* @brief View type on grid
* This is a view on a *contiguous* chunk of data defined by a grid
*/
template class Base, typename T, UInt base_dim,
UInt dim>
class GridView : public Base {
public:
using iterator = typename Base::iterator;
using const_iterator = typename Base::const_iterator;
using value_type = typename Base::value_type;
using referecence = typename Base::reference;
public:
/// Constructor
GridView(GridBase