diff --git a/src/core/fftransform.hh b/src/core/fftransform.hh index 2cf8fa0..586d0f7 100644 --- a/src/core/fftransform.hh +++ b/src/core/fftransform.hh @@ -1,115 +1,115 @@ /** * * @author Lucas Frérot * * @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 . * */ /* -------------------------------------------------------------------------- */ #ifndef FFTRANSFORM_HH #define FFTRANSFORM_HH /* -------------------------------------------------------------------------- */ #include "tamaas.hh" #include "grid.hh" #include "grid_hermitian.hh" #include "types.hh" #include /* -------------------------------------------------------------------------- */ __BEGIN_TAMAAS__ /* -------------------------------------------------------------------------- */ template class FFTransform { public: /// Constructor FFTransform(Grid & real, GridHermitian & spectral); /// Destructor virtual ~FFTransform(); public: /// Perform FFT virtual void forwardTransform() = 0; /// Perform IFFT virtual void backwardTransform() = 0; /// Normalize the real surface after backward transform virtual void normalize(); public: /// Fill a grid with modevectors: boolean if grid has hermitian dimensions template static void computeFrequencies(Grid & freq); protected: Grid & real; GridHermitian & spectral; }; /* -------------------------------------------------------------------------- */ template template void FFTransform::computeFrequencies(Grid & freq) { // If hermitian is true, we suppose the dimensions of freq are // reduced based on hermitian symetry and that it has dim components auto n = freq.sizes(); #pragma omp parallel for for (auto it = freq.begin(dim) ; it < freq.end(dim) ; ++it) { - VectorProxy wavevector(it, dim); + VectorProxy wavevector(&(*it), dim); UInt index = freq.getNbPoints() - (freq.end(dim) - it); std::array tuple = {0}; /// Computing tuple from index for (UInt d = dim-1 ; d > 0 ; d--) { tuple[d] = index % n[d]; index -= tuple[d]; index /= n[d]; } tuple[0] = index; UInt dmax = dim; if (hermitian) { dmax = dim - 1; wavevector(dim-1) = tuple[dim-1]; } for (UInt d = 0 ; d < dmax ; d++) { T td = tuple[d]; T nd = n[d]; T q = (tuple[d] <= n[d]/2) ? td : td-nd; wavevector(d) = q; } } } /* -------------------------------------------------------------------------- */ __END_TAMAAS__ /* -------------------------------------------------------------------------- */ #endif // FFTRANSFORM_HH diff --git a/src/core/grid.cpp b/src/core/grid.cpp index 9cd3e37..f45eb26 100644 --- a/src/core/grid.cpp +++ b/src/core/grid.cpp @@ -1,167 +1,156 @@ /** * * @author Lucas Frérot * * @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 . * */ /* -------------------------------------------------------------------------- */ #include "tamaas.hh" #include "grid.hh" #include #include #include /* -------------------------------------------------------------------------- */ __BEGIN_TAMAAS__ /* -------------------------------------------------------------------------- */ template Grid::Grid(): GridBase() { this->n.fill(0); this->strides.fill(1); this->nb_components = 1; } /* -------------------------------------------------------------------------- */ template Grid::Grid(const std::array & n, UInt nb_components): GridBase() { this->init(n, nb_components); } template Grid::Grid(const std::vector & n, UInt nb_components): GridBase() { if (n.size() != dim) TAMAAS_EXCEPTION("Provided sizes for grid do not match dimension"); this->init(n, nb_components); } template Grid::Grid(const std::initializer_list & n, UInt nb_components): GridBase() { if (n.size() != dim) TAMAAS_EXCEPTION("Provided sizes for grid do not match dimension"); this->init(n, nb_components); } template Grid::Grid(const Grid & o): GridBase(o), n(o.n), strides(o.strides) {} template Grid::Grid(Grid&& o): GridBase(o), n(std::move(o.n)), strides(std::move(o.strides)) {} /* -------------------------------------------------------------------------- */ template void Grid::resize(const std::array & n) { this->n = n; UInt size = this->computeSize(); this->data.resize(size); this->data.assign(size, T(0.)); this->computeStrides(); } /* -------------------------------------------------------------------------- */ template void Grid::computeStrides() { std::copy(n.begin()+1, n.end(), strides.rbegin()+2); strides[dim] = 1; strides[dim-1] = this->nb_components; std::partial_sum(strides.rbegin(), strides.rend(), strides.rbegin(), std::multiplies()); } /* -------------------------------------------------------------------------- */ -template -void Grid::uniformSetComponents(const Grid &vec) { - TAMAAS_ASSERT(vec.dataSize() == this->nb_components, - "Cannot set grid field with values of vector"); - -#pragma omp parallel for - for (UInt i = 0 ; i < dataSize() ; i++) { - this->data[i] = vec(i % this->nb_components); - } -} - -/* -------------------------------------------------------------------------- */ - template void Grid::printself(std::ostream & str) const { str << "Grid(" << dim << ", " << this->nb_components << ") {"; - for (UInt i = 0 ; i < this->data.size() - 1 ; i++) { - str << this->data[i] << ", "; + const auto size = this->dataSize(); + for (UInt i = 0 ; i < size ; i++) { + str << (*this)(i) << ", "; } - str << this->data[this->data.size()-1] << "}"; + str << "\b\b}"; } /* -------------------------------------------------------------------------- */ #define GRID_SCALAR_OPERATOR_IMPL(op) \ template \ inline Grid & Grid::operator op(const T & e) { \ + const auto size = this->dataSize(); \ _Pragma("omp parallel for") \ - for (UInt i = 0 ; i < this->data.size() ; i++) { \ - this->data[i] op e; \ + for (UInt i = 0 ; i < size ; i++) { \ + (*this)(i) op e; \ } \ return *this; \ } \ GRID_SCALAR_OPERATOR_IMPL(+=); GRID_SCALAR_OPERATOR_IMPL(*=); GRID_SCALAR_OPERATOR_IMPL(-=); GRID_SCALAR_OPERATOR_IMPL(/=); GRID_SCALAR_OPERATOR_IMPL(=); #undef GRID_SCALAR_OPERATOR_IMPL /* -------------------------------------------------------------------------- */ /// Class instanciation #define GRID_INSTANCIATE_TYPE(type) template class Grid; \ template class Grid; \ template class Grid GRID_INSTANCIATE_TYPE(Real); GRID_INSTANCIATE_TYPE(UInt); GRID_INSTANCIATE_TYPE(Complex); GRID_INSTANCIATE_TYPE(int); GRID_INSTANCIATE_TYPE(bool); GRID_INSTANCIATE_TYPE(unsigned long); #undef GRID_INSTANCIATE_TYPE __END_TAMAAS__ diff --git a/src/core/grid.hh b/src/core/grid.hh index 0d08194..9c17edb 100644 --- a/src/core/grid.hh +++ b/src/core/grid.hh @@ -1,195 +1,212 @@ /** * * @author Lucas Frérot * * @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 . * */ /* -------------------------------------------------------------------------- */ #ifndef __GRID_HH__ #define __GRID_HH__ /* -------------------------------------------------------------------------- */ #include "tamaas.hh" #include "grid_base.hh" #include #include #include #include /* -------------------------------------------------------------------------- */ __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 class Grid: public GridBase { public: /* -------------------------------------------------------------------------- */ /* Types */ /* -------------------------------------------------------------------------- */ typedef T value_type; typedef typename GridBase::iterator iterator; static constexpr UInt dimension = dim; /* -------------------------------------------------------------------------- */ /* Constructors */ /* -------------------------------------------------------------------------- */ public: /// Constructor by default (empty array) Grid(); /// Constructor Grid(const std::array & n, UInt nb_components); /// Constructor with vectors Grid(const std::vector & n, UInt nb_components); /// Constructor with initializer list Grid(const std::initializer_list & n, UInt nb_components); /// Copy constructor Grid(const Grid & o); /// Move constructor (transfers data ownership) Grid(Grid&& o); private: /// Init from standard container - template - void init(const T1 & n, UInt nb_components); + template + void init(const Container & n, UInt nb_components); public: /* -------------------------------------------------------------------------- */ /* Common operations */ /* -------------------------------------------------------------------------- */ /// Resize array void resize(const std::array & n); /// Compute size inline UInt computeSize() const; /// Compute strides virtual void computeStrides(); /// Print void printself(std::ostream & str) const; /// Get sizes const std::array & sizes() const {return n;} - /// Get total size - UInt dataSize() const {return this->data.size(); } - /// Get number of points - UInt getNbPoints() const {return dataSize() / this->getNbComponents(); } - /// Set components - void uniformSetComponents(const Grid & vec); /// Get strides const std::array & getStrides() const {return this->strides;} /* -------------------------------------------------------------------------- */ /* Access operators */ /* -------------------------------------------------------------------------- */ /// Variadic access operator (non-const) template inline T & operator()(T1... args); /// Variadic access operator template inline const T & operator()(T1... args) const; /// Tuple index access operator inline T & operator()(std::array tuple); 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 inline UInt computeOffset(std::array tuple) const; /* -------------------------------------------------------------------------- */ /* Arithmetic operators */ /* -------------------------------------------------------------------------- */ public: #define GRID_VEC_OPERATOR(op) \ template \ - Grid & operator op(const Grid & other) \ + Grid & operator op(const Grid & other) \ GRID_VEC_OPERATOR(+=); GRID_VEC_OPERATOR(*=); GRID_VEC_OPERATOR(-=); GRID_VEC_OPERATOR(/=); #undef GRID_VEC_OPERATOR -#define GRID_SCALAR_OPERATOR(op) \ +#define GRID_SCALAR_OPERATOR(op) \ Grid & operator op(const T & e) \ GRID_SCALAR_OPERATOR(+=); GRID_SCALAR_OPERATOR(*=); GRID_SCALAR_OPERATOR(-=); GRID_SCALAR_OPERATOR(/=); GRID_SCALAR_OPERATOR(=); #undef GRID_SCALAR_OPERATOR +/* -------------------------------------------------------------------------- */ +/* Move/Copy operators */ +/* -------------------------------------------------------------------------- */ + // = operator Grid & operator=(const Grid & other); // = operator (not const input, otherwise gcc is confused) Grid & operator=(Grid & other); // = operator (move) Grid & operator=(Grid && other); // Copy data from another grid template void copy(const Grid & other); // Move data from another grid template void move(Grid && other); +/* -------------------------------------------------------------------------- */ +/* Iterators */ +/* -------------------------------------------------------------------------- */ +public: + virtual iterator begin(UInt n = 1) { + std::vector strides(this->strides.begin(), this->strides.end()); + std::vector sizes(this->n.begin(), this->n.end()); + sizes.push_back(this->nb_components); + return iterator(this->getInternalData()+this->offset, 0, n, strides, sizes); + } + + virtual iterator end(UInt n = 1) { + std::vector strides(this->strides.begin(), this->strides.end()); + std::vector sizes(this->n.begin(), this->n.end()); + sizes.push_back(this->nb_components); + return iterator(this->getInternalData()+this->offset, + this->dataSize(), n, strides, sizes); + } + /* -------------------------------------------------------------------------- */ /* Member variables */ /* -------------------------------------------------------------------------- */ protected: std::array n; std::array strides; }; __END_TAMAAS__ /* -------------------------------------------------------------------------- */ /* Inline/template function definitions */ /* -------------------------------------------------------------------------- */ #include "grid_tmpl.hh" /* -------------------------------------------------------------------------- */ #endif // __GRID_HH__ diff --git a/src/core/grid_base.hh b/src/core/grid_base.hh index 9a04077..9b77ec6 100644 --- a/src/core/grid_base.hh +++ b/src/core/grid_base.hh @@ -1,144 +1,198 @@ /** * * @author Lucas Frérot * * @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 . * */ /* -------------------------------------------------------------------------- */ #ifndef __GRID_BASE_HH__ #define __GRID_BASE_HH__ /* -------------------------------------------------------------------------- */ #include "tamaas.hh" #include "array.hh" #include +#include +#include /* -------------------------------------------------------------------------- */ __BEGIN_TAMAAS__ /// Namespace for view index classes namespace view { /// Index class struct index { index(Int i = -1):i(i){} Int i; }; /// Blocked index class struct blocked : public index { blocked(Int i):index(i){} }; /// Free index class struct free : public index { free():index(-1){} }; } template class GridBase { public: /// Constructor by default GridBase() = default; /// Copy constructor GridBase(const GridBase & o): data(o.data), nb_components(o.nb_components), offset(o.offset) {} /// Move constructor (transfers data ownership) GridBase(GridBase&& o): data(std::move(o.data)), - nb_components(o.nb_components), - offset(o.offset) {} + nb_components(std::exchange(o.nb_components, 1)), + offset(std::exchange(o.offset, 0)) {} /// Destructor virtual ~GridBase() = default; /* -------------------------------------------------------------------------- */ /* Iterator class */ /* -------------------------------------------------------------------------- */ protected: struct iterator { /// constructor - iterator(T * start, ptrdiff_t offset):ptr(start), offset(offset) {} + iterator(T * data, + std::size_t start, + ptrdiff_t offset, + std::vector strides, + std::vector sizes): + data(data), + index(start), + offset(offset), + strides(strides), + n(sizes) {} /// destructor ~iterator() {} /// pre-increment - inline iterator & operator++() { ptr += offset; return *this; } + inline iterator & operator++() { index += offset; return *this; } /// comparison - inline bool operator<(const iterator & a) { return ptr < a.ptr; } + inline bool operator<(const iterator & a) { return index < a.index; } /// inequality - inline bool operator!=(const iterator & a) { return ptr != a.ptr; } + inline bool operator!=(const iterator & a) { return index != a.index; } /// increment with given offset - inline iterator & operator+=(ptrdiff_t a) { ptr += a*offset; return *this;} + inline iterator & operator+=(ptrdiff_t a) { index += a*offset; return *this;} /// dereference iterator - inline T & operator*() { return *ptr; } + inline T & operator*() { + std::vector tuple(strides.size()); + UInt index_copy = index; + + tuple.back() = index_copy % n.back(); + index_copy -= tuple.back(); + index_copy /= n.back(); + + /// Computing tuple from index + for (UInt d = tuple.size()-2 ; d > 0 ; d--) { + tuple[d] = index_copy % this->n[d]; + index_copy -= tuple[d]; + index_copy /= this->n[d]; + } + + tuple[0] = index_copy; + + UInt access_offset = 0; + for (UInt i = 0 ; i < tuple.size() ; ++i) + access_offset += tuple[i] * strides[i]; + return data[access_offset]; + } /// needed for OpenMP range calculations - inline Int operator-(const iterator & a) const { return (ptr - a.ptr)/offset; } + inline Int operator-(const iterator & a) const { return (index - a.index)/offset; } protected: - T * ptr; + T * data; + std::size_t index; ptrdiff_t offset; + std::vector strides; + std::vector n; }; -public: - /// Begin iterator with n components - virtual iterator begin(UInt n = 1) { return iterator(data.getPtr(), n); } - /// End iterator with n components - virtual iterator end(UInt n = 1) { return iterator(data.getPtr()+computeSize(), n); } - public: /// Compute size virtual UInt computeSize() const = 0; /// Get internal data pointer (const) const T * getInternalData() const {return this->data.getPtr();} /// Get internal data pointer (non-const) T * getInternalData() {return this->data.getPtr();} /// Get number of components UInt getNbComponents() const {return nb_components;} /// Set number of components void setNbComponents(UInt n) {nb_components = n;} /// Get offset UInt getOffset() const {return offset;} + /// Get total size + virtual UInt dataSize() const {return this->data.size(); } + /// Get number of points + UInt getNbPoints() const {return this->computeSize() / this->getNbComponents(); } + /// Set components + void uniformSetComponents(const GridBase & vec); + +/* -------------------------------------------------------------------------- */ +/* Iterator methods */ +/* -------------------------------------------------------------------------- */ + virtual iterator begin(UInt n = 1) = 0; + virtual iterator end(UInt n = 1) = 0; public: template void copy(const GridBase & other) { data = other.data; nb_components = other.nb_components; offset = other.offset; } template void move(GridBase && other) { data = std::move(other.data); - nb_components = other.nb_components; - offset = other.offset; - other.nb_components = 0; - other.offset = 0; + nb_components = std::exchange(other.nb_components, 1); + offset = std::exchange(other.offset, 0); } protected: Array data; UInt nb_components = 1; UInt offset = 0; }; +/* -------------------------------------------------------------------------- */ + +template +void GridBase::uniformSetComponents(const GridBase &vec) { + if (!(vec.dataSize() == this->nb_components)) + TAMAAS_EXCEPTION("Cannot set grid field with values of vector"); + + const auto size = this->dataSize(); +#pragma omp parallel for + for (UInt i = 0 ; i < size ; i++) { + this->data[i] = vec.data[i % this->nb_components]; + } +} + + __END_TAMAAS__ #endif // __GRID_BASE_HH__ diff --git a/src/core/grid_view.hh b/src/core/grid_view.hh index 9c53643..9b1b1d7 100644 --- a/src/core/grid_view.hh +++ b/src/core/grid_view.hh @@ -1,156 +1,158 @@ /** * * @author Lucas Frérot * * @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 . * */ /* -------------------------------------------------------------------------- */ #ifndef __GRID_VIEW_HH__ #define __GRID_VIEW_HH__ /* -------------------------------------------------------------------------- */ #include "tamaas.hh" #include "grid.hh" #include /* -------------------------------------------------------------------------- */ __BEGIN_TAMAAS__ /* -------------------------------------------------------------------------- */ -template -class GridView : public Base { - typedef typename Base::value_type T; - typedef typename Base::iterator iterator; +template