Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F82144761
units.hh
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Mon, Sep 9, 20:21
Size
12 KB
Mime Type
text/x-c
Expires
Wed, Sep 11, 20:21 (2 d)
Engine
blob
Format
Raw Data
Handle
20654024
Attached To
rMSPPROTO µSpectre prototype implementation
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 */
Event Timeline
Log In to Comment