Homec4science
Diffusion tamaas (maxwell_viscoelastic-branch) Browse

tamaas/src/core/grid_base.hh
393b17a3c8bbmaxwell_viscoelastic…

Last Change
Raw File

grid_base.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_BASE_HH__
#define __GRID_BASE_HH__
/* -------------------------------------------------------------------------- */
#include "tamaas.hh"
#include "array.hh"
#include "iterator.hh"
#include <cstddef>
#include <utility>
#include <vector>
#include <limits>
/* -------------------------------------------------------------------------- */
__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 <typename T>
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(std::exchange(o.nb_components, 1)),
offset(std::exchange(o.offset, 0)) {}
/// Destructor
virtual ~GridBase() = default;
/* -------------------------------------------------------------------------- */
/* Iterator class */
/* -------------------------------------------------------------------------- */
public:
using iterator = iterator_::iterator<T>;
using const_iterator = iterator_::iterator<const T>;
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<T> & vec);
/* -------------------------------------------------------------------------- */
/* Iterator methods */
/* -------------------------------------------------------------------------- */
virtual iterator begin(UInt n = 1) = 0;
virtual iterator end(UInt n = 1) = 0;
virtual const_iterator begin(UInt n = 1) const = 0;
virtual const_iterator end(UInt n = 1) const = 0;
/* -------------------------------------------------------------------------- */
/* Operators */
/* -------------------------------------------------------------------------- */
#define VEC_OPERATOR(op) \
template <typename T1> \
GridBase & operator op(const GridBase<T1> & other) \
VEC_OPERATOR(+=);
VEC_OPERATOR(*=);
VEC_OPERATOR(-=);
VEC_OPERATOR(/=);
#undef VEC_OPERATOR
template <typename T1>
T dot(const GridBase<T1> & other);
#define SCALAR_OPERATOR(op) \
GridBase & operator op(const T & e) \
SCALAR_OPERATOR(+=);
SCALAR_OPERATOR(*=);
SCALAR_OPERATOR(-=);
SCALAR_OPERATOR(/=);
SCALAR_OPERATOR(=);
#undef SCALAR_OPERATOR
virtual GridBase & operator=(const GridBase & o);
/* -------------------------------------------------------------------------- */
/* Min/max */
/* -------------------------------------------------------------------------- */
T min() const;
T max() const;
T sum() const;
T mean() const { return sum() / dataSize(); }
T var() const;
/* -------------------------------------------------------------------------- */
/* Move/Copy */
/* -------------------------------------------------------------------------- */
public:
template <typename T1>
void copy(const GridBase<T1> & other) {
data = other.data;
nb_components = other.nb_components;
offset = other.offset;
}
template <typename T1>
void move(GridBase<T1> && other) {
data = std::move(other.data);
nb_components = std::exchange(other.nb_components, 1);
offset = std::exchange(other.offset, 0);
}
protected:
Array<T> data;
UInt nb_components = 1;
UInt offset = 0;
};
/* -------------------------------------------------------------------------- */
template <typename T>
void GridBase<T>::uniformSetComponents(const GridBase<T> &vec) {
if (!(vec.dataSize() == this->nb_components))
TAMAAS_EXCEPTION("Cannot set grid field with values of vector");
auto begin_it = begin(), end_it = end();
#pragma omp parallel for
for (auto it = begin() ; it < end_it ; ++it) {
UInt i = it - begin_it;
*it = vec.data[i % this->nb_components];
}
}
/* -------------------------------------------------------------------------- */
#define SCALAR_OPERATOR_IMPL(op) \
template <typename T> \
inline GridBase<T> & GridBase<T>::operator op(const T & e) { \
_Pragma("omp parallel for") \
for (auto it = this->begin() ; it < this->end() ; ++it) { \
*it op e; \
} \
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 <typename T> \
template <typename T1> \
inline GridBase<T> & GridBase<T>::operator op(const GridBase<T1> & other) { \
TAMAAS_ASSERT(other.computeSize() == this->computeSize(), \
"surface size does not match"); \
auto begin_it = this->begin(), end_it = this->end(); \
auto other_begin = other.begin(), other_it = other.begin(); \
_Pragma("omp parallel for firstprivate(other_it)") \
for (auto it = begin() ; it < end_it ; ++it) { \
other_it = other_begin; \
other_it += it - begin_it; \
*it op *other_it; \
} \
return *this; \
} \
VEC_OPERATOR_IMPL(+=);
VEC_OPERATOR_IMPL(-=);
VEC_OPERATOR_IMPL(*=);
VEC_OPERATOR_IMPL(/=);
#undef VEC_OPERATOR
/* -------------------------------------------------------------------------- */
template <typename T>
GridBase<T> & GridBase<T>::operator=(const GridBase & o) {
this->copy(o);
return *this;
}
/* -------------------------------------------------------------------------- */
template <typename T>
inline T GridBase<T>::min() const {
T val = std::numeric_limits<T>::max();
auto begin_it = begin(), end_it = end();
#pragma omp parallel for reduction(min:val)
for (auto it = begin_it ; it < end_it ; ++it) {
val = std::min(*it, val);
}
return val;
}
/* -------------------------------------------------------------------------- */
template <typename T>
inline T GridBase<T>::max() const {
T val = std::numeric_limits<T>::lowest();
auto begin_it = begin(), end_it = end();
#pragma omp parallel for reduction(max:val)
for (auto it = begin_it ; it < end_it ; ++it) {
val = std::max(*it, val);
}
return val;
}
/* -------------------------------------------------------------------------- */
template <typename T>
inline T GridBase<T>::sum() const {
T val = 0;
auto begin_it = begin(), end_it = end();
#pragma omp parallel for reduction(+:val)
for (auto it = begin_it ; it < end_it ; ++it) {
val += *it;
}
return val;
}
/* -------------------------------------------------------------------------- */
template <typename T>
inline T GridBase<T>::var() const {
T var = 0;
const T mu = mean();
auto begin_it = begin(), end_it = end();
#pragma omp parallel for reduction(+:var)
for (auto it = begin_it ; it < end_it ; ++it) {
var += (*it - mu) * (*it - mu);
}
return var / (dataSize() - 1);
}
/* -------------------------------------------------------------------------- */
template <typename T>
template <typename T1>
T GridBase<T>::dot(const GridBase<T1> & other) {
T res = 0;
auto begin_it = begin(), end_it = end();
auto other_b = other.begin(), other_it = other.begin();
#pragma omp parallel for firstprivate(other_it) reduction(+: res)
for (auto it = begin_it ; it < end_it ; ++it) {
other_it = other_b;
other_it += it - begin_it;
res += *it * *other_it;
}
return res;
}
__END_TAMAAS__
#endif // __GRID_BASE_HH__
c4science · Help