Homec4science
Diffusion LibMultiScale (features/new_parralelism) Browse

libmultiscale/src/math/tensor.hh
8ae23a91b711features/new_parralelism

Last Change
Raw File

tensor.hh
View Options

/**
* @file tensor.hh
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
*
* @date Wed May 01 17:45:41 2013
*
* @brief This is a tensor representation
*
* @section LICENSE
*
* Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* LibMultiScale 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.
*
* LibMultiScale 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 LibMultiScale. If not, see <http://www.gnu.org/licenses/>.
*
*/
#ifndef __LIBMULTISCALE_TENSOR_HH__
#define __LIBMULTISCALE_TENSOR_HH__
/* -------------------------------------------------------------------------- */
#include "lm_common.hh"
#include <cmath>
/* -------------------------------------------------------------------------- */
__BEGIN_LIBMULTISCALE__
template <typename T, UInt Dim1, UInt Dim2, typename Daughter> class TensorTmp {
public:
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
TensorTmp() : TensorTmp<T, Dim1, Dim2, Daughter>(0.0){};
TensorTmp(Real val) {
for (UInt i = 0; i < Dim1 * Dim2; ++i)
vec[i] = val;
};
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
inline Daughter pow(UInt p) {
Daughter res;
for (UInt i = 0; i < Dim1 * Dim2; ++i)
res[i] = std::pow(vec[i], p);
return res;
}
inline Daughter fabs() {
Daughter res;
for (UInt i = 0; i < Dim1 * Dim2; ++i)
res[i] = std::fabs(vec[i]);
return res;
}
inline T sum() {
T res = T();
for (UInt i = 0; i < Dim1 * Dim2; ++i)
res += vec[i];
return res;
}
inline Daughter operator*(const Daughter &q) const {
Daughter tmp;
for (UInt i = 0; i < Dim1 * Dim2; ++i) {
tmp[i] = vec[i] * q[i];
}
return tmp;
}
inline Daughter operator+(const Daughter &q) const {
Daughter tmp;
for (UInt i = 0; i < Dim1 * Dim2; ++i) {
tmp[i] = vec[i] + q[i];
}
return tmp;
}
inline Daughter operator*(T v) const {
Daughter tmp;
for (UInt i = 0; i < Dim1 * Dim2; ++i) {
tmp[i] = vec[i] * v;
}
return tmp;
}
inline void operator+=(Daughter &t) {
for (UInt i = 0; i < Dim1 * Dim2; ++i)
vec[i] += t[i];
}
inline void operator-=(Daughter &t) {
for (UInt i = 0; i < Dim1 * Dim2; ++i)
vec[i] = t[i];
}
inline Daughter operator-(Daughter &t) {
Daughter tmp;
for (UInt i = 0; i < Dim1 * Dim2; ++i)
tmp[i] = vec[i] - t[i];
return tmp;
}
inline void operator/=(T v) {
for (UInt i = 0; i < Dim1 * Dim2; ++i)
vec[i] /= v;
}
inline void operator*=(T v) {
for (UInt i = 0; i < Dim1 * Dim2; ++i)
vec[i] *= v;
}
inline void operator*=(Daughter v) {
for (UInt i = 0; i < Dim1 * Dim2; ++i)
vec[i] *= v[i];
}
inline void square() {
for (UInt i = 0; i < Dim1 * Dim2; ++i)
vec[i] = vec[i] * vec[i];
}
inline Real &operator()(UInt i, UInt j) { return vec[i * Dim1 + j]; }
inline Real operator()(UInt i, UInt j) const { return vec[i * Dim1 + j]; }
inline Real &operator[](UInt i) { return vec[i]; }
inline Real operator[](UInt i) const { return vec[i]; }
inline Daughter &operator=(Real v) {
for (UInt i = 0; i < Dim1 * Dim2; ++i)
vec[i] = v;
return *this;
}
inline Real *getVec() { return vec; };
inline Real size() { return Dim1 * Dim2; };
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
protected:
T vec[Dim1 * Dim2];
};
/* -------------------------------------------------------------------------- */
template <typename T, UInt Dim1, UInt Dim2, typename Daughter>
inline std::ostream &operator<<(std::ostream &os,
TensorTmp<T, Dim1, Dim2, Daughter> &t) {
for (UInt i = 0; i < 9; ++i) {
os << t[i] << " ";
}
return os;
}
/* -------------------------------------------------------------------------- */
class Tensor : public TensorTmp<Real, 3, 3, Tensor> {
using TensorTmp<Real, 3, 3, Tensor>::TensorTmp;
};
/* -------------------------------------------------------------------------- */
template <typename T, UInt Dim1, UInt Dim2, typename Daughter>
inline Daughter operator*(const Real &scalar,
TensorTmp<T, Dim1, Dim2, Daughter> &q) {
Daughter tmp(q);
return tmp * scalar;
}
__END_LIBMULTISCALE__
namespace std {
template <typename T, ::libmultiscale::UInt Dim1, ::libmultiscale::UInt Dim2,
typename Daughter>
const inline ::libmultiscale::TensorTmp<T, Dim1, Dim2, Daughter> &
max(const ::libmultiscale::TensorTmp<T, Dim1, Dim2, Daughter> &a,
const ::libmultiscale::TensorTmp<T, Dim1, Dim2, Daughter> &b) {
::libmultiscale::Real na = 0., nb = 0.;
for (::libmultiscale::UInt i = 0; i < Dim1 * Dim2; ++i) {
na += a[i] * a[i];
nb += b[i] * b[i];
}
if (na > nb)
return a;
else
return b;
}
template <typename T, ::libmultiscale::UInt Dim1, ::libmultiscale::UInt Dim2,
typename Daughter>
const inline ::libmultiscale::TensorTmp<T, Dim1, Dim2, Daughter> &
min(const ::libmultiscale::TensorTmp<T, Dim1, Dim2, Daughter> &a,
const ::libmultiscale::TensorTmp<T, Dim1, Dim2, Daughter> &b) {
::libmultiscale::Real na = 0., nb = 0.;
for (::libmultiscale::UInt i = 0; i < Dim1 * Dim2; ++i) {
na += a[i] * a[i];
nb += b[i] * b[i];
}
if (na > nb)
return b;
else
return a;
}
}
/* -------------------------------------------------------------------------- */
#endif /* __LIBMULTISCALE_TENSOR_HH__ */
c4science · Help