Homec4science
Diffusion µSpectre prototype implementation (add_eigenstrain) Browse

muSpectre_proto/src/common/units.hh
dca58adc2332add_eigenstrain

Last Change
Raw File

units.hh
View Options

/**
* file units.hh
*
* @author Till Junge <till.junge@epfl.ch>
*
* @date 02 May 2017
*
* @brief Implements physical units without runtime overhead. This is inspi-
* red by section 28.7 of Stroustrup's "The C++ Programming Language",
* but extended to rational powers of basic units (fracture mechanics!)
*
* @section LICENCE
*
* Copyright (C) 2017 Till Junge
*
* µSpectre is free software; you can redistribute it and/or
* modify it under the terms of the GNU General Public License as
* published by the Free Software Foundation, either version 3, or (at
* your option) any later version.
*
* µSpectre 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
* General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with GNU Emacs; see the file COPYING. If not, write to the
* Free Software Foundation, Inc., 59 Temple Place - Suite 330,
* Boston, MA 02111-1307, USA.
*/
#include <ratio>
#include <limits>
#include <iostream>
#include <boost/preprocessor/seq/for_each.hpp>
#include "common.hh"
#include <unsupported/Eigen/CXX11/Tensor>
#include <Eigen/src/Core/functors/BinaryFunctors.h>
#include <Eigen/Dense>
#ifndef UNITS_H
#define UNITS_H
namespace muSpectre {
namespace units {
//----------------------------------------------------------------------------//
template <typename m__ratio, typename kg_ratio, typename s__ratio>
struct unit {
using m = m__ratio;
using kg = kg_ratio;
using s = s__ratio;
};
//----------------------------------------------------------------------------//
template <typename unit1, typename unit2>
struct uplus {
using type = unit<std::ratio_add<typename unit1::m, typename unit2::m >,
std::ratio_add<typename unit1::kg, typename unit2::kg>,
std::ratio_add<typename unit1::s, typename unit2::s >>;
};
//----------------------------------------------------------------------------//
template<typename unit1, typename unit2>
using unit_plus = typename uplus<unit1, unit2>::type;
//----------------------------------------------------------------------------//
template <typename unit1, typename unit2>
struct uminus {
using type = unit<std::ratio_subtract<typename unit1::m, typename unit2::m >,
std::ratio_subtract<typename unit1::kg, typename unit2::kg>,
std::ratio_subtract<typename unit1::s, typename unit2::s >>;
};
//----------------------------------------------------------------------------//
template<typename unit1, typename unit2>
using unit_minus = typename uminus<unit1, unit2>::type;
using ratio1 = std::ratio<1>;
using ratio0 = std::ratio<0>;
using nodim = unit<ratio0, ratio0, ratio0>;
using x = unit<ratio0, ratio0, ratio0>;
using m = unit<ratio1, ratio0, ratio0>;
using kg = unit<ratio0, ratio1, ratio0>;
using s = unit<ratio0, ratio0, ratio1>;
using si_speed = unit_minus<m, s>;
using si_acceleration = unit_minus<si_speed, s>;
using si_area = unit_plus<m, m>;
using N = unit_plus<kg, si_acceleration>;
using Pa = unit_minus<N, si_area>;
using J = unit_plus<N,m>;
} // units
//----------------------------------------------------------------------------//
template <typename T, typename U>
struct Quantity {
T val;
explicit constexpr Quantity(T r): val{r} {}
explicit constexpr Quantity(): val(){}
template<typename T2>
inline Quantity& operator=(const Quantity<T2, U> & other) {
this->val = other.val; return *this;}
inline auto eval(){return this->val.eval();}
};
// //----------------------------------------------------------------------------//
// template<typename T, typename U, size_t Rows>
// struct Quantity<Eigen::Matrix<T, Rows, Eigen::Dynamic>, U> {
// Eigen::Matrix<T, Rows, Eigen::Dynamic> val;
// Quantity(int dynrows, int dyncols): val(dynrows, dyncols) {}
// constexpr Quantity() = default;
// };
//----------------------------------------------------------------------------//
template<>
template<typename T2, typename U, int Rows, int Cols>
struct Quantity<Eigen::Matrix<T2, Rows, Cols>, U> {
using val_t = Eigen::Matrix<T2, Rows, Cols>;
val_t val;
explicit Quantity(int dynrows, int dyncols): val(dynrows, dyncols) {}
explicit Quantity(const val_t& m): val(m) {}
explicit Quantity(const val_t&& m): val(m) {}
constexpr Quantity():val(){};
template<typename Tother>
inline Quantity& operator=(const Quantity<Tother, U> & other) {
this->val = other.val; return *this;}
};
//----------------------------------------------------------------------------//
template<typename T, typename U>
std::ostream & operator<<(std::ostream & os, const Quantity<T, U> Q) {
os << Q.val
<< "m^" << U::m ::num/U::m ::den
<< "kg^" << U::kg::num/U::kg::den
<< "s^" << U::s ::num/U::s ::den;
return os;
}
//----------------------------------------------------------------------------//
template <typename T1, typename T2, typename U>
bool operator==(const Quantity<T1,U>& x, const Quantity<T2,U>& y) {
return x.val == y.val;
}
//----------------------------------------------------------------------------//
template <typename T1, typename T2, typename U>
auto operator+(const Quantity<T1,U>& x, const Quantity<T2,U> & y) {
return Quantity<decltype(x.val+y.val),U>(x.val + y.val);
}
//----------------------------------------------------------------------------//
template <typename T, typename U>
Quantity<T,U>& operator+=(Quantity<T,U>& x, const Quantity<T,U>& y) {
x.val += y.val;
return x;
}
//----------------------------------------------------------------------------//
template <typename T, typename U>
Quantity<T,U> operator-(const Quantity<T,U> & x, const Quantity<T,U> & y) {
return Quantity<T,U>(x.val - y.val);
}
//----------------------------------------------------------------------------//
template <typename T, typename U>
Quantity<T,U>& operator-=(Quantity<T,U> & x, const Quantity<T,U> & y) {
x.val -= y.val;
return x;
}
//----------------------------------------------------------------------------//
template <typename T, typename U, typename T2>
Quantity<T,U> operator*(Quantity<T,U> x, T2 y) {
return Quantity<T,U>(x.val * y);
}
//----------------------------------------------------------------------------//
template <typename T1, typename T2, typename U1, typename U2>
auto operator*(const Quantity<T1, U1> & x,
const Quantity<T2, U2> & y) {
return Quantity<decltype(x.val*y.val), units::unit_plus<U1, U2>>{x.val*y.val};
}
//----------------------------------------------------------------------------//
template <typename T, typename U, typename T2>
Quantity<T,U> operator*(T2 x, Quantity<T,U> y) {
return Quantity<T,U>{x * y.val};
}
//----------------------------------------------------------------------------//
template <typename T, typename U1, typename U2>
Quantity<T, units::unit_minus<U1, U2>> operator/(const Quantity<T, U1> & x,
const Quantity<T, U2> & y) {
return Quantity<T, units::unit_minus<U1, U2>>{x.val/y.val};
}
//----------------------------------------------------------------------------//
template <typename T, typename U>
auto operator/(Quantity<T,U> x, Real y) {
return Quantity<decltype(x.val/y),U>{x.val/y};
}
//----------------------------------------------------------------------------//
template <typename T, typename U>
Quantity<T,U> operator/(Real x, Quantity<T, U> y) {
return Quantity<T, units::nodim>{x}/y;
}
//----------------------------------------------------------------------------//
//! string literals for convenient use of units
#define LITERAL_declaration(literal_string) \
Quantity<Real, units::literal_string> operator""_##literal_string(long double r);\
Quantity<Real, units::literal_string> operator""_##literal_string(unsigned long long r);\
// Warning: no literals seem possible for complex scalars
#define LITERALS (m)(kg)(s)(N)(Pa)(J)(x)(nodim)
#define MACRO_declaration(r, dummy, literal_string) LITERAL_declaration(literal_string)
BOOST_PP_SEQ_FOR_EACH(MACRO_declaration, _, LITERALS)
} // muSpectre
/** Add num_traits for use in eigen arrays
*/
namespace Eigen {
/* ---------------------------------------------------------------------- */
template<>
template<typename T, typename U>
struct NumTraits<muSpectre::Quantity<T, U>>
: NumTraits<T>
{};
/* ---------------------------------------------------------------------- */
template<>
template<typename T, typename U>
struct ScalarBinaryOpTraits<muSpectre::Quantity<T, U>, muSpectre::Real, internal::scalar_product_op<muSpectre::Quantity<T, U>, muSpectre::Real>>
{
typedef muSpectre::Quantity<T,U> ReturnType;
};
/* ---------------------------------------------------------------------- */
template<>
template<typename T, typename U>
struct ScalarBinaryOpTraits<muSpectre::Real, muSpectre::Quantity<T, U>, internal::scalar_product_op<muSpectre::Real, muSpectre::Quantity<T, U>>>
{
typedef muSpectre::Quantity<T,U> ReturnType;
};
//----------------------------------------------------------------------------//
//! required for matrix multiplications with same units
template<>
template<typename U1>
struct ScalarBinaryOpTraits<muSpectre::Quantity<muSpectre::Real, U1>, muSpectre::Quantity<muSpectre::Real, U1>, internal::scalar_product_op<muSpectre::Quantity<muSpectre::Real, U1>, muSpectre::Quantity<muSpectre::Real, U1>>>
{
typedef muSpectre::Quantity<muSpectre::Real,muSpectre::units::unit_plus<U1, U1>> ReturnType;
};
//----------------------------------------------------------------------------//
//! required for matrix multiplications with differing units
template<>
template<typename T, typename U1, typename U2>
struct ScalarBinaryOpTraits<muSpectre::Quantity<T, U1>, muSpectre::Quantity<T, U2>, internal::scalar_product_op<muSpectre::Quantity<T, U1>, muSpectre::Quantity<T, U2>>>
{
typedef muSpectre::Quantity<T,muSpectre::units::unit_plus<U1, U2>> ReturnType;
};
/* ---------------------------------------------------------------------- */
template<>
template<typename T, typename U>
struct ScalarBinaryOpTraits<muSpectre::Quantity<T, U>, muSpectre::Real,
internal::scalar_quotient_op<muSpectre::Quantity<T, U>, muSpectre::Real>>
{
typedef muSpectre::Quantity<T,U> ReturnType;
};
/* ---------------------------------------------------------------------- */
template<>
template<typename T, typename U>
struct ScalarBinaryOpTraits<muSpectre::Real, muSpectre::Quantity<T, U>,
internal::scalar_quotient_op<muSpectre::Real, muSpectre::Quantity<T, U>>>
{
typedef muSpectre::Quantity<T,U> ReturnType;
};
//----------------------------------------------------------------------------//
//! required for matrix multiplications with same units
template<>
template<typename U1>
struct ScalarBinaryOpTraits<muSpectre::Quantity<muSpectre::Real, U1>, muSpectre::Quantity<muSpectre::Real, U1>,
internal::scalar_quotient_op<muSpectre::Quantity<muSpectre::Real, U1>, muSpectre::Quantity<muSpectre::Real, U1>>>
{
typedef muSpectre::Quantity<muSpectre::Real,muSpectre::units::nodim> ReturnType;
};
//----------------------------------------------------------------------------//
//! required for matrix multiplications with differing units
template<>
template<typename T, typename U1, typename U2>
struct ScalarBinaryOpTraits<muSpectre::Quantity<T, U1>, muSpectre::Quantity<T, U2>,
internal::scalar_quotient_op<muSpectre::Quantity<T, U1>, muSpectre::Quantity<T, U2>>>
{
typedef muSpectre::Quantity<T,muSpectre::units::unit_minus<U1, U2>> ReturnType;
};
} // Eigen
#endif /* UNITS_H */
c4science · Help